[R-sig-ME] mtable (package: memsic) with lmer objects
David Atkins
datkins at u.washington.edu
Mon Mar 16 19:31:56 CET 2009
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
More information about the R-sig-mixed-models
mailing list