时间戳

首页 > 解决方案 > 将 jagam 代码插入 runjags (JAGS) 模型

jags - 将 jagam 代码插入 runjags (JAGS) 模型

问题描述

我一直在尝试将平滑融入我创建的 runjags 模型中，以模拟海鸟洞穴数量和整个岛屿的分布。通过从模型输出中提取计数数据以及 x 和 y 坐标并使用此页面http://www.petrkeil.com/?p=2385上的 JAGAM 教程，我设法生成了一些平滑代码

我想我可以通过将平滑合并到 jags 模型中来提高模型性能，但我不知道如何做到这一点。你能给我一些关于如何实现这一目标的指示吗？

我在下面附上了一段 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

标签： jagsrunjags

解决方案

这是一个非常好的问题，也是使用 jagam 的（非常有用的）输出将 GAM 项添加到模型中的好主意。在您的情况下，我建议使用 jagam 仅生成 GAM 术语而不生成其他任何内容（甚至不是截距），然后将 jagam 模型输出的相关部分复制/粘贴到您现有的模型代码中以及从 jagam 获取 X 数据变量并将其用作您的数据。这最容易用一个例子来演示：

首先用一个线性项 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    
}'

请注意在 GLM 行中添加了 + eta[i] 以考虑 GAM 项，以及在监视器中添加 b 和 rho。这应该是您需要为模型做的所有事情（除了按照建议检查先验的平滑参数等）。

然后我们需要提取新的 X 数据变量以用于 JAGS：

X <- jd$jags.data$X

如果需要，您还可以提取 b 和 lambda 的初始值。最后，我们可以使用 runjags 运行模型：

results <- run.jags(model, n.chains=2, data=dataset)
results

当然，通过将 jagam 代码放入更简单的模型中，这个愚蠢的示例一无所获——jagam 可以为我们创建整个模型（包括截距和线性预测器）。但是，当将一个相对较小的 GAM 组件添加到一个较大且预先存在的模型中时，这种方法可能具有价值，该模型已被编写为使用 runjags 中的某些功能......

如果我们想使用 sim2jam 返回并在拟合的 runjags 对象上使用来自 mgcv 的相关诊断/帮助函数，目前需要直接调用 rjags 以获取更多样本：

library('rjags')
sam <- jags.samples(as.jags(results), c('b','rho'), n.iter=10000)
jam <- sim2jam(sam,jd$pregam)
plot(jam)

这里缺少两件事：

1) 无需在 rjags 中做更多样本即可使用 sim2jam 的能力。这需要对我目前正在处理的 rjags 包中的 mcarray 类进行一些添加。

2) template.jags() 自动为您完成所有这些工作的能力——这是我将来要实现的事情的清单。

希望对您有所帮助-我很想听听您的进展情况。

马特

推荐阅读