ASR008. LMM with fixed effect interaction and heterogeneous residual variances - Rat pups

The complete script for this example can be downloaded here:

Dataset

The dataset we will use is D005, and the first few rows are presented below:

weight	sex	litter	lsize	treatment
6.60	Male	1	12	Control
7.40	Male	1	12	Control
7.15	Male	1	12	Control
7.24	Male	1	12	Control
7.10	Male	1	12	Control
6.04	Male	1	12	Control

Model

The LMM that we will fit here is an extension of the example ASR007 which compares the effect of independent variables on the birth weights of rat pups. Here, a residual variance will be estimated for each level of the factor \(treatment\) as detailed below. The model can be expressed as follows:

\[ y = \mu + lsize + treatment + sex + treatment:sex + litter + e\ \] where,

\(y\) is the birth weight of rat pups,

\(\mu\) is the overall mean,

\(lsize\) is the fixed effect of litter size,

\(treatment\) is the fixed effect of treatment,

\(sex\) is the fixed effect of pup sex,

\(treatment:sex\) is the fixed effect of the interaction between treatment and sex,

\(litter\) is the random effect of litter, with \(litter \sim \mathcal{N}(0,\,\sigma^{2}_{l})\),

\(e\)

\(e \sim \mathcal{N}(0,\,\sigma^{2}_{e_t})\)

treatment

\(t=\{1,2,3\}\)

Now, let’s take a look at how to write the model with ASReml-R. Note that before fitting the model sex, litter, and treatment need to be set as factors.

asr008 <- asreml(
  fixed = weight ~ lsize + treatment + sex + treatment:sex,
  random = ~litter,
  residual = ~dsum(~units|treatment),   
  data = d005
)

The dsum() function is used for structuring the residual term. This corresponds to the direct sum or a block diagonal, where a residual variance component will be estimated for each level of the conditioning factor (e.g., treatment).

Exploring output

After fitting the mode, the statistical significance of the fixed effect terms can be evaluated with:

wald(asr008, denDF = 'numeric')$Wald

              Df denDF     F.inc           Pr
(Intercept)    1  23.1 9064.0000 0.000000e+00
lsize          1  28.5   29.6300 7.817949e-06
treatment      2  23.9   11.8800 2.611454e-04
sex            1 249.8   67.1400 1.321165e-14
treatment:sex  2 155.6    0.3174 7.285403e-01

The BLUEs for each of the fixed effects are:

summary(asr008, coef = TRUE)$coef.fixed

                                solution  std error    z.ratio
(Intercept)                   7.93716099 0.27736597 28.6162031
lsize                        -0.13000719 0.01848733 -7.0322302
treatment_Control             0.00000000         NA         NA
treatment_High               -0.80860961 0.19628409 -4.1195881
treatment_Low                -0.39027870 0.16301767 -2.3940884
sex_Female                    0.00000000         NA         NA
sex_Male                      0.40813188 0.09303422  4.3869006
treatment_Control:sex_Female  0.00000000         NA         NA
treatment_Control:sex_Male    0.00000000         NA         NA
treatment_High:sex_Female     0.00000000         NA         NA
treatment_High:sex_Male      -0.09466631 0.12919552 -0.7327368
treatment_Low:sex_Female      0.00000000         NA         NA
treatment_Low:sex_Male       -0.07601342 0.10811789 -0.7030605

The variance components estimated from this model are:

summary(asr008)$varcomp

                     component  std.error  z.ratio bound %ch
litter              0.09827257 0.03302492 2.975709     P   0
treatment_Control!R 0.26501379 0.03414137 7.762248     P   0
treatment_Low!R     0.08459389 0.01114410 7.590910     P   0
treatment_High!R    0.10835920 0.02031969 5.332719     P   0

Note that a different residual variance is estimated for each treatment level, and these appear very different.

The random effects (BLUPs) are:

summary(asr008, coef = TRUE)$coef.random
head(BLUP)

            solution std.error    z.ratio
litter_1  0.15990442 0.1631703  0.9799850
litter_2 -0.06959703 0.1564612 -0.4448196
litter_3 -0.15900548 0.2341198 -0.6791630
litter_4 -0.05114224 0.1564071 -0.3269816
litter_5  0.32064736 0.1588866  2.0180898
litter_6 -0.05665401 0.1840404 -0.3078346

And these are only for the factor litter.