3

我正在 STAN（rstan 库）中拟合逻辑模型。我的响应变量没有任何缺失值，但是我的协变量“HB”之一是二进制的并且缺少条目。

因此，目标是在每次迭代时使用伯努利先验（参数为 0.5）来估算二进制向量中缺失的条目。

但是，我遇到了问题：

  • 丢失的数据需要在参数或转换参数块中声明为实数或向量；
  • 模型块中伯努利分布的实现需要是整数；
  • 据我所知，STAN 中没有将实数或向量转换为整数的功能。

我使用了STAN 用户指南第 3.3 节中提供的指南。对于下面的模型，解析器在伯努利赋值行（模型块中的倒数第二行）给我一个错误，说它需要整数。注意：我还尝试在参数块中将 HB_miss 定义为实数并得到相同的错误。 

m2 <- '
data {                          
int<lower=0> N;                // total number of observations
int<lower=0,upper=1> y[N];     // setting the dependent variable y as binary
vector[N] X;                   // independent variable 1

int<lower=0> N_obs; 
int<lower=0> N_miss; 
int<lower=1, upper=N> ii_obs[N_obs]; 
int<lower=1, upper=N> ii_miss[N_miss]; 

vector[N_obs] HB_obs;         // independent variable 2 (observed) 

}
parameters {
real b_0;                      // intercept
real b_X;                      // beta 1,2, ...
real b_HB;
vector[N_miss] HB_miss;
}
transformed parameters {
vector[N] HB;
HB[ii_obs] = HB_obs;
HB[ii_miss] = HB_miss;
}
model {
b_0 ~ normal(0,100);           
b_X ~ normal(0,100);           
b_HB ~ normal(0,100); 
HB_miss ~ bernoulli(0.5); // This is where the parser gives me an error
y ~ bernoulli_logit(b_0 + b_X * X + b_HB * HB); // model
}

有什么想法可以在STAN中有效地在HB_miss之前分配伯努利吗？

4

2 回答 2

3

由于您提到的原因，在 Stan 程序中不可能将缺失的离散值视为未知数。Stan 中的所有算法都使用梯度，并且没有为离散未知数定义导数。

相反，您需要边缘化未知值，当一切都是二进制时，这并不太乏味。本质上，您可以使用log_mix 参数为：

  1. 缺失值的概率为 1，您所说的概率为 0.5
  2. 如果缺失值为 1，则相关观察的对数似然贡献
  3. 如果缺失值为 0，则相关观察的对数似然贡献

所以，它会是这样的

for (n in 1:N)
  target += log_mix(0.5, bernoulli_logit_lpmf(y[n] | b_0 + b_X * X[i] + b_HB),
                         bernoulli_logit_lpmf(y[n] | b_0 + b_X * X[i]));

有关更多详细信息，您可以阅读此博客文章。

于 2019-08-15T05:45:32.863 回答
2

感谢上面 Ben 的回答，这里是上述模型的完整解决方案/工作版本（添加了对混合概率的随机影响，而不是原来的 0.5 信念）：

data {                          
  int<lower=0> N;                  // total number of observations
  int<lower=0,upper=1> y[N];       // setting the dependent variable y as binary
  vector[N] X;                     // independent variable 1 (no intercept in the data section)
  int HB[N];                       // dummy coded HB with: '1-2'=0, '3-14'=1, 'Missing'=-1
}
parameters {
  real b_0;                      // intercept
  real b_X;                      // beta 1,2, ...
  real b_HB;
  real<lower=0,upper=1> lambda;  // mixture probability: lambda for HB_miss=1, and (1-lambda) for HB_miss=0 
}
model {
  b_0 ~ normal(0,100);           // priors
  b_X ~ normal(0,100);           
  b_HB ~ normal(0,100); 
  lambda ~ uniform(0,1);

  for (i in 1:N) {
    if (HB[i] == -1) {
      target += log_mix(lambda, bernoulli_logit_lpmf(y[i]| b_0 + b_X * X[i] + b_HB), bernoulli_logit_lpmf(y[i]| b_0 + b_X * X[i]));
    } else {
      HB[i] ~ bernoulli(lambda);
      y[i] ~ bernoulli_logit(b_0 + b_X * X[i] + b_HB * HB[i]); 
    }
  }   
}
'
于 2019-08-15T20:46:47.570 回答