使用 ocME 在 R 中计算有序 logit 模型的边际效应的问题

Question

我正在尝试估计一个有序的 logit 模型，包括。通过遵循代码from this tutorial 在 R 中的边际效应。我使用MASS 包中的polr 来估计模型，并使用erer 包中的ocME 来尝试计算边际效应。

估计模型没问题。

logitModelSentiment90 <- polr(availability_90_ord ~ mean_sentiment, data = data, Hess = T,
                           method = "logistic")

但是，我遇到了ocME 的问题，它会生成以下错误消息：

ocME(logitModelSentiment90)

Error in eval(predvars, data, env) : 
numeric 'envir' arg not of length one

ocME 下面的文档指出，应该使用的对象需要来自 polr 函数，这似乎正是我正在做的。

ocME(w, rev.dum = TRUE, digits = 3)
w = an ordered probit or logit model object estimated by polr from the MASS library.

那么任何人都可以帮助我了解我做错了什么吗？我已经发布了包含模型here 的两个变量的数据子集。在 R 中，我将 DV 设置为因子变量，IV 是连续的。

旁注：

我可以使用RStata 将计算从 R 传递给 Stata，以毫无问题地计算边际效应。但我不想定期这样做，所以我想了解导致 R 和 ocME 出现问题的原因。

stata("ologit availability_90_ord mean_sentiment
  mfx", data.in = data)
. ologit availability_90_ord mean_sentiment

Iteration 0:   log likelihood = -15379.121  
Iteration 1:   log likelihood = -15378.742  
Iteration 2:   log likelihood = -15378.742  

Ordered logistic regression                     Number of obs     =     11,901
                                                LR chi2(1)        =       0.76
                                                Prob > chi2       =     0.3835
Log likelihood = -15378.742                     Pseudo R2         =     0.0000

------------------------------------------------------------------------------
avail~90_ord |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
mean_senti~t |   .0044728   .0051353     0.87   0.384    -.0055922    .0145379
-------------+----------------------------------------------------------------
       /cut1 |   -1.14947   .0441059                     -1.235916   -1.063024
       /cut2 |  -.5286239    .042808                     -.6125261   -.4447217
       /cut3 |   .3127556   .0426782                      .2291079    .3964034
------------------------------------------------------------------------------
.       mfx

Marginal effects after ologit
      y  = Pr(availability_90_ord==1) (predict)
         =  .23446398
------------------------------------------------------------------------------
variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X
---------+--------------------------------------------------------------------
mean_s~t |  -.0008028      .00092   -0.87   0.384  -.002609  .001004   7.55768
------------------------------------------------------------------------------

Answer 1

您的模型只有一个解释变量 (mean_sentiment)，这对于 ocME 来说似乎是个问题。例如，尝试向模型添加第二个变量：

logitModelSentiment90 <- polr(availability_90_ord ~ mean_sentiment + I(mean_sentiment^2),
                              data = data, Hess = T,  method = "logistic")
ocME(logitModelSentiment90)

#                     effect.0 effect.1 effect.2 effect.3
# mean_sentiment        -0.004   -0.001        0    0.006
# I(mean_sentiment^2)    0.000    0.000        0    0.000

稍加修改ocME 也可以使用一个自变量正确运行。
试试下面的myocME函数

myocME <- function (w, rev.dum = TRUE, digits = 3) 
{
    if (!inherits(w, "polr")) {
        stop("Need an ordered choice model from 'polr()'.\n")
    }
    if (w$method != "probit" & w$method != "logistic") {
        stop("Need a probit or logit model.\n")
    }
    lev <- w$lev
    J <- length(lev)
    x.name <- attr(x = w$terms, which = "term.labels")
    x2 <- w$model[, x.name, drop=FALSE]
    ww <- paste("~ 1", paste("+", x.name, collapse = " "), collapse = " ")
    x <- model.matrix(as.formula(ww), data = x2)[, -1, drop=FALSE]
    x.bar <- as.matrix(colMeans(x))
    b.est <- as.matrix(coef(w))
    K <- nrow(b.est)
    xb <- t(x.bar) %*% b.est
    z <- c(-10^6, w$zeta, 10^6)
    pfun <- switch(w$method, probit = pnorm, logistic = plogis)
    dfun <- switch(w$method, probit = dnorm, logistic = dlogis)
    V2 <- vcov(w)
    V3 <- rbind(cbind(V2, 0, 0), 0, 0)
    ind <- c(1:K, nrow(V3) - 1, (K + 1):(K + J - 1), nrow(V3))
    V4 <- V3[ind, ]
    V5 <- V4[, ind]
    f.xb <- dfun(z[1:J] - c(xb)) - dfun(z[2:(J + 1)] - c(xb))
    me <- b.est %*% matrix(data = f.xb, nrow = 1)
    colnames(me) <- paste("effect", lev, sep = ".")
    se <- matrix(0, nrow = K, ncol = J)
    for (j in 1:J) {
        u1 <- c(z[j] - xb)
        u2 <- c(z[j + 1] - xb)
        if (w$method == "probit") {
            s1 <- -u1
            s2 <- -u2
        }
        else {
            s1 <- 1 - 2 * pfun(u1)
            s2 <- 1 - 2 * pfun(u2)
        }
        d1 <- dfun(u1) * (diag(1, K, K) - s1 * (b.est %*% t(x.bar)))
        d2 <- -1 * dfun(u2) * (diag(1, K, K) - s2 * (b.est %*% 
            t(x.bar)))
        q1 <- dfun(u1) * s1 * b.est
        q2 <- -1 * dfun(u2) * s2 * b.est
        dr <- cbind(d1 + d2, q1, q2)
        V <- V5[c(1:K, K + j, K + j + 1), c(1:K, K + j, K + j + 
            1)]
        cova <- dr %*% V %*% t(dr)
        se[, j] <- sqrt(diag(cova))
    }
    colnames(se) <- paste("SE", lev, sep = ".")
    rownames(se) <- colnames(x)
    if (rev.dum) {
        for (k in 1:K) {
            if (identical(sort(unique(x[, k])), c(0, 1))) {
                for (j in 1:J) {
                  x.d1 <- x.bar
                  x.d1[k, 1] <- 1
                  x.d0 <- x.bar
                  x.d0[k, 1] <- 0
                  ua1 <- z[j] - t(x.d1) %*% b.est
                  ub1 <- z[j + 1] - t(x.d1) %*% b.est
                  ua0 <- z[j] - t(x.d0) %*% b.est
                  ub0 <- z[j + 1] - t(x.d0) %*% b.est
                  me[k, j] <- pfun(ub1) - pfun(ua1) - (pfun(ub0) - 
                    pfun(ua0))
                  d1 <- (dfun(ua1) - dfun(ub1)) %*% t(x.d1) - 
                    (dfun(ua0) - dfun(ub0)) %*% t(x.d0)
                  q1 <- -dfun(ua1) + dfun(ua0)
                  q2 <- dfun(ub1) - dfun(ub0)
                  dr <- cbind(d1, q1, q2)
                  V <- V5[c(1:K, K + j, K + j + 1), c(1:K, K + 
                    j, K + j + 1)]
                  se[k, j] <- sqrt(c(dr %*% V %*% t(dr)))
                }
            }
        }
    }
    t.value <- me/se
    p.value <- 2 * (1 - pt(abs(t.value), w$df.residual))
    out <- list()
    for (j in 1:J) {
        out[[j]] <- round(cbind(effect = me[, j], error = se[, 
            j], t.value = t.value[, j], p.value = p.value[, j]), 
            digits)
    }
    out[[J + 1]] <- round(me, digits)
    names(out) <- paste("ME", c(lev, "all"), sep = ".")
    result <- listn(w, out)
    class(result) <- "ocME"
    return(result)
}

并运行以下代码：

logitModelSentiment90 <- polr(availability_90_ord ~ mean_sentiment, 
                              data = data, Hess = T,  method = "logistic")   
myocME(logitModelSentiment90)

#                effect.0 effect.1 effect.2 effect.3
# mean_sentiment   -0.001        0        0    0.001