\documentclass[12pt]{article}
\usepackage{chicago,graphicx}
\renewcommand{\baselinestretch}{1.4}
\begin{document}
\title{A Graphical Method of Exploring the Mean Structure
 in  Longitudinal Data Analysis}

\author{Julian J. Faraway \\
%  Technical Report \#281 of the Department of Statistics, \\
%  University of Michigan \\
  {\small email: \texttt{faraway@umich.edu}}}

%\author{Anonymous}

\date{\today}

\maketitle

\begin{abstract}
  In a longitudinal study, individuals are observed over some
  period of time. The investigator wishes to model the responses
  over this time as a function of various covariates measured on
  these individuals. The times of measurement may be sparse and
  not coincident across individuals. When the covariate values
  are not extensively replicated, it is very difficult to 
  propose a parametric model linking the response to the covariates
  as plots of the raw data are of little help.
  
  Although the response curve may only be observed at a few points, we
  consider the underlying curve $y(t)$. We fit a regression model
  $y(t)=x^T\beta(t)+\epsilon(t)$ and use the coefficient functions
  $\beta(t)$ to suggest a suitable parametric form. Estimates of
  $y(t)$ are constructed by simple interpolation and appropriate
  weighting is used in the regression. We demonstrate the method
  on simulated data to show its ability to recover the true structure
  and illustrate its application to some longitudinal data from the
  Panel Study of Income Dynamics.

\end{abstract}

\textbf{Keywords:} Functional Data Analysis, Nonparametric Regression,
Repeated Measures, Curve Estimation, Exploratory Data Analysis.


\section{Introduction}
\label{sec:intro}

Longitudinal data arises when a subject is observed over a period of
time, ${\cal T}$, or over some other continuous variable. Suppose the
response of the $i^{th}$ individual is $y_i(t)$ where $t \in {\cal
  T}$. Now in practice one never observes the whole continuous
function $y_i(t)$ but only at a finite number of points $t_{ij} \in
{\cal T}$ where $i=1, \dots ,n$ and $j=1, \dots , m_i$.  In many
applications, $m_i$ is quite small and sometimes the points of
observation, $t_{ij}$, differ from individual to individual. Even when
an experiment has been designed to make the points of observation the
same, data is often missing and the pattern will be broken. Now even
though the $m_i$ may be relatively small and there might seem to be
little hope of reconstructing $y_i(t)$, we take the point of view that
it is worthwhile considering the underlying function rather than the
vector $y_i(t_{ij})$ , $ j=1, \dots , m_i$.

Suppose that we observe covariates $x_i$, a vector of length $p$,
associated with each individual $i$. We wish to develop a model of
the form
\begin{displaymath}
  y_i(t) = f(t,x_i) + \epsilon_i(t).
\end{displaymath}

Our main objective is to determine an appropriate form for the mean
structure $f()$. With a sufficiently large amount of data (in the
sense that the range ${\cal T}$ is more densely covered) we might
attempt to formulate an at least partially nonparametric model. When
there is less data or a compact formulation of the model is required,
then a fully parametric form for $f()$ will be desirable. The primary
difficulty is in choosing a suitable parametric form. This is the
problem this article sets out to solve.

The error structure is contained in $\epsilon_i(t)$. I assume that
each $\epsilon_i(t)$ is an independent realization of a stochastic
process with mean 0 and covariance function $\gamma(s,t),\quad s,t \in
{\cal T}$.  So I allow for correlated errors within individuals, but
not between individuals. The errors do not depend on the covariates.
However, the purpose of this article is not to come up with novel
methods for determining the form of $\gamma()$ as this problem has
been investigated by other authors --- see for example \citeN{izenman:89}
and \citeN{grambsch:95} among others.

Little work exists on exploratory methods for finding an appropriate
mean structure for longitudinal models. When there are no covariates
or sufficiently large numbers of individuals have the same values of
the covariates, then it is possible to plot
the data to obtain a suggestion for the appropriate form for $f()$.
This is described in \citeN{rice:91} and \citeN{diggle:95}.  When the
covariates are continuous in nature, these methods will not be
effective.  When the points of measurement, $t_{ij}$ are in common
across individuals, then it may be possible to avoid specifying a
functional form for the mean structure by saturating the model with
enough parameters to model the response at each time point, but this also
requires some replication of covariate values. The alternative is to
use contextual information as much as possible, but in practice,
often the parametric form must be guessed and experimented with which is
less than satisfactory. Other semiparametric approaches to related
problems can be found in \citeN{hoover:96}, \citeN{brumbach:96}
and \citeN{zhang:98}. A fully functional approach is described
in \citeN{ramsay:97}.

The method we propose here is exploratory and graphical in nature.
Parametric assumptions are avoided because the method is intended to
suggest a suitable parameterization.  The spirit of the method is
similar to the ACE method of \citeN{breiman:85} and the generalized
additive models of \citeN{hastie:90} which although nonparametric in
nature can be used to suggest suitable parametric forms for functions
of covariates in regression problems.  The method can also be used as
a diagnostic as well as an exploratory method when the appropriate
mean structure is thought to be known, but when some caution is
required.

Methodology is laid out in Section~\ref{sec:method}.  We check the
efficacy of the method for some simulated data in Section~\ref{sec:simul}
where we know the true model and demonstrate the method on some real
data in Section~\ref{sec:rde}.  A discussion follows in
Section~\ref{sec:disc}.


\section{Method}
\label{sec:method}


\subsection{Functional Regression Analysis}
\label{sec:fra}

Suppose that we may observe the whole functional
responses, $y_i(t)$, which are assumed to arise from the model:
\begin{equation}
  \label{eq:fm}
  y = X \beta + \epsilon
\end{equation}
where $\beta$ is a vector of functions $(\beta_1(t), \dots
\beta_p(t))^T$ and $X$ is the familiar $n \times p$ design matrix
formed from the $p$-vector valued covariates $x_i$, $i=1, \dots n$. 
As in scalar (i.e.\ scalar $y$)
regression, the first column of $X$ will usually all be ones and
categorical predictors can be handled by assigning appropriate dummy
variables.  Various transformations of the predictors can be
incorporated in $X$ as in scalar regression.  Also $y$ is a vector
of response functions $(y_1(t), \dots y_n(t))^T$ and $\epsilon$ is a
vector of error functions $(\epsilon_1(t), \dots , \epsilon_n(t))^T$.

Suppose we choose $\hat\beta$ to minimize $\sum_{i=1}^n || y_i - x^T_i
\beta ||^2$ where $||\cdot ||$ is the $L_2$ norm on ${\cal T}$. In
parametric longitudinal analysis, distributional assumptions are made
and the model can be fit using maximum likelihood methods. Of course,
this fit will be more efficient than using the least squares criterion,
but at the price of making some possibly unjustifiable assumptions.
By considering each $t \in {\cal T}$ separately and provided that $X$
has full rank, as in the usual regression situation, it is clear that
the solution is
\begin{displaymath}
  \hat\beta = (X^T X)^{-1} X^T y .
\end{displaymath}

The optimality of this estimator is discussed in \citeN{ramsay:91}.
Let $\hat y = X\hat\beta$ and $\hat\epsilon = y - \hat y$. We could
also use a weight matrix $W$ to get $\hat\beta = (X^T WX)^{-1} X^T Wy $.
Since the solution is computed pointwise, the weight matrix can change
as t varies. We will need this feature below. See \citeN{faraway:97} or
\citeN{ramsay:97}
for a development of this model when the $y(t)$ are observed more densely.

\subsection{Range of Representation}
\label{sec:range}

Not all longitudinal models of interest may be represented in the form
(\ref{eq:fm}).  Clearly covariates that are explicitly a function of
time cannot be represented within this model without some further
generalization. For simplicity consider a situation where there is
only one (univariate) covariate $z$. The model with an intercept will
take the form
\begin{displaymath}
  y(t) = \beta_0 (t) + z \beta_1 (t) + \epsilon(t).
\end{displaymath}

The functional form of $\beta_0(t)$ and $\beta_1(t)$ are estimated
which may then suggest a suitable parametric model. A wide range of
parametric forms fall within this setup. Further generalization is
possible by allowing for transformations of $z$ as in

\begin{equation}
  \label{eq:addform}
  y(t) = \beta_0 (t) + \sum_l g_l(z) \beta_l(t) + \epsilon(t)
\end{equation}
where we the $g_l()$'s must be specified.  We provide no way of
estimating the correct forms for the $g_l()$ and in practice some
experimentation and perhaps luck would be needed to find these forms.
It is conceivable that with a sufficiently large number of individuals
we might attempt some nonparametric estimation of $g_l(z)$, but we do
not investigate that in this article.

Models such as
\begin{displaymath}
  y(t) = \log(z t + z^2 t^2) + \epsilon(t)
\end{displaymath}
do not fall within the factorized form (\ref{eq:addform}) and lie
outside the range of representation. We would have no chance of
discovering such a model using our method, although prior substantive
knowledge would likely be necessary to choose such a model anyway. In
some cases, an appropriate transformation of the response might allow
the model to be expressed in a representable form.

\subsection{Interpolation and Weighting}
\label{sec:weight}

In practice, we do not observe $y_i(t),  i=1, \dots ,n$, only
$y_{ij},  j=1, \dots m_i$ at times $t_{ij}$ where 
\begin{displaymath}
  y_{ij} = x_i^T \beta(t_{ij}) + \epsilon_{ij}
\end{displaymath}
where $\epsilon_{ij} = \epsilon_i(t_{ij})$.

To estimate $\beta$, we will first estimate $y_i(t)$ on a fine grid of
points $t_j, j=1, \dots m$ where the grid of points is in common for
all $i$. Call these estimates $y^l_i(t)$.

We will use linear interpolation to construct the $y^l_i(t)$ using a
constant for the estimate before the first and after the last
observation in the range ${\cal T}$, i.e. $y^l_i(t) = y_{i1}$ for $t <
t_{i1}$ and $y^l_i(t) = y_{im_i}$ for $ t > t_{im_i}$. This is not the
standard nonparametric regression problem. We wish to estimate $y(t)$,
not $Ey(t)$. Certainly, linear interpolation is a poor method of
estimating $Ey(t)$, but that is not our objective. We wish to estimate
$\beta(t)$. One might contemplate some smoothing, but there are
several reasons for preferring linear interpolation:

\begin{enumerate}
\item The $m_i$ might be very small in which case more sophisticated
  methods of nonparametric regression will not be appropriate --- our
  method of linear interpolation will work with only one point.
  
\item Not smoothing the estimates at all means that the bias will be
  minimized. At the points of observation $t_{ij}$ there is no bias
  at all.  When $\hat\beta (t)$ is constructed, the bias will be further
  reduced by the averaging. We can always smooth the estimates of
  $\beta(t)$, but if we smooth the $y_{ij}$ too much, we may lose
  features that cannot later be recovered. One exception to this idea
  of avoiding smoothing the $y_{ij}$ is that we should be beware of
  outliers. We might seek to eliminate these by a pre-screening.
  
\item In the usual nonparametric regression problem, the errors at each
  time point are assumed to be independent. In that case, linear
  interpolation is a very poor method asymptotically. However, for
  that problem, the objective is to estimate the mean curve $Ey(t)$.
  Here, we are more interested in the observed curve $y(t)$
  and the appropriate asymptotics will be different.  Smoothing might
  filter out some measurement error, but we would rather retain the
  random effects and serial correlation parts of $y(t)$ so that
  the error structure of $\epsilon (t)$ can be analyzed.

\item Linear interpolation is simple and easily implemented without
  special software.
\end{enumerate}

At a particular point $t$, we have
\begin{displaymath}
  y^l_i = x_i^T \beta + b_i + \epsilon_i , \quad i=1, \dots n,
\end{displaymath}
where $\mathrm{Var}(\epsilon) = \Omega$ (which is a diagonal matrix by
the assumption of independence between individuals) and biases $ b_i$.
The presence of the biases $b_i$ is what distinguishes this from the
usual regression set-up. At the points $t_{ij}$, the bias is zero.
How should the weight matrix $W$ be chosen in this case?

The estimator is 
\begin{displaymath}
  \hat\beta = (X^T W X)^{-1} W y^l = Ay^l
\end{displaymath}
so that $\mathrm{Bias}(\hat\beta) = Ab $ and
  $\mathrm{Var}(\hat\beta) = A\Omega A^T$
where $A = (X^T W X)^{-1} W$. Notice that the bias of $\hat\beta$
is a linear combination of the biases $b_i$. Since some cancellation
is likely to occur, the bias of $\hat\beta$ is likely to be reduced. Now
$ \mathrm{MSE}(\hat\beta) = A[bb^T + \Omega] A^T$
which suggests choosing $W = [bb^T + \Omega]^{-1}$ which gives
$  \mathrm{MSE}(\hat\beta) = [X^T[bb^T+\Omega]^{-1}X]^{-1}$.

It would now be nice to have estimates of the bias and variance.
However, given the lack of information about the true mean response
and the error structure combined with possibly very few observations
per individual it is very difficult to find good estimates of the bias
and variance. Various more ambitious methods, than the one I will
recommend below, failed because of a lack of stability in these
estimates.

Note that $\mathrm{MSE}(\hat\beta)$ has non-zero off-diagonal elements
so a separate estimate of the bias is required to approximate the
optimal weight matrix. However, since it is difficult to estimate
the bias we propose to ignore the off-diagonal elements --- we will
use weighted least squares rather than generalized least squares.

If we assume that $y(t)$ is locally quadratic, this suggests that
the bias at $t$ might be estimated by a term of the form
\begin{displaymath}
  b_i(t) = y_i(t)-y^l_i(t) = \pm \gamma \{ t^2+t_ut_l-t(t_u+t_l) \}
\end{displaymath}
for $t \in [t_l, t_u]$ where the closest bracketing timepoints are
at $t_l$ and $t_u$. This expression may be derived
by computing the distance between a line and a parabola and
by observing that $b_i(t_l)=b_i(t_u)=0$. $\gamma > 0$ is
a tuning parameter which should be chosen with respect to the
speculated amount of curvature in $y(t)$. We cannot determine
the sign of $b_i(t)$ because we do not know whether $y_i(t)$
has positive or negative curvature in the local interval.
This also illustrates why we chose to ignore the off-diagonal
elements of $W$, whereas the diagonal elements of bias are squared
and so the sign is irrelevant.
Bias is not necessarily greatest in the middle of a large gap
between observations, but without any further information it
is safest to downweight any $y^l_i(t)$ that is estimated using
relatively distant observations.
For $t < t_{i1}$, we suggest $\gamma (t_{i1}-t)^2$
and for $t > t_{im_i}$ , $\gamma (t-t_{im_i})^2$. 

Now suppose we assume that the error process, $\epsilon(t)$
has variance, $\sigma^2$, which is constant in time, then we
will choose the weights as
\begin{displaymath}
  w_i(t) = {1 \over \sigma^2 + b_i^2(t)} \equiv {1 \over 1 + b_i^2(t)}
\end{displaymath}
The last equivalency arises by observing that the multiplication
of weights by a constant term has no effect on the calculation of
the regression coefficients and that we may absorb constant multipliers
of the tuning parameter $\gamma$. If we know that the variance is
not constant in time and takes the form $\sigma^2(t) \equiv 
\sigma^2 v(t)$, for known $v(t)$, then an appropriate form for
the weights will be
\begin{displaymath}
  w_i(t) =  {1 \over v(t) + b_i^2(t)}
\end{displaymath}

The $w_i(t)$ may be plotted to guide the selection of $\gamma$.  With
larger datasets, $\gamma$ can be made bigger to localize the effect of
observations.  Since the aim of the analysis is exploratory, manual
adjustment of $\gamma$ is adequate --- there is no necessity for
smoothing parameter selection methods.  Other weighting systems can
easily be implemented if more information is available.

\subsection{Constructing the parametric model}
\label{sec:para}

Pointwise standard errors may easily be constructed using the standard
regression formula: $ se (\hat\beta_k(t)) = \sqrt{(X^T W X)^{-1}_{kk}}
\hat\sigma (t)$. We may then plot $\hat\beta_k(t) \pm 2
\mathrm{se} (\hat\beta_k(t))$ against $t$ for each $k=1, \dots , p$.
We can now
\begin{enumerate}
\item Gain insight on a suitable parametric form for $\beta (t)$, using
  the pointwise confidence bands as guide to what the reasonable forms
  are.

\item Check these estimated $\hat\beta(t)$ against an already proposed
  parametric model looking for possible inadequacies.
\end{enumerate}

It would be better to use simultaneous rather than pointwise confidence
bands, but this would require more knowledge of the error structure
than we care to presume.

\section{Simulation}
\label{sec:simul}

To see how well our method could recover the true $\beta(t)$ we 
 generated data from a known model:
\begin{equation}
  \label{simod}
  Y_{ij} =  \beta_0(t_{ij}) + z_{1i} \beta_1(t_{ij}) + z_{2i} \beta_2(t_{ij}) 
  + \epsilon_{ij}
\end{equation}
where
\begin{eqnarray*}
  \beta_0(t) &=& t/2 \\
  \beta_1(t) &=& \left\{
    \begin{array}{ll}
      0.5 - t & \mathrm{if} \quad t<0.5 \\
      0 & \mathrm{otherwise}
    \end{array}
    \right. \\
    \beta_2(t) &=& \left\{
    \begin{array}{ll}
      t - 0.5 & \mathrm{if} \quad  t>0.5 \\
      0 & \mathrm{otherwise}
    \end{array}
    \right. 
\end{eqnarray*}
and $z_{1i}$ and $z_{2i}$ are i.i.d.\ $U[-1,1]$ and $\epsilon_{ij}$ is
i.i.d.\ $N(0,0.1^2)$, $i=1, \dots 20$ and $j=1, \dots 5$. The plot
of the raw data is shown in the first panel of Figure~\ref{fig:sim}.

\begin{figure}[htbp]
  \begin{center}
    \leavevmode
    \includegraphics*{sim.ps}
    \caption{The first panel shows the data and the remaining three panels show the estimated $\beta(t)$ (solid), two standard error bands (dotted) and
      the true $\beta(t)$ (dashed)}
    \label{fig:sim}
  \end{center}
\end{figure}

Notice that it would be almost impossible to guess the correct form of
the model by looking at the plot of the raw data. There is no
suggestion of a change in trend at $t=0.5$.  However, our method
recovers the true functional form very well. We see there is linear
increase in the response over time unrelated to either covariate,
$\hat\beta_0(t)$, that the effect of $z_1$ decreases linearly over
time to zero at $t=0.5$ and is constant after that, $\hat\beta_1(t)$,
and the complimentary effect for $z_2$, $\hat\beta_2(t)$. We should
not be too concerned that the true $\beta(t)$ ventures slightly
outside the bands as these are pointwise, not simultaneous,
confidence bands.

Obviously one example is not completely convincing, but we have tried
many more.  I invite the reader to explore their own simulated models
using S+ software available from
\texttt{www.stat.lsa.umich.edu/{\char'176}faraway}.
%\texttt{www.anonymous.edu}.

\section{Example}
\label{sec:rde}

The Panel Study of Income Dynamics (PSID), begun in 1968, is a
longitudinal study of a representative sample of U.S. individuals
(\citeN{hill:92}). The study is conducted at the Survey Research
Center, Institute for Social Research, University of Michigan and is
still continuing. There are currently 8,700 households in the study
and many variables are measured. I chose to analyze an random subset
of this data, consisting of 85 heads of household who were aged
between 25-39 in 1968 and had complete data for at least 11 of the
years between 1968 and 1990. The variables included were annual
income, gender, years of education and age in 1968. The analysis here
is intended to demonstrate the methodology.  Since we have selected a
subsample and have excluded various important variables, this should
not be construed as a complete attempt to draw conclusions about the
subject matter. We will attempt to model the annual income over the
period of interest as a function of the covariates mentioned above.
We will consider a model for individual $i$ of the form:

\begin{equation}
  \label{eq:psidm}
  \log \textrm{(income)}_i(t) = \beta_0(t) + \beta_g (t) \textrm{gender}_i +
  \beta_e \textrm{education}_i + \beta_a \textrm{age}_i \epsilon_i (t)
\end{equation}

Since in this example, the response, income, is measured at yearly
intervals, it is convenient (although not necessary) to make the grid
of timepoints at which the functions will be estimated coincident with
the years. For this data it also convenient to make $\gamma$ very
large so that missing observations, although interpolated are
effectively not used.  A plot of the estimated coefficient functions
is shown in Figure~\ref{fig:psid} with a lowess-smoothed fit to aid
interpretation.

\begin{figure}[htbp]
  \begin{center}
    \leavevmode
    \includegraphics*[width=6in]{psid.ps}
    \caption{The coefficient functions $\hat\beta(t)$ (solid line) for
      the model (\ref{eq:psidm}) with $\pm 2$ standard errors (dotted
      lines) and Lowess smooth of $\hat\beta(t)$ (dashed line). 
      For gender, the effect of being female is observed.}
    \label{fig:psid}
  \end{center}
\end{figure}

The constant function $\hat\beta_0 (t)$ is not the mean income at each
timepoint so we should not be disturbed that this does not increase
monotonely with time. Age may well not be important at all given that
the zero line lies within the two standard error bands, but there is
some indication of a quadratic effect. The effect of Education seems
approximately linear with the relative effect of Education increasing
over time. The coding of the gender variable meant that the effect of
being female is represented.  We can see that women earn substantially
less although the difference is decreasing approximately linearly with
time.

This analysis gives an idea of what functional forms for a parametric
longitudinal model might be appropriate. Some more work
needs to be done in specifying a model, but at least we now have
good indication of how to do this.

\section{Discussion}
\label{sec:disc}

Our method is not appropriate for all types of longitudinal data. For
example, where each individual is measured at the same small number of
time points, then there would be little value in using our method. However,
although data may seem to have been collected at the same time points,
the concordance is lost when a different origin is used --- for
example, patients may be measured at equal time points from the start
of treatment, but if chronological age is used as the origin, the times
of measurement are no longer the same for each individual. We also have
not described how to handle time-varying covariates, although a simple
interpolation approach for this may also be feasible. 
Nevertheless, the method provides a simple way to explore certain types
of longitudinal data.


\bibliographystyle{chicago}
\bibliography{fda}


\end{document}
