QIC.binom.geeglm<-function(model.geeglm,model.independence) { #calculates binomial QAIC of Pan (2001) #obtain trace term of QAIC AIinverse<-solve(model.independence$geese$vbeta.naiv) V.msR<-model.geeglm$geese$vbeta trace.term<-sum(diag(AIinverse%*%V.msR)) #estimated mean and observed values mu.R<-model.geeglm$fitted.values y<-model.geeglm$y #scale for binary data scale <- 1 #quasilikelihood for binomial model quasi.R<-sum(y*log(mu.R/(1-mu.R))+log(1-mu.R))/scale QIC<-(-2)*quasi.R + 2*trace.term output <- c(QIC,trace.term) names(output)<-c('QIC','CIC') output }