How to perform model selection in GEE in R 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
 A: 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.  
A: [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
}

A: 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!
