Nothing
HE plot MMRA Examples"
In heplots: Visualizing Hypothesis Tests in Multivariate Linear Models
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE,
fig.height=5,
fig.width=5,
fig.align = "center",
# results='hide',
# fig.keep='none',
fig.path='fig/mmra-',
echo=TRUE,
collapse = TRUE,
comment = "#>"
)
set.seed(1071)
options(width=80, digits=4, continue=" ")
library(heplots)
library(candisc)
library(car)
library(dplyr)
library(tidyr)
library(ggplot2)
Vignette built using heplots, version r packageDescription("heplots")[["Version"]]
and candisc, version r packageDescription("candisc")[["Version"]].
Multivariate Multiple Regression Designs {-}
The ideas behind HE plots extend naturally to multivariate multiple
regression (MMRA) and multivariate analysis of covariance (MANCOVA).
In MMRA designs, the $\mathbf{X}$ matrix contains only quantitative predictors, while
in MANCOVA designs, it contains a mixture of
factors and quantitative predictors (covariates), but typically there is just
one "group" factor.
In the MANCOVA case,
there is often a subtle difference in emphasis: true MANCOVA analyses focus on
the differences among groups defined by the factors, adjusting for
(or controlling for) the quantitative covariates. Analyses concerned
with homogeneity of regression focus on quantitative predictors
and attempt to test whether the regression relations are the same for
all groups defined by the factors.
Rohwer data: Aptitude and achievement
To illustrate the homogeneity of regression flavor,
we use data
from a study by Rohwer (given in @Timm:75, 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)
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]
[^1]: This solution was suggested in an answer to a Stackoverflow question, https://stackoverflow.com/questions/73859139/how-to-make-a-scatterplot-rectangular-matrix-y1-y2-x1-x2-in-r
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
Then, 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().
#| fig.width=8,
#| fig.heigth=10,
#| fig.cap="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."
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")
Such plots form a framework for understanding model fits and statistical tests we turn to now.
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)
rohwer.ses2 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer,
subset=SES=="Lo")
Anova(rohwer.ses2)
This fits separate slopes and intercepts for each of the two groups, but it is difficult
to compare the coefficients numerically.
coef(rohwer.ses1)
coef(rohwer.ses2)
The 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).
#| echo=-1,
#| out.width="47%",
#| out.height="50%",
#| fig.show="hold",
#| collapse=FALSE,
#| fig.cap="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."
par(mar=c(4,4,1,1)+.1)
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)
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 \@ref(fig:rohwer-HE1).
#| echo=-1,
#| out.width="60%",
#| fig.cap="HE plot for `SAT` and `PPVT`, showing the effects for the PA predictors for the High and Low SES groups separately"
par(mar=c(4,4,1,1)+.1)
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")
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.
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)
Note 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.
covariates <- c("n", "s", "ns", "na", "ss")
# or: covariates <- rownames(coef(rohwer.mod))[-(1:2)]
Regr <- linearHypothesis(rohwer.mod, covariates)
print(Regr, digits=4, SSP=FALSE)
Then 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 \@ref(fig:rohwer-HE2).
#| echo=-1,
#| out.width="47%",
#| out.height="50%",
#| fig.show="hold",
#| collapse=FALSE,
#| fig.cap="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."
par(mar=c(4,4,3,1)+.1)
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")
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,
#| fig.width=6,
#| fig.height=6
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 \@ref(fig:rohwer-HE3D).
colors <- c("pink", "blue", rep("black",5), "#969696")
heplot3d(rohwer.mod, col=colors,
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))
#| echo=FALSE,
#| out.width="60%",
#| fig.cap="3D HE plot for the MANCOVA model fit to the Rohwer data."
knitr::include_graphics("fig/mmra-rohwer-HE3D.png")
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)
It 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)))
print(linearHypothesis(rohwer.mod2, coefs[grep(":", coefs)]), SSP=FALSE)
This 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 \@ref(fig:rohwer-HE1),
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 \@ref(fig:rohwer-HE4).
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.
#| echo=-1,
#| out.width="70%",
#| fig.cap="HE plot for `SAT` and `PPVT`, fitting the model `rohwer.mod2` that allows unequal slopes for the covariates."
par(mar=c(4,4,1,1)+.1)
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)]))
Comparing Figure \@ref(fig:rohwer-HE4) for the heterogeneous slopes model
with Figure \@ref(fig:rohwer-HE2) (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 \@ref(fig:rohwer-HE4), reflecting the additional variation accounted
for by differences in slopes.
Recovery from hernia repair
This example uses the Hernior data (from @MostellerTukey:1977, Data Exhibit 8, p. 567-568),
comprising
data on measures of post-operative recovery of 32 patients undergoing an
elective herniorrhaphy operation, in relation to pre-operative measures.
The outcome measures are:
leave, the patient's condition upon leaving the recovery room (a 1-4 scale, 1=best),nurse, level of nursing required one week after operation (a 1-5 scale, 1=worst) and los, length of stay in hospital after operation (in days)
The predictor variables are:
- patient
age, sex,
- `pstat, physical status ( a 1-5 scale, with 1=perfect health, ..., 5=very poor health),
build, body build (a 1-5 scale, with 1=emaciated, ..., 5=obese), and- preoperative complications with (
cardiac) heart and respiration (resp), 1-4 scales, 1=none, ..., 4=severe.
data(Hernior)
str(Hernior)
All predictors model
We begin with a model fitting all predictors. Note that the ordinal predictors,
pstat, build, cardiac and resp could arguably be treated as
factors, rather than linear, regression terms. Doing so would give tests for
nonlinear effects of their relations with the responses.
We ignore this possibility in this example.
Hern.mod <- lm(cbind(leave, nurse, los) ~ age + sex + pstat + build + cardiac + resp,
data=Hernior)
Anova(Hern.mod)
The results of the multivariate tests above are somewhat disappointing. Only the physical status
predictor (pstat) appears to be significant at conventional levels.
The univariate models for each response are implicit in the MLM Hern.mod.
These can be printed using summary(), or we can use summary() to extract
certain statistics for each univariate response model, as we do here.
Hern.summary <- summary(Hern.mod)
unlist(lapply(Hern.summary, function(x) x$r.squared))
More conveniently, the function heplots::glance.mlm() extends broom::glance.lm()
to give a one-line summary of statistics for each response variable in a MLM.
The $R^2$ and $F$ statistics are those for each overall model assessing the impact
of all predictors.
glance.mlm(Hern.mod)
Univariate tests for predictors in each of these models (not shown here)
are hard to interpret, and largely show only significant effects for
the leave variable. Yet, the $R^2$ values for the other responses
are slightly promising. We proceed to a multivariate overall test of $\mathbf{B} = 0$
for all predictors, whose term names can be easily extracted from the
rownames of the coefficients.
# test overall regression
(predictors <- rownames(coef(Hern.mod))[-1])
Regr <- linearHypothesis(Hern.mod, predictors)
print(Regr, digits=5, SSP=FALSE)
#| fig.width=7,
#| fig.height=7,
#| fig.cap="HE pairs plot for Hernior data"
clr <- c("red", "darkgray", "blue", "darkgreen", "magenta", "brown", "black")
vlab <- c("LeaveCondition\n(leave)",
"NursingCare\n(nurse)",
"LengthOfStay\n(los)")
hyp <- list("Regr" = predictors)
pairs(Hern.mod,
hypotheses=hyp,
col=clr, var.labels=vlab,
fill=c(TRUE,FALSE), fill.alpha = 0.1,
cex=1.25)
A pairs() plot for the MLM gives the set of plots shown in Figure \@ref(fig:hern-pairs)
helps to interpret the relations among the predictors which lead to the highly significant
overall test.
Among the predictors, age and sex have small and insignificant effects on the outcome measures
jointly. The other predictors all contribute to the overall test of $\mathbf{B} = 0$,
though in different ways for the various responses.
For example, in the panel for (leave, los) in Figure \@ref(fig:hern-pairs),
it can be seen that while only pstat individually is outside the
$\mathbf{E}$ ellipse, build and resp contribute to the overall test in
an opposite direction.
Canonical analysis
In this multivariate regression example, all of the terms in the model Hern.mod
have 1 df, and so plot as lines in HE plots. An alternative view of these effects
can be seen in canonical discriminant space, which, for each predictor shows the
scores on the linear combination of the responses that contributes most to
the multivariate test of that effect, together with the weights for the responses.
We use candiscList() to calculate the canonical analyses for all predictor terms in
Hern.mod.
Hern.canL <- candiscList(Hern.mod)
1D canonical discriminant plots for all terms can be obtained interactively
with a menu, simply by plotting the Hern.canL object.
plot(Hern.canL)
Plots for separate terms are produced by the lines below, and shown in
Figure \@ref(fig:hern-can1) and Figure \@ref(fig:hern-can2).
For pstat and build:
#| out.width="47%",
#| out.height="50%",
#| fig.show="hold",
#| collapse=FALSE,
#| fig.cap="1D Canonical discriminant plots for physical status (`pstat`) and `build`. The canonical scores are such that better outcomes are associated with smaller scores. Arrows show the correlations of the responses with the 1D canonical scores."
plot(Hern.canL, term="pstat")
plot(Hern.canL, term="build")
For age and cardiac:
#| out.width="47%",
#| out.height="50%",
#| fig.show="hold",
#| collapse=FALSE,
#| fig.cap="1D Canonical discriminant plots for `age` and `cardiac`. The canonical scores are such that better outcomes are associated with smaller scores."
plot(Hern.canL, term="age")
plot(Hern.canL, term="cardiac")
In these plots, the canonical scores panel shows
the linear combinations of the
response variables which have the largest possible $R^2$.
Better outcomes correspond to numerically smaller canonical scores.
The arrows in the structure panel are proportional to
the canonical weights.
These plots provide simple interpretations of the results for the
canonical combinations of the responses.
Better physical status, smaller body build, lower age and absence
of cardiac complications are all positively related to better outcomes.
Grades in a Sociology Course
The data set SocGrades contains four outcome measures on student performance
in an introductory sociology course together with six potential predictors.
These data were used by @MarascuiloLevin:1983 for an example of
canonical correlation analysis, but are also suitable as examples of
multivariate multiple regression, MANOVA, MANCOVA
and step-down analysis in multivariate linear models.
The outcome measures used here are three test scores during the course,
midterm1, midterm2, final,
and a course evaluation (eval).[^2]
[^2]:It is arguable that the students' course evaluation should not
be considered a response variable here. It could be used as a
predictor in a follow-up, step-down analysis, which would address the separate
question of whether the effects on exam grades remain, when eval is controlled for.
Predictor variables are:
class, the student's social class (an ordered factor with levels 1 > 2 > 3)sex,gpa, grade point average,boards, College Board test scores,hssoc, previous high school unit in sociology? (with values no, yes), andpretest, score on a course pretest.
str(SocGrades)
Models
The basic MLM is fit below as grades.mod with all predictor variables.
data(SocGrades)
grades.mod <- lm(cbind(midterm1, midterm2, final, eval) ~
class + sex + gpa + boards + hssoc + pretest,
data=SocGrades)
Anova(grades.mod, test="Roy")
Screening for interactions
In both univariate and multivariate
response models, it is often useful to screen for higher-order
terms (interactions, non-linear predictors). This can most easily be
done using update(), as we do below. First, try the extended
model with all pairwise interactions of the predictors.
In the model formula, ~.^2, the . represents all terms in the model,
and the ^2 generates all products of terms, such as class:sex, class:gpa, and so forth.
grades.mod2 <- update(grades.mod, . ~ .^2)
Anova(grades.mod2, test="Roy")
In the results above, only the interaction of class:sex is significant,
and the main effects of hssoc and pretest remain insignificant.
A revised model to explore is grades.mod3,
grades.mod3 <- update(grades.mod, . ~ . + class:sex - hssoc - pretest)
Anova(grades.mod3, test="Roy")
A pairwise HE plot for all responses (Figure \@ref(fig:grades-pairs) shows that nearly all
effects are in the expected directions: higher gpa, boards, class
leads to better performance on most outcomes. The interaction of
class:sex seems to be confined largely to midterm1.
#| out.width="90%",
#| fig.cap="HE pairs plot for SocGrade, model `grades.mod3`."
pairs(grades.mod3)
These effects are easier to appreciate for the three exam grades jointly in a
3D HE plot when you can rotate it interactively. A snapshot is shown in Figure \@ref(fig:grades-HE3D).
heplot3d(grades.mod3, wire=FALSE)
#| echo=FALSE,
#| out.width="70%",
#| fig.cap="3D HE plot for SocGrades, model `grades.mod3`"
knitr::include_graphics("fig/grades-HE3D.png")
Interactive
rotation of this plot shows that the effect of class is only two dimensional,
and of these, one dimension is very small. The major axis of the class ellipsoid
is aligned with increasing performance on all three grades, with the expected
ordering of the three social classes.
Canonical analysis
The representation of these effects in canonical space is particularly useful here.
Again, use candiscList() to compute the canonical decompositions for all terms
in the model, and extract the canonical $R^2$ from the terms in the result.
# calculate canonical results for all terms
grades.can <- candiscList(grades.mod3)
# extract canonical R^2s
unlist(lapply(grades.can, function(x) x$canrsq))
We use heplot() on the "candiscList" object to show
the effects of class in canonical space, giving Figure \@ref(fig:grades-can-class).
#| echo=-1,
#| out.width="80%",
#| fig.cap="Canonical HE plot for `class` effect in `grades.mod3`"
par(xpd=TRUE, mar=c(4,4,1,1)+.1)
# plot class effect in canonical space
heplot(grades.can, term="class",
scale=4, fill=TRUE, var.col="black", var.lwd=2)
It can be seen in Figure \@ref(fig:grades-can-class) that nearly all variation
in exam performance due to class is aligned with the first canonical dimension.
The three tests and course evaluation all have similar weights on this dimension,
but the course evaluation differs from the rest along a second, very small
dimension.
1D plots of the canonical scores for other effects in the model are also of
interest, and provide simple interpretations of these effects on the response
variables. The statements below produce the plots shown in Figure \@ref(fig:grades-can-all).
#| out.width="47%",
#| out.height="50%",
#| fig.show="hold",
#| collapse=FALSE,
#| fig.cap="1D Canonical discriminant plots for `sex` and `gpa`. Higher canonical scores reflect better course performance."
plot(grades.can, term="sex")
plot(grades.can, term="gpa")
It is readily seen that males perform better overall, but the effect of
sex is strongest for the midterm2.
As well, increasing course performance on tests is strongly associated with
gpa.
References
Try the heplots package in your browser
Any scripts or data that you put into this service are public.
heplots documentation built on Sept. 8, 2023, 5:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( message = FALSE, warning = FALSE, fig.height=5, fig.width=5, fig.align = "center", # results='hide', # fig.keep='none', fig.path='fig/mmra-', echo=TRUE, collapse = TRUE, comment = "#>" )
set.seed(1071) options(width=80, digits=4, continue=" ") library(heplots) library(candisc) library(car) library(dplyr) library(tidyr) library(ggplot2)
Vignette built using heplots, version r packageDescription("heplots")[["Version"]]
and candisc, version r packageDescription("candisc")[["Version"]].
Multivariate Multiple Regression Designs {-}
The ideas behind HE plots extend naturally to multivariate multiple regression (MMRA) and multivariate analysis of covariance (MANCOVA). In MMRA designs, the $\mathbf{X}$ matrix contains only quantitative predictors, while in MANCOVA designs, it contains a mixture of factors and quantitative predictors (covariates), but typically there is just one "group" factor.
In the MANCOVA case, there is often a subtle difference in emphasis: true MANCOVA analyses focus on the differences among groups defined by the factors, adjusting for (or controlling for) the quantitative covariates. Analyses concerned with homogeneity of regression focus on quantitative predictors and attempt to test whether the regression relations are the same for all groups defined by the factors.
Rohwer data: Aptitude and achievement
To illustrate the homogeneity of regression flavor, we use data from a study by Rohwer (given in @Timm:75, 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)
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]
[^1]: This solution was suggested in an answer to a Stackoverflow question, https://stackoverflow.com/questions/73859139/how-to-make-a-scatterplot-rectangular-matrix-y1-y2-x1-x2-in-r
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
Then, 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().
#| fig.width=8, #| fig.heigth=10, #| fig.cap="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." 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")
Such plots form a framework for understanding model fits and statistical tests we turn to now.
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) rohwer.ses2 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer, subset=SES=="Lo") Anova(rohwer.ses2)
This fits separate slopes and intercepts for each of the two groups, but it is difficult to compare the coefficients numerically.
coef(rohwer.ses1) coef(rohwer.ses2)
The 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).
#| echo=-1, #| out.width="47%", #| out.height="50%", #| fig.show="hold", #| collapse=FALSE, #| fig.cap="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." par(mar=c(4,4,1,1)+.1) 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)
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 \@ref(fig:rohwer-HE1).
#| echo=-1, #| out.width="60%", #| fig.cap="HE plot for `SAT` and `PPVT`, showing the effects for the PA predictors for the High and Low SES groups separately" par(mar=c(4,4,1,1)+.1) 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")
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.
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)
Note 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.
covariates <- c("n", "s", "ns", "na", "ss") # or: covariates <- rownames(coef(rohwer.mod))[-(1:2)]
Regr <- linearHypothesis(rohwer.mod, covariates) print(Regr, digits=4, SSP=FALSE)
Then 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 \@ref(fig:rohwer-HE2).
#| echo=-1, #| out.width="47%", #| out.height="50%", #| fig.show="hold", #| collapse=FALSE, #| fig.cap="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." par(mar=c(4,4,3,1)+.1) 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")
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,
#| fig.width=6, #| fig.height=6 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 \@ref(fig:rohwer-HE3D).
colors <- c("pink", "blue", rep("black",5), "#969696") heplot3d(rohwer.mod, col=colors, hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))
#| echo=FALSE, #| out.width="60%", #| fig.cap="3D HE plot for the MANCOVA model fit to the Rohwer data." knitr::include_graphics("fig/mmra-rohwer-HE3D.png")
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)
It 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))) print(linearHypothesis(rohwer.mod2, coefs[grep(":", coefs)]), SSP=FALSE)
This 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 \@ref(fig:rohwer-HE1),
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 \@ref(fig:rohwer-HE4).
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.
#| echo=-1, #| out.width="70%", #| fig.cap="HE plot for `SAT` and `PPVT`, fitting the model `rohwer.mod2` that allows unequal slopes for the covariates." par(mar=c(4,4,1,1)+.1) 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)]))
Comparing Figure \@ref(fig:rohwer-HE4) for the heterogeneous slopes model
with Figure \@ref(fig:rohwer-HE2) (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 \@ref(fig:rohwer-HE4), reflecting the additional variation accounted
for by differences in slopes.
Recovery from hernia repair
This example uses the Hernior data (from @MostellerTukey:1977, Data Exhibit 8, p. 567-568),
comprising
data on measures of post-operative recovery of 32 patients undergoing an
elective herniorrhaphy operation, in relation to pre-operative measures.
The outcome measures are:
leave, the patient's condition upon leaving the recovery room (a 1-4 scale, 1=best),nurse, level of nursing required one week after operation (a 1-5 scale, 1=worst) andlos, length of stay in hospital after operation (in days)
The predictor variables are:
- patient
age,sex, - `pstat, physical status ( a 1-5 scale, with 1=perfect health, ..., 5=very poor health),
build, body build (a 1-5 scale, with 1=emaciated, ..., 5=obese), and- preoperative complications with (
cardiac) heart and respiration (resp), 1-4 scales, 1=none, ..., 4=severe.
data(Hernior) str(Hernior)
All predictors model
We begin with a model fitting all predictors. Note that the ordinal predictors,
pstat, build, cardiac and resp could arguably be treated as
factors, rather than linear, regression terms. Doing so would give tests for
nonlinear effects of their relations with the responses.
We ignore this possibility in this example.
Hern.mod <- lm(cbind(leave, nurse, los) ~ age + sex + pstat + build + cardiac + resp, data=Hernior) Anova(Hern.mod)
The results of the multivariate tests above are somewhat disappointing. Only the physical status
predictor (pstat) appears to be significant at conventional levels.
The univariate models for each response are implicit in the MLM Hern.mod.
These can be printed using summary(), or we can use summary() to extract
certain statistics for each univariate response model, as we do here.
Hern.summary <- summary(Hern.mod) unlist(lapply(Hern.summary, function(x) x$r.squared))
More conveniently, the function heplots::glance.mlm() extends broom::glance.lm()
to give a one-line summary of statistics for each response variable in a MLM.
The $R^2$ and $F$ statistics are those for each overall model assessing the impact
of all predictors.
glance.mlm(Hern.mod)
Univariate tests for predictors in each of these models (not shown here)
are hard to interpret, and largely show only significant effects for
the leave variable. Yet, the $R^2$ values for the other responses
are slightly promising. We proceed to a multivariate overall test of $\mathbf{B} = 0$
for all predictors, whose term names can be easily extracted from the
rownames of the coefficients.
# test overall regression (predictors <- rownames(coef(Hern.mod))[-1]) Regr <- linearHypothesis(Hern.mod, predictors) print(Regr, digits=5, SSP=FALSE)
#| fig.width=7, #| fig.height=7, #| fig.cap="HE pairs plot for Hernior data" clr <- c("red", "darkgray", "blue", "darkgreen", "magenta", "brown", "black") vlab <- c("LeaveCondition\n(leave)", "NursingCare\n(nurse)", "LengthOfStay\n(los)") hyp <- list("Regr" = predictors) pairs(Hern.mod, hypotheses=hyp, col=clr, var.labels=vlab, fill=c(TRUE,FALSE), fill.alpha = 0.1, cex=1.25)
A pairs() plot for the MLM gives the set of plots shown in Figure \@ref(fig:hern-pairs)
helps to interpret the relations among the predictors which lead to the highly significant
overall test.
Among the predictors, age and sex have small and insignificant effects on the outcome measures
jointly. The other predictors all contribute to the overall test of $\mathbf{B} = 0$,
though in different ways for the various responses.
For example, in the panel for (leave, los) in Figure \@ref(fig:hern-pairs),
it can be seen that while only pstat individually is outside the
$\mathbf{E}$ ellipse, build and resp contribute to the overall test in
an opposite direction.
Canonical analysis
In this multivariate regression example, all of the terms in the model Hern.mod
have 1 df, and so plot as lines in HE plots. An alternative view of these effects
can be seen in canonical discriminant space, which, for each predictor shows the
scores on the linear combination of the responses that contributes most to
the multivariate test of that effect, together with the weights for the responses.
We use candiscList() to calculate the canonical analyses for all predictor terms in
Hern.mod.
Hern.canL <- candiscList(Hern.mod)
1D canonical discriminant plots for all terms can be obtained interactively
with a menu, simply by plotting the Hern.canL object.
plot(Hern.canL)
Plots for separate terms are produced by the lines below, and shown in Figure \@ref(fig:hern-can1) and Figure \@ref(fig:hern-can2).
For pstat and build:
#| out.width="47%", #| out.height="50%", #| fig.show="hold", #| collapse=FALSE, #| fig.cap="1D Canonical discriminant plots for physical status (`pstat`) and `build`. The canonical scores are such that better outcomes are associated with smaller scores. Arrows show the correlations of the responses with the 1D canonical scores." plot(Hern.canL, term="pstat") plot(Hern.canL, term="build")
For age and cardiac:
#| out.width="47%", #| out.height="50%", #| fig.show="hold", #| collapse=FALSE, #| fig.cap="1D Canonical discriminant plots for `age` and `cardiac`. The canonical scores are such that better outcomes are associated with smaller scores." plot(Hern.canL, term="age") plot(Hern.canL, term="cardiac")
In these plots, the canonical scores panel shows the linear combinations of the response variables which have the largest possible $R^2$. Better outcomes correspond to numerically smaller canonical scores. The arrows in the structure panel are proportional to the canonical weights.
These plots provide simple interpretations of the results for the canonical combinations of the responses. Better physical status, smaller body build, lower age and absence of cardiac complications are all positively related to better outcomes.
Grades in a Sociology Course
The data set SocGrades contains four outcome measures on student performance
in an introductory sociology course together with six potential predictors.
These data were used by @MarascuiloLevin:1983 for an example of
canonical correlation analysis, but are also suitable as examples of
multivariate multiple regression, MANOVA, MANCOVA
and step-down analysis in multivariate linear models.
The outcome measures used here are three test scores during the course,
midterm1, midterm2, final,
and a course evaluation (eval).[^2]
[^2]:It is arguable that the students' course evaluation should not
be considered a response variable here. It could be used as a
predictor in a follow-up, step-down analysis, which would address the separate
question of whether the effects on exam grades remain, when eval is controlled for.
Predictor variables are:
class, the student's social class (an ordered factor with levels1>2>3)sex,gpa, grade point average,boards, College Board test scores,hssoc, previous high school unit in sociology? (with valuesno,yes), andpretest, score on a course pretest.
str(SocGrades)
Models
The basic MLM is fit below as grades.mod with all predictor variables.
data(SocGrades) grades.mod <- lm(cbind(midterm1, midterm2, final, eval) ~ class + sex + gpa + boards + hssoc + pretest, data=SocGrades) Anova(grades.mod, test="Roy")
Screening for interactions
In both univariate and multivariate
response models, it is often useful to screen for higher-order
terms (interactions, non-linear predictors). This can most easily be
done using update(), as we do below. First, try the extended
model with all pairwise interactions of the predictors.
In the model formula, ~.^2, the . represents all terms in the model,
and the ^2 generates all products of terms, such as class:sex, class:gpa, and so forth.
grades.mod2 <- update(grades.mod, . ~ .^2) Anova(grades.mod2, test="Roy")
In the results above, only the interaction of class:sex is significant,
and the main effects of hssoc and pretest remain insignificant.
A revised model to explore is grades.mod3,
grades.mod3 <- update(grades.mod, . ~ . + class:sex - hssoc - pretest) Anova(grades.mod3, test="Roy")
A pairwise HE plot for all responses (Figure \@ref(fig:grades-pairs) shows that nearly all
effects are in the expected directions: higher gpa, boards, class
leads to better performance on most outcomes. The interaction of
class:sex seems to be confined largely to midterm1.
#| out.width="90%", #| fig.cap="HE pairs plot for SocGrade, model `grades.mod3`." pairs(grades.mod3)
These effects are easier to appreciate for the three exam grades jointly in a 3D HE plot when you can rotate it interactively. A snapshot is shown in Figure \@ref(fig:grades-HE3D).
heplot3d(grades.mod3, wire=FALSE)
#| echo=FALSE, #| out.width="70%", #| fig.cap="3D HE plot for SocGrades, model `grades.mod3`" knitr::include_graphics("fig/grades-HE3D.png")
Interactive
rotation of this plot shows that the effect of class is only two dimensional,
and of these, one dimension is very small. The major axis of the class ellipsoid
is aligned with increasing performance on all three grades, with the expected
ordering of the three social classes.
Canonical analysis
The representation of these effects in canonical space is particularly useful here.
Again, use candiscList() to compute the canonical decompositions for all terms
in the model, and extract the canonical $R^2$ from the terms in the result.
# calculate canonical results for all terms grades.can <- candiscList(grades.mod3) # extract canonical R^2s unlist(lapply(grades.can, function(x) x$canrsq))
We use heplot() on the "candiscList" object to show
the effects of class in canonical space, giving Figure \@ref(fig:grades-can-class).
#| echo=-1, #| out.width="80%", #| fig.cap="Canonical HE plot for `class` effect in `grades.mod3`" par(xpd=TRUE, mar=c(4,4,1,1)+.1) # plot class effect in canonical space heplot(grades.can, term="class", scale=4, fill=TRUE, var.col="black", var.lwd=2)
It can be seen in Figure \@ref(fig:grades-can-class) that nearly all variation in exam performance due to class is aligned with the first canonical dimension. The three tests and course evaluation all have similar weights on this dimension, but the course evaluation differs from the rest along a second, very small dimension.
1D plots of the canonical scores for other effects in the model are also of interest, and provide simple interpretations of these effects on the response variables. The statements below produce the plots shown in Figure \@ref(fig:grades-can-all).
#| out.width="47%", #| out.height="50%", #| fig.show="hold", #| collapse=FALSE, #| fig.cap="1D Canonical discriminant plots for `sex` and `gpa`. Higher canonical scores reflect better course performance." plot(grades.can, term="sex") plot(grades.can, term="gpa")
It is readily seen that males perform better overall, but the effect of
sex is strongest for the midterm2.
As well, increasing course performance on tests is strongly associated with
gpa.
References
Try the heplots package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.