[R-sig-ME] mtable (package: memsic) with lmer objects

[David Atkins]

lmer useRs--

Not too long ago, I discovered the mtable() function in the memisc 
package, which formats tables of regression results for one or more 
modeles in "friendly" ways.  At least, it formats them in ways that I 
often see in psychology, psychiatry, and sociology journals (where I 
usually read/publish).  There is an example below.

With some assistance from Martin Elff, the package author, I now have a 
summary method that works with lmer() and glmer(), though admittedly, it 
is pretty bare bones and could be extended in useful ways.

That is part of the reason why I'm posting this.  It could be 
interesting to have options for including random-effects variances, or 
group sizes, but these will vary depending on the model and gets a bit 
beyond my primitive coding skills.

So, here's a plain-jane version for those who might find it useful, and 
for those who might want to take a crack at extending (though could you 
let me know if you do?  I think Martin may role this into the next 
memisc update).  [BTW, this is based off the getSummary.glm code, which 
explains most of the commented out code.]

cheers, Dave

library(lme4)
library(memisc)

### create three models
fm1 <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy)
fm1.1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1.2 <- lmer(Reaction ~ as.factor(Days) + (Days|Subject), sleepstudy)

### note: need to run the code below fro setCoefTemplate and
### getSummary.lmer first

mtable("Model 1"=fm1, "Model 2"=fm1.1, "Model 3"=fm1.2,
     			coef.style = "est.ci", # using "homegrown" est.ci, specified above
     			summary.stats=c("AIC","BIC"),
     			getSummary = "getSummary.lmer")#,

setCoefTemplate(
  est.ci=c(
	est = "($est:#)($p:*)",
	ci = "[($lwr:#),($upr:#)]"))

getSummary.lmer <- function (obj, alpha = 0.05, ...)
{
     require(lme4)
     smry <- summary(obj)
     #N <- if (length(weights(obj))) ### NOTE: how to deal with 
groups/samp size?
     #    sum(weights(obj))
     #else sum(smry$df[1:2])
     coef <- smry at coefs
     lower <- qnorm(p = alpha/2, mean = coef[, 1], sd = coef[,
         2])
     upper <- qnorm(p = 1 - alpha/2, mean = coef[, 1], sd = coef[,
         2])
     if (ncol(smry at coefs) == 3) {
     	p <- (1 - pnorm(smry at coefs[,3]))*2 # NOTE: no p-values for lmer() 
due to
     											# unclear dfs; calculate p-values based on z
     	coef <- cbind(coef, p, lower, upper)
    	} else {
    			coef <- cbind(coef, lower, upper) # glmer will have 4 columns 
with p-values
    			}
     colnames(coef) <- c("est", "se", "stat", "p", "lwr", "upr")
     #phi <- smry$dispersion
     #LR <- smry$null.deviance - smry$deviance
     #df <- smry$df.null - smry$df.residual
     ll <- smry at AICtab[3][,1]
     deviance <- smry at AICtab[4][,1]
     #if (df > 0) {
     #    p <- pchisq(LR, df, lower.tail = FALSE)
     #    L0.pwr <- exp(-smry$null.deviance/N)
     #    McFadden <- 1 - smry$deviance/smry$null.deviance
     #    Cox.Snell <- 1 - exp(-LR/N)
     #    Nagelkerke <- Cox.Snell/(1 - L0.pwr)
     #}
     #else {
     #    LR <- NA
     #    df <- NA
     #    p <- NA
     #    McFadden <- NA
     #    Cox.Snell <- NA
     #    Nagelkerke <- NA
     #}
     AIC <- smry at AICtab[1][,1] # NOTE: these are both data.frames? not 
sure why...
     BIC <- smry at AICtab[2][,1]
     ### NOTE: don't see a similar slot for "xlevels" to get levels of
     ###		factor variables used as predictors; for time being, force
     ###		user to specify explicitly; nope that didn't work...
     #if (fac != NULL) {
     # 	n <- length(fac)
     #	xlevels <- vector(n, mode = "list")
     #	for (i in 1:n) {
     #		xlevels[i] <- levels(obj at frame[,fac[i]])
     #		}
     #	}
     #sumstat <- c(phi = phi, LR = LR, df = df, p = p, logLik = ll,
     #    deviance = deviance, McFadden = McFadden, Cox.Snell = Cox.Snell,
     #    Nagelkerke = Nagelkerke, AIC = AIC, BIC = BIC, N = N)
     sumstat <- c(logLik = ll, deviance = deviance, AIC = AIC, BIC = BIC)
     list(coef = coef, sumstat = sumstat,
     	 contrasts = attr(model.matrix(obj), "contrasts"),
         xlevels = NULL, call = obj at call)
}

-- 
Dave Atkins, PhD
Research Associate Professor
Center for the Study of Health and Risk Behaviors
Department of  Psychiatry and Behavioral Science
1100 NE 45th Street, Suite 300
Seattle, WA  98105
206-616-3879
datkins at u.washington.edu