Longitudinal data and conditioning

Suppose you have a predictor \(X_{it}\) on person \(i\) at time \(t\), and an outcome \(Y_{it}\), where the outcome is generated by \[Y_{it}\mid X_{it}=x = a_i +1\times x+\epsilon\] where \(\epsilon\) are iid \(N(0,1)\) errors and \(a_i\) are iid \(N(0,1)\) random intercepts. I will show you the data generation code later, but this is honestly

Further suppose that you fit a linear model by least squares

lm(y~x, data=the_data)

## 
## Call:
## lm(formula = y ~ x, data = the_data)
## 
## Coefficients:
## (Intercept)            x  
##    -0.01108      0.99969

and also by maximum likelihood¹ with a linear mixed model

library(lme4)

## Loading required package: Matrix

lmer(y~x+(1|id),data=the_data, REML=FALSE)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: y ~ x + (1 | id)
##    Data: the_data
##       AIC       BIC    logLik  deviance  df.resid 
##  9844.550  9868.576 -4918.275  9836.550      2996 
## Random effects:
##  Groups   Name        Std.Dev.
##  id       (Intercept) 1.0586  
##  Residual             0.9668  
## Number of obs: 3000, groups:  id, 1000
## Fixed Effects:
## (Intercept)            x  
##   -0.004785     0.818309

The results are different. The least-squares coefficient clearly estimates 1 and the maximum-likelihood coefficient clearly estimates something that isn’t 1.

In this setting I’ve told you enough for you to verify that the least squares estimate should be consistent for 1, the value in the generating equation: the Normal equations will be unbiased at an intercept of zero and a slope of 1. What’s wrong with the maximum likelihood estimator? How did I rig it to fail?

The journal series Communications in Statistics does not have a terribly high reputation in the field of statistical methods, and it tends not to get the most interesting research. There are, of course, important exceptions. One such exception was by Margaret Pepe and Garnet Anderson in 1994, after they couldn’t get it published in various other venues. They were writing about “working correlations in GEE.² The same result was rediscovered by Pan, Louis, and Connett six years later. These papers were about GEE, and so naturally weren’t read by the sort of people who only do likelihood inference. But we nearly digress.

The profile score equations for the regression parameters \(\beta\) in a linear mixed model \[Y=X\beta+Zb+\epsilon\] with \(b\) and \(\epsilon\) zero-mean Normally distributed are \[\sum_i \sum_{s,t} x_{is}w_{ist}(y_{it}-\mu_{it})=0\] where \(\mu_{it}\) is the linear predictor \(x_{it}\beta\) and \(w_{ist}\) is the \((s,t)\) element of the modelled precision matrix \(\mathrm{var}[Y_{i\cdot}]^{-1}\) for person \(i\) at the current values of the variance parameters.

What does it take for the profile score equations to be unbiased? If \(w_{ist}= 0\) for \(s\neq t\) we only need \[E[Y_{it}|X_{it}=x]=x_{it}\beta.\] For general \(w\), though, we need \[E[Y_{it}|X_{it}=x,X_{is}=x']=x_{it}\beta.\] That is, we need the residual \(Y_{it}-x_{it}\beta\) to be uncorrelated with past and future \(x_{is}\) for person \(i\). We can, in principle, always arrange for past \(x\) to be uncorrelated with the residuals by including lagged \(x\) in the predictors. Conditioning on future \(x\) in a similar way can get metaphysically expensive.

So have we broken maximum likelihood? Not really. The likelihood we have written down is for the distribution of \(Y_{i\cdot}|X_{i\cdot}\). If we think it’s for the distribution of \(Y_{it}|X_{it}\) we only have ourselves to blame.³

Here is the data generation code

a<-rnorm(1000)
x0<-rnorm(1000)
y0<-rnorm(1000)+x0
x1<-rnorm(1000)+y0/sqrt(2)
y1<-rnorm(1000)+x1
x2<-rnorm(1000)+y1/2
y2<-rnorm(1000)+x2

the_data<-data.frame(x=c(x0,x1,x2),
                     y=c(y0,y1,y2)+rep(a,3),
                    id=rep(1:1000,3))

As I promised, every \(Y\) is generated only from the \(X\) at the corresponding time and iid errors. What I didn’t say was that \(Y\) affects future \(X\).

The relevance of GEE to all this is that GEE comes from using these profile score equations (or analogous equations in non-Normal models) to estimate regression relationships for \(Y_{i\cdot}|X_{i\cdot}\) or \(Y_{it}|X_{it}\). In GEE, the weight matrix \(w_{i\cdot\cdot}\) is a choice by the analyst, so a reasonable number of people are aware of the Pepe & Anderson and Pan, Louis & Connett results that inform the choice. The analyst might choose a diagonal weight matrix so as not to have to condition on future predictors. In a linear mixed model the choice of weight matrix happens under the hood, a consequence of the covariance model. The analyst is thus less likely to realise the implications. But the results on cross-time correlations apply to linear mixed models and Bayesian linear models, and mutatis mutandis to generalised linear mixed models. “No ad-hocery” will only get you so far; sometimes you actually want estimators set up for the purpose at hand.