Вставка кода Jagam в модель Runjags (JAGS)
Я пытался включить сглаживание в модель runjags, которую я создал, чтобы смоделировать номера нор морских птиц и распределение по острову. Мне удалось сгенерировать некоторый код сглаживания, извлекая данные подсчета и координаты x и y из выходных данных модели и используя руководство JAGAM на этой странице http://www.petrkeil.com/?p=2385
Я думаю, что смогу улучшить производительность модели, включив сглаживание в модель зазубрин, но не знаю, как это сделать. Можете ли вы предложить мне какие-либо советы о том, как этого добиться?
Я приложил раздел кода runjags и вывода JAGAM ниже.
код runjags:
for(i in 1:K) {
S1[i]~dpois(lambda1[i])
SS1[i]~dpois(lambda1[i])
lambda1[i]<-exp(a0+
a1*Tussac[i]+
a2*normalise_DEM_aspect[i]+
a3*normalise_DEM_slope[i]+
a4*Tussac[i]*normalise_DEM_aspect[i]+
a5*Tussac[i]*normalise_DEM_slope[i]+
a6*normalise_sentinel1[i]+
a7*normalise_setinel3[i]+
a8*normalise_sentinel4[i]+
a9*normalise_sentinel5[i]+
a10*normalise_sentinel8[i]+
a11*normalise_sentinel10[i]+
a12*S2[i])
}
Выход JAGAM:
readLines("jagam.bug")
"model {"
" eta <- X %*% b ## linear predictor"
" for (i in 1:n) { mu[i] <- exp(eta[i]) } ## expected response"
" for (i in 1:n) { y[i] ~ dpois(mu[i]) } ## response "
" ## Parametric effect priors CHECK tau=1/35^2 is appropriate!"
" for (i in 1:1) { b[i] ~ dnorm(0,0.00083) }"
" ## prior for s(x,y)... "
" K1 <- S1[1:29,1:29] * lambda[1] + S1[1:29,30:58] * lambda[2]"
" b[2:30] ~ dmnorm(zero[2:30],K1) "
" ## smoothing parameter priors CHECK..."
" for (i in 1:2) {"
" lambda[i] ~ dgamma(.05,.005)"
" rho[i] <- log(lambda[i])"
" }"
"}"
Пример данных:
S1 Logit_tussac soil_moisture DEM_slope DEM_aspect DEM_elevation sentinel1 sentinel2 sentinel3 sentinel4 sentinel5 sentinel6 sentinel7 sentinel8 sentinel9 sentinel10
NA NA NA 14.917334 256.1612 12.24432 0.0513 0.0588 0.0541 0.1145 0.1676 0.1988 0.1977 0.1658 0.1566 0.0770
0 -9.210240 1 23.803741 225.1231 16.88028 0.1058 0.1370 0.2139 0.2387 0.2654 0.2933 0.3235 0.2928 0.3093 0.1601
NA NA NA 20.789165 306.0945 18.52480 0.0287 0.0279 0.0271 0.0276 0.0290 0.0321 0.0346 0.0452 0.0475 0.0219
NA -9.210240 1 6.689442 287.9641 36.08975 0.0462 0.0679 0.1274 0.1535 0.1797 0.2201 0.2982 0.2545 0.4170 0.2252
0 -9.210240 1 25.476444 203.0659 23.59964 0.0758 0.1041 0.1326 0.1571 0.2143 0.2486 0.2939 0.2536 0.3336 0.1937
1 -1.385919 3 1.672511 270.0000 39.55215 0.0466 0.0716 0.1227 0.1482 0.2215 0.2715 0.3334 0.2903 0.3577 0.1957
1 ответ
Это очень хороший вопрос и хорошая идея использовать (чрезвычайно полезный) вывод jagam, чтобы добавить термин GAM в вашу модель. В вашем случае я бы порекомендовал использовать jagam для генерации только термина (ов) GAM и ничего больше (даже перехвата), а затем скопировать / вставить соответствующие разделы вывода модели jagam в существующий код модели как а также взять переменную данных X из jagam и использовать ее в качестве своих данных. Это проще всего продемонстрировать на примере:
Сначала смоделируйте некоторые данные с помощью одного линейного члена X1 и одного нелинейного члена X2 (в данном случае полином, но это не имеет значения):
library('runjags')
library('mgcv')
set.seed(2018-09-06)
N <- 100
dataset <- data.frame(X1 = runif(N,-1,1), X2 = runif(N,-1,1))
dataset$ll <- with(dataset, 1 + 0.15*X1 + 0.25*X2 - 0.2*X2^2 + 0.15*X2^3 + rnorm(N,0,0.1))
dataset$Y <- rpois(N, exp(dataset$ll))
# Non-linear relationship with log lambda:
with(dataset, plot(X2, ll))
Затем запустите jagam, НО не забудьте исключить член перехвата, указав 0 + в правой части:
# Get the JAGAM stuff excluding intercept:
jd <- jagam(Y ~ 0 + s(X2), data=dataset, file='jagam.txt',
sp.prior="gamma",diagonalize=TRUE,family='poisson')
В качестве альтернативы вы можете оставить здесь термин перехват и удалить его из вашей модели. Это дает нам файл jagam.txt, который выглядит следующим образом:
model {
eta <- X %*% b ## linear predictor
for (i in 1:n) { mu[i] <- exp(eta[i]) } ## expected response
for (i in 1:n) { y[i] ~ dpois(mu[i]) } ## response
## prior for s(X2)...
for (i in 1:8) { b[i] ~ dnorm(0, lambda[1]) }
for (i in 9:9) { b[i] ~ dnorm(0, lambda[2]) }
## smoothing parameter priors CHECK...
for (i in 1:2) {
lambda[i] ~ dgamma(.05,.005)
rho[i] <- log(lambda[i])
}
}
Вы можете удалить первую и последнюю строки, а также две строки, начинающиеся с for (i in 1:n), так как мы будем копировать их сами. Теперь скопируйте все оставшееся содержимое файла и перейдите к вашей (не GAM) модели только с линейными предикторами (и / или случайными эффектами или чем-либо еще) - например:
model <- 'model{
for(i in 1:N){
log(mean[i]) <- intercept + coef*X1[i]
Y[i] ~ dpois(mean[i])
}
# Our priors:
intercept ~ dnorm(0, 10^-6)
coef ~ dnorm(0, 10^-6)
#data# N, X1, Y
#monitor# intercept, coef
}'
Затем вставьте бит GAM, который вы скопировали, в конец, чтобы получить:
model <- 'model{
for(i in 1:N){
log(mean[i]) <- intercept + coef*X1[i] + eta[i]
Y[i] ~ dpois(mean[i])
}
# Our priors:
intercept ~ dnorm(0, 10^-6)
coef ~ dnorm(0, 10^-6)
#data# N, X1, Y, X
#monitor# intercept, coef, b, rho
## JAGAM
eta <- X %*% b ## linear predictor
## prior for s(X2)...
for (i in 1:8) { b[i] ~ dnorm(0, lambda[1]) }
for (i in 9:9) { b[i] ~ dnorm(0, lambda[2]) }
## smoothing parameter priors CHECK...
for (i in 1:2) {
lambda[i] ~ dgamma(.05,.005)
rho[i] <- log(lambda[i])
}
## END JAGAM
}'
Обратите внимание на добавление + eta[i] к строке GLM для учета термина (ов) GAM, а также добавление b и rho к мониторам. Это должно быть все, что вам нужно сделать для модели (кроме проверки априорных параметров на параметры сглаживания и т. Д., Как предлагается).
Затем нам нужно извлечь новую переменную данных X для использования с JAGS:
X <- jd$jags.data$X
Вы также можете извлечь начальные значения для b и лямбда при необходимости. Наконец, мы можем запустить модель, используя runjags:
results <- run.jags(model, n.chains=2, data=dataset)
results
Конечно, этот глупый пример ничего не выигрывает, помещая код jagam в более простую модель - jagam мог бы создать для нас всю модель (включая перехват и линейный предиктор). Но этот подход может иметь значение при добавлении относительно небольшого компонента GAM к более крупной и уже существующей модели, которая была написана для использования некоторых функций в runjags...
Если мы хотим использовать sim2jam, чтобы вернуться назад и использовать соответствующие диагностические / вспомогательные функции из mgcv для встроенного объекта runjags, в настоящее время необходимо вызвать rjags напрямую, чтобы получить больше образцов:
library('rjags')
sam <- jags.samples(as.jags(results), c('b','rho'), n.iter=10000)
jam <- sim2jam(sam,jd$pregam)
plot(jam)
Здесь не хватает двух вещей:
1) Возможность использовать sim2jam без необходимости делать больше сэмплов в rjags. Для этого нужны некоторые дополнения к классу mcarray в пакете rjags, над которым я сейчас работаю.
2) Возможность для template.jags() сделать все это автоматически для вас - это в моем списке вещей, которые нужно реализовать в будущем.
Надеюсь, это поможет - мне было бы интересно услышать, как вы поживаете.
Matt