具有三个层次结构的多级 stan 模型

如何解决具有三个层次结构的多级 stan 模型

假设我有一个多级数据结构。使用全局分布，我从中绘制高级分布，从中绘制低级分布，从中绘制响应变量。我将如何在 stan 模型中实现这样的事情。

下面是一个最小的例子，我希望可以说明问题。在斯坦代码中有

• 一个评论过的“模型”部分正在工作，但忽略了多层次方面并将每个较低级别视为平等，无论高级起源如何，并为此提供不按高级顺序收缩（参见图片）。
• 一个带有 forloop 的“模型”部分，我虽然会做我想做的事，但需要很长时间才能完成，并且有很多警告（Rhat、树深度、贝叶斯分数、低 ESS）

我对建模非常缺乏经验，所有关于 ML 建模的教程都没有循环方法我虽然在这里很有意义，所以我怀疑我完全朝着错误的方向前进。因此，我们将不胜感激任何帮助。

用于生成和运行模型的 R 代码

```{r}
library(ggplot2)
library(rstan)
library(bayesplot)
library(loo)


set.seed(0)
numberOFTest <- 50
resposeList <- c()
highLevelIndexList <- c()
lowLevelIndexList  <- c()
for (highLevelIndex in seq(3)) {
  highLevelMu <- rnorm(1,5) #High Level avarage drawn drom normal dist
  for (lowLevelIndex in seq(6)) {
    lowLevelMu <- rnorm(1,highLevelMu,3) #low Level avarage drawn drom normal dist
    resposeList <- c(resposeList,rnorm(numberOFTest,lowLevelMu,1))
    highLevelIndexList <- c(highLevelIndexList,rep(highLevelIndex,numberOFTest) )
    lowLevelIndexList <- c(lowLevelIndexList,rep(paste0(highLevelIndex,lowLevelIndex),numberOFTest) )
  
  }
}
lowLevelIndexList <- as.integer( as.factor(lowLevelIndexList))

dat <- list(Nr=length(resposeList),Nh=length(unique(highLevelIndexList)),Nl=length(unique(lowLevelIndexList)),#############
              r=seq(length(resposeList)),response=resposeList,highLevelIndex=(highLevelIndexList),lowLevelIndex=(lowLevelIndexList)
              )


model <- stan(data=dat,file="modle.stan",chains = 4,cores = 8,control=list(adapt_delta=0.95),iter = 2000)

标准代码

data {
  int<lower=1> Nr; // number of rows (entries)
  int<lower=1> Nh; // number of highlevels 
  int<lower=1> Nl; // number of low levels 
  int<lower=1> r[Nr]; // rownumber
  real response[Nr]; // 
  int<lower=1> highLevelIndex[Nr]; // 
  int<lower=1> lowLevelIndex[Nr]; // 
}


parameters {
  real<lower=0> sigma; //
  real<lower=0> sigma_global; //
  real<lower=0> sigma_hL; //
  real<lower=0> sigma_lL; //
  real a_global; //
  real a_hL[Nh]; //
  real a_lL[Nl]; //
}

transformed parameters {
  vector[Nr] mu; 
  for (rowIndex in 1:Nr) { //
    mu[rowIndex] = a_lL[lowLevelIndex[rowIndex]];
  }
}


//Idea of what I want

model {
  sigma_global ~ exponential(0.1);  //hyper
  sigma_hL ~ exponential(0.1);  //hyper
  sigma_lL ~ exponential(0.1);  //hyper
  a_global ~ normal(0,sigma_global);
  //-----------
  for (rowIndex in 1:Nr) { //
    a_hL[highLevelIndex[rowIndex]] ~   normal(a_global,sigma_hL);  //hyper
    a_lL[lowLevelIndex[rowIndex]] ~   normal(a_hL[highLevelIndex[rowIndex]],sigma_lL);  
  }
    ////
  sigma ~ exponential(0.1);  
  response ~ normal(mu,sigma);
}


/*
model {///// Working but not what I want
  a_lL ~   normal(0,10); 
  ////
  sigma ~ exponential(0.1);  
  response ~ normal(mu,sigma);
}
*/

generated quantities { // for the PPC and loo
  real n_rep[Nr];
  vector[Nr] log_lik;
  
  for (i in 1:Nr){
    n_rep[i] = normal_rng(mu[i],sigma);
    log_lik[i] = normal_lpdf(response[i] | mu[i],sigma);
  }

}

每块(1-6)、(7-12)和(12-18)内无收缩的工作模型的mcmc_areas

解决方法

发现错误：我需要通过查找表将低级值映射到高级值。下面是一个工作版本，它也只需要一秒钟即可完成。

data {
  int<lower=1> Nr; // number of rows (entries)
  int<lower=1> Nh; // number of highlevels 
  int<lower=1> Nl; // number of low levels 
  real response[Nr]; // 
  int<lower=1> highLevelIndex[Nr]; // 
  int<lower=1> lowLevelIndex[Nr]; // 
  int<lower=1> lookUp[Nl];//lookuptable which highlevel value a lowlevel value fits to
}

parameters {
  real<lower=0> sigma; //
  real<lower=0> sigma_hL; //
  real<lower=0> sigma_lL; //
  real a_hL[Nh]; //
  real a_lL[Nl]; //
}

transformed parameters {
  vector[Nr] mu; 
  for (rowIndex in 1:Nr) { //
    mu[rowIndex] = a_lL[lowLevelIndex[rowIndex]];
  }
}

model {
  sigma_hL ~ exponential(0.01);
  a_hL ~          normal(0,sigma_hL); 
  sigma_lL ~ exponential(0.01);
  for (lLIdx in 1:Nl) { //
    a_lL[lLIdx] ~  normal(a_hL[lookUp[lLIdx]],sigma_lL);  
  }
  ////
  sigma ~ exponential(0.1);  
  response ~ normal(mu,sigma);
}

generated quantities { // for the PPC and loo
  real n_rep[Nr];
  vector[Nr] log_lik;
  
  for (i in 1:Nr){
    n_rep[i] = normal_rng(mu[i],sigma);
    log_lik[i] = normal_lpdf(response[i] | mu[i],sigma);
  }
}