#Adjusted R-square after McElree, using fixed effects and random effects (level 1)
FRsquared <- function(model, observed) {
SSR <- sum(residuals(model,level=1)^2)
SS <- sum((observed - mean(observed))^2)
N <- model$dims$N
fixedparms <- length(model$coefficients$fixed)
randparms <- model$dims$qvec[1]
parms = fixedparms + randparms + randparms*(randparms-1)/2 + 1
#means + variances + covariances
rsq = 1 - (SSR/(N-parms))/(SS/(N-1)) 
}