Nothing
## ---- include = FALSE---------------------------------------------------------
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 = "#>"
)
## ----setup, echo=FALSE--------------------------------------------------------
set.seed(1071)
options(width=80, digits=4, continue=" ")
library(heplots)
library(candisc)
library(car)
library(dplyr)
library(tidyr)
library(ggplot2)
## ---- rohwer-some-------------------------------------------------------------
data(Rohwer)
Rohwer |> dplyr::sample_n(6)
## ----rohwer-long--------------------------------------------------------------
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
## ----rohwer-long-ggplot-------------------------------------------------------
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")
## ---- rohwer-separate---------------------------------------------------------
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)
## ---- rohwer-coef-------------------------------------------------------------
coef(rohwer.ses1)
coef(rohwer.ses2)
## ----rohwer-coefplot----------------------------------------------------------
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)
## ---- rohwer-HE1--------------------------------------------------------------
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")
## ---- rohwer-mod--------------------------------------------------------------
# MANCOVA, assuming equal slopes
rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ SES + n + s + ns + na + ss,
data=Rohwer)
Anova(rohwer.mod)
## ---- rohwer-cov-names--------------------------------------------------------
covariates <- c("n", "s", "ns", "na", "ss")
# or: covariates <- rownames(coef(rohwer.mod))[-(1:2)]
## ---- rohwer-mod-test---------------------------------------------------------
Regr <- linearHypothesis(rohwer.mod, covariates)
print(Regr, digits=4, SSP=FALSE)
## ---- rohwer-HE2--------------------------------------------------------------
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")
## ---- rohwer-HE3--------------------------------------------------------------
pairs(rohwer.mod, col=colors,
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")),
cex=1.3, lwd=c(2, rep(3,5), 4))
## ---- rohwer-HE3D-code, eval=FALSE--------------------------------------------
# colors <- c("pink", "blue", rep("black",5), "#969696")
# heplot3d(rohwer.mod, col=colors,
# hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))
## ----rohwer-HE3D--------------------------------------------------------------
knitr::include_graphics("fig/mmra-rohwer-HE3D.png")
## ---- rohwer-mod2-------------------------------------------------------------
rohwer.mod2 <- lm(cbind(SAT, PPVT, Raven) ~ SES * (n + s + ns + na + ss),
data=Rohwer)
Anova(rohwer.mod2)
## ---- rohwer-mod2-test--------------------------------------------------------
(coefs <- rownames(coef(rohwer.mod2)))
print(linearHypothesis(rohwer.mod2, coefs[grep(":", coefs)]), SSP=FALSE)
## ---- rohwer-HE4--------------------------------------------------------------
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)]))
## ----hern-str-----------------------------------------------------------------
data(Hernior)
str(Hernior)
## ---- hern1-------------------------------------------------------------------
Hern.mod <- lm(cbind(leave, nurse, los) ~ age + sex + pstat + build + cardiac + resp,
data=Hernior)
Anova(Hern.mod)
## ---- hern.rsq----------------------------------------------------------------
Hern.summary <- summary(Hern.mod)
unlist(lapply(Hern.summary, function(x) x$r.squared))
## ----hern-glance--------------------------------------------------------------
glance.mlm(Hern.mod)
## ---- hern2-------------------------------------------------------------------
# test overall regression
(predictors <- rownames(coef(Hern.mod))[-1])
Regr <- linearHypothesis(Hern.mod, predictors)
print(Regr, digits=5, SSP=FALSE)
## ---- hern-pairs--------------------------------------------------------------
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)
## ---- hern-canL---------------------------------------------------------------
Hern.canL <- candiscList(Hern.mod)
## ---- hern.canL2, eval=FALSE--------------------------------------------------
# plot(Hern.canL)
## ---- hern-can1---------------------------------------------------------------
plot(Hern.canL, term="pstat")
plot(Hern.canL, term="build")
## ---- hern-can2---------------------------------------------------------------
plot(Hern.canL, term="age")
plot(Hern.canL, term="cardiac")
## -----------------------------------------------------------------------------
str(SocGrades)
## ---- grades1-----------------------------------------------------------------
data(SocGrades)
grades.mod <- lm(cbind(midterm1, midterm2, final, eval) ~
class + sex + gpa + boards + hssoc + pretest,
data=SocGrades)
Anova(grades.mod, test="Roy")
## ---- grades2-----------------------------------------------------------------
grades.mod2 <- update(grades.mod, . ~ .^2)
Anova(grades.mod2, test="Roy")
## ---- grades3-----------------------------------------------------------------
grades.mod3 <- update(grades.mod, . ~ . + class:sex - hssoc - pretest)
Anova(grades.mod3, test="Roy")
## ---- grades-pairs------------------------------------------------------------
pairs(grades.mod3)
## ---- grades-HE3D-code, eval=FALSE--------------------------------------------
# heplot3d(grades.mod3, wire=FALSE)
## ---- grades-HE3D-------------------------------------------------------------
knitr::include_graphics("fig/grades-HE3D.png")
## ---- grades4-----------------------------------------------------------------
# calculate canonical results for all terms
grades.can <- candiscList(grades.mod3)
# extract canonical R^2s
unlist(lapply(grades.can, function(x) x$canrsq))
## ---- grades-can-class--------------------------------------------------------
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)
## ---- grades-can-all----------------------------------------------------------
plot(grades.can, term="sex")
plot(grades.can, term="gpa")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.