rm numDeriv dependency (closes #41)

Author: leeper

DESCRIPTION

@@ -4,17 +4,16 @@ 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.2.9
-Date: 2016-08-30
+Version: 0.2.10
+Date: 2016-08-31
 Authors@R: c(person("Thomas J.", "Leeper",
                     role = c("aut", "cre"),
                     email = "thosjleeper@gmail.com"))
 Imports:
     stats,
     graphics,
     compiler,
-    MASS,
-    numDeriv
+    MASS
 Suggests:
     sandwich,
     webuse,

NAMESPACE

@@ -51,4 +51,4 @@ importFrom(graphics,points)
 importFrom(graphics,polygon)
 importFrom(graphics,rug)
 importFrom(graphics,segments)
-importFrom(numDeriv,grad)
+importFrom(utils,head)

NEWS.md

@@ -1,5 +1,9 @@
 # CHANGES TO margins 0.3.0 #
 
+## margins 0.2.10
+
+* Replaced `numDeriv::jacobian()` with an internal alternative. (#41)
+
 ## margins 0.2.8
 
 * Added a `prediction()` method for "ivreg" objects (from `AER::ivreg()`). (#3)

R/build_margins.R

@@ -7,8 +7,7 @@
 #' @param vce A character string indicating the type of estimation procedure to use for estimating variances. The default (\dQuote{delta}) uses the delta method. Alternatives are \dQuote{bootstrap}, which uses bootstrap estimation, or \dQuote{simulation}, which averages across simulations drawn from the joint sampling distribution of model coefficients. The latter two are extremely time intensive.
 #' @param iterations If \code{vce = "bootstrap"}, the number of bootstrap iterations. If \code{vce = "simulation"}, the number of simulated effects to draw. Ignored otherwise.
 #' @param unit_ses If \code{vce = "delta"}, a logical specifying whether to calculate and return unit-specific marginal effect variances. This calculation is time consuming and the information is often not needed, so this is set to \code{FALSE} by default.
-#' @param method A character string indicating the numeric derivative method to use when estimating marginal effects variances. \dQuote{simple} optimizes for speed; \dQuote{Richardson} optimizes for accuracy. See \code{\link[numDeriv]{jacobian}} for details. Currently this is only used when calculating marginal effect variances.
-#' @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.
+#' @param eps A numeric value specifying the \dQuote{step} to use when calculating numerical derivatives.
 #' @param \dots Ignored.
 #' @details Generally, it is not necessary to call this function directly because \code{\link{margins}} provides a simpler interface. To just get marginal effects without building a \dQuote{margins} object, call \code{\link{marginal_effects}} instead, which handles the effect estimation of a model object without building a \dQuote{margins} object.
 #' 
@@ -21,7 +20,6 @@
 #' @keywords models
 #' @import stats
 #' @importFrom compiler cmpfun
-#' @importFrom numDeriv grad
 #' @importFrom MASS mvrnorm
 #' @export
 build_margins <- 
@@ -32,16 +30,14 @@ function(model,
          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()
-         eps = 1e-7,
+         eps = 1e-4,
          ...) {
 
     # variables in the model
     allvars <- all.vars(model[["terms"]])[-1]
 
     # march.arg() for arguments
     type <- match.arg(type)
-    method <- match.arg(method)
     vce <- match.arg(vce)
     if (is.function(vcov)) {
         vcov <- vcov(model)
@@ -53,12 +49,12 @@ function(model,
     # variance estimation technique
     variances <- get_effect_variances(data = data, model = model, allvars = names(mes), 
                                       type = type, vcov = vcov, vce = vce, 
-                                      iterations = iterations, method = method)
+                                      iterations = iterations, eps = eps)
 
     # get unit-specific effect variances (take derivative of `.build_grad_fun()` for every row separately)
     if ((vce == "delta") && (isTRUE(unit_ses))) {
         vmat <- do.call("rbind", lapply(seq_len(nrow(data)), function(datarow) {
-            delta_once(data = data[datarow,], model = model, type = type, vcov = vcov, method = method)
+            delta_once(data = data[datarow,], model = model, type = type, vcov = vcov, eps = eps)
         }))
         colnames(vmat) <- paste0("se.", names(mes))
         vmat <- as.data.frame(vmat)

R/delta.R

@@ -3,29 +3,47 @@ function(data,
          model, 
          type = c("response", "link", "terms"), 
          vcov = vcov(model),
-         method = c("simple", "Richardson", "complex")) {
+         eps = 1e-4) {
     # take the derivative of each marginal effect from a model with respect to each model coefficient
 
     type <- match.arg(type)
-    method <- match.arg(method)
     if (is.function(vcov)) {
         vcov <- vcov(model)
     }
 
-    # express each marginal effect as a function of all coefficients
-    # holding data constant
-    # this is what .build_grad_fun() will do
-    # then:  numDeriv::grad(.build_grad_fun(), model$coef)
-    # gives `gradmat`, such that v %*% V %*% t(v)
+    # express each marginal effect as a function of estimated coefficients
+    # holding data constant (using `.build_grad_fun()`)
+    # use `jacobian(.build_grad_fun(), model$coef)`
+    # to get `jacobian`, an ME-by-beta matrix,
+    # such that jacobian %*% V %*% t(jacobian)
     # gives the variance of each marginal effect
-    # `gradmat` should be an ME-by-beta matrix
-    
-    # Some references:
     # http://www.soderbom.net/lecture10notes.pdf
     # http://stats.stackexchange.com/questions/122066/how-to-use-delta-method-for-standard-errors-of-marginal-effects
 
-    FUN <- .build_grad_fun(data = data, model = model, type = type, method = method)
-    gradmat <- numDeriv::jacobian(FUN, model[["coefficients"]], method = method)
-    vout <- diag(gradmat %*% vcov %*% t(gradmat))
+    FUN <- .build_grad_fun(data = data, model = model, type = type)
+    #jacobian <- numDeriv::jacobian(FUN, model[["coefficients"]], method = "simple")
+    jacobian <- jacobian(FUN, model[["coefficients"]], eps = eps)
+    vout <- diag(jacobian %*% vcov %*% t(jacobian))
     return(vout)
 }
+
+.build_grad_fun <- function(data, model, type = "response", eps = 1e-4) {
+    
+    # factory function to return prediction holding data constant but varying coefficients
+    FUN <- function(coefs) {
+        model[["coefficients"]] <- coefs
+        colMeans(marginal_effects(model = model, data = data, type = type, eps = eps), na.rm = TRUE)
+    }
+    return(compiler::cmpfun(FUN))
+}
+
+jacobian <- function(FUN, coefficients, eps = 1e-4) {
+    F0 <- FUN(coefficients)
+    out <- matrix(NA_real_, nrow = length(F0), ncol = length(coefficients))
+    for (i in seq_along(coefficients)) {
+        coeftemp <- coefficients
+        coeftemp[i] <- coeftemp[i] + eps
+        out[, i] <- (FUN(coeftemp) - F0) / eps
+    }
+    out
+}

R/factories.R

@@ -6,12 +6,11 @@ function(data = data,
          vcov = vcov(model),
          vce = c("delta", "simulation", "bootstrap", "none"),
          iterations = 50L, # if vce == "bootstrap" or "simulation"
-         method = c("simple", "Richardson", "complex"), # passed to marginal_effects()
+         eps = 1e-4,
          ...) {
 
     # march.arg() for arguments
     type <- match.arg(type)
-    method <- match.arg(method)
     vce <- match.arg(vce)
     if (is.function(vcov)) {
         vcov <- vcov(model)
@@ -24,7 +23,7 @@ function(data = data,
     } else if (vce == "delta") {
 
         # default method
-        variances <- delta_once(data = data, model = model, type = type, vcov = vcov, method = method)
+        variances <- delta_once(data = data, model = model, type = type, vcov = vcov, eps = eps)
 
     } else if (vce == "simulation") {
 

R/get_effect_variances.R

@@ -1,4 +1,4 @@
-get_instant_pdiff <- function(data, model, variable, type = c("response", "link"), eps = 1e-7) {
+get_instant_pdiff <- function(data, model, variable, type = c("response", "link"), eps = 1e-4) {
     # @title Instantaneous change in fitted values (numerical derivative)
     # @description This is an internal function used to calculate instantaneous change (numerical derivative) in y-hat between observed values in `data` and the smallest machine-precise change in the value of `data`. This is used by \code{marginal_effects} for numeric variables. It currently only uses the "simple" derivative method. This might change in the future
     # @param data The dataset on which to to calculate `predict(model)` (and the slope thereof)

R/marginal_effects_internal.R

@@ -1,6 +1,6 @@
 #' @rdname marginal_effects
 #' @export
-marginal_effects.lm <- function(model, data, type = c("response", "link"), eps = 1e-7, ...) {
+marginal_effects.lm <- function(model, data, type = c("response", "link"), eps = 1e-4, ...) {
 
     # setup data, if missing
     if (missing(data)) {

R/marginal_effects_methods.R

@@ -31,6 +31,7 @@ prediction <- function(model, data, ...) {
     UseMethod("prediction")
 }
 
+#' @importFrom utils head
 #' @export
 print.prediction <- function(x, digits = 4, ...) {
     f <- x[["fitted"]]
@@ -39,7 +40,7 @@ print.prediction <- function(x, digits = 4, ...) {
         m <- sprintf(paste0("%0.", digits, "f"), m)
         message(paste0("Average prediction: ", m, ", for ", length(f), " ", ngettext(length(f), "observation", "observations")))
     } else if (is.factor(f)) {
-        m <- sort(table(p$fitted), decreasing = TRUE)[1]
+        m <- sort(table(x[["fitted"]]), decreasing = TRUE)[1]
         message(paste0("Modal prediction: ", shQuote(names(m)), " for ", m, " of ", length(f), " ", 
                 ngettext(length(f), "observation", "observations"),
                 " with total ", nlevels(f), " ", ngettext(nlevels(f), "level", "levels") ))

R/prediction.R

@@ -218,5 +218,5 @@ prediction.polr <- function(model, data, ...) {
               class = c("prediction", "data.frame"), 
               row.names = seq_len(length(pred[["fit"]])),
               model.class = class(model),
-              type = type)
+              type = NULL)
 }

R/prediction_methods.R

@@ -29,7 +29,7 @@ By comparison, R has no robust functionality in the base tools for drawing out m
 
 Given the challenges of interpreting the contribution of a given regressor in any model that includes quadratic terms, multiplicative interactions, a non-linear transformation, or other complexities, there is a clear need for a simple, consistent way to estimate marginal effects for popular statistical models. This package aims to correctly calculate marginal effects that include complex terms and provide a uniform interface for doing those calculations. Thus, the package implements a single S3 generic method (`margins()`) that can be easily generalized for any type of model implemented in R.
 
-Some technical details of the package are worth briefly noting. The estimation of marginal effects relies on numeric derivatives produced using `predict()` and [`numDeriv::grad()`](https://cran.r-project.org/package=numDeriv). While symbolic differentiation of some models (e.g., linear models) is possible using `D()` and `deriv()`, R's modelling language (the "formula" class) is sufficiently general to enable the construction of model formulae that contain terms that fall outside of R's symbolic differentiation rule table (e.g., `y ~ factor(x)` or `y ~ I(FUN(x))` for any arbitrary `FUN()`). By relying on numeric differentiation, `margins()` supports *any* model that can be expressed in R formula syntax. Even Stata's `margins` command is limited in its ability to handle variable transformations (e.g., including `x` and `log(x)` as predictors) and quadratic terms (e.g., `x^3`); these scenarios are easily expressed in an R formula and easily handled, correctly, by `margins()`.
+Some technical details of the package are worth briefly noting. The estimation of marginal effects relies on numeric derivatives produced using `predict()` and a numerical approximation of [the Jacobian matrix](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant). While symbolic differentiation of some models (e.g., linear models) is possible using `D()` and `deriv()`, R's modelling language (the "formula" class) is sufficiently general to enable the construction of model formulae that contain terms that fall outside of R's symbolic differentiation rule table (e.g., `y ~ factor(x)` or `y ~ I(FUN(x))` for any arbitrary `FUN()`). By relying on numeric differentiation, `margins()` supports *any* model that can be expressed in R formula syntax. Even Stata's `margins` command is limited in its ability to handle variable transformations (e.g., including `x` and `log(x)` as predictors) and quadratic terms (e.g., `x^3`); these scenarios are easily expressed in an R formula and easily handled, correctly, by `margins()`.
 
 ## Simple code examples ##
 

README.Rmd

@@ -29,7 +29,7 @@ By comparison, R has no robust functionality in the base tools for drawing out m
 
 Given the challenges of interpreting the contribution of a given regressor in any model that includes quadratic terms, multiplicative interactions, a non-linear transformation, or other complexities, there is a clear need for a simple, consistent way to estimate marginal effects for popular statistical models. This package aims to correctly calculate marginal effects that include complex terms and provide a uniform interface for doing those calculations. Thus, the package implements a single S3 generic method (`margins()`) that can be easily generalized for any type of model implemented in R.
 
-Some technical details of the package are worth briefly noting. The estimation of marginal effects relies on numeric derivatives produced using `predict()` and [`numDeriv::grad()`](https://cran.r-project.org/package=numDeriv). While symbolic differentiation of some models (e.g., linear models) is possible using `D()` and `deriv()`, R's modelling language (the "formula" class) is sufficiently general to enable the construction of model formulae that contain terms that fall outside of R's symbolic differentiation rule table (e.g., `y ~ factor(x)` or `y ~ I(FUN(x))` for any arbitrary `FUN()`). By relying on numeric differentiation, `margins()` supports *any* model that can be expressed in R formula syntax. Even Stata's `margins` command is limited in its ability to handle variable transformations (e.g., including `x` and `log(x)` as predictors) and quadratic terms (e.g., `x^3`); these scenarios are easily expressed in an R formula and easily handled, correctly, by `margins()`.
+Some technical details of the package are worth briefly noting. The estimation of marginal effects relies on numeric derivatives produced using `predict()` and a numerical approximation of [the Jacobian matrix](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant). While symbolic differentiation of some models (e.g., linear models) is possible using `D()` and `deriv()`, R's modelling language (the "formula" class) is sufficiently general to enable the construction of model formulae that contain terms that fall outside of R's symbolic differentiation rule table (e.g., `y ~ factor(x)` or `y ~ I(FUN(x))` for any arbitrary `FUN()`). By relying on numeric differentiation, `margins()` supports *any* model that can be expressed in R formula syntax. Even Stata's `margins` command is limited in its ability to handle variable transformations (e.g., including `x` and `log(x)` as predictors) and quadratic terms (e.g., `x^3`); these scenarios are easily expressed in an R formula and easily handled, correctly, by `margins()`.
 
 ## Simple code examples ##
 
@@ -85,7 +85,7 @@ microbenchmark(marginal_effects(x))
 ```
 ## Unit: milliseconds
 ##                 expr      min       lq     mean   median       uq      max neval
-##  marginal_effects(x) 7.942455 8.409444 9.199376 9.115538 9.792552 13.37642   100
+##  marginal_effects(x) 8.047778 8.400351 9.020149 8.620828 9.537136 12.67893   100
 ```
 
 ```r
@@ -94,8 +94,8 @@ microbenchmark(margins(x))
 
 ```
 ## Unit: milliseconds
-##        expr      min       lq     mean   median      uq      max neval
-##  margins(x) 63.40545 69.59018 74.86074 72.59378 75.9127 169.0655   100
+##        expr      min       lq     mean   median       uq      max neval
+##  margins(x) 63.93785 67.05056 74.14901 70.09276 76.57553 180.7601   100
 ```
 
 In addition to the estimation procedures and `plot()` generic, **margins** offers several plotting methods for model objects. First, there is a new generic `cplot()` that displays predictions or marginal effects (from an "lm" or "glm" model) of a variable conditional across values of third variable (or itself). For example, here is a graph of predicted probabilities from a logit model: