% 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}
%\bibliographystyle{abbrvnat}
%\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*{\given}{\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},\ldots,{#1}_{#2}}}
\newcommand*{\Var}{\ensuremath{\mathsf{Var}}}

\newcommand{\sizedmat}[2]{%
  \mathord{\mathop{\mat{#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]{{\texttt{#1}}}
\newcommand{\func}[1]{{\texttt{#1}()}}

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


\title{HE Plots for Repeated Measures Designs}
\author{Michael Friendly \\ York University, Toronto}
%\date{For \Rpackage{heplots} version \Sexpr{packageDescription("heplots")[["Version"]]} , \Sexpr{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 \citet{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 \pkg{heplots} for \proglang{R} \citep{Fox-etal:07:heplots} and \proglang{SAS}
	\citep{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 \pkg{heplots} for \proglang{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: \url{http://datavis.ca}
}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/1/31336-5053
%% Fax: +43/1/31336-734

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


\begin{document}



%\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}

<<preliminaries,echo=FALSE,results=hide>>=
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

@

\section{Introduction}
	Hypothesis Error (HE) plots, introduced in \citet{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 \pkg{heplots} \citep{Fox-etal:07:heplots} for \proglang{R} \citep{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 \emph{Journal of Statistical Software}
	\citep{Friendly:10:JSS}.

\subsection{Multivariate linear models: Notation}
To set notation, we express the \MLM\ as
\begin{equation}\label{eq:mglm}
\sizedmat{Y}{n \times p} = \sizedmat{X}{n \times q} \sizedmat{\Beta}{q \times p} +
\sizedmat{U}{n \times p}
\comma
\end{equation}
where, $\mat{Y} \equiv (\nvec{\vec{y}}{p})$ is the matrix of responses for $n$ subjects on $p$ variables,
$\mat{X}$ is the design matrix for $q$ regressors, $\mat{\Beta}$ is the $q \times p$ matrix
of regression coefficients or model parameters and $\mat{U}$ is the $n \times p$
matrix of errors,
with $\textrm{vec}(\mat{U}) \sim \mathcal{N}_p ( \mat{0}, \mat{I}_n \otimes \mat{\Sigma} )$,
where $\otimes$ is the Kronecker product.

A convenient feature of the \MLM\ for general multivariate responses is that
\emph{all} tests of linear hypotheses (for null effects) can be represented in the form of a general
linear test,
\begin{equation}\label{eq:mglt}
H_0 : \sizedmat{L}{h \times q}
\sizedmat{\Beta}{q \times p} =
\sizedmat{0}{h \times p}
\comma
\end{equation}
where $\mat{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 \proglang{R} all such tests can be carried out using
the functions \func{Anova} and \func{linearHypothesis} in the
\pkg{car}.%
\footnote{
Both the \pkg{car} and the \pkg{heplots} are being actively developed.
Except where noted, all results in this paper were produced using the old-stable versions on CRAN,
\textbf{car} 1.2-16 (2009/10/10)
% \Sexpr{packageDescription("car")[["Version"]]} (\Sexpr{(packageDescription("car")[["Date"]])})
and
\textbf{heplots} 0.8-11 (2009-12-08)
% \Sexpr{packageDescription("heplots")[["Version"]]}  (\Sexpr{(packageDescription("heplots")[["Date"]])})
running under \Sexpr{R.version$version.string}.
}

For \emph{any} such hypothesis of the form \eqref{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}:
\begin{equation} \label{eq:hmat}
\mat{H}  \equiv \mat{SSP}_H =
 (\mat{L} \widehat{\mat{B}})\trans \,
 [\mat{L} (\mat{X}\trans \mat{X} )^{-} \mat{L}\trans]^{-1} \,
 (\mat{L} \widehat{\mat{B}})
 \comma
\end{equation}
and
\begin{equation} \label{eq:emat}
\mat{E}  \equiv \mat{SSP}_E =
 \mat{Y}\trans \mat{Y} -
 \widehat{\mat{B}}\trans (\mat{X}\trans \mat{X}) \widehat{\mat{B}}
 =
  \widehat{\mat{U}}\trans  \widehat{\mat{U}}
 \comma
\end{equation}
where $\widehat{\mat{U}} = \mat{Y} - \mat{X} \widehat{\mat{B}}$ is the matrix of residuals.
Multivariate test statistics (Wilks' $\Lambda$, Pillai trace, Hotelling-Lawley trace, Roy's maximum root)
for testing \eqref{eq:mglt} are based on the $s = \min(p, h)$ non-zero latent roots of
$\mat{H} \mat{E}^{-1}$ and attempt to capture how ``large'' $\mat{H}$ is, relative to
$\mat{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 \eqref{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 (\mat{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.

\subsection{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 \emph{data ellipse} (or covariance ellipse),
 described by \citet{Dempster:69} and
\citet{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\mat{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}$:%
\begin{equation}
\mathcal{E}_{c}(\mathbf{y;S,}\overline{\vec{y}})\equiv\left\{ y\text{: }(%
\vec{y}-\overline{\vec{y}})\trans\mat{S}^{-1}(\vec{y}-\overline{%
\vec{y}})\leq c^{2}\right\}
\label{eqn-data-ellipse}
\end{equation}
Here, $\mat{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
\mat{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.
}


\subsection{HE plots}
The essential idea behind HE plots is that any multivariate hypothesis
test \eqref{eq:mglt} can be represented visually by ellipses (or ellipsoids in 3D) which express
the size  of co-variation against a multivariate null hypothesis
($\mat{H}$) relative to error covariation ($\mat{E}$).
The multivariate tests, based on the latent roots of $\mat{H} \mat{E}^{-1}$,
are thus translated directly to the sizes of the \mat{H} ellipses for
various hypotheses, relative to the size of the \mat{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 \mat{E} matrix is first scaled to a covariance matrix
($\mat{E}/df_e = \widehat{\Var}(\mat{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 \mat{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 \mat{E} ellipse, the \mat{H} ellipses
can be scaled to show either the \emph{effect size} or strength of
\emph{evidence} against $H_0$ (significance).

For effect size scaling, each \mat{H} is divided by $df_e$ to conform
to \mat{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 \mat{H} ellipses
relative to that of the \mat{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 \mat{H} ellipse is scaled to $\mat{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 \mat{H} ellipse extends anywhere outside the \mat{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 \mat{H} and \mat{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 $\mat{H} = \min(p, h) \le 2$, all significant effects can be
observed directly in 2D HE plots; when rank$(\mat{H})=3$, some rotation
of a 3D plot will
reveal each significant effect as extending somewhere outside the \mat{E}
ellipsoid.

In our \proglang{R} implementation, the basic plotting functions in the \pkg{heplots} are \func{heplot}
and \func{heplot3d} for \code{mlm} objects.
These rely heavily on the \func{Anova} and other functions from the \pkg{car}
\citep{car} for computation.
For more than three response variables, all pairwise HE plots can be
shown using a \func{pairs} function for \code{mlm} objects. Alternatively, the related
\pkg{candisc} \citep{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 \citet{Friendly:07:manova,FoxFriendlyMonette:09:compstat} for details
and examples for between-S MANOVA designs, MMREG and MANCOVA models.

\section{Repeated measures designs}\label{sec:repmeas}

The framework for the \MLM\ described above
pertains to the situation in which the response vectors (rows,
$\vec{y}_i\trans$ of $\mat{Y}_{n \times p})$ are \emph{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 \eqref{eq:mglm} applies, but hypotheses
about between- and within-individual variation are tested by an extended form
of the general linear test \eqref{eq:mglt}, which becomes
\begin{equation}\label{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
\end{equation}
where $\mat{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 \mat{H} and \mat{E} matrices for testing \eqref{eq:mglt-rep} become
\begin{equation} \label{eq:hmat-rep}
\mat{H}  =
 (\mat{L} \widehat{\mat{B}} \mat{M})\trans \,
 [\mat{L} (\mat{X}\trans \mat{X} )^{-} \mat{L}\trans]^{-1} \,
 (\mat{L} \widehat{\mat{B}} \mat{M})
 \comma
\end{equation}
and
\begin{equation} \label{eq:emat-rep}
\mat{E}  =
 (\mat{Y} \mat{M})\trans
 [\mat{I} - (\mat{X}\trans \mat{X} )^{-}  \mat{X}\trans] \,
 (\mat{Y} \mat{M})
 \period
\end{equation}
This may be easily seen to be just the ordinary \MLM\ applied to the
transformed responses $\mat{Y} \mat{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
\citet{Hsu:38} and is discussed further by \citet{Rencher:95} and \citet{Timm:80}.
In what follows, we refer to hypotheses pertaining to between-individual
effects (specified by \mat{L}) as ``between-S'' and hypotheses pertaining
to within-individual effects (\mat{M}) as ``within-S.''

\begin{table}[htb]
\renewcommand{\arraystretch}{1.5}
\renewcommand{\tabcolsep}{0.5cm}
%	\vspace{6pt}
	\label{tab:BWtab}
\setlength{\extrarowheight}{4pt}
\begin{tabular}{|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}
  $\mat{M}$ for Within-S effects & Intercept & $\mat{L} = \mat{L}_A$ & $\mat{L} = \mat{L}_B$ & $\mat{L} = \mat{L}_{AB}$ \\[1ex] \hline
  $\mat{M} = \mat{M}_1 = \left(
                     \begin{array}{rrr}
                       1 & 1 & 1 \\
                     \end{array}
                    \right)\trans
  $  & $\vec{\mu}_{..}$ & \code{A} & \code{B} & \code{A:B}   \\[.5ex] 
%  \hline
%\rule{0pt}{5pt}
  $\mat{M} = \mat{M}_C = \left(
                     \begin{array}{rrr}
                       1 & -1 & 0 \\
                       0 & 1 & -1 \\
                     \end{array}
                   \right)\trans
  $  & \code{C} & \code{A:C} & \code{B:C} & \code{A:B:C} \\
  \hline
\end{tabular}
\caption{Three-way design: Tests for between- (A, B) and within-S (C) effects
	are constructed using
	various \mat{L} and \mat{M} matrices. Table entries give the term actually tested
	via the general linear test in \eqref{eq:mglt-rep}.}
\end{table}

In the general case, various \mat{L} matrices provide contrasts or
select the particular coefficients
tested for between-S effects, while various \mat{M} matrices specify
linear combinations of responses for the within-S effects. This is
illustrated in \tabref{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
$\mat{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 $\mat{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, $\mat{L B M}=\mat{0}$
tests between-S \emph{differences} among the responses transformed by \mat{M}.

For more than one within-S factor, the full
\mat{M} matrices for various within-S terms are generated as Kronecker products
of the one-way \mat{M} contrasts
with the unit vector (\vec{1}) of appropriate size.
For example, with $c$ levels of factor C and $d$ of factor D,
\begin{equation}\label{eq:Mcd}
	\begin{split}
	\mat{M}_{C\otimes D}  & =  (\vec{1}_c, \mat{M}_C ) \otimes (\vec{1}_d, \mat{M}_D ) \\
		& =  ( \vec{1}_c \otimes \vec{1}_d, \vec{1}_c \otimes \mat{M}_D,
		  \mat{M}_C \otimes \vec{1}_d, \mat{M}_C \otimes \mat{M}_D ) \\
		& = ( \mat{M}_1, \mat{M}_D, \mat{M}_C, \mat{M}_{CD} ) \period
  \end{split}
\end{equation}
Each of the within-S terms combine with any between-S terms in an obvious way to
give an extended version of \tabref{tab:BWtab} with additional rows for
$\mat{M}_D$  and  $\mat{M}_{CD}$.

In passing, we note that \emph{all} software (\proglang{SAS}, \proglang{SPSS}, \proglang{R}, etc.)
that handles repeated measure designs through this extension of the MLM effectively
works via the general linear test \eqref{eq:mglt}, with either implicit
or explicit specifications for the \mat{L} and \mat{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 \mat{U} is unstructured.  Alternative analysis
based on mixed (or hierarchical) models, e.g.
\citep{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 (\secref{sec:profile}),
designs with multiple between- and within-S effects (\secref{sec:BW}),
and doubly-multivariate designs
(\secref{sec:doubly}), where two or more separate responses (e.g., weight loss
and self esteem) are \emph{each} observed in a factorial structure over
multiple within-S occasions. In \secref{sec:addendum} we describe
a simplified interface for these plots in the development versions of
the \textbf{heplots} and \textbf{car} packages. Finally (\secref{sec:compare})
we compare these methods with visualizations based on the mixed model.

\section{One sample profile analysis}\label{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 \code{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).

<<voc0, results=verbatim>>=
some(VocabGrowth,5)
@
A boxplot of these scores (\figref{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 \code{grade}).
We can see that vocabulary
scores increase with age, but the trend of means appears non-linear.

%<<voc1,fig=TRUE, include=FALSE>>=
%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
% <<voc1a, fig=TRUE, include=FALSE>>=
% xyplot(Vocabulary ~ Grade, groups=id, data=voc, type='b', pch=16:25)
% @
%

<<voc1,fig=TRUE, include=FALSE>>=
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}
\begin{figure}[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).}
\label{fig:voc-boxplot}
\end{figure}


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}$.

<<voc2, results=verbatim>>=
(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 \eqref{eq:mglt}, with $\mat{L} = \mat{I}$.
<<voc3, results=verbatim>>=
(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
\mat{H} and \mat{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 \mat{H} and \mat{E} matrices can be printed with \code{summary(Vocab.aov0)},
or extracted from the \code{Anova.mlm} object.  In this case, \mat{H}
is simply $n \overline{\vec{y}} \overline{\vec{y}}\trans$ and \mat{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}})$.
<<voc-Anova, results=verbatim>>=
Vocab.aov0$SSP     # H matrix
Vocab.aov0$SSPE    # E matrix
@

The HE plot for the \code{Vocab.mod} model
shows the test for the \code{(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)$.
<<markH0, keep.source=TRUE>>=
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 \code{grade8} and \code{grade9}
variables in \figref{fig:voc4}.  The \mat{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.
	
<<voc4, fig=TRUE, include=FALSE, width=6, height=6>>=
heplot(Vocab.mod, terms="(Intercept)", type="III")
mark.H0(0,0)
title(expression(paste("Multivariate test of ", H[0], " : ", bold(mu)==0)))
@

\begin{figure}[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 \mat{H} ellipse outside the \mat{E} ellipse
	signals that this $H_0$ is clearly rejected.}
\label{fig:voc4}
\end{figure}

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

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

Thus, to test the within-S effect of grade, we need to construct a \code{grade}
variable for the levels of grade and use this as a model formula,
 \verb|idesign=~grade| to specify the within-S design
 in the call to \code{Anova}.
<<voc5, results=verbatim>>=
idata <-data.frame(grade=ordered(8:11))
(Vocab.aov <- Anova(Vocab.mod, idata=idata, idesign=~grade, type="III"))
@

As shown in \tabref{tab:BWtab}, any such within-S test is effectively carried out using
a transformation \mat{Y} to \mat{Y M}, where the columns of
\mat{M} provide contrasts among the grades.  For the \emph{overall} test of \code{grade},
\emph{any} set of 3 linearly independent contrasts will give the same
test statistics, though, of course the interpretation of the parameters will differ.
Specifying \code{grade} as an ordered factor
(\code{grade=ordered(8:11)}) will cause \func{Anova} to use the polynomial contrasts
shown in $\mat{M}_{poly}$ below.
\[
\mat{M}_{poly} = \left(
\begin{array}{rrr}
-3 &  1  & -1 \\
-1 & -1  &  3 \\
 1 & -1  & -3 \\
 3 &  1  &  1 \\
\end{array}
 \right)
\quad\quad
\mat{M}_{first} = \left(
\begin{array}{rrr}
-1 & -1  & -1 \\
 1 &  0  &  0 \\
 0 &  1  &  0 \\
 0 &  0  &  1 \\
\end{array}
 \right)
\quad\quad
\mat{M}_{last} = \left(
\begin{array}{rrr}
 1 &  0  &  0 \\
 0 &  1  &  0 \\
 0 &  0  &  1 \\
-1 & -1  & -1 \\
\end{array}
 \right)
\]
Alternatively,
$\mat{M}_{first}$ would test the
gains in vocabulary
between grade 8 (baseline) and each of grades 9--11, while
$\mat{M}_{last}$ would test the
difference between each of grades 8--10 from grade 11. (In \R, these contrasts
are constructed with
\code{M.first <- contr.sum(factor(11:8))[4:1,3:1]}, and
\code{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,
$\mat{M} \vec{\mu}_y =\vec{0}$, such as we visualized in \figref{fig:voc4}.


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

%\[
%\mat{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 $\mat{M}_{poly}$, the MANOVA test for the \code{grade} effect is then testing $H_0 :
\mat{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.
<<voc6, results=verbatim>>=
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 \code{(Intercept)} term in \code{within.mod}
is identical to the test of \code{grade} in \code{Vocab.mod} at the beginning of
this subsection.

We can show this test visually as follows. \figref{fig:voc8}(a)
shows a scatterplot of the transformed linear and quadratic trend scores,
overlayed with a 68\% data ellipse.
 \figref{fig:voc8}(b) is the corresponding HE plot for these two variables.
Thus, we can see that the \mat{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).
 \figref{fig:voc8} is produced as follows:
<<voc8, echo=TRUE, fig=TRUE, include=FALSE, width=10, height=5>>=
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)
@

\begin{figure}[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 \figref{fig:voc4}, the size of the \mat{H} ellipse for Grade relative
to the \mat{E} ellipse shows the strength of evidence against $H_0$.}
\label{fig:voc8}
\end{figure}

We interpret \figref{fig:voc8}(b) as follows, bearing in mind that we are
looking at the data in the transformed space ($\mat{Y M}$) of the linear (slopes) and
quadratic (curvatures) of the original data (\mat{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 \mat{H} ellipse plots as a degenerate line because the \mat{H} matrix has rank 1
(1 df for the MANOVA test of the Intercept).
Its projection outside the \mat{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.

\section{Between- and within-S effects}\label{sec:BW}
When there are both within- and between-S effects, the multivariate and
univariate hypotheses tests can all be obtained together using \func{Anova}
with the \code{idata} and \code{idesign} specifying the within-S levels
and the within-S design, as shown above.
\func{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 \func{heplot} and related methods
it is again convenient
 to explicitly transform $\mat{Y} \mapsto \mat{Y M}$
for a given set of within-S contrasts, in the same way as done for
the \code{VocabGrowth} data. See \secref{sec:addendum} for simplified
code producing these HE plots directly, without the need for explicit
transformation.

To illustrate, we use the data from \citet{OBrienKaiser:85} contained in the
data frame \code{OBrienKaiser} in the \pkg{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 \code{session} in the analysis below.


<<obk1>>=
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).
<<obk2, results=verbatim>>=
table(OBK$gender, OBK$treatment)
@
The factors \code{gender} and \code{treatment} were specified with
the following contrasts, $\mat{L}_{gender}$, and $\mat{L}_{treatment}$, shown below.
The contrasts for \code{treatment} are nested (Helmert) contrasts
testing  a comparison of the
control group with the average of treatments A and B (\code{treatment1})
and the difference between treatments A and B (\code{treatment2}).
<<obk3, results=verbatim>>=
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.
<<obk-MANOVA, results=verbatim>>=
# MANOVA model
mod.OBK <- lm(cbind(pre, post, fup) ~  treatment*gender,  data=OBK)
mod.OBK
@

If we regarded the repeated measure effect of \code{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 (\code{post - pre}) and (\code{fup - post}).
<<obk4, results=verbatim, keep.source=TRUE>>=
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:
<<obk4, results=verbatim>>=
# 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 \code{print} methods
for \code{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 \code{summary}
method.%
\footnote{
In contrast, \proglang{SAS} \code{proc glm} and \proglang{SPSS} General Linear Model provide only
more complete, but often bewildering outputs that still recall the days
of \proglang{Fortran} coding in spite of more modern look and feel.
}
For example, using \code{summary(aov.OBK)} as shown below,	we can display all the multivariate
tests together with the \mat{H} and \mat{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.
<<obk5, results=hide>>=
# 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 \emph{nature} of these effects.
Ordinarily, from a data-analytic point of view, I would show traditional
plots of means and other measures (as in \figref{fig:voc-boxplot})
or their generalizations in effect plots
\citep{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 \emph{visual} summary of an MLM,
mirroring the tabular presentation from \code{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 (\mat{Y M}).

HE plots for between-S effects can be displayed for any pair of
responses with \func{heplot}.  \figref{fig:obk-HE1} shows this
for \code{pre} and \code{post}.  By default,
\mat{H} ellipses for all model terms (excluding the intercept) are displayed.
Additional MLM tests can also be displayed using the \code{hypotheses} argument;
here we specify the two contrasts for the \code{treatment} effect
shown above as \code{contrasts(OBK$treatment)}.
<<obk-HE1, fig=TRUE, include=FALSE>>=
# 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")
@

\begin{figure}[htb]
\centering
\includegraphics[width=.65\textwidth]{fig/plot-obk-HE1}
\caption{HE plot for the \code{mod.OBK} repeated measures model, showing between-S
effects and contrasts in the space of the \code{pre} and \code{post} variables.
Main effect means of the \code{treatment} (A, B, Control)
and \code{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 \code{pre = post}.}
\label{fig:obk-HE1}
\end{figure}

In \figref{fig:obk-HE1}, we see that the treatment effect is significant,
and the large vertical extent of this \mat{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 \code{gender} is not significant, but the HE plot shows that
that males are higher than females at both Pre and Post tests.
Likewise, the \code{treatment:gender} interaction fails significance,
but
the orientation of the \mat{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 \code{post} than \code{pre}.

Using \code{heplot3d(mod.OBK, ...)} gives an interactive 3D version of \figref{fig:obk-HE1}
for \code{pre}, \code{post}, and \code{fup},
that can be rotated and zoomed, or played as an animation.
<<obk-3d, eval=FALSE, echo=TRUE, keep.source=TRUE>>=
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 \func{pairs}
for the \code{mlm} object \code{mod.OBK}, with the result shown in \figref{fig:obk-HE2}.
<<obk-HE2, fig=TRUE, include=FALSE>>=
pairs(mod.OBK, col=c("red", "black", "blue", "brown"))
@
\begin{figure}[htb]
\centering
\includegraphics[width=.75\textwidth]{fig/plot-obk-HE2}
\caption{HE plot for the \code{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 \figref{fig:obk-HE1}.}
\label{fig:obk-HE2}
\end{figure}
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
$\mat{Y} \mapsto \mat{Y M}$, here using the (profile) contrasts for
\code{session}.
<<obk6, echo=TRUE, results=verbatim>>=
# 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  \tabref{tab:BWtab}, with these (or any other)
contrasts as $\mat{M}_{session}$,
the tests of the effects contained in \code{treatment*gender}
in \code{mod1.OBK} are identical to the interactions of these terms with \code{session}, as shown above
for the full model in \code{aov.OBK}.
The \code{(Intercept)} term in this model represents the
\code{session} effect.


The HE plot for within-S effects (\figref{fig:obk-HE3}) is constructed from the \code{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 \code{session}, as just described.
%\footnote{For technical reasons, it is necessary to calculate
%the \mat{H} matrices using Type III tests, and explicitly retain the
%\code{(Intercept)} term which is ignored by default in \func{heplot}.
%}
<<obk-HE3, fig=TRUE, include=FALSE>>=
# 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)
@

\begin{figure}[htb]
\centering
\includegraphics[width=.65\textwidth]{fig/plot-obk-HE3}
\caption{HE plot for the \code{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.}
\label{fig:obk-HE3}
\end{figure}

\figref{fig:obk-HE3} provides an interpretation of the within-S effects
shown in the MANOVA table shown
above for \code{Anova(mod.OBK)}.
We can see that the effects of
\code{session} and \code{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 \code{post-pre} difference
than to \code{fup-post}.

\subsection{Higher-order designs}
The scheme described above and the obvious generalization of \tabref{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 \func{lm}
call, e.g., \verb|~ treatment * gender|.

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


<<obk2-1, echo=TRUE>>=
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 \code{within} data frame looks like this:
<<obk2-2, results=verbatim>>=
str(within)
within
@

The repeated measures MANOVA analysis can then be carried out as follows:
<<obk2-3, results=verbatim>>=
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 \code{treatment}, \code{gender},
\code{session} and all interactions among them are identical to
what was found in the simplified analysis above.
Among the effects including \code{hour}, only the main effect is
significant here.

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

\begin{equation*}
\mat{M}_{session} = \left(
	\begin{array}{rr}
		-1 & 0 \\ 1 & -1 \\ 0 & 1 \\
   \end{array}
\right)	
\quad\quad
\mat{M}_{hour} = \left(
	\begin{array}{rrrr}
		-2 & 2  & -1 & 1 \\
		-1 & -1 &  2 & -4 \\
		 0 & -2 &  0 & 6 \\
		 1 & -1 & -2 & -4 \\
		 2 &  2 &  1 & 1 \\
   \end{array}
\right)	
\end{equation*}
Tests involving the interaction of \code{session:hour} use the Kronecker
product, $\mat{M}_{session} \otimes \mat{M}_{hour}$.

For HE plots, it is necessary to explicitly carry out the transformation
of $\mat{Y} \mapsto \mat{Y} \mat{M}_w$, where $ \mat{M}_w$
conforms to \mat{Y} and represents the contrasts for the within-S effect.
In the present example, this means that
$\mat{M}_{session}$ and $\mat{M}_{hour}$ are both expanded
as Kronecker products with the unit vector,
\begin{eqnarray*}
	\mat{M}_s & = & \mat{M}_{session} \otimes \vec{1}_5 \comma \\
	\mat{M}_h & = & \vec{1}_3 \otimes \mat{M}_{hour} \period  \\
\end{eqnarray*}

These calculations in \proglang{R} are shown below:
<<obk2-4, results=verbatim, keep.source=TRUE>>=
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 \code{M.h}, we can construct the within-model for all terms involving
the \code{hour} effect,

<<obk2-5>>=
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 \code{hour},
representing within-session slope and curvature.  \figref{fig:obk2-HE1}
is produced as shown below.  As shown by the \code{Anova(mod.OBK2, ...)} above,
all interactions with \code{hour} are small, and so these appear wholly contained
within the \mat{E} ellipse. In particular, there are no differences among groups
(treatment $\times$ gender) in the slopes or curvatures over \code{hour}.
For the main effect of \code{hour},
the linear effect is almost exactly zero, while the quadratic
effect is huge.
<<obk2-HE1, fig=TRUE, include=FALSE>>=
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()
@
\begin{figure}[htb]
\centering
\includegraphics[width=.6\textwidth]{fig/plot-obk2-HE1}
\caption{HE plot for the \code{mod.OBK2} repeated measures model, showing within-S
effects for linear and quadratic contrasts in \code{hour}.
As in \figref{fig:obk-HE3}, we are viewing hypothesis and error variation in the
transformed space of the repeated measures contrasts, here given by $\mat{M}_h$.}
\label{fig:obk2-HE1}
\end{figure}


The \func{pairs} function shows these effects (\figref{fig:obk2-HE2}) for all contrasts in \code{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 \code{hour}
are more complex than just the large quadratic component we saw in \figref{fig:obk2-HE1}.
<<obk2-HE2, fig=TRUE, include=FALSE>>=
pairs(mod.OBK2.hour, type="III", remove.intercept=FALSE, term.labels="hour", col=colors)
@
\begin{figure}[htb]
\centering
\includegraphics[width=.6\textwidth]{fig/plot-obk2-HE2}
\caption{HE plot for the \code{mod.OBK2} repeated measures model, showing within-S
effects for all pairs of contrasts in \code{hour}.}
\label{fig:obk2-HE2}
\end{figure}

\section{Doubly-multivariate designs}\label{sec:doubly}
In the designs discussed above the same measure is observed on all occasions.
Sometimes, there are two or more \emph{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 \mat{M} matrix in forming
the linear hypotheses to be tested via \eqref{eq:mglt-rep}.  For example, with two measures,
$Y_1, Y_2$ observed on three repeated sessions, the full \mat{M} matrix for the design is generated
as in \eqref{eq:Mcd} as

\begin{equation}\label{eq:Mmc}
	\mat{M}_{CM}
		= (\vec{1}, \mat{M}_{session} ) \otimes \mat{I}_2
    = \left(
	\begin{array}{rrr}
		1 & -1 & 0 \\ 1 & 1 & -1 \\ 1 & 0 & 1 \\
   \end{array}
  \right)	
  \otimes
  \left(
	\begin{array}{rr}
		1 & 0 \\ 0 & 1 \\
   \end{array}
\right) \period		
\end{equation}
	
In \R, we can express this as follows, using \code{M.measure} to represent $Y_1, Y_2$.

<<doubly, results=verbatim>>=
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 \func{Anova} does not (yet)%
\footnote{The new version of the \pkg{car} (2.0-0) on CRAN now includes
	enhanced \func{Anova} and \func{linearHypothesis} which perform these analyses.}
 allow such a doubly-multivariate design to be specified
directly.

\subsection{Example: Weight loss and self esteem}
To illustrate, we consider the data frame \code{WeightLoss} originally
from Andrew Ainsworth (\url{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 (\code{wl, se}) and month (\code{1:3}).

<<wl1, results=verbatim>>=
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 \figref{fig:wl-means},%
\footnote{
These plots were drawn using \func{plotmeans} in the \pkg{gplots} \citep{gplots}.
The code is not shown, but is available in the \proglang{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.
<<wl-means, echo=FALSE, fig=TRUE, include=FALSE, width=12, height=4>>=
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)
@
\begin{figure}[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.}
\label{fig:wl-means}
\end{figure}


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 \code{group},
which are labeled \code{group1} and \code{group1} respectively.

<<wl2, results=verbatim>>=
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 \code{group} effect.

<<wl3, results=verbatim>>=
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, \figref{fig:wl-HE1} shows the test of \code{group} and the two contrasts
for weight loss and for self esteem at months 1 and 2.
%<<wl-HE1, fig=TRUE, include=FALSE, keep.source=TRUE>>=
%heplot(wl.mod, hypotheses=c("group1", "group2"),
%	xlab="Weight Loss, month 1", ylab="Weight Loss, month 2")
%@

<<wl-HE1, fig=TRUE, include=FALSE, keep.source=TRUE, width=8, height=5>>=
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)
@


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

This is helpful, but doesn't illuminate the \emph{overall} \code{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 \mat{M} matrix for the measures, used in the
between-S analysis.  We use $\mat{M} = \mat{I}_2 \otimes \vec{1}/3$
so that the resulting scores are the means (not sums) for weight
loss and self esteem.

<<wl4, results=verbatim>>=
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 (\figref{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 \figref{fig:wl-means}.
<<wl-HE2, fig=TRUE, include=FALSE, keep.source=TRUE>>=
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")
@
\begin{figure}[htb]
\centering
\includegraphics[width=.6\textwidth]{fig/plot-wl-HE2}
\caption{HE plot for the \code{between.mod} doubly-multivariate design, showing
	overall between-S
effects for weight loss and self esteem.}
\label{fig:wl-HE2}
\end{figure}

Next, for the within-S analysis, we define the \mat{M} matrix for months,
using orthogonal polynomials representing linear and quadratic trends.
As before, the test of the \code{(Intercept)} term in these trend scores
corresponds to the month effect in the doubly-multivariate model,
and the \code{group} effect tests the group $\times$ month interaction.
<<wl5, results=verbatim>>=
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 (\figref{fig:wl-HE3}) can be produced as follows.
The \mat{H} and \mat{E} matrices are all $4 \times 4$, but the \mat{H}
matrices for the \code{month} and \code{group:month} effects are
rank 1 and 2 respectively.
<<wl-HE3, fig=TRUE, include=FALSE, keep.source=TRUE, width=10, height=5>>=
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()

\begin{figure}[htb]
\centering
\includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-wl-HE3}
\caption{HE plots for the \code{within.mod} doubly-multivariate design, showing
	the effects of \code{month} and the interaction \code{group:month}
for weight loss vs.\ self esteem. (a) Linear effects; (b) Quadratic effects.}
\label{fig:wl-HE3}
\end{figure}
\figref{fig:wl-HE3} shows the plots for
the linear and quadratic effects separately for weight loss vs.\ self esteem.
The plot of linear effects (\figref{fig:wl-HE3}(a)) shows that the effect of \code{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 \code{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
(\figref{fig:wl-HE3}(b)) is similar, except here,  differences in curvature
over months are driven largely by the difference between the \code{DietEx}
group from the others on weight loss.

The interested reader might wish to
compare the standard univariate plots of means in \figref{fig:wl-means} with the
HE plots in \figref{fig:wl-HE2} and \figref{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.

\section{Simplified interface: heplots 0.9 and car 2.0}\label{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 \textbf{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
\textbf{heplots} could only repeated measures by explicitly transforming
$\mat{Y} \mapsto \mat{Y M}$ and re-fitting submodels in terms of the transformed
responses.

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

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

For the \code{VocabularyGrowth} data, \figref{fig:voc8}(b) can be produced by

<<voc-new1,fig=FALSE, echo=TRUE, results=hide>>=
(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 \code{OBrienKaiser} data, the code for plots of between-S effects is the same
as shown above for \figref{fig:obk-HE1} and \figref{fig:obk-HE2}.
The HE plot for within-S effects involving \code{session} (\figref{fig:obk-HE3})
can be produced using \code{iterm="session"}:
<<obk-new0, fig=FALSE, echo=FALSE, results=hide>>=
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)
@
<<obk-new1,fig=FALSE, echo=TRUE, results=hide, keep.source=TRUE>>=
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 \code{hour} can be obtained using the
expanded model (\code{mod.OBK2}) for the 15 combinations of \code{hour} and \code{session}:
<<obk-new2,fig=FALSE, echo=TRUE, results=hide, keep.source=TRUE, eval=FALSE>>=
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}

\section{Comparison with other approaches}\label{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 \citep{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?

\subsection{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 \code{Sex} (\figref{fig:ortho-xyplot} (left)),

<<ortho-xyplot1-code, keep.source=TRUE, eval=TRUE>>=
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")
@
<<ortho-xyplot1,fig=TRUE, include=FALSE, echo=FALSE, width=6, height=6>>=
print(trellis.last.object())
@

%<<ortho-xyplot1,fig=TRUE, include=FALSE, echo=FALSE, width=6, height=6>>=
%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 \code{distance} on \code{age} (\figref{fig:ortho-xyplot} (right)).
<<ortho-xyplot2-code, keep.source=TRUE, eval=TRUE>>=
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')
                })
@
<<ortho-xyplot2,fig=TRUE, include=FALSE, echo=FALSE, width=6, height=6>>=
plot(trellis.last.object())
@


\begin{figure}[htb]
%\centering
 \begin{minipage}[c]{.49\textwidth}
  \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-xyplot1}
 \end{minipage}
 \begin{minipage}[c]{.49\textwidth}
  \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-xyplot2}
 \end{minipage}
\caption{Profile plot of \code{Orthodont} data, by sex (left); Pooled OLS and individual linear regressions on age, by sex (right)}
\label{fig:ortho-xyplot}
\end{figure}

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:

\begin{eqnarray}
 \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
\end{eqnarray}
where the vector of residuals for subject $i$ is $\vec{e}_{i \bullet} \sim N ( \vec{0}, \mat{R}_i )$.
(For this example, we take $\mat{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,
\begin{eqnarray}
  \beta_{0i} & = & \gamma_{00} + \gamma_{01} \mathrm{Sex}_{i} + u_{0i} \comma \\
  \beta_{1i} & = & \gamma_{10} + \gamma_{11} \mathrm{Sex}_{i} + u_{1i} \period
\end{eqnarray}
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 \code{lme} as follows:

<<Ortho.mix1, results=verbatim>>=
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, results=verbatim>>=
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,
<<Ortho.anova, results=verbatim>>=
anova(Ortho.mix1, Ortho.mix3)
@

To aid interpretation, we can plot the estimated fixed effects from the linear
model (m1) as follows, using the \code{predict} method for \code{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 \figref{fig:Ortho-fm}(left).
Similar code produces a plot of the cubic model \figref{fig:Ortho-fm}(right).

<<Ortho-fm1-code, keep.source=TRUE, eval=TRUE>>=
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")
@
<<Ortho-fm1,fig=TRUE, include=FALSE, echo=FALSE, width=6, height=6>>=
plot(trellis.last.object())
@

<<Ortho-fm3,fig=TRUE, include=FALSE, echo=FALSE, width=6, height=6>>=
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"))
@

\begin{figure}[htb]
%\centering
 \begin{minipage}[c]{.49\textwidth}
  \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-Ortho-fm1}
 \end{minipage}
 \begin{minipage}[c]{.49\textwidth}
  \includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-Ortho-fm3}
 \end{minipage}
\caption{Fitted values showing the estimated fixed effects of \code{age}, \code{Sex}, and their interaction
in the linear growth model (left) and cubic growth model (right).}
\label{fig:Ortho-fm}
\end{figure}


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.

\subsection{MVLM approach}

For the multivariate approach, the \code{Orthodont} data must first be reshaped to wide format with the
distance values as separate columns.
<<ortho-prep, results=verbatim>>=
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 \code{Sex} as the between-S factor.
Age is quantitative, so the intra-subject data frame (\code{idata}) is created
with \code{age} as an ordered factor.
<<ortho-mlm, results=verbatim>>=
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 \code{Sex} and \code{age} are highly significant and their interaction is nearly
significant.


\figref{fig:ortho-HE} shows HE plots for the between- and within-S effects, produced as shown
below.
The left panel plots the effect of \code{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 \code{age} main effect and the \code{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 \mat{E}
ellipse, we see that the main effect of \code{age} is essentially linear, and the overall
\code{Sex:age} effect is nearly significant due to a difference in slopes, but not curvature.

<<ortho-HE, fig=TRUE, include=FALSE, keep.source=TRUE, width=8, height=4>>=
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)
@
\begin{figure}[htb]
\centering
\includegraphics[width=\textwidth, trim=0 30 0 30]{fig/plot-ortho-HE}
\caption{HE plots for the \code{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 \code{age} and \code{Sex:age} (right).}
\label{fig:ortho-HE}
\end{figure}

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 \func{pairs} plot (not shown here),
<<ortho-pairs, fig=FALSE, include=FALSE, keep.source=TRUE, width=8, height=8>>=
pairs(Ortho.mod, idata=idata, idesign=~age, iterm="age", col=c("red", "blue", "brown"))
@
%\begin{figure}[htb]
%\centering
%\includegraphics[width=.8\textwidth, trim=0 30 0 30]{fig/plot-ortho-pairs}
%\caption{}
%\label{fig:ortho-pairs}
%\end{figure}
and, in particular, the panel corresponding to the nonlinear (quadratic and cubic)
components of trend, shown in \figref{fig:ortho-nonlin-HE}.
<<ortho-nonlin-HE, fig=TRUE, include=FALSE, keep.source=TRUE, width=8, height=8>>=
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")
@

\begin{figure}[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 \code{age} and \code{Sex:age},
showing that they are quite small in relation to within-S error.}
\label{fig:ortho-nonlin-HE}
\end{figure}

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 \code{Sex:age}
effect, as tested with \func{linearHypothesis}. The ``response transformation matrix''
shown below is equivalent to $\mat{M}_{poly}$ described earlier (\secref{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 \mat{H} and \mat{E} matrices.

<<ortho-lh1, results=verbatim>>=
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 $\mat{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),

<<ortho-lh2, results=verbatim>>=
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 \code{age} and the \code{Sex:age} interaction can be tested as follows,
using the contrasts for quadratic and cubic trends in age.

<<ortho-lh3, results=verbatim>>=
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 (\figref{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 \emph{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. 

\section{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 \citep{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 (\pkg{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.

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

%\bibliography{graphics,Rpackages,timeref}
\bibliography{repeated}

%\include{notes}

<<tidyup,echo=FALSE,results=hide>>=
options(old.opts )
@


\end{document}
