statistics - 如何通过多元回归获得 WinBUGS 中的多项概率

Question

在 WinBUGS 中，我指定一个具有多项似然函数的模型，并且我需要确保多项概率均在 0 到 1 之间且总和为 1。

这是指定可能性的代码部分:

e[k,i,1:9] ~ dmulti(P[k,i,1:9],n[i,k]) 

这里，数组 P[] 指定多项分布的概率。

这些概率将根据我的数据(矩阵 e[])使用一系列固定和随机效应的多元线性回归来估计。例如，以下是用于预测 P[] 元素之一的多元线性回归:

P[k,1,2] <- intercept[1,2]  +  Slope1[1,2]*Covariate1[k]  +
Slope2[1,2]*Covariate2[k]  +  Slope3[1,2]*Covariate3[k]  
+ Slope4[1,2]*Covariate4[k] +  RandomEffect1[group[k]]  +  
RandomEffect2[k]

编译时，模型产生错误:

elements of proportion vector of multinomial e[1,1,1] must be between zero and one

如果我理解正确的话，这意味着向量 P[k,i,1:9](上面多项式似然函数中的概率向量)的元素可能是非常大(或小)的数字。实际上，它们都需要在 0 和 1 之间，并且总和为 1。

我是 WinBUGS 的新手，但通过阅读周围的内容，似乎使用 beta 回归而不是多重线性回归可能是前进的方向。然而，虽然这允许每个元素介于 0 和 1 之间，但这似乎并没有触及问题的核心，即 P[k,i,1:9] 的所有元素都必须为正且总和为 1。

响应变量可能可以非常简单地转换为比例。我尝试过将每个元素除以 P[k,i,1:9] 的总和，但到目前为止没有成功。

如有任何提示，我们将不胜感激!

(我已经提供了模型的有问题的部分；整个事情相当长。)

Answer 1

执行此操作的通常方法是使用 logit 链接的多项式等价项将转换后的概率限制在区间 (0,1) 内。例如(对于单个预测器，但对于您需要的任意多个预测器来说，原理相同):

Response[i, 1:Categories] ~ dmulti(prob[i, 1:Categories], Trials[i])

phi[i,1] <- 1
prob[i,1] <- 1 / sum(phi[i, 1:Categories])
for(c in 2:Categories){
    log(phi[i,c]) <- intercept[c] + slope1[c] * Covariate1[i]
    prob[i,c] <- phi[i,c] / sum(phi[i, 1:Categories])
}

为了可识别性，phi[1] 的值设置为 1，但截距和斜率 1 的其他值是独立估计的。当类别数等于 2 时，这会崩溃为通常的逻辑回归，但编码为多项式响应:

log(phi[i,2]) <- intercept[2] + slope1[2] * Covariate1[i]
prob[i,2] <- phi[i, 2] / (1 + phi[i, 2])
prob[i,1] <- 1 / (1 + phi[i, 2])

即:

logit(prob[i,2]) <- intercept[2] + slope1[2] * Covariate1[i]
prob[i,1] <- 1 - prob[i,2]

在此模型中，我按类别对斜率 1 进行了索引，这意味着结果的每个级别都与预测变量具有独立的关系。如果您有一个序数响应，并且想要假设与协变量相关的优势比在响应的连续水平之间是一致的，那么您可以将索引放在斜率 1 上(并稍微重新编写代码，以便 phi 是累积的)以获得比例优势逻辑回归 (POLR)。

---

编辑

这里是一些示例代码的链接，涵盖我教授的类(class)中的逻辑回归、多项式回归和 POLR:

http://runjags.sourceforge.net/examples/squirrels.R

请注意，它使用 JAGS(而不是 WinBUGS)，但据我所知，这些类型的模型的模型语法没有差异。如果您想从 WinBUGS 背景快速开始使用 runjags 和 JAGS，那么您可以按照以下小插图进行操作:

http://runjags.sourceforge.net/quickjags.html