2

我正在尝试估计一个有序的 logit 模型，包括。通过遵循本教程中的代码，R 中的边际效应。我正在使用polrfrom the MASSpackage 来估计模型并ocME从erer package 来尝试计算边际效应。

估计模型没有问题。

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.

那么有人可以帮助我理解我做错了什么吗？我在这里发布了包含模型的两个变量的数据子集发布了包含模型的两个变量的数据子集。在 RI 中将 DV 设置为因子变量，IV 是连续的。

边注：

我可以将计算从 R 传递给 Stata，RStata以计算边际效应而没有任何问题。但我不想定期这样做，所以我想了解是什么导致了 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
------------------------------------------------------------------------------
4

1 回答 1

2

您的模型只有一个解释变量 ( 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
于 2018-12-05T13:20:27.133 回答