Nothing
#############################################################################
# Copyright (c) 2013-2020 Alexandra Kuznetsova, Per Bruun Brockhoff, and
# Rune Haubo Bojesen Christensen
#
# This file is part of the lmerTest package for R (*lmerTest*)
#
# *lmerTest* is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# *lmerTest* is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# <https://www.r-project.org/Licenses/> and/or
# <http://www.gnu.org/licenses/>.
#############################################################################
#
# step.R - implementation of backward elimination for lmerModLmerTest objects
# ------- Contents: --------
#
# --- Generics: ---
#
# step
# get_model
#
# --- methods: ---
#
# step.lmerModLmerTest
# step.default
# get_model.step_list
# print.step_list
# plot.step_list
#
# --- other exported function: ---
#
# --- utility functions: ---
#
# ran_redTable
# fix_redTable
# reduce_random
# ranova_lm
# reduce_fixed
#
##############################################
######## step()
##############################################
#' Generic Step Function
#'
#' Generic step function with default method \code{stats::step}. This
#' construction ensures that \code{stats::step} still works on \code{lm}
#' objects etc. after loading the \pkg{lmerTest} package.
#'
#' @param object a model object.
#' @param ... currently not used.
#'
#' @author Rune Haubo B. Christensen
#' @seealso \code{\link[=step.lmerModLmerTest]{step}}
#' @export
#' @keywords internal
step <- function(object, ...) UseMethod("step")
##############################################
######## step.default()
##############################################
#' @rdname step
#' @export
#' @keywords internal
step.default <- function(object, ...) stats::step(object, ...)
##############################################
######## step.lmerModLmerTest()
##############################################
#' Backward Elimination for Linear Mixed Models
#'
#' Backward elimination of random-effect terms followed by backward elimination
#' of fixed-effect terms in linear mixed models.
#'
#' Tests of random-effects are performed using \code{\link{ranova}} (using
#' \code{reduce.terms = TRUE}) and tests of fixed-effects are performed using
#' \code{\link[=drop1.lmerModLmerTest]{drop1}}.
#'
#' The step method for \code{\link{lmer}} fits has a print method.
#'
#' @param object a fitted model object. For the \code{lmerModLmerTest} method
#' an \code{\link{lmer}} model fit (of class \code{"lmerModLmerTest"}.)
#' @param ddf the method for computing the denominator degrees of freedom and
#' F-statistics. \code{ddf="Satterthwaite"} (default) uses Satterthwaite's method;
#' \code{ddf="Kenward-Roger"} uses Kenward-Roger's method.
#' @param alpha.random alpha for random effects elimination
#' @param alpha.fixed alpha for fixed effects elimination
#' @param reduce.fixed reduce fixed effect structure? \code{TRUE} by default.
#' @param reduce.random reduce random effect structure? \code{TRUE} by default.
#' @param keep an optional character vector of fixed effect terms which should
#' not be considered for eliminated. Valid terms are given by
#' \code{attr(terms(object), "term.labels")}. Terms that are marginal to terms
#' in keep will also not be considered for eliminations.
#' @param ... currently not used.
#'
#' @return \code{step} returns a list with elements \code{"random"} and
#' \code{"fixed"} each
#' containing anova-like elimination tables. The \code{"fixed"} table is
#' based on \code{drop1} and the \code{"random"} table is
#' based on \code{ranova} (a \code{drop1}-like table for random effects). Both
#' tables have a column \code{"Eliminated"} indicating the order in which terms
#' are eliminated from the model with zero (\code{0}) indicating that the term
#' is not eliminated from the model.
#'
#' The \code{step} object also contains the final model as an attribute which
#' is extractable with \code{get_model(<step_object>)}.
#' @seealso \code{\link[=drop1.lmerModLmerTest]{drop1}} for tests of marginal
#' fixed-effect terms and \code{\link{ranova}} for a
#' \code{\link[=drop1.lmerModLmerTest]{drop1}}-like table of reduction of
#' random-effect terms.
#' @author Rune Haubo B. Christensen and Alexandra Kuznetsova
#' @export
#' @examples
#'
#' # Fit a model to the ham dataset:
#' fm <- lmer(Informed.liking ~ Product*Information+
#' (1|Consumer) + (1|Product:Consumer)
#' + (1|Information:Consumer), data=ham)
#'
#' # Backward elimination using terms with default alpha-levels:
#' (step_res <- step(fm))
#' final <- get_model(step_res)
#' anova(final)
#'
#' \dontrun{
#' # Fit 'big' model:
#' fm <- lmer(Informed.liking ~ Product*Information*Gender*Age +
#' + (1|Consumer) + (1|Consumer:Product) +
#' (1|Consumer:Information), data=ham)
#' step_fm <- step(fm)
#' step_fm # Display elimination results
#' final_fm <- get_model(step_fm)
#' }
#'
step.lmerModLmerTest <- function(object, ddf=c("Satterthwaite", "Kenward-Roger"),
alpha.random=0.1, alpha.fixed=0.05,
reduce.fixed=TRUE, reduce.random=TRUE,
keep, ...) {
# Check for and warn about deprecated arguments:
ignored <- c("type", "fixed.calc", "lsmeans.calc", "difflsmeans.calc",
"test.effs")
dots <- list(...)
for(nm in ignored) if(any(pmatch(names(dots), nm, nomatch = 0)))
warning(paste0("Argument '", nm, "' is deprecated and ignored."))
if(any(pmatch(names(dots), "keep.effs", nomatch = 0)))
warning("Argument 'keep.effs' is deprecated: use 'keep' instead")
# reduce random and fixed parts?
if(!reduce.random) alpha.random <- 1
if(!reduce.fixed) alpha.fixed <- 1
if(missing(keep)) keep <- character(0L)
# Reduce random and fixed parts:
red_random <- eval.parent(reduce_random(object, alpha=alpha.random))
model <- attr(red_random, "model")
# 'model' may be 'lmerMod' rather than 'lmerModLmerTest', so we coerce to
# 'lmerModLmerTest' if required:
if(inherits(model, "lmerMod") && !inherits(model, "lmerModLmerTest"))
model <- as_lmerModLmerTest(model)
stopifnot(inherits(model, "lmerModLmerTest") || inherits(model, "lm"))
red_fixed <- eval.parent(reduce_fixed(model, ddf=ddf,
alpha=alpha.fixed, keep=keep))
# get 'reduction' tables:
step_random <- ran_redTable(red_random)
step_fixed <- fix_redTable(red_fixed)
# organize results and return:
step_list <- list(random=step_random, fixed=step_fixed)
class(step_list) <- "step_list"
attr(step_list, "model") <- attr(red_fixed, "model")
attr(step_list, "drop1") <- attr(red_fixed, "drop1")
step_list
}
##############################################
######## get_model()
##############################################
#' Extract Model from an Object
#'
#' @param x an object.
#' @param ... currently not used.
#'
#' @seealso \code{\link{get_model.step_list}}
#' @export
#' @keywords internal
get_model <- function(x, ...) UseMethod("get_model")
##############################################
######## get_model.step_list()
##############################################
#' @rdname step.lmerModLmerTest
#' @param x a step object.
#' @export
get_model.step_list <- function(x, ...) {
attr(x, "model")
}
##############################################
######## print.step_list()
##############################################
#' @importFrom stats formula
#' @export
#' @keywords internal
print.step_list <- function(x, digits = max(getOption("digits") - 2L, 3L),
signif.stars = getOption("show.signif.stars"),
...) {
print(x[["random"]])
cat("\n")
print(x[["fixed"]])
cat("\nModel found:", deparse2(formula(attr(x, "model"))), sep="\n")
invisible(x)
}
##############################################
######## plot.step_list()
##############################################
#' Plot LS-means for Backward Reduced Model
#'
#' Computes the LS-means for the final backward reduced model and passes these
#' to \code{\link{plot.ls_means}}.
#'
#' Error bars are confidence intervals - the default is 95% CI but the confidence
#' level can be changed.
#'
#' @param x a \code{step_list} object; the result of running
#' \code{\link[=step.lmerModLmerTest]{step}}.
#' @param y not used and ignored with a warning.
#' @param which optional character vector naming factors for which LS-means should
#' be plotted. If \code{NULL} (default) plots for all LS-means are generated.
#' @param mult if \code{TRUE} and there is more than one term for which to plot
#' LS-means the plots are organized in panels with \code{facet_wrap}.
#' @param pairwise pairwise differences of LS-means?
#' @param level confidence level.
#' @param ddf denominator degree of freedom method.
#' @param ... currently not used.
#'
#' @export
#' @author Rune Haubo B. Christensen and Alexandra Kuznetsova
#' @seealso \code{\link[=ls_means.lmerModLmerTest]{ls_means}} and
#' \code{\link{plot.ls_means}}
#' @keywords internal
#' @examples
#'
#' \dontrun{
#' # Fit example model:
#' tv <- lmer(Sharpnessofmovement ~ TVset * Picture +
#' (1 | Assessor:TVset) + (1 | Assessor:Picture) +
#' (1 | Assessor:Picture:TVset) + (1 | Repeat) + (1 | Repeat:Picture) +
#' (1 | Repeat:TVset) + (1 | Repeat:TVset:Picture) + (1 | Assessor),
#' data = TVbo)
#'
#' # Backward reduce the model:
#' (st <- step(tv)) # takes ~10 sec to run
#'
#' # Pairwise comparisons of LS-means for Picture and TVset:
#' plot(st, which=c("Picture", "TVset"), pairwise = TRUE)
#' }
#'
plot.step_list <- function(x, y=NULL, which=NULL, pairwise=FALSE, mult=TRUE,
level=0.95, ddf=c("Satterthwaite", "Kenward-Roger"),
...) {
plot(ls_means(get_model(x), pairwise=pairwise, level=level, ddf=ddf),
y=y, which=which, mult=mult)
}
##############################################
######## step utility functions below
##############################################
ran_redTable <- function(table) {
aov <- attr(table, "ranova")[-1, , drop=FALSE]
stopifnot(nrow(table) >= 1)
tab <- rbind(cbind("Eliminated"=c(NA_real_, seq_len(nrow(table)-1)), table),
cbind("Eliminated"=rep(0, nrow(aov)), aov))
class(tab) <- c("anova", "data.frame")
attr(tab, "heading") <- "Backward reduced random-effect table:\n"
tab
}
fix_redTable <- function(table) {
aov <- attr(table, "drop1")
tab <- rbind(cbind("Eliminated"=seq_len(nrow(table)), table),
cbind("Eliminated"=rep(0, nrow(aov)), aov))
class(tab) <- c("anova", "data.frame")
attr(tab, "heading") <- "Backward reduced fixed-effect table:"
if(!is.null(ddf <- attr(table, "ddf"))) {
ddf <- switch(ddf, "Satterthwaite" = "Satterthwaite",
"Kenward-Roger" = "Kenward-Roger")
attr(tab, "heading") <-
c(attr(tab, "heading"), paste("Degrees of freedom method:", ddf, "\n"))
}
tab
}
#' @importFrom stats formula update
#' @importFrom lme4 getME
reduce_random <- function(model, alpha=0.1) {
ran <- ranova(model)
reduced <- ran[1L, ]
newfit <- model
newform <- formula(model)
forms <- attr(ran, "formulae")
pvals <- ran[-1, "Pr(>Chisq)"]
above <- (!is.na(pvals) & pvals > alpha)
while(any(above)) {
remove <- which.max(pvals)
newform <- forms[[remove]]
reduced <- rbind(reduced, ran[1 + remove, ])
if(!has_ranef(newform)) { # If no random effects: fit with lm
reml <- getME(newfit, "is_REML")
lm_call <- get_lm_call(newfit, formula=newform)
newfit <- eval.parent(as.call(lm_call))
ran <- ranova_lm(newfit, REML=reml)
break
}
newfit <- eval.parent(update(newfit, formula. = newform))
# newfit <- update(newfit, formula = newform)
ran <- ranova(newfit)
forms <- attr(ran, "formulae")
pvals <- ran[-1, "Pr(>Chisq)"]
above <- (!is.na(pvals) & pvals > alpha)
}
attr(reduced, "model") <- newfit
attr(reduced, "formula") <- newform
attr(reduced, "ranova") <- ran
reduced
}
ranova_lm <- function(model, REML=TRUE) {
# Compute a ranova table for an lm-object only containing a '<none>' row
# and the right header.
aov <- mk_LRtab(get_logLik(model, REML=REML))
rownames(aov) <- "<none>"
head <- c("ANOVA-like table for random-effects: Single term deletions",
"\nModel:", deparse2(formula(model)))
# attr(aov, "formulae") <- new_forms
structure(aov, heading = head, class = c("anova", "data.frame"))
}
#' @importFrom stats nobs formula
reduce_fixed <- function(model, ddf=c("Satterthwaite", "Kenward-Roger"), alpha=0.05,
keep) {
if(missing(keep)) keep <- character(0L)
stopifnot(is.character(keep))
term_names <- attr(terms(model), "term.labels")
# Test validity of
if(!all(keep %in% term_names)) {
offending <- paste(setdiff(keep, term_names), collapse = " ")
txt1 <- sprintf("Invalid 'keep' ignored: %s.", offending)
txt2 <- sprintf("Valid terms are: %s.", paste(term_names, collapse = " "))
warning(paste(txt1, txt2, sep="\n"), call. = FALSE)
}
ddf <- match.arg(ddf)
aov <- if(inherits(model, "lmerMod")) drop1.lmerModLmerTest(model, ddf=ddf) else
drop1(model, test="F")[-1L, , drop=FALSE]
reduced <- aov[0L, ]
newfit <- model
newform <- orig_form <- formula(model)
nobs_model <- nobs(model)
terms <- rownames(aov)
consider <- setdiff(terms, keep)
pvals <- aov[consider, "Pr(>F)"]
above <- (!is.na(pvals) & pvals > alpha)
if(any(above)) while(any(above)) {
remove <- consider[which.max(pvals)]
newform <- rm_complete_terms(remove, orig_form, random = FALSE)[[1L]]
reduced <- rbind(reduced, aov[remove, ])
newfit <- eval.parent(update(newfit, formula = newform))
# newfit <- update(newfit, formula = newform)
nobs_newfit <- nobs(newfit)
if(all(is.finite(c(nobs_model, nobs_newfit))) && nobs_newfit != nobs_model)
stop("number of rows in use has changed: remove missing values?",
call.=FALSE)
aov <- if(inherits(newfit, "lmerMod")) drop1.lmerModLmerTest(newfit, ddf=ddf) else
drop1(newfit, test="F")[-1L, , drop=FALSE]
# aov <- drop1(newfit)
orig_form <- formula(newfit)
terms <- rownames(aov)
consider <- setdiff(terms, keep)
pvals <- aov[consider, "Pr(>F)"]
above <- (!is.na(pvals) & pvals > alpha)
}
attr(reduced, "model") <- newfit
attr(reduced, "formula") <- newform
attr(reduced, "drop1") <- aov
attr(reduced, "ddf") <- if(inherits(model, "lmerMod")) ddf else NULL
reduced
}
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.