[R] Using Conditional AIC with lmer

[Kyle Edwards]

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