## Analysis of variance with random effects

Back to course page

Some functions from the packages listed below are used inside this tutorial:

library("car")
library("nlme")


## Data

Data used in this tutorial come from a text by Luc Duchateau, Paul Janssen, and John Rowlands that is available at a webpage of International Livestock Research Institute in Nairobi, Kenya.

The following dataset corresponds to Table 1.2.2 from the above mentioned text and comes from a randomised block design study that aimed to find out whether low (L) or high (H) dose of a drug Berenil leads to different results in a treatment of trypanosomosis (parasite in cattle transmitted by a tse-tse fly). Three ($$I = 3$$) herds, each consisting of 10 cows were divided (each of them) into two parts (each consisting of 5 cows). One half of each herd was treated by a low dose of Berenil, the second half of the herd was treated by a high dose of Berenil. Severity of the disease was measured by so called packed cell volume (PCV, proportion of the volume of blood serum consisting of red blood cells). The primary outcome of the study was then increase of PCV (PCV.diff) over a specified time period.

dt122 <- data.frame(herd = gl(3, 10),
dose = rep(gl(2, 5, labels = c("L", "H")), 3),
PCV.diff = c(0.9, 2.0, 2.0, 2.2, 2.0,  7.4, 6.8, 7.1, 6.7, 7.9,
2.3, 2.7, 1.3, 1.6, 2.1,  6.6, 7.1, 6.2, 7.8, 7.2,
3.4, 2.7, 2.8, 3.0, 3.4,  8.8, 8.3, 7.9, 8.2, 8.1))
print(dt122)

##    herd dose PCV.diff
## 1     1    L      0.9
## 2     1    L      2.0
## 3     1    L      2.0
## 4     1    L      2.2
## 5     1    L      2.0
## 6     1    H      7.4
## 7     1    H      6.8
## 8     1    H      7.1
## 9     1    H      6.7
## 10    1    H      7.9
## 11    2    L      2.3
## 12    2    L      2.7
## 13    2    L      1.3
## 14    2    L      1.6
## 15    2    L      2.1
## 16    2    H      6.6
## 17    2    H      7.1
## 18    2    H      6.2
## 19    2    H      7.8
## 20    2    H      7.2
## 21    3    L      3.4
## 22    3    L      2.7
## 23    3    L      2.8
## 24    3    L      3.0
## 25    3    L      3.4
## 26    3    H      8.8
## 27    3    H      8.3
## 28    3    H      7.9
## 29    3    H      8.2
## 30    3    H      8.1


Object dt122 is a classical data.frame which can be worked with in a standard way.

## Exploratory plots

Some basic exploratory plots useful in a context of comparison of the effect of low/high dose of Berenil on PCV.diff, while taking into account a grouped structure of the data are the following:

COL <- rep(rep(c("red", "darkblue"), each = 5), 3)
BG <- rep(rep(c("pink", "skyblue"), each = 5), 3)
COL2 <- rep(c("pink", "skyblue"), 3)

par(mfrow = c(1, 2), bty = "n")
plot(PCV.diff ~ as.numeric(herd:dose), data = dt122, xaxt = "n", xlab = "", col = COL, bg = BG, pch = 23)
axis(1, at = 1:6, labels = paste("H", rep(1:3, each = 2), "/", rep(c("L", "H"), 3), sep = ""))
plot(PCV.diff ~ as.factor(herd:dose), data = dt122, xaxt = "n", xlab = "", col = COL2)
axis(1, at = 1:6, labels = paste("H", rep(1:3, each = 2), "/", rep(c("L", "H"), 3), sep = ""))


par(mfrow = c(1, 1))


What can be concluded from here?

## Grouped data

More possibilities to calculate descriptive statistics and draw plots are obtained if data are stored as a groupedData structure (capabilities of the nlme package):

gdt122 <- groupedData(formula = PCV.diff ~ 1 | herd / dose, data = dt122)
print(gdt122)

## Grouped Data: PCV.diff ~ 1 | herd/dose
## <environment: 0x98ecba4>
##    herd dose PCV.diff
## 1     1    L      0.9
## 2     1    L      2.0
## 3     1    L      2.0
## 4     1    L      2.2
## 5     1    L      2.0
## 6     1    H      7.4
## 7     1    H      6.8
## 8     1    H      7.1
## 9     1    H      6.7
## 10    1    H      7.9
## 11    2    L      2.3
## 12    2    L      2.7
## 13    2    L      1.3
## 14    2    L      1.6
## 15    2    L      2.1
## 16    2    H      6.6
## 17    2    H      7.1
## 18    2    H      6.2
## 19    2    H      7.8
## 20    2    H      7.2
## 21    3    L      3.4
## 22    3    L      2.7
## 23    3    L      2.8
## 24    3    L      3.0
## 25    3    L      3.4
## 26    3    H      8.8
## 27    3    H      8.3
## 28    3    H      7.9
## 29    3    H      8.2
## 30    3    H      8.1

summary(gdt122)

##  herd   dose      PCV.diff
##  1:10   L:15   Min.   :0.900
##  2:10   H:15   1st Qu.:2.225
##  3:10          Median :4.800
##                Mean   :4.883
##                3rd Qu.:7.350
##                Max.   :8.800

gsummary(gdt122)

##     herd dose PCV.diff
## 1/L    1    L     1.82
## 1/H    1    H     7.18
## 2/L    2    L     2.00
## 2/H    2    H     6.98
## 3/L    3    L     3.06
## 3/H    3    H     8.26

gsummary(gdt122, FUN = mean)

##     herd dose PCV.diff
## 1/L    1    L     1.82
## 1/H    1    H     7.18
## 2/L    2    L     2.00
## 2/H    2    H     6.98
## 3/L    3    L     3.06
## 3/H    3    H     8.26

gsummary(gdt122, FUN = sd)

##     herd dose  PCV.diff
## 1/L    1    L 0.5215362
## 1/H    1    H 0.4868265
## 2/L    2    L 0.5567764
## 2/H    2    H 0.6099180
## 3/L    3    L 0.3286335
## 3/H    3    H 0.3361547

plot(gdt122)


## Classical two-way ANOVA

### Two-way ANOVA with interaction

summary(m0 <- aov(PCV.diff ~ herd + dose + herd:dose, data = dt122))

##             Df Sum Sq Mean Sq F value   Pr(>F)
## herd         2   9.05    4.52  19.225 1.04e-05 ***
## dose         1 201.24  201.24 855.140  < 2e-16 ***
## herd:dose    2   0.18    0.09   0.387    0.683
## Residuals   24   5.65    0.24
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Do we need interaction?
• Provided interaction was significant, would it be possible to evaluate an effect of dose?

### Two-way ANOVA without interaction

summary(m0.noInt <- aov(PCV.diff ~ herd + dose, data = dt122))

##             Df Sum Sq Mean Sq F value   Pr(>F)
## herd         2   9.05    4.52   20.18 5.13e-06 ***
## dose         1 201.24  201.24  897.48  < 2e-16 ***
## Residuals   26   5.83    0.22
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(m0.noInt, type = "II")

## Anova Table (Type II tests)
##
## Response: PCV.diff
##            Sum Sq Df F value    Pr(>F)
## herd        9.049  2  20.177 5.133e-06 ***
## dose      201.243  1 897.482 < 2.2e-16 ***
## Residuals   5.830 26
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• We have balanced design, so type I and type II ANOVA tables are the same.
• Are there differences between doses?
• Are there differences between herds?
• How to interprete results of this model?
• Estimated effect of dose?
lm0.noInt <- lm(PCV.diff ~ herd + dose, data = dt122)
summary(lm0.noInt)

##
## Call:
## lm(formula = PCV.diff ~ herd + dose, data = dt122)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -1.0100 -0.2975  0.0350  0.3050  0.8100
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.9100     0.1729  11.046 2.57e-11 ***
## herd2        -0.0100     0.2118  -0.047    0.963
## herd3         1.1600     0.2118   5.478 9.57e-06 ***
## doseH         5.1800     0.1729  29.958  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4735 on 26 degrees of freedom
## Multiple R-squared:  0.973,  Adjusted R-squared:  0.9699
## F-statistic: 312.6 on 3 and 26 DF,  p-value: < 2.2e-16

confint(lm0.noInt)["doseH",]

##    2.5 %   97.5 %
## 4.824581 5.535419


## Mixed-effects two-way ANOVA, classical analysis using the decomposition of the sum of squares

We will now consider herds as a random sample (of size 3) from a population of herds and will treat corresponding parameters as random variables (random effects) rather than fixed quantities (fixed effects). In the following, let $$Y_{i,j,k}$$ denote PCV.diff for cow $$k$$ ($$k = 1,…,5 = K$$) in part of the herd $$i$$ ($$i = 1, 2, 3 = I$$) being treated by dose $$j$$ ($$j = 1, 2 = J$$). Fully, we could write $$k = k(i,j)$$ and $$j = j(i)$$ since with respect to the design of the study, cows are nested within a dose which is nested within a herd.

Classical mixed-effects ANOVA model corresponding to the design of this study is the following $Y_{i,j,k} = \mu + b_i^A + \beta_j^B + b_{i,j}^{AB} + \varepsilon_{i,j,k},$ where

• $$\beta_j^B$$: (fixed) effect of dose $$j$$ (some identif. constr. needed, e.g., with contr.treatment: $$\beta_1^B = 0$$, or with contr.sum: $$\beta_1^B + \beta_2^B = 0$$);
• $$b_i^A$$: (random) effects of herd, $$\mbox{E}(b_i) = 0$$, $$\mbox{var}(b_i) = d_1^2$$, $$i = 1, 2, 3$$;
• $$b_{i,j}^{AB}$$: (random) effects of dose (within a herd), $$\mbox{E}(b_{i,j}) = 0$$, $$\mbox{var}(b_{i,j}) = d_2^2$$, $$i = 1, 2, 3$$, $$j = j(i) = 1, 2$$;
• $$\varepsilon_{i,j,k}$$: random error = (random) effect of a cow, $$\mbox{E}(\varepsilon_{i,j,k}) = 0$$, $$\mbox{var}(\varepsilon_{i,j,k}) = \sigma^2$$, $$i = 1, 2, 3$$, $$j = j(i) = 1, 2$$, $$k = k(i,j) = 1, 2, 3, 4, 5$$.

All random terms are assumed to be mutually uncorrelated.

That is, variability of response is decomposed into three sources ($$\sigma^2$$ $$=$$ variability between cows coming from the same herd and being treated by the same dose, $$d_2^2$$ $$=$$ variability between doses inside a herd, $$d_1^2$$ $$=$$ variability between herds.

ANOVA table for the above model:

summary(m1 <- aov(PCV.diff ~ dose + Error(herd/dose), data = dt122))

##
## Error: herd
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  2  9.049   4.524
##
## Error: herd:dose
##           Df Sum Sq Mean Sq F value   Pr(>F)
## dose       1 201.24  201.24    2211 0.000452 ***
## Residuals  2   0.18    0.09
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24  5.648  0.2353


In the output:

• Error: herd, Mean Sq (Residuals) $$= \mbox{MS}_A$$;
• Error: herd:dose, Mean Sq (dose) $$= \mbox{MS}_B$$;
• Error: herd:dose, Mean Sq (Residuals) $$= \mbox{MS}_{AB}$$;
• Error: Within, Mean Sq (Residuals) $$= \mbox{MS}_e$$.

Can you find a test for significance of the fixed effect of dose in the above output? How to interprete the results of this test?

The mean squares that correspond to the random components of the model satisfy the following (try to derive as an optional homework): $\begin{array}{rcl} \mbox{E}(\mbox{MS}_A) & = & \sigma^2 \;+\; K\,d_2^2 \;+\; J\cdot K\,d_1^2, \ \mbox{E}(\mbox{MS}_{AB}) & = & \sigma^2 \;+\; K\,d_2^2, \ \mbox{E}(\mbox{MS}_{e}) & = & \sigma^2. \ \end{array}$ Using above expressions, unbiased estimators of all variance components can easily be derived.

Here is unbiased estimate of $$\sigma^2$$ and related estimate of $$\sigma$$:

MSe <- deviance(m1$Within) / m1$Within$df.residual sigmaHat <- sqrt(MSe) data.frame(MSe = MSe, sigmaHat = sigmaHat)  ## MSe sigmaHat ## 1 0.2353333 0.4851117  Here is unbiased estimate of $$d_1^2$$ and related estimate of $$d_1$$: I <- 3 J <- 2 K <- 5 MSA <- deviance(m1$herd) / m1$herd$df.residual
MSAB <- deviance(m1$"herd:dose") / m1$"herd:dose"\$df.residual
d1sqHat <- (MSA - MSAB) / (J * K)
d1Hat <- sqrt(d1sqHat)
data.frame(MSA = MSA, MSAB = MSAB, d1sqHat = d1sqHat, d1Hat = d1Hat)

##        MSA  MSAB   d1sqHat     d1Hat
## 1 4.524333 0.091 0.4433333 0.6658328


Here is unbiased estimate of $$d_2^2$$:

(d2sqHat <- (MSAB - MSe) / K)

## [1] -0.02886667

• Negative estimate of the variance component. How does it come?

## Linear mixed model, (RE)ML analysis

To fit a linear mixed model by (RE)ML method, capabilities of the package nlme can be used. In a context of the current study, be aware of the fact that asymptotic is somehow questionable here since only $$I = 3$$ herds are available.

The same model as above, now fitted by the REML method (default to lme function):

m2lme <- lme(PCV.diff ~ dose, random = ~1 | herd/dose, data = dt122)
summary(m2lme)

## Linear mixed-effects model fit by REML
##  Data: dt122
##        AIC      BIC    logLik
##   59.02353 65.68455 -24.51177
##
## Random effects:
##  Formula: ~1 | herd
##         (Intercept)
## StdDev:   0.6557517
##
##  Formula: ~1 | dose %in% herd
##          (Intercept)  Residual
## StdDev: 1.591771e-05 0.4735301
##
## Fixed effects: PCV.diff ~ dose
##                Value Std.Error DF   t-value p-value
## (Intercept) 2.293333 0.3978511 24  5.764301  0.0000
## doseH       5.180000 0.1729088  2 29.957998  0.0011
##  Correlation:
##       (Intr)
## doseH -0.217
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.17303689 -0.65651639  0.07550478  0.60397762  1.67043567
##
## Number of Observations: 30
## Number of Groups:
##           herd dose %in% herd
##              3              6

• $$\widehat{d}_1 = 0.6557517$$;
• $$\widehat{d}_2 = 1.591771e-05$$;
• $$\widehat{\sigma} = 0.4735301$$;
• Estimated effect of dose?
• P-value to evaluate the effect of dose?

#### Estimated fixed effects

fixef(m2lme)

## (Intercept)       doseH
##    2.293333    5.180000


#### Empirical Bayes estimates of random effects (BLUP's)

ranef(m2lme)

## Level: herd
##   (Intercept)
## 1  -0.3643349
## 2  -0.3738393
## 3   0.7381742
##
## Level: dose %in% herd
##       (Intercept)
## 1/L -6.158230e-10
## 1/H  4.011472e-10
## 2/L  4.548455e-10
## 2/H -6.751214e-10
## 3/L  1.609775e-10
## 3/H  2.739742e-10


#### (Approximate) confidence intervals for model parameters

intervals(m2lme)

## Approximate 95% confidence intervals
##
##  Fixed effects:
##                lower     est.    upper
## (Intercept) 1.472209 2.293333 3.114458
## doseH       4.436034 5.180000 5.923966
## attr(,"label")
## [1] "Fixed effects:"
##
##  Random Effects:
##   Level: herd
##                     lower      est.    upper
## sd((Intercept)) 0.2338472 0.6557517 1.838852
##   Level: dose
##                 lower         est. upper
## sd((Intercept))     0 1.591771e-05   Inf
##
##  Within-group standard error:
##     lower      est.     upper
## 0.3608343 0.4735301 0.6214231


Apparently, there are some problems in estimation of $$d_2^2$$. We are approaching a boundary of the parameter space. See also a negative value of its unbiased estimate.

#### ML method

Also the ML method can be used:

m2MLlme <- lme(PCV.diff ~ dose, random = ~1 | herd/dose, data = dt122, method = "ML")
summary(m2MLlme)

## Linear mixed-effects model fit by maximum likelihood
##  Data: dt122
##        AIC      BIC   logLik
##   57.06219 64.06818 -23.5311
##
## Random effects:
##  Formula: ~1 | herd
##         (Intercept)
## StdDev:   0.5291783
##
##  Formula: ~1 | dose %in% herd
##          (Intercept)  Residual
## StdDev: 1.113661e-05 0.4646783
##
## Fixed effects: PCV.diff ~ dose
##                Value Std.Error DF  t-value p-value
## (Intercept) 2.293333 0.3397556 24  6.74995  0.0000
## doseH       5.180000 0.1756319  2 29.49350  0.0011
##  Correlation:
##       (Intr)
## doseH -0.258
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.23260293 -0.66583662  0.07333894  0.59731185  1.68408520
##
## Number of Observations: 30
## Number of Groups:
##           herd dose %in% herd
##              3              6

intervals(m2MLlme)

## Approximate 95% confidence intervals
##
##  Fixed effects:
##                lower     est.    upper
## (Intercept) 1.615889 2.293333 2.970777
## doseH       4.449941 5.180000 5.910059
## attr(,"label")
## [1] "Fixed effects:"
##
##  Random Effects:
##   Level: herd
##                     lower      est.    upper
## sd((Intercept)) 0.2233861 0.5291783 1.253568
##   Level: dose
##                 lower         est. upper
## sd((Intercept))     0 1.113661e-05   Inf
##
##  Within-group standard error:
##     lower      est.     upper
## 0.3558910 0.4646783 0.6067192


#### Slightly different mixed models

The three models below are almost the same. Nevertheless, they (slightly) differ. How?

gdt122.2 <- groupedData(formula = PCV.diff ~ 1 | herd/dose, data = dt122)
gdt122.34 <- groupedData(formula = PCV.diff ~ 1 | herd, data = dt122)
m2 <- lme(PCV.diff ~ dose, random = ~1 | herd/dose,   data = gdt122.2)
m3 <- lme(PCV.diff ~ dose, random = ~1 + dose | herd, data = gdt122.34)
m4 <- lme(PCV.diff ~ dose, random = pdDiag(~ 1 + dose), data = gdt122.34)

summary(m2)

## Linear mixed-effects model fit by REML
##  Data: gdt122.2
##        AIC      BIC    logLik
##   59.02353 65.68455 -24.51177
##
## Random effects:
##  Formula: ~1 | herd
##         (Intercept)
## StdDev:   0.6557517
##
##  Formula: ~1 | dose %in% herd
##          (Intercept)  Residual
## StdDev: 1.591771e-05 0.4735301
##
## Fixed effects: PCV.diff ~ dose
##                Value Std.Error DF   t-value p-value
## (Intercept) 2.293333 0.3978511 24  5.764301  0.0000
## doseH       5.180000 0.1729088  2 29.957998  0.0011
##  Correlation:
##       (Intr)
## doseH -0.217
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.17303689 -0.65651639  0.07550478  0.60397762  1.67043567
##
## Number of Observations: 30
## Number of Groups:
##           herd dose %in% herd
##              3              6

summary(m3)

## Linear mixed-effects model fit by REML
##  Data: gdt122.34
##        AIC      BIC    logLik
##   61.02353 69.01676 -24.51177
##
## Random effects:
##  Formula: ~1 + dose | herd
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev       Corr
## (Intercept) 6.557517e-01 (Intr)
## doseH       2.439662e-05 0
## Residual    4.735301e-01
##
## Fixed effects: PCV.diff ~ dose
##                Value Std.Error DF   t-value p-value
## (Intercept) 2.293333 0.3978511 26  5.764301       0
## doseH       5.180000 0.1729088 26 29.957998       0
##  Correlation:
##       (Intr)
## doseH -0.217
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.17303689 -0.65651640  0.07550478  0.60397763  1.67043567
##
## Number of Observations: 30
## Number of Groups: 3

summary(m4)

## Linear mixed-effects model fit by REML
##  Data: gdt122.34
##        AIC      BIC    logLik
##   59.02353 65.68455 -24.51177
##
## Random effects:
##  Formula: ~1 + dose | herd
##  Structure: Diagonal
##         (Intercept)      doseH  Residual
## StdDev:   0.6557517 2.5681e-05 0.4735301
##
## Fixed effects: PCV.diff ~ dose
##                Value Std.Error DF   t-value p-value
## (Intercept) 2.293333 0.3978511 26  5.764301       0
## doseH       5.180000 0.1729088 26 29.957998       0
##  Correlation:
##       (Intr)
## doseH -0.217
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.17303689 -0.65651639  0.07550478  0.60397762  1.67043567
##
## Number of Observations: 30
## Number of Groups: 3

ranef(m2)

## Level: herd
##   (Intercept)
## 1  -0.3643349
## 2  -0.3738393
## 3   0.7381742
##
## Level: dose %in% herd
##       (Intercept)
## 1/L -6.158230e-10
## 1/H  4.011472e-10
## 2/L  4.548455e-10
## 2/H -6.751214e-10
## 3/L  1.609775e-10
## 3/H  2.739742e-10

ranef(m3)

##   (Intercept)         doseH
## 2  -0.3738393 -2.431589e-09
## 1  -0.3643349  1.181571e-10
## 3   0.7381742  2.313431e-09

ranef(m4)

##   (Intercept)         doseH
## 2  -0.3738393 -1.757296e-09
## 1  -0.3643349  1.044159e-09
## 3   0.7381742  7.131364e-10


m2 is, in fact, a model with two-levels of nesting, m3 and m4 are specified as classical (single-level) linear mixed models. Can you write mathematically the three models? Models m2 and m3 or m4 have (among other things) somehow different interpretation of the “interaction” random effects “$$b_{i,j}^{AB}$$”.

Back to course page