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`
         diff         lwr        upr     p adj
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

(PVPS11.png here)

>

> ### Random-Effect Model

> require(lme4)

Loading required package: lme4
Loading required package: lattice
Loading required package: Matrix

> 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