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

Question

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

下面是一个最小的例子，我希望它能说明问题。在 stan 代码中有

有人评论说“模型”部分正在运行，但忽略了多层次方面并平等对待每个较低层次，而不管高层来源如何，因此不会按高层顺序收缩（见图）。
带有 forloop 的“模型”部分，我虽然会做我想做的事，但需要很长时间才能完成，并且有很多警告（Rhat、treedepth、贝叶斯分数、低 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,0,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);
  }

}

mcmc_areas(1-6),(7-12),(12-18)[=13各块无收缩的工作模型=]

Answer 1

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

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

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

multilevel stan model with three hierarchies

bayesian

multi-level

stan