1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
#' Display model fit summary of x or x like objects, fast
#'
#' Faster than the implementation in the arm package because it avoids refitting
#'
#' @title fastdisp: faster display of model summaries
#' @param x a model object
#' @param ... additional arguments to pass to \code{arm::\link[arm]{display}}
#' including number of digits
#' @details The time saving is only noticeable for large, time-consuming (g)lmer
#' fits.
#' @import arm
#' @return A printed summary of a x object
#' @examples
#' \donttest{
#' #Compare the time for displaying this modest model
#' require(arm)
#' m1 <- lmer(y ~ lectage + studage + (1|d) + (1|s), data=InstEval)
#' system.time(display(m1))
#' system.time(fastdisp(m1))
#' }
#' @seealso \code{\link[arm]{display}}
#' @rdname fastdisp
#' @export fastdisp
fastdisp <- function (x, ...) {
UseMethod("fastdisp", x)
}
#' @rdname fastdisp
#' @importFrom stats df pt
#' @export
fastdisp.merMod <- function (x, ...)
{
.local <- function (x, digits = 2, detail = FALSE)
{
out <- NULL
out$call <- x@call
print(out$call)
fcoef <- fixef(x)
useScale <- getME(x, "devcomp")$dims["useSc"]
corF <- lme4::vcov.merMod(x)@factors$correlation
coefs <- cbind(fcoef, corF@sd)
if (length(fcoef) > 0) {
if (!useScale) {
coefs <- coefs[, 1:2, drop = FALSE]
out$z.value <- coefs[, 1]/coefs[, 2]
out$p.value <- 2 * pnorm(abs(out$z.value), lower.tail = 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"]
cat("\nError terms:\n")
vc <- easyVarCorr(VarCorr(x), useScale = useScale,
digits)
print(vc[, c(1:2, 4:ncol(vc))], quote = FALSE)
out$ngrps <- lapply(x@flist, function(x) length(levels(x)))
is_REML <- isREML(x)
llik <- logLik(x, REML = is_REML)
out$AIC <- AIC(llik)
# out$deviance <- deviance(refitML(x))
out$n <- getME(x, "devcomp")$dims["n"]
# Dhat <- -2 * (llik)
# pD <- out$deviance - Dhat
# out$DIC <- out$deviance + pD
cat("---\n")
cat(sprintf("number of obs: %d, groups: ", out$n))
cat(paste(paste(names(out$ngrps), out$ngrps, sep = ", "),
collapse = "; "))
cat(sprintf("\nAIC = %g", round(out$AIC, 1)))
# cat(round(out$DIC, 1))
# cat("\ndeviance =", fround(out$deviance, 1), "\n")
if (useScale < 0) {
out$sigma.hat <- sigma(x)
cat("overdispersion parameter =", fround(out$sigma.hat,
1), "\n")
}
return(invisible(out))
}
.local(x, ...)
}
#' @rdname fastdisp
#' @export
fastdisp.merModList <- function(x, ...){
.local <- function (x, digits = 2, detail = FALSE)
{
out <- NULL
useScale <- getME(x[[1]], "devcomp")$dims["useSc"]
#useScale <- TRUE
out$call <- x[[1]]@call
listFE <- modelFixedEff(x)
row.names(listFE) <- listFE$term
out$t.value <- listFE$statistic
out$coef <- listFE$estimate
out$se <- listFE$std.error
listRE <- modelRandEffStats(x)
out$ngrps <- lapply(x[[1]]@flist, function(x) length(levels(x)))
is_REML <- isREML(x[[1]])
llik <- lapply(x, logLik, REML = is_REML)
out$AIC <- mean(unlist(lapply(llik, AIC)))
out$n <- round(mean(unlist(lapply(lapply(lapply(x, getME, "devcomp"),
"[[", "dims"), "[", 2))), 0) # round to nearest integer
print(out$call)
if (!detail) {
pfround(listFE[, 2:3], digits)
}
else {
listFE$p.value <- 2 * pt(abs(listFE$statistic), listFE$df, lower.tail = FALSE)
pfround(listFE[, 2:6], digits)
}
cat("\nError terms:\n")
vc <- easyVarCorr(VarCorr(x[[1]]), useScale = useScale,
digits)
# Resort the output of the random effect summary
listRE <- listRE[grep("cor_", listRE$term, invert=TRUE), ]
resid <- listRE[listRE$group == "Residual", ]
listRE <- listRE[listRE$group != "Residual", ]
listRE <- rbind(listRE, resid)
#
vc[, 3] <- as.character(round(listRE$estimate^2, digits = digits))
vc[, 4] <- as.character(round(listRE$estimate, digits = digits))
print(vc[, c(1:2, 4:ncol(vc))], quote = FALSE)
cat("---\n")
cat(sprintf("number of obs: %d, groups: ", out$n))
cat(paste(paste(names(out$ngrps), out$ngrps, sep = ", "),
collapse = "; "))
cat(sprintf("\nAIC = %g", round(out$AIC, 1)))
cat("---\n")
# cat(round(out$DIC, 1))
# cat("\ndeviance =", fround(out$deviance, 1), "\n")
if (useScale < 0) {
out$sigma.hat <- sigma(x)
cat("overdispersion parameter =", fround(out$sigma.hat,
1), "\n")
cat("---\n")
}
return(invisible(out))
}
.local(x, ...)
}
|
