Fish net data from Dr. Longnecker’s Notes

Strength of fish net material made from four randomly chosen machines.

y <- c(128,127,129,126,128,121,120,123,122,125,
       126,125,127,125,124,125,126,129,128,127)
machine <- as.factor(c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4))

Is there significant machine-to-machine variability?

One-way random effects model

Suppose

\[ Y_{ij} = \mu + A_i + \varepsilon_{ij}, \quad i = 1,\dots,a, \quad j = 1,\dots,n_{i} \] where

• \(Y_{ij}\) is the reponse of EU \(j\) in treatment group \(i\).
• the \(A_i\) are independent \(\mathcal{N}(0,\sigma_A^2)\) rvs called random effects.
• the \(\varepsilon_{ij}\) are independent \(\mathcal{N}(0,\sigma^2_\varepsilon)\).
• the \(A_i\) and the \(\varepsilon_{ij}\) are independent of each other.
• \(\mu\) is the overall mean
• we call \(\sigma^2_A\) and \(\sigma_\varepsilon^2\) variance components.

Until now we have studied fixed-effects models.

Assume a balanced design, i.e. \(n_1 = \dots = n_a = n\).

Goals in the one-way random effects model

In the one-way random effects model we wish to

1. Visualize the data.
2. Decompose the variability in the \(Y_{ijk}\) into its sources.
3. Test \(H_0\): \(\sigma_A^2 = 0\) versus \(H_1\): \(\sigma_A^2 > 0\).
4. Estimate the variance components \(\sigma_A^2\) and \(\sigma_\varepsilon^2\).
5. Estimate the overall mean \(\mu\).
6. “Predict” the realized values of \(A_1,\dots,A_a\).
7. Check whether the model assumptions are satisfied.

Fish net material

Side-by-side boxplots are a natural choice if the \(n_i\) are not too small.

boxplot(y~machine)

Fish net material

Dotplots are a better when the \(n_i\) are small.

stripchart(y~machine, vertical=T, method="jitter", xlab = "machine")

Sums of squares for one-way random effects model

We have the decomposition \(\operatorname{SS}_{\operatorname{Tot}}= \SStrt + \operatorname{SS}_{\operatorname{Error}}\).

Fish net material

Easiest way to obtain sums of squares:

lm_out <- lm(y ~ machine)

y.. <- mean(y)
SStot <- sum((y - y..)**2)
SSerror <- sum(lm_out$residuals**2)
SStrt <- SStot - SSerror

We obtain \(\operatorname{SS}_{\operatorname{Tot}}= 122.95\), \(\SStrt = 87.75\), and \(\operatorname{SS}_{\operatorname{Error}}= 35.2\).

ANOVA table for one-way random effects model

Obtain the MS values by dividing the SS values by the Df values.

We reject \(H_0\): \(\sigma_A^2 = 0\) if \(\operatorname{F}_{\operatorname{test}}> F_{a-1,a(n-1),\alpha}\).

The corresponding p-value is \(P(F > \operatorname{F}_{\operatorname{test}})\), where \(F \sim F_{a-1,a(n-1),\alpha}\).

Discuss: It can be shown that \(\SStrt\) and \(\operatorname{SS}_{\operatorname{Error}}\) are independent and \[ \frac{\operatorname{SS}_{\operatorname{Trt}}}{n\sigma_A^2 + \sigma_\varepsilon^2} \sim \chi^2_{a-1} \quad \text{ and } \quad \frac{\operatorname{SS}_{\operatorname{Error}}}{\sigma_\varepsilon^2} \sim \chi_{a(n-1)}^2. \] Show that \(\displaystyle \operatorname{F}_{\operatorname{test}}= \frac{\MStrt}{\MSerror} \sim F_{a-1,a(n-1)}\) under \(H_0\): \(\sigma_A^2 = 0\).

Fish net material

We obtain the ANOVA table just as we did in the fixed-effects case.

anova(lm_out)

Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value  Pr(>F)    
machine    3   87.7    29.2    13.3 0.00013 ***
Residuals 16   35.2     2.2                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The small p-value indicates strong evidence against \(H_0\): \(\sigma^2_A = 0\).

Expected MS values in one-way random-effects model

In models with random effects, we like to track how the variance components \(\sigma_A^2\) and \(\sigma_\varepsilon^2\) affect the expected values of the mean squares.

Method of moments for estimating \(\sigma_A^2\) and \(\sigma_\varepsilon^2\)

The method of moments (MoM) is one way to estimate the variance components.

1. Set the observed mean squares equal to the expected mean squares: \[ \MStrt \overset{\operatorname{set}}{=} n\sigma^2_A + \sigma_\varepsilon^2 \quad \text{ and } \quad \MSerror \overset{\operatorname{set}}{=} \sigma_\varepsilon^2 \]
2. Solve these equations for \(\sigma_A^2\) and \(\sigma_\varepsilon^2\), which gives \[ \dot \sigma_A^2 = \frac{\MStrt - \MSerror}{n} \quad \text{ and } \quad \dot \sigma_\varepsilon^2 = \MSerror. \]

MoM method is deprecated because it is possible to get \(\dot \sigma_A^2 < 0\).

Fish net material

Compute the MoM estimators \(\dot \sigma_A^2\) and \(\dot \sigma_\varepsilon^2\) for the fish net data.

n <- 5
a <- 4
MStrt <- SStrt / (a - 1)
MSerror <- SSerror / (a*(n-1))

sigma_e_mom <- sqrt(MSerror)
sigma_A_mom <- sqrt((MStrt - MSerror)/n)

We obtain \(\dot \sigma_A^2 = 5.41\) and \(\dot \sigma_\varepsilon^2 = 2.2\).

REML for estimating the variance components

• REML stands for restricted maximum likelihood.
• It is a recipe (beyond course) for estimating variance components.
• In the one-way random-effects model it results in \[ \hat \sigma_A^2 = \left\{ \begin{array}{ll} \dot \sigma_A^2& \text{if } \dot \sigma_A^2 \geq 0\\ 0 & \text{if } \dot \sigma_A^2 < 0 \end{array}\right. \text{ and } ~ \hat \sigma_\varepsilon^2 = \left\{ \begin{array}{ll} \dot \sigma_\varepsilon^2& \text{if } \dot \sigma_A^2 \geq 0\\ \dfrac{\operatorname{SS}_{\operatorname{Tot}}}{an - 1} & \text{if } \dot \sigma_A^2 < 0 \end{array}\right. \]
• Will always return a nonnegative estimate of \(\sigma_A^2\).
• Same as MoM estimator if \(\dot \sigma_A^2 > 0\).
• The “state of the art” is to use REML estimation.

Fish net material

Use R package lmerTest to compute the REML estimates of \(\sigma_A^2\) and \(\sigma_\varepsilon^2\).

# first time run install.packages("lmerTest")
library(lmerTest) # may need to run install.packages("Matrix") for this to work
lmer_out <- lmer(y ~ 1 + (1 | machine))
summary(lmer_out)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ 1 + (1 | machine)

REML criterion at convergence: 79.7

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.653 -0.688 -0.102  0.485  1.718 

Random effects:
 Groups   Name        Variance Std.Dev.
 machine  (Intercept) 5.41     2.33    
 Residual             2.20     1.48    
Number of obs: 20, groups:  machine, 4

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   125.55       1.21   3.00     104    2e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimating \(\mu\) and predicting \(A_1,\dots,A_a\).

• We estimate the overall mean \(\mu\) with \(\bar Y_{..}\).
• We also wish to guess at the realized values of \(A_1,\dots,A_a\).
• Call this prediction instead of estimation, because \(A_1,\dots,A_a\) are random variables rather than parameters with fixed values.
• A recipe (beyond course) for the best linear unbiased predictor yields \[ \hat A_i = \left(\frac{n\hat \sigma_A^2}{n\hat \sigma_A^2 + \hat \sigma_\varepsilon^2}\right) (\bar Y_{i.} - \bar Y_{..}) \quad \text{ for } i =1,\dots,a. \]

Fish net material

Obtain the predicted values of \(A_1\), \(A_2\), \(A_3\), and \(A_4\) for the fish net data.

yi. <- aggregate(y, by = list(machine = machine), mean)$x
# use MoMs, since sA2_dot > 0
Ahat <- (n * sigma_A_mom**2)/(n * sigma_A_mom**2 + sigma_e_mom**2) * (yi. - y..)
Ahat

[1]  1.896 -3.098 -0.139  1.341

Fish net material

Can use predict() on the lmer() output to obtain the \(\hat A_i\) values:

predict(lmer_out,random.only = TRUE) # gives estimates of the A_i for each obs

     1      2      3      4      5      6      7      8      9     10     11 
 1.896  1.896  1.896  1.896  1.896 -3.098 -3.098 -3.098 -3.098 -3.098 -0.139 
    12     13     14     15     16     17     18     19     20 
-0.139 -0.139 -0.139 -0.139  1.341  1.341  1.341  1.341  1.341 

Fish net material

boxplot(y~machine)
points(yi.~as.factor(levels(machine)),col="red",pch=19)
points((y.. + Ahat)~as.factor(levels(machine)),col="blue",pch=18)
abline(h = y..,col = "red")

Fish net material

stripchart(y~machine, vertical = T, method = "jitter", xlab = "machine")
points(yi.~as.factor(levels(machine)),col="red",pch=19)
points((y.. + Ahat)~as.factor(levels(machine)),col="blue",pch=18)
abline(h = y..,col = "red")

Confidence interval for the overall mean

• Under the one-way random effects model, we can show \[ \bar Y_{..} \sim \mathcal{N}\Big(\mu,\frac{n\sigma_A^2 + \sigma_\varepsilon^2}{na}\Big). \]
• From here, using \(\displaystyle \frac{\SStrt}{n\sigma_A^2 + \sigma_\varepsilon^2}\sim \chi^2_{a-1}\), we have \[ \frac{\bar Y_{..} - \mu}{\sqrt{\MStrt / (na)}} \sim t_{a-1}. \]
• Thus a \((1-\alpha)100\%\) CI for \(\mu\) can be obtained as \[ \bar Y_{..} \pm t_{a-1,\alpha/2}\sqrt{\frac{\MStrt}{an}}. \]
• confint() function on lmer() output uses another method.

Fish net material

Construct a \(95\%\) confidence interval for \(\mu\) in fish net example.

alpha <- 0.05
tval <- qt(1-alpha/2,a-1)
me <- tval * sqrt(MStrt/(a*n))
lo <- y.. - me
up <- y.. + me

We obtain the interval \([121.701,129.399]\).

Fish net material

Default confidence intervals from confint() function.

confint(lmer_out)

              2.5 % 97.5 %
.sig01        0.968   5.09
.sigma        1.087   2.19
(Intercept) 122.890 128.21

Fish net material

Wald-type CI \(Y_{..} \pm z_{\alpha/2}\sqrt{\MStrt/(na)}\) on fish net data.

zval <- qnorm(1-alpha/2)
me <- zval * sqrt(MStrt/(a*n))
lo <- y.. - me
up <- y.. + me

We obtain the interval \([123.18,127.92]\).

Compare to:

confint(lmer_out,method="Wald")

            2.5 % 97.5 %
.sig01         NA     NA
.sigma         NA     NA
(Intercept)   123    128

Checking model assumptions

Validity of the methods in these slides depends on these assumptions:

1. The responses are normally distributed around the treatment means.
   (Check QQ plot of residuals)

2. The treatment means are normally distributed.
   (The number of groups \(a\) is often too small for this to be checked)

3. The response has the same variance in all treatment groups.
   (Check residuals vs fitted values plot)

4. The response values in each treatment group are independent.
   (No way to check; must trust experimental design)

Fitted values and residuals

• Build fitted values from \(\hat \mu\) and the predictions \(\hat A_1,\dots,\hat A_a\) as \[ \hat Y_{ij} = \hat \mu + \hat A_i, \quad \text{ for } i =1 ,\dots, a, \quad j=1,\dots,n. \]

• Our residuals then become \[ \hat \varepsilon_{ij} = Y_{ij} - \hat Y_{ij}, \quad \text{ for } i =1 ,\dots, a, \quad j=1,\dots,n. \]

• These residuals are not simply the deviations of the responses around their treatment group means (as in fixed effects models)!

• So \(\operatorname{SS}_{\operatorname{Error}}\) is not the sum of squared residuals here.

Fish net material

Use predict() on the lmer() output to obtain the \(\hat Y_{ij} = \hat \mu + \hat A_i\).

predict(lmer_out)

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
127 127 127 127 127 122 122 122 122 122 125 125 125 125 125 127 127 127 127 127 

Fish net material

Check the residuals vs fitted values plot: Use plot() on lmer() output.

plot(lmer_out)

Fish net material

Check the normal q-q plot of the residuals:

yhat <- predict(lmer_out)
ehat <- y - yhat
qqnorm(ehat)

Triglyceride level data from Kuehl (2000)

Triglyceride levels in serum samples measured with 4 machines on 4 days.

tg <- c(142.3,144.0,148.6,146.9,142.9,147.4,133.8,133.2,
        134.9,146.3,145.2,146.3,125.9,127.6,108.9,107.5,
        148.6,156.5,148.6,153.1,135.5,138.9,132.1,149.7,
        152.0,151.4,149.7,152.0,142.9,142.3,141.7,141.2)
days <- as.factor(c(rep(1,8),rep(2,8),rep(3,8),rep(4,8)))
machine <- as.factor(c(1,1,2,2,3,3,4,4,1,1,2,2,3,3,4,4,
                       1,1,2,2,3,3,4,4,1,1,2,2,3,3,4,4))

Machines randomly sampled. Regard days as randomly sampled.

Two-way random effects model

Assume \[ Y_{ijk} = \mu + A_i + B_j + (AB)_{ij} + \varepsilon_{ijk} \] for \(i=1,\dots,a\), \(j=1,\dots,b\), and \(n=1,\dots,n\), where

• \(A_i\overset{\operatorname{ind}}{\sim}\mathcal{N}(0,\sigma_A^2)\)
• \(B_j \overset{\operatorname{ind}}{\sim}\mathcal{N}(0,\sigma_B^2)\)
• \((AB)_{ij} \overset{\operatorname{ind}}{\sim}\mathcal{N}(0,\sigma_{AB}^2)\)
• \(\varepsilon_{ijk} \overset{\operatorname{ind}}{\sim}\mathcal{N}(0,\sigma_\varepsilon^2)\)
• the \(A_i\), \(B_j\), and \(\varepsilon_{ijk}\) are all independent.
• \(\mu\) is the overall mean.
• we call \(\sigma_A^2\), \(\sigma_B^2\), \(\sigma_{AB}^2\), and \(\sigma_\varepsilon^2\) variance components.

Still assume a balanced design, that is \(n_{ij} = n\) for all \(i,j\).