# Getting the variance-covariance matrix of regression coefficients in GEE

I fitted a GEE model using the function genZcor with user defined correlation matrix. I want to get the var-cov matrix of the regression coefficients. But the output provides only limited information.

I would be very much thankful if you could kindly let me know how to get it since I am struggling lot getting this.

The naive VC matrix of the parameter estimates is stored in the vbeta.naiv component of the GEE fit.

Here is an illustration with the Ohio data:

library(geepack)
data(ohio)
fm <- resp ~ age*smoke
gee.fit <- geese(fm, id=id, data=ohio, family=binomial,
corstr="exch", scale.fix=TRUE)
summary(gee.fit)


This is a basic model where the working correlation is considered symmetric with correlation $\rho$. Here, the within-cluster correlation is estimated at $\hat\rho=0.355$. Estimates with robust standard error (computed from a sandwich estimator, as detailed in the JSS paper, pp. 4-5) are shown below:

            estimate  san.se     wald       p
(Intercept) -1.90050 0.11909 254.6860 0.00000
age         -0.14124 0.05820   5.8889 0.01524
smoke        0.31383 0.18784   2.7912 0.09478
age:smoke    0.07083 0.08828   0.6438 0.42234


The robust and naive VC matrices are obtained as follows:

> gee.fit$vbeta [,1] [,2] [,3] [,4] [1,] 0.014182 0.002677 -0.014182 -0.002677 [2,] 0.002677 0.003387 -0.002677 -0.003387 [3,] -0.014182 -0.002677 0.035285 0.005348 [4,] -0.002677 -0.003387 0.005348 0.007793 > gee.fit$vbeta.naiv
[,1]      [,2]      [,3]      [,4]
[1,]  0.01407  0.002400 -0.014072 -0.002400
[2,]  0.00240  0.003139 -0.002400 -0.003139
[3,] -0.01407 -0.002400  0.034991  0.005373
[4,] -0.00240 -0.003139  0.005373  0.007938


We can check that the Wald statistics computed using those values (as the ratio of the estimates to their standard errors, which are the diagonal entries of the VC matrix) match the ones displayed in the summary table:

> (gee.fit$beta/sqrt(diag(gee.fit$vbeta)))^2
(Intercept)         age       smoke   age:smoke
254.6860      5.8889      2.7912      0.6438


(If you use geeglm instead, coefficients are available through the accessor coef(), and the robust VC matrix is stored in gee.fit$geese$vbeta where gee.fit now holds the results of the call to geeglm.)

A more detailed account of GEE computing is available in this excellent tutorial, Generalized Estimating Equations (GEE), by Søren Højsgaard and Ulrich Halekoh.

• I am trying to perform Tukey-type comparison using the glht package after fitting a GEE model. To do this I have to specify the vcov matrix in the glht function, which lead me to this post... Now, the problem is that the robust elements in the vcov matrix are smaller than the naive elements. Do you have inkling or former experience as to how this can be? – Notquitesure Aug 19 '14 at 13:49