У меня есть набор данных из 200 цыплят, отслеживаемых в течение 53 недель. Каждую курицу проверяют каждую неделю на наличие заболеваний, и они классифицируются как «1», если они не заражены, или «2», если они заражены. Таким образом, у меня есть матрица 200 на 53, S, 1 и 2, обозначающая статус инфекции каждой курицы на каждую неделю.
Я построил ряд байесовских моделей для расчета 2-х -2 матрица перехода, обозначающая вероятность перехода точки данных S [t-1, c] из состояния '1' или '2' в новое состояние в точке S [t, c] (Здесь t - наш индекс для времени, c - это наш индекс для курицы).
Один конкретный вариант модели предполагает, что эти вероятности перехода могут незначительно изменяться со временем. Мы предполагаем, что существует некоторая «средняя» вероятность перехода, альфа, а затем есть небольшой поправочный член для каждой из 53 недель.
В выражении с использованием JAGS в R мы пишем:
for(t in 1:Time){
# P12
pi[1,2,t] <- ilogit(alpha1 + chick_corr[1,t]) #ilogit forces this value between 0 and 1
# P11
pi[1,1,t] <- 1 - pi[1,2,t]
# P21
pi[2,1,t] <- ilogit(alpha2 + chick_corr[2,t])
# P22
pi[2,2,t] <- 1 - pi[2,1,t]
}
Тогда у меня есть два разных способа составления поправочных терминов.
В модели 0 мы считаем, что 53 члена chick_corr [1, t] взяты из нормального распределения с центром в 0. Аналогично для chick_corr [2, t].
в модели 1 вместо этого мы извлекаем оба из многомерного нормального распределения для поиска любой ковариации.
Обе модели выдают практически идентичный выходной сигнал с некоторыми хорошими результатами. Так что я тогда спросил бы, «какой из них лучше?». Нужна ли многомерная норма?
Конечно, практически не слишком важный вопрос для запроса, но я использую его в качестве инструмента для практики выбора модели с помощью байесовского фактора, используя многомерный подход MCM C, с использованием переменной `m ' в качестве переменной переключения для выбора между моделями.
Я записал свой код JAGS для выбора модели, однако при запуске апостериор для 'm' всегда будет полностью выбирать один вариант. Он никогда не выбирает другую модель. Я предоставил разумные псевдоприоры, чтобы сделать каждый выбор более «заманчивым». Любопытно также, что если я инициализирую цепочку в определенном состоянии переключателя c, он никогда не выйдет из этого состояния, предполагая, что он просто остается постоянно застрявшим в определенном состоянии.
Учитывая, что модели очень похожи, я предположил, что увижу апостериор для m со средним значением «0,5». Я также попытался заменить модель 1 на реплику модели 0, то есть сравнить две идентичные модели, которые наверняка имели бы равномерное распределение m, но наблюдалось то же явление.
Скорее всего, этот момент у меня есть это неправильно запрограммировано, так как я новичок в выборе байесовской модели. Ниже приведен код JAGS и некоторые дополнительные пояснения. Большое спасибо за любую помощь, которую вы можете предложить.
jags_model <- "model{
# Define how data is used
for(c in Chick){
for(t in 2:Time){
S[c, t] ~ dcat(pi_choice[S[c, t - 1], 1:2, t, (m+1)])
}
}
#pi_choice switches between the transition matrices for our two models
for(t in 1:Time){
pi_choice[1:2,1:2, t, 1] <- pi_M0[1:2,1:2,t]
pi_choice[1:2,1:2, t, 2] <- pi_M1[1:2,1:2,t]
}
# Inits
for(c in Chick){
S[c, 1] ~ dcat(Strt_Pr)
}
#MODEL 1
###############
for(t in 1:Time){
# P12
pi_M0[1,2,t] <- ilogit(alpha1_M0[m+1] + chick_corr_M0[1,t,m+1])
# P11
pi_M0[1,1,t] <- 1 - pi_M0[1,2,t]
# P21
pi_M0[2,1,t] <- ilogit(alpha2_M0[m+1] + chick_corr_M0[2,t,m+1])
# P22
pi_M0[2,2,t] <- 1 - pi_M0[2,1,t]
}
for(t in 1:Time){
chick_corr_M0[1,t,1] ~ dnorm(chick1mean[1], 1/chick1sd[1]) #prior
chick_corr_M0[2,t,1] ~ dnorm(chick2mean[1], 1/chick2sd[1]) #prior
chick_corr_M0[1,t,2] ~ dnorm(chick1mean[2], 1/chick1sd[2]) #pseudo
chick_corr_M0[2,t,2] ~ dnorm(chick2mean[2], 1/chick2sd[2]) #pseudo
}
chick1mean[1] ~ dnorm(0, 0.001) #prior
chick1sd[1] ~ dunif(0, 500) #prior
chick2mean[1] ~ dnorm(0, 0.001) #prior
chick2sd[1] ~ dunif(0, 500) #prior
chick1mean[2] ~ dnorm(-3.4, 1/(0.7*(3.4-0.6))) #PSEUDO
chick1sd[2] ~ dunif(0.05,0.3) #PSEUDO
chick2mean[2] ~ dnorm(-11.4, 1/(0.7*(11.4-6.6))) #PSEUDO
chick2sd[2] ~ dunif(0.1,0.7) #PSEUDO
alpha1_M0[1] ~ dunif(-50,50)
alpha2_M0[1] ~ dunif(-50,50)
alpha1_M0[2] ~ dnorm(2.8,1/(0.7*(5.4-2.8))) #PSEUDO
alpha2_M0[2] ~ dnorm(11.6,1/(0.7*(21.8-11.6))) #PSEUDO
#Now for Model 1 , using the MVN
######################
for(t in 1:Time){
# P12
pi_M1[1,2,t] <- ilogit(alpha1_M1[m+1] + chick_corr_M1[1,t,m+1])
# P11
pi_M1[1,1,t] <- 1 - pi_M1[1,2,t]
# P21
pi_M1[2,1,t] <- ilogit(alpha2_M1[m+1] + chick_corr_M1[2,t,m+1])
# P22
pi_M1[2,2,t] <- 1 - pi_M1[2,1,t]
}
for(t in 1:Time){
chick_corr_M1[1:2,t,1] ~ dmnorm(Mu1, CovMat1[1:2,1:2]) #PSEUDO
chick_corr_M1[1:2,t,2] ~ dmnorm(Mu2, CovMat2[1:2,1:2])
}
CovMat1 ~ dwish((0.5*CovPseudo[1:2,1:2]) , 2) #PSEUDO
CovPseudo[1,1] <- 6.18
CovPseudo[1,2]<- 1.83
CovPseudo[2,1]<- 1.83
CovPseudo[2,2] <- 3.67
CovMat2 ~ dwish(WishPr[1:2,1:2] , 2) #Prior
alpha1_M1[1] ~ dnorm(-0.63, 1/(0.7*(0.63-0.45))) #PSEUDO
alpha2_M1[1] ~ dnorm(0.23,1/(0.7*(0.48-0.23))) #PSEUDO
alpha1_M1[2] ~ dunif(-50,50) #PRIOR
alpha2_M1[2] ~ dunif(-50,50) #PRIOR
m ~ dbern(0.5)
}"
Strt_Pr <- rep(1/2, 2)
model_output <- run.jags(model = jags_model, monitor = c("pi_M0", "pi_M1", "alpha1_M0", "alpha2_M0",
"chick_corr_M0", "chick_corr_M1", "alpha1_M1",
"alpha2_M1", "m"),
data = list(S = S,
Chick = 1:nrow(S),
Time = ncol(S),
Strt_Pr = Strt_Pr,WishPr = diag(x=1,nrow=2),
Mu1 = c(-0.003,0.0037), Mu2 = c(0,0)),
sample = 100000)
Чтобы подвести итог, проблема не переключается. Модель остается с любым состоянием, в котором инициализировано m.