Nothing
#' @title lmerExp: transform the unit of coefficients (internal function)
#' @description Transform the unit of coefficients to "Coeff", "OR" or "RR"
#' @param lmer.coef lmer coefficients.
#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc..., Default: 'binomial'
#' @param dec Decimal point
#' @return The transforemed coefficients(95% CI), p-value
#' @details DETAILS
#' @examples
#' # EXAMPLE1
#' @rdname lmerExp
#' @importFrom stats pnorm
lmerExp <- function(lmer.coef, family = "binomial", dec) {
if (class(lmer.coef)[1] == "numeric") {
lmer.coef <- t(data.frame(lmer.coef))
}
pv <- 2 * (1 - pnorm(abs(lmer.coef[, 3])))
if (family == "binomial") {
OR <- paste(round(exp(lmer.coef[, 1]), dec), " (", round(exp(lmer.coef[, 1] - 1.96 * lmer.coef[, 2]), dec), ",", round(exp(lmer.coef[, 1] + 1.96 * lmer.coef[, 2]), dec), ")", sep = "")
return(cbind(OR, pv))
} else if (family %in% c(NULL, "gaussian")) {
coeff <- paste(round(lmer.coef[, 1], dec), " (", round(lmer.coef[, 1] - 1.96 * lmer.coef[, 2], dec), ",", round(lmer.coef[, 1] + 1.96 * lmer.coef[, 2], dec), ")", sep = "")
return(cbind(coeff, pv))
} else if (family %in% c("poisson", "quasipoisson")) {
RR <- paste(round(exp(lmer.coef[, 1]), dec), " (", round(exp(lmer.coef[, 1] - 1.96 * lmer.coef[, 2]), dec), ",", round(exp(lmer.coef[, 1] + 1.96 * lmer.coef[, 2]), dec), ")", sep = "")
return(cbind(RR, pv))
}
}
#' @title lmer.display: table for "lmerMod" or "glmerMod" object (lme4 package)
#' @description Make mixed effect model results from "lmerMod" or "glmerMod" object (lme4 package)
#' @param lmerMod.obj "lmerMod" or "glmerMod" object
#' @param dec Decimal, Default: 2
#' @param ci.ranef Show confidence interval of random effects?, Default: F
#' @return Table: fixed & random effect
#' @details DETAILS
#' @examples
#' library(geepack)
#' data(dietox)
#' dietox$Cu <- as.factor(dietox$Cu)
#' l1 <- lme4::lmer(Weight ~ Time + Cu + (1 | Pig) + (1 | Evit), data = dietox)
#' lmer.display(l1)
#' @rdname lmer.display
#' @export
#' @importFrom lme4 lmer glmer confint.merMod
#' @importFrom stats update formula
lmer.display <- function(lmerMod.obj, dec = 2, ci.ranef = F) {
sl <- summary(lmerMod.obj)
fixef <- sl$coefficients[-1, ]
mdata <- lmerMod.obj@frame
forms <- as.character(sl$call[[2]])
y <- forms[2]
xfr <- strsplit(forms[3], " \\+ ")[[1]]
xr <- xfr[grepl("\\|", xfr)]
xf <- xfr[!grepl("\\|", xfr)]
family.lmer <- ifelse(is.null(sl$family), "gaussian", sl$family)
uni.res <- ""
basemodel <- update(lmerMod.obj, stats::formula(paste(c(". ~ .", xf), collapse = " - ")), data = mdata)
unis <- lapply(xf, function(x) {
summary(stats::update(basemodel, stats::formula(paste0(". ~ . +", x)), data = mdata))$coefficients
})
unis2 <- Reduce(rbind, unis)
# uni.res <- unis2[rownames(unis2) != "(Intercept)", ]
adj.res <- sl$coefficients
adj.res <- adj.res[rownames(adj.res) != "(Intercept)", , drop = FALSE]
uni.res <- matrix(nrow = dim(adj.res)[1], ncol = dim(adj.res)[2])
rownames(uni.res) <- rownames(adj.res)
colnames(uni.res) <- colnames(adj.res)
for (i in 1:nrow(adj.res)) {
try(uni.res[rownames(adj.res)[i], ] <- unis2[rownames(adj.res)[i], ],
silent = TRUE
)
}
if (length(xf) == 1) {
fix.all <- lmerExp(uni.res, family = family.lmer, dec = dec)
family.label <- colnames(fix.all)[1]
colnames(fix.all) <- c(paste(family.label, "(95%CI)", sep = ""), "P value")
rownames(fix.all) <- rownames(sl$coefficients)[-1]
} else {
fix.all <- cbind(lmerExp(uni.res, family = family.lmer, dec = dec), lmerExp(fixef, family = family.lmer, dec = dec))
family.label <- colnames(fix.all)[1]
colnames(fix.all) <- c(paste("crude ", family.label, "(95%CI)", sep = ""), "crude P value", paste("adj. ", family.label, "(95%CI)", sep = ""), "adj. P value")
rownames(fix.all) <- rownames(sl$coefficients)[-1]
}
ranef <- data.frame(sl$varcor)[, c(1, 4)]
ranef.out <- cbind(as.character(round(ranef[, 2], dec)), matrix(NA, nrow(ranef), ncol(fix.all) - 1))
rownames(ranef.out) <- ranef[, 1]
if (ci.ranef) {
ranef.ci <- confint.merMod(lmerMod.obj, parm = 1:nrow(ranef), oldNames = F)^2
ranef.paste <- paste(round(ranef$vcov, dec), " (", round(ranef.ci[, 1], dec), ",", round(ranef.ci[, 2], dec), ")", sep = "")
ranef.out <- cbind(ranef.paste, matrix(NA, nrow(ranef), 3))
rownames(ranef.out) <- ranef[, 1]
}
ranef.na <- rbind(rep(NA, ncol(fix.all)), ranef.out)
rownames(ranef.na)[1] <- "Random effects"
colnames(ranef.na) <- colnames(fix.all)
## rownames
rn.uni <- lapply(unis, rownames)
fix.all.list <- lapply(1:length(xf), function(x) {
if (grepl(":", xf[x])) {
a <- unlist(strsplit(xf[x], ":"))[1]
b <- unlist(strsplit(xf[x], ":"))[2]
fix.all[grepl(a, rownames(fix.all)) & grepl(b, rownames(fix.all)), ]
} else {
fix.all[rownames(fix.all) %in% rn.uni[[x]], ]
}
# fix.all[grepl(x, rownames(fix.all)), ]
})
varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all))
lapply(varnum.mfac, function(x) {
fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])
})
fix.all.unlist <- Reduce(rbind, fix.all.list)
rn.list <- lapply(1:length(xf), function(x) {
if (grepl(":", xf[x])) {
a <- unlist(strsplit(xf[x], ":"))[1]
b <- unlist(strsplit(xf[x], ":"))[2]
rownames(fix.all)[grepl(a, rownames(fix.all)) & grepl(b, rownames(fix.all))]
} else {
rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]
}
# rownames(fix.all)[grepl(x, rownames(fix.all))]
})
varnum.2fac <- which(lapply(xf, function(x) {
length(sapply(mdata, levels)[[x]])
}) == 2)
lapply(varnum.2fac, function(x) {
rn.list[[x]] <<- paste(xf[x], ": ", levels(mdata[, xf[x]])[2], " vs ", levels(mdata[, xf[x]])[1], sep = "")
})
lapply(varnum.mfac, function(x) {
if (grepl(":", xf[x])) {
a <- unlist(strsplit(xf[x], ":"))[1]
b <- unlist(strsplit(xf[x], ":"))[2]
if (a %in% xf && b %in% xf) {
ref <- paste0(a, levels(mdata[, a])[1], ":", b, levels(mdata[, b])[1])
rn.list[[x]] <<- c(paste(xf[x], ": ref.=", ref, sep = ""), gsub(xf[x], " ", rn.list[[x]]))
} else {
rn.list[[x]] <<- c(paste(xf[x], ": ref.=NA", sep = ""), gsub(xf[x], " ", rn.list[[x]]))
}
} else {
rn.list[[x]] <<- c(paste(xf[x], ": ref.=", levels(mdata[, xf[x]])[1], sep = ""), gsub(xf[x], " ", rn.list[[x]]))
}
# rn.list[[x]] <<- c(paste(xf[x], ": ref.=", levels(mdata[, xf[x]])[1], sep = ""), gsub(xf[x], " ", rn.list[[x]]))
})
if (class(fix.all.unlist)[1] == "character") {
fix.all.unlist <- t(data.frame(fix.all.unlist))
}
rownames(fix.all.unlist) <- unlist(rn.list)
tb.df <- as.data.frame(rbind(fix.all.unlist, ranef.na))
## metric
no.grp <- sl$ngrps
no.obs <- nrow(mdata)
ll <- round(sl$logLik[[1]], dec)
aic <- round(sl$AICtab, dec)[1]
met <- c(NA, no.grp, no.obs, ll, aic)
met.mat <- cbind(met, matrix(NA, length(met), ncol(fix.all) - 1))
rownames(met.mat) <- c("Metrics", paste("No. of groups (", rownames(met.mat)[2:(length(no.grp) + 1)], ")", sep = ""), "No. of observations", "Log-likelihood", "AIC value")
met.mat[, 1] <- as.character(met.mat[, 1])
colnames(met.mat) <- colnames(tb.df)
tb.lmerMod <- rbind(tb.df, met.mat)
lapply(seq(2, ncol(fix.all), by = 2), function(x) {
tb.lmerMod[, x] <<- as.numeric(as.vector(tb.lmerMod[, x]))
})
## caption
cap.lmerMod <- paste(sl$methTitle, " : ", y, " ~ ", forms[3], sep = "")
return(list(table = tb.lmerMod, caption = cap.lmerMod))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.