如何使用多元回归在 WinBUGS 中获得多项式概率

How to obtain multinomial probabilities in WinBUGS with multiple regression

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

任何提示将不胜感激!

(我已经提供了模型中有问题的部分;整个过程相当长。)

执行此操作的通常方法是使用 logit link 的多项式等效项将转换后的概率限制在区间 (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]

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


编辑

这里是 link 一些示例代码,涵盖我教授的课程中的逻辑回归、多项回归和 POLR:

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

请注意,它使用 JAGS(而不是 WinBUGS),但据我所知,这些类型的模型在模型语法上没有差异。如果您想从 WinBUGS 背景快速开始使用 runjags 和 JAGS,那么您可以遵循这个小插图:

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