How to perform model selection in GEE in R

Question

I would like to do model selection for generalized estimating equations (GEE). Pan (2001) is most frequently cited for developing a method using QIC. I am wondering if anyone knows of a way to do this in R? I am currently using the package 'geepack' for my GEE analysis. I heard about 'yags' but it is not available for R version 2.13 or 2.14 and because of other things I'm doing, I'd rather not try to revert to an older version of R.

Any Help would be appreciated. Thanks, Dan

Automatic model selection is generally a bad idea. I wrote a lengthy answer about this to the following question: http://stats.stackexchange.com/questions/20836/algorithms-for-automatic-model-selection/20856#20856 which you may want to read. Using the AIC (or its variants, e.g. QIC) *is* better than traditional approaches, but is *still* subject to all of the same problems. — gung - Reinstate Monica, Jan 26 '12 at 20:11
Thanks for the thoughts gung. I appreciate your lengthy discussion and examples that you link to. Let's say I did what to compare models using QIC (say to test models specifying specific hypotheses), any idea how I'd go about it? — djhocking, Jan 27 '12 at 21:29

score 3 · Answer 1 · answered Jan 28 '12 at 05:48

If you want to select amongst pre-specified models, this should work the same with GEE as elsewhere. For example, if you were comparing a nested model to a full model, you could test that. If the models weren't nested, you could use an informational criterion (such as the QIC) to help adjudicate between them. Another approach is to use the Parametric Bootstrap Cross-fitting Method; this is a very solid approach, but computationally very expensive.

djhocking · Accepted Answer · 2014-01-14T03:08:33.100

[UPDATE: I improved on the code below and made a small R package hosted on GitHub: https://github.com/djhocking/qicpack]

I figured out a solution for calculating QIC from geepack package output. My code is below. This is one of the first functions I've ever written, so I apologize if it's messy but hopefully others find it useful. I definitely recommend reading gung's thoughts on model selection (linked above) before using this or any other information criterion model selection techniques (e.g. AIC, BIC, DIC). Also much of this code was pieced together from other sources, which I tried to reference at the start. I also received valuable input from Jun Yan, the geepack author.

######################################################################################
# QIC for GEE models
# Daniel J. Hocking
# 07 February 2012
# Refs:
  # Pan (2001)
  # Liang and Zeger (1986)
  # Zeger and Liang (1986)
  # Hardin and Hilbe (2003)
  # Dornmann et al 2007
  # # http://www.unc.edu/courses/2010spring/ecol/562/001/docs/lectures/lecture14.htm
######################################################################################

# Poisson QIC for geeglm{geepack} output
# Ref: Pan (2001)
QIC.pois.geese <- function(model.R, model.indep) {
  library(MASS)
  # Fitted and observed values for quasi likelihood
  mu.R <- model.R$fitted.values
      # alt: X <- model.matrix(model.R)
          #  names(model.R$coefficients) <- NULL
      #  beta.R <- model.R$coefficients
          #  mu.R <- exp(X %*% beta.R)
      y <- model.R$y

  # Quasi Likelihood for Poisson
  quasi.R <- sum((y*log(mu.R)) - mu.R) # poisson()$dev.resids - scale and weights = 1

  # Trace Term (penalty for model complexity)
  AIinverse <- ginv(model.indep$geese$vbeta.naiv) # Omega-hat(I) via Moore-Penrose generalized inverse of a matrix in MASS package
  # Alt: AIinverse <- solve(model.indep$geese$vbeta.naiv) # solve via indenity
  Vr <- model.R$geese$vbeta
  trace.R <- sum(diag(AIinverse %*% Vr))
  px <- length(mu.R) # number non-redunant columns in design matrix

  # QIC
  QIC <- (-2)*quasi.R + 2*trace.R
  QICu <- (-2)*quasi.R + 2*px    # Approximation assuming model structured correctly
  output <- c(QIC, QICu, quasi.R, trace.R, px)
  names(output) <- c('QIC', 'QICu', 'Quasi Lik', 'Trace', 'px')
  output
}

It would help if you documented this a little better. What are the inputs to this function? model.R and model.indep? — Ryan Simmons, Apr 28 '15 at 20:07
Any idea why the QIC calculation in your mini-package gives a very different value compared to the QIC() function in the MuMIn package? — commscho, Apr 26 '19 at 17:28

score 2 · Answer 3 · edited Feb 17 '15 at 03:05

2

You can use the model.sel command from the MuMIn package:

library(MuMIn)
model.sel(gee.0, gee.1, gee.2, gee.3, rank = QIC)

This uses MSE of prediction for model selection (Mean square error of prediction)--The smaller the better!

edited Feb 17 '15 at 03:05

Patrick Coulombe

2,536
6
20
28

answered Feb 17 '15 at 02:13

Keiv

21
1

I don't think this package works. When I try to use this command, it gives me the same QIC for every model, regardless of the correlation structure used. – Ryan Simmons Apr 28 '15 at 20:01
The `MuMIn` package seems to work fine (in 2019). But where can I find information on your claim about it using MSE for model selection? – Frederick Oct 16 '19 at 10:04

How to perform model selection in GEE in R

3 Answers3