Abstract

We explore, review and analyze the ergoStool data in the first chapter of Mixed Effects Models in S and S-Plus by Pinheiro and Bates (2000). The ergoStool data provide an example of a randomized block design that may be analyzed within the framework of mixed effects models, those which contain both fixed and random factors. The data are available from the nlme package. We comment on some observations/insights regarding the use of R Markdown and on some differences found in generating PDF and HTML files.

Opening Observations

This document reflects initial attempts using RMarkdown for Pinheiro and Bates’ ergoStool example. Much of the text is directly copied from Pinheiro and Bates (2000), designated as PB or P&B throughout. Commentary is provided in connection with certain RMarkdown commands, etc. The RMarkdown (.Rmd) files for this example (both HTML and PDF) are available on request from W. Greg Alvord, greg.alvord@nih.gov.

Randomized Blocks Design Example, Pinheiro & Bates, page 12

Introduction

We address the second example discussed by Pinheiro & Bates (2000). From P&B, page 12, “A randomized block design is a type of experiment in which there are two classification factors: an experimental factor for which we use fixed effects and a blocking factor for which we use random effects.” The data are contained in the object ergoStool. We show a summary of the data and the default plot.

Here is a summary of the data.

##      effort      Type      Subject  
##  Min.   : 7.00   T1:9   8      : 4  
##  1st Qu.: 8.00   T2:9   5      : 4  
##  Median :10.00   T3:9   4      : 4  
##  Mean   :10.25   T4:9   9      : 4  
##  3rd Qu.:12.00          6      : 4  
##  Max.   :15.00          3      : 4  
##                         (Other):12

FIGURE 1.5 (P&B, p. 13) displays the default plot of the data.

FIGURE 1.5 (p.13) in P&B. Effort required (Borg scale) to arise from a stool for nine different subjects each using four different types of stools. Different colors, show in the key at the top of the plot, are used for different types of stools.

FIGURE 1.5 (p.13) in P&B. Effort required (Borg scale) to arise from a stool for nine different subjects each using four different types of stools. Different colors, show in the key at the top of the plot, are used for different types of stools.

‘The experimenters recorded the effort required by each of nine different subjects to arise from four different types of stools, which we will label as T1, T2, T3, and T4.’ We will want to compare these four particular types of stools; hence, we will use fixed effects for the Type factor. The nine different subjects represent a sample from the population about which we wish to make inferences so we use random effects to model the Subject factor.’ (P&B, p. 13)

Comment: Pinheiro and Bates do a nice job of introducing the reader to the idea of a randomized block (or randomized blocks) design. They emphasize that the blocking factor is a random variable. An investigation of statistics texts and experimental design texts from 40 or 50 years ago reveals that there was (and still is) some confusion about what a randomized blocks design is. In the older texts, the idea of a blocking factor was sometimes introduced as representing a random variable, but not always. In older texts, it was rarely a introduced as a subject. Frequently, it was designated as a partition of land (a block or plot, say), a classroom (with one or more students), a litter (with one or more animals), a type of fertilizer, a blend of some crop, or contiguous periods of time.

Analysis

From FIGURE 1.5 (P&B, p. 13) we can see systematic types of differences among the different stool types. For example, stool T2 required the greatest effort across Subjects. On the other hand, stool T1 consistently required a low amount of energy across Subjects. Subject is “a blocking factor because it represents a known source of variability in the experiment.” (p. 13)

We next compare (visually) the magnitudes of the effects of Type and Subject factors.

FIGURE 1.6 (p.14) in P&B. Design plot for the data in the stool ergometric experiment. The mean of the response (effort) is plotted for each level of each of the factors Type and Subject.

FIGURE 1.6 (p.14) in P&B. Design plot for the data in the stool ergometric experiment. The mean of the response (effort) is plotted for each level of each of the factors Type and Subject.

“The variability associated with the Type factor (a fixed factor) is comparable to the variability associated with the Subject factor (a random factor). We also see [FIGURE 1.6 in P&B, p. 14] that the average effort according to stool type is in the order T1 \(\leq\) T4 \(\leq\) T3 \(\leq\) T2.” (p. 13)

1.2.1 Choosing Contrasts for Fixed-Effects Terms

[We quote directly from Pinheiro & Bates, p. 14.] A model with fixed effects \(\beta_{j}\) for the Type factor and random effects \(b_{i}\) for the Subject factor could be written

\[ \hspace{33 mm} y_{ij} = \beta_{j} + b_{i} + \epsilon_{ij}, \hspace{5 mm} i = 1,...,9, \hspace{5 mm} j = 1, ..., 4, \hspace{33 mm} \mbox{PB (1.6)} \]

with random parameters distributed as

\[ \hspace{49 mm} b_{i} \sim \mathcal{N}(0, \sigma^{2}_{b}), \hspace{4 mm} \epsilon_{ij} \sim \mathcal{N} (0, \sigma^{2}) \hspace{49 mm} \mbox{PB (1.6)} \]

In matrix notation, this can be written, equivalently, as

\[ \boldsymbol{y_{i}} = \boldsymbol{X_{i}}\boldsymbol{\beta} + \boldsymbol{Z_{i}}b_{i} + \boldsymbol{\epsilon_{i}}, \hspace{5 mm} i = 1,...,9, \]

with random parameters distributed as

\[ b_{i} \sim \mathcal{N}(0, \sigma^{2}_{b}), \hspace{4 mm} \boldsymbol{\epsilon}_{i} \sim \mathcal{N} (\boldsymbol{0}, \sigma^{2}\boldsymbol{I}) \]

where, for \(i = 1,...,9\)

\[ \boldsymbol{y}_{i} = \left[ \begin{array}{c} y_{i1} \\ y_{i2} \\ y_{i3} \\ y_{i4} \end{array} \right] , \hspace{5 mm} \boldsymbol{X}_{i} = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right] , \hspace{5 mm} \boldsymbol{Z}_{i} = \boldsymbol{1} = \left[ \begin{array}{c} 1 \\ 1 \\ 1 \\ 1 \end{array} \right] , \hspace{5 mm} \boldsymbol{y}_{i} = \left[ \begin{array}{c} \epsilon_{i1} \\ \epsilon_{i2} \\ \epsilon_{i3} \\ \epsilon_{i4} \end{array} \right] \]

[Quoting from P&B] “This form of fixed-effects matrix \(\boldsymbol{X}_{i}\) is sometimes called the cell means form because the \(j\)th component of \(\boldsymbol{\beta}\) represents what would be the mean effort to arise from the \(j\)th type of stool if the whole population were tested.”

The \(\beta_{j}\) have a simple interpretation, but are not convenient to use when assessing differences between stool types, Type [recall: Type is a fixed factor]. To make it easier to assess these differences we use an alternative form of the \(\boldsymbol{X}_{i}\) matrices with one column expressing the “overall mean” or reference level and three columns representing changes among the types of stools. The three columns representing the changes are called contrasts. There are several different choices available for these contrasts (Venables and Ripley, Sect. 6.2). In R, the default contrasts for an unordered factor are the “treatment” contrasts.

##    T2 T3 T4
## T1  0  0  0
## T2  1  0  0
## T3  0  1  0
## T4  0  0  1

The \(\boldsymbol{X}_{i}\) matrices can be displayed with the model.matrix() function. To save space we show the \(\boldsymbol{X}_{1}\) matrix only (i.e., for the first subject). (We use option echo = TRUE to see the logic of how this is done.)

tv.Subject.1 <- with(ergoStool, Subject == '1') # create a 'truth' vector
ergoStool1 <- ergoStool[tv.Subject.1,]
model.matrix( effort ~ Type, data = ergoStool1) # X matrix for Subject 1
##   (Intercept) TypeT2 TypeT3 TypeT4
## 1           1      0      0      0
## 2           1      1      0      0
## 3           1      0      1      0
## 4           1      0      0      1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Type
## [1] "contr.treatment"

Fitting the model in this form with lme produces

fm1Stool <- lme( effort ~ Type, data = ergoStool, random = ~ 1 | Subject )
summary(fm1Stool)
## Linear mixed-effects model fit by REML
##  Data: ergoStool 
##        AIC      BIC    logLik
##   133.1308 141.9252 -60.56539
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:    1.332465 1.100295
## 
## Fixed effects: effort ~ Type 
##                Value Std.Error DF   t-value p-value
## (Intercept) 8.555556 0.5760123 24 14.853079  0.0000
## TypeT2      3.888889 0.5186838 24  7.497610  0.0000
## TypeT3      2.222222 0.5186838 24  4.284348  0.0003
## TypeT4      0.666667 0.5186838 24  1.285304  0.2110
##  Correlation: 
##        (Intr) TypeT2 TypeT3
## TypeT2 -0.45               
## TypeT3 -0.45   0.50        
## TypeT4 -0.45   0.50   0.50 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.80200345 -0.64316591  0.05783115  0.70099706  1.63142054 
## 
## Number of Observations: 36
## Number of Groups: 9

[Quoting P&B, p. 16] “By convention, the coefficient corresponding to the first column in the \(\boldsymbol{X}_{i}\), is called the intercept, (Intercept). Since we are using treatment contrasts here (the default in R as opposed to Helmert contrasts in S), the (Intercept) is actually the mean of all responses for Type T1. The other effects are quantities that need to be added to, or subtracted from, the Intercept.” We confirm this with some simple R code.

fixef(fm1Stool) # fixed effects for object fm1Stool
## (Intercept)      TypeT2      TypeT3      TypeT4 
##   8.5555556   3.8888889   2.2222222   0.6666667
with(ergoStool, tapply(effort, Type, mean)) # calculate means for Type factor
##        T1        T2        T3        T4 
##  8.555556 12.444444 10.777778  9.222222

We can see that if we add the ‘effect’ for T2 to the Intercept, i.e., 8.5555556 + 3.8888889, or 12.4444444 \(\equiv\) 12.4444444. Similarly, if we add the ‘effect’ for T3 to the Intercept, i.e., 8.5555556 + 2.2222222, or 10.7777778 \(\equiv\) 10.7777778. The same operation holds true for T4 as well.

The significance of the Type term is assessed with anova function.

anova( fm1Stool )
##             numDF denDF  F-value p-value
## (Intercept)     1    24 455.0075  <.0001
## Type            3    24  22.3556  <.0001

The row labelled Type in the anova table represents a test of the hypothesis

\[ \mbox{H}_{0}: \beta_{2} = \beta_{3} = \beta_{4} = 0, \]

which is equivalent to reducing model (P&B (1.6)) to

\[ \boldsymbol{y}_{i} = \boldsymbol{1}\beta + \boldsymbol{Z}_{i}b_{i} + \boldsymbol{\epsilon}_{i}, \hspace{5 mm} i = 1,...,9, \hspace{5 mm} b_{i} \sim \mathcal{N}(0, \sigma^{2}_{b}), \hspace{5 mm} \epsilon_{i} \sim \mathcal{N} (\boldsymbol{0}, \sigma^{2}\boldsymbol{I}) \]

Reiterating some of the points made on page 19,

  • The overall effect of the factor should be assessed with anova, not by examining the t-value’s or p-values’s associated with the fixed effects parameters. The anova does not depend on the choice of contrasts as long as the intercept term is retained in the model.

  • Interpretation of the parameter estimates for a fixed-effects term depends on the contrasts being used.

  • For REML estimation, likelihood-ratio tests or comparisons of AIC or BIC require the same fixed-effects structure and the same choice of contrasts in all models.

  • The “cell means” parameters can be estimated by adding -1 to amodel formula but this will usually make the results from anova meaningless.

1.2.2 Examining the Model

P&B recommend examining the fitted model, both numerically and graphically.

The intervals function provides an indication of the precision of the estimates of the variance components.

intervals( fm1Stool )
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                  lower      est.    upper
## (Intercept)  7.3667247 8.5555556 9.744386
## TypeT2       2.8183781 3.8888889 4.959400
## TypeT3       1.1517114 2.2222222 3.292733
## TypeT4      -0.4038442 0.6666667 1.737177
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Subject 
##                    lower     est.    upper
## sd((Intercept)) 0.749509 1.332465 2.368835
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.8292494 1.1002946 1.4599324

The estimate of \(\sigma\) is \(\hat\sigma\) = 1.1002946. The estimate of \(\sigma_{b}\) is \(\hat\sigma_{b}\) = 1.332465.

P&B (p. 20) state that “\(\sigma\) is estimated relatively precisely, whereas [the estimate of] \(\sigma_{b}\), or \(\hat\sigma_{b}\), can vary by a factor of about 3, which is a factor of 9 if we express the estimates as variances.”

Another useful function to produce both standard deviations and variances of the random effects in a model isVarCorr. Here is the VarCorr output for the fm1Stool object.

VarCorr(fm1Stool)
## Subject = pdLogChol(1) 
##             Variance StdDev  
## (Intercept) 1.775463 1.332465
## Residual    1.210648 1.100295

FIGURE 1.7 in P&B shows a plot of the standardized residuals for the fit fm1Stool.

FIGURE 1.7. Standardized residuals versus the fitted values for the REML fit of a random-effects model to the data in the ergometrix experiment on types of stools.

FIGURE 1.7. Standardized residuals versus the fitted values for the REML fit of a random-effects model to the data in the ergometrix experiment on types of stools.

[P&B, p 20] “The plot of the standardized residuals vs. the fitted values, shown in Figure 1.7, does not indicate a violation of the assumption of constant variance for the \(\epsilon_{ij}\) terms.” [Note data in ergoStool ‘groupedData’ data frame.]

[Quoting P&B, pp 21-22] “Figure 1.7 shows the overall behavior of the residuals relative to the fitted values. It may be more informative to examine this behavior according to the Subject factor or according to the Type factor, which can be done by providing an explicit formula to the plot method for the fitted lme model. The formula can use functions such as resid or fitted applied to the fitted model. As a shortcut, a”." appearing in the formula is interpreted as the fitted model object itself. Thus, to plot the standardized or “Pearson,” residuals versus the fitted values by Subject we use" [. . . the code in the code chunk below, to produce Figure 4, or FIGURE 1.8 in P&B.]

plot( fm1Stool,
      form = resid(., type = 'p') ~ fitted(.) | Subject, 
      between = list(x = .5, y = .3), pch = 16, abline = 0 ) # pch works here
P&B FIGURE 1.8. Standardized residuals versus the fitted values by Subject for the REML fit of a random effects model to the dat in the ergometric experiment on types of stools.

P&B FIGURE 1.8. Standardized residuals versus the fitted values by Subject for the REML fit of a random effects model to the dat in the ergometric experiment on types of stools.

The residuals in both plots are homogeneously distributed, indicating a nice fit produced from the fm1Stool model.

Comments for RMarkdown Users

Many \(\LaTeX\) \(\hspace{.1 mm}\) commands work when the .Rmd file is generating a PDF document, but not when it is generating an HTML document. Thus, we have chosen to use the R Markdown choice whenever possible. As an example, the syntax $$ y = \beta_{0} + \beta_{1}x $$ produces an equation, here, … \[ y = \beta_{0} + \beta_{1}x \] as expected. The \(\LaTeX\) \(\hspace{.1 mm}\) syntax \begin{equation} y = \beta_{0} + \beta_{1}x \end{equation}, here

numbers the equation whereas the syntax $$ y = \beta_{0} + \beta_{1}x $$ does not. If the .Rmd file is used to generate an HTML file, the syntax \[ y = \beta_{0} + \beta_{1}x \] works but the “\begin{equation} y = \beta_{0} + \beta_{1}x \end{equation}” syntax, here . . . \begin{equation} y = \beta_{0} + \beta_{1}x \end{equation} is ignored.

One of the most useful commands we have discovered is ‘boldsymbol’, preceded with a \ which allows variables, such as $\boldsymbol{y}_i$ or $\boldsymbol{X}_{i}$ to be typeset in R Markdown. For example the command ‘boldsymbol’ between dollar signs, (i.e., $’s), produces \(\boldsymbol{X}_{i}\).

In this document we have used several “Inline R Code” calculations such as 'r fm1Stool$sigma' where the quotes we see here, in plain text, are ‘backquotes’, above the tab key on the keyboard.

References

Pinheiro JC and Bates DM. (2000) Mixed-Effects Models in S and S-Plus, Springer

Venables WN and Ripley BD. (2002) Modern Applied Statistics with S 4th ed., Springer