From JAGS to Stan - 2 parameter probit IRT with covariates

[Jake Johnson]

Greetings Stan users,

I recently began experimenting with Stan as a tool to simulate two parameter IRT models (probit link) with congressional voting data.  Originally I started with the pscl package in R and then moved on to JAGS via rjags.  The current model I am interested in adds additional covariates to the standard model as in:

P(y_ij = 1) = Phi( beta_j * x_i - alpha_j + gamma * z_ij)

where z_ij represents data (legislator and constituent characteristics).  I have written the code for this model (attached) but after running it a number of times the model doesn't mix well.  It moves to the prior means for the x_i and stays there while all the other parameters remain at their initial values (at least up to 500 iterations which took around 3 hours).  

I am certain the code works since I have run it with z_ij equal to all zeros and the output is comparable to using ideal() in pscl which is the exact same model.  I was also successful with a canned example I made for this post (see attached).  So I am thinking this is related to my data.  I have had success using JAGS with this data (code attached) and it does return sensible results.  However it is very slow and so I would really like to get this working in Stan to compare.  I am hoping someone might be able to provide some advice or tips to get things moving.  Better starting values?  Rescaling the data?

As an aside, can someone provide some intuition for the "vanishing density" error?  Either of the following triggers it:

-  Passing initial values for alpha and beta using the coefficients from probit regressions of y_j on the initial x_i and an intercept (one regression per j).  This is a default in ideal().

-  Passing initial values for gamma using the coefficients from a probit regression of y on z (no intercept) with all other initial values at zero

Thanks in advance!

Jake

[Bob Carpenter]

Those didn't look like at all the same models in Stan and JAGS
given that you give beta a direct multinormal prior in the
JAGS approach.

For speed (and robustness) in Stan, you're better off using

y ~ bernoulli_logit(beta[rc[r]]*xi[id[r]] - alpha[rc[r]] + zg[r]);

That assumes the parameter's on the logit scale. It usually doesn't make
much difference in modeling other than scale of coefficients.

I didn't trace through the values for your mu_* and sigma_* constants
or the form of your predictors, but you have to be very careful about
identifiability with these models with independent additive factors.

If the predictors are not standardized, doing that could help. But
the problem in these models is usually the non-identifiability.

Whatever you do, it would help to generate fake data based on your model
and see if the program can fit it. That'll test that everything is OK before
moving to real data.

More below.

> On Jun 8, 2015, at 10:11 PM, Jake Johnson <jak...@gmail.com> wrote:
>
> Greetings Stan users,
>
> I recently began experimenting with Stan as a tool to simulate two parameter IRT models (probit link) with congressional voting data. Originally I started with the pscl package in R and then moved on to JAGS via rjags. The current model I am interested in adds additional covariates to the standard model as in:
>
> P(y_ij = 1) = Phi( beta_j * x_i - alpha_j + gamma * z_ij)
>
> where z_ij represents data (legislator and constituent characteristics). I have written the code for this model (attached) but after running it a number of times the model doesn't mix well. It moves to the prior means for the x_i and stays there while all the other parameters remain at their initial values (at least up to 500 iterations which took around 3 hours).
>
> I am certain the code works since I have run it with z_ij equal to all zeros and the output is comparable to using ideal() in pscl which is the exact same model. I was also successful with a canned example I made for this post (see attached). So I am thinking this is related to my data. I have had success using JAGS with this data (code attached) and it does return sensible results. However it is very slow and so I would really like to get this working in Stan to compare. I am hoping someone might be able to provide some advice or tips to get things moving. Better starting values? Rescaling the data?
>
> As an aside, can someone provide some intuition for the "vanishing density" error?

It's triggered when the density is zero (log density negative infinity).
This usually results from a constraint being violated. A culprit might
be Phi_approx underflowing to 0 or overflowing to 1 with mismatched
data.

> Either of the following triggers it:
> - Passing initial values for alpha and beta using the coefficients from probit regressions of y_j on the initial x_i and an intercept (one regression per j). This is a default in ideal().
> - Passing initial values for gamma using the coefficients from a probit regression of y on z (no intercept) with all other initial values at zero

You can use print statements to track down where the
density goes to zero. You can use the constant lp__ (not in the
doc and soon to be deprecated and replaced with a function) to inspect
the current value at any point using:

print("lp==", lp__);

- Bob

[Jake Johnson]

Thanks for the the fast reply.  It isn't evident in the code I sent but the sigma's in the dmnrom() in JAGS are diagonal.  I read that this could improve performance versus a loop statement.  My impression is the same is not true for Stan - i.e. the loop is better than using multi_normal().  

I am not too concerned with identifiability per se since I am mainly interested in the rankings of the x_i.  I will also be placing a restriction on the distribution of the posterior means of the x_i and then transforming alpha_j and beta_j.  In fact, the magnitudes of gamma are not that important either since I am mainly interested in their sign.

I'll try standardizing the data and see what happens.  From there I will try fitting some simulated data.  

Thanks again,

Jake 

[Bob Carpenter]

> On Jun 8, 2015, at 11:41 PM, Jake Johnson <jak...@gmail.com> wrote:
>
> Thanks for the the fast reply. It isn't evident in the code I sent but the sigma's in the dmnrom() in JAGS are diagonal. I read that this could improve performance versus a loop statement. My impression is the same is not true for Stan - i.e. the loop is better than using multi_normal().

Absolutely --- and vectorizing is better still.

> I am not too concerned with identifiability per se since I am mainly interested in the rankings of the x_i.

It's important for Stan --- if the model's not identifiable,
Stan will run all over the posterior and the results won't make
sense. There's a discussion in the manual chapter on problematic
posteriors.

Specifically, you can't define an improper posterior and then
hope to standardize it back.

> I will also be placing a restriction on the distribution of the posterior means of the x_i and then
> transforming alpha_j and beta_j. In fact, the magnitudes of gamma are not that important either since I am mainly interested in their sign.

You need to make it at least weakly identifiable as the
density is defined in the parameters. If it isn't,
there will be problems for sampling with both HMC and Gibbs.

> I'll try standardizing the data and see what happens. From there I will try fitting some simulated data.

Let us know how it goes. A model along the lines of what you're
doing should be estimable.

- Bob

[Jake Johnson]

Bob: Just wanted to provide an update.  

Things are mixing well now on account of standardizing the data I think.  I was able to track down the "vanishing density" error as well using print().  The issue was related to the initial values I was trying to use.  Specifically, they were pushing the value of the link (xb) to such extremes that Phi_approx() was essentially returning zero's and one's.  Same error occurred using Phi().  Hence the "vanishing density."

Regarding identification, I realized my code did not include any of my priors.  I think identification won't be an issue by specifying priors for the xi of N(-1,1) and N(1,1) with the means being based on political party (democrats are -1, republicans are 1).  The xi are most definitely bimodal so this prevents sampling from values that are equivalent via rotation.  This should also provide the scale for the other parameters, etc.  Do you agree with that?

In any case, I am now seeing comparable results with JAGS and really nice performance running parallel in Windows 7.

Thanks again for your help.

Jake