Understanding the Limitations of GEE

Question

Review of GEE: Generalized Estimating Equations (GEE) are used to estimate the parameters of a generalized linear model. The main advantage of GEE is that they essentially allow you to model a multivariate correlated random variable, without knowing the exact form of this random variable (this is really useful because its really difficult to write a multivariate correlated non-random variable). GEE works through the Quasi-Likelihood function - by only knowing some relative moment based properties between the response and the covariates, we can still estimate model parameters, and these model parameters will still be consistent, unbiased and asymptotic normal (even we pick the wrong correlation structure!). The only drawback is however that GEE can only estimate population level parameters and not individual level parameters. I am trying to understand why this is.

The quasi-likelihood function for GEE is given by:

$$ Q(\beta; y) = \sum_{i=1}^{n} (y_i - \mu_i(\beta))^T V_i^{-1}(\beta, \alpha) (y_i - \mu_i(\beta))$$

where:

• $y_i$ is the response for the $i$-th subject,
• $\mu_i(\beta)$ is the mean response for the $i$-th subject,
• $V_i(\beta, \alpha)$ is the variance function for the $i$-th subject,
• $\beta$ are the fixed effects parameters, and
• $\alpha$ are the dispersion parameters.

The GEE are obtained by solving (where $D_i$ is the derivative of $\mu_i(\beta)$ with respect to $\beta$):

$$\frac{d}{d\beta} Q(\beta; y) = U(\beta) = \sum_{i=1}^{n} D_i^T V_i^{-1} (y_i - \mu_i(\beta)) = 0$$

Question: However, reading the above equations, it is unclear to me what happens if I would trick/force GEE to estimate individual level parameters instead of population level parameters. I tried to create a proof-of-concept problem to see what can go wrong if I use GEE to estimate individual level parameters.

Suppose we start with a Simple Longitudinal Model for the $i$-th person in the $j$-th group ( $y_{ij}$ is the response, $x_{ij}$ is the covariate, $b_{0ij}$ and $b_{1ij}$ are the individual-specific intercept and slope, and $\epsilon_{ij}$ is the error term):

$$y_{ij} = b_{0ij} + x_{ij}b_{1ij} + \epsilon_{ij}$$ $$b_{0i} = b_{0} + u_{0i}$$

If we write the GEE :

$$\mu_{ij}(\beta) = b_{0ij} + x_{ij}b_{1ij}$$

$$D_{ij} = \frac{\partial \mu_{ij}(\beta)}{\partial \beta} = (1, x_{ij})^T$$

$$U(\beta) = \sum_{i=1}^{n} D_i^T V_i^{-1} (y_i - \mu_i(\beta)) = 0$$ $$U(\beta) = \sum_{i=1}^{n} \sum_{j=1}^{m} (1, x_{ij})^T V_{ij}^{-1} (y_{ij} - b_{0ij} - x_{ij}b_{1ij}) = 0$$

Problem: Perhaps I can show that GEE provides unbiased Population estimates but biased Individual estimates?

• Population Level Bias: As bias is linked to the Expected Value, we can write:

$$U(\beta) = \sum_{i=1}^{n} D_i^T V_i^{-1} (y_i - \mu_i(\beta)) = 0$$

$$E[U(\beta)] = E\left[\sum_{i=1}^{n} D_i^T V_i^{-1} (y_i - \mu_i(\beta))\right]$$

Substituting $y_i = \mu_i(\beta) + \epsilon_i$, where $\epsilon_i$ is the error term, we get:

$$E[U(\beta)] = E\left[\sum_{i=1}^{n} D_i^T V_i^{-1} \epsilon_i\right]$$ If $\epsilon_i$ is independent of $D_i^T$, then $E[U(\beta)] = 0$ and the estimates are unbiased.

• Individual Level Bias:

$$E[U(\beta)] = E\left[\sum_{i=1}^{n} D_i^T V_i^{-1} (y_i - \mu_i(\beta))\right]$$

Substituting $y_{ij} = b_{0ij} + x_{ij}b_{1ij} + \epsilon_{ij}$ and $b_{0i} = b_{0} + u_{0i}$, we get:

$$E[U(\beta)] = E\left[\sum_{i=1}^{n} D_i^T V_i^{-1} \epsilon_{ij}\right]$$

If $\epsilon_{ij}$ is independent of $D_i^T$, then $E[U(\beta)] = 0$ and the estimates are STILL unbiased???

Thus it seems like GEE is just as suitable for estimating Individual parameters compared to Population parameters (contradiction). Where have I made a mistake?

• Note: Independence of $\epsilon_i$ and $D_i^T$ justification:

$$E[U(\beta)] = E\left[\sum_{i=1}^{n} D_i^T V_i^{-1} \epsilon_i\right]$$

$$E[U(\beta)] = \sum_{i=1}^{n} E[D_i^T V_i^{-1} \epsilon_i]$$

Since $\epsilon_i$ is independent of $D_i^T$, we can say ($D_i$ and $V_i$ are fixed quantities): $$E[D_i^T V_i^{-1} \epsilon_i] = D_i^T V_i^{-1} E[\epsilon_i]$$

Assuming that $E[\epsilon_i] = 0$, we get:

$$E[U(\beta)] = \sum_{i=1}^{n} D_i^T V_i^{-1} \cdot 0 = 0$$

Answer 1

There are two sorts of problems that will arise: general ones and boring technical ones for specific cases. I will start with the general one, where we can simplify to the independence working model and the identity link -- linear regression.

Suppose you have longitudinal data with $t=1,\dots,T_n$ observations on person $n$, and people $n=1,\dots,N$. You have an outcome $Y_{nt}$ on person $n$ at time $t$ and a single predictor $X_{nt}$. You might fit a model A $$E[Y_{nt}|X_{nt}]=\alpha+\beta X_{nt}$$ In this model, $\beta$ is an average difference in $Y$ per unit difference in $X$, averaged over times and people, and $\alpha$ is an intercept averaged over people.

You could also fit a model B $$E[Y_{nt}|X_{nt}]=\alpha_n+\beta X_{nt}$$ that has a separate intercept for each person (still just linear regression; no random effects or anything). Or you could even fit a model C $$E[Y_{nt}|X_{nt}]=\alpha_n+\beta_n X_{nt}$$ that has separate slope and intercept for each person.

As long as you have at least two distinct $X$ values for each person all these models can be fitted and will give unbiased estimates of the parameters they estimate (because it's just linear regression).

However, the last model will give absolutely terrible unbiased estimates in most settings, because each $\beta_n$ is estimated only from $T_n$ points. Unless $T_n$ is large, the variance of $\hat\beta_n$ will be large and poorly estimated, so you won't learn much by fitting the model. Your estimates of $\beta_n$ won't improve as $N$ increases, only as $T_n$ increases

If you want a good estimate of $\beta_n$, and $T_n$ isn't big enough to get one from model C, you need to accept some bias in your estimates. A reasonable approach is to take a weighted combination of $\hat\beta$ from model A and $\hat\beta_n$ from model C.

What weighted combination do we want? Well, one approach is to minimise the expected mean squared error, which depends on the bias of $\beta$ and on the variances of $\hat\beta$ and $\hat\beta_n$. We can estimate the variances (it's just linear regression) and we can estimate $E[(\beta-\beta_n)^2]$, so that's all feasible. We'll give more weight to $\hat\beta$ when $T_n$ is small (so $\hat\beta_n$ is bad) or when $E[(\beta-\beta_n)^2]$ is small, so that $\hat\beta$ is good.

After all this we have good (though not unbiased) estimates of individual-level slopes, done with just linear regression and no random effects. Except, it turns out that the estimates we get this way are exactly the linear mixed model estimates! So, it works, but it doesn't change anything.

The same is basically true for log links: Ulrike Grömping wrote a paper on using Poisson GEE to fit Poisson mixed models. However, the $\hat\beta_n$ are no longer exactly unbiased -- they can be infinite -- and the pooling approach doesn't work very well. Having unbiased estimating equations (as you still do) isn't enough on its own to guarantee unbiased estimators (you need information going to infinity to get consistent estimators, and even they aren't unbiased outside linear models)

For logistic regression the whole thing goes horribly wrong. Even in model B, with separate intercepts and a common slope, $\hat\beta$ is badly biased for $\beta$ (conditional logistic regression was invented to work around this problem).

Getting back to the original question: yes, you can fit individual parameters by ordinary regression or GEE. You probably won't have much information about each parameter, and this will be a problem. The type of problem it is varies depending on which model you're fitting. You're still probably going to need the regularisation that a random-effects model gives you if you want individual-level parameter estimates.