The dataset for this research was collected by Ferreira et al. (2018) with the aim of assessing the biotherapeutic effect of chicken serum on rodents experimentally infected by Trypanosoma cruzi, considering clinical, immunological and parasitological parameters. The project for the execution of the experiments was approved by the Ethics Committee on Research Involving Animals of the State University of Maringá, Paraná, Brazil, under legal opinion CEUA 2401220716/2016.

Based on Ferreira et al. (2018), the experiments were run through blind, controlled and randomized assays. They used 57 male Swiss mice, aged 56 days, sourced from the central vivarium of the State University of Maringá. The mice were divided into treatment groups in a way that the initial mean weight of the animals in each group was approximately equal, and kept in cages with a maximum of 5 individuals; all groups were maintained under the same experimental conditions. With n being the number of individuals per group, and g the number of cages, the experimental groups (treatments) were characterized as follows:

· G1: (Non-infected control). Non-infected and non-treated animals (n=5)(g=1)

· G2: Animals treated with biotherapeutic chicken serum with 13cH dilution (n=13)^[1]

· G3: Animals treated with biotherapeutic chicken serum with 6cH dilution (n=13)(g=3)

· G4: (Infection control). The animals were infected and received no treatment (n=13)(g=3).

· G5: Animals treated with biotherapeutic chicken serum with 3cH dilution (n=13)(g=4).

Their body weight was monitored over 12 weeks, measured on a semi-analytical scale (Balance BEL®)

According to Currie and Durbán (2002), a P-spline is a combination of B-splines and a difference penalty function of the L order associated with the estimated coefficients of the B-splines, and the use of this tool helps reduce the flexibility of the B-splines, thus preventing an overfitting of the curve.

Authors Eilers and Marx (1996) present in their article an interesting result concerning modeling through the combination of B-splines; they state that the integral of the second derivative squared can be expressed as a quadratic function in the coefficients associated with the sum. Thus, the roughness of a curve is expressed by taking a base with several B-splines associated with a penalty, as follows:

in which [a,b] is the study interval, d is the order of the L difference operator for, B is the regression basis function, λ is the control parameter of the smoothing, and is a function that composes B-spline sums.

Generalized additive mixed models using P-splines

Over the last decade, advancements have been observed in the application of P-splines to the context of generalized additive models. This is due to the possibility of rewriting the linear term by summing unknown smooth functions, not parametrically specified. In this sphere, it is possible to assign peculiar random effects to subjects, generating a new class, called Generalized Additive Mixed Models.

Currie and Durbán (2002), Durbán, Harezlak, Wand, and Carroll (2005) and Durbán (2014) describe that the appropriate model for this situation is that which allows estimating the curves of specific differences between individuals using non-parametric smooth terms. In this model, a smoothing non-parametric gi function is incorporated, thus enabling, according to Ruppert, Wand, and Carroll (2003), an interaction between treatment groups with the continuous predictor. This model can be described as:

in which yij indicates the response variable, () = 1+2+, in which 1, 2 denote the random effects, and () is a smoothing function that explains the trajectory of each individual, and indicates the i-th sampling unit observed at the j-th instant, and the k index indicates the k-ith treatment group. Besides, the fxij function is estimated using P-splines, which can be rewritten in the context of mixed model; besides:

in which the specific curves of the ( ) sampling units have parametric components 1+2 and non-parametric components, with both being random.

Thus, model 3 can be rewritten in its matrixial form as:

The X matrix of fixed effects is expressed by:

in which Ti is a matrix that indicates if the i-th subject receives the k-th treatment.

The β parameter vector is defined by:

The Z matrix of random effects is defined by:

The vector of random effects is expressed by:

And the covariance matrix is written as:

According to Diggle, Heagerty, Liang, and Zeger (2002), correlation structures are used for modeling covariance matrices, referring to the dependence between observations; when it comes to mixed-effect models, they are used for modeling the dependence between errors within the group. Because the variance is related to components Z_jGZ_iT and R_i, one can describe the structures of G or R, or both. Littell, Pendergast, and Natarajan (2000) and Wolfinger (1993) mention that the covariance between two observations for the same individual depends on the length of the time interval between measurements, and that the variance is considered constant over time. Covariance can be described in terms of the variance, and correlations can be expressed as a function of the time interval between measurements.

Model selection

Choosing or selecting the model that best explains data is undoubtedly fundamental in any statistical analysis in the modeling context. There are numerous tests and criteria for model selection in the literature, with one of them being the Akaike information criterion (AIC). According to Akaike (1974), the idea behind this criterion is to compare models, expressed by the following equation:

with l() being the logarithm of the likelihood function applied to the estimator of maximum likelihood of θ, () and p the number of explanatory variables. Thus, the satisfactory model will the one with the lowest AIC statistical value.

To conduct the statistical analyses, only the 'Treatment groups' and 'Experiment follow-up time' variables were considered, that is, the Cage variable was not treated due to the study objective. Table 1 contains summarized statistical data on the 'weight' response variable for each one of the treatment types:

Table 1.
Summarized measures for each treatment.

Observing Table 1 allows concluding that the mean of the weights in relation to each one of the treatments is 45.49 g for the 'Non-infected control' treatment group, and 41.44 g for the individuals allocated in the homeopathic treatment with 3cH dilution.

It is also possible to see a difference as to the maximum and minimum values for the weights in each one of the treatment groups; the minimum weights of the individuals allocated in the homeopathic treatment group with 3cH dilution (G5) is 31.00 g, whereas the minimum weight of the individuals in the Control treatment group (G1) is 39.20 g; as for the maximum weight, it ranges from 55.00 g, representing the 'Infection control' group of individuals (G4), to 49.20 g, representing the group under homeopathic treatment with 6cH dilution (G3). Finally, about the coefficient of variation, there is a small oscillation, with values ranging from 0.076 g and 0.106 g, approximately.

Figure 1 displays the behavior of the response variable for each one of the treatments; the group under homeopathic treatment with 3cH dilution (G5) had an average weight level lower than that of the other treatment groups, while the non-infected control group had the highest average weight level compared to the other treatment groups.

Figure 1.
Boxplot for the response variable by treatment.

Figure 2 presents the behavior of the response variable as a function of time, with the weight showing an increase behavior over time and oscillating as of the sixth week. Another remarkable fact was to find potentially outlier measures in the first week and in the fourth week of monitoring.

Figure 2.
Boxplot for the response variable over time.

About the profile chart, Figure 3 shows an increase behavior as of the second week as for weight over time. An oscillating behavior is seen for the weights soon in the first week, a factor that can be explained by this being the beginning of the treatment; besides, the response variable seems to behave linearly from the second to the sixth week for each subject in relation to the treatment groups as time varies.

Figure 3.
Average profile for each subject observed over time.

Figure 4 shows dispersion charts by treatment; there may be some variability between these observations, justifying the possibility of using some type of smoothing that takes into account the variability found.

It is also possible to notice that the weight gain of the mice in each treatment group varies greatly between the intercept and the slope of the lines that have been fitted; this allows considering mixed models with random effects, be they independent or correlated between intercept and slope in relation to time – in this case, for the Time covariable.

Figure 4.
Dispersion chart by treatment.

Considering the smoothing function over time, the model considered is described by:

in which Weightij means the weight of the i-th subject in the j-th week, for i = 1,...,57, j = 1,..., 12. The homeopathic treatment groups range from k = 1,..., 5. Moreover, f₁ refers to the mean curve for the individuals allocated in the Non-infected control group, f₂ refers to mean curve for the subjects that received the homeopathic treatment with 13cH dilution, f₃ refers to the mean curve for the individuals that received the treatment with 6cH dilution, f₄ refers to the mean curve for the individuals allocated in the Infection control group, and f5 refers to the mean curve for the individuals that received the homeopathic treatment with 3cH dilution, while a_i1 and a_i2 refer to the random effects corresponding to the intercepts and slopes of the curves.

Also, f(Time) is the unknown smoothing function that describes the trajectory of this variable in relation to Weight evolution, in which f( ) is estimated by means of P-splines, with third-polynomial-degree B-splines basis function, and a second-order discrete penalty. The variability between subjects is determined through the ai random effect and a random effect in the slope for the Time variable, in which (a_ia, a_i2)’~N(0, E) and e_ij ~N(0, sigma_e), assuming the dependence of the random terms.

Table 2 displayed the fitted models, considering different structures of variances and covariances, defined by: a model without considering a correlation structure in the errors, and three other models using the AR(1), CS and ARMA(1,1) correlation structures. For model selection, the Akaike criterion was used.

Table 2.
Comparing models with different correlation structure

Observing Table 2 allows concluding that the model to be selected is the one with the AR(1) correlation structure, since it provides the lowest value for the AIC statistics; in addition, this structure allows improving the fit of the curves both for the independent variables and the dependent variable. Figure 5 contains the chart for the smoothed curve of the Time covariable after the determination of the correlation structure; it is possible to see that the Weight evolution presented an increase behavior until the 7th week approximately and, after this period, it decreased until the 12th week; as a consequence, the effect of this covariable follows a non-linear pattern.

Figure 5.
Smooth curve with 5.89 effective degrees of freedom (e.d.f), considering the AR(1) correlation structure.

Figure (6) displays the correlogram with the AR(1) structure in the residuals; it is possible to see a significant portion of the autocorrelation values contained within the significance levels, that is, there is a decrease in the issue found as to the assumption of independence in the residuals of the model without considering a correlation structure.

Figure 6.
Correlogram considering the AR(1) correlation structure.

Table 3 contains the estimates for the λ parameter, with their respective effective degrees of freedom (e.d.f), F test, and p-value. It can be concluded that the smoothing of the Time covariable for each homeopathic treatment group was highly significant in relation to Weight evolution.

Table 3.
Summarized measures for each treatment.

Figure 7 shows the behavior of the estimated smoothing non-parametric functions of model 11 for each treatment group; the continuous line is the estimated smoothing function, and the gray line represents the ranges with 95% confidence. Because the effect of the Time variable for each treatment group is not linear, the estimated smoothing functions are significant with 3.751, 4.763, 2.726, 6.864 and 2.000 effective degrees of freedom. Additionally, the evolution of the weights over time has a period of increase until the 7th week approximately in the groups under homeopathic treatment with 13cH dilution (Figure 7b), 6cH dilution (Figure 7c), and 3cH dilution (Figure 7e) with the passing of time; as for the Non-infected control group (Figure 7a) and the Infection control group (Figure 7d), there is a slight rise in the beginning of the treatment, with an oscillation as of the 2nd week. Figures 7b and 7d show a significant decrease, thus evidencing that the animals lost weight after said period.

Figure 8 contains the specific estimated curves for the subjects allocated in their treatment groups; model 11 describe said curves, taking into account the non-linearity and the associated random effects. It is possible to observe the fitting capacity of the model using 10 nodes, which presents itself satisfactorily as to interpolation.

Thus, with the model defined in 11, it is possible to describe the mean trajectories by treatment group and verify the objective proposed in the beginning of this application, in which Figure 9 presents the estimated curves by treatment group:

Figure 7.
Smoother for Weekly Time in each Treatment Group.

Figure 8.
Individual curves for the subjects, fitted by treatment (G1: Group – Control; G2: Group – 13cH; G3: Group – 6cH; G4: Group – Infection Control, and G5: Group – 3cH).

Observing Figure 9, it is possible to verify that the animals allocated in the control group had, during a significant portion of the study, an average weight superior to that of the other treatment groups. When it comes to the group under homeopathic treatment with 13cH dilution, none of the animals reached the end of the study and, in said group, the greatest dilution of the biotherapeutic serum was used for treating the Chagas disease. Finally, one can further observe that, on average, the subjects allocated in the groups under homeopathic treatment with 13cH, 6cH and 3cH dilutions and the Infection control group had a decrease in weight as of the 7th week.

Figure 9.
Estimated curves by treatment.