expand and test support for lme4 'merMod' models (#56)

leeper · May 11, 2018 · 3b01749fc1f78ba94f3489303cf0a9ee09ba1da4 · 3b01749
1 parent d3eeca4
commit 3b01749fc1f78ba94f3489303cf0a9ee09ba1da4
diff --git a/.travis.yml b/.travis.yml
@@ -21,9 +21,11 @@ r_packages:
 - rmarkdown
 - testthat
 - covr
+- AER
 - betareg
 - ggplot2
 - gapminder
+- lme4
 - nnet
 - ordinal
 - sandwich

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -4,7 +4,7 @@ Title: Marginal Effects for Model Objects
 Description: An R port of Stata's 'margins' command, which can be used to
     calculate marginal (or partial) effects from model objects.
 License: MIT + file LICENSE
-Version: 0.3.21
+Version: 0.3.22
 Date: 2018-05-11
 Authors@R: c(person("Thomas J.", "Leeper",
                     role = c("aut", "cre"), 
@@ -27,17 +27,18 @@ Imports:
     grDevices,
     MASS
 Suggests:
+    methods,
     knitr,
     rmarkdown,
     testthat,
     ggplot2,
     gapminder,
     sandwich,
-    stargazer
+    stargazer,
+    lme4
 Enhances:
     AER,
     betareg,
-    lme4,
     nnet,
     ordinal,
     survey

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+# margins 0.3.22
+
+* Expanded support for objects of class "merMod" from **lme4**, including support for variance estimation and an expanded test suite. (#56)
+
 # margins 0.3.21
 
 * Modified the internals of `gradient_factory()` to be more robust to an expanded set of model classes through the introduction of an internal function `reset_coefs()`. A test suite for this function has been added.

diff --git a/R/find_terms_in_model.R b/R/find_terms_in_model.R
@@ -12,14 +12,22 @@ find_terms_in_model.default <- function(model, variables = NULL) {
     # identify classes of terms in `model`
     if (!is.null(attributes(terms(model))[["dataClasses"]])) {
         ## first look in the `terms(model)`
-        classes <- attributes(terms(model))[["dataClasses"]][-1]
-    } else if (!is.null(attributes(terms(model$model))[["dataClasses"]])) {
+        classes <- attributes(terms(model))[["dataClasses"]][-1L]
+    } else if ("model" %in% names(model) && !is.null(attributes(terms(model$model))[["dataClasses"]])) {
         ## then look in the `terms(model$model)`
-        classes <- attributes(terms(model$model))[["dataClasses"]][-1]
+        classes <- attributes(terms(model$model))[["dataClasses"]][-1L]
     } else {
-        ## probably should try something else before giving up but we'll leave it like this for now
-        ## ^ famous last words
-        stop("No variable classes found in model.")
+        ## then get desperate and use `find_data()`
+        ## this is used for "merMod" objects
+        att <- attributes(terms(find_data(model)))
+        if ("dataClasses" %in% names(att)) {
+            classes <- att[["dataClasses"]][-1L]
+            rm(att)
+        } else {
+            ## probably should try something else before giving up but we'll leave it like this for now
+            ## ^ famous last words
+            stop("No variable classes found in model.")
+        }
     }
 
     # drop specially named "(weights)" variables
@@ -65,6 +73,22 @@ find_terms_in_model.default <- function(model, variables = NULL) {
     return(vars)
 }
 
+find_terms_in_model.merMod <- function(model, variables = NULL) {
+
+    # require lme4 package in order to identify random effects terms
+    requireNamespace("lme4")
+
+    # call default method
+    varslist <- find_terms_in_model.default(model, variables)
+
+    # sanitize `varlist` to remove random effects terms
+    fixed <- as.character(formula(model, fixed.only = TRUE))[-c(1:2)]
+    varslist[] <- lapply(varslist, function(vartype) vartype[vartype %in% fixed])
+    varslist
+}
+
+find_terms_in_model.lmerMod <- find_terms_in_model.merMod
+
 # call gsub_bracket on all common formula operations
 clean_terms <- function(terms) {
     # the use of paste("`", x, "`") is a hack to deal with variables that have spaces in their names

diff --git a/R/get_effect_variances.R b/R/get_effect_variances.R
@@ -54,9 +54,16 @@ function(data,
                                 varslist = varslist,
                                 ...)
         # get jacobian
-        jacobian <- jacobian(FUN, coef(model)[names(coef(model)) %in% c("(Intercept)", colnames(vcov))], weights = weights, eps = eps)
-        # sandwich
-        vc <- jacobian %*% vcov %*% t(jacobian)
+        if (inherits(model, "merMod")) {
+            requireNamespace("lme4")
+            jacobian <- jacobian(FUN, lme4::fixef(model)[names(lme4::fixef(model)) %in% c("(Intercept)", colnames(vcov))], weights = weights, eps = eps)
+            # sandwich
+            vc <- as.matrix(jacobian %*% vcov %*% t(jacobian))
+        } else {
+            jacobian <- jacobian(FUN, coef(model)[names(coef(model)) %in% c("(Intercept)", colnames(vcov))], weights = weights, eps = eps)
+            # sandwich
+            vc <- jacobian %*% vcov %*% t(jacobian)
+        }
         # extract variances from diagonal
         variances <- diag(vc)
 

diff --git a/R/marginal_effects.R b/R/marginal_effects.R
@@ -1,5 +1,5 @@
 #' @rdname marginal_effects
-#' @title Differentiate a Model Object with Respect to All Variables
+#' @title Differentiate a Model Object with Respect to All (or Specified) Variables
 #' @description Extract marginal effects from a model object, conditional on data, using \code{\link{dydx}}.
 #' @param data A data.frame over which to calculate marginal effects. This is optional, but may be required when the underlying modelling function sets \code{model = FALSE}.
 #' @param variables A character vector with the names of variables for which to compute the marginal effects. The default (\code{NULL}) returns marginal effects for all variables.
@@ -9,7 +9,7 @@
 #' @param as.data.frame A logical indicating whether to return a data frame (the default) or a matrix.
 #' @param varslist A list structure used internally by \code{\link{margins}}. Users should not set this.
 #' @param \dots Arguments passed to methods, and onward to \code{\link{dydx}} methods and possibly further to \code{\link[prediction]{prediction}} methods. This can be useful, for example, for setting \code{type} (predicted value type), \code{eps} (precision), or \code{category} (category for multi-category outcome models), etc.
-#' @details This function extracts unit-specific marginal effects from an estimated model with respect to \emph{all} variables specified in \code{data} and returns a data.frame. (Note that this is not each \emph{coefficient}.) See \code{\link{dydx}} for computational details, or to extract the marginal effect for only one variable. Note that for factor and logical class variables, discrete changes in the outcome are reported rather than instantaneous marginal effects.
+#' @details Users likely want to use the fully featured \code{\link{margins}} function rather than \code{marginal_effects}, which merely performs estimation of the marginal effects but simply returns a data frame. \code{\link{margins}}, by contrast, does some convenient packaging around these results and supports additional functionality, like variance estimation and counterfactual estimation procedures. The methods for this function provide lower-level functionality that extracts unit-specific marginal effects from an estimated model with respect to \emph{all} variables specified in \code{data} (or the subset specified in \code{variables}) and returns a data frame. See \code{\link{dydx}} for computational details. Note that for factor and logical class variables, discrete changes in the outcome are reported rather than instantaneous marginal effects.
 #' 
 #' Methods are currently implemented for the following object classes:
 #' \itemize{
@@ -18,13 +18,13 @@
 #'   \item \dQuote{ivreg}, see \code{\link[AER]{ivreg}}
 #'   \item \dQuote{lm}, see \code{\link[stats]{lm}}
 #'   \item \dQuote{loess}, see \code{\link[stats]{loess}}
-#'   \item \dQuote{merMod}, see \code{\link[lme4]{lmer}}
+#'   \item \dQuote{merMod}, see \code{\link[lme4]{lmer}}, \code{\link[lme4]{glmer}}
 #'   \item \dQuote{nnet}, see \code{\link[nnet]{nnet}}
 #'   \item \dQuote{polr}, see \code{\link[MASS]{polr}}
 #'   \item \dQuote{svyglm}, see \code{\link[survey]{svyglm}}
 #' }
 #'
-#' A methods is also provided for the object classes \dQuote{margins} to return a simplified data frame from complete \dQuote{margins} objects.
+#' A method is also provided for the object classes \dQuote{margins} to return a simplified data frame from complete \dQuote{margins} objects.
 #' 
 #' @return An data frame with number of rows equal to \code{nrow(data)}, where each row is an observation and each column is the marginal effect of a variable used in the model formula.
 #' @examples

diff --git a/R/marginal_effects_default.R b/R/marginal_effects_default.R
@@ -30,7 +30,6 @@ function(model,
     }
 
     out <- c(out1, out2)
-
     if (isTRUE(as.data.frame)) {
         out <- do.call("cbind.data.frame", out[vapply(out, function(x) length(x) > 0, FUN.VALUE = logical(1))])
     } else {

diff --git a/R/marginal_effects_merMod.R b/R/marginal_effects_merMod.R
@@ -1,44 +1,32 @@
 #' @rdname marginal_effects
 #' @importFrom prediction find_data
 #' @export
-marginal_effects.merMod <- 
-function(model, 
-         data = find_data(model), 
+marginal_effects.merMod <-
+function(model,
+         data = find_data(model),
          variables = NULL,
-         type = c("response", "link"), 
+         type = c("response", "link"),
          eps = 1e-7,
          as.data.frame = TRUE,
+         varslist = NULL,
          ...) {
 
     type <- match.arg(type)
 
-    # extract term names from `model`
-    term_names <- all.vars(terms(model))
-
     # identify classes of terms in `model`
-    nnames <- colnames(attributes(terms(model))[["factors"]])
-    lnames <- NULL
-    fnames <- NULL
-    classes <- rep("numeric", length(nnames))
-    warning("factor variables are not handled as factor for models of class 'merMod'")
-
-    # subset of variables for which to compute the marginal effects
-    if (!is.null(variables)) {
-        if (any(!variables %in% nnames)) {
-            stop("Some values in 'variables' are not in the model terms.")
-        }
-        nnames <- nnames[nnames %in% variables]
+    if (is.null(varslist)) {
+        varslist <- find_terms_in_model(model, variables = variables)
     }
 
     # estimate numerical derivatives with respect to each variable (for numeric terms in the model)
     # add discrete differences for logical terms
-    out1 <- lapply(c(nnames, lnames), dydx, data = data, model = model, type = type, eps = eps, as.data.frame = as.data.frame, ...)
+    out1 <- lapply(c(varslist$nnames, varslist$lnames), dydx, data = data, model = model, type = type, eps = eps, as.data.frame = as.data.frame, ...)
 
     # add discrete differences for factor terms
     ## exact number depends on number of factor levels
     out2 <- list()
-    for (i in seq_along(fnames)) {
-        out2[[i]] <- dydx.factor(data = data, model = model, fnames[i], type = type, fwrap = FALSE, as.data.frame = as.data.frame, ...)
+    for (i in seq_along(varslist$fnames)) {
+        out2[[i]] <- dydx.factor(data = data, model = model, varslist$fnames[i], type = type, fwrap = FALSE, as.data.frame = as.data.frame, ...)
     }
 
     out <- c(out1, out2)