Identification of multinomial (by means of Poisson) IRT/ideal point model

[Thomas]

I’m fitting a multinomial (by means of Poisson) IRT/ideal point model for category counts in policy documents with

Count(document, category) ~ Poisson(exp(alpha_doc + psi_cat + beta_cat*theta_doc))

As constraints, one psi and one beta are set to zero, and the ideal points theta receive a standard normal prior. A conditional statement based on the order of two betas is used to reverse the scale when necessary (I know that half-normal priors might be better for that, but I had mixed experiences when working with those).

It seems that the constraints are sufficient to identify the model, at least I get very similar results across different runs. It seems to work since the standard normal prior on the ideal points theta apparently leads to a standard normal posterior distribution of the ideal points. Is there a simple answer why this happens? I would not necessarily have expected that, and when I code the same model in Jags, theta does not come out standardized in the posterior and the model lacks identification. Is it ok to rely on the prior doing the job in Stan or would it still be better to standardize the ideal points deterministically (probably at the expense of speed?)? Any hints appreciated.

Code is as follows:

data {
int<lower=1> dir[2]; // doc id pair to set direction
int<lower=1> D; // number of docs
int<lower=1> T; // number of terms/categories
int<lower=1> N; // number of obs = number of docs * number of cats
int<lower=1> doc[N]; // doc id
int<lower=1> term[N]; // cat id
int<lower=0> y[N]; // counts as one long vector
}

parameters {
vector<lower=-20, upper=20>[D] alpha;
vector[T-1] psi_rest; // unconstrained elements of psi
vector[D] theta;
vector[T-1] beta_raw; // unconstrained elements of beta
real<lower=0> tau_psi;
real<lower=0> tau_beta;
}

transformed parameters {
vector[T] psi;
vector[T] beta;
real<lower=0> sigma_psi;
real<lower=0> sigma_beta;
psi[1] <- 0; // first cat. as reference cat.
beta[1] <- 0;
for (t in 2:T){
psi[t] <- psi_rest[t-1];
if (beta_raw[dir[2]-1]-beta_raw[dir[1]-1] > 0) {
beta[t] <- beta_raw[t-1];
}
else {
beta[t] <- -beta_raw[t-1];
}
}
sigma_beta <- 1/sqrt(tau_beta);
sigma_psi <- 1/sqrt(tau_psi);
}

model {
vector[N] logmu;
// alpha ~ uniform(-20,20);  implicit above
psi_rest ~ normal(0,sigma_psi);
beta_raw ~ normal(0,sigma_beta);
theta ~ normal(0,1);
tau_beta ~ gamma(.01,.01);
tau_psi ~ gamma(.01,.01);
for(n in 1:N){
logmu[n] <- alpha[doc[n]] + psi[term[n]] + theta[doc[n]] * beta[term[n]];
}
y ~ poisson_log(logmu) ;
}

[Bob Carpenter]

On Oct 9, 2014, at 2:42 PM, Thomas <daub...@tcd.ie> wrote:

> I’m fitting a multinomial (by means of Poisson) IRT/ideal point model for category counts in policy documents with
>
> Count(document, category) ~ Poisson(exp(alpha_doc + psi_cat + beta_cat*theta_doc))
>
> As constraints, one psi and one beta are set to zero, and the ideal points theta receive a standard normal prior. A conditional statement based on the order of two betas is used to reverse the scale when necessary (I know that half-normal priors might be better for that, but I had mixed experiences when working with those).

I'm curious what the mixed experience was. Did you
constrain the parameters to be positive? The log transform on the
variables does change the posterior geometry and hence affects
sampling.

> It seems that the constraints are sufficient to identify the model, at least I get very similar results across different runs.

That does sound sufficient to identify the model.

> It seems to work since the standard normal prior on the ideal points theta apparently leads to a standard normal posterior distribution of the ideal points. Is there a simple answer why this happens?

Yes. When you have too many degrees of freedom in the parameterization (i.e.,
non-identifiability), what happens is that one parameter can go up
or down in scale or shift in location, with a corresponding shift or scaling
of another parameter group without affecting the posterior
density. By assigning one group of variables a normal(0,1) prior, that
identifies the scale and location that's going to fit best, hence you
get normal(0,1) posteriors for the constrained variable and everything else
floats.

There's a discussion in the Gelman and Hill's book section on (multilevel)
IRT and ideal-point models.

> I would not necessarily have expected that, and when I code the same model in Jags, theta does not come out standardized in the posterior and the model lacks identification.

My guess would be that either (a) it's not exactly the same model in
JAGS, or (b) you haven't run JAGS long enough to converge.
Gibbs sampling is slow to converge if you have correlations in the posterior,
as you do in these essentially unidentified models.

> Is it ok to rely on the prior doing the job in Stan or would it still be better to standardize the ideal points deterministically (probably at the expense of speed?)? Any hints appreciated.

It is critical to rely on the prior or another means of
identification of the model in Stan. Just standardizing the
resulting parameter draws won't work. There's a discussion
of the issue for a very simple example with posterior draws
in the problematic posteriors chapter of the manual, with
more explanation about what's going on and why the post-standardization
seems to work in Gibbs (as far as I know, nobody's shown it gives
the right answer formally, though Andrew's convinced it does), but
not in Hamiltonian Monte Carlo.

- Bob

[Thomas Daeubler]

Thanks very much for the rapid reply, Bob!

Regarding the half-normal priors: Yes, I constrained one beta to be positive in some earlier versions of the model. What sometimes happened was that this beta then was estimated as barely larger than zero (although it was one of the more discriminating ones in the unconstrained version), and the scale was still in the wrong direction.

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

--
Dr. Thomas Däubler
Postdoc Fellow
Mannheim Centre for European Social Research
68131 Mannheim
Germany
Phone: +49 (0)621 181 2809
Email: thomas....@mzes.uni-mannheim.de

[Thomas]

After experimenting a good deal with the above model, I have for now settled on the following specification (see code below):

- pin down one of the term/item parameters each, draw the rest from a bivariate normal, with mean zero, using the Cholesky factor re-parameterization, using a non-hierarchical prior on the covariance matrix (3s on the diagonal and zeros elsewhere)
- keep the soft constraint on the ideal points using the standard normal prior
- resolve the reflection invariance only in the generated quantities block

This works well and the Rhats generally look good after a few hundred iterations (from random inits), and I can replicate the means I get from an ML version of the model. However, there is high autocorrelation in the chains, the effective sample sizes are well below a third of the number of iterations for many parameters, and the warmup seems slow (for a 350x55 cells table, 200 warmup iterations take something between 2 and 4 mins, sampling itself is quicker, that's on a 3.6 Ghz machine with 16 GB RAM).
Ok, I understand that the merely soft constraint on the ideal points leads to shifts in other parameters across draws, but e.g. using an additional hard sum-to-zero constraint on the ideal points only somewhat alleviates the problem (I've also read the post on STZ constraints, where it was argued that hard constraints do not essentially change the model.)

Is there anything wrong with the code, is it just the nature of the model, or is there anything else I can try to do? Does the thinning interval choice actually affect the warmup behavior - I guess not?
Please let me know if you need further information about the model or its behavior.

Thanks for any hints!

Thomas

data {

int<lower=1> D; // number of docs

int<lower=1> T; // number of terms/cats

int<lower=1> N; // number of obs = number of docs * number of cats
int<lower=1> doc[N]; // doc id
int<lower=1> term[N]; // cat id
int<lower=0> y[N]; // counts as one long vector

cov_matrix[2] Sigma_termpar; // prior for co-variance matrix

transformed data {
matrix[2,2] L;
L <- cholesky_decompose(Sigma_termpar);

}

parameters {
vector<lower=-20, upper=20>[D] alpha;

matrix[T-1,2] termpar_e;
vector[D] theta_raw;
}

transformed parameters {
matrix[T,2] termpar_raw;

// constrain first termpar pair to 0; and remaining termpar ~ multi_normal(0,Sigma_termpar)
termpar_raw <- append_row(rep_row_vector(0,2), termpar_e*L);
}

model {

vector[N] logmu;
// note: alpha is uniform
to_vector(termpar_e) ~ normal(0,1);
theta_raw ~ normal(0,1);

for(n in 1:N){
logmu[n] <- alpha[doc[n]] + termpar_raw[term[n],1] + theta_raw[doc[n]] * termpar_raw[term[n],2];
}
y ~ poisson_log(logmu);
}

generated quantities {
vector[D] theta;
matrix[T,2] termpar;
real reflect;

reflect <- if_else(theta_raw[1] > theta_raw[2],-1,1);
for(t in 1:T) {
termpar[t] <- to_row_vector(termpar_raw[t]);
termpar[t,2] <- reflect*termpar[t,2];
}
theta <- reflect*theta_raw;

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Andrew Gelman]

The hard constraint on alpha does not look good to me.  I prefer a prior such as normal.  Also I’m guessing that instead of that “reflect” thing, you could use an ordered data type in Stan and that would work better.

A

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.

[Bob Carpenter]

What's an "ML version of the model"? You mean software
to fit a point estimate?

If you only want a posterior mode, you can use Stan's
optimization functionality.

And you absolutely don't want to do a signum-type
reflection if it leads to different posterior modes. For
Stan to be efficient, the parameters themselves should be
identified.

- Bob

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.
>
>
>

> --
> 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.

[Thomas]

Thank you so much for the suggestions, Andrew and Bob!

(Yes, I used some ML code outside of stan to compare it to my posterior means - but I'd like to get the Bayesian version working.)

I replaced the uniform on alpha with a normal(mu_alpha, sigma_alpha) and that makes it run somewhat faster.

I tried the ordered vector suggestion for the reflection invariance, but that doesn't work since the two constrained ideal points may end up being placed very close together (that's similar to what happens when using two half-normals).

I am aware that the above version does not actually solve the bimodality in the model. Alternatively, I have also tried pinning down two ideal points to fix the direction, but I am struggling with achieving that the constraints bite, i.e. that the remaining ideal points are placed relative to the two fixed ones.
If I fix two positions known to be quite extreme to -1 and 1, and I use a normal(0,5) prior or a hierarchical normal prior (mu_theta,sigma_theta) for the others, then the polarity of the remaining ideal points can still vary across chains: as compensation, the rows with the two constrained thetas just get two alphas that are "off" (when I compare them across chains).  I may be wrong, but I guess another problem in choosing the prior for the unconstrained ideal points is the multiplicative indeterminacy between the ideal points and the discrimination parameters. In my original model, the standard normal prior on the ideal points pretty much fixes their variance and the discrimnations adjust (but the direction is not fixed). When pinning down two ideal points, the direction is fixed, but it is a rather weak constraint on the variance of the ideal points (which still needs to be inferred), and that makes it also harder to infer the variance of the discriminations at the same time.
Any suggestions on what priors to use on the remaining ideal points if I fix two? What looks most promising so far is a normal(0, sigma_theta), but I need to test it further.

Thanks!

[Andrew Gelman]

> On Mar 11, 2015, at 2:13 PM, Thomas <daub...@tcd.ie> wrote:
>
> I tried the ordered vector suggestion for the reflection invariance, but that doesn’t work since the two constrained ideal points may end up being placed very close together (that's similar to what happens when using two half-normals).

The ordering should work even if the two points are close. But if that is a concern (that you don’t want to constrain them to be ordered), you can simply give them different priors with equal variances but with one prior having a mean lower than the other. This should break the aliasing and identify the model. In general we much prefer this sort of identification strategy, (rather than having an unidentified model and then summing to average over aliases).

cc-ing Bob because he likes to collect this sort of advice

[Thomas]

Thank you very much for the suggestions - with your help I've found an elegant solution:

What eventually worked for me is to use an ordered[2] type on two discrimination parameters (whose variance is unconstrained relative to the ideal points with the standard normal prior). And I can confirm that this works well even if applied to two items with weak/similar discriminations. However, the same strategy with the ordered vector type did not work for me when applied to the ideal points. The constraint does not fix the direction, as even ideal points that should be far away from each other can end up at the same location in the posterior (with some tiny difference  to just satisfy the constraint), with the scale being reversed. I also tried the suggestion with the different prior means - but without success.

I have a follow-up question. I am now running a hierarchical multi-group version of the model - what Jean-Paul Fox in his Bayesian IRT Modelling book refers to as a 'random item' model. Item parameters, the mean of the latent variable, and the variance of the latent variable/discriminations are allowed to differ across groups.

I tried various different ways of fixing the direction in this hierarchical model, but without success. For instance, using the ordered vector type on the same two discriminations in each group does not solve the problem. I also experimented with combinations of constraints on the discriminations and two group means of the ideal points. The result is always that, when comparing across chains, the direction is the same in some groups, but not in others. I'd be grateful if you have any recommendations.

Code for the model can be found below (this is a version without any constraints seeking to fix reflection invariance). Except for the sign indeterminacy, the results look ok when comparing posterior means across chains.

Thomas
 
data {
int<lower=1> D; // number of docs
int<lower=1> T; // number of terms/cats

int<lower=1> G; // number of groups

int<lower=1> N; // number of obs = number of docs * number of cats
int<lower=1> doc[N]; // doc id
int<lower=1> term[N]; // cat id

int<lower=1> gr[N]; // group id
int<lower=1> grdoc[D]; // group id at doc level

int<lower=0> y[N]; // counts as one long vector
}

parameters {

real mu_alpha;
real<lower=0> sigma_alpha;
vector[T-1] mu_psi;
vector<lower=0>[T-1] sigma_psi;
vector[G] mu_theta; // group means
vector<lower=0>[G] sigma_beta; // group-specific sd

vector[D] alpha;
matrix[T-1,G] psi_rest;
matrix[G,T] beta;
vector[D] theta;
}

transformed parameters {

matrix[T,G] psi;
psi <- append_row(psi_rest,rep_row_vector(0,G)); // set one difficulty to zero in each group
}


model {
// local variables
vector[N] logmu;
vector[D] thetamean;

mu_alpha ~ normal(0,5);
sigma_alpha ~ cauchy(0,5);
alpha ~ normal(mu_alpha,sigma_alpha);

mu_psi ~ normal(0,5);
sigma_psi ~ cauchy(0,5);
for(t in 2:T){
psi_rest[t-1] ~ normal(mu_psi[t-1],sigma_psi[t-1]);
}

sigma_beta ~ cauchy(0,5); // vector - different for each group
for(g in 1:G){
beta[g] ~ normal(0,sigma_beta[g]);
}

mu_theta ~ normal(0,1); // means of latent variable in group
for(d in 1:D){
thetamean[d] <- mu_theta[grdoc[d]];
}
theta ~ normal(thetamean,1); // (group-level sd is one as well)

for(n in 1:N){
logmu[n] <- alpha[doc[n]] + psi[term[n],gr[n]] + theta[doc[n]] * beta[gr[n],term[n]];
}
y ~ poisson_log(logmu);
}
'

[Bob Carpenter]

> On Mar 25, 2015, at 4:07 AM, Thomas <daub...@tcd.ie> wrote:
>
> Thank you very much for the suggestions - with your help I've found an elegant solution:
>
> What eventually worked for me is to use an ordered[2] type on two discrimination parameters (whose variance is unconstrained relative to the ideal points with the standard normal prior).

The worry Andrew brought up with this solution is that if you have

ordered[2] a;

then if the posterior mean and sd of a[1] is 0,1 and of a[2] is 1,1, then
there are probably cases where you should have a[2] < a[1], but won't with
the costraint.

I have no idea if this is relevant for your model.

> And I can confirm that this works well even if applied to two items with weak/similar discriminations. However, the same strategy with the ordered vector type did not work for me when applied to the ideal points. The constraint does not fix the direction, as even ideal points that should be far away from each other can end up at the same location in the posterior (with some tiny difference to just satisfy the constraint), with the scale being reversed.

I'm not sure what you mean by scale being reversed. The
posterior uncertainty or some scale parameter?

> I also tried the suggestion with the different prior means - but without success.
>
> I have a follow-up question. I am now running a hierarchical multi-group version of the model - what Jean-Paul Fox in his Bayesian IRT Modelling book refers to as a 'random item' model. Item parameters, the mean of the latent variable, and the variance of the latent variable/discriminations are allowed to differ across groups.

How did the model in the book resolve the multimodalities?

- Bob

[Andrew Gelman]

Yes, as far as I know this is still an open problem.

[Thomas]

I made some changes to the group-partitioned model.

Regarding reflection invariance:
The model in the book did not face this issue, as it was a latent ability type model with a unidirectional latent variable. Since I wasn't happy with any other solution, I now address that issue by initializing from near the posterior means of results from an unpartitioned model.

I switched to drawing the discrimination parameters from a group-level mean, having a reference category, as for the difficulty parameters.

I also use a negative binomial model now to address outliers, and I have a follow-up question concerning that:
I get quite different results depending on the level at which I let the overdispersion parameter (inverse phi) vary. I use either group (in my case country) or item.

I have an external measure I can validate the ideal points with. When I compare the results from the validation exercise across specifications, there seem to be three types of countries:

a) countries where both use of country-level or item-level phi works well (i.e. ideal points are highly correlated with the validation criterion)
b) countries where item-level phi gives much better results than country-level phi. It seems that in these cases
there are outlier cells with large counts not due to ideology and the item-level phi seems good at fixing such measurement errors.
c) countries where country-level phi gives much better results than item-level phi. It seems that for these cases the country-level phi somehow makes it easier that items receive country-specific discrimination parameters that are quite different from their group-level/international means. E.g. what is internationally a fairly neutral policy category may in one country be stongly discriminating, and the model is much better at using the information in those items with country-level phi.

Of course I'd like to have the best of both worlds - any ideas how to get there? Is there a way to let phi draw on country-level and item-level information (ideally without resorting to a model with cell-level cross-classified random effects)?

Letting phi vary at the country-item-combination seems not to be the solution, since the results then are very similar to the version with item-level phi.

Thanks a lot in advance!

Thomas



data {
int<lower=1> D; // number of docs
int<lower=1> T; // number of terms/cats
int<lower=1> G; // number of groups
int<lower=1> N; // number of obs = number of docs * number of cats
int<lower=1> doc[N]; // doc id
int<lower=1> term[N]; // cat id
int<lower=1> gr[N]; // group id
int<lower=1> grdoc[D]; // group id at doc level

// int<lower=1> grsize[G]; // number of docs in each group

int<lower=0> y[N]; // counts as one long vector
}

parameters {

vector<lower=0>[G] sigma_alpha;
vector[G] mu_alpha;
vector[D] alpha;

vector[T-1] mu_psi;
real<lower=0,upper=10> cauchy_scale_psi;
vector<lower=0>[T-1] sigma_psi;
matrix[T-1,G] psi_e;

vector[T-1] mu_beta;
real<lower=0,upper=10> cauchy_scale_beta;
vector<lower=0>[T-1] sigma_beta;
matrix[T-1,G] beta_e;

row_vector[G] mu_theta;
vector[D] theta;

vector<lower=0,upper=100>[T] phi_inv; // term/item-level phi
// vector<lower=0,upper=100>[G] phi_inv; // alternatively: group-level phi

}

transformed parameters {

matrix[T,G] psi;

matrix[T,G] beta;

matrix[T-1,G] psi_rest;
matrix[T-1,G] beta_rest;

for(t in 2:T){
psi_rest[t-1] <- mu_psi[t-1] + psi_e[t-1] * sigma_psi[t-1];
beta_rest[t-1] <- mu_beta[t-1] + beta_e[t-1] * sigma_beta[t-1];
}
psi <- append_row(psi_rest,rep_row_vector(0,G));
beta <- append_row(beta_rest,rep_row_vector(0,G));

}


model {
// local variables
vector[N] logmu;
vector[D] thetamean;

vector[D] alphamean;
vector[D] alphasd;
vector[N] philong;

mu_alpha ~ normal(0,5);
sigma_alpha ~ cauchy(0,5);

mu_psi ~ normal(0,5);
sigma_psi ~ cauchy(0,cauchy_scale_psi);
to_vector(psi_e) ~ normal(0,1);

mu_beta ~ normal(0,5);
sigma_beta ~ cauchy(0,cauchy_scale_beta);
to_vector(beta_e) ~ normal(0,1);

mu_theta ~ normal(0,1);

for(d in 1:D){
thetamean[d] <- mu_theta[grdoc[d]];

alphamean[d] <- mu_alpha[grdoc[d]];
alphasd[d] <- sigma_alpha[grdoc[d]];
}
theta ~ normal(thetamean,1);
alpha ~ normal(alphamean,alphasd);

for(n in 1:N){
logmu[n] <- alpha[doc[n]] + psi[term[n],gr[n]] + theta[doc[n]] * beta[term[n],gr[n]];
philong[n] <- 1/phi_inv[term[n]];
// philong[n] <- 1/phi_inv[gr[n]];
}
y ~ neg_binomial_2_log(logmu,philong);