5

次のデータ フレームを検討してください。

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))

オブジェクトを返す次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[,,]要約すると、私の次のステップは、事後分布の平均と 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を使用して出力を単一の結果データ フレームにコンパイルするにはどうすればよいでしょうか? の各レベルで、結果のオブジェクトを に出力し、 をデータ フレームに出力する必要があります。forloopapplyclustersizemc.arraysamplesresult_listcoeff_outputresult_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)

これが実際のデータフレームへのリンクです。このソリューションは、クラスターサイズが最大 600 の大きなデータフレームを処理できる必要があります。

download.file("https://drive.google.com/file/d/1ZYIQtb_QHbYsInDGkta-5P2EJrFRDf22/view?usp=sharing",temp)
4

2 に答える 2