Nominal Response Model from BUGS - (slow code)

[Jason Engle]

Hello!

I am relatively new to STAN, so bear with me.

I'm trying to translate a Nominal Response model (Bock 1972) from a working BUGS code that runs interminably slowly.  I've tried STAN with the same code, but the performance fares equally as poorly, if not worse, which leads me to suspect I'm doing something suboptimal in the world of STAN, but have had no luck in figuring it out over the course of a day of putzing.

My primary thought is that (from forums and such):

resp[k,j]'~categorical(p[j,k]);

should be done all at once as vectors.  That's fine, but I run across an error when I form it as resp'<-categorical(p):

no matches for function name="transpose"
    arg 0 type=int[2]
available function signatures for transpose:
0.  transpose(vector) : row vector
1.  transpose(row vector) : vector
2.  transpose(matrix) : matrix
expression is ill formed

Probably something painfully simple that I've been unable to uncover or a gap in my understanding.  Thanks for any input.

Best,
Par Jason Engle

---------------------------------

data{
int N;
int T;
int NC;
int resp[T,N];
real lam[T,NC];
}

parameters{
real theta[N];
real thetaN[T];
matrix[T,NC] z;
}

model
{
vector[NC] tt[N,T];
vector[NC] p[N,T];
matrix[N,T] sumtt;

// Nominal model  - Bock 1972

// iterate over respondents
for (j in 1:N) {
//Iterate over items
for (k in 1:T) {
// Iterate over categories
for (l in 1:10) {
tt[j,k,l]<- exp(z[k,l]+lam[k,l]*(fabs(theta[j]-thetaN[k])));
}
sumtt[j,k]<-tt[j,k,1]+tt[j,k,2]+tt[j,k,3]+tt[j,k,4]+tt[j,k,5]+tt[j,k,6]+tt[j,k,7]+tt[j,k,8]+tt[j,k,9]+tt[j,k,10];
}
}

// iterate over respondents
for (j in 1:N) {
//Iterate over items
for (k in 1:T) {
// Iterate over categories
for (l in 1:10) {
p[j,k,l]<-tt[j,k,l]/sumtt[j,k];
}}}


//Iterate over items
for (j in 1:N) {
for (k in 1:T) {
resp[k,j]'~categorical(p[j,k]);
}
}

//Iterate over items
for (k in 1:T) {
// Iterate over categories
for (l in 1:10) {
z[k,l]~normal(0,1);
}
}

theta~normal(0,1);
thetaN~normal(0,1);

[Bob Carpenter]

On 1/21/14, 3:33 PM, Jason Engle wrote:

> I'm trying to translate a Nominal Response model (Bock 1972) from a working BUGS code that runs interminably slowly.
> I've tried STAN with the same code, but the performance fares equally as poorly, if not worse,

Is it mixing slowly or just running slowly per iteration?
Speeding up the implementation won't help with the mixing.
Often poor mixing arises because a model's not well identified.

What you want to compare is time per effective sample, not
time per iteration. The number of effective samples you get out
in the end is all that matters in terms of posterior inference.
The residual error in estimating posterior means is reported
as MCMC std error in Stan output, and that's just (estimated)
posterior variance divided by the sqrt of the number of effective
samples.

> which leads me to suspect
> I'm doing something suboptimal in the world of STAN, but have had no luck in figuring it out over the course of a day of
> putzing.

> My primary thought is that (from forums and such):
>
> resp[k,j]'~categorical(p[j,k]);
>
> should be done all at once as vectors. That's fine, but I run across an error when I form it as resp'<-categorical(p):
>
> no matches for function name="transpose"
> arg 0 type=int[2]
> available function signatures for transpose:
> 0. transpose(vector) : row vector
> 1. transpose(row vector) : vector
> 2. transpose(matrix) : matrix
> expression is ill formed

Stan's strongly typed and each function is very picky about
the type of arguments it works on. Your variable "resp"
is a 2D integer array (i.e., int[2]), whereas the transpose
operator (') only applies to vectors, row vectors, or matrices.

I rewrote the model for you more efficiently. To understand
the changes, note that if x is a matrix then x[n] is a row
vector, and look up the definition of categorical_logit, which
will lead you to softmax. You were implementing softmax yourself
in code. But using the correct parameterization of categorical
will be much faster. Also note that I vectorized the prior for z.

And please check my work on the transforms
to make sure that both I did it right and you understand what I did.
(We should have a version of categorical_logit that takes a row vector
as a parameter, and then you wouldn't need that final transposition.)

You could definite temporaries for the argument to categorical_logit,
but it's not necessary in Stan.

- Bob

data {
int N;
int T;
int NC;
int resp[T,N];

matrix[T,NC] lam;

}
parameters {
real theta[N];
real thetaN[T];
matrix[T,NC] z;
}
model {

theta ~ normal(0,1);
thetaN ~ normal(0,1);

for (t in 1:T)
z[t] ~ normal(0,1);

for (n in 1:N)
for (t in 1:T)
resp[t,n] ~ categorical_logit((z[t]
+ lam[t] * fabs(theta[n] - thetaN[t]))');
}

[Jason Engle]

Hi Bob,

Thanks for the rewrite.

I'm implementing your rewrite now, which seems a more elegant way to the same end, so I do appreciate it.  It is still on the slow side however still - in spite of your elegance here.

The difficulty is, unfortunately, both in mixing and slow sampling.  I expected the mixing, which is the same in WINBUGS, just a lot of parameters with high auto-correlation.  However, slow sampling is what really concerns me in order to make running this kind of model viable.  My understanding is, to compensate for the proliferation of parameters in these types of models, I'll need to compensate with a larger burnin and sample (which seems to be a hallmark of most IRT models I've tried to run).

Thanks again.  I'll keep plugging/experimenting with different things.  I've tried to "vectorize" the response data to no avail thus far, and perhaps "categorical logit" may help me here?

Best,
Jason

[Andrew Gelman]

No burnin.  We only do warmup now.

--
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/groups/opt_out.

[Bob Carpenter]

Andrew's always hung up the lexicon, even though
I'm the linguist in the group.

We haven't had any problem getting IRT models to mix
and sample well, even at large scale.

The best thing you can do if the model's not well
identified (often the problem with correlated parameters)
is come up with some kind of reparameterization that is
well identified. But keep in mind it has to be well
identified on the unconstrained scale for really
optimal behavior, though you can often get better mixing
in practice if you just define some transformed parameters
or generated quantities of interest that are identified.

Our multivariate distributions like categorical_logit
aren't vectorized, so you're not going to have luck
there. The categorical_logit will should make things a lot
faster compared to your computing softmax yourself in
Stan code --- it precomputes all the derivatives rather
than figuring them all out with the chain rule at run time.
The other vectorization shouldn't make much of a difference
other than simplifying the code.

- Bob

>> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.

[Jason Engle]

Hi Bob,

Excellent!  In other words, I'm mentally barking up the wrong tree on vectors - and that my own resignation on poor mixing and slow sampling is likely ill-founded.

It'll take some thought when I get the mental space to return in a day or two, but it's encouraging that it may eventually be worth the effort if I think through it and try alternate parameterizations.

Thanks for the thoughts and help.

Best,
Jason

>> stan-users+...@googlegroups.com <mailto:stan-users+unsub...@googlegroups.com>.

[Jason Engle]

Hi Bob,

Since no good deed goes unpunished, I thought I'd run this by you.

The preceding model is, alas, unidentified.  Both good news and bad news, but it brings up a corollary of difficulty. 

I want (need) to constrain certain parameters in this model such that my "z" parameters (matrix[T,NC] z;) sum to 0 for every T - a k-1 categorical constraint (duh).  Or set one to 0 so my k-1 parameters are set.  My difficulty is where to do such an artful thing exactly? 

I'd also like to peg a number of my thetaN[T] parameters so that they are constrained positive or negative to prevent scale reversal which I get in my chains (e.g. real<upper=0> thetaN[1]; AND real<lower=0> thetaN[2];).  This should "nail my scale to the wall" so to speak.  I'm simply unsure how to nail it down.

It has to be easy, I'm just too dense at this moment for that specific light to come on after trying to emulate 9.2 in the Reference Manual.

Best,
Jason

[Bob Carpenter]

There are several ways you can identify a model like this.
But they'll all be pretty ugly to code in Stan as is.

What we really need is a mean_zero vector data type (a K-vector
gets K-1 free parameters). And a first_element_zero vector data type
(also K-1 free params). In the meantime...

One approach is to define each row as a scaled and offset
simplex:

parameters {
simplex[NC] z_raw[T];
real<lower=0> z_scale[T];
...
}
transformed parameters {
matrix[T,NC] z;
...

for (t in 1:T)

z[t] <- z_scale[t] * (z_raw[t] - 1.0 / T);

...
}

This lets the z be anything but insures each z[t] sums to 0.
You then want to define priors on the z_raw (perhaps uniform)
and on the z_scale (whatever makes sense).

The other standard approach is to fix one of the values in each
row, such as z[t,1] == 0. To do this, you need to define (NC-1)-vectors
that are free, then put them back together in z. Not very pretty.
And also asymmetric in terms of the priors.

To constrain the sign of just one entry in an M-vector, you have
to define it like this:

parameters {
vector[M-1] a2_M;
real<lower=0> a1;
}
transformed parameters {
vector[M] a;
a[1] <- a1;
for (m in 2:M)
a[m] <- a2_M[m-1];
}

File this under "stupid Stan tricks" :-)

- Bob

> >> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.
> >> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

> >
> > --
> > 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/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Andrew Gelman]

There's also the post processing trick, which goes as follows:

Suppose you have parameters theta_1,…,theta_J whose overall level is nonidentified jointly with some other parameter. For example, there is a term in the model, mu + theta_j, and the likelihood is unchanged if you add a const to mu and subtract that const from the vector theta.

You get identification in two steps:
(1) Assign weak priors to the theta_j's, just so they have a proper posterior. For example, independent normal priors with mean 0 and scale 10.
(2) Create derived parameters:
mu.adj <- mu + mean(theta)
theta.adj_j <- theta_j - mean(theta).
You run the HMC on (mu,theta), but you save and monitor convergence of (mu.adj, theta.adj).
This is what we call the additive redundant parameterization in the ARM book.

Ideally this could be done in a more automatic way (I think this might be what Bob is saying when he is talking about a mean_zero vector data type). As it is, it's a bit of a hack because you have to define these additive parameters manually. But it can work.

[Bob Carpenter]

Michael should weigh in here, but I don't think the post-processing
trick is a great idea with our adaptive HMC, because it's not
where the sampling will happen. We'd rather have the parameters
defined in the parameters block be well identified.

If someone wants to actually evaluate both options empirically,
that'd be great (and it'd be a paper, as Andrew likes to say).

- Bob

[Andrew Gelman]

Step (1) below assures that the parameters will be well identified.

[Michael Betancourt]

Right — ideally we would have a data type that transforms to an unconstrained space of
smaller dimension as we do for simplices and covariance matrices. Actually, if all you
need is for the sum to be zero then why not just use a simplex (n variables constrained
to sum to 1) and then subtract one?

[Bob Carpenter]

A twofer followup; see below.

On 1/23/14, 11:27 PM, Michael Betancourt wrote:
> Right -- ideally we would have a data type that transforms to an unconstrained space of

> smaller dimension as we do for simplices and covariance matrices. Actually, if all you
> need is for the sum to be zero then why not just use a simplex (n variables constrained
> to sum to 1) and then subtract one?

You need to subtract 1 total, so to do that symmetrically, you'd subtract 1.0 / N
from each of the N entries. (Be sure not to use 1/N, as that'll round to 0
given the way integer arithmetic works.)

On 1/23/14, 11:18 PM, Andrew Gelman wrote:
> Step (1) below assures that the parameters will be well identified.

Are they well identified or just weakly identified? I found
the IRT models we were playing with would take forever to
converge with only weakly identified parameters.

- Bob

>> On 1/23/14, 10:36 PM, Andrew Gelman wrote:
>>> There's also the post processing trick, which goes as follows:
>>>

>>> Suppose you have parameters theta_1,...,theta_J whose overall level is nonidentified jointly with some other parameter. For example, there is a term in the model, mu + theta_j, and the likelihood is unchanged if you add a const to mu and subtract that const from the vector theta.

[Andrew Gelman]

All that matters is that the postprocessed quantities be well identified.

[Marcus Brubaker]

On Thu, Jan 23, 2014 at 5:27 PM, Michael Betancourt <betan...@gmail.com> wrote:

Actually, if all you
need is for the sum to be zero then why not just use a simplex (n variables constrained
to sum to 1) and then subtract one?

That's not quite the same because simplex values are also constrained to be positive.

[Marcus Brubaker]

On Thu, Jan 23, 2014 at 5:36 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

All that matters is that the postprocessed quantities be well identified.

Non-identifiable in the parameters can, at least in theory, play havoc with the adaptation and sampling efficiency.

Cheers,

Marcus

 

On Jan 23, 2014, at 11:32 PM, Bob Carpenter wrote:

On 1/23/14, 11:18 PM, Andrew Gelman wrote:
> Step (1) below assures that the parameters will be well identified.

Are they well identified or just weakly identified?  I found
the IRT models we were playing with would take forever to
converge with only weakly identified parameters.

[Bob Carpenter]

On 1/23/14, 11:44 PM, Marcus Brubaker wrote:

> On Thu, Jan 23, 2014 at 5:27 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>
> Actually, if all you
> need is for the sum to be zero then why not just use a simplex (n variables constrained
> to sum to 1) and then subtract one?
>
>
> That's not quite the same because simplex values are also constrained to be positive.

But if you subtract 1.0/N from each, they won't all be positive. You need
to then scale by a positive-constrained variable to get an arbitrary vector that
sums to 1. So to get a transformed parameter theta that sums to 1 in a way
that's well identified, you can do this:

parameters {
simplex[N] theta_raw;
real<lower=0> theta_scale;
...
}
transformed parameters {
vector[N] theta;
theta <- theta_scale * (theta_raw - 1.0 / N);
...

I'm having a deja vu moment, having just written this in response
to another message today or yesterday.

You then need a prior on theta_scale.

I'm still a little confused here, because this doesn't seem
over-parameterized, yet there are N free parameters, whereas
a vector that sums to 0 should really only have (N-1) free
parameters. Am I missing something here?

- Bob

[Jason Engle]

Hi Bob (and all),

Just an update. The workaround:

parameters {
   simplex[N] theta_raw;
   real<lower=0> theta_scale;
   ...
}
transformed parameters {
   vector[N] theta;
   theta <- theta_scale * (theta_raw - 1.0 / N);
   ...

it works. 

That is, it (a) accomplishes sampling speed that this model was unable to attain in BUGS and (b) ceases to reverse scale intermittently. 

Yet, exactly as it does in BUGS, the samples meander.  The parameters are plausible after a hefty sample, but, as I understand, it may require a much larger sample to produce a high enough "effective sample".  Being able to be confident of convergence may be another matter as well.  Maybe the correlation of "z" may not be something that's easily overcome.

As for setting multiple category (NC-1) vectors and sett the first of the series to "0", I have not been successful - but given the volume of code, it may be my own error.  "Warmup" seems to sample correctly, and then when it gets to "sampling", it samples the same values over and over.

So, I'll keep working that option as well intermittently. 

Again, thanks.  All this is not to say that this isn't possible.  By all accounts it should be with some coaxing.

Best,
Jason

Best,
Jason

[Andrew Gelman]

Hi, to check if you're really having convergence, you might consider saving some predictive quantities and checking the convergence of them.  Sometimes it's possible to get good convergence on predictive quantities even if some intermediate parameters are only weakly identified.  Also in your model below perhaps it will help to use a proper prior on theta_raw.

[Bob Carpenter]

Because the simplex is compact, you get a proper uniform prior
on it. Same as for correlation matrices and unit-length vectors.

I'm pretty sure there's a proper prior on the unconstrained space
through the Jacobian. If not, we've (or at least I've) been
giving people bad advice on the list.

Practically speaking, there's no issue at all just sampling
from a simplex with no priors declared in Stan. It mixes
very well.

- Bob

>> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.

[Bob Carpenter]

On 1/24/14, 5:10 PM, Jason Engle wrote:
> Hi Bob (and all),
>
> Just an update. The workaround:
>
> parameters {
> simplex[N] theta_raw;
> real<lower=0> theta_scale;
> ...
> }
> transformed parameters {
> vector[N] theta;
> theta <- theta_scale * (theta_raw - 1.0 / N);
> ...
>
> it works.
>
> That is, it (a) accomplishes sampling speed that this model was unable to attain in BUGS and (b) ceases to reverse scale
> intermittently.
>
> Yet, exactly as it does in BUGS, the samples meander. The parameters are plausible after a hefty sample, but, as I
> understand, it may require a much larger sample to produce a high enough "effective sample". Being able to be confident
> of convergence may be another matter as well. Maybe the correlation of "z" may not be something that's easily overcome.

The only way to compare sampling speeds is to compare
effective samples per second. Time per iteration is
irrelevant.

> As for setting multiple category (NC-1) vectors and sett the first of the series to "0", I have not been successful -
> but given the volume of code, it may be my own error. "Warmup" seems to sample correctly, and then when it gets to
> "sampling", it samples the same values over and over.

That's very strange. Do you recall what the tree depth is?
Or even better, do you have a smallish model that demonstrates
this behavior?

I've gotten this technique (pinning a value to 0) to work before,
so I know it's possible.

- Bob