Error when fitting the spike and slab priors

Hi everyone,
I am trying to fit a spike and slab regularization prior for logistic regression model. For that I am going to use following mixture of priors.
image


My code is as follows:

data {
  int<lower=1> N;
  int<lower=1> K1; 
  int<lower=0,upper=1> y1[N];
  matrix[N,K1] x1;
  
}

parameters {
   real alpha;
   real<lower=0,upper=1> theta;
  vector<lower=0>[K1] tauj_sqr;
  int <lower=0,upper=1> gammaj[K1] ;
    vector[K1] spike_j;
    vector[K1] slab_j;
      }

transformed parameters {
  vector[K1] beta_j = gammaj .* slab_j + (1-gammaj) .* spike_j;
}

model {
     
    slab_j ~ normal(0, sqrt(tauj_sqr)) ;
    spike_j ~ normal(0, sqrt(0.001)) ;
    gammaj ~ bernoulli(theta);
   theta ~ beta(0.5,0.5);
   tauj_sqr ~ inv_gamma(0.5,0.5);
  
     alpha ~ normal(0, 5);
  
  y1 ~ bernoulli_logit_glm(x1, alpha, beta_j);
  
}


However, I am getting an error regarding to the following line of code :
int <lower=0,upper=1> gammaj[K1] ;

Can anyone guide me how to fix this ?

Also will there be any efficient parameterization that can be used for spike and slab priors ?

Thank you

1 Like

Stan doesn’t support discrete parameters. It’s possible to marginalize out discrete parameters; however, it’s generally recommended to use regularized horseshoe priors instead. The easiest way to do this would probably be to fit the model in brms instead of Stan directly.

I would also recommend taking a look at the Prior Choice recommendations for some of your other priors.

3 Likes

@Christopher-Peterson Thank you for your reply. I will go through the resources you provided. They will be very useful. Unfortunately I cannot use any package to fulfill my task as fitting a univariate logistic model is not my final task. So I may need to use Stan to code these priors. Also I am comparing different regularization priors including spike and slab prior and horseshoe.

The brms package is quite flexible in the models it can fit, but if it doesn’t work for you, you can also use the stan code it creates as a starting point for a more custom model.

1 Like

You can use a marginalized approached (treat the spike and slab a mixture distribution, and marginalize the mixture). You can use a horseshoe (or its variants). Or you can use an approximate spike and slab, which swaps out the point mass at zero with a very small scale normal. The latter two will be your best bet. The first option still requires a point-mass, which is not continuous; it can be ‘simulated’ using a gigantic increment, but it’s going to be less stable than just approximating the spike or using another regularization approach.

1 Like

@Stephen_Martin Thank you for the response. Did you see any example code anywhere that I can follow? I tried to fit a spike and slab prior like follows. But I believe this is not the right way of doing this.

data {
int<lower=1> N;
int<lower=1> K1;
int<lower=0,upper=1> y1[N];
matrix[N,K1] x1;

}

parameters {
real alpha;
real<lower=0,upper=1> theta;
vector<lower=0,upper=1>[K1] gammaj;
vector[K1] spike_j;
vector<lower=0>[K1] xa;
vector<lower=0>[K1] xb;

}

transformed parameters {
vector[K1] slab_j = xa .* sqrt(xb);
vector[K1] beta_j = gammaj .* slab_j + (1-gammaj) .* spike_j;
}

model {

xa ~ std_normal();

xb ~ inv_gamma(0.5, 0.5);
spike_j ~ normal(0, sqrt(0.001)) ;
gammaj ~ beta(1,1);

alpha ~ normal(0, 5);
y1 ~ bernoulli_logit_glm(x1, alpha, beta_j);

}