Introduction

There has been a significant increase worldwide in the capacity to produce, store and transmit information; the latter is characterized as data and demands more and more advances from statistics, both as to the development of methodologies and as to new, ever-complex indicators that require modern equipment, statistical software and trained technicians (Ignácio, 2010). Mixed Models were widely studied by Fisher in 1918, with major impacts on quantitative genetics studies, and referred to by the author as components of variance models (Scheffe, 1999). The development of linear mixed models combined in one single equation is a result of primordial investigations (Harville, 1976; 1977) that facilitated this achievement; later, their use would be discussed in (Laird & Ware, 1982) for longitudinal data, which are data characterized by a time sequence of two or more observations on each individual.

The applications of Hierarchical Linear Models have been growing due to the great extension of problems that have hierarchically structured data, just as in this study, which analyzed body weight data of individuals infected by Trypanosoma cruzi, the agent of Chagas disease. The latter is one of the most widely distributed pathologies in the American continent. Vectors of the disease can be found from the south of the United States to Argentina. There are more than one hundred species responsible for the natural transmission of the infection by Trypanosoma cruzi, directly helping it spread in the home environment or participating in the maintenance of chagasic enzooty. It is estimated that 16 to 18 million individuals are infected, and that approximately 80 million people are at risk of contamination in Latin America (Schmunis, 1997; WHO, 1991 apud Vinhaes & Dias, 2000). Hence the importance of the several studies conducted in this field.

Material and methods

The experiment described in this research was carried out by (Ferreira et al., 2018) through blind, controlled and randomized assays. The objective was to verify the efficiency of using variations of a biotherapeutics produced from a serum of Gallus gallus domesticus (chicken), under parasitological, clinical and immunological parameters, in mice experimentally infected with Trypanosoma cruzi. The expected efficiency of each biotherapeutics would be observed in the subjects' body weight.

Also, according to the authors, the experiment used 57 Swiss mice, aged 56 days old, all male and sourced from the Central Vivarium of the State University of Maringá. The animals were distributed into treatment groups and housed in cages with a maximum of 5 animals. With polysulfone (ALESCO^®), 20 × 32 × 21 cm³ in dimension, controlled temperature (22 ± 2ºC), a 12-hour light/dark cycle, and water and feed being supplied ad libitum, the cages were transformed into a microenvironment. All groups were subjected to the same experimental conditions.

The animals were divided into 5 groups and subdivided into 14 cages: group 1 was composed of 5 animals allocated in cage 1; groups 2, 3 and 4 were composed of 13 animals and divided into 3 cages each, respectively in cages 2, 3, 4, 5, 6, 7, 8, 9 and 10; and group 5 was composed of 13 animals as well, but they were divided into 4 cages – 11, 12, 13 and 14, respectively. The infected animals were inoculated intraperitoneally with 1,400 blood trypomastigotes (infectious forms) of Trypanosoma cruzi (Y strain) (Nussenzweig et al., 1953).

Here is a description of the experimental groups (treatments):

1: NIC - (Non-infected control) Non-infected and non-treated animals (n = 5);

2: G13cH - Animals treated with 13cH chicken serum biotherapeutics (n = 13);

3: G6cH - Animals treated with 6cH chicken serum biotherapeutics (n = 13);

4: ICG - (Infection control) The animals were infected and received no treatment (n = 13);

5: G3cH - Animals treated with 3cH chicken serum biotherapeutics (n = 13).

The medicine was diluted in water (1mL 10mL^-1) and supplied ad libitum in a sterile amber bottle, in accordance with Aleixo et al. (2013), for 16 consecutive hours (medicine available for the animals from 16:00 to 8:00), on the 4^th, 7^th and 10^th days after infection (totaling 3 doses). The treatment scheme is based on the action of the drug, which is linked to its immunological effects and to the specific evolution of the Y strain of Trypanosoma cruzi in Swiss mice (Aleixo et al., 2013; Ferraz et al., 2016). The project from which the experiments were run was approved by the Ethics Committee on Research Involving Animals of the State University of Maringá, Paraná, CEUA Opinion 2401220716/2016.

Body weight evolution was monitored over 12 weeks using a semi-analytical balance (Balance BEL^®). In this study, the assessment performed at the beginning of the treatment was used to analyze the initial weight of the mice. Due to a great loss of information as of the 7^th study week, the modeling was based on the longitudinal measures observed up to the 6^th week.

Hierarchical Linear Models

For modeling the weight of the mice (y) as a function of the predictor variables (Z) at the individual level and/or at a higher level (W), the dataset is assumed to be multilevel; theoretically, the model can be presented as a hierarchical system of regression equations.

The Hierarchical Linear Model, in its general form, can be formulated through two equations:

in which (1) represents the model equation within the group – level 1 –, and (2) represents the model between groups – level 2. The explanatory variables at the subject level are represented by index (q),while the explanatory variables at the group level are represented by index (p), so (q) and (p) represent the number of variables. Thus:

y_i is a result vector.

is a matrix of explanatory variables;

i s one of the components of the vector of unknown fixed parameters;

is the vector of level-1;

is a matrix of level-2 explanatory variables;

y is a vector of fixed effects;

Is a vector of random effects.

This model can be combined to result in a single Linear Mixed Model, with and

Generally referred to as Hierarchical Model in (Verbeke & Molenberghs, 2000). corresponds to the fixed part, with representing the vector of components , while corresponds to the random part of the Model.

From the model given in (3), errors will be assumed as independent between groups, and different between levels, that is, it is assumed that random effects and errors follow a normal distribution, with mean equal to 0, and have correlated residuals, where Phi is the covariance between b_iand e_i, and (D) and (R) are the covariance matrices, with both matrices being positive and defined by hypothesis, therefore not singular:

where

and

with the following assumptions, , ; can be zero for .

Frequently, it is assumed that . Thus, the total variance of the model (3) for response vector is given as:

where such considerations mean that:

Parameter Estimation

When there is a Hierarchical Linear Model, in the form given in (3), with matrix of variances and covariances such as those presented in (7), there is usually interest in predicting fixed and random effects and in estimating the components of variance. Commonly, as a procedure for estimating Hierarchical Linear Models, Maximum Likelihood Estimation (MLE) is used. The Maximum Likelihood Estimation method consists of maximizing the likelihood function of observations in relation to the fixed effects and to the components of variance, requiring data normality assumption. Assuming that the vector of the observations has as mean, and as matrix of variances and covariances, the likelihood function (L) of is:

where represents the determinant of the matrix.

Model diagnoses

The Hierarchical Linear Model, just as in Ordinary Linear Regression, has distributive assumptions that may or may not be valid when used in practice. However, the diagnoses to assess these assumptions and consequent alternatives, for assumption violation suspicions, have not been fully developed for this model, mainly because the analysis tool is relatively recent.

The diagnoses were divided into two stages, residual analysis and influential point analysis; for the residual analysis, the model given in (3) encompasses the uncertainty both at the individual level ei and and at the group level b_i. A residual term should allow assessing the distributional assumptions of the model. Thus, level-1 (individual) residuals and level-2 (group) residuals are deemed essential; therefore, regardless of the form of the covariance matrix, these residuals are of interest. All upper-level residuals are the best linear unbiased predictors (BLUPs) of random effects. Following the convention established in (Goldstein, 2011) and (Raudenbush & Bryk, 2002) for level specification, the residuals are titled by the level at which they were introduced in the Model.

Checking if their definitions are interrelated, level-1 and level-2 residuals are fundamental for modeling, since a deficiency at one level of the model can be perceived in the residual analysis at another level. An ascending residual analysis is recommended for verifying the validation of level-1 residuals; this way, having concluded the appropriate model for this level, move on to level 2. By doing so, the impacts on the residual analysis, resulting from a possible confusion when modeling, will be minimized (Loy & Hofmann, 2013).

One way of indicating whether the residual variation within the group is constant between groups is to use the programming protocol presented by (Buja et al., 2009), which employs several sets of simulated data from which null graphs are built.

In model fitting and parameter estimation, not all observations, or groups, have the same effect. Some observations or groups stand out from the others, and the fit of the model detects these differences. These observations or groups are referred to as influential or leverage points. We are especially interested in the points of influence on fitted values, fixed effects (estimates) and components of variance.

Considering that, in addition to fixed effects, there are random effects influencing the result of a studied phenomenon, Hierarchical Linear Models can be used to study the best covariance structure. For this case, we assume that the covariance structure, , is fixed, which is a generalization of the linear regression, with denoting leveraging at level i.

The multiple statistics that can define ‘leverage points’ for fitted values in a Hierarchical Linear Model are described by general leverage points (H), leverage points in fixed effects (H₁), leverage points in random effects (H₂), and leverage points in non-confounding random effects (H₂^*). Following the definition provided by (Demidenko & Stukel, 2005), the leveraging of group i is the sum of the leverages for fixed effects H₁i and random effects H_2i , where

Here is some confusion as to the diagnosis of influential points, which occurs between levels. Because the leveraging of the random effects (11) results from the leveraging of the fixed effects (10). Optionally, according to (Nobre & Singer, 2011), the leveraging for the random effects can be defined as

which solves such a confusion.

Influential observations in fixed-effect estimations can be made using Cook's distance, or also MDFFITS statistics, which is a multivariate version of DFFITS statistics (Belsley, Kuh, & Welsch, 2004). Both statistics determine the distance between fixed-effect estimations deriving from complete data, and those from reduced data, and are generalized by (Christensen, Pearson, & Johnson, 1992) and (Schabenberger, 2005) for the Hierarchical Linear Model as follows

These statistics present great values for influential observation and, because is used, there is no exact reference distribution for it.

Components of variance estimation is another attribute that allows analyzing the influence of possible model observations. Even though components of variance estimations are not of primary interest for researchers, it is fundamental to investigate this part of the model because components of variance impact its fixed part (Loy & Hofmann, 2013). For the Variance Components, one can directly compare the relative change for each variance component, theta₁

Note that the relative variance change (RVC) will be close to zero when the i-th unit does not influence the component of variance in question.

All decision measures for choosing the best proposed fit, estimating parameters and diagnosing the model were obtained from software R version 3.4.4 (R Core Team, 2019) with the following packages (ggplot2, HLMdiag, lme4, nlme, mlmRev, nullabor, fitdistrplus, dbplyr , plyr).

The model selection was based on the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC) and the likelihood-ratio test.

Results

Initially, a descriptive analysis was performed to ascertain the behavior of the subjects' weight during the experiment; this analysis was divided in two parts – subjects' initial weight (1^st week), and subjects' weight in the following experiment weeks (2^nd to 6^th week). The dataset used for both the descriptive analysis and the fittings is made up of 342 observations, of which 16 were lost (not recorded).

Mice weight descriptive analysis

To study the initial weight of the mice, a statistical summary was produced, observing that the means of the groups range from 39.24 g for group 4 to 42.19 g for group 2, with standard deviations equal to 3.68 and 4.04 g, respectively.

Figure 1 shows the scatter of the weight of each subject at the beginning of the experiment and how discrepant the weights of the mice are. The treatment group to which each subject belongs is identified with different colors.

Figure 2 shows the weight evolution of each subject as a function of time, evidencing a similar behavior as the weeks go by. In addition, Figure 3 shows that all groups had an average weight increase as a function of time.

The highest mean between groups was found in the 6^th week, for group 2, while the lowest one was found in the 2^{nd week, for group 3. The largest deviation between groups also occurred in the 6th week of group 2, while the lowest} one occurred in the 5^th week of group 3.

Figure 4 displays a histogram, a Cumulative Distribution Function (CDF) graph, a QQ-Plot and a PP-Plot, all as a function of weight observations for all treatments from the 2^nd to the 6^th week, considering normal distribution as reference.

The Kolmogorov-Smirnov test presented D = 0.077755 and p-value = 0.077 as statistics, indicating that the subjects' weights follow normal distribution.

Model fit

The analysis of the models sought to ascertain the behavior of the subjects infected by Trypanosoma cruzi by relating the response variable 'weight' to possible influencing factors, cage and time, which belong to the fixed part of the model. The treatment group and the subjects themselves belong to the random part.

The models were fitted for response variable , in which i = 1, 2,. . . , 5 (group number), j = 1, 2, . . . , ni (number of subjects, where if i = 1, ni = 5, and if i = 1, ni = 13), k = 1, . . . , mi (cage number, where if i = 1 mi = 1; i = 2 mi = 2, 3, 4; i = 3 mi = 5, 6, 7; i = 4 mi = 8, 9, 10; i = 5 mi = 11, 12, 13, 14) and t = 2, . . . , 6 (treatment week).

The proposed models are described in Table 1; the first one is given only by the intercept of the fixed part, while the effect of the subjects is given by means of a random intercept for each one of them. In the second, the effect of time is added to the fixed part; in the third model, the effect of the groups is added to the fixed part; in the fourth model, in addition to time and group, the effect of the cages is added to the fixed part.

All proposed models presented violation in the residual normality assumption. The graphs in Figure 5 show how observation 16 stands out in 3 of the 4 fittings, so it is an observation worth investigating.

Observation 16 corresponds to the body weight of subject 4 in the second week; said subject belongs to group G1, which is the control with non-infected animals. Table 2 specifically shows this subject's data; it is possible to see that, in the second week, its weight had a sharp rise in relation to its initial weight, which was 41.00 g, then a drop in the 3rd week, presenting small variations until the 6th week.

This way, there is evidence that this observation has gone through a collection error, so the fittings were redone without observation 16.

Through the Shapiro-Wilk residual normality test, Model 2, Model 3 and Model 4 obtained, as p-value, 0.267, 0.266 and 0.265, respectively, thus presenting residual normality.

The values of all ANOVA decision measures for each selected model are displayed in Table 3.

Models 2, 3 and 4 are statistically equivalent but, considering the model selection criteria, AIC and BIC, and the likelihood ratio test, Model 2 was chosen as the one the best explains the data. Another factor that led it to be chosen is that Model 2 is more parsimonious than the other equivalent ones.

Table 4 presents estimations of fixed-effect parameters for Model 2 with their respective confidence intervals

The validity of Model 2 can be confirmed by verifying its assumptions, which were assessed through diagnosis analysis.

Residual analysis

The level-1 and level-2 residuals of Model 2 were analyzed in an ascending manner. This analysis mode is necessary because the residuals are inter-related, therefore confounding, and can make the model diagnosis difficult if not analyzed correctly.

Residuals resulting from level 1 of least squares (LS) are not confounded with the residuals at level 2. By fitting the LS regression models separately, random effects are treated as fixed. Figure 6 presents a graph for the level-1 residuals of the LS, Time. It is suggested that time may not be linearly related to weight. In order to address the assumption of level-1 homoscedastic residuals, we used semi-standardized residuals, just as shown in (Snijders & Berkhof, 2008). Figure 7 display this graph, which indicates no linearity violation by said assumption.

Figure 8 shows the normal quantile graph of the level-1 semi-standardized residuals; visually, the semi-standardized residuals seem normal, thus showing no evidence against their normality assumption. The normality assumption of the residuals was assessed by the normal quantile graph presented in Figure 9.

The programming protocol proposed by (Buja et al., 2009) was also used to corroborate with the verification of the assumption of level-1 homoscedastic residuals. It is displayed in Figure 10. It is not possible to distinguish the real-data graph from the simulated graph. Therefore, the analysis can move on without need for corrective measures.

The random effects, commonly referred to as level-2 effects, are defined by Zb, or b only. As previously said, an ascending residual analysis is performed, so residuals EB at level 2 are used. This choice was made because least square residuals are more variable than EB residuals, and EB estimates of b were obtained directly.

Analysis of influential observations

This sub-section will present the use of diagnoses to assess changes in the components of variance estimation using RVC, the fixed-effect estimation using Cook's distance, and fitted values using leveraging. These quantities are used to assess influences in level-1 and level-2 units.

For fixed effects, we have two statistics commonly used to measure whether there has been any changes, namely Cook's distance and MDFFITS, already mentioned. Figure 11 shows that no subject was identified by both statistics (which present the same values) as influential in the fixed effects, considering an internal scale. Figure 11 displays graphs for the abovementioned statistics in the form of dot plot and modified dot plot.

For the components of variance, we used the RVC, which was presented in equation 15; it measures changes in the estimations of the A-th component of variance, θA, with and without unit i.

Table 5 illustrates the RVC, which presents as output a matrix where each column represents a component of variance, σ2, which is the residual variance, and D11, which is the variance associated with the random intercept of the subjects. Note that the value farthest from zero in D11 for n=95, characterizing this unit as influential in the component of variance.

Figure 12 shows the modified dot plot of the level-2 RVC for the random intercept of the subjects. By using the internal scale, n=95 was identified as influential unit, just as shown in Table 8, and is worth attention.

Subject n=95 corresponds to the weight of subject 38 at the beginning of the treatment; as shown in Table 6, said subject is the heaviest among the animals. This justifies the fact that it was identified as influential in the graph of Figure 12.

In the fitted values, in addition to checking how the observations of one same subject directly impact the parameters of the adopted model, it was interesting to explore whether these observations are atypical in relation to the fitted values and to the explanatory variables of such model. This exploration was done through multiple statistics.

These statistics are displayed in Table 7, which reveals that subjects 9, 27, 44 and 56 have high leveraging in the fixed effects of the fitted values. Likewise, these subjects are identified as having high leveraging in the random effects, which is natural, since H2 depends on H1. Considering H2, all subjects present the same metrics, so all of them have the same influence under the random effects of the fitted values of Model 2.

Because the intercept of subject # 38 was identified as being influential in the components of variance, and the observations for subjects # 9, 27, 44 and 56 proved to be influential in the fitted values, both were removed, and the fitting of the chosen model was repeated in order to indicate whether there was any significant change in the estimations of the parameters.

The estimations of the parameters without influential observations are displayed in Table 8.

Note that, in comparison with Table 4, the parameters went through small changes. Therefore, since the influential observations have little impact on the selected model, we chose to keep them in order to preserve the data. Thus, we consider the initial estimations for Model 2.

Conclusion

The proposed Linear Mixed Model proved to be satisfactory for the dataset with hierarchical structure, besides managing to describe the behavior of the subjects' weight during the experiment weeks. The Group and Cage variables showed no significance for the Model; only the Time variable and the intercept (initial weight) were capable of describing the weight of the individuals in the Model. Such an importance of the intercept in the Model is due to it being considered as random, since each subject presented a different initial weight. Thus, the Model with only intercept and time was capable of describing the subjects' weight, with acceptable estimations.

Taking the analysis results into account, we suggest, as recommendation, that researchers should start their experiments with the subjects' initial weight as similar as possible and divide the subjects equally into groups. It is also advisable to test new biotherapeutics for the treatment, as the ones analyzed in the present study (G13cH, G6cH and G3cH) did not prove to be efficient.

References

Aleixo, D. L., Veiga, F. K., Ciupa, L., Portocarrero, A. R., Sandri, P. F., Ferraz, F. N., ... Araújo, S. M. (2013). Highly diluted medication modifies parasitological, immunological and histopathological parameters of Trypanosoma cruzi infection. International Journal of High Dilution Research, 12(44), 107-108.

Belsley, D. A., Kuh, E., & Welsch, R. E. (2004). Regression diagnostics: Identifying influential data and sources of collinearity (Wiley series in probability and statistics). Hoboken, NJ: John Wiley & Sons, Inc.

Buja, A., Cook, D., Hofmann, H., Lawrence, M., Lee, E. K., Swayne, D. F., & Wickham, H. (2009). Statistical inference for exploratory data analysis and model diagnostics. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1906), 4361-4383. doi: 10.1098/rsta.2009.0120

Christensen, R., Pearson, L. M., & Johnson, W. (1992). Case-deletion diagnostics for mixed models. Technometrics, 34(1), 38-45. doi: 10.2307/1269550

Demidenko, E., & Stukel, T. A. (2005). Influence analysis for linear mixed‐effects models. Statistics in Medicine, 24(6), 893-909. doi:10.1002/sim.1974

Ferraz, F. N., Aleixo, D. L., Gruendling, A. P., Gomes, M. L., Toledo, M. J. O., & Araújo, S. M. (2016). Trypanosoma cruzi: Evaluation of PCR as a laboratory tool to follow up the evolution of parasite load. Iranian Journal of Parasitology, 11(3), 389-395. PMID: 28127346

Ferreira, É. C., Ciupa, L., Portocarrero, A. R., Brustolin, C. F., Massini, P. F., Aleixo, D. L., & Araújo, S. M. (2018). Phosphorus protects cardiac tissue by modifying the immune response in rats infected by Trypanosoma cruzi. Cytokine, 102, 102-106. doi: 10.1016/j.cyto.2017.07.017

Goldstein, H. (2011). Multilevel statistical models (4th ed.). Bristol, UK: John Wiley & Sons

Harville, D. A. (1976). Extension of the Gauss-Markov theorem to include the estimation of random effects. The Annals of Statistics, 4(2), 384-395. Retrieved from www.jstor.org/stable/2958209

Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358), 320-338. doi: 10.2307/2286796

Ignácio, S. A. (2010). Importância da estatística para o processo de conhecimento e tomada de decisão. Revista Paranaense de Desenvolvimento, (118), 175-192.

Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963-974. PMID: 7168798

Loy, A., & Hofmann, H. (2013). Diagnostic tools for hierarchical linear models. WIREs Computational Statistics, 5(1), 48-61. Doi 10.1002/wics.1238

Nobre, J. S., & Singer, J. M. (2011). Leverage analysis for linear mixed models. Journal of Applied Statistics, 38(5), 1063-1072. doi: 10.1080/02664761003759016

Nussenzweig, V., Sonntag, R., Biancalana, A., Freitas, J. D., Amato Neto, V., & Kloetzel, J. (1953). Ação de corantes trifenil-metânicos sobre o Trypanosoma cruzi “in vitro”. Emprego da violeta de genciana na profilaxia da transmissão da moléstia de Chagas por transfusão de sangue. O hospital, 44(6), 731-744.

R Core Team. (2019). R: A Language and Environment for Statistical Computing - Version 3.4.4. Vienna, AU: R Foundation for Statistical Computing. Retrieved from https://www.R-project.org/

Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (Advanced quantitative Techniques in the Social Sciences, Vol 1, 2nd ed.). Thousand Oaks, CA: Sage Publications

Schabenberger, O. (2005). Mixed model influence diagnostics. In SAS Users Group International Conference (Vol. 29, paper 189). Retrieved from https://support.sas.com/resources/papers/proceedings/proceedings/ sugi29/189-29.pdf

Scheffe, H. (1999). The analysis of variance. New York: Wiley-Interscience Publication.

Snijders, T. A., & Berkhof, J. (2008). Diagnostic checks for multilevel models. In J. Leeuw, & E. Meijer (Eds.), Handbook of multilevel analysis (p. 141-175). New York, NY: Springer. Doi: 10.1007/978-0-387-73186-5_3

Verbeke, G., & Molenberghs, G. (2009). Linear mixed models for longitudinal data. New York, NY: Springer.

Vinhaes, M. C., & Dias, J. C. P. (2000). Doença de chagas no Brasil. Cadernos de Saúde Pública, 16(Suppl. 2), S7-S12. Doi: 10.1590/S0102-311X2000000800002

Author notes

roney_peterson@hotmail.com