questions around hierarchical ordered logit

[Tobias Konitzer]

Hi All,

I have been trying to build a hierarchical ordered logit (below), and am hoping for feedback. This is in the framework for MRP (multi-level modeling and post-stratification to adjust for non-response bias in convenient samples), so I am trying to put out the linear predictor per row in my population-level data set. Please ignore the empty loops in the model, they stem from a script automating the model itself (I hope it does not unduly influence stan).

#1) does the model look fine to you?

#2) For confidence intervals, I currently compute three dimensional arrays and assess uncertainty from there by simulating N_pops(100)*N_category distributions over all responses.  

 

holder<-array(NA, c(X[['N_pops']],data_structure$responses$n_levels,100))
  for (p in 1:data_structure$responses$n_levels){
    if (p==1)
    {holder[,p,]<-(1/(1+exp(params[["eta"]][,w,]-params[[paste0("c_",data_structure$responses$response_name)]][,p])))}
    else if (p>1 & p<data_structure$responses$n_levels)
    {holder[,p,]<-(1/(1+exp(params[["eta"]][,w,]-params[[paste0("c_",data_structure$responses$response_name)]][,p])))-
      (1/(1+exp(params[["eta"]][,w,]-params[[paste0("c_",data_structure$responses$response_name)]][,p-1])))}
    else {holder[,p,]<-(1 - (1/(1+exp( params[["eta"]][,w,]- params[[paste0("c_",data_structure$responses$response_name)]][,p-1] ))))}
  }

for( j in 1:data_structure$responses$n_level){

        y.pred[,j] <- sapply(1:100, function(x) sum( input_data$pops_N * holder[,j,][,x] ) / sum( input_data$pops_N ))

      }

I then awkwardly assess confidence by 

for( j in 1:data_structure$responses$n_level){
        y.pred[,j] <- sapply(1:100, function(x) sum( input_data$pops_N * holder[,j,][,x] ) / sum( input_data$pops_N ))
      }
      #adjust cis
      pred_temp_national<-rbind(apply(y.pred,2,quantile,prob=0.025), y.pred[which.max(params$lp__),],apply(y.pred,2,quantile,prob=0.975))

 

Is there a better way? I really have a hard time of thinking about this in a >2 choice setting.

#3) I would like to make credible intervals comparable to margin of error (currently, uncertainty bands in my model are much wider). Is it fair to scale somewhere by sqrt(N) or something akin?

Many thanks,

Tobi

[Tobias Konitzer]

And, I would like to ask one follow up. Do globally fixed cut-points produce some weird constraints here, and should they always vary by group? In any level with many groups, that can get very untractable, no? My worry is that globally fixed cut-points fixes the rank-order of class probabilities by group, if I understand correctly...

Thanks again!

[Bob Carpenter]

We don't have good prior recommendations for these as far as I know. I think global cut points are more common than varying ones, but we're big fans of hierarchical modeling, so if you can figure it out, that'd be great. I didn't see any hierarchical structure in your model or a prior for c_vote.

Empty loops won't do anything---the C++ compiler is usually good enough to ignore all this.

For standard errors on the posterior means, we divide the posterior standard deviation by the sqrt of the number of effective samples. It's all estimated.

- Bob

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Tobias Konitzer]

Thanks for the last reply, Bob. I have been torturing myself to getting to this, but that makes sense. My intuition is that hierarchies on the cut points should be dependent on groups for which you would want to enforce rank in class probabilities to be different. For example, if we want to predict classes globally, the non-hierachical model would fits best. If we want to predict classes for every party, we would need to let cut-points vary by party... Seem correct?

Tobi

[Tobias Konitzer]

Hi All,


One quick update: Any reason against the attached hierarchical ordered logit? It assumes uniform priors on the cutpoints, which I am fine with. With sample data, it converges currently at Max Rhat = 1.094, but this goes drastically up if we increase the number of groups...

Thanks again for all the awesome work you guys are doing. I am (and have been) a big fan of stan, and trying to promote it heavily at my (and others) research team at Facebook.

Best,
Tobi

[Andrew Gelman]

Hi, Tobias, my only thought is on this:

approval[i] ~ ordered_logistic(eta, c_approval[group[i]]);

Is the ordered_logistic distribution vectorized?  If so, I think the program will run faster to take this statement outside the loop.

A

<ordinal_dynamic_generic_custom_dynamics_now_temp.stan>

[Tobias Konitzer]

HI Andrew,

thanks for the quick answer! But, I can't see a way of indexing the cut-points by group-membership of individual i without the more costly loop. Maybe I am missing something obvious...?

Thanks again,

Tobi

[Andrew Gelman]

You can define the argument to the distribution in a loop, then put the ~ statement outside the loop.

[Tobias Konitzer]

Of course! Really missed that; thanks!

Tobi

[Bob Carpenter]

No --- we need to vectorize it. The code Sean's
working on now should make that much easier.

- Bob

[Bob Carpenter]

I also realized we did say something on our priors page. I should look
at our own Wiki:

https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations


• Prior for cutpoints in ordered logit or probit regression
• For full Bayes, uniform priors typically should be ok, I think.
• For modal estimation, put in some pseudodata in each category to prevent "cutpoint collapse."
• Cutpoints are ordered (by definition). Put the prior on the differences between the cutpoints rather than the cutpoints themselves.
• Alternatively, put a prior on the cutpoints and partially pool them, not to a constant, but to a linear function. That is, if the cutpoints are c[1],...,c[K], don't do c[k] ~ normal(mu, sigma); instead do c[k] ~ normal(a + b*k, sigma), estimating the hyperparameters a and b from data. It will probably make sense to put informative priors on a, b, and sigma too.

- Bob

[Andrew Gelman]

Hey, I think I wrote that!

[Tobias Konitzer]

Thanks for all the answers so far! I never liked putting group hierarchies on the cut-points, that never seemed logical to me. So I have gone back and put the hierarchies on the linear predictor, roughly modeled after page 143 in stan reference 2.14 and gelman/hill 2007, chapter 13 (model attached).

The idea is the following: We have a group-specific vector of random effects, that we draw from a multivariate normal, modeling the across-group-within-random-effects variations and correlations in the Varcov matrix of this multivariate normal. Then, we draw class of outcome from an ordered logit where the combination of group-specific random effects is the linear predictor.

But, the model does not converge at all. Is there something obvious that I am missing?

Thanks again,

Tobi

[Daniel Lee]

Hi Tobias,

When you say it's not converging at all, it actually helps to give more info:

- does it run quickly?

- does it not initialize?

- are the rhats ridiculous?

For me, convergence has a specific meaning and it's hard to tell if you also mean the same thing. And there are ways of dealing with it, depending on how you got there. 

So... more info if you want help. We're happy to help if we know how. 

Daniel

[Tobias Konitzer]

Hi Daniel,

Sorry for the lack of info provided. I was hoping there was something obviously wrong in the model. The rhats are all over the place. The max rhat is 812.51, and the model runs very fast. I am thinking there is something in modeling the Covar of within/across group-random effects where things go awry, but I am not sure. Any help is greatly appreciated, I have been relying on this forum multiple times, and the work you guys are doing is unbelievable.

Many thanks,

Tobi

[Bob Carpenter]

As always, I'd suggest starting with a simpler model. Lots
can go wrong in those ordinal models with cutpoints collapsing.
You don't have any priors to keep them from doing that, and you're
also letting them run between -infinity and infinity. There's
often a translational non-identifiability in cutpoints when you
do that if you have an intercept in your model.

Did your program get garbled in transmission? I'm seeing things
like this in it:

for (k in 1:(n_waves-1)) {

}

and then very oddly spaced code that I can't believe anyone intended:

for(i in 1:N){
real eta;

eta =

beta_stat_race[race[i]]+
beta_stat_age[age[i]]+
beta_stat_education[education[i]]+
beta_stat_gender[gender[i]]
;

approval[i] ~ ordered_logistic(eta, c_approval[group[i]]);
}

It might help to organize the code a bit better.

You can define

vector[N] eta = beta_stat_race[race] + beta_stat_age[age] + ...;

and then use eta[n] in the loop. We haven't gotten around to
vectorizing ordered_logistic yet.

You can declare/define now all at once so this:

transformed data {
real scale_const;
real scale_dyn;

scale_dyn =0.05;
scale_const =5;
}

can be just

transformed data {
real scale_const = 5;
real scale_dyn = 0.05;
}

- Bob

[Tobias Konitzer]

Hi Bob,

Thanks for the comment, and a very late apologies for the code. It stems from a text-parsing script, hence the odd spacing etc. but I should have checked before submitting! 

I am still not getting a hierarchical ordered logit to converge. Again, the model is roughly modeled after page 143 in stan reference 2.14 and gelman/hill 2007, chapter 13. The idea is the following: We have a group-level vector of means for every random effect (i.e. a length-7 vector for 7 groups, and then draw a random-effect-specific (say for age_bracket=1) vector of group effects from a multivariate normal with zeros in the off-diagonals to reduce the number of parameters to be estimated . Then, we draw class of outcome from an ordered logit where the combination of group-specific random effects is the linear predictor, and every group has a separate ordered vector of cutpoints. I tried a number of fixes ("break" meaning bad convergence, as in R Hat):

- change the order of hierarchies: draw a global vector of random effects (i.e. a length-4 vector for 4 age-brackets , and then draw a group-specific vector of random effects from a multivariate normal 

- let the model try to estimate off-diagonals of the correlation matrix of the multivariate priors (breaks)

- run simple models with one class of random effects only (breaks)

- run as Bernoulli to avoid "cut-point collapse" (breaks)

I must be missing something obvious here! The cleaned up model is attached (I am not concerned about the efficiency of the loop right now). 

Any thoughts?

Thanks again,

Tobi