ASR007. LMM with fixed effect interaction - Rat pups

The complete script for this example can be downloaded here:

Dataset

In this example we will use the D005 dataset, 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

We will fit a LMM for comparing the effect of the independent variables on the birth weights of rat pups. The model can be expressed as:

\[ 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\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e})\).


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.

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


Exploring output

The visual diagnostics of the residuals can be done with:

plot(asr007)

Note that one data point clearly deviates from the assumed Normal distribution and its expected performance. Outliers can be identified for removal using the residuals() function as shown below:

res <- residuals(asr007)
d005$weight[which.min(res)] <- NA


Below is the command in ASReml-R to re-fit (update) the model without the outlier:

asr007 <- update(asr007)


The statistical significance of fixed effects can be evaluated with:

wald(asr007, denDF = 'numeric')$Wald
              Df denDF     F.inc           Pr
(Intercept)    1  23.1 9.013e+03 0.000000e+00
lsize          1  28.5 3.157e+01 4.845266e-06
treatment      2  23.8 1.388e+01 1.013476e-04
sex            1 297.4 6.093e+01 1.018075e-13
treatment:sex  2 299.1 2.975e-02 9.706930e-01


The BLUEs for each of the fixed effects are:

summary(asr007, coef = TRUE)$coef.fixed
                                 solution  std error     z.ratio
(Intercept)                   8.057416225 0.27214843 29.60669782
lsize                        -0.133471396 0.01858875 -7.18022294
treatment_Control             0.000000000         NA          NA
treatment_High               -0.896290954 0.19257886 -4.65415016
treatment_Low                -0.464506379 0.15913967 -2.91885983
sex_Female                    0.000000000         NA          NA
sex_Male                      0.340039711 0.06543112  5.19691103
treatment_Control:sex_Female  0.000000000         NA          NA
treatment_Control:sex_Male    0.000000000         NA          NA
treatment_High:sex_Female     0.000000000         NA          NA
treatment_High:sex_Male      -0.028641244 0.11744535 -0.24386869
treatment_Low:sex_Female      0.000000000         NA          NA
treatment_Low:sex_Male       -0.009322738 0.09421320 -0.09895362


The random effects (BLUPs) are:

BLUP <- summary(asr007, coef = TRUE)$coef.random
head(BLUP)
           solution std.error    z.ratio
litter_1  0.1470404 0.1394820  1.0541894
litter_2 -0.1017843 0.1344999 -0.7567609
litter_3 -0.2571826 0.2196628 -1.1708067
litter_4 -0.0860600 0.1344756 -0.6399673
litter_5  0.3224513 0.1357976  2.3744988
litter_6  0.2288453 0.1648274  1.3883935


More information

Note that an extension of this model, where the residual variance is estimated under each level of treatment, is available here: ASR008