# Do GEE and GLM estimate the same coefficients?

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?

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<-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

• And does this remain the case when we choose a non-linear link function? Jan 10, 2019 at 0:34
• The OP asked about a GLM, however, not a mixed GLM, correct? So it should be glm(fev ~ week) vs the geeglm. In that case, the GEE and GLM will agree, ut but differ in standard errors. This isn't necessarily the case for a multilevel model; see an example in a question I posted when the parameters of a multilevel model don't have anywhere near the same coefficients as a GEE: stats.stackexchange.com/questions/358231/… Jan 10, 2019 at 2:30
• @MarkWhite nice catch, I didn't even notice he was fitting a mixed effects model. It is curious that they have the same coefficients, since I was under the impression that GLMM and GEE do not produce equivalent estimates, as you said. Jan 10, 2019 at 2:34
• @Marcel, I've updated the question to address your comments as well as Mark's. I've also included Mark White's excellent post on multilevel models in the body of the answer for future users who might stumble upon the post. Jan 10, 2019 at 14:27
• @StatsStudent in your updated answer, you are using lm, which is not a generalized linear model, it is ordinary least squares. OLS regression's coefficients have a population averaged interpretation, which implies that they will have the same coefficients as a model estimated with generalized estimating equations - in fact it reduces to the same case under constant variance. My question was regarding the coefficients for a generalized linear model fit with ML and that of a generalized linear model fit with GEE - and if THOSE two models have the same estimated coefficients. E.g., a logit link Jan 10, 2019 at 19:17

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

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