R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化
最近我们被客户要求撰写关于生态学的研究报告,包括一些图形和统计输出。
本文,我通过两个种群生态学家可能感兴趣的例子来说明使用“JAGS”来模拟数据:首先是线性回归,其次是估计动物存活率(公式化为状态空间模型)。
最近,我一直在努力模拟来自复杂分层模型的数据。我现在正在使用 JAGS
。
模拟数据 JAGS
很方便,因为你可以使用(几乎)相同的代码进行模拟和推理,并且你可以在相同的环境(即JAGS
)中进行模拟研究(偏差、精度、区间覆盖 )。
线性回归示例
我们首先加载本教程所需的包:
library(R2jags)
然后直接切入正题,让我们从线性回归模型生成数据。使用一个 data
块,并将参数作为数据传递。
data{
# 似然函数:
for (i in 1:N){
y[i] ~ # tau是精度(1/方差)。}
这里, alpha
和 beta
是截距和斜率、 tau
方差的精度或倒数、 y
因变量和 x
解释变量。
我们为用作数据的模型参数选择一些值:
# 模拟的参数
N # 样本
x <- 1:N # 预测因子
alpha # 截距
beta # 斜率
sigma# 残差sd1/(sigma*sigma) # 精度
# 在模拟步骤中,参数被当作数据处理
现在运行 JAGS
; 请注意,我们监控因变量而不是参数,就像我们在进行标准推理时所做的那样:
# 运行结果
out
输出有点乱,需要适当格式化:
# 重新格式化输出
mcmc(out)
dim
dat
现在让我们将我们用来模拟的模型拟合到我们刚刚生成的数据中。不再赘述,假设读者熟悉 JAGS
线性回归。
# 用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()
让我们看看结果并与我们用来模拟数据的参数进行比较(见上文):
# 总结后验
print(res)
检查收敛:
# 追踪图
plot(res)
绘制回归参数和残差标准差的后验分布:
# 后验分布
plot(res)
模拟示例
我现在说明如何使用 JAGS
来模拟来自具有恒定生存和重新捕获概率的模型的数据。我假设读者熟悉这个模型及其作为状态空间模型的公式。
让我们模拟一下!
# 恒定的生存和重新捕获概率
for (i in 1:nd){for (t in f:(on-1)){#概率
for (i in 1:nid){# 定义潜伏状态和第一次捕获时的观察值z[i,f[i] <- 1mu2[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
让我们为参数选择一些值并将它们存储在数据列表中:
# 用于模拟的参数
n = 100 # 个体的数量
meanhi <- 0.8 # 存活率
meap <- 0.6 # 重捕率
data<-list
现在运行 JAGS
:
out
格式化输出:
as.mcmc(out)
head(dat)
我只监测了检测和非检测,但也可以获得状态的模拟值,即个人在每种情况下是生是死。你只需要修改对JAGS
的调用 monitor=c("y","x")
并相应地修改输出。
现在我们将 Cormack-Jolly-Seber (CJS) 模型拟合到我们刚刚模拟的数据中,假设参数不变:
# 倾向性和约束
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]] <- 1for (t in (f[i]+1):nions){# 状态过程z[i,t] ~ dbern(mu1[i,t])# 观察过程y[i,t] ~ dbern(mu2[i,t])
准备数据:
# 标记的场合的向量
gerst <- function(x) min(which(x!=0))
# 数据
jagta
# 初始值
for (i in 1:dim]){min(which(ch[i,]==1))max(which(ch[i,]==1))function(){list(meaphi, mep , z ) }
我们想对生存和重新捕获的概率进行推断:
标准 MCMC 设置:
ni <- 10000
准备运行 JAGS
!
# 从R中调用JAGS
jags(nin = nb, woy = getwd() )
总结后验并与我们用来模拟数据的值进行比较:
print(cj3)
非常接近!
跟踪图
trplot
后验分布图
denplot
R语言生态学JAGS模拟数据、线性回归、Cormack-Jolly-Seber (CJS) 模型拟合MCMC 估计动物存活率和可视化相关推荐
- 【视频】线性回归中的贝叶斯推断与R语言预测工人工资数据|数据分享
最近我们被客户要求撰写关于线性回归的研究报告,包括一些图形和统计输出. 在这个视频中,我们转向简单线性回归中的贝叶斯推断. 我们将使用一个参照先验分布,它提供了频率主义解决方案和贝叶斯答案之间的联系. ...
- R语言淮河流域水库水质数据相关性分析、地理可视化、广义相加模型GAM调查报告...
采样地点:淮河流域一带,昭平台水库.白龟山水库.燕山水库.石漫滩水库.板桥水库.宿鸭湖水库.博山水库.南湾水库.石山口水库.五岳水库.泼河水库.鲶鱼山水库(点击文末"阅读原文"获取 ...
- 价值1143元的《R语言统计分析微生物组数据(Statistical Analysis of Microbiome Data with R)》系列图书
文章目录 <R语言统计分析微生物组数据> 本书简介 作者简介 章节简介 猜你喜欢 写在后面 <R语言统计分析微生物组数据> Statistical Analysis of Mi ...
- R语言第六讲 数据的统计分析
基本命令练习 下面的代码涵盖了一些分析数据常用的一些R语言的命令: #基本向量.矩阵的一般操作 x <- c(1,3,2,5) x x = c(1,6,2) x y = c(1,4,3) len ...
- 《R语言统计分析微生物组数据》图书简介
<R语言统计分析微生物组数据> Statistical Analysis of Microbiome Data with R 京东上原版图书售价1143元 https://item.jd. ...
- R语言data.table导入数据实战:data.table中编写函数并使用SD数据对象
R语言data.table导入数据实战:data.table中编写函数并使用SD数据对象 目录 R语言data.table导入数据实战:data.table中编写函数并使用SD数据对象 #data.t ...
- R语言data.table导入数据实战:把data.frame数据转化为data.table数据
R语言data.table导入数据实战:把data.frame数据转化为data.table数据 目录 R语言data.table导入数据实战:把data.frame数据转化为data.table数据 ...
- R语言data.table导入数据实战:data.table使用by函数进行数据分组(aggregate)
R语言data.table导入数据实战:data.table使用by函数进行数据分组(aggregate) 目录 R语言data.table导入数据实战:data.table使用by函数进行数据分组( ...
- R语言dplyr包通过数据列的索引重命名数据列实战(Rename Column by Index Position)
R语言dplyr包通过数据列的索引重命名数据列实战(Rename Column by Index Position) 目录 R语言dplyr包通过数据列的索引重命名数据列实战(Rename Colum ...
最新文章
- 连连看html游戏全代码js、jquery操作
- c语言utc时间转换北京时间_C/C++标准库之转换UTC时间到local本地时间详解
- matlab如何配置weka,matlab调用weka
- 对于高并发短连接造成Cannot assign requested address解决方法
- python核心编程——python对象
- oracle12c之 控制pdb中sga 与 pga 内存使用
- python编程实现撤销上一步操作_78行Python代码实现现微信撤回消息功能
- 设置mysql8的root可以远程访问
- 技术文档(3)--查看和修改Linux服务器的时区和时间
- Python 学习记录(1)对象命名导致的问题
- Mybatis的两种分页方式:RowBounds和PageHelper
- DOTween的Sequence图例说明
- 游戏开发之C++引用(C++基础)
- Java使用qq邮箱发送email
- 法向量与切向量的转化
- 微型计算机的地址加法器,地址加法器
- C#对Redis的读写与使用
- openbravo erp介绍(一)
- Ubuntu系统如何屏幕截图
- 【20保研】热忱欢迎全国2020届优秀本科毕业生免试攻读重庆大学研究生
热门文章
- 《Adobe Illustrator CC经典教程》—第0课0.19节使用图像描摹
- pycharm 激活码及使用方式
- 美国登月技术退步了?50年前就能载人着陆,怎么现在只能带着史努比绕一圈...
- Nevercenter CameraBag Pro照片滤镜软件 v2023.2.0
- python循环体结束标志_python如何结束循环
- 在asp.net中如何播放MP3的歌曲?
- Spket 破解方法很好的javascript脚本编译器
- 怎样打开微信定位服务器地址,微信位置服务功能,能知道对方位置,你们知道怎么用吗!...
- 移动端H5下载后端文件
- C 宏定义及函数宏定义