如何获取多项式logit的置信区间？

Question

让我以UCLA的多项式逻辑回归为例子——

library(nnet)
library(foreign)

ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")

我想知道如何获取95％的置信区间？

- user2966726

---

我的意思是预测均值周围的95%置信区间。 - user2966726

Answer 1

这可以通过使用“effects”软件包来实现，我在Cross Validated的另一个问题中展示了这个软件包，链接在这里。

让我们看一下你的例子。

library(nnet)
library(foreign)

ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

我们使用 effects 库的 Effect() 方法代替了 base 库的 predict() 方法。

require(effects)

fit.eff <- Effect("ses", test, given.values = c("write" = mean(ml$write)))

data.frame(fit.eff$prob, fit.eff$lower.prob, fit.eff$upper.prob)

  prob.academic prob.general prob.vocation L.prob.academic L.prob.general L.prob.vocation U.prob.academic
1     0.4396845    0.3581917     0.2021238       0.2967292     0.23102295      0.10891758       0.5933996
2     0.4777488    0.2283353     0.2939159       0.3721163     0.15192359      0.20553211       0.5854098
3     0.7009007    0.1784939     0.1206054       0.5576661     0.09543391      0.05495437       0.8132831
  U.prob.general U.prob.vocation
1      0.5090244       0.3442749
2      0.3283014       0.4011175
3      0.3091388       0.2444031

如果需要的话，我们也可以使用effects中的工具绘制预测概率及其相应的置信区间。

plot(fit.eff)

- Johan Larsson

---

Johan，当你只想为截距模型（例如，公式为：ses〜1）获取置信区间时，是否有一种使用effect函数的方法？似乎仅在effect函数中写入“1”或“Intercept”并不起作用。 - Jamie

Answer 2

只需在您的模型对象上使用confint函数即可。

ci <- confint(test, level=0.95)

请注意，confint是一个通用函数，在multinom中运行的是特定版本，可以通过运行来查看。

> methods(confint)
[1] confint.default   confint.glm*      confint.lm*       confint.multinom*
[5] confint.nls* 

编辑：

关于计算预测概率的置信区间一事，我引用自：https://stat.ethz.ch/pipermail/r-help/2004-April/048917.html

有没有可能使用multinom函数估计概率的置信区间？

不行，因为置信区间（sic）只适用于单个参数而非概率（sic）。预测是一个概率分布，因此不确定性应该是Kd空间中的某个区域，而不是一个区间。你为什么需要关于预测的不确定性陈述（通常称为容限区间/区域）？在这种情况下，你有一个发生或不发生的事件，有意义的不确定性是概率分布。如果你真的需要置信区间，可以从拟合参数的不确定性中进行模拟，预测并总结所得到的经验分布。

- nico

---

我需要的不是多项式系数的置信区间。相反，我需要预测概率的置信区间 - predict（test，newdata = dses，“probs”）。 - user2966726

不理解为什么 Ripley 教授（上述引语的作者）对此大惊小怪。在不深入思考的情况下，似乎他的反对意见适用于任何“置信区间”（或者你想称其为什么都可以），但这是一种相当标准的做法。我调查了 nnet:::predict.multinom、nnet:::predict.nnet 和 nnet:::vcov.multinom。如果你能构造出模型矩阵 X，使得预测的概率是 X %*% coef(model)，那么方差（在 logit 尺度上）就是 X %*% vcov(model) %*% t(X)，然后你可以在 logit 尺度上构建置信区间并进行反变换。 - Ben Bolker

然而，上述函数并不够透明，让我认为我无法在短时间内对它们进行黑客攻击。如果您可以证明一个序数响应，ordinal::predict.clm()确实允许使用se.fit选项。 - Ben Bolker

@BenBolker，emmeans和marginaleffects软件包的开发版本支持nnet::multinom和mclogit::mblogit多项式模型，并且可以通过Delta方法在响应比例上计算95％的置信区间。从logit或潜在（中心化logit）比例尺通过softmax变换反向转换原则上应该更好（多项式的反向链接函数是softmax，而不是逆logit），但仍需在marginaleffects软件包中实现此功能... - Tom Wenseleers