Use apsrtable with `merMod` objects

[jknowles]

I'm attempting to update this package to construct tables for mixed effect models fit of the class merMod. Below is a simple example and the corresponding error message. I am wondering if there is some way I can work around the fundamental reliance on the $ operator so that S4 objects like merMod objects can be easily adapted by creating a summary function.

(fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))

# adapted from code in package already
apsrtableSummary.merMod <- function (object, ...) {
  obj <- summary(object)
  fcoef <- coef(obj)
  out <- list()
  useScale <- obj$useScale
  corF <- vcov(object)@factors$correlation
  coefs <- cbind(fcoef[, 1:2])
  if (length (fcoef) > 0){
      if (obj$useScale == FALSE) {
        coefs <- coefs[, 1:2, drop = FALSE]
        out$z.value <- coefs[, 1]/coefs[, 2]
        out$p.value <- 2 * pnorm(abs(out$z.value), lower = FALSE)
        coefs <- cbind(coefs,
                       `z value` = out$z.value,
                       `Pr(>|z|)` = out$p.value)
      }
      else {
        out$t.value <- coefs[, 1]/coefs[, 2]
        coefs <- cbind(coefs, `t value` = out$t.value)
      }
    dimnames(coefs)[[2]][1:2] <- c("coef.est", "coef.se")
#       if(detail){
#         pfround (coefs, digits)
#       }
#       else{
#         pfround(coefs[,1:2], digits)
#     }
  }
  out$coef <- coefs[,"coef.est"]
  out$se <- coefs[,"coef.se"]
#   vc <- as.matrix(VarCorr(object, digits))
#   vc[,1] <-
#   print (vc[,c(1:2,4:ncol(vc))], quote=FALSE)

  out$ngrps <- lapply(object@flist, function(x) length(levels(x)))
  ## Model fit statistics.
  ll <- logLik(object)[1]
  deviance <- deviance(object)
  AIC <- AIC(object)
  BIC <- BIC(object)
  N <- as.numeric(length(obj$residuals))
  G <- as.numeric(obj$ngrps)
  sumstat <- c(logLik = ll, deviance = deviance, AIC = AIC,
               BIC = BIC, N = N, Groups = G)

  ## Return model summary.
  list(coef = obj$coefficients, sumstat = sumstat,
       contrasts = attr(model.matrix(object), "contrasts"),
       xlevels = NULL, call = object@call)
}

apsrtableSummary(fm1)

Results in :

$coef
             Estimate Std. Error   t value
(Intercept) 251.40510   6.824556 36.838311
Days         10.46729   1.545789  6.771485

$sumstat
   logLik  deviance       AIC       BIC         N    Groups 
-871.8141 1743.6283 1755.6283 1774.7860  180.0000   18.0000 

$contrasts
NULL

$xlevels
NULL

$call
lmer(formula = Reaction ~ Days + (Days | Subject), data = sleepstudy)

Which I think is OK, but when I move on to apsrtable the object with apsrtable(fm1), I get:

Error in x$se : $ operator not defined for this S4 class

The traceback sheds some light on where things are not going well:

traceback()
3: FUN(X[[1L]], ...)
2: lapply(models, function(x) {
       s <- try(apsrtableSummary(x), silent = TRUE)
       if (inherits(s, "try-error")) {
           s <- summary(x)
       }
       if (!is.null(x$se) && se != "vcov") {
           est <- coef(x)
           if (class(x$se) == "matrix") {
               x$se <- sqrt(diag(x$se))
           }
           s$coefficients[, 3] <- tval <- est/x$se
           e <- try(s$coefficients[, 4] <- 2 * pt(abs(tval), length(x$residuals) - 
               x$rank, lower.tail = FALSE), silent = TRUE)
           if (inherits(e, "try-error")) {
               s$coefficients[, 4] <- 2 * pnorm(abs(tval), lower.tail = FALSE)
           }
           s$se <- x$se
       }
       if (se == "pval") {
           s$coefficients[, 2] <- s$coefficients[, 4]
       }
       return(s)
   })
1: apsrtable(fm1)

Somehow, the reliance on x$se is a problem. What am I doing wrong to allow getting around this step and moving to creating the summary.

[malecki]

I will look in more depth someday this week.

Yeah, using $ has a bad smell to it. It was convenient for lm & glm, and I think that's something worth preserving.

I can think of a couple ways around this. The internal models and summaries should each be a class that can then just have a se getter that 'does the right thing' for s3, custom-se, or any other objects.

Thanks for the good work! And like I said I'll try to look at it myself in the coming week.

On Aug 11, 2014, at 23:38, Jared Knowles notifications@github.com wrote:

I'm attempting to update this package to construct tables for mixed effect models fit of the class merMod. Below is a simple example and the corresponding error message. I am wondering if there is some way I can work around the fundamental reliance on the $ operator so that S4 objects like merMod objects can be easily adapted by creating a summary function.

(fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))

adapted from code in package already

apsrtableSummary.merMod <- function (object, ...) {
obj <- summary(object)
fcoef <- coef(obj)
out <- list()
useScale <- obj$useScale
corF <- vcov(object)@factors$correlation
coefs <- cbind(fcoef[, 1:2])
if (length (fcoef) > 0){
if (obj$useScale == FALSE) {
coefs <- coefs[, 1:2, drop = FALSE]
out$z.value <- coefs[, 1]/coefs[, 2]
out$p.value <- 2 * pnorm(abs(out$z.value), lower = FALSE)
coefs <- cbind(coefs,
z value = out$z.value,
Pr(>|z|) = out$p.value)
}
else {
out$t.value <- coefs[, 1]/coefs[, 2]
coefs <- cbind(coefs, t value = out$t.value)
}
dimnames(coefs)[[2]][1:2] <- c("coef.est", "coef.se")

if(detail){

pfround (coefs, digits)

}

else{

pfround(coefs[,1:2], digits)

}

}
out$coef <- coefs[,"coef.est"]
out$se <- coefs[,"coef.se"]

vc <- as.matrix(VarCorr(object, digits))

vc[,1] <-

print (vc[,c(1:2,4:ncol(vc))], quote=FALSE)

out$ngrps <- lapply(object@flist, function(x) length(levels(x)))

Model fit statistics.

ll <- logLik(object)[1]
deviance <- deviance(object)
AIC <- AIC(object)
BIC <- BIC(object)
N <- as.numeric(length(obj$residuals))
G <- as.numeric(obj$ngrps)
sumstat <- c(logLik = ll, deviance = deviance, AIC = AIC,
BIC = BIC, N = N, Groups = G)

Return model summary.

list(coef = obj$coefficients, sumstat = sumstat,
contrasts = attr(model.matrix(object), "contrasts"),
xlevels = NULL, call = object@call)
}

apsrtableSummary(fm1)
Results in :

$coef
Estimate Std. Error t value
(Intercept) 251.40510 6.824556 36.838311
Days 10.46729 1.545789 6.771485

$sumstat
logLik deviance AIC BIC N Groups
-871.8141 1743.6283 1755.6283 1774.7860 180.0000 18.0000

$contrasts
NULL

$xlevels
NULL

$call
lmer(formula = Reaction ~ Days + (Days | Subject), data = sleepstudy)

Which I think is OK, but when I move on to apsrtable the object with apsrtable(fm1), I get:

Error in x$se : $ operator not defined for this S4 class
The traceback sheds some light on where things are not going well:

traceback()
3: FUN(X[[1L]], ...)
2: lapply(models, function(x) {
s <- try(apsrtableSummary(x), silent = TRUE)
if (inherits(s, "try-error")) {
s <- summary(x)
}
if (!is.null(x$se) && se != "vcov") {
est <- coef(x)
if (class(x$se) == "matrix") {
x$se <- sqrt(diag(x$se))
}
s$coefficients[, 3] <- tval <- est/x$se
e <- try(s$coefficients[, 4] <- 2 * pt(abs(tval), length(x$residuals) -
x$rank, lower.tail = FALSE), silent = TRUE)
if (inherits(e, "try-error")) {
s$coefficients[, 4] <- 2 * pnorm(abs(tval), lower.tail = FALSE)
}
s$se <- x$se
}
if (se == "pval") {
s$coefficients[, 2] <- s$coefficients[, 4]
}
return(s)
})
1: apsrtable(fm1)
Somehow, the reliance on x$se is a problem. What am I doing wrong to allow getting around this step and moving to creating the summary.

—
Reply to this email directly or view it on GitHub.

[jknowles]

@malecki Any update on this? No worries if there is not. I'm happy to help. I have been using xtable as an alternative, but it only gives the coefficients and model reported standard errors in its summary, and then doesn't play nice combining this with other model types.

[malecki]

I broke something else, it sorta works for merMods now... to be continued. (checkout dev branch)

[jknowles]

@malecki I'll check it out and see if I can help contribute. I also wrote an extension that can summarize regression objects from the rms package. Maybe I'll contribute that as well. Those are easy since they are just lm objects with additional slots.

[jknowles]

@malecki -- Just curious, do you think PR #2 is worth pulling in?

[malecki]

yes, just haven't gotten a chance to look in more detail (sorry, been
really busy)

On Fri, Sep 12, 2014 at 5:02 PM, Jared Knowles notifications@github.com
wrote:

@malecki https://github.com/malecki -- Just curious, do you think PR #2
#2 is worth pulling in?

—
Reply to this email directly or view it on GitHub
#1 (comment).

[jknowles]

No problem, just checking in!

[jknowles]

I'd say this is fixed on dev' now, but nowdev' needs to be rebuilt to pass R CMD CHECK. I'll work on this slowly, but I'm going to open a new issue to indicate that is the next step.