前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >生态学JAGS模拟数据、回归、CORMACK-JOLLY-SEBER (CJS) 模型拟合MCMC 估计动物存活率

生态学JAGS模拟数据、回归、CORMACK-JOLLY-SEBER (CJS) 模型拟合MCMC 估计动物存活率

作者头像
拓端
发布2021-12-16 12:34:57
6150
发布2021-12-16 12:34:57
举报
文章被收录于专栏:拓端tecdat拓端tecdat

原文链接:http://tecdat.cn/?p=24721

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

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

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

线性回归示例

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

代码语言:javascript
复制
library(R2jags)

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

代码语言:javascript
复制
data{
# 似然函数:
for (i in 1:N){
y\[i\] ~  # tau是精度(1/方差)。

}

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

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

代码语言:javascript
复制
# 模拟的参数 
N  # 样本
x <- 1:N # 预测因子
alpha # 截距
beta  # 斜率
sigma# 残差sd
 1/(sigma*sigma) # 精度
# 在模拟步骤中,参数被当作数据处理

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

代码语言:javascript
复制
# 运行结果
out

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

代码语言:javascript
复制
# 重新格式化输出
mcmc(out)
代码语言:javascript
复制
dim
代码语言:javascript
复制
dat

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

代码语言:javascript
复制
# 用BUGS语言指定模型
model <-     

for (i in 1:N){
y\[i\] ~ dnorm(mu\[i\], tau) # tau是精度(1/方差)


alpha  截距
beta # 斜率
sigma  # 标准差


# 数据
dta <- list(y = dt, N = length(at), x = x)

# 初始值
inits 


# MCMC设置
ni <- 10000


# 从R中调用JAGS
jags()

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

代码语言:javascript
复制
# 总结后验
print(res)

检查收敛:

代码语言:javascript
复制
# 追踪图
plot(res)

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

代码语言:javascript
复制
# 后验分布
plot(res)

模拟示例

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

让我们模拟一下!

代码语言:javascript
复制
# 恒定的生存和重新捕获概率
for (i in 1:nd){
   for (t in f:(on-1)){

#概率
for (i in 1:nid){
   # 定义潜伏状态和第一次捕获时的观察值
   z\[i,f\[i\] <- 1
   mu2\[i,1\] <- 1 * z\[i,f\[i\] # 在第一次捕获时检测为1("以第一次捕获为条件")。
   # 然后处理以后的情况
   for (t in (f\[i\]+1):non){
      # 状态进程
      mu1\[i,t\] <- phi * z 
      # 观察过程
      mu2\[i,t\] <- p * z

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

代码语言:javascript
复制
# 用于模拟的参数 
n = 100 # 个体的数量
meanhi <- 0.8 # 存活率
meap <- 0.6 # 重捕率
data<-list

现在运行 JAGS

代码语言:javascript
复制
out

格式化输出:

代码语言:javascript
复制
as.mcmc(out)
代码语言:javascript
复制
head(dat)

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

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

代码语言:javascript
复制
# 倾向性和约束
for (i in 1:nd){
   for (t in f\[i\]:(nn-1)){


mehi ~ dunif(0, 1) # 平均生存率的先验值
Me ~ dunif(0, 1) # 平均重捕的先验值
# 概率
for (i in 1:nd){
   # 定义第一次捕获时的潜伏状态
   z\[i\]\] <- 1
   for (t in (f\[i\]+1):nions){
      # 状态过程
      z\[i,t\] ~ dbern(mu1\[i,t\])
      # 观察过程
      y\[i,t\] ~ dbern(mu2\[i,t\])

准备数据:

代码语言:javascript
复制
# 标记的场合的向量
gerst <- function(x) min(which(x!=0))
# 数据
jagta
代码语言:javascript
复制
# 初始值
for (i in 1:dim\]){
       min(which(ch\[i,\]==1))
      max(which(ch\[i,\]==1))

function(){list(meaphi, mep , z ) }

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

标准 MCMC 设置:

代码语言:javascript
复制
ni <- 10000

准备运行 JAGS

代码语言:javascript
复制
# 从R中调用JAGS
jags(nin = nb, woy = getwd() )

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

代码语言:javascript
复制
print(cj3)

非常接近!

跟踪图

代码语言:javascript
复制
trplot

后验分布图

代码语言:javascript
复制
denplot
本文摘选《R语言生态学JAGS模拟数据、线性回归、CORMACK-JOLLY-SEBER (CJS) 模型拟合MCMC 估计动物存活率和可视化》
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-12-15,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 拓端数据部落 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 线性回归示例
  • 模拟示例
  • 标准 MCMC 设置:
  • 跟踪图
  • 后验分布图
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档