Логистическая регрессия с категориальными предикторами с использованием JAGS

Question

Я новичок в JAGS и пытаюсь предсказать двоичный результат (0/1), используя 9 не непрерывных предикторов. Значения предиктора могут быть 0, 1 или 2. Это мой первый раз, и хотя я могу запустить модель, я на 100% уверен, что здесь определенно есть ряд проблем.

Пример файла данных (список)

$y
[1] 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0

$N
[1] 50

$oAnt
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1
[29] 1 1 1 1 1 1 2 1 1 0 1 1 1 1 1 2 1 1 1 1 1 1

$nAnt
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

$cAnt
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

$oPen
[1] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 0
[29] 1 0 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 0 2 1 1

$nPen
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1

$cPen
[1] 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0
[29] 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0

$oFin
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

$nFin
[1] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1
[29] 1 1 1 1 2 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 3

$cFin
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1
[29] 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0

модель

model {
    for( i in 1 : N ){
        y[i] ~ dbern(mu[i])
        mu[i] <- 1/(1+exp(-(b0 + b1*oAnt[i] + b2*nAnt[i] + b3*cAnt[i] + b4*oPen[i] + b5*nPen[i] + b6*cPen[i] + b7*oFin[i] + b8*nFin[i] + b9*cFin[i])))
}
b0 ~ dnorm(0, 1.0e-12)
b1 ~ dnorm(0, 1.0e-12)
b2 ~ dnorm(0, 1.0e-12)
b3 ~ dnorm(0, 1.0e-12)
b4 ~ dnorm(0, 1.0e-12)
b5 ~ dnorm(0, 1.0e-12)
b6 ~ dnorm(0, 1.0e-12)
b7 ~ dnorm(0, 1.0e-12)
b8 ~ dnorm(0, 1.0e-12)
b9 ~ dnorm(0, 1.0e-12)
}

Я использовал оценки из glm() модель в виде единиц (как предполагает А. Гельман) - но для простоты, давайте просто предположим, что я позволю JAGS выбрать начальные значения для цепочек.

Бегущая модель и т. Д.

jagsModel = jags.model(file = "antPenFin.txt", data = dataList, n.chains = 2, n.adapt = 500)

update(jagsModel, n.iter = 500)

codaSamples = coda.samples(jagsModel,
variable.names = names(dataList)[3:11], n.iter = 5000)

Проблемы

Вывод моей модели выглядит совершенно нечетким (что становится ясно, когда я пытаюсь построить ее). Я уверен, что здесь есть какая-то очень простая проблема. Может ли кто-нибудь помочь?

Большое спасибо.

Answer 1

Ваша проблема здесь в том, что вы не можете предоставить диапазон чисел (от 0 до 3 в вашем случае) в качестве категориальных ковариат. В настоящее время ваша модель интерпретирует эти числа как непрерывные. Что вам нужно сделать, это преобразовать их в фиктивные переменные. Вы можете сделать это легко используя model.matrix функция в R. Я собираюсь сгенерировать некоторые данные в качестве примера, которые вы могли бы затем применить к своим данным.

# generate y
y <- sample(0:1, 30, replace = TRUE)
# add three different categorical covariates. Note here that these are all
# factors.
oAnt <- factor(sample(0:2, 30, replace = TRUE))
cAnt <- factor(sample(0:1, 30, replace = TRUE))
nFin <- factor(sample(0:3, 30, replace = TRUE))

# create your model matrix
my_matrix <- model.matrix(y ~ oAnt + cAnt + nFin)

head(my_matrix)

   (Intercept) oAnt1 oAnt2 cAnt1 nFin1 nFin2 nFin3
1           1     1     0     1     0     1     0
2           1     1     0     1     1     0     0
3           1     1     0     0     0     0     1
4           1     1     0     1     0     0     1
5           1     0     0     1     0     1     0
6           1     1     0     1     0     0     1

Для каждого уровня фактора для предиктора вам нужно будет создать n-1 столбцов (например, nFin колеблется от 0 до 3, вы создадите 3 столбца фиктивных переменных). Таким образом, у вас будет еще много коэффициентов регрессии для включения в вашу модель. Как только вы сделаете свою матрицу модели, вы можете преобразовать ее в список как таковой.

# remove the intercept from the list as you already have it in your model
matrix_list <- as.list(data.frame(my_matrix[,-1]))

Оттуда все, что вам нужно сделать, это создать еще несколько предикторов для модели. Кроме того, если nAnt в вашей модели действительно все, идите и удалите это, так как вы просто включаете два перехвата.