Do GEE and GLM estimate the same coefficients?

Question

In a GLM, the likelihood equations depend on the assumed distribution only through the mean and the variance. The likelihood equations are

$$\sum_i^n (\frac{\partial \mu_i}{\partial \eta_i}) \frac{y_i - \mu_i}{Var(Y_i)}x_{ij} = 0, \quad (j = 1, ..., p)$$

and in the quasi-likelihood case, we just let $Var(Y_i) = v(\mu_i)$ be some function of the mean. For GEE, the response is extended to be multivariate, with an assumed correlation structure, with the quasi-likelihood equations.

Does this imply that GEE and GLM will have the same parameter (say $\beta$) estimates (population averaged) with the only difference being correct standard errors in the GEE case (Assuming clustered data?)

If the estimated coefficients are not the same, then what is the difference?

Answer 1

Yes. GEE and GLM will indeed have the same coefficients, but different standard errors. To check, run an example in R. I've taken this example from Chapter 25 of Applied Regression Analysis and Other Multivariable Methods, 5th by Kleinbaum, et. al (just because it's on my desk and references GEE and GLM):

library(geepack)
library(lme4)

#get book data from 
mydf<-read.table("http://www.hmwu.idv.tw/web/bigdata/rstudio-readData/tab/ch25q04.txt", header=TRUE)
mydf<-data.frame(subj=mydf$subj, week=as.factor(mydf$week), fev=mydf$fev)
#Make 5th level the reference level to match book results
mydf$week<-relevel(mydf$week, ref="5")

#Fit GLM Mixed Model
mixed.model<-summary(lme4::lmer(fev~week+(1|subj),data=mydf))
mixed.model$coefficients

                Estimate Std. Error     t value
(Intercept)  6.99850  0.2590243 27.01870247
week1        2.81525  0.2439374 11.54087244
week2       -0.15025  0.2439374 -0.61593680
week3        0.00325  0.2439374  0.01332309
week4       -0.04700  0.2439374 -0.19267241

#Fit a gee model with any correlation structure.  In this case AR1
gee.model<-summary(geeglm(fev~week, id=subj, waves=week, corstr="ar1", data=mydf))
gee.model$coefficients

            [Estimate   Std.err         Wald  Pr(>|W|)
(Intercept)  6.99850 0.2418413 8.374312e+02 0.0000000
week1        2.81525 0.2514376 1.253642e+02 0.0000000
week2       -0.15025 0.2051973 5.361492e-01 0.4640330
week3        0.00325 0.2075914 2.451027e-04 0.9875090
week4       -0.04700 0.2388983 3.870522e-02 0.8440338][1]

UPDATE

As Mark White pointed out in his comment, I did indeed previously fit a "single-level" Mixed Effects GLM. Since you didn't specify whether you wanted a "fixed effects" or "random" effects GLM model, I just picked "random" since that's the model fit in the book I selected from. But indeed, Mark is right that the coefficients do not necessarily agree in multilevel models, and someone provided a nice answer about that question previously. For your reference, I've added a "fixed" effects GLM model below using lm.

#Fit Traditional GLM Fixed Effect Model (i.e. not Random effects)
glm.fixed<-summary(lm(fev~week, data=mydf))
glm.fixed$coefficients
            Estimate Std. Error     t value     Pr(>|t|)
(Intercept)  6.99850  0.2590243 27.01870247 7.696137e-68
week1        2.81525  0.3663157  7.68531179 7.287752e-13
week2       -0.15025  0.3663157 -0.41016538 6.821349e-01
week3        0.00325  0.3663157  0.00887213 9.929302e-01
week4       -0.04700  0.3663157 -0.12830465 8.980401e-01

Note the first and second columns of the output in each model. They coefficients are identity, but standard errors differ.

You also added a comment which asked, "And does this remain the case when we choose a non-linear link function?" Note first that this is a different question since non-linear link functions generally aren't General Linear Models but Generalized Linear models. In this case, the coefficients do not necessarily match. Here's an example again in R:

#Fit Generalized Linear Mixed Effects Model with, say, Binomail Link
nlmixed.model<-summary(lme4::glmer(I(mydf$fev>mean(mydf$fev))~week+(1|subj), family="binomial", data=mydf))
nlmixed.model$coefficients

#Fit GEE model with, say, Binomial Link
nlgee.model<-summary(geeglm(I(mydf$fev>mean(mydf$fev))~week, id=subj, waves=week, family="binomial", data=mydf))
nlgee.model$coefficients

Answer 2

It depends on exactly what you mean and what you're assuming.

If you use the independence working correlation, the parameter estimates $\hat\beta$ in glm and GEE will be identical, with only the standard errors being potentially different

If you use another working correlation, the parameter estimates $\hat\beta$ will not be identical. They will estimate the same underlying parameter $\beta$ if $$E[Y_{it}|{\mathbf{X_{i}=x_i}]}=x_{it}\beta$$ for a linear model (or the equivalent with a link function for a generalised linear model). That is, they will estimate the same parameter if the marginal model is correctly specified with respect to past and future $x$, not just current $x$. (see this question).

If it's only true that $$E[Y_{it}|X_{it}=x_{it}]=x_{it}\beta$$ then the underlying parameters will differ (except with a diagonal working correlation matrix). What happens is that $Y_{it}-\mu_{it}$ is correlated with $X_{is}$ for $s\neq t$, and using a non-diagonal working correlation matrix brings terms like $(Y_{it}-\mu_{it})X_isV^{st}$ into the expectation of the estimating equations, where $V^{st}$ is the $(s,t)$ element of the inverse of the working covariance matrix.

It's also quite possible that even $$E[Y_{it}|X_{it}=x_{it}]=x_{it}\beta$$ isn't true, in which case the parameters and estimates will again depend on the working correlation matrix. For example, it could be that the difference in $Y$ for a one-unit difference in $X$ between people and within people might not be the same. The difference between smokers and non-smokers might be larger than that effect of starting or stopping smoking, either because duration matters or because there are other differences between smokers and non-smokers. Different working correlation structures will weight between-person and within-person contrasts differently. This paper talks about examples

Answer 3

I think it may not.

The estimated equations, in their formulas, depend on the inverse of the working covariance matrix. If we change it, the beta coefficients will change too, because the entire equation will change.

While in the GLM the working correlation is not applicable - it's fixed, independent.

And it can be shown:

For the independence, the beta are the same. Notice, that GEE reports the original ones taken from the GLM, so we can compare the output from the GLM and GEE instatntly:

> coef(summary(gee(Reaction ~ Days, id = Subject,data = sleepstudy,corstr = "independence")))
Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
running glm to get initial regression estimate
(Intercept)        Days 
      251.4        10.5 
            Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept)    251.4       6.61   38.03        6.63    37.91
Days            10.5       1.24    8.45        1.50     6.97

But not, when we change the structure. They are close, but not identical.

> coef(summary(gee(Reaction ~ Days, id = Subject,data = sleepstudy,corstr = "AR-M", Mv=1)))
Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
running glm to get initial regression estimate
(Intercept)        Days 
      251.4        10.5 
            Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept)    253.5      10.71   23.67        6.36    39.88
Days            10.5       1.67    6.25        1.44     7.27

or

> coef(summary(gee(Reaction ~ Days, id = Subject,data = sleepstudy,corstr = "AR-M", Mv=2)))
Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
running glm to get initial regression estimate
(Intercept)        Days 
      251.4        10.5 
            Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept)    253.8      10.83   23.45        6.45    39.34
Days            10.4       1.57    6.65        1.46     7.18
> 

or

> coef(summary(gee(Reaction ~ Days, id = Subject,data = sleepstudy,corstr = "unstructured")))
Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
running glm to get initial regression estimate
(Intercept)        Days 
      251.4        10.5 
            Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept)   252.04       7.33   34.40       10.09    24.97
Days            9.78       1.02    9.62        2.64     3.71