Error when fitting the spike and slab priors

Modeling
, ,
1

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

image
image1041×276 11.3 KB

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
2

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
3

@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.

4

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
5

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
6

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

}