Serious Coding Error in QIC function in geepack? [closed]

I believe the QIC function in geepack has a significant error. The function appears to incorrectly specify the independence model, which is needed to calculate QIC. The function will therefore often if not typically produce an incorrect QIC score. In many cases, the discrepancy between the real QIC score and the QIC score generated by the QIC function in geepack will be statistically significant. Am I missing something?

One can view the source code for QIC() by typing geepack:::QIC.geeglm into the R console. I have reproduced that code below. 

function (object, tol = .Machine$double.eps, ...) { if (!("geeglm" %in% class(object))) { stop("QIC requires a geeglm object as input") } invert <- if ("MASS" %in% loadedNamespaces()) { MASS::ginv } else { solve } computeqic <- function(object) { mu <- object$fitted.values y <- object$y type <- family(object)$family quasi <- switch(type, poisson = sum((y * log(mu)) - mu), gaussian = sum(((y - mu)^2)/-2), binomial = sum(y * log(mu/(1 - mu)) + log(1 - mu)), Gamma = sum(-y/(mu - log(mu))), stop("Error: distribution not recognized")) object$call$corstr <- "independence" object$call$zcor <- NULL model.indep <- eval(object, parent.frame()) AIinverse <- invert(model.indep$geese$vbeta.naiv, tol = tol) Vr <- object$geese$vbeta trace <- sum(diag(AIinverse %*% Vr)) params <- length(coef(object)) kpm <- params + length(object$geese$alpha) QIC <- -2 * (quasi - trace) QICu <- -2 * (quasi - params) QICC <- QIC + (2 * kpm * (kpm + 1))/(length(object$residuals) - kpm - 1) output <- c(QIC, QICu, quasi, trace, params, QICC) names(output) <- c("QIC", "QICu", "Quasi Lik", "CIC", "params", "QICC") output } if (length(list(...))) { results <- lapply(list(object, ...), computeqic) check <- sapply(list(object, ...), function(x) { length(x$y) }) if (any(check != check[1])) warning("models are not all fitted to the same number of observations") res <- do.call("rbind", results) Call <- match.call() Call$k <- NULL row.names(res) <- as.character(Call[-1L]) res } else { computeqic(object) } } <bytecode: 0x16818dc10> <environment: namespace:geepack> 

In the above function, the input (i.e., "object") is the model whose QIC you wish to calculate. Call this the "target model." Line 29 is where the target model's QIC is calculated: 

QIC <- -2 * (quasi - trace).

"quasi" is the quasi-likelihood of the target model (line 17). "trace" is the sum of the diagonal elements of a matrix product (line 25). That matrix product is the result of multiplying the inverse of the naive variance-covariance matrix of the independence model (line 24) by the robust variance-covariance matrix of the target model (line 25). The independence model is created in lines 21-22. The independence model is supposed to be the analog of the target model, except with an independence correlation structure. 

The QIC function is written in such a way that it uses the target model as the basis to create the independence model, but the code does not accomplish this task. The error occurs on line 21. It appears the intent was to take the target model (i.e., "object") and change it into an independence model by changing the correlation structure to "independence" (line 21):

object$call$corstr <- "independence")

The code then evaluates this object and calls it "model.indep" (line 23): 

model.indep <- eval(object, parent.frame()))

However, line 21 does not in fact change the correlation structure of the target model to an independence correlation structure; line 21 simply changes the name of the original correlation structure to "independence." The result is that model.indep is not an independence model, and the QIC calculation in line 29 is therefore (probably) incorrect.

As proof of this, we can use the data provided from the fixed2Zcor documentation in geepack. (While I will use this data in what's below, the errors do not apply only to models with fixed correlation structures. Indeed, I won't run a model with a fixed correlation structure.):

timeorder <- rep(1:5, 6) tvar <- timeorder + rnorm(length(timeorder)) idvar <- rep(1:6, each=5) uuu <- rep(rnorm(6), each=5) yvar <- 1 + 2*tvar + uuu + rnorm(length(tvar)) simdat <- data.frame(idvar, timeorder, tvar, yvar) head(simdat,12) simdatPerm <- simdat[sample(nrow(simdat)),] simdatPerm <- simdatPerm[order(simdatPerm$idvar),] head(simdatPerm) 

We can fit the following model to this data, which has an exchangeable correlation structure:

mod_exch <- geeglm (yvar~tvar, id=idvar, data=simdatPerm, corstr="exchangeable")

Here's the model summary: 

> print(summary(mod_exch),10) Call: geeglm(formula = yvar ~ tvar, data = simdatPerm, id = idvar,     corstr = "independence")  Coefficients:                 Estimate      Std.err      Wald Pr(>|W|)     (Intercept) 0.9007246379 0.4303760528   4.38014  0.03636 *   tvar        2.0537574814 0.1185602549 300.06827  < 2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation structure = exchangeable Estimated Scale Parameters:                Estimate      Std.err (Intercept) 2.541291877 0.6195246062   Link = identity Estimated Correlation Parameters:           Estimate       Std.err alpha 0.7091991363 0.08091013334 Number of clusters:   6  Maximum cluster size: 5  

Running geepack's QIC() function, we get mod_exch's QIC score:

> QIC(mod_exch)[1] QIC 79.91151102

Below is what leads me to conclude the above QIC score is incorrect. We can extract the relevant part of the QIC function from geepack and run it directly in R, with mod_exch taking the place of "object." To do so, we first have to load library(MASS) and set the value of tol to .Machine$double.eps. We also need to change "invert" (from their code) to ginv (see line 11 of the below code).

mu <- mod_exch$fitted.values y <- mod_exch$y type <- family(mod_exch)$family quasi <- switch(type, poisson = sum((y * log(mu)) - mu), gaussian = sum(((y - mu)^2)/-2), binomial = sum(y * log(mu/(1 - mu)) + log(1 - mu)), Gamma = sum(-y/(mu - log(mu))), stop("Error: distribution not recognized")) mod_exch$call$corstr <- "independence" mod_exch$call$zcor <- NULL model.indep <- eval(mod_exch, parent.frame()) AIinverse <- ginv(model.indep$geese$vbeta.naiv, tol = tol) Vr <- mod_exch$geese$vbeta trace <- sum(diag(AIinverse %*% Vr)) params <- length(coef(mod_exch)) kpm <- params + length(mod_exch$geese$alpha) QIC <- -2 * (quasi - trace) 

Running this code then typing in "QIC" in the command line, we see that the QIC is equal to the QIC score we got from the QIC() function in geepack.

> QIC[1] 79.91151102

In other words, the code we extracted and ran directly in R is doing the same thing as the code in the QIC function in geepack. Good.

Here comes the problem. Let's run a summary of model.indep, which, recall, is supposed to be an independence model. 

> print(summary(model.indep),10) Call: geeglm(formula = yvar ~ tvar, data = simdatPerm, id = idvar, corstr = "independence") Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) 0.9007246379 0.4303760528 4.38014 0.03636 * tvar 2.0537574814 0.1185602549 300.06827 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation structure = exchangeable Estimated Scale Parameters: Estimate Std.err (Intercept) 2.541291877 0.6195246062 Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.7091991363 0.08091013334 Number of clusters: 6 Maximum cluster size: 5 

The call says model.indep has an "independence" correlation structure, but that is simply because we had changed the name of the correlation structure of mod_exch to "independence." In fact, model.indep is not actually a model with an independence correlation structure; it is just our original exchangeable model, mod_exch. How do we know this? You can see at the bottom of the summary of model.indep that it fits a correlation estimate of alpha; independence models don't do this in geeglm. Moreover, the value of alpha is .709, which is the same value of alpha for the exchangeable model. Also, about halfway down the model summary, you'll see that it says ``Correlation structure = exchangeable".

So I'm fairly confident the QIC function in geepack is wrong. It does not correctly specify an independence model, and therefore the QIC score will generally be incorrect. In the above example, I believe the correct QIC score is 84.438, not 79.912. This is a non-trivial difference if we adopt the convention from AIC that a 2-point difference in AIC scores is statistically significant.  

Thoughts?