1 Rohwer data: Aptitude and achievement
To illustrate the homogeneity of regression flavor, we use data from a study by Rohwer (given in Timm (1975), Ex. 4.3, 4.7, and 4.23) on kindergarten children, designed to determine how well a set of paired-associate (PA) tasks predicted performance on measures of achievement:
- the Peabody Picture Vocabulary test (
PPVT), - a student achievement test (
SAT), and - the Raven Progressive matrices test (
Raven).
The PA tasks were considered measures of learning aptitude and
varied in how the stimuli were presented, and are called named (n), still (s), named still (ns),
named action (na), and sentence still (ss).
Two groups were tested: a group of \(n=37\) children from a low socioeconomic
status (SES) school, and a group of \(n=32\) high SES children from an
upper-class, white residential school. The data are in the data frame
Rohwer in the heplots package:
data(Rohwer)
Rohwer |> dplyr::sample_n(6)
#> group SES SAT PPVT Raven n s ns na ss
#> 31 1 Lo 9 63 11 2 12 5 25 14
#> 7 1 Lo 6 71 21 0 1 20 23 18
#> 18 1 Lo 45 54 10 0 6 6 14 16
#> 17 1 Lo 19 66 13 7 12 21 35 27
#> 50 2 Hi 4 87 14 1 4 14 25 19
#> 59 2 Hi 36 89 15 1 6 15 23 281.1 Preliminary plots
Before fitting models, it is usually useful to do some data exploration and graphing. With multivariate multiple regression data, among the most helpful plots are scatterplots of each response variable, Y, against each predictor, X, and we can get a better sense of the relationships by adding linear regression lines, loess smooths or other enhancements.
A scatterplot matrix, using graphics::pairs() or GGally::ggpairs() is easy to do.
However, with 3 response variables, 4 predictors and a group factor (SES), this can
be overwhelming. An alternative is to compose a rectangular matrix of plots for
only the Y variables against the Xs.
This turned out to be not as easy as it might seem, because none of the pairs() methods
allow for this possibility. The trick is to reshape the data from wide to long format
and use facets in ggplot2 to compose the pairwise scatterplots into the desired
rectangular matrix format.1
library(tidyr)
library(dplyr)
library(ggplot2)
yvars <- c("SAT", "PPVT", "Raven" ) # outcome variables
xvars <- c("n", "s", "ns", "na", "ss") # predictors
xvars <- c("n", "s", "ns") # make a smaller example
Rohwer_long <- Rohwer %>%
dplyr::select(-group, -na, -ss) |>
tidyr::pivot_longer(cols = all_of(xvars),
names_to = "xvar", values_to = "x") |>
tidyr::pivot_longer(cols = all_of(yvars),
names_to = "yvar", values_to = "y") |>
dplyr::mutate(xvar = factor(xvar, levels = xvars),
yvar = factor(yvar, levels = yvars))
Rohwer_long
#> # A tibble: 621 × 5
#> SES xvar x yvar y
#> <fct> <fct> <int> <fct> <int>
#> 1 Lo n 1 SAT 49
#> 2 Lo n 1 PPVT 48
#> 3 Lo n 1 Raven 8
#> 4 Lo s 2 SAT 49
#> 5 Lo s 2 PPVT 48
#> 6 Lo s 2 Raven 8
#> 7 Lo ns 6 SAT 49
#> 8 Lo ns 6 PPVT 48
#> 9 Lo ns 6 Raven 8
#> 10 Lo n 5 SAT 47
#> # ℹ 611 more rowsThen, we can use ggplot2 to make produce the pairwise plots for each combination of x and y variables. Using color=SES in the aesthetic results in a separate regression line for
the two SES groups produced by geom_smooth().
ggplot(Rohwer_long, aes(x, y, color = SES, shape = SES)) +
geom_jitter(size=1.5) +
geom_smooth(method = "lm",
se = FALSE,
formula = y ~ x,
size=1.5) +
facet_grid(yvar ~ xvar, # plot matrix of Y by X
scales = "free") +
theme_bw(base_size = 16) +
theme(legend.position = "bottom")
Figure 1.1: Plot matrix of the response variables (SAT, PPVT and Raven) against each of three predictors (n, s, ns). The lines show separate linear regressions for each SES group.
Such plots form a framework for understanding model fits and statistical tests we turn to now.
1.2 Separate models
As one approach, we might be tempted to fit separate regression models for each of the High and Low SES groups. This approach is not generally recommended because it lacks power (smaller sample sizes in each group than a combined analysis) and does not allow hypotheses about equality of slopes and intercepts to be tested directly.
rohwer.ses1 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer,
subset=SES=="Hi")
Anova(rohwer.ses1)
#>
#> Type II MANOVA Tests: Pillai test statistic
#> Df test stat approx F num Df den Df Pr(>F)
#> n 1 0.202 2.02 3 24 0.1376
#> s 1 0.310 3.59 3 24 0.0284 *
#> ns 1 0.358 4.46 3 24 0.0126 *
#> na 1 0.465 6.96 3 24 0.0016 **
#> ss 1 0.089 0.78 3 24 0.5173
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rohwer.ses2 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer,
subset=SES=="Lo")
Anova(rohwer.ses2)
#>
#> Type II MANOVA Tests: Pillai test statistic
#> Df test stat approx F num Df den Df Pr(>F)
#> n 1 0.0384 0.39 3 29 0.764
#> s 1 0.1118 1.22 3 29 0.321
#> ns 1 0.2252 2.81 3 29 0.057 .
#> na 1 0.2675 3.53 3 29 0.027 *
#> ss 1 0.1390 1.56 3 29 0.220
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This fits separate slopes and intercepts for each of the two groups, but it is difficult to compare the coefficients numerically.
coef(rohwer.ses1)
#> SAT PPVT Raven
#> (Intercept) -28.4675 39.69709 13.24384
#> n 3.2571 0.06728 0.05935
#> s 2.9966 0.36998 0.49244
#> ns -5.8591 -0.37438 -0.16402
#> na 5.6662 1.52301 0.11898
#> ss -0.6227 0.41016 -0.12116
coef(rohwer.ses2)
#> SAT PPVT Raven
#> (Intercept) 4.15106 33.00577 11.17338
#> n -0.60887 -0.08057 0.21100
#> s -0.05016 -0.72105 0.06457
#> ns -1.73240 -0.29830 0.21358
#> na 0.49456 1.47042 -0.03732
#> ss 2.24772 0.32396 -0.05214The function heplots::coefplot() makes this a bit easier, by plotting bivariate confidence ellipses for the
coefficients in a multivariate linear model. In this problem, with three response variables,
the 95% confidence regions are 3D ellipsoids, but we only plot them in 2D.
The 3D versions have the property that a given predictor is significant by a multivariate test if
the ellipsoid excludes the point (0, 0, 0).
coefplot(rohwer.ses1, fill=TRUE, cex.label=1.5, cex.lab=1.5)
text(-10, 3, "High SES group", pos=4, cex=1.4)
coefplot(rohwer.ses2, fill=TRUE, cex.label=1.5, cex.lab=1.5)
text(-4.7, 2.5, "Low SES group", pos=4, cex=1.4)
Figure 1.2: Coefficient plots for the separate models for the High and Low SES groups in the Rohwer data. The ellipses are 95% confidence regions for the pairs of regression coefficients of (SAT, PPVT) for each predictor in the model.
Alternatively, we can visualize the results of the multivariate tests for the predictors
with HE plots. Here we make use of the fact
that several HE plots can be overlaid using the option add=TRUE
as shown in Figure 1.3.
heplot(rohwer.ses1,
ylim=c(40,110), # allow more room for 2nd plot
col=c("red", "black"),
fill = TRUE, fill.alpha = 0.1,
lwd=2, cex=1.2)
heplot(rohwer.ses2,
add=TRUE,
col=c("brown", "black"),
grand.mean=TRUE, error.ellipse=TRUE, # not shown by default when add=TRUE
fill = TRUE, fill.alpha = 0.1,
lwd=2, cex=1.2)
# label the groups at their centroid
means <- aggregate(cbind(SAT,PPVT)~SES, data=Rohwer, mean)
text(means[,2], means[,3], labels=means[,1], pos=3, cex=2, col="black")
Figure 1.3: HE plot for SAT and PPVT, showing the effects for the PA predictors for the High and Low SES groups separately
We can readily see the difference in means for the two SES groups (Hi has greater scores on both variables)
and it also appears that the slopes of the predictor ellipses are shallower for the High
than the Low group, indicating greater relation with the SAT score.
1.3 MANCOVA model
Alternatively (and optimistically), we can fit a MANCOVA model that allows different means for the two SES groups on the responses, but constrains the slopes for the PA covariates to be equal.
# MANCOVA, assuming equal slopes
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 ' ' 1Note that, although the
multivariate tests for two of the covariates (ns and na)
are highly significant, univariate multiple regression tests for the
separate responses [from summary(rohwer.mod)] are relatively weak.
We can also test the global 5 df hypothesis, \(\mathbf{B}=\mathbf{0}\),
that all covariates have null effects
for all responses as a linear hypothesis.
First, extract the names
of the PA tests predictors from the model.
car::linearHypothesis() takes a vector of the names coefficients to be tested simultaneously.
Regr <- linearHypothesis(rohwer.mod, covariates)
print(Regr, digits=4, SSP=FALSE)
#>
#> Multivariate Tests:
#> Df test stat approx F num Df den Df Pr(>F)
#> Pillai 5 0.6658 3.537 15 186 2.31e-05 ***
#> Wilks 5 0.4418 3.812 15 166 8.28e-06 ***
#> Hotelling-Lawley 5 1.0309 4.032 15 176 2.79e-06 ***
#> Roy 5 0.7574 9.392 5 62 1.06e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Then 2D views of the
additive MANCOVA model rohwer.mod and the overall test for all covariates
are produced as follows, giving the plots in Figure 1.4.
colors <- c("red", "blue", rep("black",5), "#969696")
heplot(rohwer.mod,
col=colors, variables=c(1,2),
hypotheses=list("Regr" = covariates),
fill = TRUE, fill.alpha = 0.1,
cex=1.5, lwd=c(2, rep(3,5), 4),
main="(SAT, PPVT) in Rohwer MANCOVA model")
heplot(rohwer.mod,
col=colors, variables=c(1,3),
hypotheses=list("Regr" = covariates),
fill = TRUE, fill.alpha = 0.1,
cex=1.5, lwd=c(2, rep(3,5), 4),
main="(SAT, Raven) in Rohwer MANCOVA model")
Figure 1.4: HE plot for SAT and PPVT (left) and for SAT and Raven (right) using the MANCOVA model. The ellipses labeled ‘Regr’ show the test of the overall model, including all predictors.
The positive orientation of the Regr ellipses shows that the predicted
values for all three responses are positively correlated (more so for
SAT and PPVT). As well, the High SES group
is higher on all responses than the Low SES group.
Alternatively, all pairwise plots among these responses could be drawn using
the pairs.mlm() function,
pairs(rohwer.mod, col=colors,
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")),
cex=1.3, lwd=c(2, rep(3,5), 4))
or as a 3D plot, using heplot3d() as shown in Figure 1.5.
colors <- c("pink", "blue", rep("black",5), "#969696")
heplot3d(rohwer.mod, col=colors,
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))
Figure 1.5: 3D HE plot for the MANCOVA model fit to the Rohwer data.
1.4 Testing homogeneity of regression
The MANCOVA model, rohwer.mod, has relatively simple interpretations
(large effect of SES, with ns and na as the major predictors)
but the test of relies on the assumption of homogeneity of slopes for the predictors.
We can test this assumption as follows, by adding interactions of SES
with each of the covariates:
rohwer.mod2 <- lm(cbind(SAT, PPVT, Raven) ~ SES * (n + s + ns + na + ss),
data=Rohwer)
Anova(rohwer.mod2)
#>
#> 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 ' ' 1It appears from the above that there is only weak evidence of unequal slopes
from the separate SES: terms. The evidence for heterogeneity is
stronger, however, when these terms are tested collectively using the
linearHypothesis() function:
(coefs <- rownames(coef(rohwer.mod2)))
#> [1] "(Intercept)" "SESLo" "n" "s" "ns"
#> [6] "na" "ss" "SESLo:n" "SESLo:s" "SESLo:ns"
#> [11] "SESLo:na" "SESLo:ss"
print(linearHypothesis(rohwer.mod2, 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 ' ' 1This model (rohwer.mod2) is similar in spirit to the two models
(rohwer.ses1 and rohwer.ses2)
fit for the two SES groups separately and show in Figure 1.3,
except that model rohwer.mod2 assumes a common within-groups error covariance matrix
and allows overall tests.
To illustrate model rohwer.mod2, we construct an HE plot for
SAT and PPVT shown in Figure 1.6.
To simplify this display, we show the hypothesis ellipses
for the overall effects of the PA tests in the baseline high-SES group, and
a single combined ellipse for all the SESLo: interaction terms that
we tested previously, representing differences in slopes between the low and
high-SES groups.
Because SES is “treatment-coded” in this model, the ellipse for each
covariate represents the hypothesis that the slopes for that covariate are
zero in the high-SES baseline category. With this parameterization, the ellipse for
Slopes represents the joint hypothesis that slopes for the covariates do not differ
in the low-SES group.
colors <- c("red", "blue", rep("black",5), "#969696")
heplot(rohwer.mod2, col=c(colors, "brown"),
terms=c("SES", "n", "s", "ns", "na", "ss"),
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss"),
"Slopes" = coefs[grep(":", coefs)]))
Figure 1.6: HE plot for SAT and PPVT, fitting the model rohwer.mod2 that allows unequal slopes for the covariates.
Comparing Figure 1.6 for the heterogeneous slopes model
with Figure 1.4 (left) for the homogeneous slopes model,
it can be seen that most of the covariates have ellipses of similar
size and orientation, reflecting similar evidence against the respective
null hypotheses, as does the effect of SES, showing the
greater performance of the high-SES group on all response measures.
Somewhat more subtle, the error ellipse is noticeably smaller in
Figure 1.6, reflecting the additional variation accounted
for by differences in slopes.
