提高层次 logit/probit 模型的 Stan 程序的速度

Question

我使用下面的代码来估计一个分层的 logit 模型。

data {
  int<lower=1> D;   // number of covariates
  int<lower=0> N;   // number of observations
  int<lower=1> L;   // number of groups (levels)
  int<lower=0,upper=1> y[N];  // binary response
  int<lower=1,upper=L> ll[N];    // group ids
  row_vector[D] x[N];   //covariates
}
parameters {
  real mu[D];
  real<lower=0> sigma[D];
  vector[D] beta[L];
}
model {
  
  mu ~ normal(0, 100);  // prior for mu
//  sigma ~ uniform(0, 100) ; // prior for sigma
  for (l in 1:L) beta[l] ~ normal(mu, sigma);
  
  {
    vector[N] x_beta_ll;
    for (n in 1:N) x_beta_ll[n] = x[n] * beta[ll[n]];
    y ~ bernoulli_logit(x_beta_ll);
  }
}

对于我的数据，N = 264,713（总观察数），L = 10,000（组数），D = 26（协变量数）。然后，我有这个消息，需要很长时间才能采样：

Chain 1: Gradient evaluation took 0.131 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1310 seconds.
Chain 1: Adjust your expectations accordingly!

不知道有没有什么提高速度的好方法。顺便说一句，这是我使用的 R 代码：

fit <- stan(
  file = "hierarchical_logit.stan",  # Stan program
  data = data,    # named list of data
  chains = 1,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 5,              # number of cores (could use one per chain)
  verbose = TRUE
)

我会感谢专家的提示。

谢谢。

Answer 1

首先，在 logit 转换的上下文中，您的先验看起来可能非常糟糕。这取决于您提供的变量是什么，但我的猜测是，它将把这些变量的几乎所有密度都放在 0 和 1。对于连续变量，如果您在一开始就提供标准化变量，这些将意味着概率完全没有空间从 0 步长到 1，同样，对于截距项，大部分密度也将在 0 和 1 处。值得考虑这些 - 并尝试模拟它，例如在 R 中，任何截距项的先验都是hist(boot::inv.logit(rnorm(1000, 0, 100)), breaks = 100))。例如，在截距项上更健康、信息量不足的先验将类似于normal(0, 1.5)，它在概率尺度上大致均匀地分布在 [0,1] 之间（如 McElreath 所建议的那样）。类似的思考和模拟将帮助您找出其他变量的先验（我不知道它们是什么）。提供一些常识性先验可能会很快，具体取决于您的变量/数据是什么样的，并且无论如何都很重要。

除非有特殊原因不这样做，否则我要做的第一件事就是使用 brms 或 rstanarm 运行模型。这是一个标准模型，它们将自动运行一个精心设计的模型，您可以只使用 lmer 类型的语法。如果这仍然非常缓慢，那么您将需要做一些额外的工作。

如果您想自己做（/brms 太慢了），请确保您传入标准化变量，或者在您的模型中标准化它们 (https://mc-stan.org/docs/2_25/stan-users-guide/standardizing-predictors-and-outputs.html)。

目前你有一个居中的参数化 - 一个非居中的参数化可以更快，所以如果它很慢，那么值得尝试一下（在 Stan 手册中有一节）。它看起来像这样：

parameters {
  real mu[D];
  real<lower=0> sigma[D];
  vector[D] beta_z[L];
}
transformed parameters {
  vector[D] beta[L];
  for(l in 1:L){
    beta[l] = mu + sigma * beta_z[l];
  }
}
model {
...as before, but replace beta[l] ~ normal(mu, sigma) with...
for (l in 1:L) beta_z[l] ~ normal(0, 1);
}

最后，您可以编写一个使用链内并行化来实现主要加速的模型。然而，这需要付出一些努力，除非您发现更简单的步骤仍然太慢，否则我不会这样做。请参阅 https://mc-stan.org/users/documentation/case-studies/reduce_sum_tutorial.html 以获取工作示例。

另外，关于 Stan 的问题，Stan discourse 非常好 - 你可能会做其他我没有想到的事情..