rename vc argument to vcov (closes #34)

NEWS.md

@@ -3,6 +3,7 @@
 ## margins 0.2.2
 
 * Removed support for "marginal effects at means" (MEMs) and the `atmeans` argument throughout package. (#35)
+* Renamed the `vc` argument to `vcov` for consistency with other packages. (#34)
 
 ## margins 0.2.1
 

R/build_margins.R

@@ -3,7 +3,7 @@
 #' @param model A model object.
 #' @param data A data.frame over which to calculate marginal effects.
 #' @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 vc A matrix containing the variance-covariance matrix for estimated model coefficients.
+#' @param vcov A matrix containing the variance-covariance matrix for estimated model coefficients, a function to perform the estimation with \code{model} as its only argument
 #' @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 method A character string indicating the numeric derivative method to use when estimating marginal effects. \dQuote{simple} optimizes for speed; \dQuote{Richardson} optimizes for accuracy. See \code{\link[numDeriv]{grad}} for details.
@@ -26,7 +26,7 @@ build_margins <-
 function(model, 
          data,
          type = c("response", "link", "terms"),
-         vc = vcov(model),
+         vcov = vcov(model),
          vce = c("delta", "simulation", "bootstrap"),
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          method = c("simple", "Richardson", "complex"), # passed to marginal_effects()
@@ -39,18 +39,21 @@ function(model,
     type <- match.arg(type)
     method <- match.arg(method)
     vce <- match.arg(vce)
+    if (is.function(vcov)) {
+        vcov <- vcov(model)
+    }
 
     # obtain gradient with respect to each variable in data
     mes <- marginal_effects(model = model, data = data, type = type, method = method)
 
     # variance estimation technique
     variances <- get_effect_variances(data = data, model = model, allvars = names(mes), 
-                                      type = type, vc = vc, vce = vce, 
+                                      type = type, vcov = vcov, vce = vce, 
                                       iterations = iterations, method = method)
     # get unit-specific effect variances (take derivative of `.build_grad_fun()` for every row separately)
     if (vce == "delta") {
         vmat <- do.call("rbind", lapply(seq_len(nrow(data)), function(datarow) {
-            delta_once(data = data[datarow,], model = model, type = type, vc = vc, method = method)
+            delta_once(data = data[datarow,], model = model, type = type, vcov = vcov, method = method)
         }))
         colnames(vmat) <- paste0("se.", names(mes))
         vmat <- as.data.frame(vmat)

R/delta.R

@@ -2,12 +2,15 @@ delta_once <-
 function(data, 
          model, 
          type = c("response", "link", "terms"), 
-         vc = vcov(model),
+         vcov = vcov(model),
          method = c("simple", "Richardson", "complex")) {
     # 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
@@ -23,6 +26,6 @@ function(data,
 
     FUN <- .build_grad_fun(data = data, model = model, type = type, method = method)
     gradmat <- numDeriv::jacobian(FUN, model[["coefficients"]], method = method)
-    vout <- diag(gradmat %*% vc %*% t(gradmat))
+    vout <- diag(gradmat %*% vcov %*% t(gradmat))
     return(vout)
 }

R/get_effect_variances.R

@@ -3,7 +3,7 @@ function(data = data,
          model = model, 
          which = all.vars(model[["terms"]])[-1], # which mes do we need variances of
          type = c("response", "link", "terms"),
-         vc = vcov(model),
+         vcov = vcov(model),
          vce = c("delta", "simulation", "bootstrap"),
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          method = c("simple", "Richardson", "complex"), # passed to marginal_effects()
@@ -13,11 +13,14 @@ function(data = data,
     type <- match.arg(type)
     method <- match.arg(method)
     vce <- match.arg(vce)
+    if (is.function(vcov)) {
+        vcov <- vcov(model)
+    }
 
     if (vce == "delta") {
 
         # default method
-        variances <- delta_once(data = data, model = model, type = type, vc = vc, method = method)
+        variances <- delta_once(data = data, model = model, type = type, vcov = vcov, method = method)
 
     } else if (vce == "simulation") {
 
@@ -26,7 +29,7 @@ function(data = data,
         tmpmodel$model <- NULL # remove data from model for memory
 
         # simulate from multivariate normal
-        coefmat <- MASS::mvrnorm(iterations, coef(model), vc)
+        coefmat <- MASS::mvrnorm(iterations, coef(model), vcov)
 
         # estimate AME from from each simulated coefficient vector
         effectmat <- apply(coefmat, 1, function(coefrow) {