This script reproduces all of the analysis and graphs for the MANCOVA of the Rohwer data in the paper and also includes other analyses not described there. It is set up as an R script that can be “compiled” to HTML, Word, or PDF using knitr::knit(). This is most convenient within R Studio via the File -> Compile Notebook option.

Load packages and the data

library(heplots)
data(Rohwer, package="heplots")

Fit the MANCOVA model, test hypotheses:

SES is the group factor, representing Hi and Lo socioeconomic status

rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ SES + n + s + ns + na + ss,
                 data=Rohwer)
Anova(rohwer.mod)
## 
## Type II MANOVA Tests: Pillai test statistic
##     Df test stat approx F num Df den Df  Pr(>F)    
## SES  1     0.379    12.18      3     60 2.5e-06 ***
## n    1     0.040     0.84      3     60  0.4773    
## s    1     0.093     2.04      3     60  0.1173    
## ns   1     0.193     4.78      3     60  0.0047 ** 
## na   1     0.231     6.02      3     60  0.0012 ** 
## ss   1     0.050     1.05      3     60  0.3770    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pairwise HE plots

pairs(rohwer.mod, 
      hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")), 
      fill=TRUE, fill.alpha=0.1)

Fit heterogeneous regression model with SES interactions:

rohwer.mod1 <- lm(cbind(SAT, PPVT, Raven) ~ SES * (n + s + ns + na + ss), 
                  data=Rohwer)
Anova(rohwer.mod1)
## 
## Type II MANOVA Tests: Pillai test statistic
##        Df test stat approx F num Df den Df  Pr(>F)    
## SES     1     0.391    11.78      3     55 4.5e-06 ***
## n       1     0.079     1.57      3     55 0.20638    
## s       1     0.125     2.62      3     55 0.05952 .  
## ns      1     0.254     6.25      3     55 0.00100 ***
## na      1     0.307     8.11      3     55 0.00015 ***
## ss      1     0.060     1.17      3     55 0.32813    
## SES:n   1     0.072     1.43      3     55 0.24417    
## SES:s   1     0.099     2.02      3     55 0.12117    
## SES:ns  1     0.118     2.44      3     55 0.07383 .  
## SES:na  1     0.148     3.18      3     55 0.03081 *  
## SES:ss  1     0.057     1.12      3     55 0.35094    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Assess significance of interactions

linearHypothesis() takes a vector of terms to test simulataneously. Here we use grep() to find all interaction terms, which have a “:” in their name.

coefs <- rownames(coef(rohwer.mod1))   # store coefficient names in a vector
print(linearHypothesis(rohwer.mod1,    # only test for interaction effects
                       coefs[grep(":", coefs)]), SSP=FALSE) 
## 
## Multivariate Tests: 
##                  Df test stat approx F num Df den Df  Pr(>F)   
## Pillai            5    0.4179    1.845     15  171.0 0.03209 * 
## Wilks             5    0.6236    1.894     15  152.2 0.02769 * 
## Hotelling-Lawley  5    0.5387    1.927     15  161.0 0.02396 * 
## Roy               5    0.3846    4.385      5   57.0 0.00191 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Alternatively, fit models for each group separately:

This also allows us to see possible differences in within-group error variances and covariances

rohwer.ses1 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, 
                  data = Rohwer, subset = SES == "Hi")
rohwer.ses2 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, 
                  data = Rohwer, subset = SES == "Lo")

Visualize results with overlaid HE plots:

For this to work correctly, it is necessary to set the x,y axis limits the explicitly to be the same in both plots

## Low SES students:
heplot(rohwer.ses2, col = c("red", rep("black",5), "blue"), 
       hypotheses = list("B=0, Low SES" = c("n", "s", "ns", "na", "ss")), 
       level = 0.5, cex = 1.25, 
       fill = c(TRUE, FALSE), fill.alpha = 0.05, 
       xlim = c(-15, 110), ylim = c(40, 110),
       xlab = "Student Achievement Test", 
       ylab = "Peabody Picture Vocabulary Test", 
       label.pos = c(1, rep(NULL, 5), 1))

## High SES students:
heplot(rohwer.ses1, col = c("red", rep("black", 5), "blue"), 
       hypotheses = list("B=0, High SES" = c("n", "s", "ns", "na", "ss")), 
       level = 0.5, cex = 1.25, 
       add = TRUE,   # place both plots on same graphic
       error = TRUE, # error ellipse is not drawn by default with add = TRUE
       fill = c(TRUE, FALSE), fill.alpha = 0.05, 
       xlim = c(-15, 110), ylim = c(40, 110))