我使用以下代码来估计分层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=264713(总观察数(，L=10000(组数(，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
)
如果能得到专家的建议，我将不胜感激。
谢谢。
第一件事是，在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例如。
此外，对于关于斯坦的问题，斯坦的论述非常好——你可能会做其他我没有想到的事情