2

在 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] 的总和，但到目前为止还没有成功。

任何提示将不胜感激！

（我已经提供了模型有问题的部分；整个内容相当长。）

4

1 回答 1

2

执行此操作的常用方法是使用 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])

IE：

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

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

编辑

这是我教授的课程中包含逻辑回归、多项回归和 POLR 的示例代码的链接：

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

请注意，它使用 JAGS（而不是 WinBUGS），但据我所知，这些类型的模型的模型语法没有区别。如果你想从 WinBUGS 背景快速开始使用 runjags 和 JAGS，那么你可以遵循这个小插曲：

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

于 2018-02-16T09:49:31.057 回答