add loess support (#3)

NAMESPACE

@@ -3,17 +3,22 @@
 S3method(confint,margins)
 S3method(cplot,glm)
 S3method(cplot,lm)
+S3method(cplot,loess)
 S3method(extract_marginal_effects,margins)
 S3method(extract_marginal_effects,marginslist)
 S3method(marginal_effects,glm)
 S3method(marginal_effects,lm)
+S3method(marginal_effects,loess)
 S3method(margins,lm)
+S3method(margins,loess)
 S3method(persp,glm)
 S3method(persp,lm)
+S3method(persp,loess)
 S3method(plot,margins)
 S3method(plot,marginslist)
 S3method(prediction,glm)
 S3method(prediction,lm)
+S3method(prediction,loess)
 S3method(print,margins)
 S3method(print,marginslist)
 S3method(summary,margins)

NEWS.md

@@ -3,6 +3,9 @@
 ## margins 0.2.7
 
 * `marginal_effects()` and `prediction()` are now S3 generics, with methods for "lm" and "glm" objects, improving extensability. (#39, #40)
+* Added preliminary support for "loess" objects, including methods for `prediction()`, `marginal_effects()`, `cplot()`, and `persp()`. No effect variances are currently calculated.
+* `prediction()` returns a new class ("prediction") and gains a `print()` method.
+* Internal function `get_effect_variances()` gains a "none" option for the `vce` argument, to skip calculation of ME variances.
 
 ## margins 0.2.7
 

R/build_margins.R

@@ -29,7 +29,7 @@ function(model,
          data,
          type = c("response", "link", "terms"),
          vcov = stats::vcov(model),
-         vce = c("delta", "simulation", "bootstrap"),
+         vce = c("delta", "simulation", "bootstrap", "none"),
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          method = c("simple", "Richardson", "complex"), # passed to marginal_effects()
@@ -77,7 +77,7 @@ function(model,
                   cbind(data, pred, mes)
               }, 
               class = c("margins", "data.frame"), 
-              Variances = setNames(variances, names(mes)),
+              Variances = if (is.null(variances)) variances else setNames(variances, names(mes)),
               type = type,
               call = model[["call"]],
               df.residual = model[["df.residual"]],

R/confint.R

@@ -15,7 +15,8 @@ function(object, parm, level = 0.95, ...) {
     a <- c(a, 1 - a)
     fac <- qnorm(a)
     ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, a))
-    ses <- sqrt(attributes(object)[["Variances"]])
+    variances <- attributes(object)[["Variances"]]
+    ses <- if (is.null(variances)) rep(NA_real_, length(cf)) else sqrt(variances)
     ci[] <- cf + ses %o% fac
     colnames(ci) <- sprintf("%0.2f%%", 100 * a)
     ci

R/cplot.R

@@ -254,3 +254,7 @@ function(object,
 #' @rdname cplot
 #' @export
 cplot.glm <- cplot.lm
+
+#' @rdname cplot
+#' @export
+cplot.loess <- cplot.lm

R/get_effect_variances.R

@@ -4,7 +4,7 @@ function(data = data,
          which = all.vars(model[["terms"]])[-1], # which mes do we need variances of
          type = c("response", "link", "terms"),
          vcov = vcov(model),
-         vce = c("delta", "simulation", "bootstrap"),
+         vce = c("delta", "simulation", "bootstrap", "none"),
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          method = c("simple", "Richardson", "complex"), # passed to marginal_effects()
          ...) {
@@ -17,7 +17,11 @@ function(data = data,
         vcov <- vcov(model)
     }
 
-    if (vce == "delta") {
+    if (vce == "none") {
+
+        return(NULL)
+
+    } else if (vce == "delta") {
 
         # default method
         variances <- delta_once(data = data, model = model, type = type, vcov = vcov, method = method)

R/marginal_effects.R

@@ -1,7 +1,7 @@
 #' @rdname marginal_effects
 #' @title Differentiate a Model Object
 #' @description Extract marginal effects (via numerical differentiation) and predicted differences in factor changes from a model object, conditional on data
-#' @param data A data.frame over which to calculate marginal effects.
+#' @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 model A model object, perhaps returned by \code{\link[stats]{lm}} or \code{\link[stats]{glm}}
 #' @param type A character string indicating the type of marginal effects to estimate. Mostly relevant for non-linear models, where the reasonable options are \dQuote{response} (the default) or \dQuote{link} (i.e., on the scale of the linear predictor in a GLM).
 #' @param eps A numeric value specifying the \dQuote{step} to use when calculating numerical derivatives. By default this is the smallest floating point value that can be represented on the present architecture.
@@ -82,3 +82,7 @@ marginal_effects.lm <- function(model, data, type = c("response", "link"), eps =
 #' @rdname marginal_effects
 #' @export
 marginal_effects.glm <- marginal_effects.lm
+
+#' @rdname marginal_effects
+#' @export
+marginal_effects.loess <- marginal_effects.lm

R/margins.R

@@ -2,11 +2,13 @@
 #' @title Marginal Effects Estimation
 #' @description \code{margins} is an S3 generic function for building a \dQuote{margins} object from a model object. Methods are currently implemented for \dQuote{lm} (and, implicitly, \dQuote{glm}) class objects.
 #' @param model A model object of class \dQuote{lm}.
-#' @param data A data.frame containing the data at which to evaluate the marginal effects, as in \code{\link[stats]{predict}}.
+#' @param data A data.frame containing the data at which to evaluate the marginal effects, as in \code{\link[stats]{predict}}. This is optional, but may be required when the underlying modelling function sets \code{model = FALSE}.
 #' @param at A list of one or more named vectors, specifically values at which to calculate the marginal effects. See \code{\link{build_datalist}} for details on use.
 #' @param \dots Arguments passed to methods, and in turn to \code{\link{build_margins}}.
 #' @details Methods for this generic return a \dQuote{marginslist} list of one or more \dQuote{margins} objects. A \dQuote{margins} object is a data.frame consisting of the original data, predicted values and standard errors thereof, estimated marginal effects from the model \code{model}, with attributes describing various features of the marginal effects estimates.
 #' 
+#' Some modelling functions set \code{model = FALSE} by default. For \code{margins} to work best, this should be set to \code{TRUE}. Otherwise the \code{data} argument to \code{margins} is probably required.
+#' 
 #' See \code{\link{marginal_effects}} for details on estimation of marginal effects.
 #'
 #' The \code{margins} method for objects of class \dQuote{lm} or \dQuote{glm} simply constructs a list of data.frames (using \code{\link{build_datalist}}) and calls \code{\link{build_margins}} on each. You can call \code{build_margins} directly, but it requires the explicit specification of a dataset over which to estimate the quantities of interest.
@@ -77,3 +79,35 @@ function(model,
     # return value
     structure(out, class = "marginslist")
 }
+
+#' @rdname margins
+#' @export
+margins.loess <- function(model, 
+         data, 
+         at = NULL, 
+         ...){
+
+    # setup data
+    if (missing(data)) {
+        if (!is.null(model[["call"]][["data"]])) {
+            data <- eval(model[["call"]][["data"]], parent.frame()) 
+        } else { 
+            data <- get_all_vars(model[["terms"]], data = model[["model"]])
+        }
+    }
+    data_list <- build_datalist(data, at = at)
+
+    # reduce memory profile
+    model[["model"]] <- NULL
+    attr(model[["terms"]], ".Environment") <- NULL
+
+    # calculate marginal effects
+    out <- lapply(data_list, function(thisdata) {
+        m <- build_margins(model = model, data = thisdata, vcov = NULL, vce = "none", ...)
+        attr(m, "at") <- attributes(thisdata)[["at"]]
+        m
+    })
+
+    # return value
+    structure(out, class = "marginslist")
+}

R/persp.R

@@ -1,3 +1,4 @@
+#' @rdname persp
 #' @title Perspective plots for models
 #' @description Draw one or more perspectives plots reflecting predictions or marginal effects from a model. Currently methods exist for \dQuote{lm} and \dQuote{glm} models.
 #' @param x A model object.
@@ -114,5 +115,10 @@ function(x,
     invisible(out)
 }
 
+#' @rdname persp
 #' @export
 persp.glm <- persp.lm
+
+#' @rdname persp
+#' @export
+persp.loess <- persp.lm

R/prediction.R

@@ -6,7 +6,7 @@
 #' @param type A character string indicating the type of marginal effects to estimate. Mostly relevant for non-linear models, where the reasonable options are \dQuote{response} (the default) or \dQuote{link} (i.e., on the scale of the linear predictor in a GLM).
 #' @param \dots Additional arguments passed to \code{\link[stats]{predict}} methods.
 #' @details This function is simply a wrapper around \code{\link[stats]{predict}} that returns a data.frame containing predicted values with respect to all variables specified in \code{data}. It is used internally by \code{\link{build_margins}}.
-#' @return An data.frame with number of rows equal to number of rows in \code{data}, where each row is an observation and the two columns represent fitted/predicted values and the standard errors thereof.
+#' @return A data.frame with class \dQuote{prediction} that has a number of rows equal to number of rows in \code{data}, where each row is an observation and the two columns represent fitted/predicted values (\code{fitted}) and the standard errors thereof (\code{se.fitted}).
 #' require("datasets")
 #' x <- lm(mpg ~ cyl * hp + wt, data = mtcars)
 #' prediction(x)
@@ -40,7 +40,7 @@ prediction.lm <- function(model, data, type = "response", ...) {
     # obs-x-2 data.frame of predictions
     structure(list(fitted = pred[["fit"]], 
                    se.fitted = pred[["se.fit"]]), 
-              class = "data.frame", 
+              class = c("prediction", "data.frame"), 
               row.names = seq_len(length(pred[["fit"]])))
 }
 
@@ -66,6 +66,31 @@ prediction.glm <- function(model, data, type = c("response", "link"), ...) {
     # obs-x-2 data.frame of predictions
     structure(list(fitted = pred[["fit"]], 
                    se.fitted = pred[["se.fit"]]), 
-              class = "data.frame", 
+              class = c("prediction", "data.frame"), 
+              row.names = seq_len(length(pred[["fit"]])))
+}
+#' @rdname prediction
+#' @export
+prediction.loess <- function(model, data, type = "response", ...) {
+    # setup data
+    if (missing(data)) {
+        if (!is.null(model[["call"]][["data"]])) {
+            data <- eval(model[["call"]][["data"]], parent.frame()) 
+        } else { 
+            data <- get_all_vars(model[["terms"]], data = model[["model"]])
+        }
+    }
+
+    type <- match.arg(type)
+
+    # extract predicted value at input value (value can only be 1 number)
+    pred <- predict(model, newdata = data, type = type, se = TRUE, ...)
+    class(pred[["fit"]]) <- c("fit", "numeric")
+    class(pred[["se.fit"]]) <- c("se.fit", "numeric")
+
+    # obs-x-2 data.frame of predictions
+    structure(list(fitted = pred[["fit"]], 
+                   se.fitted = pred[["se.fit"]]), 
+              class = c("prediction", "data.frame"), 
               row.names = seq_len(length(pred[["fit"]])))
 }

R/print.R

@@ -1,7 +1,7 @@
 #' @export
 print.margins <- 
 function(x, digits = 4, ...) {
-    print(colMeans(extract_marginal_effects(x)), digits = digits, ...)
+    print(colMeans(extract_marginal_effects(x), na.rm = TRUE), digits = digits, ...)
     invisible(x)
 }
 

R/summary.R

@@ -12,9 +12,10 @@ function(object, digits = 4, level = 0.95, ...) {
     }
     fmt <- paste0("%0.", ifelse(digits > 7, 7, digits), "f")
     mes <- extract_marginal_effects(object)
+    variances <- attributes(mes)[["Variances"]]
     tab <- structure(list(Factor = names(mes), 
-                          "dy/dx" = colMeans(mes),
-                          "Std.Err." = sqrt(attributes(mes)[["Variances"]])
+                          "dy/dx" = colMeans(mes, na.rm = TRUE),
+                          "Std.Err." = if (is.null(variances)) rep(NA_real_, ncol(mes)) else sqrt(variances)
                           ),
                      class = "data.frame", row.names = names(mes))
     tab[["z value"]] <- tab[,"dy/dx"]/tab[,"Std.Err."]
@@ -29,5 +30,5 @@ function(object, digits = 4, level = 0.95, ...) {
 #' @export
 summary.marginslist <- 
 function(object, row.names = FALSE, ...) {
-    invisible(lapply(object, summary, ...))
+    lapply(object, summary, ...)
 }