R中多项式回归的拟合值：参考类别的系数是什么？

Question

我正在使用nnet包中的multinom函数运行多项式逻辑回归。

在多项式逻辑回归中，据我所知，系数是响应概率与参考响应概率比值的对数变化量（即，ln(P(i)/P(r))=B₁+B₂*X... 其中i是一个响应类别，r是参考类别，X是一些预测因素）。

然而，fitted(multinom(...))为每个类别（包括参考类别r）生成估计值。 编辑：例如：

set.seed(1)
library(nnet)
DF <- data.frame(X = as.numeric(rnorm(30)), 
                 Y = factor(sample(letters[1:5],30, replace=TRUE)))
DF$Y<-relevel(DF$Y, ref="a") #ensure a is the reference category
model <- multinom(Y ~ X, data = DF)
coef(model)
#  (Intercept)           X
#b   0.1756835  0.55915795
#c  -0.2513414 -0.31274745
#d   0.1389806 -0.12257963
#e  -0.4034968  0.06814379

head(fitted(model))
#          a         b          c         d         e
#1 0.2125982 0.2110692 0.18316042 0.2542913 0.1388810
#2 0.2101165 0.1041655 0.26694618 0.2926508 0.1261210
#3 0.2129182 0.2066711 0.18576567 0.2559369 0.1387081
#4 0.1733332 0.4431170 0.08798363 0.1685015 0.1270647
#5 0.2126573 0.2102819 0.18362323 0.2545859 0.1388516
#6 0.1935449 0.3475526 0.11970164 0.2032974 0.1359035

head(DF)
#           X Y
#1 -0.3271010 a

为了计算第一行响应 b 和响应 a 之间的预测概率比率，我们计算 exp(0.1756835+0.55915795*(-0.3271010))=0.9928084。我看到这对应于第一行拟合的 P(b)/P(a) (0.2110692/0.2125982=0.9928084)。

参考类别的拟合概率是通过代数方式计算得出的吗 (例如， 0.2110692/exp(0.1756835+0.55915795*(-0.3271010)))？

有没有办法获得参考类别预测概率的方程式？

- jflournoy

---

我在你的问题中添加了一些数据。你能否使用这个例子澄清你的问题？期望的输出是什么？ - agstudy

除了缺乏数据示例或从中提取函数“multinom”的包规范外，我认为您错误地指定了如何表示赔率比。赔率是“P（i）/ P（！i）”，因此赔率比具有4个概率表达式的比率形式。 - IRTFM

DWin，关于二项式逻辑回归，您是正确的。然而，我认为在多项式逻辑回归中，预测的赔率比实际上是某个响应概率与参考响应概率之间的比率--请参阅适当的Wikipedia页面。 - jflournoy

赔率是一个事件的概率（可能取决于其他特征）除以非事件的概率（在相同特征条件下）。维基百科的引用并不改变这个事实，尤其是当维基百科关于“赔率比”和“赔率”的部分有正确的公式时。 - IRTFM

DWin，感谢您对我的术语进行纠正。我应该说“概率比”而不是“赔率”。 - jflournoy

Answer 1

我有同样的问题，在查找后我认为解决方案是： 给定三个类：a，b，c和算法拟合（model）的概率pa，pb，pc输出，您可以从以下三个方程式中重构这些概率：

log(pb/pa) = beta1*X

log(pc/pa) = beta2*X

pa+pb+pc=1

其中beta1,beta2是coef(model)的输出行，X是您的输入数据。

通过操作这些方程式，您可以得到：

pb = exp(beta1*X)/(1+exp(beta1*X)+exp(beta2*X))

pc = exp(beta2*X)/(1+exp(beta1*X)+exp(beta2*X))

pa = 1 - pb - pc

- Gabriela

---

我认为上面答案中的这些方程式缺少了截距系数，需要加以考虑。请查看下面我的回答以获取更多详细信息。 - Nicholas G Reich

Answer 2

重点在于multinom()的帮助文件中指出，“拟合对数线性模型，第一类系数为零。”

这意味着可以直接计算参考类别的预测值，并假定“a”类的系数为零。例如，对于给定的样本行，可以使用softmax变换计算类别“a”的预测概率： exp(0+0)/(exp(0+0) + exp(0.1756835 + 0.55915795*(-0.3271010)) + exp(-0.2513414 + (-0.31274745)*(-0.3271010)) + exp(0.1389806 + (-0.12257963)*(-0.3271010)) + exp(-0.4034968 + 0.06814379*(-0.3271010))) 或者更简单地说，使用非硬编码数值，我们可以计算整个数据的第一行的所有概率：

softMax <- function(x){
    expx <- exp(x)
    return(expx/sum(expx))
}
coefs <- rbind(c(0,0), coef(model))
linear.predictor <- as.vector(coefs%*%c(1,-0.3271010))
softMax(linear.predictor)

值得一提的是：原问题中的示例对我来说无法完全复制，我的随机种子会给出不同的随机数。因此，我已经重新生成了示例，并在下面的计算中进行了说明。

library(nnet)
set.seed(1)
DF <- data.frame(
    X = as.numeric(rnorm(30)), 
    Y = factor(sample(letters[1:5],30, replace=TRUE)))
DF$Y<-relevel(DF$Y, ref="a") #ensure a is the reference category
model <- multinom(Y ~ X, data = DF)
coef(model)
##   (Intercept)             X
## b -0.33646439  1.200191e-05
## c -0.36390688 -1.773889e-01
## d -0.45197598  1.049034e+00
## e -0.01418543  3.076309e-01

DF[1,]
##            X Y
## 1 -0.6264538 c

fitted.values(model)[1,]
##          a          b          c          d          e 
## 0.27518921 0.19656378 0.21372240 0.09076844 0.22375617 

coefs <- rbind(c(0,0), coef(model))
linear.predictor <- as.vector(coefs%*%c(1,DF[1,"X"]))
softMax(linear.predictor)
## [1] 0.27518921 0.19656378 0.21372240 0.09076844 0.22375617

- Nicholas G Reich