# How to get an overall p-value for an independent categorical variable using generalized estimating equations (geeglm) in R

I would like to know if there is a way to get an overall p-value for an independent categorical variable using generalized estimating equations in R. When I run the analysis, the output provides p-values for each level of the categorical variable, instead of an overall p-value for the variable itself.

One example of an output is this:

install.packages("geepack", dependencies=TRUE)
library(geepack
#some code using geeglm with binomial errors
# Results

Call:
geeglm(formula = modality ~ factor(smoking_status), family = "binomial",
data = airleak_cleaned_included, id = as.factor(mrn))

Coefficients:
Estimate Std.err  Wald Pr(>|W|)
(Intercept)                0.875   0.266 10.82    0.001 **
factor(smoking_status)1   -0.359   0.338  1.13    0.287
factor(smoking_status)2   -1.099   0.722  2.32    0.128
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = independence
Estimated Scale Parameters:

Estimate Std.err
(Intercept)        1  0.0597
Number of clusters:   176  Maximum cluster size: 1

• Edit to SO standards. I think this is really a stats question.
– DWin
Jan 28, 2022 at 18:07
• Yes, but its easier to do using score tests rather than Wald tests. In principle: fit an intercept-only model and save the results in, say, fit0. Then fit your model above above and save the results in, say, fit1. Then call anova(fit0, fit1). For more details, please ask on StackExchange. I am voting for migration. Jan 28, 2022 at 18:30

I'm posting an "answer", but I don't really want it considered as such because I don't have much statistical theory behind it. With Thomas Lumley's endorsement of this as acceptable practice, I'll consider it an answer. I can show the mechanics of getting a result using R. I seconded Limey's vote to migrate to stats.stackexchange.com after I searched that site for an answer to the general question of "how to do a principled model comparison" on GEE modeling results to get "overall significance" estimates from those GEEs. I didn't find what I considered to be a satisfying answer on either SO or CV.com.

This code installs and loads the R library 'geepack' and then builds 2 models (with one that sets the RHS of a formula to 1 to build a "Null model") based on the example code in ?geeglm and then runs anova on the two models. That package does have an anova.geeglm function which I take to imply that the authors think this is valid in at least some situations.

install.packages('geepack')
library(geepack)
data(dietox)
dietox$$Cu <- as.factor(dietox$$Cu)
mf <- formula(Weight ~ Cu * (Time + I(Time^2) + I(Time^3)))
gee1 <- geeglm(mf, data=dietox, id=Pig, family=poisson("identity"), corstr="ar1")
# warnings noted
gee1
coef(gee1)
vcov(gee1)
summary(gee1)
coef(summary(gee1))

mf2 <- formula(Weight ~ 1)
gee2 <- geeglm(mf2, data=dietox, id=Pig, family=poisson("identity"), corstr="ar1")
# more warnings

methods("anova") # to get available class-specific function names
[1] anova.coxph*       anova.coxphlist*   anova.geeglm*      anova.glm*         anova.glmlist*     anova.glmmPQL*
[7] anova.lm*          anova.lmlist*      anova.loess*       anova.loglm*       anova.mlm*         anova.negbin*
[13] anova.nls*         anova.polr*        anova.survreg*     anova.survreglist*
see '?methods' for accessing help and source code

anova(gee1,gee2)
#---------------
Analysis of 'Wald statistic' Table

Model 1 Weight ~ Cu * (Time + I(Time^2) + I(Time^3))
Model 2 Weight ~ 1
Df   X2 P(>|Chi|)
1 11 9327    <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


If this were "regular" regression, I would have confidence that this was a valid approach to model comparison and that a X^2 difference of 9327 between nested models differing by 11 degrees of freedom pointed to a high "statistical significance" and a corresponding low p-value. However, I'm being appropriately modest because this is a statistical area about which I have no experience.

• Yes, that's the right idea. This isn't an anova in the traditional sense, because it's not based on a likelihood, but it does the tests that the original question wanted. In principle a score test might be easier because you don't need to fit the more-complex model but in practice it actually isn't easier and fitting the model is not a problem. The warnings are because the outcome variable isn't integer, which a Poisson variable would be; you can use quasipoisson to silence them. Jul 26, 2023 at 7:32
• @ThomasLumley : Thanks, Thomas. I wasn't familiar with the expression Analysis of 'Wald statistic' Table but I guess its being done with inference based on the usual difference in chisquare-distributed ( X^2) statistic with the difference in adjusted df.
– DWin
Jul 27, 2023 at 0:42
• No, it's not being done with the usual difference in chisquared statistics, which doesn't work in this context. It's based on $(\hat\beta-\beta_0)^TV^{-1}(\hat\beta-\beta_0)$ for the subspace of $\beta$ that is in the larger model but not the smaller one. It doesn't have the nice property of traditional anova that the test statistics for A vs B and B vs C add to give the statistic for A vs C. Jul 27, 2023 at 1:21
• @ThomasLumley So those beta-hats are vectors and V^-1 is the full variance-covariance matrix? Shouldn’t the the displayed results say something along those lines?
– DWin
Jul 27, 2023 at 4:18
• Not the full covariance matrix, the covariance matrix for the subset being tested. That's the coded message of "Wald Statistic" rather than likelihood ratio Jul 27, 2023 at 7:22