Multilevel IRT/Ideal Point Models: Additive Redundant Parameterization Problem

[Rick Farouni]

Hi,

I am trying to implement the rescaling method discussed in Gelman and Hill (2007) and  Bafumi(2005) http://www.stat.columbia.edu/~gelman/research/published/171.pdf, but I keep getting this error message

" Error in function stan::prob::bernoulli_logit_log(N4stan5agrad3varE): Logit transformed probability parameter[0]  is 1.#QNAN:0, but must not be nan!Rejecting proposed initial value with zero density.

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

error occurred during calling the sampler; sampling not done "

Here is the model

data {

  int<lower=1> J;                // number of students

  int<lower=1> K;                // number of questions

  int<lower=1> N;                // number of observations

  int<lower=1,upper=J> jj[N];    // student for observation n

  int<lower=1,upper=K> kk[N];    // question for observation n

  int<lower=0,upper=1> y[N];     // correctness of observation n

}

parameters {

  vector[J] a_raw;

  vector[K] b_raw;

  vector<lower=0>[K] g_raw;

  real<lower=0> mu_g_raw;

  real mu_b_raw;

  real mu_a_raw;

  real<lower=0> sigma_a_raw;

  real<lower=0> sigma_b_raw;

  real<lower=0> sigma_g_raw;

} 

transformed parameters {

  vector[J] a;

  vector[K] b;

  vector[K] g;

  vector[N] p;

  real scale;

  real shift;

  shift <- mean(a);

  scale <- sd(a);

  a <- (a_raw - shift)/scale;

  b <- (b_raw - shift)/scale;

  g <- g_raw*scale;

  for (n in 1:N)

  p[n]<-g[kk[n]]*(a[jj[n]] - b[kk[n]]) ;

  

}

model {

  mu_a_raw ~ normal(0, 10);

  sigma_a_raw ~ cauchy(0,3);

  a_raw ~ normal(mu_a_raw, sigma_a_raw);

  

  mu_b_raw ~ normal(0, 10);

  sigma_b_raw ~ cauchy(0,3);

  b_raw ~ normal(mu_b_raw, sigma_b_raw);

  mu_g_raw ~ uniform(0, 10);

  sigma_g_raw ~ cauchy(0,3);

  g_raw ~ lognormal(mu_g_raw, sigma_g_raw);

  

  y ~ bernoulli_logit(p);

}

Thanks

Rick,

[Daniel Lee]

Hi Rick,

Can you try printing the variable 'p' either in the model or within the transformed parameters block?

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Rick Farouni]

Hi Daniel,

I get a long a  line of  "1.#QNAN"

[Daniel Lee]

The next step is to find out which part of the expression is causing p to be NaN:

- kk[n]

- g[kk[n]]

- jj[n]

- a[jj[n]]

- b[kk[n]]

Try printing each of those in the transformed parameters block after they get set. See if they're all sensible.

[Rick Farouni]

I tried printing the " p" and all the variables that precede it. I get ""1.#QNAN"". But  I think I understand why I get this error. Both "shift" and"scale" are a function of "a", but "a" is declared after "shift" and"scale", which in turn depends on  them. 

On Thursday, June 26, 2014 1:21:25 PM UTC-4, Daniel Lee wrote:

[Daniel Lee]

Yes, I think you've identified the problem. In the original model, you had:

  shift <- mean(a);

  scale <- sd(a);

  a <- (a_raw - shift)/scale;

  b <- (b_raw - shift)/scale;

The first two lines should come after the last two. Both shift and scale will be NaN at the time of evaluation in the original model.

[Rick Farouni]

Thanks Daniel,

The model I have is a slight modification of the code here https://github.com/stan-dev/stan/blob/feature/ARM/src/models/ARM/Ch.19/item_response.stan , which has a typo that someone else has pointed out

https://groups.google.com/forum/#!searchin/stan-users/rescaling/stan-users/KUGdsg6AFVI/l3SUA68p9sMJ


Rick

[Bob Carpenter]

[cc-ed to stan-dev]

Peter --- is this one of the ones you fixed?

- Bob

[Peter Li]

I didn't see this before since the model doesn't have any data so I haven't compiled/run it. In any case, it's fixed now.

[Rick Farouni]

Hi,

So I manged to fit the model but I always seem to get a very low number of effective samples for the scaled variance and location parameters.Is there a reason to be concerned here? I have attached the model and here is a part of  the output.  

                   mean se_mean  sd 2.5%  25%  50%  75% 97.5% n_eff Rhat

gamma[47]    3.0     0.0 0.4  2.3  2.8  3.0  3.3   3.9  8000  1.0

gamma[48]    2.9     0.0 0.4  2.2  2.6  2.8  3.1   3.6  8000  1.0

gamma[49]    3.3     0.0 0.4  2.6  3.0  3.2  3.5   4.1  8000  1.0

gamma[50]    3.6     0.0 1.2  1.9  2.7  3.4  4.2   6.4  2581  1.0

mu_alpha     -0.8     0.9 1.9 -4.0 -2.6 -0.6  0.7   2.5      4    2.0

mu_beta       2.9     1.8 2.8 -1.2  0.7  1.7  5.5   8.1      2     2.7

mu_gamma     0.8     0.3 0.5  0.3  0.5  0.7  1.2   1.9    3  2.5

sigma_alpha  3.4     1.1 1.9  1.1  1.8  2.9  4.6   7.9     3   1.9

sigma_beta   5.1     1.6 2.8  1.6  2.7  4.4  6.9  11.7     3   1.9

sigma_gamma  0.5     0.2 0.3  0.1  0.3  0.4  0.6   1.1   3  2.1  

Thanks,

Rick

[Bob Carpenter]

OK, I see what's going on. This is NOT a good way to code
these models. This model relies on Andrew's trick for BUGS of
specifying an ill-defined posterior to sample from, then
attempting to pull it all back into something with reasonable
n_eff and Rhat through standardization. Here's the offending
code in our example:

> transformed parameters {
> vector[J] a;
> vector[K] b;
> vector[K] g;
> vector[N] p;
> real scale;
> real shift;
>
> shift <- mean(a);
> scale <- sd(a);

> a <- (a_raw - shift) / scale;
> b <- (b_raw - shift) / scale;
> g <- g_raw * scale;
>
> for (n in 1:N)
> p[n]<-g[kk[n]]*(a[jj[n]] - b[kk[n]]);

The problem here is that the raw values are not identified.

The problem is that the posterior for the parameters is improper.
And what you see is that Stan's HMC is effcient enough to race all over
the posterior and at least realize it's undefined.

There's a chapter on problematic posteriors that goes through
an example of the effects of such non-identified models. (Look at the
plot of the long ridge, which is what happens for this kind of example.)

This model needs to be recoded so that it's properly identified.

We're working on rewriting all the ARM models, which are a pretty direct translation
of what's in Gelman and Hill, and thus not guaranteed to work.

This looks like a good place to start!

For models that don't have data we can release, we can simulate
data according to the model. I want to bring up the question of how
to simulate predictors for regressions --- it seems like for proper testing
we shouldn't make them all independent.

- Bob

On Jun 27, 2014, at 12:28 AM, Rick Farouni <rfar...@gmail.com> wrote:

> Hi,
>

> So I manged to fit the model but I always see to get a very low number of effective samples for the scaled variance and location parameters.Is that a reason to be concerned here? I have attached the model and here is a part of the output.

> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> gamma[47] 3.0 0.0 0.4 2.3 2.8 3.0 3.3 3.9 8000 1.0
> gamma[48] 2.9 0.0 0.4 2.2 2.6 2.8 3.1 3.6 8000 1.0
> gamma[49] 3.3 0.0 0.4 2.6 3.0 3.2 3.5 4.1 8000 1.0
> gamma[50] 3.6 0.0 1.2 1.9 2.7 3.4 4.2 6.4 2581 1.0
> mu_alpha -0.8 0.9 1.9 -4.0 -2.6 -0.6 0.7 2.5 4 2.0
> mu_beta 2.9 1.8 2.8 -1.2 0.7 1.7 5.5 8.1 2 2.7
> mu_gamma 0.8 0.3 0.5 0.3 0.5 0.7 1.2 1.9 3 2.5
> sigma_alpha 3.4 1.1 1.9 1.1 1.8 2.9 4.6 7.9 3 1.9
> sigma_beta 5.1 1.6 2.8 1.6 2.7 4.4 6.9 11.7 3 1.9
> sigma_gamma 0.5 0.2 0.3 0.1 0.3 0.4 0.6 1.1 3 2.1
>
>
> Thanks,
>
> Rick
>
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <irt.stan>

[Rick Farouni]

I see. So the method mentioned in Gelman and Hill(2007)  is better suited to work with BUGS, not STAN.

[Bob Carpenter]

The problem isn't the software but rather than the posterior
is improper in the parameters over which it's defined.

It doesn't really work with BUGS, either, when you look
at convergence of the raw parameters. They'll drift like
they do in Stan.

In both BUGS and Stan, the transformed parameters may look
OK, but they're still problematically defined as transforms
of variables with improper posteriors.

The manual chapter on problematic posteriors goes over the
simplest possible example of where this can happen.

- Bob

[Peter Li]

I'm confused how you would fix this. If you look at how it all gets computed in p[i] for example, g(a-b) leaves you with just (a_raw - b_raw)*g_raw so why even bother shifting and scaling? And what's the point of d? Or is this all some kind of reparameterization trick?

[Bob Carpenter]

It's a reparameterization. You have an
additive non-identifiability, which you can fix
by fixing one of the coefficients. And you have
a multiplicative non-identifiability that can get
fixed by fixing one of the scales, say the student
abilities, to normal(0,1). And then there's a sign
non-identifiability, which you can fix by constraining
the discrimination (g) to b positive. (Unless it's
an item-response model, then I'm not sure what to do
about the sign --- Andrew and Jennifer discuss in their
book.)

- Bob

[Rick Farouni]

For ideal point models, Gelman and Hill discuss an approach to deal with reflection invariance using group level regression predictors. Simon Jackman has an example ("legislators") on his website that implements the approach in Winbugs. http://jackman.stanford.edu/mcmc/index.php. I would be very interested to see how to code a similar model in STAN but for a two-dimensional case. In two dimensions, two more indeterminacy arise if I am not mistaken: label switching and rotational indeterminacy. 

[Peter Li]

Ok, but then how do you fix the non-identifiability of the raw values mentioned in the previous email? Would it be just to fix one of the a_raw and b_raw to N(0,1)?

[Bob Carpenter]

We don't have a solution for label switching, I'm afraid.
There's a discussion in the manual, but all the possible solutions
we've come up with or seen from others are flawed.

I don't know about the rotational indeterminancy. Is there
a mathematical constraint on the params that solve it?

Also, using group-level predictors doesn't solve the basic
indeterminacy problem.

- Bob

[Rick Farouni]

If I am not mistaken, for a d-dimensional model, if we specify a normal(0,1) prior on  each on the ability parameters (factors), we would still need d*(d-1)/2 additional constraints on the discrimination parameters. So for a two dimensional model, we would need to impose one restriction on the discrimination(loading) matrix. But these restrictions won't take care of reflective invariance and label switching. In Bayesian Factor Analysis, the positive lower triangular identification scheme is usually used to to mitigate rotational invariance, sign flips, and label switching. Has anyone fit a model in stan using the positive lower triangular identification scheme?

[Bob Carpenter]

This thread's drifted from typos in ARM models, right?

But the answer is that I don' tknow if anyone's done that,
but I'd be interested in seeing how it works --- these problems
with identification are really important and I'd like to include
some more realistic examples in the manual.

Is there a relatively simple reference to these models I could
start with? Or if someone sends me a working Stan model, I can
include that (with attribution, of course).

- Bob

[Rick Farouni]

Hi Bob,

For a factor model y=Lf +e, where f~ N(0,Sigma) , I think there are three important cases to consider:

1. The positive lower triangular (PLT) identification scheme mentioned here http://apps.olin.wustl.edu/faculty/zhou/GZ-RFS96.pdf , where Sigma=I (identity matrix) and the the diagonal entries of L are polarity constrained.
2.The modified PLT where the factors,f, are still orthogonal but the elements of diagonal of Sigma Aguilar, O. & West, M.(2000) Bayesian dynamic factor models and portfolio allocation.
3.. Sigma is a correlation matrix, but now instead of requiring D(D-1)/2 constraints on L, we need D(D-1) to ensure local uniqueness under rotation and the polarity constraint to ensure global rotational uniqueness.This paper lists the conditions http://link.springer.com/article/10.1007/s11336-012-9259-3

But as this paper shows https://www.stat.washington.edu/research/reports/2011/tr589.pdf , imposing constraints is problematic.

I'll try work on something and if I get a working Stan model, I'll send it to you.

Rick,