ASR003. Simple LMM - Twins’ IQ

The complete script for this example can be downloaded here:

Dataset

In this example we will use the D002 dataset, and the first few rows are presented below:

pair	twin	iq
112	A	113
112	B	109
114	A	94
114	B	100
126	A	99
126	B	86

Model

The basic LMM that we will fit in this example is:

\[ y = \mu + twin + pair + e\ \] where,

\(y\) is the IQ score,

\(\mu\) is the overall mean,

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

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

\(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, twin and pair need to be set as factors.

asr003 <- asreml(
  fixed = iq ~ twin,
  random = ~pair,
  residual = ~units,
  data = d002
)

Exploring output

Evaluate the significance of the fixed effects:

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

            Df denDF    F.inc         Pr
(Intercept)  1    31 1499.000 0.00000000
twin         1    31    3.416 0.07412364

Having the factor twin as not significant is relevant, as we are not expecting differences between them in terms of IQ.

Lets obtain the predicted means, based on the above model with:

predict(asr003, classify = 'twin')$pvals

  twin predicted.value std.error    status
1    A        93.21875  2.568317 Estimable
2    B        96.12500  2.568317 Estimable

We can also check the variance components:

summary(asr003)$varcomp

        component std.error  z.ratio bound %ch
pair    171.53572  48.84192 3.512059     P   0
units!R  39.56357  10.04961 3.936825     P   0

We can estimate the intra-correlation between twins (\(r^2\)), that, if close to one implies that we have a strong genetic component in this response variable. This can be done using:

vpredict(asr003, r2 ~ V1/(V1+V2))

    Estimate         SE
r2 0.8125831 0.06100241

Finally, we can report the random effects (BLUP) associated with each twin pair.

BLUP <- summary(asr003, coef = TRUE)$coef.random
head(BLUP)

           solution std.error    z.ratio
pair_112  14.639529  4.747755  3.0834633
pair_114   2.087359  4.747755  0.4396517
pair_126  -1.947267  4.747755 -0.4101449
pair_132 -14.499437  4.747755 -3.0539564
pair_136  -5.981893  4.747755 -1.2599414
pair_148   3.432234  4.747755  0.7229172

More information

Note that an alternative of this model, using a more complex residual structure, can be found here: ASR025.