#' Get p-values for a multilevel regression model
#'
#' Takes a model fit with lmer and (optionally) an estimation method, returns a
#' data frame containing fixed effect coefficient estimates, SE, t-values, and
#' p-values. Can use normal approximation, Kenward-Roger approximation for
#' degrees of freedom, or both. If Kenward-Roger approximation is used, result
#' also contains approximated degrees of freedom.
#'
#' @param model a merMod model object fit by lmer
#' @param method p-value estimation method, defaults to "normal", can also be
#' "KR" or "all"
#' @return data frame containing fixed effect estimates, SE, t-values, (degrees
#' of freedom if appropriate) and p-values.
#' @export
#' @examples
#' require(lme4)
#' model <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#' get_pvalues(model)
#' get_pvalues(model, "KR")
get_pvalues <- function(model, method = "normal") {
# Extract fixed effect coefficients from model
coefs <- data.frame(coef(summary(model)))
# Check/warn if bad method name used
supported <- c("normal", "KR", "all")
if (!is.element(method, supported)) {
warning("Invalid p-value estimation method. Use 'normal', 'KR'",
"or 'all'. Defaulting to normal approximation.")
method <- "normal"
}
# Apply normal approx if asked for
if (method %in% c("normal", "all")) {
coefs$p.normal <- get_pvalues_normal(coefs$t.value)
coefs$p.normal.star <- get_p_stars(coefs$p.normal)
}
# Apply KR if asked for
if (method %in% c("KR", "all")) {
#run a minimal KRmodcomp
restrictionMatrix <- cbind(t(rep(0, length(fixef(model))-1)), 1)
kr <- pbkrtest::KRmodcomp(model, restrictionMatrix)
#get ddf
coefs$df.KR <- pbkrtest::getKR(kr, "ddf")
coefs$p.KR <- get_pvalues_kr(coefs$t.value, coefs$df.KR)
coefs$p.KR.star <- get_p_stars(coefs$p.KR)
}
coefs
}
get_pvalues_normal <- function(ts) {
2 * (1 - pnorm(abs(ts)))
}
get_pvalues_kr <- function(ts, dfs) {
2 * (1 - pt(abs(ts), dfs))
}
get_p_stars <- function(ps) {
symnum(ps, na = FALSE, cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", " "), legend = FALSE)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.