Как получить полиномиальные вероятности в 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], но пока безуспешно.

Любые советы будут очень благодарны!

(Я предоставил проблемные разделы модели; вся она довольно длинная.)


regression statistics multinomial winbugs winbugs14
person Sprog    schedule 13.02.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Обычный способ сделать это — использовать полиномиальный эквивалент логит-ссылки, чтобы ограничить преобразованные вероятности интервалом (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 (и немного переформулировать код, чтобы фи было кумулятивным), чтобы получить логистическая регрессия пропорциональных шансов (POLR).

Редактировать

Вот ссылка на пример кода, охватывающего логистическую регрессию, полиномиальную регрессию и POLR из курса, который я преподаю:

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

Обратите внимание, что он использует JAGS (а не WinBUGS), но, насколько я знаю, нет различий в синтаксисе моделей для этих типов моделей. Если вы хотите быстро начать работу с ранджагами и JAGS на фоне WinBUGS, вы можете следовать этой виньетке:

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

person Matt Denwood    schedule 16.02.2018
comment
Действительно полезный ответ. Спасибо - было бы здорово взглянуть на пример кода, который вы упомянули! - person Sprog; 18.02.2018