Spike and Slab redux

Modeling
1

Dear all,

I’m trying to set up a spike and slab prior in Stan. First off, it’s unclear given the back and forth on the forum as to how to set this up. In any case, this is what I tried

modelString = "
data {
int<lower=0> n;
vector [n] readscore;
vector [n] Female; vector [n] ESCS;
vector [n] METASUM; vector [n] PERFEED;
vector [n] JOYREAD; vector [n] MASTGOAL;
vector [n] ADAPTIVITY; vector [n] TEACHINT;
vector [n] SCREADDIFF; vector [n] SCREADCOMP;
}

parameters {
real alpha;
real beta1; real beta6;
real beta2; real beta7;
real beta3; real beta8;
real beta4; real beta9;
real beta5; real beta10;
real<lower=0> sigma;
real<lower=0> pi;
real<lower=0,upper=1> lambda; // Needed for Spike and Slab
real<lower=0> tau; // Needed for Spike and Slab
real<lower=0> c;

}

model {
real mu[n];
for (i in 1:n)
mu[i] = alpha + beta1Female[i] + beta2ESCS[i] + beta3METASUM[i]
+ beta4PERFEED[i] + beta5JOYREAD[i] + beta6MASTGOAL[i]
+ beta7ADAPTIVITY[i] + beta8TEACHINT[i]
+ beta9SCREADDIFF[i] + beta10SCREADCOMP[i] ;

// Spike and Slab Priors
sigma ~ cauchy(0,1);
lambda ~ bernoulli(pi);
tau ~ cauchy(0,1);
c ~ cauchy(0,10);
alpha ~ uniform(-3,3);
beta1 ~ normal(0, clambda); beta6 ~ normal(0, clambda);
beta2 ~ normal(0, clambda); beta7 ~ normal(0, clambda);
beta3 ~ normal(0, clambda); beta8 ~ normal(0, clambda);
beta4 ~ normal(0, clambda); beta9 ~ normal(0, clambda);
beta5 ~ normal(0, clambda); beta10 ~ normal(0, clambda);

// Likelihood
readscore ~ normal(mu, sigma);

}

// For posterior predictive checking and loo cross-validation
generated quantities {
vector[n] readscore_rep;
vector[n] log_lik;
for (i in 1:n) {
readscore_rep[i] = normal_rng(alpha + beta1Female[i] + beta2ESCS[i] + beta3METASUM[i]
+ beta4PERFEED[i] + beta5JOYREAD[i] + beta6MASTGOAL[i]+ beta7ADAPTIVITY[i] + beta8TEACHINT[i]
+ beta9SCREADDIFF[i] + beta10SCREADCOMP[i], sigma);

log_lik[i] = normal_lpdf(readscore[i] | alpha + beta1*Female[i] + beta2*ESCS[i] + beta3*METASUM[i]
        + beta4*PERFEED[i] + beta5*JOYREAD[i] + beta6*MASTGOAL[i]+ beta7*ADAPTIVITY[i] + beta8*TEACHINT[i]
        + beta9*SCREADDIFF[i] + beta10*SCREADCOMP[i], sigma);

}
}

"

I’m getting the error message, which I am not sure I understand.

int ~ bernoulli(real)
int ~ bernoulli(real[ ])
int ~ bernoulli(vector)
int ~ bernoulli(row_vector)
int[ ] ~ bernoulli(real)
int[ ] ~ bernoulli(real[ ])
int[ ] ~ bernoulli(vector)
int[ ] ~ bernoulli(row_vector)

Real return type required for probability function.
error in ‘model7a11534daddd_57a8adcf55dc2a129787780e00b7e4db’ at line 38, column 29

36:   // Spike and Slab Priors
37:       sigma ~ cauchy(0,1);
38:       lambda ~ bernoulli(pi);
                                ^
39:       tau ~ cauchy(0,1);

A couple of related questions:

I know that one can basically get the same results from the horseshoe. What priors are recommended for horseshoe get essentially spike and slab?

Is there a simple example of Stan code for the “Finnish” horseshoe for the Gaussian linear model?

(I realize there are simpler ways to set up the model using matrix notation, but this analysis is for demonstration purposes and I want the students to see the detail) . :-)

2

Have a look at @betanalpha’s excellent case study on different approaches to sparsity-inducing priors, all with Stan examples

1 Like
3

Forgot the link! Here it is: Sparsity Blues

4

Hi Andrew,

Thanks for the quick reply. I’ve seen this (and another similar link) from Betancourt and they have been really helpful. My issue is two fold. First, from my post, I’m getting error message I don’t quite understand. Second, I’m not sure that this is the way one directly sets up a spike and slab prior, so if you (or others) have a better suggestion, please share. Similarly, perhaps there is a choice of hyperparameter in the horseshoe prior which would essentially yield the results from a spike and slab prior.

Thanks

5

The issue is that you can’t directly construct a spike-and-slab prior in Stan since it requires a discrete parameter. This is why you’re getting the error with the bernoulli prior, because you’re attempting to use it a with a continuous parameter when it’s only defined for integer data.

I’m not sure of a horsehoe equivalent for the spike-and-slab, @avehtari would you know?

6

lambda ~ bernoulli(pi);

lambda must be integer for this to be true.

7

See Figures 1 and 2 in Sparsity information and regularization in the horseshoe and other shrinkage priors

See Appendix C in the same paper.

1 Like