work on clm support (#60)

NAMESPACE

@@ -14,6 +14,7 @@ S3method(dydx,ordered)
 S3method(image,glm)
 S3method(image,lm)
 S3method(image,loess)
+S3method(marginal_effects,clm)
 S3method(marginal_effects,default)
 S3method(marginal_effects,glm)
 S3method(marginal_effects,lm)
@@ -24,6 +25,7 @@ S3method(marginal_effects,merMod)
 S3method(marginal_effects,nnet)
 S3method(marginal_effects,polr)
 S3method(margins,betareg)
+S3method(margins,clm)
 S3method(margins,default)
 S3method(margins,glm)
 S3method(margins,lm)

R/marginal_effects_clm.R

@@ -0,0 +1,41 @@
+#' @rdname marginal_effects
+#' @importFrom prediction find_data
+#' @export
+marginal_effects.clm <- 
+function(model, 
+         data = find_data(model, parent.frame()), 
+         variables = NULL,
+         type = NULL,
+         eps = 1e-7, 
+         varslist = NULL,
+         as.data.frame = TRUE,
+         ...) {
+
+    if (!is.null(type)) {
+        warning(sprintf("'type' is ignored for models of class '%s'", class(model)))
+    }
+
+    # identify classes of terms in `model`
+    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(varslist$nnames, varslist$lnames), dydx, data = data, model = model, type = NULL, 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(varslist$fnames)) {
+        out2[[i]] <- dydx.factor(data = data, model = model, varslist$fnames[i], fwrap = FALSE, type = NULL, as.data.frame = as.data.frame, ...)
+    }
+
+    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 {
+        out <- do.call("cbind", out[vapply(out, function(x) length(x) > 0, FUN.VALUE = logical(1))])
+    }
+    return(out)
+}

R/marginal_effects_glm.R

@@ -7,8 +7,8 @@ function(model,
          variables = NULL,
          type = c("response", "link"), 
          eps = 1e-7, 
-         as.data.frame = TRUE,
          varslist = NULL,
+         as.data.frame = TRUE,
          ...) {
 
     type <- match.arg(type)
@@ -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 {

R/margins_clm.R

@@ -0,0 +1,58 @@
+#' @rdname margins
+#' @export
+margins.clm <- 
+function(model, 
+         data = find_data(model, parent.frame()), 
+         variables = NULL,
+         at = NULL, 
+         type = c("response", "link"),
+         vce = "none",
+         eps = 1e-7,
+         ...) {
+
+    # match.arg()
+    type <- match.arg(type)
+    if (type == "link") {
+        type <- "linear.predictor"
+    }
+    vce <- match.arg(vce)
+
+    # setup data
+    data_list <- build_datalist(data, at = at)
+    if (is.null(names(data_list))) {
+        names(data_list) <- NA_character_
+    }
+    at_specification <- attr(data_list, "at_specification")
+
+    # identify classes of terms in `model`
+    varslist <- find_terms_in_model(model, variables = variables)
+
+    # calculate marginal effects
+    out <- list()
+    for (i in seq_along(data_list)) {
+        out[[i]] <- build_margins(model = model, 
+                                  data = data_list[[i]], 
+                                  variables = variables,
+                                  type = type,
+                                  vce = vce, 
+                                  vcov = NULL,
+                                  unit_ses = FALSE, 
+                                  eps = eps,
+                                  varslist = varslist,
+                                  ...)
+        out[[i]][["_at_number"]] <- i
+    }
+
+    # return value
+    structure(do.call("rbind", out), 
+              class = c("margins", "data.frame"), 
+              at = if (is.null(at)) NULL else at_specification,
+              type = NULL,
+              call = if ("call" %in% names(model)) model[["call"]] else NULL,
+              model_class = class(model),
+              vce = vce, 
+              vcov = NULL,
+              jacobian = NULL,
+              weighted = FALSE,
+              iterations = NULL)
+}

tests/testthat/tests-reset_coefs.R

@@ -117,5 +117,20 @@ if (requireNamespace("lme4")) {
 #if (requireNamespace("nnet")) {
 #}
 
-#if (requireNamespace("ordinal")) {
-#}
+# if (requireNamespace("ordinal")) {
+    # test_that("reset_coefs() works for 'clm' objects", {
+        # data("wine", package = "ordinal")
+        # # base object
+        # mod1 <- ordinal::clm(rating ~ temp + contact, data = wine)
+        # # modified object
+        # mod2 <- reset_coefs(mod1, c(tempwarm = 3, contactyes = 2))
+        # # expect coefs to have been changed
+        # expect_true(!isTRUE(all.equal(coef(mod1), coef(mod2))), label = "coefficients reset in 'clm' object")
+        # # expect prediction from modified object to be correct
+        # expect_true(!isTRUE(all.equal(predict(mod1, newdata = wine), predict(mod2, newdata = wine))),
+                    # label = "predictions differ from original 'clm' object")
+        # expect_true(isTRUE(all.equal(predict(mod2, newdata = wine)$fit,
+                                     # 3*(as.integer(wine$temp)-1L) + 2*(as.integer(wine$contact)-1L), check.attributes = FALSE)),
+                    # label = "predictions correct from reset 'clm' object")
+    # })
+# }