0
点赞
收藏
分享

微信扫一扫

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化

老榆 2022-11-10 阅读 53

本文,我通过两个种群生态学家可能感兴趣的例子来说明使用“JAGS”来模拟数据:首先是线性回归,其次是估计动物存活率(公式化为状态空间模型)。

最近,我一直在努力模拟来自复杂分层模型的数据。我现在正在使用 ​​JAGS​​。

模拟数据 ​​JAGS​​​ 很方便,因为你可以使用(几乎)相同的代码进行模拟和推理,并且你可以在相同的环境(即​​JAGS​​)中进行模拟研究(偏差、精度、区间覆盖 )。

线性回归示例

我们首先加载本教程所需的包:

  1.  library(R2jags)
  2.   
  3.   

然后直接切入正题,让我们从线性回归模型生成数据。使用一个 ​​data​​ 块,并将参数作为数据传递。


1.  data{
2. # 似然函数:
3. for (i in 1:N){
4. y[i] ~ # tau是精度(1/方差)。
5.
6. }
7.

这里, ​​alpha​​​ 和 ​​beta​​​ 是截距和斜率、 ​​tau​​​ 方差的精度或倒数、 ​​y​​​ 因变量和 ​​x​​ 解释变量。

我们为用作数据的模型参数选择一些值:


1.   
2. # 模拟的参数
3. N # 样本
4. x <- 1:N # 预测因子
5. alpha # 截距
6. beta # 斜率
7. sigma# 残差sd
8. 1/(sigma*sigma) # 精度
9. # 在模拟步骤中,参数被当作数据处理
10.

现在运行 ​​JAGS​​; 请注意,我们监控因变量而不是参数,就像我们在进行标准推理时所做的那样:

  1.  # 运行结果
  2.  out

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_数据

输出有点乱,需要适当格式化:

  1.  # 重新格式化输出
  2.  mcmc(out)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_02

dim

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_03

dat

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_数据_04

现在让我们将我们用来模拟的模型拟合到我们刚刚生成的数据中。不再赘述,假设读者熟悉 ​​JAGS​​ 线性回归。


1.   
2. # 用BUGS语言指定模型
3. model <-
4.
5. for (i in 1:N){
6. y[i] ~ dnorm(mu[i], tau) # tau是精度(1/方差)
7.
8.
9. alpha 截距
10. beta # 斜率
11. sigma # 标准差
12.
13.
14. # 数据
15. dta <- list(y = dt, N = length(at), x = x)
16.
17. # 初始值
18. inits
19.
20.
21. # MCMC设置
22. ni <- 10000
23.
24.
25. # 从R中调用JAGS
26. jags()
27.

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_数据_05

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_06

让我们看看结果并与我们用来模拟数据的参数进行比较(见上文):

  1.  # 总结后验
  2.  print(res)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_07

检查收敛:

  1.  # 追踪图
  2.  plot(res)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_数据_08

绘制回归参数和残差标准差的后验分布:

  1.  # 后验分布
  2.  plot(res)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_09

模拟示例

我现在说明如何使用 ​​JAGS​​ 来模拟来自具有恒定生存和重新捕获概率的模型的数据。我假设读者熟悉这个模型及其作为状态空间模型的公式。

让我们模拟一下!


1.   
2. # 恒定的生存和重新捕获概率
3. for (i in 1:nd){
4. for (t in f:(on-1)){
5.
6. #概率
7. for (i in 1:nid){
8. # 定义潜伏状态和第一次捕获时的观察值
9. z[i,f[i] <- 1
10. mu2[i,1] <- 1 * z[i,f[i] # 在第一次捕获时检测为1("以第一次捕获为条件")。
11. # 然后处理以后的情况
12. for (t in (f[i]+1):non){
13. # 状态进程
14. mu1[i,t] <- phi * z
15. # 观察过程
16. mu2[i,t] <- p * z
17.
18.

  1.  

让我们为参数选择一些值并将它们存储在数据列表中:


1.   
2. # 用于模拟的参数
3. n = 100 # 个体的数量
4. meanhi <- 0.8 # 存活率
5. meap <- 0.6 # 重捕率
6. data<-list

现在运行 ​​JAGS​​:

out

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_10

格式化输出:

as.mcmc(out)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_11

head(dat)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_12

我只监测了检测和非检测,但也可以获得状态的模拟值,即个人在每种情况下是生是死。你只需要修改对​​JAGS​​​ 的调用 ​​monitor=c("y","x")​​ 并相应地修改输出。

现在我们将 Cormack-Jolly-Seber (CJS) 模型拟合到我们刚刚模拟的数据中,假设参数不变:


1.   
2. # 倾向性和约束
3. for (i in 1:nd){
4. for (t in f[i]:(nn-1)){
5.
6.
7. mehi ~ dunif(0, 1) # 平均生存率的先验值
8. Me ~ dunif(0, 1) # 平均重捕的先验值
9. # 概率
10. for (i in 1:nd){
11. # 定义第一次捕获时的潜伏状态
12. z[i]] <- 1
13. for (t in (f[i]+1):nions){
14. # 状态过程
15. z[i,t] ~ dbern(mu1[i,t])
16. # 观察过程
17. y[i,t] ~ dbern(mu2[i,t])

准备数据:


1.   
2. # 标记的场合的向量
3. gerst <- function(x) min(which(x!=0))
4. # 数据
5. jagta
6.
7. # 初始值
8. for (i in 1:dim]){
9. min(which(ch[i,]==1))
10. max(which(ch[i,]==1))
11.
12. function(){list(meaphi, mep , z ) }
13.

我们想对生存和重新捕获的概率进行推断:

标准 MCMC 设置:

ni <- 10000

准备运行 ​​JAGS​​!

  1.  # 从R中调用JAGS
  2.  jags(nin = nb, woy = getwd() )

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_13

总结后验并与我们用来模拟数据的值进行比较:

print(cj3)

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_14

非常接近!

跟踪图

trplot

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_线性回归_15

后验分布图

denplot

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_r语言_16

R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化_数据_17



举报

相关推荐

0 条评论