Hi all,
This was posted originally on r-sig-mixed-models, but I thought I
would post here as well as it might be of more general interest.
With a colleague, I have been trying to implement the Conditional AIC
described by Vaida and Blanchard 2005 Biometrika, "Conditional Akaike
information for mixed-effects models". This quantity is derived in a
way analogous to the AIC, but is appropriate for scenarios where one
is interested in the particular coefficient estimates for individual
random effects. The formula for the asymptotic CAIC is given as
-2*log(likelihood of observed values, conditional on ML estimates of
fixed effects and empirical Bayes estimates of random effects) + 2*K
where K = rho + 1, and rho = "effective degrees of freedom" = trace
of the hat matrix mapping predicted values onto observed values.
After some thinking and some off-list advice, we have decided that
appropriate code for CAIC is
CAIC <- function(model) {
sigma <- attr(VarCorr(model), 'sc')
observed <- attr(model, 'y')
predicted <- fitted(model)
cond.loglik <- sum(dnorm(observed, predicted, sigma, log=TRUE))
rho <- hatTrace(model)
p <- length(fixef(model))
N <- nrow(attr(model, 'X'))
K.corr <- N*(N-p-1)*(rho+1)/((N-p)*(N-p-2)) + N*(p+1)/((N-p)*(N-p-2))
CAIC <- -2*cond.loglik + 2*K.corr
return(CAIC)
}
where K.corr is the finite-sample correction for K, for ML model fits.
I am posting this so that 1) This code can be of use to any other
souls in the statistical wilderness trying to do model selection with
mixed models, and 2) So that wiser minds can point out any errors in
our approach.
Thanks,
Kyle
