In friendly/heplots: Visualizing Hypothesis Tests in Multivariate Linear Models

% template for Sweave vignettes %\VignetteIndexEntry{HE Plots for Repeated Measures Designs} %\VignetteDepends{reshape, car, gplots, nlme, lattice, rgl} %\VignetteKeywords{data ellipse, HE plot, HE plot matrix, profile analysis, repeated measures, MANOVA, doubly-multivariate designs, mixed models} %\VignettePackage{heplots}

%\documentclass[12pt]{article} %\documentclass[article]{jss} \documentclass[nojss]{jss} \usepackage{float} \usepackage{amssymb, amsmath, amsfonts} %\usepackage{latexsym} %\usepackage[noae]{Sweave} % \usepackage{fancyvrb} % loaded by Sweave %\usepackage{color} % loaded by jss.cls \usepackage{url} %\usepackage[round]{natbib} % loaded by jss.cls \usepackage{upquote} % for Anova() output with `***' %\usepackage{hyperref} % %\bibpunct{(}{)}{;}{a}{,}{,} \usepackage{bm} % for bold math vectors/matrices \usepackage{array} % for \extrarowheight in Tables %\usepackage{geometry} %\geometry{left=1.1in, right=1.1in, top=1in, bottom=1in} %\usepackage{comment}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% document-specific definitions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% math stuff \newcommand{; | ;}{\ensuremath{\, | \,}} \renewcommand{\vec}[1]{\ensuremath{\bm{#1}}} \newcommand{\mat}[1]{\ensuremath{\bm{#1}}} \newcommand{\trans}{\ensuremath{^\mathsf{T}}} \newcommand{\diag}[1]{\ensuremath{\mathrm{diag} (#1)}} %\newcommand{\implies}{ \ensuremath{\mapsto} } \newcommand{\Beta}{B} \newcommand{\Epsilon}{E} \newcommand{\period}{\:\: .} \newcommand{\comma}{\:\: ,} \newcommand{\nvec}[2]{\ensuremath{{#1}{1}, {#1}{2},...,{#1}_{#2}}} \newcommand*{\Var}{\ensuremath{\mathsf{Var}}}

\newcommand{\sizedmat}[2]{% \mathord{\mathop{$\mathbf{#1}$}\limits_{(#2)}}% }

% abbreviations \newcommand{\LM}{LM} \newcommand{\MLM}{MVLM}

% provided by amsmath %\newenvironment{equation*}{\displaymath}{\enddisplaymath}%

% \ref extensions \newcommand{\tabref}[1]{Table~\@ref(#1)} \newcommand{\figref}[1]{Figure~\@ref(#1)} \newcommand{\secref}[1]{Section~\@ref(#1)} \renewcommand{\eqref}[1]{Eqn.~(\@ref(#1))}

% R stuff %\newcommand{\pkg}[1]{\textsf{#1} package} \renewcommand{\pkg}[1]{{\fontseries{b}\selectfont #1} package} \newcommand{\R}{\textsf{R}} %\newcommand{\code}[1]{{#1}} \newcommand{\func}[1]{{#1()}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% JSS preamble %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{HE Plots for Repeated Measures Designs} \author{Michael Friendly \ York University, Toronto} %\date{For \Rpackage{heplots} version r packageDescription("heplots")[["Version"]] , r Sys.Date()} \Plainauthor{Michael Friendly} %% comma-separated \Plaintitle{HE Plots for Repeated Measures Designs} %% without formatting %\Shorttitle{HE plots for Repeated Measures Designs} %% a short title (if necessary)

\definecolor{Soutput}{rgb}{0,0,0.56} \definecolor{Sinput}{rgb}{0.56,0,0} %\DefineVerbatimEnvironment{Sinput}{Verbatim}{formatcom={\color{Sinput}},fontsize=\footnotesize, baselinestretch=0.75} %\DefineVerbatimEnvironment{Soutput}{Verbatim}{formatcom={\color{Soutput}},fontsize=\footnotesize, baselinestretch=0.75} \DefineVerbatimEnvironment{Sinput}{Verbatim}{formatcom={\color{Sinput}},fontsize=\small, baselinestretch=0.75} \DefineVerbatimEnvironment{Soutput}{Verbatim}{formatcom={\color{Soutput}},fontsize=\small, baselinestretch=0.75}

\Abstract{ Hypothesis Error (HE) plots, introduced in @Friendly:07:manova, provide graphical methods to visualize hypothesis tests in multivariate linear models, by displaying hypothesis and error covariation as ellipsoids and providing visual representations of effect size and significance. These methods are implemented in the heplots for R [@Fox-etal:07:heplots] and SAS [@Friendly:06:hesoft], and apply generally to designs with fixed-effect factors (MANOVA), quantitative regressors (multivariate multiple regression) and combined cases (MANCOVA).

This paper describes the extension of these methods to repeated measures designs
in which the multivariate responses represent the outcomes on one or more
``within-subject'' factors.  This extension is illustrated using
the `heplots` for R.  Examples describe one-sample profile analysis,
designs with multiple between-S and within-S factors, and doubly-multivariate designs,
with multivariate responses observed on multiple occasions.

}

\Keywords{data ellipse, HE plot, HE plot matrix, profile analysis, repeated measures, MANOVA, doubly-multivariate designs, mixed models} \Plainkeywords{data ellipse, HE plot, HE plot matrix, profile analysis, repeated measures, MANOVA, doubly-multivariate designs, mixed models} %% without formatting %% at least one keyword must be supplied

%% publication information %% NOTE: This needs to filled out ONLY IF THE PAPER WAS ACCEPTED. %% If it was not (yet) accepted, leave them commented. \Volume{37} \Issue{4} \Month{November} \Year{2010} \Submitdate{2009-12-22} \Acceptdate{2010-01-10}

%% The address of (at least) one author should be given %% in the following format: \Address{ Michael Friendly\ Psychology Department\ York University\ Toronto, ON, M3J 1P3\ Canada \ E-mail: \email{friendly```yorku.ca}\ URL: [http://datavis.ca] } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: [phone] %% Fax: [phone]

%% for those who use Sweave please include the following line (with % symbols): %% need no \usepackage{Sweave.sty}

%\SweaveOpts{engine=R,height=6,width=6,results=hide,fig=FALSE,echo=TRUE} \SweaveOpts{engine=R,height=6,width=6,results=hide,fig=FALSE,echo=TRUE,eps=FALSE} \SweaveOpts{prefix.string=fig/plot} \setkeys{Gin}{width=0.65\textwidth}

set.seed(1071)
old.opts <- options(width=85, digits=6, useFancyQuotes = FALSE, continue="  ", prompt="R> ")
library(heplots)
library(lattice)
library(nlme)
library(car)
# avoid deprecated warning
if (car2 <- packageDescription("car")[["Version"]] >= 2) data.ellipse <- dataEllipse

Introduction

Hypothesis Error (HE) plots, introduced in @Friendly:07:manova,
provide graphical methods to visualize hypothesis tests in multivariate linear models,
by displaying hypothesis and error covariation as ellipsoids and providing  visual
representations of effect size and significance.

The heplots [@Fox-etal:07:heplots] for R [@R:2010] implements these methods for the general class of the multivariate linear model (\MLM) including fixed-effect factors (MANOVA), quantitative regressors (multivariate multiple regression (MMREG)) and combined cases (MANCOVA). Here, we describe the extension of these methods to repeated measures designs in which the multivariate responses represent the outcomes on one or more ``within-subject'' factors. This vignette also appears in the Journal of Statistical Software [@Friendly:10:JSS].

Multivariate linear models: Notation

To set notation, we express the \MLM\ as {#eq:mglm} \sizedmat{Y}{n \times p} = \sizedmat{X}{n \times q} \sizedmat{\Beta}{q \times p} + \sizedmat{U}{n \times p} \comma

where, $$\mathbf{Y}$ \equiv (\nvec{\vec{y}}{p})$ is the matrix of responses for $n$ subjects on $p$ variables, $$\mathbf{X}$$ is the design matrix for $q$ regressors, $$\mathbf{\Beta}$$ is the $q \times p$ matrix of regression coefficients or model parameters and $$\mathbf{U}$$ is the $n \times p$ matrix of errors, with $\textrm{vec}($\mathbf{U}$) \sim \mathcal{N}_p ( $\mathbf{0}$, $\mathbf{I}$_n \otimes $\mathbf{\Sigma}$ )$, where $\otimes$ is the Kronecker product.

A convenient feature of the \MLM\ for general multivariate responses is that all tests of linear hypotheses (for null effects) can be represented in the form of a general linear test, {#eq:mglt} H_0 : \sizedmat{L}{h \times q} \sizedmat{\Beta}{q \times p} = \sizedmat{0}{h \times p} \comma

where $$\mathbf{L}$$ is a matrix of constants whose rows specify $h$ linear combinations or contrasts of the parameters to be tested simultaneously by a multivariate test. In R all such tests can be carried out using the functions Anova() and linearHypothesis() in the car.% \footnote{ Both the car and the heplots are being actively developed. Except where noted, all results in this paper were produced using the old-stable versions on CRAN, car 1.2-16 (2009/10/10) % r packageDescription("car")[["Version"]] (r (packageDescription("car")[["Date"]])) and heplots 0.8-11 (2009-12-08) % r packageDescription("heplots")[["Version"]] (r (packageDescription("heplots")[["Date"]])) running under r R.version$version.string. }

For any such hypothesis of the form \Eqn. @ref(eq:mglt), the analogs of the univariate sums of squares for hypothesis ($SS_H$) and error ($SS_E$) are the $p \times p$ sum of squares and crossproducts (SSP) matrices given by \citep[Ch. 3,5]{Timm:75}: {#eq:hmat} $\mathbf{H}$ \equiv $\mathbf{SSP}$_H = ($\mathbf{L}$ \widehat{$\mathbf{B}$})\trans \, [$\mathbf{L}$ ($\mathbf{X}$\trans $\mathbf{X}$ )^{-} $\mathbf{L}$\trans]^{-1} \, ($\mathbf{L}$ \widehat{$\mathbf{B}$}) \comma

and {#eq:emat} $\mathbf{E}$ \equiv $\mathbf{SSP}$_E = $\mathbf{Y}$\trans $\mathbf{Y}$ - \widehat{$\mathbf{B}$}\trans ($\mathbf{X}$\trans $\mathbf{X}$) \widehat{$\mathbf{B}$} = \widehat{$\mathbf{U}$}\trans \widehat{$\mathbf{U}$} \comma

where $\widehat{$\mathbf{U}$} = $\mathbf{Y}$ - $\mathbf{X}$ \widehat{$\mathbf{B}$}$ is the matrix of residuals. Multivariate test statistics (Wilks' $\Lambda$, Pillai trace, Hotelling-Lawley trace, Roy's maximum root) for testing \Eqn. @ref(eq:mglt) are based on the $s = \min(p, h)$ non-zero latent roots of $$\mathbf{H}$ $\mathbf{E}$^{-1}$ and attempt to capture how ``large'' $$\mathbf{H}$$ is, relative to $$\mathbf{E}$$ in $s$ dimensions. All of these statistics have transformations to $F$ statistics giving either exact or approximate null hypothesis $F$ distributions. The corresponding latent vectors give a set of $s$ orthogonal linear combinations of the responses that produce maximal univariate $F$ statistics for the hypothesis in \Eqn. @ref(eq:mglt); we refer to these as the canonical discriminant dimensions.

In a univariate, fixed-effects linear model, it is common to provide $F$ tests for each term in the model, summarized in an analysis-of-variance (ANOVA) table. The hypothesis sums of squares, $SS_H$, for these tests can be expressed as differences in the error sums of squares, $SS_E$, for nested models. For example, dropping each term in the model in turn and contrasting the resulting residual sum of squares with that for the full model produces so-called Type-III tests; adding terms to the model sequentially produces so-called Type-I tests; and testing each term after all terms in the model with the exception of those to which it is marginal produces so-called Type-II tests. Closely analogous MANOVA tables can be formed similarly by taking differences in error sum of squares and products matrices ($\mathbf{E}$) for such nested models. Type I tests are sensible only in special circumstances; in balanced designs, Type II and Type III tests are equivalent. Regardless, the methods illustrated in this paper apply to any multivariate linear hypothesis.

Data ellipses and ellipsoids

In what follows, we make extensive use of ellipses (or ellipsoids in 3+D) to represent joint variation among two or more variables, so we define this here. The data ellipse (or covariance ellipse), described by @Dempster:69 and @Monette:90, is a device for visualizing the relationship between two variables, $Y_{1}$ and $Y_{2}$. Let $D_{M}^{2}(% \vec{y})=(\vec{y}-\overline{\vec{y}})\trans$\mathbf{S}$^{-1}(\vec{y}-% \overline{\vec{y}})$ represent the squared Mahalanobis distance of the point $\vec{y}=(y_{1},y_{2})\trans$ from the centroid of the data $% \overline{\vec{y}}=(\overline{Y}{1},\overline{Y}{2})\trans$. The data ellipse $\mathcal{E}{c}$ of size $c$ is the set of all points $\vec{y}$ with $D{M}^{2}(\vec{y})$ less than or equal to $c^{2}$:%

\mathcal{E}_{c}(\mathbf{y;S,}\overline{\vec{y}})\equiv\left{ y\text{: }(% \vec{y}-\overline{\vec{y}})\trans$\mathbf{S}$^{-1}(\vec{y}-\overline{% \vec{y}})\leq c^{2}\right} {#eqn-data-ellipse}

Here, $$\mathbf{S}$=\sum_{i=1}^{n}(\vec{y}-\overline{\vec{y}})\trans(% \vec{y}-\overline{\vec{y}})/(n-1) = \widehat{\Var}(\vec{y})$ is the sample covariance matrix.

Many properties of the data ellipse hold regardless of the joint distribution of the variables, but if the variables are bivariate normal, then the data ellipse represents a contour of constant density in their joint distribution. In this case, $D_{M}^{2}(\vec{y})$ has a large-sample $\chi^{2}$ distribution with 2 degrees of freedom, and so, for example, taking $c^{2}=\chi_{2}^{2}(0.95)=5.99\approx6$ encloses approximately 95 percent of the data. Taking $c^{2}=\chi_{2}^{2}(0.68)=2.28$ gives a bivariate analog of the univariate $\pm 1$ standard deviation interval, enclosing approximately 68% of the data.

The generalization of the data ellipse to more than two variables is immediate: Applying Equation \@ref(eqn-data-ellipse) to $\vec{y}% =(y_{1},y_{2},y_{3})\trans$, for example, produces a data ellipsoid in three dimensions. For $p$ multivariate-normal variables, selecting $c^{2}=\chi _{p}^{2}(1-\alpha)$ encloses approximately $100(1-\alpha)$ percent of the data.% \footnote{ Robust versions of data ellipses (e.g., based on minimum volume ellipsoid (MVE) or minimum covariance determinant (MCD) estimators of $\mathbf{S}$) are also available, as are small-sample approximations to the enclosing $c^2$ radii, but these refinements are outside the scope of this paper. }

HE plots

The essential idea behind HE plots is that any multivariate hypothesis test \Eqn. @ref(eq:mglt) can be represented visually by ellipses (or ellipsoids in 3D) which express the size of co-variation against a multivariate null hypothesis ($$\mathbf{H}$$) relative to error covariation ($$\mathbf{E}$$). The multivariate tests, based on the latent roots of $$\mathbf{H}$ $\mathbf{E}$^{-1}$, are thus translated directly to the sizes of the $\mathbf{H}$ ellipses for various hypotheses, relative to the size of the $\mathbf{E}$ ellipse. Moreover, the shape and orientation of these ellipses show something more-- the directions (linear combinations of the responses) that lead to various effect sizes and significance.

In these plots, the $\mathbf{E}$ matrix is first scaled to a covariance matrix ($$\mathbf{E}$/df_e = \widehat{\Var}($\mathbf{U}$_i)$). The ellipse drawn (translated to the centroid $\overline{\vec{y}}$ of the variables) is thus the data ellipse of the residuals, reflecting the size and orientation of residual variation. In what follows (by default), we always show these as ``standard'' ellipses of 68% coverage. This scaling and translation also allows the means for levels of the factors to be displayed in the same space, facilitating interpretation.

The ellipses for $\mathbf{H}$ reflect the size and orientation of covariation against the null hypothesis. They always proportional to the data ellipse of the fitted effects (predicted values) for a given hypothesized term. In relation to the $\mathbf{E}$ ellipse, the $\mathbf{H}$ ellipses can be scaled to show either the effect size or strength of evidence against $H_0$ (significance).

For effect size scaling, each $\mathbf{H}$ is divided by $df_e$ to conform to $\mathbf{E}$. The resulting ellipses are then exactly the data ellipses of the fitted values, and correspond visually to multivariate analogs of univariate effect size measures (e.g., $(\bar{y_1}-\bar{y_2})/s$ where $s$=within group standard deviation). That is, the sizes of the $\mathbf{H}$ ellipses relative to that of the $\mathbf{E}$ reflect the (squared) differences and correlation of the factor means relative to error covariation.

For significance scaling, it turns out to be most visually convenient to use Roy's largest root test as the test criterion. In this case the $\mathbf{H}$ ellipse is scaled to $$\mathbf{H}$/(\lambda_\alpha df_e)$ where $\lambda_\alpha$ is the critical value of Roy's statistic. Using this gives a simple visual test of $H_0$: Roy's test rejects $H_0$ at a given $\alpha$ level if and only if the corresponding $\alpha$-level $\mathbf{H}$ ellipse extends anywhere outside the $\mathbf{E}$ ellipse.% \footnote{Other multivariate tests (Wilks' $\Lambda$, Hotelling-Lawley trace, Pillai trace) also have geometric interpretations in HE plots (e.g., Wilks' $\Lambda$ is the ratio of areas (volumes) of the $\mathbf{H}$ and $\mathbf{E}$ ellipses (ellipsoids)), but these statistics do not provide such simple visual comparisons. All HE plots shown in this paper use significance scaling, based on Roy's test. } Consequently, when the rank of $$\mathbf{H}$ = \min(p, h) \le 2$, all significant effects can be observed directly in 2D HE plots; when rank$($\mathbf{H}$)=3$, some rotation of a 3D plot will reveal each significant effect as extending somewhere outside the $\mathbf{E}$ ellipsoid.

In our R implementation, the basic plotting functions in the heplots are heplot() and heplot3d() for mlm objects. These rely heavily on the Anova() and other functions from the car [@car] for computation. For more than three response variables, all pairwise HE plots can be shown using a pairs() function for mlm objects. Alternatively, the related candisc [@candisc] produces HE plots in canonical discriminant space. This shows a low-rank 2D (or 3D) view of the effects for a given term in the space of maximum discrimination, based on the linear combinations of responses which produce maximally significant test statistics. See @Friendly:07:manova,FoxFriendlyMonette:09:compstat for details and examples for between-S MANOVA designs, MMREG and MANCOVA models.

Repeated measures designs {#sec:repmeas}

The framework for the \MLM\ described above pertains to the situation in which the response vectors (rows, $\vec{y}i\trans$ of $$\mathbf{Y}${n \times p})$ are iid and the $p$ responses are separate, not necessarily commensurate variables observed on individual $i$.

In principle, the \MLM\ extends quite elegantly to repeated-measure (or within-subject) designs, in which the $p$ responses per individual can represent the factorial combination of one or more factors that structure the response variables in the same way that the between-individual design structures the observations. In the multivariate approach to repeated measure data, the same model \Eqn. @ref(eq:mglm) applies, but hypotheses about between- and within-individual variation are tested by an extended form of the general linear test \Eqn. @ref(eq:mglt), which becomes {#eq:mglt-rep} H_0 : \sizedmat{L}{h \times q} \sizedmat{\Beta}{q \times p} \sizedmat{M}{p \times k} = \sizedmat{0}{h \times k} \comma

where $$\mathbf{M}$$ is a matrix of constants whose columns specify $k$ linear combinations or contrasts among the responses, corresponding to a particular within-individual effect. In this case, the $\mathbf{H}$ and $\mathbf{E}$ matrices for testing \Eqn. @ref(eq:mglt-rep) become {#eq:hmat-rep} $\mathbf{H}$ = ($\mathbf{L}$ \widehat{$\mathbf{B}$} $\mathbf{M}$)\trans \, [$\mathbf{L}$ ($\mathbf{X}$\trans $\mathbf{X}$ )^{-} $\mathbf{L}$\trans]^{-1} \, ($\mathbf{L}$ \widehat{$\mathbf{B}$} $\mathbf{M}$) \comma

and {#eq:emat-rep} $\mathbf{E}$ = ($\mathbf{Y}$ $\mathbf{M}$)\trans [$\mathbf{I}$ - ($\mathbf{X}$\trans $\mathbf{X}$ )^{-} $\mathbf{X}$\trans] \, ($\mathbf{Y}$ $\mathbf{M}$) \period

This may be easily seen to be just the ordinary \MLM\ applied to the transformed responses $$\mathbf{Y}$ $\mathbf{M}$$ which form the basis for a given within-individual effect. The idea for this approach to repeated measures through a transformation of the responses was first suggested by @Hsu:38 and is discussed further by @Rencher:95 and @Timm:80. In what follows, we refer to hypotheses pertaining to between-individual effects (specified by $\mathbf{L}$) as between-S'' and hypotheses pertaining to within-individual effects ($\mathbf{M}$) aswithin-S.''

[htb] \renewcommand{\arraystretch}{1.5} \renewcommand{\tabcolsep}{0.5cm} % \vspace{6pt} {#tab:BWtab} \setlength{\extrarowheight}{4pt} {|l|l|l|l|l|} \hline % after \: \hline or \cline{col1-col2} \cline{col3-col4} ... & \multicolumn{4}{c|}{Between-S effect tested} \ \cline{2-5} $$\mathbf{M}$$ for Within-S effects & Intercept & $$\mathbf{L}$ = $\mathbf{L}$A$ & $$\mathbf{L}$ = $\mathbf{L}$_B$ & $$\mathbf{L}$ = $\mathbf{L}${AB}$ \[1ex] \hline $$\mathbf{M}$ = $\mathbf{M}$_1 = \left( {rrr} 1 & 1 & 1 \

                \right)\trans

$ & $\vec{\mu}_{..}$ & A & B & A:B \[.5ex] % \hline %\rule{0pt}{5pt} $$\mathbf{M}$ = $\mathbf{M}$_C = \left( {rrr} 1 & -1 & 0 \ 0 & 1 & -1 \

               \right)\trans

$ & C & A:C & B:C & A:B:C \ \hline

\caption{Three-way design: Tests for between- (A, B) and within-S (C) effects are constructed using various $\mathbf{L}$ and $\mathbf{M}$ matrices. Table entries give the term actually tested via the general linear test in \Eqn. @ref(eq:mglt-rep).}

In the general case, various $\mathbf{L}$ matrices provide contrasts or select the particular coefficients tested for between-S effects, while various $\mathbf{M}$ matrices specify linear combinations of responses for the within-S effects. This is illustrated in \Table @ref(tab:BWtab) for a three-way design with two between-S factors (A, B) and one within-S factor (C).

The between-S terms themselves are tested using the unit vector $$\mathbf{M}$ = ( \vec{1}_p )$, giving a test based on the sums over the within-S effects. This simply reflects the principle of marginality, by which effects for any term in a linear model are tested by averaging over all factors not included in that term. Tests using a matrix $$\mathbf{M}$$ of contrasts for a within-S effect provide tests of the interactions of that effect with each of the between-S terms. That is, $$\mathbf{L B M}$=$\mathbf{0}$$ tests between-S differences among the responses transformed by $\mathbf{M}$.

For more than one within-S factor, the full $\mathbf{M}$ matrices for various within-S terms are generated as Kronecker products of the one-way $\mathbf{M}$ contrasts with the unit vector (\vec{1}) of appropriate size. For example, with $c$ levels of factor C and $d$ of factor D, {#eq:Mcd}

$\mathbf{M}$_{C\otimes D}  & =  (\vec{1}_c, $\mathbf{M}$_C ) \otimes (\vec{1}_d, $\mathbf{M}$_D ) \\
    & =  ( \vec{1}_c \otimes \vec{1}_d, \vec{1}_c \otimes $\mathbf{M}$_D,
      $\mathbf{M}$_C \otimes \vec{1}_d, $\mathbf{M}$_C \otimes $\mathbf{M}$_D ) \\
    & = ( $\mathbf{M}$_1, $\mathbf{M}$_D, $\mathbf{M}$_C, $\mathbf{M}$_{CD} ) \period

Each of the within-S terms combine with any between-S terms in an obvious way to give an extended version of \Table @ref(tab:BWtab) with additional rows for $$\mathbf{M}$D$ and $$\mathbf{M}${CD}$.

In passing, we note that all software (SAS, SPSS, R, etc.) that handles repeated measure designs through this extension of the MLM effectively works via the general linear test \Eqn. @ref(eq:mglt), with either implicit or explicit specifications for the $\mathbf{L}$ and $\mathbf{M}$ matrices involved in testing any hypothesis for between- or within-S effects. This mathematical elegance is not without cost, however. The MLM approach does not allow for missing data (a particular problem in longitudinal designs), and the multivariate test statistics (Wilks' $\Lambda$, etc.) assume the covariance matrix of $\mathbf{U}$ is unstructured. Alternative analysis based on mixed (or hierarchical) models, e.g. [@PinheiroBates:2000,VerbekeMolenberghs:00] are more general in some ways, but to date visualization methods for this approach remain primitive and the mixed model analysis does not easily accommodate multivariate responses.

The remainder of the paper illustrates these MLM analyses, shows how they may be performed in \R, and how HE plots can be used to provide visual displays of what is summarized in the multivariate test statistics. We freely admit that these displays are somewhat novel and take some getting used to, and so this paper takes a more tutorial tone. We exemplify these methods in the context of simple, one-sample profile analysis (\@ref(sec:profile)), designs with multiple between- and within-S effects (\@ref(sec:BW)), and doubly-multivariate designs (\@ref(sec:doubly)), where two or more separate responses (e.g., weight loss and self esteem) are each observed in a factorial structure over multiple within-S occasions. In \@ref(sec:addendum) we describe a simplified interface for these plots in the development versions of the heplots and car packages. Finally (\@ref(sec:compare)) we compare these methods with visualizations based on the mixed model.

One sample profile analysis {#sec:profile}

The simplest case of a repeated-measures design is illustrated by the data on vocabulary growth of school children from grade 8 to grade 11, in the data frame VocabGrowth, recording scores on the vocabulary section of the Cooperative Reading Test for a cohort of 64 students. (The scores are scaled to a common, but arbitrary origin and unit of measurement, so as to be comparable over the four grades.) Since these data cover an age range in which physical growth is beginning to decelerate, it is of interest whether a similar effect occurs in the acquisition of new vocabulary. Thus, attention here is arguably directed to polynomial trends in grade: average rate of change (slope, or linear trend) and shape of trajectories (quadratic and cubic components).

some(VocabGrowth,5)

A boxplot of these scores (\Figure @ref(fig:voc-boxplot)) gives an initial view of the data. To do this, we first reshape the data from wide to long format (i.e., each 4-variate row becomes four rows indexed by grade). We can see that vocabulary scores increase with age, but the trend of means appears non-linear.

%r %library(reshape) %voc <- melt(VocabGrowth) %colnames(voc) <- c("Grade", "Vocabulary") %voc$Grade <- factor(as.numeric(substr(voc$Grade,6,8))) %plot(voc, main="Vocabulary Growth data") %abline(lm(Vocabulary ~ as.numeric(Grade), data=voc), col="red") %

% profile plot % r % xyplot(Vocabulary ~ Grade, groups=id, data=voc, type='b', pch=16:25) % %

voc <- reshape(VocabGrowth, direction="long", varying=list(grade=1:4), timevar="Grade", v.names="Vocabulary")
boxplot(Vocabulary ~ Grade, data=voc, col="bisque",
    ylab="Vocabulary", main="Vocabulary Growth data")
abline(lm(Vocabulary ~ as.numeric(Grade), data=voc), col="red")
means <- tapply(voc$Vocabulary, voc$Grade, mean)
points(1:4, means, pch=7, col="blue")
lines(1:4, means, col="blue", lwd=2)

%\setkeys{Gin}{width=0.7\textwidth} [htb] \centering \includegraphics[width=0.6\textwidth]{fig/plot-voc1} \caption{Boxplots of Vocabulary score by Grade, with linear regression line (red) and lines connecting grade means (blue).} {#fig:voc-boxplot}

The standard univariate and multivariate tests for the differences in vocabulary with grade can be carried out as follows. First, we fit the basic \MLM\ with an intercept only on the right-hand side of the model, since there are no between-S effects. The intercepts estimate the means at each grade level, $\mu_8 , \dots , \mu_{11}$.

(Vocab.mod <- lm(cbind(grade8,grade9,grade10,grade11) ~ 1, data=VocabGrowth))

We could test the multivariate hypothesis that all means are simultaneously zero, $\mu_8 = \mu_9 = \mu_{10} = \mu_{11} = 0$. This point hypothesis is the simplest case of a multivariate test under \Eqn. @ref(eq:mglt), with $$\mathbf{L}$ = $\mathbf{I}$$.

(Vocab.aov0 <- Anova(Vocab.mod, type="III"))

This hypothesis tests that the vocabulary means are all at the arbitrary origin for the scale. Often this test is not of direct interest, but it serves to illustrate the $\mathbf{H}$ and $\mathbf{E}$ matrices involved in any multivariate test, their representation by HE plots, and how we can extend these plots to the repeated measures case.

The $\mathbf{H}$ and $\mathbf{E}$ matrices can be printed with summary(Vocab.aov0), or extracted from the Anova.mlm object. In this case, $\mathbf{H}$ is simply $n \overline{\vec{y}} \overline{\vec{y}}\trans$ and $\mathbf{E}$ is the sum of squares and crossproducts of deviations from the column means, $ \sum_{i=1}^{n}(\vec{y}_i-\overline{\vec{y}})\trans(% \vec{y}_i-\overline{\vec{y}})$.

Vocab.aov0$SSP     # H matrix
Vocab.aov0$SSPE    # E matrix

The HE plot for the Vocab.mod model shows the test for the (Intercept) term (all means = 0). To emphasize that the test is assessing the (squared) distance of $\bar{\vec{y}}$ from $\vec{0}$, in relation to the covariation of observations around the grand mean, we define a simple function to mark the point hypothesis $H_0 = (0, 0)$.

mark.H0 <- function(x=0, y=0, cex=2, pch=19, col="green3", lty=2, pos=2) {
  points(x,y, cex=cex, col=col, pch=pch)
  text(x,y, expression(H[0]), col=col, pos=pos)
  if (lty>0) abline(h=y, col=col, lty=lty)
  if (lty>0) abline(v=x, col=col, lty=lty)
}

Here we show the HE plot for the grade8 and grade9 variables in \Figure @ref(fig:voc4). The $\mathbf{E}$ ellipse reflects the positive correlation of vocabulary scores across these two grades, but also shows that variability is greater in Grade 8 than in Grade 9. Its position relative to (0,0) indicates that both means are positive, with a larger mean at Grade 9 than Grade 8.

heplot(Vocab.mod, terms="(Intercept)", type="III")
mark.H0(0,0)
title(expression(paste("Multivariate test of ", H[0], " : ", bold(mu)==0)))

[htb] \centering \includegraphics[width=.65\textwidth]{fig/plot-voc4} \caption{HE plot for Vocabulary data, for the MANOVA test of $H_0 : \vec{\mu}_y=0$. The size of the (degenerate) ellipse for the intercept term relative to that for Error gives the strength of evidence for the difference between the sample means (marked by $+$) and the means under $H_0$ (marked by the cross-hairs and green dot). The projection of this $\mathbf{H}$ ellipse outside the $\mathbf{E}$ ellipse signals that this $H_0$ is clearly rejected.} {#fig:voc4}

The $\mathbf{H}$ ellipse plots as a degenerate line because the $\mathbf{H}$ matrix has rank 1 (1 df for the MANOVA test of the Intercept). The fact that the $\mathbf{H}$ ellipse extends outside the $\mathbf{E}$ ellipse (anywhere) signals that this $H_0$ is clearly rejected (for some linear combination of the response variables). Moreover, the projections of the $\mathbf{H}$ and $\mathbf{E}$ ellipses on the grade8 and grade9 axes, showing $\mathbf{H}$ widely outside $\mathbf{E}$, signals that the corresponding univariate hypotheses, $\mu_8=0$ and $\mu_9=0$ would also be rejected.

Testing within-S effects {#sec:within}

For the Anova() function, the model for within-S effects--- giving rise the $\mathbf{M}$ matrices in \Eqn. @ref(eq:mglt-rep)--- is specified through the arguments idata (a data frame giving the factor(s) used in the intra-subject model) and idesign (a model formula specifying the intra-subject design). That is, if $$\mathbf{Z}$ = [ \mathtt{idata} ]$, the $\mathbf{M}$ matrices for different intra-subject terms are generated from columns of $\mathbf{Z}$ indexed by the terms in idesign, with factors and interactions expanded expanded to contrasts in the same way that the design matrix $\mathbf{X}$ is generated from the between-S design formula.

Thus, to test the within-S effect of grade, we need to construct a grade variable for the levels of grade and use this as a model formula, idesign=~grade to specify the within-S design in the call to Anova.

idata <-data.frame(grade=ordered(8:11))
(Vocab.aov <- Anova(Vocab.mod, idata=idata, idesign=~grade, type="III"))

As shown in \Table @ref(tab:BWtab), any such within-S test is effectively carried out using a transformation $\mathbf{Y}$ to $\mathbf{Y M}$, where the columns of $\mathbf{M}$ provide contrasts among the grades. For the overall test of grade, any set of 3 linearly independent contrasts will give the same test statistics, though, of course the interpretation of the parameters will differ. Specifying grade as an ordered factor (grade=ordered(8:11)) will cause Anova() to use the polynomial contrasts shown in $$\mathbf{M}${poly}$ below. [ $\mathbf{M}${poly} = \left( {rrr} -3 & 1 & -1 \ -1 & -1 & 3 \ 1 & -1 & -3 \ 3 & 1 & 1 \

\right) \quad\quad $\mathbf{M}$_{first} = \left( {rrr} -1 & -1 & -1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \

\right) \quad\quad $\mathbf{M}$_{last} = \left( {rrr} 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \ -1 & -1 & -1 \

\right) ] Alternatively, $$\mathbf{M}${first}$ would test the gains in vocabulary between grade 8 (baseline) and each of grades 9--11, while $$\mathbf{M}${last}$ would test the difference between each of grades 8--10 from grade 11. (In \R, these contrasts are constructed with M.first <- contr.sum(factor(11:8))[4:1,3:1], and M.last <- contr.sum(factor(8:11)) respectively.) In all cases, the hypothesis of no difference among the means across grade is transformed to an equivalent multivariate point hypothesis, $$\mathbf{M}$ \vec{\mu}_y =\vec{0}$, such as we visualized in \Figure @ref(fig:voc4).

Correspondingly, the HE plot for the effect of grade can be constructed as follows. For expository purposes we explicitly transform $\mathbf{Y}$ to $\mathbf{Y M}$, where the columns of $\mathbf{M}$ provide contrasts among the grades reflecting linear, quadratic and cubic trends using $$\mathbf{M}$_{poly}$.

%[ %$\mathbf{M}$ = \bordermatrix{ %& \text{lin} & \text{quad} & \text{cube} \cr %& -3 & 1 & -1 \cr %& -1 & -1 & 3 \cr %& 1 & -1 & -3 \cr %& 3 & 1 & 1 %} %]

Using $$\mathbf{M}${poly}$, the MANOVA test for the grade effect is then testing $H_0 : $\mathbf{M}$ \vec{\mu}_y =\vec{0} \leftrightarrow \mu{Lin} = \mu_{Quad} = \mu_{Cubic} = 0$. That is, the means across grades 8--11 are equal if and only if their linear, quadratic and cubic trends are simultaneously zero.

trends <- as.matrix(VocabGrowth) %*% poly(8:11, degree=3)
colnames(trends)<- c("Linear", "Quad", "Cubic")

# test all trend means = 0 == Grade effect
within.mod <- lm(trends ~ 1)

Anova(within.mod, type="III")

It is easily seen that the test of the (Intercept) term in within.mod is identical to the test of grade in Vocab.mod at the beginning of this subsection.

We can show this test visually as follows. \Figure @ref(fig:voc8)(a) shows a scatterplot of the transformed linear and quadratic trend scores, overlayed with a 68% data ellipse. \Figure @ref(fig:voc8)(b) is the corresponding HE plot for these two variables. Thus, we can see that the $\mathbf{E}$ ellipse is simply the data ellipse of the transformed vocabulary scores; its orientation indicates a slight tendency for those with greater slopes (gain in vocabulary) to have greater curvatures (leveling off earlier). \Figure @ref(fig:voc8) is produced as follows:

op <- par(mfrow=c(1,2))
data.ellipse(trends[,1:2], xlim=c(-4,8), ylim=c(-3,3), levels=0.68,
    main="(a) Data ellipse ")
mark.H0(0,0)

heplot(within.mod, terms="(Intercept)", col=c("red", "blue"), type="III",
    term.labels="Grade", , xlim=c(-4,8), ylim=c(-3,3),
    main="(b) HE plot for Grade effect")
mark.H0(0,0)
par(op)

[htb] \centering \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-voc8} \caption{Plots of linear and quadratic trend scores for the Vocabulary data. (a) Scatterplot with 68% data ellipse; (b) HE plot for the effect of Grade. As in \Figure @ref(fig:voc4), the size of the $\mathbf{H}$ ellipse for Grade relative to the $\mathbf{E}$ ellipse shows the strength of evidence against $H_0$.} {#fig:voc8}

We interpret \Figure @ref(fig:voc8)(b) as follows, bearing in mind that we are looking at the data in the transformed space ($$\mathbf{Y M}$$) of the linear (slopes) and quadratic (curvatures) of the original data ($\mathbf{Y}$). The mean slope is positive while the mean quadratic trend is negative. That is, overall, vocabulary increases with Grade, but at a decreasing rate. The $\mathbf{H}$ ellipse plots as a degenerate line because the $\mathbf{H}$ matrix has rank 1 (1 df for the MANOVA test of the Intercept). Its projection outside the $\mathbf{E}$ ellipse shows a highly significant rejection of the hypothesis of equal means over Grade.

In such simple cases, traditional plots (boxplots, or plots of means with error bars) are easier to interpret. HE plots gain advantages in more complex designs (two or more between- or within-S factors, multiple responses), where they provide visual summaries of the information used in the multivariate hypothesis tests.

Between- and within-S effects {#sec:BW}

When there are both within- and between-S effects, the multivariate and univariate hypotheses tests can all be obtained together using Anova() with the idata and idesign specifying the within-S levels and the within-S design, as shown above. linearHypothesis() can be used to test arbitrary contrasts in the within- or between- effects.

However, to explain the visualization of these tests for within-S effects and their interactions using heplot() and related methods it is again convenient to explicitly transform $$\mathbf{Y}$ \mapsto $\mathbf{Y M}$$ for a given set of within-S contrasts, in the same way as done for the VocabGrowth data. See \@ref(sec:addendum) for simplified code producing these HE plots directly, without the need for explicit transformation.

To illustrate, we use the data from @OBrienKaiser:85 contained in the data frame OBrienKaiser in the car. The data are from an imaginary study in which 16 female and male subjects, who are divided into three treatments, are measured at a pretest, posttest, and a follow-up session; during each session, they are measured at five occasions at intervals of one hour. The design, therefore, has two between-subject and two within-subject factors.

For simplicity here, we initially collapse over the five occasions, and consider just the within-S effect of session, called session in the analysis below.

library("car")  # for OBrienKaiser data

# simplified analysis of OBrienKaiser data, collapsing over hour
OBK <- OBrienKaiser
OBK$pre <- rowMeans(OBK[,3:7])
OBK$post <- rowMeans(OBK[,8:12])
OBK$fup <- rowMeans(OBK[,13:17])
# remove separate hour scores
OBK <- OBK[,-(3:17)]

Note that the between-S design is unbalanced (so tests based on Type II sum of squares and crossproducts are preferred, because they conform to the principle of marginality).

table(OBK$gender, OBK$treatment)

The factors gender and treatment were specified with the following contrasts, $$\mathbf{L}${gender}$, and $$\mathbf{L}${treatment}$, shown below. The contrasts for treatment are nested (Helmert) contrasts testing a comparison of the control group with the average of treatments A and B (treatment1) and the difference between treatments A and B (treatment2).

contrasts(OBK$gender)
contrasts(OBK$treatment)

We first fit the general MANOVA model for the three repeated measures in relation to the between-S factors. As before, this just establishes the model for the between-S effects.

# MANOVA model
mod.OBK <- lm(cbind(pre, post, fup) ~  treatment*gender,  data=OBK)
mod.OBK

If we regarded the repeated measure effect of session as equally spaced, we could simply use polynomial contrasts to examine linear (slope) and quadratic (curvature) effects of time. Here, it makes more sense to use profile contrasts, testing (post - pre) and (fup - post).

session <- ordered(c("pretest", "posttest", "followup"),
    levels=c("pretest", "posttest", "followup"))
contrasts(session) <- matrix(c(-1,  1, 0,
                                0, -1, 1), ncol=2)
session
idata <- data.frame(session)

The multivariate tests for all between- and within- effects are then calculated as follows:

# Multivariate tests for repeated measures
aov.OBK <- Anova(mod.OBK, idata=idata, idesign=~session, test="Roy")
aov.OBK

It is useful to point out here that the default print methods for Anova.mlm objects, as shown above, gives an optimally compact summary for all between- and within-S effects, using a given test statistic, yet all details and other test statistics are available using the summary method.% \footnote{ In contrast, SAS proc glm and SPSS General Linear Model provide only more complete, but often bewildering outputs that still recall the days of Fortran coding in spite of more modern look and feel. } For example, using summary(aov.OBK) as shown below, we can display all the multivariate tests together with the $\mathbf{H}$ and $\mathbf{E}$ matrices, and/or all the univariate tests for the traditional univariate mixed model, under the assumption of sphericity and with Geiser-Greenhouse and Huhyn-Feldt corrected $F$ tests. To conserve space in this article the results are not shown here.

# All multivariate tests
summary(aov.OBK, univariate=FALSE)
# Univariate tests for repeated measures
summary(aov.OBK, multivariate=FALSE)

OK, now its time to understand the nature of these effects. Ordinarily, from a data-analytic point of view, I would show traditional plots of means and other measures (as in \Figure @ref(fig:voc-boxplot)) or their generalizations in effect plots [@Fox:87,effects2]. But I'm not going to do that here. Instead, I'd like for you to understand how HE plots provide a compact visual summary of an MLM, mirroring the tabular presentation from Anova(mod.OBK) above, but which also reveals the nature of effects. Here, you have to bear in mind that between-S effects are displayed most naturally in the space of the response variables, while within-S effects are most easily seen in the contrast space of transformed responses ($\mathbf{Y M}$).

HE plots for between-S effects can be displayed for any pair of responses with heplot(). \Figure @ref(fig:obk-HE1) shows this for pre and post. By default, $\mathbf{H}$ ellipses for all model terms (excluding the intercept) are displayed. Additional MLM tests can also be displayed using the hypotheses argument; here we specify the two contrasts for the treatment effect shown above as contrasts(OBK$treatment).

# HE plots for Between-S effects
heplot(mod.OBK, hypotheses=c("treatment1", "treatment2"),
  col=c("red", "black", "blue", "brown", "gray40", "gray40"),
  hyp.labels=c("(A,B)-Control", "A-B"),
  main="Between-S effects and contrasts"
    )
lines(c(3,7), c(3,7), col="green")

[htb] \centering \includegraphics[width=.65\textwidth]{fig/plot-obk-HE1} \caption{HE plot for the mod.OBK repeated measures model, showing between-S effects and contrasts in the space of the pre and post variables. Main effect means of the treatment (A, B, Control) and gender (M, F) groups are marked by points. Contrasts among the treatment groups appear as lines labeled at one extreme. The green line of unit slope shows equality of pre = post.} {#fig:obk-HE1}

In \Figure @ref(fig:obk-HE1), we see that the treatment effect is significant, and the large vertical extent of this $\mathbf{H}$ ellipse shows this is largely attributable to the differences among groups in the Post session. Moreover, the main component of the treatment effect is the contrast between the Control group and groups A \& B, which outperform the Control group at Post test. The effect of gender is not significant, but the HE plot shows that that males are higher than females at both Pre and Post tests. Likewise, the treatment:gender interaction fails significance, but the orientation of the $\mathbf{H}$ ellipse for this effect can be interpreted as showing that the differences among the treatment groups are larger for males than for females. Finally, the line of unit slope shows that for all effects, scores are greater on post than pre.

Using heplot3d(mod.OBK, ...) gives an interactive 3D version of \Figure @ref(fig:obk-HE1) for pre, post, and fup, that can be rotated and zoomed, or played as an animation.

heplot3d(mod.OBK, hypotheses=c("treatment1", "treatment2"),
    col=c("pink", "black", "blue", "brown", "gray40", "gray40"),
    hyp.labels=c("(A,B)-Control", "A-B")
    )
# rotate around y axis
play3d( rot8y <- spin3d(axis=c(0,1,0)),duration=12 )

This plot is not shown here, but an animated version can be generated from the code included in the supplementary materials.

Alternatively, all pairwise HE plots for the session means can be shown using pairs() for the mlm object mod.OBK, with the result shown in \Figure @ref(fig:obk-HE2).

pairs(mod.OBK, col=c("red", "black", "blue", "brown"))

[htb] \centering \includegraphics[width=.75\textwidth]{fig/plot-obk-HE2} \caption{HE plot for the mod.OBK repeated measures model, showing between-S effects for all pairs of sessions. The panel in row 2, column 1 is identical to that shown separately in \Figure @ref(fig:obk-HE1).} {#fig:obk-HE2}

Here we see that (a) the treatment effect is largest in the combination of post-test and follow up; (b) this 2 $df$ test is essentially 1-dimensional in this view, i.e., treatment means at post-test and follow up are nearly perfectly correlated; (c) the superior performance of males relative to females, while not significant, holds up over all three sessions.

As before, for expository purposes, HE plots for within-S effects are constructed by transforming $$\mathbf{Y}$ \mapsto $\mathbf{Y M}$$, here using the (profile) contrasts for session.

# Transform to profile contrasts for within-S effects
OBK$session.1 <- OBK$post - OBK$pre
OBK$session.2 <- OBK$fup - OBK$post
mod1.OBK <- lm(cbind(session.1, session.2) ~ treatment*gender,  data=OBK)
Anova(mod1.OBK, test="Roy", type="III")

From the schematic summary in \Table @ref(tab:BWtab), with these (or any other) contrasts as $$\mathbf{M}$_{session}$, the tests of the effects contained in treatment*gender in mod1.OBK are identical to the interactions of these terms with session, as shown above for the full model in aov.OBK. The (Intercept) term in this model represents the session effect.

The HE plot for within-S effects (\Figure @ref(fig:obk-HE3)) is constructed from the mod1.OBK object as shown below. The main manipulation done here is to re-label the terms plotted to show each of them as involving session, as just described. %\footnote{For technical reasons, it is necessary to calculate %the $\mathbf{H}$ matrices using Type III tests, and explicitly retain the %(Intercept) term which is ignored by default in heplot(). %}

# HE plots for Within-S effects
heplot(mod1.OBK,
  main="Within-S effects: Session * (Treat*Gender)",
  remove.intercept=FALSE, type="III",
  xlab="Post-Pre", ylab="Fup-Post",
  term.labels=c("session", "treatment:session", "gender:session",
                "treatment:gender:session"),
  col=c("red", "black", "blue", "brown"),
  xlim=c(-2,4), ylim=c(-2,3)
)
mark.H0(0,0)

[htb] \centering \includegraphics[width=.65\textwidth]{fig/plot-obk-HE3} \caption{HE plot for the mod1.OBK model, showing within-S effects in the space of contrasts among sessions. The point labeled $H_0$ here marks the comparison point for no difference over session in contrast space.} {#fig:obk-HE3}

\Figure @ref(fig:obk-HE3) provides an interpretation of the within-S effects shown in the MANOVA table shown above for Anova(mod.OBK). We can see that the effects of session and treatment:session are significant. More importantly, for both of these, but the interaction in particular, the significance of the effect is more attributable to the post-pre difference than to fup-post.

Higher-order designs

The scheme described above and the obvious generalization of \Table @ref(tab:BWtab) easily accommodates designs with two or more within-S factors. Any number of between-S factors are handled automatically, by the model formula for between-S effects specified in the lm() call, e.g., ~ treatment * gender.

For example, for the O'Brien-Kaiser data with session and hour as two within-S factors, first create a data frame, within specifying the values of these factors for the $3 \times 5$ combinations.

session <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
    levels=c("pretest", "posttest", "followup"))
contrasts(session) <- matrix(c(-1,  1, 0,
                                    0, -1, 1), ncol=2)
hour <- ordered(rep(1:5, 3))
within <- data.frame(session, hour)

The within data frame looks like this:

str(within)
within

The repeated measures MANOVA analysis can then be carried out as follows:

mod.OBK2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
                     post.1, post.2, post.3, post.4, post.5,
                     fup.1, fup.2, fup.3, fup.4, fup.5) ~  treatment*gender,
                data=OBrienKaiser)
(aov.OBK2 <- Anova(mod.OBK2, idata=within, idesign=~session*hour, test="Roy"))

Note that the test statistics for treatment, gender, session and all interactions among them are identical to what was found in the simplified analysis above. Among the effects including hour, only the main effect is significant here.

The following $\mathbf{M}$ matrices, corresponding to profile contrasts for session and polynomial contrasts for hour are used internally in Anova() in calculating these effects (shown here as integers, rather than in normalized form).

$\mathbf{M}$_{session} = \left( {rr} -1 & 0 \ 1 & -1 \ 0 & 1 \

\right) \quad\quad $\mathbf{M}$_{hour} = \left( {rrrr} -2 & 2 & -1 & 1 \ -1 & -1 & 2 & -4 \ 0 & -2 & 0 & 6 \ 1 & -1 & -2 & -4 \ 2 & 2 & 1 & 1 \

\right)

Tests involving the interaction of session:hour use the Kronecker product, $$\mathbf{M}${session} \otimes $\mathbf{M}${hour}$.

For HE plots, it is necessary to explicitly carry out the transformation of $$\mathbf{Y}$ \mapsto $\mathbf{Y}$ $\mathbf{M}$w$, where $ $\mathbf{M}$_w$ conforms to $\mathbf{Y}$ and represents the contrasts for the within-S effect. In the present example, this means that $$\mathbf{M}${session}$ and $$\mathbf{M}$_{hour}$ are both expanded as Kronecker products with the unit vector,

$\mathbf{M}$_s & = & $\mathbf{M}$_{session} \otimes \vec{1}_5 \comma \\
$\mathbf{M}$_h & = & \vec{1}_3 \otimes $\mathbf{M}$_{hour} \period  \\

These calculations in R are shown below:

M.session <- matrix(c(-1,  1, 0,
                       0, -1, 1), ncol=2)
rownames(M.session) <-c("pre", "post", "fup")       
colnames(M.session) <-paste("s", 1:2, sep="")       

M.hour <- matrix(c(-2, -1, 0, 1, 2,
                    2, -1, -2, -1, 1,
                    -1, 2, 0, -2, 1,
                    1, -4, 6, -4, 1), ncol=4)
rownames(M.hour) <- paste("hour", 1:5, sep="")  
colnames(M.hour) <- c("lin", "quad", "cubic", "4th")
M.hour

unit <- function(n, prefix="") {
    J <-matrix(rep(1, n), ncol=1)
    rownames(J) <- paste(prefix, 1:n, sep="")
    J
}

M.s <- kronecker( M.session, unit(5, "h"), make.dimnames=TRUE)

(M.h <- kronecker( unit(3, "s"), M.hour, make.dimnames=TRUE))

M.sh <- kronecker( M.session, M.hour, make.dimnames=TRUE)

Using M.h, we can construct the within-model for all terms involving the hour effect,

Y.hour <- as.matrix(OBrienKaiser[,3:17]) %*% M.h
mod.OBK2.hour <- lm(Y.hour ~  treatment*gender,  data=OBrienKaiser)

We can plot these effects for the linear and quadratic contrasts of hour, representing within-session slope and curvature. \Figure @ref(fig:obk2-HE1) is produced as shown below. As shown by the Anova(mod.OBK2, ...) above, all interactions with hour are small, and so these appear wholly contained within the $\mathbf{E}$ ellipse. In particular, there are no differences among groups (treatment $\times$ gender) in the slopes or curvatures over hour. For the main effect of hour, the linear effect is almost exactly zero, while the quadratic effect is huge.

labels <- c("hour", paste(c("treatment","gender","treatment:gender"),":hour", sep=""))
colors <- c("red", "black", "blue", "brown", "purple")
heplot(mod.OBK2.hour, type="III", remove.intercept=FALSE, term.labels=labels, col=colors)
mark.H0()

[htb] \centering \includegraphics[width=.6\textwidth]{fig/plot-obk2-HE1} \caption{HE plot for the mod.OBK2 repeated measures model, showing within-S effects for linear and quadratic contrasts in hour. As in \Figure @ref(fig:obk-HE3), we are viewing hypothesis and error variation in the transformed space of the repeated measures contrasts, here given by $$\mathbf{M}$_h$.} {#fig:obk2-HE1}

The pairs() function shows these effects (\Figure @ref(fig:obk2-HE2)) for all contrasts in hour. To reduce clutter, we only label the hour effect, since all interactions with hour are small and non-significant. The main additional message here is that the effects of hour are more complex than just the large quadratic component we saw in \Figure @ref(fig:obk2-HE1).

pairs(mod.OBK2.hour, type="III", remove.intercept=FALSE, term.labels="hour", col=colors)

[htb] \centering \includegraphics[width=.6\textwidth]{fig/plot-obk2-HE2} \caption{HE plot for the mod.OBK2 repeated measures model, showing within-S effects for all pairs of contrasts in hour.} {#fig:obk2-HE2}

Doubly-multivariate designs {#sec:doubly}

In the designs discussed above the same measure is observed on all occasions. Sometimes, there are two or more different measures, $Y_1, Y_2, \dots$, observed at each occasion, for example response speed and accuracy. In these cases, researchers often carry out separate repeated measures analyses for each measure. However the tests of between-S effects and each within-S effect can also be carried out as multivariate tests of $Y_1, Y_2, \dots$ simultaneously, and these tests are often more powerful, particularly when the effects for the different measures are weak, but correlated.

In the present context, such doubly-multivariate designs can be easily handled in principle by treating the multiple measures as an additional within-S factor, but using an identity matrix as the $\mathbf{M}$ matrix in forming the linear hypotheses to be tested via \Eqn. @ref(eq:mglt-rep). For example, with two measures, $Y_1, Y_2$ observed on three repeated sessions, the full $\mathbf{M}$ matrix for the design is generated as in \Eqn. @ref(eq:Mcd) as

{#eq:Mmc} $\mathbf{M}${CM} = (\vec{1}, $\mathbf{M}${session} ) \otimes $\mathbf{I}$_2 = \left( {rrr} 1 & -1 & 0 \ 1 & 1 & -1 \ 1 & 0 & 1 \

\right)
\otimes \left( {rr} 1 & 0 \ 0 & 1 \

\right) \period

In \R, we can express this as follows, using M.measure to represent $Y_1, Y_2$.

M.measure <- diag(2)
rownames(M.measure)<- c("Y1", "Y2")
colnames(M.measure)<- c("Y1", "Y2")

kronecker(cbind(1, M.session), M.measure, make.dimnames=TRUE)

In the result, the first two columns correspond to the within-S Intercept term, and are used to test all between-S terms for $Y_1, Y_2$ simultaneously. The remaining columns correspond to the session effect for both variables and all interactions with session. In practice, this analysis must be performed in stages because Anova() does not (yet)% \footnote{The new version of the car (2.0-0) on CRAN now includes enhanced Anova() and linearHypothesis() which perform these analyses.} allow such a doubly-multivariate design to be specified directly.

Example: Weight loss and self esteem

To illustrate, we consider the data frame WeightLoss originally from Andrew Ainsworth ([http://www.csun.edu/~ata20315/psy524/main.htm]), giving (contrived) data on weight loss and measures of self esteem after each of three months for 34 individuals, who were observed in one of three groups: Control, Diet group, Diet $+$ Exercise group. The within-S factors are thus measure (wl, se) and month (1:3).

table(WeightLoss$group)
some(WeightLoss)

Because this design is complex, and to facilitate interpretation of the effects we will see in HE plots, it is helpful to view traditional plots of means with standard errors for both variables. These plots, shown in \Figure @ref(fig:wl-means),% \footnote{ These plots were drawn using plotmeans() in the gplots [@gplots]. The code is not shown, but is available in the R example code for this paper. } show that, for all three groups, the amount of weight lost each month declines, but only the Diet $+$ Exercise maintains substantial weight loss through month 2. For self esteem, all three groups have a U-shaped pattern over months, and by month 3, the groups are ordered Control $<$ Diet $<$ Diet $+$ Exercise.

op <- par(mfrow=c(1,2))
library("gplots")
WLlong <-reshape(WeightLoss, varying=list(WeightLoss=2:4, SelfEsteem=5:7), direction="long",
    timevar = "month",
    v.names=c("WeightLoss", "SelfEsteem"))
WLlong <- WLlong[ order(WLlong$id),]

plotmeans(WeightLoss ~ interaction(month, group), data=WLlong,
    ylab="Weight Loss", xlab="Month / Group",
    connect=list(1:3, 4:6, 7:9), pch=7, p=0.68,
    n.label=FALSE, main="Weight Loss: Group x Month",
    legends =  paste(rep(1:3, 3),c("", "Control", "", "", "Diet", "", "", "DietEx", ""),sep="\n")
    )
abline(v=c(3.5, 6.5), col="gray")

plotmeans(SelfEsteem ~ interaction(month, group), data=WLlong,
    ylab="Self Esteem", xlab="Month / Group",
    connect=list(1:3, 4:6, 7:9), pch=7, p=0.68,
    n.label=FALSE, main="Self Esteem: Group x Month",
    legends =  paste(rep(1:3, 3),c("", "Control", "", "", "Diet", "", "", "DietEx", ""),sep="\n")
    )
abline(v=c(3.5, 6.5), col="gray")
par(op)

[htb] \centering \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-wl-means} \caption{Means for weight loss and self esteem by group and month. Error bars show $\pm 1$ standard error for each mean.} {#fig:wl-means}

Research interest in the differences among groups would likely be focused on the questions: (a) Do the two diet groups differ from the control group? (b) Is there an additional effect of exercise, given diet? These questions may be tested with the (Helmert) contrasts used below for group, which are labeled group1 and group1 respectively.

contrasts(WeightLoss$group) <- matrix(c(-2,1,1, 0, -1, 1),ncol=2)
(wl.mod<-lm(cbind(wl1,wl2,wl3,se1,se2,se3)~group, data=WeightLoss))

A standard between-S MANOVA, ignoring the within-S structure shows a highly significant group effect.

Anova(wl.mod)

As before, it is often useful to examine HE plots for pairs of variables in this analysis before proceeding to the within-S analysis. For example, \Figure @ref(fig:wl-HE1) shows the test of group and the two contrasts for weight loss and for self esteem at months 1 and 2. %r %heplot(wl.mod, hypotheses=c("group1", "group2"), % xlab="Weight Loss, month 1", ylab="Weight Loss, month 2") %

op <- par(mfrow=c(1,2))
heplot(wl.mod, hypotheses=c("group1", "group2"),
    xlab="Weight Loss, month 1", ylab="Weight Loss, month 2")
heplot(wl.mod, hypotheses=c("group1", "group2"), variables=4:5,
    xlab="Self Esteem, month 1", ylab="Self Esteem, month 2")
par(op)

[htb] \centering \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-wl-HE1} \caption{HE plot for the wl.mod MANOVA model, showing between-S effects for weight loss (left) and self esteem (right) at months 1 and 2.} {#fig:wl-HE1}

This is helpful, but doesn't illuminate the overall group effect for weight loss and self esteem for all three months, and, of course cannot shed light on any interactions of group with measure or month. In the following discussion, we will assume that the researcher is particularly interested in understanding the relation between weight loss and self esteem as it is expressed in changes over time and differences among groups. % pairs(wl.mod, variables=1:3)

To carry out the doubly-multivariate analysis, we proceed as follows. First, we define the $\mathbf{M}$ matrix for the measures, used in the between-S analysis. We use $$\mathbf{M}$ = $\mathbf{I}$_2 \otimes \vec{1}/3$ so that the resulting scores are the means (not sums) for weight loss and self esteem.

measure <- kronecker(diag(2), unit(3, 'M')/3, make.dimnames=TRUE)
colnames(measure)<- c('WL', 'SE')
measure

between <- as.matrix(WeightLoss[,-1]) %*% measure

between.mod <- lm(between ~ group, data=WeightLoss)
Anova(between.mod, test="Roy", type="III")

The HE plot for this component of the analysis (\Figure @ref(fig:wl-HE2)) shows a striking effect: Averaging over all three months, the means for the Control, Diet and DietEx group on both weight loss and self esteem are highly correlated and in the expected direction. This is something that is not at all obvious in \Figure @ref(fig:wl-means).

heplot(between.mod, hypotheses=c("group1", "group2"),
    xlab="Weight Loss", ylab="Self Esteem",
    col=c("red", "blue", "brown"),
    main="Weight Loss & Self Esteem: Group Effect")

[htb] \centering \includegraphics[width=.6\textwidth]{fig/plot-wl-HE2} \caption{HE plot for the between.mod doubly-multivariate design, showing overall between-S effects for weight loss and self esteem.} {#fig:wl-HE2}

Next, for the within-S analysis, we define the $\mathbf{M}$ matrix for months, using orthogonal polynomials representing linear and quadratic trends. As before, the test of the (Intercept) term in these trend scores corresponds to the month effect in the doubly-multivariate model, and the group effect tests the group $\times$ month interaction.

month <- kronecker(diag(2), poly(1:3, degree=2), make.dimnames=TRUE)
colnames(month)<- c('WL1', 'WL2', 'SE1', 'SE2')
round(month,digits=4)
trends <- as.matrix(WeightLoss[,-1]) %*% month
within.mod <- lm(trends ~ group, data=WeightLoss)
Anova(within.mod, test="Roy", type="III")

HE plots corresponding to this model (\Figure @ref(fig:wl-HE3)) can be produced as follows. The $\mathbf{H}$ and $\mathbf{E}$ matrices are all $4 \times 4$, but the $\mathbf{H}$ matrices for the month and group:month effects are rank 1 and 2 respectively.

op <- par(mfrow=c(1,2))
heplot(within.mod, hypotheses=c("group1", "group2"), variables=c(1,3),
    xlab="Weight Loss - Linear", ylab="Self Esteem - Linear",
    type="III", remove.intercept=FALSE,
    term.labels=c("month", "group:month"),
    main="(a) Within-S Linear Effects")
mark.H0()

heplot(within.mod, hypotheses=c("group1", "group2"), variables=c(2,4),
    xlab="Weight Loss - Quadratic", ylab="Self Esteem - Quadratic",
    type="III", remove.intercept=FALSE,
    term.labels=c("month", "group:month"),
    main="(b) Within-S Quadratic Effects")
mark.H0()
par(op)

%heplot(within.mod, hypotheses=c("group1", "group2"), % xlab="Weight Loss", ylab="Self Esteem", % type="III", remove.intercept=FALSE, % term.labels=c("month", "group:month"), % main="Weight Loss & Self Esteem: Within-S Effects") %mark.H0()

[htb] \centering \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-wl-HE3} \caption{HE plots for the within.mod doubly-multivariate design, showing the effects of month and the interaction group:month for weight loss vs.\ self esteem. (a) Linear effects; (b) Quadratic effects.} {#fig:wl-HE3}

\Figure @ref(fig:wl-HE3) shows the plots for the linear and quadratic effects separately for weight loss vs.\ self esteem. The plot of linear effects (\Figure @ref(fig:wl-HE3)(a)) shows that the effect of month can be be described as negative slopes for weight loss combined with positive slopes for self esteem-- all groups lose progressively less weight over time, but generally feel better about themselves. Differences among groups in the group:month effect are in the same direction, but with greater differences among groups in the slopes for self esteem. The interpretation of the quadratic effects (\Figure @ref(fig:wl-HE3)(b)) is similar, except here, differences in curvature over months are driven largely by the difference between the DietEx group from the others on weight loss.

The interested reader might wish to compare the standard univariate plots of means in \Figure @ref(fig:wl-means) with the HE plots in \Figure @ref(fig:wl-HE2) and \Figure @ref(fig:wl-HE3). The univariate plots have the advantage of showing the data directly, but cannot show the sources of significant effects in multivariate repeated measures models. HE plots have the advantage that they show directly what is expressed in the multivariate tests for relevant hypotheses.

Simplified interface: heplots 0.9 and car 2.0 {#sec:addendum}

It sometimes happens that the act of describing and illustrating software spurs development to make both simpler, and such is the case here. At the beginning, the stable version of car on CRAN provided the computation for multivariate linear hypotheses including repeated measures designs, but could not handle doubly-multivariate designs directly; the CRAN version of heplots could only repeated measures by explicitly transforming $$\mathbf{Y}$ \mapsto $\mathbf{Y M}$$ and re-fitting submodels in terms of the transformed responses.

The new versions of these packages on CRAN ([http://cran.us.r-project.org/]) now handle these cases directly from the basic mlm object. heplot() now provides the arguments idata, idesign, icontrasts, or, for the doubly-multivariate case, imatrix, which are passed to Anova() to calculate the appropriate $\mathbf{H}$ and $\mathbf{E}$ matrices.

Omitting these arguments in the call to heplot() gives an HE plot for all between-S effects (or the subset specified by the terms argument), just as before. For the within-S effects, $\mathbf{E}$ matrices differ for for different within-S terms, so it is necessary to specify the intra-subject term (iterm, corresponding to $\mathbf{M}$) for which HE plots are desired. Several examples are given below.

For the VocabularyGrowth data, \Figure @ref(fig:voc8)(b) can be produced by

(Vocab.mod <- lm(cbind(grade8,grade9,grade10,grade11) ~ 1, data=VocabGrowth))

idata <-data.frame(grade=ordered(8:11))

heplot(Vocab.mod, type="III", idata=idata, idesign=~grade, iterm="grade",
    main="HE plot for Grade effect")

For the OBrienKaiser data, the code for plots of between-S effects is the same as shown above for \Figure @ref(fig:obk-HE1) and \Figure @ref(fig:obk-HE2). The HE plot for within-S effects involving session (\Figure @ref(fig:obk-HE3)) can be produced using iterm="session":

mod.OBK <- lm(cbind(pre, post, fup) ~  treatment*gender,  data=OBK)
session <- ordered(c("pretest", "posttest", "followup"),
    levels=c("pretest", "posttest", "followup"))
contrasts(session) <- matrix(c(-1,  1, 0,
                                0, -1, 1), ncol=2)

idata <- data.frame(session)
heplot(mod.OBK, idata=idata, idesign=~session, iterm="session",
  col=c("red", "black", "blue", "brown"),
  main="Within-S effects: Session * (Treat*Gender)")

Similarly, HE plots for terms involving hour can be obtained using the expanded model (mod.OBK2) for the 15 combinations of hour and session:

mod.OBK2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
                     post.1, post.2, post.3, post.4, post.5,
                     fup.1, fup.2, fup.3, fup.4, fup.5) ~  treatment*gender,
                data=OBrienKaiser)
heplot(mod.OBK2, idata=within, idesign=~hour, iterm="hour")
heplot(mod.OBK2, idata=within, idesign=~session*hour, iterm="session:hour")

%\SweaveInput{compare.Rnw} %\SweaveOpts{engine=R,height=6,width=6,results=hide,fig=FALSE,echo=TRUE} %\SweaveOpts{prefix.string=fig/plot}

Comparison with other approaches {#sec:compare}

The principal goals of this paper have been (a) to describe the extension of the classical \MLM\ to repeated measures designs; (b) to explain how HE plots provide compact and understandable visual summaries of the effects shown in typical numerical tables of MANOVA tests; and (c) illustrate these in a variety of contexts ranging from single-sample designs to complex doubly-multivariate designs.

In the context of repeated measures designs, I mentioned earlier that mixed models for longitudinal data provide an attractive alternative to the \MLM\ (because the former easily accommodate missing or unbalanced data over intra-subject measurements, time-varying covariates, and often allows the residual covariation to be modeled with fewer parameters).

Here we consider a classic data set [@PotthoffRoy:64] used in the first application of the \MLM\ to growth-curve analysis. These data are often used as illustrations of longitudinal models, e.g., \citet[\S 17.4]{VerbekeMolenberghs:00}.

Investigators at the University of North Carolina Dental School followed the growth of 27 children (16 males, 11 females) from age 8 until age 14 in a study designed to establish typical patterns of jaw size useful for orthodontic practice. Every two years they measured the distance between the pituitary and the pterygomaxillary fissure, two points that are easily identified on x-ray exposures of the side of the head. The questions of interest include (a) Over this age range, can growth be adequately represented as linear in time, or is some more complex function necessary? (b) Are separate growth curves needed for boys and girls, or can both be described by the same growth curve?

Longitudinal, mixed model approach

The mixed model for longitudinal data is very general and flexible for the reasons noted above, but it is inappropriate here to relate any more than the barest of details necessary for this example. We begin with simple plots of the data: A profile plot grouped by Sex (\Figure @ref(fig:ortho-xyplot) (left)),

data("Orthodont", package="nlme")
library("lattice")
xyplot(distance ~ age|Sex, data=Orthodont, type='b', groups=Subject, pch=15:25,
    col=palette(), cex=1.3, main="Orthodont data")

print(trellis.last.object())

%r %print( %xyplot(distance ~ age|Sex, data=Orthodont, type='b', groups=Subject, pch=15:25, % col=palette(), cex=1.3, main="Orthodont data") %) %

and also a summary plot showing fitted lines for each individual, together with the pooled ordinary least squares regression of distance on age (\Figure @ref(fig:ortho-xyplot) (right)).

xyplot(distance ~ age | Sex, data = Orthodont, groups = Subject,
    main = "Pooled OLS and Individual linear regressions ~ age", type = c('g', 'r'),
    panel = function(x, y, ...) {
        panel.xyplot(x, y, ..., col = gray(0.5))
        panel.lmline(x, y, ..., lwd = 3, col = 'red')
                })

plot(trellis.last.object())

[htb] %\centering [c]{.49\textwidth} \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-xyplot1}

[c]{.49\textwidth} \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-xyplot2}

\caption{Profile plot of Orthodont data, by sex (left); Pooled OLS and individual linear regressions on age, by sex (right)} {#fig:ortho-xyplot}

From these plots, we can see that boys generally have larger jaw distances than girls, and the rate of growth (slopes) for boys is generally larger than for girls. It is difficult to discern any patterns within the sexes, except that one boy seems to stand out, with a lower intercept and steeper slope.

With the longitudinal mixed model, contemplate fitting two models describing an individual's pattern of growth: a model fitting only linear growth and a model fitting each person's trajectory exactly by including quadratic and cubic trends in time. For the sake of interpretation of coefficients in these models, it is common to recenter the time variable so that time=0 corresponds to initial status. Using Year = (age-8), we have:

\mathrm{m1:} \quad\quad y_{it} & = & \beta_{0i} + \beta_{1i} \mathrm{Year}{it} + e{it} \ \mathrm{m3:} \quad\quad y_{it} & = & \beta_{0i} + \beta_{1i} \mathrm{Year}{it} + \beta{2i} \mathrm{Year}^2_{it} + \beta_{3i} \mathrm{Year}^3_{it}+ e_{it} \comma

where the vector of residuals for subject $i$ is $\vec{e}{i \bullet} \sim N ( \vec{0}, $\mathbf{R}$_i )$. (For this example, we take $$\mathbf{R}$_i$ to be unstructured, even though other specifications require fewer parameters.) For the linear model (m1), we entertain the possibility that the person-level intercepts ($\beta{0i}$) and slopes ($\beta_{1i}$) depend on Sex, and so specify them as random coefficients,

\beta_{0i} & = & \gamma_{00} + \gamma_{01} \mathrm{Sex}{i} + u{0i} \comma \ \beta_{1i} & = & \gamma_{10} + \gamma_{11} \mathrm{Sex}{i} + u{1i} \period

In these equations the $\gamma$s are the fixed effects, while the $u$ (along with the errors $e_{it}$) are random effects. Note that Sex is coded 0=Male, 1=Female, so $\gamma_{00}$ and $\gamma_{10}$ are the intercept and slope for Males; $\gamma_{01}$ pertains to the difference in intercepts for Females relative to Males, while $\gamma_{11}$ is the difference in slopes.

The linear growth model can be fit using lme as follows:

Ortho <- Orthodont
Ortho$year <- Ortho$age - 8  # make intercept = initial status

Ortho.mix1 <- lme(distance ~ year * Sex, data=Ortho,
        random = ~ 1 + year | Subject, method="ML")
#Ortho.mix1
anova(Ortho.mix1)

Similarly, the model (m3) allowing cubic growth at level 1 can be fit using:

Ortho.mix3 <- lme(distance ~ year*Sex + I(year^2) + I(year^3), data=Ortho,
        random = ~ 1 + year | Subject, method="ML")
anova(Ortho.mix3)

A likelihood ratio test confirms that the quadratic and cubic components of Year do not improve the model,

anova(Ortho.mix1, Ortho.mix3)

To aid interpretation, we can plot the estimated fixed effects from the linear model (m1) as follows, using the predict method for lme objects to calculate the the fitted values for boys and girls over the range of years (0 to 6) corresponding to ages 8 to 14, as in \Figure @ref(fig:Ortho-fm)(left). Similar code produces a plot of the cubic model \Figure @ref(fig:Ortho-fm)(right).

grid <- expand.grid(year=0:6, Sex=c("Male", "Female"))
grid$age <- grid$year+8  # plot vs. age
fm.mix1 <-cbind(grid, distance = predict(Ortho.mix1, newdata = grid, level=0))

xyplot(distance ~ age, data=fm.mix1, groups=Sex, type="b",
    par.settings = list(superpose.symbol = list(cex = 1.2, pch=c(15,16))),
    auto.key=list(text=levels(fm.mix1$Sex), points = TRUE, x=0.05, y=0.9, corner=c(0,1)),
    main="Linear mixed model: predicted growth")

plot(trellis.last.object())

fm.mix3 <-cbind(grid, distance = predict(Ortho.mix3, newdata = grid, level=0))
print(xyplot(distance ~ age, data=fm.mix3, groups=Sex, type="b",
    par.settings = list(superpose.symbol = list(cex = 1.2, pch=c(15,16))),
    auto.key=list(text=levels(fm.mix3$Sex), points = TRUE, x=0.05, y=0.9, corner=c(0,1)),
    main="Cubic mixed model: predicted growth"))

[htb] %\centering [c]{.49\textwidth} \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-Ortho-fm1}

[c]{.49\textwidth} \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-Ortho-fm3}

\caption{Fitted values showing the estimated fixed effects of age, Sex, and their interaction in the linear growth model (left) and cubic growth model (right).} {#fig:Ortho-fm}

For this simple example with only two predictors, such plots provide a direct visual summary of the fitted fixed effects in the model, as far as they go.

MVLM approach

For the multivariate approach, the Orthodont data must first be reshaped to wide format with the distance values as separate columns.

library("nlme")
Orthowide <- reshape(Orthodont, v.names="distance", idvar=c("Subject", "Sex"),
    timevar="age", direction="wide")
some(Orthowide, 4)

The \MLM\ is then fit as follows, with Sex as the between-S factor. Age is quantitative, so the intra-subject data frame (idata) is created with age as an ordered factor.

Ortho.mod <- lm(cbind(distance.8, distance.10, distance.12, distance.14) ~ Sex, data=Orthowide)

idata <- data.frame(age=ordered(seq(8,14,2)))
Ortho.aov <- Anova(Ortho.mod, idata=idata, idesign=~age)
Ortho.aov

We see that both Sex and age are highly significant and their interaction is nearly significant.

\Figure @ref(fig:ortho-HE) shows HE plots for the between- and within-S effects, produced as shown below. The left panel plots the effect of Sex for ages 8 and 14, with a green line of unit slope. Males clearly show greater growth by age 14, and the difference between males and females is greater at at 14 than at age 8. The right panel shows the linear and quadratic trends with age, reflecting the overall age main effect and the Sex:age interaction. Recalling that the contributions of each displayed variable to each effect in an HE plot can be seen by their horizontal and vertical shadows relative to the $\mathbf{E}$ ellipse, we see that the main effect of age is essentially linear, and the overall Sex:age effect is nearly significant due to a difference in slopes, but not curvature.

op <- par(mfrow=c(1,2))
heplot(Ortho.mod, variables=c(1,4), asp=1, col=c("red", "blue"),
    xlim=c(18,30), ylim=c(18,30),
    main="Orthodont data: Sex effect")
abline(0,1, col="green")

heplot(Ortho.mod, idata=idata, idesign=~age, iterm="age", col=c("red", "blue", "brown"),
    main="Orthodont data: Within-S effects")
par(op)

[htb] \centering \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-HE} \caption{HE plots for the Ortho.mod MANOVA model, showing between-S effects for distance at ages 8 and 14 (left) and the linear and quadratic within-S effects for age and Sex:age (right).} {#fig:ortho-HE}

To examine the questions of interest here in more detail, we focus on the intra-subject design and ask if linear growth is sufficient to explain both average development over time and differences between boys and girls. This is easily answered visually from the pairs() plot (not shown here),

pairs(Ortho.mod, idata=idata, idesign=~age, iterm="age", col=c("red", "blue", "brown"))

%[htb] %\centering %\includegraphics[width=.8\textwidth, trim=0 30 0 30]{fig/plot-ortho-pairs} %\caption{} % {#fig:ortho-pairs} % and, in particular, the panel corresponding to the nonlinear (quadratic and cubic) components of trend, shown in \Figure @ref(fig:ortho-nonlin-HE).

heplot(Ortho.mod, idata=idata, idesign=~age, iterm="age", col=c("red", "blue", "brown"),
    variables=c(2,3), main="Orthodont data: Nonlinear Within-S effects")

[htb] \centering \includegraphics[width=.5\textwidth, trim=0 30 0 30]{fig/plot-ortho-nonlin-HE} \caption{Nonlinear components of the within-S effects of age and Sex:age, showing that they are quite small in relation to within-S error.} {#fig:ortho-nonlin-HE}

We can confirm the impression that no nonlinear effects are important by testing linear hypotheses. To explain this, we first show the details of the test of the overall Sex:age effect, as tested with linearHypothesis(). The ``response transformation matrix'' shown below is equivalent to $$\mathbf{M}$_{poly}$ described earlier (\@ref(sec:within)) for a 4-level factor with linear, quadratic and cubic trend components. The univariate tests for individual contrasts in age are then based on the diagonal elements of the $\mathbf{H}$ and $\mathbf{E}$ matrices.

linearHypothesis(Ortho.mod, hypothesis="SexFemale", idata=idata, idesign=~age, iterms="age", title="Sex:age effect")

From this, the linear and nonlinear terms can be tested by selecting the appropriate columns of $$\mathbf{M}$_{poly}$ supplied as the contrasts associated with the age effect. For example, for tests of the linear effect of age and the Sex:age interaction (differences in slopes),

linear <- idata
contrasts(linear$age, 1) <- contrasts(linear$age)[,1]
print(linearHypothesis(Ortho.mod, hypothesis="(Intercept)",
    idata=linear, idesign=~age, iterms="age", title="Linear age"), SSP=FALSE)
print(linearHypothesis(Ortho.mod, hypothesis="SexFemale",
    idata=linear, idesign=~age, iterms="age", title="Linear Sex:age"), SSP=FALSE)

Similarly, the nonlinear effects of age and the Sex:age interaction can be tested as follows, using the contrasts for quadratic and cubic trends in age.

nonlin <- idata
contrasts(nonlin$age, 2) <- contrasts(nonlin$age)[,2:3]
print(linearHypothesis(Ortho.mod, hypothesis="(Intercept)",
    idata=nonlin, idesign=~age, iterms="age", title="Nonlinear age"), SSP=FALSE)
print(linearHypothesis(Ortho.mod, hypothesis="SexFemale",
    idata=nonlin, idesign=~age, iterms="age", title="Nonlinear Sex:age"), SSP=FALSE)

These examples show that, in simple cases, traditional plots of fitted values from mixed models (\Figure @ref(fig:Ortho-fm)) have the advantage of simple visual interpretation in terms of slopes and intercepts, at least for linear models. But such profile plots are generic: comparable plots could equally well be drawn for the fitted values from the \MLM\ in this section. HE plots for repeated measure designs have the additional advantage of showing the nature of significance tests and linear hypotheses, though the structure of the \MLM\ requires separate plots of between-S and within-S effects. In more complex designs with multiple between-S and within-S effects, and designs with multivariate responses (where mixed models do not apply), HE plots gain greater advantage.

Discussion

Graphical methods for univariate response models are well-developed, but analogous methods for multivariate responses are still developing. Indeed, this is a fruitful area for new research [@Friendly:07:manova]. HE plots provide one new direction, providing direct visualizations of effects in \MLM s, in the space of response variables (or in the the reduced-rank canonical space (candisc) displaying maximal differences).

This paper has shown how these methods can be extended to repeated measures designs, by displaying effects in the transformed space of contrasts or linear combinations for within-S effects. As we hope to have shown, these plots can provide insights into the relations among effects and interpretations of those that are significant (or not) which go beyond what is available in traditional, univariate displays.

Acknowledgments

I am grateful to John Fox for collaboration, advice and continued development of the car which provides the computational infrastructure for the heplots. The Associate Editor and one reviewer helped considerably in clarifying the presentation of these methods.

%# References

References

%\include{notes}

options(old.opts )

friendly/heplots documentation built on May 4, 2024, 12:05 a.m.