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 |

# 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:

### 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\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e_t})\), where the variance is estimated for each level of`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.

```
<- asreml(
asr008 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`

.