Problems with IRT model

[Emmanuel Charpentier]

I discovered Stan through Andrew Gelman's blog two days ago. The description was enticing, so I installed it and started to play around.

I tried the "blocker" example of the BUGS manual vol 1 with no problems, and got a somewhat disappointing performance, when I compared it to the performance of JAGS on thsi very simple model. See attached blocker.(R|stan) file, which is annotated (and contains data).

Most of the elapsed time was spent in compilation. This did not really surprise me : I did not expect the heavy compilation process of Stan to be as fast as the light compilation offered by JAGS. I was more surprised buy the fact that, even using only one processor, sampling from an already fit model was about 30% faster in JAGS. But again, the overhead entailed by a compilation to C++ might easily explain this in the case of such a tiny model.

So I tried my hand at a problem I *knew* to be hard or JAGS : an IRT model (part of a larger model which can take *hours* to converge correctly on JAGS). This model follows largely Gelman & Hill's 2007 textbook.

My first attempt on a dataset containing missing values failed for reasons I have not yet understood (but I understand that the released stan 1.0.0 is not yet complete...). So I wrote another model dealing only with full data. This time, I got a runtime error complaining about a negative infinite log-likelihood, which I have been unable to trace. As far as I can tell, all my initialisations are reasonable.

Could you let me know where I goofed ? See ViStan.R for data and annotated source of calls and ViStan.stan and ViStanFull for the models.

Sincerely yours,
                                                                             Emmanuel Charpentier

[Bob Carpenter]

I just love the top line comment, and thought Andrew
would too, so I thought I'd call that out first:

## Un essai d'IRT avec Stan, le nouveau joujou d'Andrew Gelman.

On 9/8/12 6:25 PM, Emmanuel Charpentier wrote:
> I discovered Stan through Andrew Gelman's blog two days ago. The description was enticing, so I installed it and started
> to play around.

> I tried the "blocker" example of the BUGS manual vol 1 with no problems, and got a somewhat disappointing performance,
> when I compared it to the performance of JAGS on thsi very simple model. See attached blocker.(R|stan) file, which is
> annotated (and contains data).

Thanks for taking the time to try it and give us
feedback.

Responding to the note in the file, you're right that
Stan hasn't implemented random sample generation in the
generated quantities block. We're working on it!
For now, you can do what you did and define it in the
model, which is definitely slower.

We'll be speeding up the model considerably when we
vectorize all the distributions. This speedup will
scale with N.

> Most of the elapsed time was spent in compilation. This did not really surprise me : I did not expect the heavy
> compilation process of Stan to be as fast as the light compilation offered by JAGS. I was more surprised buy the fact
> that, even using only one processor, sampling from an already fit model was about 30% faster in JAGS. But again, the
> overhead entailed by a compilation to C++ might easily explain this in the case of such a tiny model.

When you say 30% faster, what exactly are you measuring?
Does that include the compilation time? Are both Stan
and JAGS giving roughly similar effective sample sizes?
Everyone uses their own ESS estimator.

There's a note in the Wiki on how to up the optimization
level of R's built-in compiler to -O3 (optimization level
3). That makes a big difference to Stan. It slows
down compilation considerably, while at the same time
substantially speeding up sampling.

You should only have to compile a model once.

> So I tried my hand at a problem I *knew* to be hard or JAGS : an IRT model (part of a larger model which can take
> *hours* to converge correctly on JAGS). This model follows largely Gelman & Hill's 2007 textbook.
>
> My first attempt on a dataset containing missing values failed for reasons I have not yet understood (but I understand
> that the released stan 1.0.0 is not yet complete...). So I wrote another model dealing only with full data. This time, I
> got a runtime error complaining about a negative infinite log-likelihood, which I have been unable to trace. As far as I
> can tell, all my initialisations are reasonable.

We're also planning to implement the models from
Gelman and Hill's book. The one you did is
a pretty straightforward model with a discrimination
term.

> Could you let me know where I goofed ?

I don't see anything obviously wrong in terms of
what Stan expects.

Did you try letting Stan do the initialization?
Stan initializes randomly uniformly on [-2,2]
on the unconstrained scale.

Positive-constrained parameters like your gamma
get a log transform to unconstrain them, and exp
transform to come back; so initializing around 0 on
the unconstrained scale initializes around 1.

You can set init=0 to get all-zero inits on the
unconstrained scale.

The binomial(1,inv_logit(phi)) could be replaced with
bernoulli_logit(phi); I saw you did that in some other
models. But Stan's pretty efficient at figuring out
which parts of the binomial calculation to do to get
the log probability up to a proportion.

> ## Error in function stan::prob::binomial_log<N4stan5agrad3varE>(N4stan5agrad3varE): Probability, theta, is nan:0,
> but must be finite!

The error you report is likely to arise if theta is
never initialized.

## SAMPLING FOR MODEL 'IRT1' NOW (CHAIN 1).
## Erreur dans sampler$call_sampler(args_list[[i]]) : variable ImpCp missing

I hope this is self explanatory. It's expecting to find
a variable and not finding it. If you want to do missing
data models in Stan, they have to be coded differently than
in BUGS/JAGS. There's a section in the programming part of
the manual with some examples.

I'll code up a version of the model you were trying
to debug and post it back to the list.

- Bob

[Andrew Gelman]

Bob:

So, are you saying that's ok to run Stan without initializing anything (because then Stan will choose the inits), and it's ok to initialize all the variables (then Stan will just use those inits), but it's not ok to initialize only _some_ of the variables?  If so, perhaps we should set Stan so that, if the user initializes some variables but not all, Stan will automatically initialize the variables that the user did not initialize.

Does this make sense?

A

P.S.  Joujou, indeed!

[Bob Carpenter]

On 9/8/12 7:42 PM, Andrew Gelman wrote:
> Bob:
> So, are you saying that's ok to run Stan without initializing anything (because then Stan will choose the inits), and
> it's ok to initialize all the variables (then Stan will just use those inits), but it's not ok to initialize only _some_
> of the variables? If so, perhaps we should set Stan so that, if the user initializes some variables but not all, Stan
> will automatically initialize the variables that the user did not initialize.
> Does this make sense?

Absolutely. It won't be trivial given the read
routines, which are shared between data (which must
be there) and inits (which I agree would be better
the way you describe). But it should be doable because
we know all the types and sizes of all parameters, so
we can calculate the number of unconstrained params to
fill in. I can just drop flags through the code indicating
whether to fill in randomly or not.

There are quite a few other priorities first, like making
sure matrix ops won't cause R to segfault and getting
user model types going at least in example form!

> P.S. Joujou, indeed!

Mitzi set me straight on the translation. It was
a faux ami problem, as the French say --- I
was confusing "joujou" (toy) with "juju" (witchcraft).

Programmers often use "fu" (as in Kung Fu)
or "juju" (as in witchcraft) in similar contexts,
e.g., "that's some wicked R fu, dude" and
"i'm going to need crazy template juju for this".

Not to be confused with the generic variable "foo", of course.

- Bob

[Emmanuel Charpentier]

Dear Bob,

Thhank you for your prompt answer.

Le dimanche 9 septembre 2012 01:34:23 UTC+2, Bob Carpenter a écrit :

I just love the top line comment, and thought Andrew
would too, so I thought I'd call that out first:

## Un essai d'IRT avec Stan, le nouveau joujou d'Andrew Gelman.

I meant "joujou" (toy, thingie), but with some strong connotations of "juju" (charm, black magic). This apparently didn't came across the linguistic barrier, alas.
 

On 9/8/12 6:25 PM, Emmanuel Charpentier wrote:
> I discovered Stan through Andrew Gelman's blog two days ago. The description was enticing, so I installed it and started
> to play around.

> I tried the "blocker" example of the BUGS manual vol 1 with no problems, and got a somewhat disappointing performance,
> when I compared it to the performance of JAGS on thsi very simple model. See attached blocker.(R|stan) file, which is
> annotated (and contains data).

Thanks for taking the time to try it and give us
feedback.

Responding to the note in the file, you're right that
Stan hasn't implemented random sample generation in the
generated quantities block.  We're working on it!
For now, you can do what you did and define it in the
model, which is definitely slower.

We'll be speeding up the model considerably when we
vectorize all the distributions.  This speedup will
scale with N.

> Most of the elapsed time was spent in compilation. This did not really surprise me : I did not expect the heavy
> compilation process of Stan to be as fast as the light compilation offered by JAGS. I was more surprised buy the fact
> that, even using only one processor, sampling from an already fit model was about 30% faster in JAGS. But again, the
> overhead entailed by a compilation to C++ might easily explain this in the case of such a tiny model.

When you say 30% faster, what exactly are you measuring?
Does that include the compilation time?  Are both Stan
and JAGS giving roughly similar effective sample sizes?
Everyone uses their own ESS estimator.

Compare the times elapsed for *second* run (i.e. sampling from an already compiled model).

There's a note in the Wiki on how to up the optimization
level of R's built-in compiler to -O3 (optimization level
3).  That makes a big difference to Stan.  It slows
down compilation considerably, while at the same time
substantially speeding up sampling.

I'd rather not : the "-O2" option is hardwired in the R stan() function, and, as far as I know, there is but one way to change it, which is to hack stan(). Which I'm a bit reluctant to do at the time, since, IIUC, a 1.0.1 version is close to release.
 

You should only have to compile a model once.

So I noticed, as exemplified in the annotated  R file, which contains timings for a second run (that's where my "30% slower for sampling" figure comes from).

OK. Just a little problem : THERE IS NO theta VARIABLE IN THE MODELS. I suspect that the compiler is goofing in building it's error message.

## SAMPLING FOR MODEL 'IRT1' NOW (CHAIN 1).
## Erreur dans sampler$call_sampler(args_list[[i]]) : variable ImpCp missing

I hope this is self explanatory.  It's expecting to find
a variable and not finding it.

ImpCp, representing *imputed values* for missing Cps..
 

If you want to do missing
data models in Stan, they have to be coded differently than
in BUGS/JAGS.  There's a section in the programming part of
the manual with some examples.

compare ViStan.stan and ViStanFull.stan... I got inspiration from the manual for the handling of missing variables in ViStan.stan (the one which aimed to fit a model to data with incomplete observations).

I'll code up a version of the model you were trying
to debug and post it back to the list.

I'll be very glad !

- Bob

Again, thank you for a very fast and on-the-point answer.

                                                                  Emmanuel Charpentier

[Emmanuel Charpentier]

Le dimanche 9 septembre 2012 01:51:44 UTC+2, Bob Carpenter a écrit :

On 9/8/12 7:42 PM, Andrew Gelman wrote:
> Bob:
> So, are you saying that's ok to run Stan without initializing anything (because then Stan will choose the inits), and
> it's ok to initialize all the variables (then Stan will just use those inits), but it's not ok to initialize only _some_
> of the variables?  If so, perhaps we should set Stan so that, if the user initializes some variables but not all, Stan
> will automatically initialize the variables that the user did not initialize.
> Does this make sense?

Absolutely.  It won't be trivial given the read
routines, which are shared between data (which must
be there) and inits (which I agree would be better
the way you describe).  But it should be doable because
we know all the types and sizes of all parameters, so
we can calculate the number of unconstrained params to
fill in.  I can just drop flags through the code indicating
whether to fill in randomly or not.

I second this position : initializing just the parameters that *must* have a  well-defined value is quite handy when you have a ton of parameters and only a few need precise initialization. I'd like to add that JAGS default behaviour for initialization (to 0, rather thn sampling from the prior as in OpenBUGS) seems to me more reasonable in many cases. I love the stan() R function range of options : "0", "random" or user-specified values. The best of ll worlds would be a named list of parameters, and a default, which could be either "0" or "random". But that might be treated at the R level.

Thank you again.

                                                Emmanuel Charpentier

[Emmanuel Charpentier]

Wups. Typo (mouso ?) correction

[ Bandwidth savings... ]

      ## SAMPLING FOR MODEL 'IRT1' NOW (CHAIN 1).

## Erreur dans sampler$call_sampler(args_list[[i]]) : variable ImpCp missing

I hope this is self explanatory.  It's expecting to find
a variable and not finding it.

ImpCp, representing *imputed values* for missing Cps..

I meant : ImpCp, representing *imputed values* for missing Cps, is declared as a *parameter* in the "parameters' section. That results from my understanding of the relevant section of the (impressive !) user's manual. Am I wrong ?

[ Other bandwidth savings ... ]

                                                                                  Emmanuel Charpentier

[Jiqiang Guo]

A couple of comments below.

Jiqiang

 

I'd rather not : the "-O2" option is hardwired in the R stan() function, and, as far as I know, there is but one way to change it, which is to hack stan(). Which I'm a bit reluctant to do at the time, since, IIUC, a 1.0.1 version is close to release.

stan() itself does not hardcode any optimization level for compiling c++ code: the configuration is up to the default that is used by R CMD SHLIB.  This would be the same in 1.0.1 version.  And I do not consider changing the optimization level a hack. 

ImpCp is declared as a *discrete* parameter, which is not supported in Stan for now.  This might be the reason.  

                                                                                  Emmanuel Charpentier

[Emmanuel Charpentier]

Dear Jiqiang,

Thank your for taking the time to answer me. A couple answers below.

Le dimanche 9 septembre 2012 15:36:57 UTC+2, Jiqiang Guo a écrit :

A couple of comments below.

Jiqiang

 

I'd rather not : the "-O2" option is hardwired in the R stan() function, and, as far as I know, there is but one way to change it, which is to hack stan(). Which I'm a bit reluctant to do at the time, since, IIUC, a 1.0.1 version is close to release.

stan() itself does not hardcode any optimization level for compiling c++ code: the configuration is up to the default that is used by R CMD SHLIB.  This would be the same in 1.0.1 version.  And I do not consider changing the optimization level a hack. 

OK. I'll look into changing that cleanly (maybe on a per-sesion basis, through R's options()).

On Sun, Sep 9, 2012 at 5:48 AM, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:

Wups. Typo (mouso ?) correction

[ Bandwidth savings... ]

 

I meant : ImpCp, representing *imputed values* for missing Cps, is declared as a *parameter* in the "parameters' section. That results from my understanding of the relevant section of the (impressive !) user's manual. Am I wrong ?

[ Other bandwidth savings ... ]

ImpCp is declared as a *discrete* parameter, which is not supported in Stan for now.  This might be the reason.  

Aha. I overlooked that in the manual. So currently Stan can't impute *discrete* missing data. Nor assess the distributions of *discrete* parameters (so much for model choice algorithms such as umbrella sampling or hierarchical model choice (à la Carlin & Chib (1995) or otherwise)...).

Any idea of the ETA for a discrete sampler ?

Again, thank you very much !

                                                                     Emmanuel Charpentier

[Bob Carpenter]

On 9/9/12 10:41 AM, Emmanuel Charpentier wrote:
> Dear Jiqiang,
>
> Thank your for taking the time to answer me. A couple answers below.
>
> Le dimanche 9 septembre 2012 15:36:57 UTC+2, Jiqiang Guo a écrit :
>
> A couple of comments below.
>
> Jiqiang
>
> I'd rather not : the "-O2" option is hardwired in the R stan() function, and, as far as I know, there is but one
> way to change it, which is to hack stan(). Which I'm a bit reluctant to do at the time, since, IIUC, a 1.0.1
> version is close to release.
>
>
> stan() itself does not hardcode any optimization level for compiling c++ code: the configuration is up to the
> default that is used by R CMD SHLIB. This would be the same in 1.0.1 version. And I do not consider changing
> the optimization level

> <http://code.google.com/p/stan/wiki/RStanGettingStarted#Modify_C++_compiler_optimization_level> a hack.

>
>
> OK. I'll look into changing that cleanly (maybe on a per-sesion basis, through R's options()).

I don't think that's possible. If you do figure out a way
to do it, could you let us know?

As a side benefit, -O3 will make other R apps compiled
with R CMD SHLIB faster :-)

Here's a thread on stackoverflow on the topic:

http://stackoverflow.com/questions/5789673/r-c-optimization-flag-when-using-the-inline-package

> Aha. I overlooked that in the manual. So currently Stan can't impute *discrete* missing data. Nor assess the
> distributions of *discrete* parameters (so much for model choice algorithms such as umbrella sampling or hierarchical
> model choice (à la Carlin & Chib (1995) or otherwise)...).

Right. Unless you can marginalize out the
discrete parameter, Stan can't handle it.

I'm about to put a compile-time check to flag attempted
uses of discrete params in the parser, which I forgot to
do the first time around.

> Any idea of the ETA for a discrete sampler ?

We'll probably write a simple Gibbs for finite
and slice for infinite sampler. Probably some
time in the next six or nine months.

- Bob

[Bob Carpenter]

n 9/9/12 5:26 AM, Emmanuel Charpentier wrote:
> Dear Bob,
>
> Thhank you for your prompt answer.
>
> Le dimanche 9 septembre 2012 01:34:23 UTC+2, Bob Carpenter a écrit :

> > ## Error in function stan::prob::binomial_log<N4stan5agrad3varE>(N4stan5agrad3varE): Probability, theta, is
> nan:0,
> > but must be finite!
>
> The error you report is likely to arise if theta is
> never initialized.
>
> OK. Just a little problem : THERE IS NO theta VARIABLE IN THE MODELS. I suspect that the compiler is goofing in building
> it's error message.

I see the confusion here.

We should make that error message clearer.
That's just a generic report from binomial_log, where
the probability-of-success parameter is called theta.
We should just get rid of the names altogether and refer
to the names of the variable, like chance-of-success
or scale or location.

The problem for us is at the point that warning (and it's
just a warning) gets reported, we've lost the link to the
underlying variable in the Stan model. It would be possible
to fix this by passing arguments representing the names of
all the variables to the distributions or by trying to catch
errors in the Stan code, but both are a lot of work and
would make the code very ugly. Probably worth it in the
longer term, as this is going to be very frustrating for
users.

- Bob

[Bob Carpenter]

That's what we're aiming for. It'll take a bit of
work and there are other fires to put out at the
moment on basic error checking that we're going to
address first.

The biggest hurdle for the implementation is that the
inits are on the unconstrained scale.

I don't think JAGS always initializes to 0. It says
in the manual that it takes reasonable values.
People like to use very fat priors even when they
have a good sense of what inits should be. Maybe
Andrew has something to say about that -- he likes
dispersed initializations.

RStan and Stan use the same basic C++ code, but
they have different command-line controls for marshaling
the arguments. The C++ needs the data to be given
to the model's class constructor and the constructor
computes the transformed data, all of which then remains
fixed throughout sampling.

So far, we haven't had problems with Stan converging,
even from pretty far-out inits. If they do get so
bad that they result in NaN or infinite values, we
resample.

- Bob

[Andrew Gelman]

On Sep 9, 2012, at 11:15 PM, Bob Carpenter wrote:

> I don't think JAGS always initializes to 0. It says
> in the manual that it takes reasonable values.
> People like to use very fat priors even when they
> have a good sense of what inits should be. Maybe
> Andrew has something to say about that -- he likes
> dispersed initializations.

If you're doing Gibbs, it's a good idea to avoid deterministic inits. For example, if you set theta_1,…,theta_J=0, then you might estimate tau=0, etc. This is probably less of an issue with Nuts or HMC. Still, I'd generally prefer some randomness in the inits (except when we want to do some pure comparison as in our speed comparisons with Jags).

[Emmanuel Charpentier]

Dear Bob,

Le lundi 10 septembre 2012 05:11:30 UTC+2, Bob Carpenter a écrit :

 

> OK. Just a little problem : THERE IS NO theta VARIABLE IN THE MODELS. I suspect that the compiler is goofing in building
> it's error message.

I see the confusion here.

We should make that error message clearer.
That's just a generic report from binomial_log, where
the probability-of-success parameter is called theta.
We should just get rid of the names altogether and refer
to the names of the variable, like chance-of-success
or scale or location.

The problem for us is at the point that warning (and it's
just a warning) gets reported, we've lost the link to the
underlying variable in the Stan model.  It would be possible
to fix this by passing arguments representing the names of
all the variables to the distributions or by trying to catch
errors in the Stan code, but both are a lot of work and
would make the code very ugly.  Probably worth it in the
longer term, as this is going to be very frustrating for
users.

Indeed...

Re-reading my code (both R and Stan), I have been unable to identify the uninitialized parameters. I'm probably becoming myopic with this code, so I'll give it a rest, and try to debug it a few days later.

Thank you again !

                                                                             Emmanuel Charpentier

[Matt Hoffman]

>> Any idea of the ETA for a discrete sampler ?
>
>
> We'll probably write a simple Gibbs for finite
> and slice for infinite sampler. Probably some
> time in the next six or nine months.
>
> - Bob

Note that there are no plans to write an *efficient* Gibbs or slice
sampler in the immediate future. The first version that goes out will
not be able to exploit any independence structure that may be present
in the model, which means it'll be dramatically slower than a
custom-built Gibbs sampler.

We were talking at one point about making it possible for users to
write their own Gibbs samplers for interleaving with NUTS samples,
though.

Matt