# Nested Designs

## PVPS Example

> ### Data entry

> pvps <- data.frame(score=c(12,13,11,12,15,19,17,17,10,12,11,17,

+   18,20,21,16,20,14,17,15,19,22,21,21),

+   school=gl(2,12), class=gl(6,4))

> save.image()

> ### Fixed-Effect Model

> ## Homogeneous Variance Tests

> car::leveneTest(score~class, data=pvps)

```Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group  5  0.9493 0.4737
18
```

> bartlett.test(score~class, data=pvps)

```Bartlett test of homogeneity of variances

data:  score by class
Bartlett's K-squared = 5.4979, df = 5, p-value = 0.3582
```

> fligner.test(score~class, data=pvps)

```Fligner-Killeen test of homogeneity of variances

data:  score by class
Fligner-Killeen:med chi-squared = 4.9199, df = 5, p-value =
0.4257
```

> # All three test conclude consistently that the variances

> # appear to be homogeneous

> ## Fit AOV

> pvps.aov <- aov(score ~ school + class %in% school, data=pvps)

> summary(pvps.aov)

```             Df  Sum Sq Mean Sq F value    Pr(>F)
school        1 140.167 140.167 31.7358 2.408e-05 ***
school:class  4  96.833  24.208  5.4811  0.004574 **
Residuals    18  79.500   4.417
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

> # Note that coding classes of school 2 as 1, 2, and 3 will

> # give the same results as demonstrated below. However, this

> # will give correct Tukey's all pairwise comparisons.

> pvps\$classes <- gl(3,4,24)

> pvps.aov2 <- aov(score ~ school + classes %in% school, data=pvps)

> summary(pvps.aov2)

```               Df  Sum Sq Mean Sq F value    Pr(>F)
school          1 140.167 140.167 31.7358 2.408e-05 ***
school:classes  4  96.833  24.208  5.4811  0.004574 **
Residuals      18  79.500   4.417
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

> pvps.tk <- TukeyHSD(pvps.aov, "school:classes")

> pvps.tk

```  Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = score ~ school + classes %in% school, data = pvps)

\$`school:classes`
2:1-1:1  6.75   2.0272933 11.4727067 0.0029272
1:2-1:1  5.00   0.2772933  9.7227067 0.0344276
2:2-1:1  4.50  -0.2227067  9.2227067 0.0670139
1:3-1:1  0.50  -4.2227067  5.2227067 0.9993228
2:3-1:1  8.75   4.0272933 13.4727067 0.0001776
1:2-2:1 -1.75  -6.4727067  2.9727067 0.8415845
2:2-2:1 -2.25  -6.9727067  2.4727067 0.6602357
1:3-2:1 -6.25 -10.9727067 -1.5272933 0.0059730
2:3-2:1  2.00  -2.7227067  6.7227067 0.7567426
2:2-1:2 -0.50  -5.2227067  4.2227067 0.9993228
1:3-1:2 -4.50  -9.2227067  0.2227067 0.0670139
2:3-1:2  3.75  -0.9727067  8.4727067 0.1688999
1:3-2:2 -4.00  -8.7227067  0.7227067 0.1256704
2:3-2:2  4.25  -0.4727067  8.9727067 0.0922817
2:3-1:3  8.25   3.5272933 12.9727067 0.0003530
```

> # The session below is used to validate the first c.i. above.

> MSE <- sum(resid(pvps.aov2)^2)/pvps.aov\$df.resid

> MSE

```[1] 4.416667
```

> 6.75+c(-1,1)*qtukey(.05,6,pvps.aov2\$df.resid,low=F)*sqrt(MSE/4)

```[1]  2.027293 11.472707
```

> # Now produce graphical representation of Tukey's comparisons.

> windows(width=4,height=4,pointsize=10)

> plot(pvps.tk, las=1)

> title(main="PVPS Example", cex.main=1, line=1)

> dev.off()

```null device
1
```

>

> ### Random-Effect Model

> require(lme4)

```Loading required package: lme4
```

> pvps.lmer <- lmer(score ~ 1 + (1|school) + (1|school:class),

+   data=pvps)

> pvps.lmer

```Linear mixed-effects model fit by REML ['lmerMod']
Formula: score ~ 1 + (1 | school) + (1 | school:class)
Data: pvps
REML criterion at convergence: 112.8758
Random effects:
Groups       Name        Std.Dev.
school:class (Intercept) 2.224
school       (Intercept) 3.109
Residual                 2.102
Number of obs: 24, groups: school:class, 6; school, 2
Fixed Effects:
(Intercept)
16.25
```

> # Inspect the variance components!

>

> ### Mixed-Effect Model

> pvps.lmer.2 <- lmer(score ~ school + (1|school:class), data=pvps)

> pvps.lmer.2

```Linear mixed model fit by REML ['lmerMod']
Formula: score ~ school + (1 | school:class)
Data: pvps
REML criterion at convergence: 106.8868
Random effects:
Groups       Name        Std.Dev.
school:class (Intercept) 2.224
Residual                 2.102
Number of obs: 24, groups: school:class, 6
Fixed Effects:
(Intercept)      school2
13.833        4.833
```

> # Again, inspect the variance components.

> summary(pvps.lmer.2)

```Linear mixed model fit by REML ['lmerMod']
Formula: score ~ school + (1 | school:class)
Data: pvps

REML criterion at convergence: 106.8868

Random effects:
Groups       Name        Variance Std.Dev.
school:class (Intercept) 4.948    2.224
Residual                 4.417    2.102
Number of obs: 24, groups: school:class, 6

Fixed effects:
Estimate Std. Error t value
(Intercept)   13.833      1.420   9.740
school2        4.833      2.009   2.406

Correlation of Fixed Effects:
(Intr)
school2 -0.707
```