IRT model help

[Chris Weber]

Hello,

I am fitting a two parameter logistic IRT model on a relatively large dataset of congressional roll call votes (about 4000). In the data, I also include respondents from a public opinion survey, who have answered questions about these votes. This results in a relatively large dataset (about 56K respondents and 4000 roll call votes), but also large amount of missing data (survey respondents only answer questions about 9 of the votes; senators don't vote on house measures, etc.). I'm trying to obtain the posterior ability distributions for each respondent. The problem is that the model is incredibly slow. I've let it run for a day, with very little progress. Is there anything I can do to speed things up? Am I asking too much by analyzing such a large dataset?

Any advice would be greatly appreciated. Thanks in advance.

Chris 

[Bob Carpenter]

That amount of data's totally manageable. It looks like you're not
doing anything with the missingness yet --- that's a separate issue.

When you say "very little progress", what happens when you run it?

I'm not that familiar with the ideal-point model applications;
in standard IRT models, we can usually safely assume alpha[k] > 0.
My understanding is that they can be hard to identify in terms
of sign of alpha.

What is the mean of y? If it's not close to 0, you may want to
add a prior mean for beta (centered parameterization) or an
intercept for the model as a whole (same model, non-centered
parameterization) --- it's a single scalar unconstrained parameter
either way.

One thing that might be a little more efficient is to define


eta[n] <- alpha[kk[n]]*theta[jj[n]]- beta[kk[n]];

and then use

y ~ bernoulli_logit(eta);

instead of the loop. But that'll only give you a small speedup if
it's even noticeable.

- Bob

P.S. I'm thinking I should just append this to all my messages, but
I'll say it again --- we strongly recommend indenting so programs
are easier to read.

> --
> 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 roll call.R><final.data.Rdata>

[Chris Weber]

I apologize -- I've indented the code, which is also attached. You're correct, I haven't yet tried to tackle the missing data issue. 

Regarding the alpha parameter, that sounds correct. Yet, in this application it probably shouldn't be constrained to be greater than 0, as some votes are more appealing to conservatives, others to liberals -- here, the latent variable is ideology. There are far too many roll call votes to manually recode the responses, so that 1 is a conservative (or liberal) response.  This is the reason why I think alpha[k] should be allowed to be either positive or negative. 

I've added the model intercept (and taking y out of the loop didn't do much), but I still haven't seen much improvement in speed. By "very little progress" I mean 100 iterations in 12 hours or so. After the first few iterations I get the "informational message" 

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Exception thrown at line 22:

stan::math::normal_log: Scale parameter is 0, but must be > 0!

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

which I often get when I run IRT models in stan, but other than that, I'm not sure what/if there is a problem.  Any other suggestions on what else I might change? Thanks in advance.

Chris

[Bob Carpenter]

> On Aug 13, 2015, at 11:51 PM, Chris Weber <christoph...@gmail.com> wrote:
>
> I apologize -- I've indented the code, which is also attached. You're correct, I haven't yet tried to tackle the missing data issue.
>
> Regarding the alpha parameter, that sounds correct. Yet, in this application it probably shouldn't be constrained to be greater than 0, as some votes are more appealing to conservatives, others to liberals -- here, the latent variable is ideology. There are far too many roll call votes to manually recode the responses, so that 1 is a conservative (or liberal) response. This is the reason why I think alpha[k] should be allowed to be either positive or negative.

Right --- that's typical in ideal-point models. I was just trying to say
that I don't have much experience with them, so can't give you good coding
advice. The problem is that you get a sign invariance problem in that you can swap
all the signs and get the same answer.

> I've added the model intercept (and taking y out of the loop didn't do much), but I still haven't seen much improvement in speed. By "very little progress" I mean 100 iterations in 12 hours or so. After the first few iterations I get the "informational message"
>
> Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
> Exception thrown at line 22:
> stan::math::normal_log: Scale parameter is 0, but must be > 0!
> If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
> but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
>
> which I often get when I run IRT models in stan, but other than that, I'm not sure what/if there is a problem. Any other suggestions on what else I might change? Thanks in advance.

No idea offhand. Our IRT models with sigma > 0 used for education
applications typically run much much faster than this. The problem
is that I don't know the effect of allowing alpha to be positive or
negative.

I would suggest trying to fit simulated data just to make sure
your program is working.

You can start with a smaller stepsize (say 0.1 or 0.01) and
a higher target acceptance rate (adapt_delta of 0.9 or even 0.99),
which often gets rid of the numerical instability.

You can also try reasonable initializations, but I don't know how
to do that if you don't know the signs of the deltas.

- Bob

[Gregory Eady]

As I understand it, reflection invariance can be solved by constraining the sign of one of the betas, alphas, or thetas (see Bafumi, Gelman, Park, and Kaplan 2005, 176-179). Probably the most practical solution here is to constrain the sign on an alpha for a bill that you know to be conservative (or liberal). This takes a tiny bit of work in Stan, but it's straightforward enough.

I've simulated a binary item-response model (hopefully this was done correctly), and then coded up the item-response model where the sign on one of the alphas is constrained. Because this is simulated data, I took as the "known" conservative bill the one with the minimum alpha from the simulation. This would have to be set manually in the real-data case, and would be the bill that is most clearly a conservative vote in the data.

One of the problems I see after running this, however, is that some chains get stuck sampling around the "wrong" mode (which has a lower log posterior due to constraining the sign of alpha as it should), but I'm not sure how this can be prevented.

Bob, do you have any suggestions on this?

library(rstan)

library(reshape2) # install.packages("reshape2")

library(shinystan)

n_items <- 30 # Items

n_respondents <- 1000 # Respondents

n_obs <- n_items * n_respondents # Total number of observations

theta  <- rnorm(n_respondents, 0, 1) # Abilities

alpha <- rnorm(n_items, 0, 1) # Discrimination

beta <- rnorm(n_items, 0, 1) # Difficulty

slope.ability <- outer(theta, alpha) # Multiply slope vector by ability vector

intercept <- matrix(rep(beta, n_respondents), nrow = n_respondents, byrow = TRUE)

prob <- plogis(slope.ability - intercept)

D <- data.frame(t(apply(prob, 1, function(x) rbinom(length(x), 1, x))))

names(D) <- 1:ncol(D)

D <- data.frame(respondent = 1:nrow(D), D)

D <- melt(D, id = "respondent", value.name = "y")

D$issue <- as.numeric(gsub("X", "", D$variable))

# SET CONSTRAINT ON ALPHA

constraint <- rep(0, n_items)

# Here I just use the minimum alpha as the one I constrain,

# but any, say, clear conservative roll call could be used here.

constraint[alpha == min(alpha)] <- 1

constraint[(which(constraint == 1) + 1):length(constraint)] <- 2

# This constraint vector will read something like 0, 0, 0, 0, 1, 2, 2, 2, 2, 2 ... 2,

# where 1 indicates the constrained alpha, and the 2s are for bookkeeping in the

# transformed parameters block as you can see below

model_data <- list(n_items = length(unique(D$issue)),

                   n_respondents = length(unique(D$respondent)),

                   n_obs = length(D$y),

                   issue = D$issue,

                   respondent = D$respondent,

                   y = D$y,

                   constraint = constraint)

model_code <- "data {

                    int n_items;       // number of items

                    int n_respondents; // number of persons

                    int n_obs;         // number of responses

                    int respondent[n_obs];   // respondent ID

                    int issue[n_obs];        // issue ID

                    int y[n_obs];            // respondent response

                    int constraint[n_items]; // constraint on difficulty

                } 

               parameters {

                    vector[n_respondents] theta; // ability

                    vector[n_items-1] alpha_raw;           // discrimination (-1 because excludes constrained alpha)

                    real<upper = 0> alpha_raw_constrained; // constrain sign of the assigned alpha 

                    real<lower = 0> mu_alpha;              // discrimination mean

                    real<lower = 0> sigma_alpha;           // discrimination scale

                    vector[n_items] beta;       // difficulty

                    real mu_beta;               // difficulty mean

                    real<lower = 0> sigma_beta; // difficulty scale

                } 

                transformed parameters {

                    vector[n_items] alpha;

                    for(i in 1:n_items) {

                        if(constraint[i] == 0) alpha[i] <- alpha_raw[i];

                        if(constraint[i] == 1) alpha[i] <- alpha_raw_constrained;

                        if(constraint[i] == 2) alpha[i] <- alpha_raw[i-1];

                    }

                }

                model {

                    // Hyperparameters

                    mu_alpha ~ cauchy(0, 3);

                    sigma_alpha ~ cauchy(0, 3);

                    mu_beta ~ cauchy(0, 3);

                    sigma_beta ~ cauchy(0, 3);

                    // Priors

                    alpha_raw ~ normal(mu_alpha, sigma_alpha);  

                    beta ~ normal(mu_beta, sigma_beta);

                    theta ~ normal(0, 1);

                    // Likelihood

                    for(i in 1:n_obs) {

                       y[i] ~ bernoulli_logit(alpha[issue[i]] * theta[respondent[i]] - beta[issue[i]]); 

                    }

                }

                "

# Compile to C

model_irt <- stan_model(model_code = model_code)

# Sample from posterior

options(mc.cores = 6)

# Sample from posterior

posterior_irt <- sampling(model_irt, data = model_data, warmup = 500, iter = 1000,

                          open_progress = FALSE, chains = 6)

launch_shinystan(posterior_irt)

# Plot estimates versus true values of theta (ability)

cor(sapply(extract(posterior_irt, paste0("theta[", 1:n_respondents, "]")), mean), theta)

plot(sapply(extract(posterior_irt, paste0("theta[", 1:n_respondents, "]")), mean), theta)

# Plot estimates versus true values of alpha (discrimination)

cor(sapply(extract(posterior_irt, paste0("alpha[", 1:n_items, "]")), mean), alpha)

plot(sapply(extract(posterior_irt, paste0("alpha[", 1:n_items, "]")), mean), alpha)

# Plot estimates versus true values of beta (difficulty)

cor(sapply(extract(posterior_irt, paste0("beta[", 1:n_items, "]")), mean), beta)

plot(sapply(extract(posterior_irt, paste0("beta[", 1:n_items, "]")), mean), beta)

[Bob Carpenter]

Multimodality is still an open research issue for both optimization
and MCMC. So no, I don't have a general solution. In fact, to have
a proper sample from the posterior, you need to be able to visit all
the modes. Of course, if one is much lower mass than the others (note
this is not lower density at its local optimum, but lower mass in measure
theoretic terms), then there's usually not too much bias in estimation from
ignoring it. IF that's the case, then you can initialize around the good
mode or just toss out the chains that are from the very-low-mass mode.
But overall, you need to proceed with extreme caution if the posterior
is multimodal.

Also, this is why we don't believe the people who don't want to run
multiple chains to diagnose convergence because they think one long chain
is enough, because in theory it will explore all the modes. Another case
where something works in theory, but not in practice.

- Bob

> On Aug 14, 2015, at 6:54 PM, Gregory Eady <gregor...@gmail.com> wrote:
>
> As I understand it, reflection invariance can be solved by constraining the sign of one of the betas, alphas, or thetas (see Bafumi, Gelman, Park, and Kaplan 2005, 176-179). Probably the most practical solution here is to constrain the sign on an alpha for a bill that you know to be conservative (or liberal). This takes a tiny bit of work in Stan, but it's straightforward enough.
>
> I've simulated a binary item-response model (hopefully this was done correctly), and then coded up the item-response model where the sign on one of the alphas is constrained. Because this is simulated data, I took as the "known" conservative bill the one with the minimum alpha from the simulation. This would have to be set manually in the real-data case, and would be the bill that is most clearly a conservative vote in the data.
>
> One of the problems I see after running this, however, is that some chains get stuck sampling around the "wrong" mode (which has a lower log posterior due to constraining the sign of alpha as it should), but I'm not sure how this can be prevented.
>
> Bob, do you have any suggestions on this?
>
> library(rstan)
> library(reshape2) # install.packages("reshape2")
>

> n_items <- 30 # Items
> n_respondents <- 1000 # Respondents
> n_obs <- n_items * n_respondents # Total number of observations
>
> theta <- rnorm(n_respondents, 0, 1) # Abilities
> alpha <- rnorm(n_items, 0, 1) # Discrimination
> beta <- rnorm(n_items, 0, 1) # Difficulty
>
> slope.ability <- outer(theta, alpha) # Multiply slope vector by ability vector
> intercept <- matrix(rep(beta, n_respondents), nrow = n_respondents, byrow = TRUE)
> prob <- plogis(slope.ability - intercept)
>
> D <- data.frame(t(apply(prob, 1, function(x) rbinom(length(x), 1, x))))
> names(D) <- 1:ncol(D)
> D <- data.frame(respondent = 1:nrow(D), D)
>
> D <- melt(D, id = "respondent", value.name = "y")
> D$issue <- as.numeric(gsub("X", "", D$variable))
>
>

> # SET CONSTRAIN ON ALPHA

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

> <Convergence diagnostics.png>