Как получить полиномиальные вероятности в WinBUGS с множественной регрессией
В 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, но из прочтения кажется, что каким-то образом использование бета-регрессии, а не множественных линейных регрессий, могло бы стать шагом вперед. Однако, хотя это позволило бы каждому элементу быть между 0 и 1, это, похоже, не доходит до сути проблемы, которая заключается в том, что все элементы P [k, i, 1: 9] должны быть положительными и сумма до 1.
Может случиться так, что переменная отклика может быть очень просто преобразована в пропорцию. Я попробовал это, пытаясь разделить каждый элемент на сумму P[k,i,1:9], но пока безуспешно.
Любые советы будут очень благодарны!
(Я предоставил проблемные разделы модели; все это довольно долго.)
1 ответ
Обычный способ сделать это состоит в том, чтобы использовать полиномиальный эквивалент логит-ссылки, чтобы ограничить преобразованные вероятности интервалом (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, но другие значения intercept и slope1 оцениваются независимо. Когда количество категорий равно 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).
редактировать
Вот ссылка на пример кода, охватывающего логистическую регрессию, полиномиальную регрессию и POLR из курса, который я преподаю:
http://runjags.sourceforge.net/examples/squirrels.R
Обратите внимание, что он использует JAGS (а не WinBUGS), но, насколько мне известно, нет различий в синтаксисе модели для этих типов моделей. Если вы хотите быстро начать работу с runjags & JAGS из фона WinBUGS, вы можете использовать следующую виньетку: