ch10: linear hypotheses

10-mlm-review.qmd

@@ -50,28 +50,37 @@ where
 The structure of the model matrix $\mathbf{X}$ is the same as the univariate linear model, and may contain, therefore, 
 
 * quantitative predictors, such as `age`, `income`, years of `education`
-* transformed predictors like $\sqrt{\text{age}}$ or $\log{\text{income}}$
+* transformed predictors like $\sqrt{\text{age}}$ or $\log{(\text{income})}$
 * polynomial terms: $\text{age}^2$, $\text{age}^3, \dots$ (using `poly(age, k)` in R)
 * categorical predictors ("factors"), such as treatment (Control, Drug A, drug B), or sex; internally a factor with `k` levels is transformed to `k-1` dummy (0, 1) variables, representing comparisons with a reference level, typically the first.
 * interaction terms, involving either quantitative or categorical predictors, e.g., `age * sex`, `treatment * sex`.
 
 ### Assumptions
-Just as in univariate models, the assumptions of the multivariate linear model entirely concern the behavior of the errors (residuals). 
+Just as in univariate models, the assumptions of the multivariate linear model almost entirely concern the behavior of the errors (residuals). 
 Let $\mathbf{\epsilon}_{i}^{\prime}$ represent the $i$th row of $\Epsilon$. Then it is assumed that:
 
-* The residuals, $\mathbf{\epsilon}_{i}^{\prime}$ are distributed as multivariate normal, 
+* **Normality**: The residuals, $\mathbf{\epsilon}_{i}^{\prime}$ are distributed as multivariate normal, 
 $\mathcal{N}_{p}(\mathbf{0},\boldsymbol{\Sigma})$, 
 where $\mathbf{\Sigma}$ 
 is a non-singular error-covariance matrix. This can be assessed visually using a $\chi^2$ QQ plot
 of Mahalanobis squared distance against their corresponding $\chi^2_p$ values.
-* The error-covariance matrix $\mathbf{\Sigma}$ is constant across all observations and grouping factors.
+* **Homoscedasticity**: The error-covariance matrix $\mathbf{\Sigma}$ is constant across all observations and grouping factors.
 Graphical methods to determine if this assumption is met are described in @sec-eqcov.
-* $\mathbf{\epsilon}_{i}^{\prime}$ and $\mathbf{\epsilon}_{j}^{\prime}$ are independent for $i\neq j$; and 
+* **Independence**: $\mathbf{\epsilon}_{i}^{\prime}$ and $\mathbf{\epsilon}_{j}^{\prime}$ are independent for $i\neq j$; and 
 * The predictors, $\mathbf{X}$, are fixed or independent of $\Epsilon$. 
 
 These statements are simply the
 multivariate analogs of the assumptions of normality, constant variance and independence
-of the errors in univariate models.
+of the errors in univariate models. Note that it is unnecessary to assume that the predictors 
+(regressors, columns of $\mathbf{X}$) are normally distributed. 
+
+Implicit in the above is perhaps the most important assumption---that the model has been correctly specified. This means:
+
+* **Linearity**: The form of the relations between each $\mathbf{y}$ and the $\mathbf{x}$s is correct. Typically this means that the relations are _linear_, but if not, we have specified a correct transformation of $\mathbf{y}$ and/or $\mathbf{x}$.
+
+* **Completeness**: No relevant predictors have been omitted from the model.
+
+* **Additive effects**:  The combined effect of different predictors is the sum of their individual effects.
 
 ## Fitting the model
 
@@ -165,7 +174,8 @@ The same idea can be applied to test the difference between any pair of _nested_
 
 The import of the above is that when there are $p > 1$ response variables, and we are testing
 a hypothesis comprising $\text{df}_h >1$ coefficients or degrees of freedom, there are $s > 1$
-possible dimensions in which $\H$ can be large relative to $\E$. There are several test
+possible dimensions in which $\H$ can be large relative to $\E$, each measured by the
+eigenvalue $\lambda_i$. There are several test
 statistics that combine these into a single measure. 
 
 $$
@@ -177,28 +187,31 @@ $$
 \end{align*}
 $$
 
-Each of these statistics have different distributions under the null hypothesis, but they
+These correspond to different kinds of "means" of the $\lambda_i$: arithmetic, geometric, harmonic
+and supremum. See @Friendly-etal:ellipses:2013 for the geometry behind these measures.
+
+Each of these statistics have different sampling distributions under the null hypothesis, but they
 can all be converted to $F$ statistics, which are exact when $s \le 2$, and approximations otherwise.
 As well, each has an analog of the $R^2$-like partial $\eta^2$ measure, giving the partial association
 accounted for by each term in the MLM.
 
-<!-- Does the `parse-latex` filter work? -->
-
-<!-- ```{=latex} -->
-<!-- \begin{center} -->
-<!-- \begin{tabular}{|l|l|l|l|} -->
-<!--   \hline -->
-<!--   % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ... -->
-<!--   Criterion & Formula &  Partial $\eta^2$   \\ -->
-<!--   \hline -->
-<!--   Wilks's $\Lambda$ & $\Lambda = \prod^s_i \frac{1}{1+\lambda_i}$ &  $\eta^2 = 1-\Lambda^{1/s}$   \\ -->
-<!--   Pillai trace & $V = \sum^s_i \frac{\lambda_i}{1+\lambda_i}$ &  $\eta^2 = \frac{V}{s} $   \\ -->
-<!--   Hotelling-Lawley trace & $H = \sum^s_i \lambda_i$ & $\eta^2 = \frac{H}{H+s}$   \\ -->
-<!--   Roy maximum root & $R = \lambda_1$  &  $ \eta^2 = \frac{\lambda_1}{1+\lambda_1}  \\ -->
-<!--   \hline -->
-<!-- \end{tabular} -->
-<!-- \end{center} -->
-<!-- ``` -->
+Does the `parse-latex` filter work?
+
+```{=latex}
+\begin{center}
+\begin{tabular}{|l|l|l|l|}
+  \hline
+  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
+  Criterion & Formula &  Partial $\eta^2$   \\
+  \hline
+  Wilks's $\Lambda$ & $\Lambda = \prod^s_i \frac{1}{1+\lambda_i}$ &  $\eta^2 = 1-\Lambda^{1/s}$   \\
+  Pillai trace & $V = \sum^s_i \frac{\lambda_i}{1+\lambda_i}$ &  $\eta^2 = \frac{V}{s} $   \\
+  Hotelling-Lawley trace & $H = \sum^s_i \lambda_i$ & $\eta^2 = \frac{H}{H+s}$   \\
+  Roy maximum root & $R = \lambda_1$  &  $ \eta^2 = \frac{\lambda_1}{1+\lambda_1}$  \\
+  \hline
+\end{tabular}
+\end{center}
+```
 
 ### Testing contrasts and linear hypotheses
 
@@ -207,7 +220,7 @@ in $\mathbf{B}$.  Suppose we want to test the hypothesis that a subset of rows (
 and/or columns (responses) simultaneously have null effects.
 This can be expressed in the general linear test,
 $$
-H_0 : \mathbf{C}_{h \times q} \, \vec{B}_{q \times p} = \mat{0}_{h \times p} \comma
+H_0 : \mathbf{C}_{h \times q} \, \mathbf{B}_{q \times p} = \mathbf{0}_{h \times p} \comma
 $$
 where $\mathbf{C}$ is a  full rank $h \le q$ hypothesis matrix of constants, that selects 
 subsets or linear combinations (contrasts) of the coefficients in $\mathbf{B}$ to be tested
@@ -222,10 +235,10 @@ $$
 $${#eq-hmat}
 
 where there are $s = \min(h, p)$ non-zero eigenvalues of $\mat{H}\mat{E}^{-1}$.
-In @eq-hmat $\mathbd{H}$ measures the (Mahalanobis) squared distances (and cross products) among
+In @eq-hmat $\mathbf{H}$ measures the (Mahalanobis) squared distances (and cross products) among
 the linear combinations $\mathbf{C} \widehat{\mathbf{B}}$ from the origin under the null hypothesis.
 An animated display of these ideas and the relations between data ellipses and HE plots can be seen at
-http://www.datavis.ca/gallery/animation/manova/.
+[https://www.datavis.ca/gallery/animation/manova/](https://www.datavis.ca/gallery/animation/manova/).
 
 
 

R/parenting-ex.R

@@ -0,0 +1,118 @@
+#' ---
+#' title: "HE Plot Examples using Parenting Data"
+#' author: "Michael Friendly and Matthew Sigal"
+#' date: "18 Aug 2016"
+#' ---
+
+#' This script reproduces all of the analysis and graphs for the MANOVA of the parenting 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.
+
+#+ echo=FALSE
+library(knitr)
+library(rgl)
+knitr::opts_chunk$set(warning=FALSE, message=FALSE, R.options=list(digits=4))
+knitr::knit_hooks$set(webgl = hook_webgl)
+if(packageVersion("heplots") < "1.3.2") {
+    warning("Installing package: heplots 1.3.2+ from R-Forge")
+    install.packages("heplots", repos="http://R-Forge.R-project.org")
+}
+
+#' ## Load packages and the data
+library(car)
+library(heplots)
+data(Parenting, package="heplots")
+
+#' ## Initial view: side-by-side boxplots for a multivariate response
+library(reshape2)
+library(ggplot2)
+
+parenting.long <- melt(Parenting)
+
+ggplot(parenting.long, aes(x=group, y=value, fill=group)) +
+		geom_boxplot(outlier.size=2.5, alpha=.5) + 
+		stat_summary(fun.y=mean, colour="white", geom="point", size=2) +
+		facet_wrap(~ variable) +
+    theme_bw() + 
+    theme(legend.position="top") +
+    theme(axis.text.x = element_text(angle = 15, hjust = 1)) 
+
+#' ## Run the MANOVA
+
+# NB: order responses so caring, play are first
+parenting.mod <- lm(cbind(caring, play, emotion) ~ group, data=Parenting)
+Anova(parenting.mod)
+
+#' test linear hypotheses (contrasts)
+contrasts(Parenting$group)   # display the contrasts
+print(linearHypothesis(parenting.mod, "group1"), SSP=FALSE)
+print(linearHypothesis(parenting.mod, "group2"), SSP=FALSE)
+
+
+
+#' ## Figure 4: compare effect and significance scaling
+
+op <- par(mar=c(4,4,1,1)+0.1)
+res <- heplot(parenting.mod,
+		fill=TRUE, fill.alpha=c(0.3, 0.1),
+		lty = c(0,1),
+		cex=1.3, cex.lab=1.5)
+label <- expression(paste("Significance scaling:", H / lambda[alpha], df[e]))
+text(8.5, 11.5, label, cex=1.6)
+par(op)
+#dev.copy2pdf(file="parenting-HE2.pdf")
+
+op <- par(mar=c(4,4,1,1)+0.1)
+res <- heplot(parenting.mod, size="effect",
+              fill=TRUE, fill.alpha=c(0.3, 0.1), 
+              lty = c(0,1),
+              cex=1.3, cex.lab=1.5, label.pos=c(1,2),
+              xlim=res$xlim, ylim=res$ylim)
+label <- expression(paste("Effect size scaling:", H / df[e]))
+text(8.5, 11.5, label, cex=1.6)
+par(op)
+#dev.copy2pdf(file="parenting-HE1.pdf")
+
+#' ## Figure 5: showing contrasts
+# display tests of contrasts
+hyp <- list("N:MP" = "group1", "M:P" = "group2")
+
+# Fig 5: make a prettier heplot plot
+op <- par(mar=c(4,4,1,1)+0.1)
+heplot(parenting.mod, hypotheses=hyp, asp=1, 
+       fill=TRUE, fill.alpha=c(0.3,0.1), 
+       col=c("red", "blue"), 
+       lty=c(0,0,1,1), label.pos=c(1,1,3,2),
+       cex=1.4, cex.lab=1.4, lwd=3)
+par(op)
+#dev.copy2pdf(file="parenting-HE3.pdf")
+
+#' ## other HE plots not shown in paper
+pairs(parenting.mod, fill=TRUE, fill.alpha=c(0.3, 0.1))
+
+# This 3D plot should be interactive: zoom and rotate under mouse control
+
+#+ webgl=TRUE
+heplot3d(parenting.mod, wire=FALSE)
+
+#' ## Canonical discriminant analysis
+
+library(candisc)
+parenting.can <- candisc(parenting.mod)
+heplot(parenting.can)
+
+#' ## Figure 6
+op <- par(mar=c(5,4,1,1)+.1)
+pos <- c(4, 3, 3)
+
+heplot(parenting.can, 
+	fill=TRUE, fill.alpha=.1,  scale=6.5,
+	cex.lab=1.5, var.cex=1.2, 
+	var.lwd=2, var.col="black",  var.pos=pos,
+	label.pos=c(4, 3, 3), 
+	cex=1.2, 
+	xlim=c(-6,6), ylim=c(-6,6),
+	prefix="Canonical dimension ")
+par(op)
+#dev.copy2pdf(file="parenting-hecan.pdf")