Bernoulli Prior in R STAN

Question

I am fitting a logistic model in STAN (rstan library). My response variable does not have any missing values, however one of my covariates "HB" is binary and has missing entries.

Thus the goal is to impute at each iteration the missing entries in the binary vector, using a Bernoulli prior (with parameter, say, 0.5).

However, I am running into issues:

The missing data needs to be declared either as real or vector in the parameters or transformed parameters block;
Realizations from a Bernoulli distribution in the model block need to be integers;
As far as I am aware there is no function in STAN to convert a real or vector to integer.

I used the guidelines provided in section 3.3 of the STAN user guide. With the model below, the parser gives me an error at the bernoulli assignment line (penultimate line in the model block), saying it needs integers. Note: I also tried defining HB_miss as a real in the parameters block and getting the same error.

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
}

Any ideas how I can assign a bernoulli prior to HB_miss effectively in STAN?

score 3 · Accepted Answer · answered Aug 15 '19 at 05:45

For the reasons you mentioned, it is not possible to treat missing discrete values as unknowns in a Stan program. All of the algorithms in Stan utilize gradients, and derivatives are not defined for discrete unknowns.

Instead, you need to marginalize over the unknown values, which is not too tedious when everything is binary. Essentially, you can use the log_mix function whose arguments are:

The probability the missing value is 1, which you say is 0.5 in your case
The log-likelihood contribution for the observation in question if the missing value were 1
The log-likelihood contribution for the observation in question if the missing value were 0

So, it would be something like

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

For more details, you could read this blog post.

Got it - it works! Will post full working stan code in case it helps others. — Neodyme, Aug 15 '19 at 11:59

score 2 · Answer 2 · edited Dec 02 '20 at 00:07

Thanks to Ben's answer above, here is the full solution / working version of the model above (added a random effect on the mixture probability instead of the original 0.5 belief):

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

not desperately important, but could you indent this code ... ? — Ben Bolker, Aug 15 '19 at 20:49

Bernoulli Prior in R STAN

2 Answers2