应用贝叶斯模型 (JAGS) 进行各种迭代
Apply a Bayesian model (JAGS) for various iterations
考虑以下数据框：
set.seed(5678)
sub_df<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
假设我想要 运行 以下贝叶斯线性模型，其中 returns samples，一个 mc.array 对象：
library("rjags")
library("coda")
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
# identity
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
##compile
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
##samples returns a list of mcarray objects
samples<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
鉴于 samples$beta1[,,] 表示来自 jags 模型参数后验分布的随机样本，那么总结一下，我的下一步是计算后验分布的均值和 95% 可信区间.所以我会使用：
coeff_output<- round(quantile(samples$beta1[,,],probs=c(0.5,0.025,0.975)),3)
现在，假设我的实际数据框有多个级别 clustersize。
set.seed(5672)
df<- data.frame(clustersize= c(rep(1, 4), rep(2,4), rep(3, 3)),
lepsp= c("A", "B", "C", "D", "B", "C", "D", "E", "A", "D", "F"),
dens= round(runif(11, c(0, 1)), 3),
db= sample(1:10, 11, replace=TRUE))
我如何 运行 这个模型分别用于 clustersize 的每个级别，并使用 forloop 或 apply 函数将输出编译成单个结果数据框？对于 clustersize 的每个级别，生成的 mc.array 对象 samples 应该输出到 result_list 并且 coeff_output 应该输出到数据框 result_coeff.
下面我分别计算每个 clustersize 的输出，以生成预期的结果列表和数据框。
#clustersize==1
sub_df1<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples1<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output1<-
data.frame(as.list(round(quantile(samples1$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
#clustersize==2
sub_df2<- data.frame(clustersize= rep(2,4),
lepsp= c( "B", "C", "D", "E"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples2<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output2<-
data.frame(as.list(round(quantile(samples2$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
#clustersize==3
sub_df3<- data.frame(clustersize= rep(3, 3),
lepsp= c("A", "D", "F"),
dens= round(runif(3, c(0, 1)), 3),
db= sample(1:10, 3, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples3<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output3<-
data.frame(as.list(round(quantile(samples3$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
期望的最终输出：
result_list<- list(samples1, samples2, samples3)
result_coeff<-rbind(coeff_output1, coeff_output2, coeff_output3)
这里是一个link到实际的数据框。该解决方案应该能够处理集群大小高达 600 的大型数据帧。
download.file("https://drive.google.com/file/d/1ZYIQtb_QHbYsInDGkta-5P2EJrFRDf22/view?usp=sharing",temp)
您可以使用 purrr 包中的 map 和 split 不同的 clustersize:
library(rjags)
library(coda)
library(purrr)
set.seed(5678)
set.seed(5672)
df<- data.frame(clustersize= c(rep(1, 4), rep(2,4), rep(3, 3)),
lepsp= c("A", "B", "C", "D", "B", "C", "D", "E", "A", "D", "F"),
dens= round(runif(11, c(0, 1)), 3),
db= sample(1:10, 11, replace=TRUE))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
# identity
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
# split data for different clustersize and calculate result
result <- df %>% split(.$clustersize) %>% map(~{
dataForJags <- list(dens=.x$dens, db=.x$db, N=length(.x$dens))
##compile
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
##samples returns a list of mcarray objects
samples<-jags.samples(model= mod1,variable.names=c("beta1","int","mu","tau"),n.iter=100000)
coeff_output<- data.frame(as.list(round(quantile(samples$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
list(samples = samples, coeff_output = coeff_output)
}) %>% transpose()
result$samples
result$coeff_output
请注意，根据您的要求，使用 purrr::transpose 将最终结果转换为 list 样本和 list 系数。
这里有几个问题需要考虑，这些问题是由您尝试做的事情的规模引起的。您正在创建 550 多个不同的 jags.sample 对象，每个对象迭代 100000 次，然后尝试将所有对象存储在一个列表中。在大多数机器上，这会导致内存问题：输出太大。
至少有两种方法可以解决这个问题：
- 采取措施尽可能减少我们输入数据的内存使用量。
- 调整我们的 JAGS 输出，使其不会从每个链中保存太多迭代。
我已经对您的代码进行了一些修改，以使其能够与您的实际数据集一起使用。
创建输入数据：
在您的原始代码中，clustersize 和 db 都具有数据类型 numeric，即使它们只需要是整数。 numeric 类型占用 8 个字节，而 integer 类型只占用 4 个字节。如果我们将这两列强制转换为 integer 类型，我们实际上可以在下一步中将数据帧列表的内存大小减少大约 30%。
library("tidyverse")
#### Load Raw Data ####
df <- read_csv("example.csv") %>%
select(-1) %>%
mutate(clustersize = as.integer(clustersize),
db = as.integer(db))
初始 JAGS 调整
您对每个链使用的迭代次数过多； niter = 100000 非常高。您还应该使用 n.burn 指定老化期，使用 n.adapt、 指定适应期，并使用 thin 指定细化参数。细化参数在这里尤为重要——这直接减少了我们从每个链中保存的迭代次数。细化参数 50 意味着我们只保存每 50 个结果。
有 post-hoc 方法用于 select 细化参数、老化和适应期，但该讨论超出了 SO 的范围。有关所有这些参数的作用的一些基本信息，这里有一个很好的答案：。目前，我提供的值将允许此代码在您的整个数据集上 运行，但我建议您仔细 select 您用于最终分析的值。
使用 tidybayes
以下解决方案使用 tidybayes 包。这提供了一个干净的输出，并允许我们将所有系数汇总整齐地行绑定到一个数据帧中。请注意，我们使用 coda.samples() 而不是 jags.samples()，因为这提供了一个更通用的 MCMC 对象，我们可以将其传递给 spread_draws()。我们还使用 dplyr::group_split()，它比 split().
的计算效率略高
library("rjags")
library("coda")
library("tidybayes")
set.seed(5672)
result <- df %>% group_split(clustersize) %>% map(~{
dataForJags <- list(dens=.x$dens, db=.x$db, N=length(.x$dens))
# Declare model structure
mod1 <- jags.model(textConnection(model),
data=dataForJags,
n.chains=2)
# samples returns a list of mcmc objects
samples<-coda.samples(model=mod1,
variable.names=c("beta1","int","mu","tau"),
n.burn=10000,
n.adapt=5000,
n.iter=25000,
thin=50
)
# Extract individual draws
samp <- spread_draws(samples, beta1)
# Summarize 95% credible intervals
coeff_output <- spread_draws(samples, beta1) %>%
median_qi(beta1)
list(samples = samp, coeff_output = coeff_output)
}) %>% transpose()
# List of sample objects
result$samples
# Dataframe of coefficient estimates and 95% credible intervals
result_coeff <- bind_rows(result$coeff_output, .id = "clustersize")
考虑以下数据框：
set.seed(5678)
sub_df<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
假设我想要 运行 以下贝叶斯线性模型，其中 returns samples，一个 mc.array 对象：
library("rjags")
library("coda")
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
# identity
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
##compile
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
##samples returns a list of mcarray objects
samples<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
鉴于 samples$beta1[,,] 表示来自 jags 模型参数后验分布的随机样本，那么总结一下，我的下一步是计算后验分布的均值和 95% 可信区间.所以我会使用：
coeff_output<- round(quantile(samples$beta1[,,],probs=c(0.5,0.025,0.975)),3)
现在，假设我的实际数据框有多个级别 clustersize。
set.seed(5672)
df<- data.frame(clustersize= c(rep(1, 4), rep(2,4), rep(3, 3)),
lepsp= c("A", "B", "C", "D", "B", "C", "D", "E", "A", "D", "F"),
dens= round(runif(11, c(0, 1)), 3),
db= sample(1:10, 11, replace=TRUE))
我如何 运行 这个模型分别用于 clustersize 的每个级别，并使用 forloop 或 apply 函数将输出编译成单个结果数据框？对于 clustersize 的每个级别，生成的 mc.array 对象 samples 应该输出到 result_list 并且 coeff_output 应该输出到数据框 result_coeff.
下面我分别计算每个 clustersize 的输出，以生成预期的结果列表和数据框。
#clustersize==1
sub_df1<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples1<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output1<-
data.frame(as.list(round(quantile(samples1$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
#clustersize==2
sub_df2<- data.frame(clustersize= rep(2,4),
lepsp= c( "B", "C", "D", "E"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples2<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output2<-
data.frame(as.list(round(quantile(samples2$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
#clustersize==3
sub_df3<- data.frame(clustersize= rep(3, 3),
lepsp= c("A", "D", "F"),
dens= round(runif(3, c(0, 1)), 3),
db= sample(1:10, 3, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
samples3<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
coeff_output3<-
data.frame(as.list(round(quantile(samples3$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
期望的最终输出：
result_list<- list(samples1, samples2, samples3)
result_coeff<-rbind(coeff_output1, coeff_output2, coeff_output3)
这里是一个link到实际的数据框。该解决方案应该能够处理集群大小高达 600 的大型数据帧。
download.file("https://drive.google.com/file/d/1ZYIQtb_QHbYsInDGkta-5P2EJrFRDf22/view?usp=sharing",temp)
您可以使用 purrr 包中的 map 和 split 不同的 clustersize:
library(rjags)
library(coda)
library(purrr)
set.seed(5678)
set.seed(5672)
df<- data.frame(clustersize= c(rep(1, 4), rep(2,4), rep(3, 3)),
lepsp= c("A", "B", "C", "D", "B", "C", "D", "E", "A", "D", "F"),
dens= round(runif(11, c(0, 1)), 3),
db= sample(1:10, 11, replace=TRUE))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
# identity
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"
# split data for different clustersize and calculate result
result <- df %>% split(.$clustersize) %>% map(~{
dataForJags <- list(dens=.x$dens, db=.x$db, N=length(.x$dens))
##compile
mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)
##samples returns a list of mcarray objects
samples<-jags.samples(model= mod1,variable.names=c("beta1","int","mu","tau"),n.iter=100000)
coeff_output<- data.frame(as.list(round(quantile(samples$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
list(samples = samples, coeff_output = coeff_output)
}) %>% transpose()
result$samples
result$coeff_output
请注意，根据您的要求，使用 purrr::transpose 将最终结果转换为 list 样本和 list 系数。
这里有几个问题需要考虑，这些问题是由您尝试做的事情的规模引起的。您正在创建 550 多个不同的 jags.sample 对象，每个对象迭代 100000 次，然后尝试将所有对象存储在一个列表中。在大多数机器上，这会导致内存问题：输出太大。
至少有两种方法可以解决这个问题：
- 采取措施尽可能减少我们输入数据的内存使用量。
- 调整我们的 JAGS 输出，使其不会从每个链中保存太多迭代。
我已经对您的代码进行了一些修改，以使其能够与您的实际数据集一起使用。
创建输入数据：
在您的原始代码中，clustersize 和 db 都具有数据类型 numeric，即使它们只需要是整数。 numeric 类型占用 8 个字节，而 integer 类型只占用 4 个字节。如果我们将这两列强制转换为 integer 类型，我们实际上可以在下一步中将数据帧列表的内存大小减少大约 30%。
library("tidyverse")
#### Load Raw Data ####
df <- read_csv("example.csv") %>%
select(-1) %>%
mutate(clustersize = as.integer(clustersize),
db = as.integer(db))
初始 JAGS 调整
您对每个链使用的迭代次数过多； niter = 100000 非常高。您还应该使用 n.burn 指定老化期，使用 n.adapt、 指定适应期，并使用 thin 指定细化参数。细化参数在这里尤为重要——这直接减少了我们从每个链中保存的迭代次数。细化参数 50 意味着我们只保存每 50 个结果。
有 post-hoc 方法用于 select 细化参数、老化和适应期，但该讨论超出了 SO 的范围。有关所有这些参数的作用的一些基本信息，这里有一个很好的答案：
使用 tidybayes
以下解决方案使用 tidybayes 包。这提供了一个干净的输出，并允许我们将所有系数汇总整齐地行绑定到一个数据帧中。请注意，我们使用 coda.samples() 而不是 jags.samples()，因为这提供了一个更通用的 MCMC 对象，我们可以将其传递给 spread_draws()。我们还使用 dplyr::group_split()，它比 split().
library("rjags")
library("coda")
library("tidybayes")
set.seed(5672)
result <- df %>% group_split(clustersize) %>% map(~{
dataForJags <- list(dens=.x$dens, db=.x$db, N=length(.x$dens))
# Declare model structure
mod1 <- jags.model(textConnection(model),
data=dataForJags,
n.chains=2)
# samples returns a list of mcmc objects
samples<-coda.samples(model=mod1,
variable.names=c("beta1","int","mu","tau"),
n.burn=10000,
n.adapt=5000,
n.iter=25000,
thin=50
)
# Extract individual draws
samp <- spread_draws(samples, beta1)
# Summarize 95% credible intervals
coeff_output <- spread_draws(samples, beta1) %>%
median_qi(beta1)
list(samples = samp, coeff_output = coeff_output)
}) %>% transpose()
# List of sample objects
result$samples
# Dataframe of coefficient estimates and 95% credible intervals
result_coeff <- bind_rows(result$coeff_output, .id = "clustersize")