ref <-
vlm_reference(
source = "Example 3"
, evaluate = TRUE
, level = 1
)\[ \newcommand{\vy}{\mathbf{y}} \newcommand{\vx}{\mathbf{x}} \newcommand{\vb}{\mathbf{b}} \newcommand{\ve}{\mathbf{e}} \newcommand{\vone}{\mathbf{1}} \newcommand{\vbeta}{\boldsymbol{\beta}} \newcommand{\vmu}{\boldsymbol{\mu}} \newcommand{\mX}{\mathbf{X}} \newcommand{\mZ}{\mathbf{Z}} \newcommand{\mV}{\mathbf{V}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mJ}{\mathbf{J}} \newcommand{\mSigma}{\boldsymbol{\Sigma}} \newcommand{\vg}{\mathbf{g}} \newcommand{\vc}{\mathbf{c}} \newcommand{\vz}{\mathbf{z}} \newcommand{\vlambda}{\boldsymbol{\lambda}} \newcommand{\mA}{\mathbf{A}} \newcommand{\mB}{\mathbf{B}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mD}{\mathbf{D}} \newcommand{\mL}{\mathbf{L}} \newcommand{\mQ}{\mathbf{Q}} \newcommand{\vu}{\mathbf{u}} \newcommand{\vr}{\mathbf{r}} \newcommand{\vl}{\mathbf{l}} \newcommand{\vk}{\mathbf{k}} \newcommand{\vtheta}{\boldsymbol{\theta}} \newcommand{\mG}{\mathbf{G}} \newcommand{\mR}{\mathbf{R}} \newcommand{\mP}{\mathbf{P}} \newcommand{\mC}{\mathbf{C}} \newcommand{\mH}{\mathbf{H}} \newcommand{\mW}{\mathbf{W}} \newcommand{\mPhi}{\boldsymbol{\Phi}} \newcommand{\E}{\mathrm{E}} \newcommand{\Var}{\mathrm{Var}} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\tr}{\operatorname{tr}} \]
1 Abstract
Chapter 3 of the book applies the mixed model to a designed experiment with a covariate (Example 3.1), to the estimation of heritability and breeding values (Example 3.2) and to repeated measures with several covariance structures (Example 3.3). The book prints only part of the data of Examples 3.1 and 3.2, so their analyses run but cannot reproduce the printed numbers; the complete data of Example 3.3 reproduce every printed number.
2 Example 3.1: doses nested in drug, with a covariate
m1 <-
vlm_fit(
fixed = PCV2 ~ drug + drug:dose
, data = ex31
, random = ~ drug:herd
, repeated = NULL
, method = c("REML", "ML", "Type3")
, ddfm = "satterth"
, type = 3
, class = c("herd", "dose")
, bound = TRUE
, control = list()
)
m1$tests effect num_df den_df F p
<char> <int> <num> <num> <num>
1: drug 1 10.08194 0.001544608 0.9694174
2: drug:dose 1 31.51959 2.600111652 0.1168256
vlm_estimate(
fit = m1
, L = list(
"Mean Samorin" = list(Intercept = 1, drug = c(0, 1), "drug:dose" = c(0, 0, 1))
, "Berenil 2 doses" = list("drug:dose" = c(1, -1, 0))
, "Ber vs Sam at dose 1" = list(drug = c(1, -1), "drug:dose" = c(1, 0, -1))
)
, random = NULL
, level = 0.95
)Estimates (denominator df: satterth)
Estimates
label estimate se df t p lower upper space
Mean Samorin 24.17 1.038 15.36 23.28 <1e-04 21.96 26.38 broad
Berenil 2 doses -2.102 1.303 31.52 -1.612 0.117 -4.758 0.5548 broad
Ber vs Sam at dose 1 -1.000 1.468 15.36 -0.6811 0.506 -4.123 2.123 broad
Sources: Duchateau, Janssen and Rowlands (1998); McLean, Sanders and Stroup
(1991); Satterthwaite (1946); Giesbrecht and Burns (1985); Fai and Cornelius
(1996). Full references: vlm_cite(x).Model 3 lets each treatment have its own slope on the PCV at treatment, and compares the treatments for an animal with PCV1 = 18:
m3 <-
vlm_fit(
fixed = PCV2 ~ drug + drug:dose + PCV1:drug:dose
, data = ex31
, random = ~ drug:herd
, repeated = NULL
, method = c("REML", "ML", "Type3")
, ddfm = "kr"
, type = 3
, class = c("herd", "dose")
, bound = TRUE
, control = list()
)
vlm_lsmeans(
fit = m3
, specs = ~ drug:dose
, at = list(PCV1 = 18)
, diff = FALSE
, adjust = c("none", "tukey", "bonferroni", "sidak")
, level = 0.95
)Least-squares means (denominator df: kr; 95% limits)
Least squares means
drug dose estimate se df t p lower upper
BERENIL 1 24.20 1.987 29.16 12.18 <1e-04 20.14 28.27
BERENIL 2 24.79 1.174 13.62 21.11 <1e-04 22.27 27.32
SAMORIN 1 25.26 2.169 26.49 11.64 <1e-04 20.80 29.71
Sources: Duchateau, Janssen and Rowlands (1998); Searle, Speed and Milliken
(1980); Kenward and Roger (1997). Full references: vlm_cite(x).3 Example 3.2: heritability and sire predictions
f32 <-
vlm_fit(
fixed = Ww ~ sex + agew + breed
, data = ex32
, random = ~ breed:sire_id
, repeated = NULL
, method = c("REML", "ML", "Type3")
, ddfm = "satterth"
, type = 3
, class = "breed"
, bound = TRUE
, control = list()
)
f32$covariance parameter subject estimate se z p
<char> <char> <num> <num> <num> <num>
1: breed:sire_id 0.000000 NA NA NA
2: Residual 5.447347 1.020382 5.338539 4.684922e-08
vlm_lsmeans(
fit = f32
, specs = ~ breed
, at = NULL
, diff = TRUE
, adjust = "tukey"
, level = 0.95
)$lsmeans breed estimate se df t p lower upper
<char> <num> <num> <num> <num> <num> <num> <num>
1: 1 13.123981 0.6055657 57 21.67227 3.343968e-29 11.911357 14.33661
2: 2 12.396055 0.7201282 57 17.21368 3.012048e-24 10.954023 13.83809
3: 3 10.598945 0.8773799 57 12.08022 2.364957e-17 8.842022 12.35587
4: 4 12.740250 0.6569679 57 19.39250 8.957490e-27 11.424695 14.05581
5: 5 9.722305 0.8625235 57 11.27193 3.896394e-16 7.995131 11.44948
6: 6 10.050959 0.9686565 57 10.37618 9.546083e-15 8.111258 11.99066On this part of the data the sire variance is estimated as zero, so the heritability \(h^2 = 4\sigma_s^2/ (\sigma_s^2 + \sigma^2)\) is zero; with the complete data the book finds 0.49.
4 Example 3.3: covariance structures for repeated measures
Model 1 (a random animal effect, equivalently compound symmetry), Model 2 (first-order autoregressive errors) and Model 3 (a random intercept and slope per animal):
cs <-
vlm_fit(
fixed = PCV ~ breed + breed:time
, data = ex33
, random = NULL
, repeated = vlm_repeated(
type = c("vc", "cs", "ar1", "sp_pow", "un")[2]
, subject = ~ animal_id:breed
, order = NULL
)
, method = c("REML", "ML", "Type3")
, ddfm = c("default", "contain", "residual", "betwithin", "satterth", "kr")
, type = 3
, class = NULL
, bound = TRUE
, control = list()
)
ar <-
vlm_fit(
fixed = PCV ~ breed + breed:time
, data = ex33
, random = NULL
, repeated = vlm_repeated(
type = c("vc", "cs", "ar1", "sp_pow", "un")[3]
, subject = ~ animal_id:breed
, order = NULL
)
, method = c("REML", "ML", "Type3")
, ddfm = c("default", "contain", "residual", "betwithin", "satterth", "kr")
, type = 3
, class = NULL
, bound = TRUE
, control = list()
)
rc <-
vlm_fit(
fixed = PCV ~ breed + breed:time
, data = ex33
, random = vlm_random(
effects = ~ 1 + time
, subject = ~ animal_id:breed
, type = c("vc", "un")[2]
)
, repeated = NULL
, method = c("REML", "ML", "Type3")
, ddfm = c("default", "contain", "residual", "betwithin", "satterth", "kr")
, type = 3
, class = NULL
, bound = TRUE
, control = list()
)
rc$covariance parameter subject estimate se z p
<char> <char> <num> <num> <num> <num>
1: UN(1,1) animal_id:breed 12.464572630 6.049874059 2.060303 1.968480e-02
2: UN(2,1) animal_id:breed -0.255121593 0.142828341 -1.786211 7.406506e-02
3: UN(2,2) animal_id:breed 0.005758072 0.003777797 1.524188 6.373091e-02
4: Residual 4.353042151 0.513010930 8.485282 1.075986e-17
cmp <-
vlm_compare(
CS = cs
, AR1 = ar
, RC = rc
, test = TRUE
)
cmpComparison of 3 models; smallest AIC: RC
Models
model method cov_parms m2ll AIC AICC BIC chisq df p p_boundary
CS REML 2 783.45 787.45 787.53 788.42 . . .
AR1 REML 2 793.14 797.14 797.21 798.11 . . .
RC REML 4 767.56 775.56 775.81 777.50 25.583 2 <1e-04 .
Sources: Duchateau, Janssen and Rowlands (1998); Akaike (1974); Schwarz (1978).
Full references: vlm_cite(x).
plot(cmp)
The book compared Models 1 and 2 only (the restricted log-likelihoods \(-391.73\) and \(-396.57\) favour compound symmetry); the random coefficients model has a much smaller AIC still, as the negative covariance of the intercepts and slopes (animals that start higher fall faster) suggests.
knitr::kable(x = ref[, c("source", "quantity", "value", "computed", "agrees")], digits = 4,
col.names = c("Source", "Quantity", "Book", "Computed", "Agrees"))| Source | Quantity | Book | Computed | Agrees |
|---|---|---|---|---|
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | animal variance | 5.4780 | 5.4785 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | residual variance | 4.9650 | 4.9649 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | fixed effect, row 1 | 35.0600 | 35.0602 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | standard error, row 1 | 1.0560 | 1.0558 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | containment df, row 1 | 10.0000 | 10.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | fixed effect, row 2 | -0.8420 | -0.8420 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | standard error, row 2 | 1.4930 | 1.4931 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | containment df, row 2 | 10.0000 | 10.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | fixed effect, row 4 | -0.4130 | -0.4134 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | standard error, row 4 | 0.0225 | 0.0225 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | containment df, row 4 | 154.0000 | 154.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | fixed effect, row 5 | -0.2760 | -0.2761 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | standard error, row 5 | 0.0225 | 0.0225 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | containment df, row 5 | 154.0000 | 154.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 3 | -0.6720 | -0.6720 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 4 | -2.4970 | -2.4967 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 5 | -2.3830 | -2.3827 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 2 | 3.0380 | 3.0379 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 6 | 0.0324 | 0.0324 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 1 | 2.4810 | 2.4811 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 7 | -1.6670 | -1.6671 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 8 | 2.6260 | 2.6264 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 9 | -1.6670 | -1.6671 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 10 | -0.3720 | -0.3723 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 11 | -1.7540 | -1.7543 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | prediction of animal 12 | 2.8340 | 2.8344 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | standard error of prediction | 1.0910 | 1.0912 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | df of the predictions | 154.0000 | 154.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | -2 restricted log-likelihood, compound symmetry | 783.4510 | 783.4514 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | Akaike’s criterion (SAS 6) | -393.7300 | -393.7257 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | Schwarz’s criterion (SAS 6) | -396.8300 | -396.8255 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | null model likelihood ratio chi-square | 80.8450 | 80.8447 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | Type 3 F, row 1 | 470.0500 | 470.0533 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | denominator df, row 1 | 154.0000 | 154.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | Type 3 F, row 2 | 0.3200 | 0.3180 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | denominator df, row 2 | 10.0000 | 10.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | Type 3 F, row 3 | 18.6200 | 18.6168 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 1 | denominator df, row 3 | 154.0000 | 154.0000 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | AR(1) correlation | 0.6210 | 0.6214 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | residual variance | 10.1900 | 10.1898 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | -2 restricted log-likelihood | 793.1390 | 793.1394 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | Akaike’s criterion (SAS 6) | -398.5700 | -398.5697 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | Schwarz’s criterion (SAS 6) | -401.6700 | -401.6696 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 2 | null model likelihood ratio chi-square | 71.1560 | 71.1566 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | UN(1,1) | 12.4600 | 12.4646 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | UN(2,1) | -0.2550 | -0.2551 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | UN(2,2) | 0.0058 | 0.0058 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | residual variance | 4.3530 | 4.3530 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | standard error of the intercept | 1.5010 | 1.5014 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | standard error of breed BO | 2.1230 | 2.1233 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | standard error of a slope | 0.0370 | 0.0375 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | random intercept of BO241 | -1.8360 | -1.8358 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | random slope of BO241 | 0.0563 | 0.0563 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | standard error of a random intercept | 1.6400 | 1.6376 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | standard error of a random slope | 0.0417 | 0.0417 | TRUE |
| Duchateau, Janssen and Rowlands (1998), Example 3.3, Model 3 | df of the predictions | 144.0000 | 144.0000 | TRUE |