Hierarchical Dynamic Ordered Probit Regression

[Tobias Konitzer]

Hi All,

I have been trying to write a Hierarchical Dynamic Ordered Probit Regression. The model features individual- and state-level predictors (all random effects), and time-interactions to add dynamics. Time interactions are drawn from normal priors sharing a variance-term, and all other predictors are drawn from normal priors sharing a variance term. The model does not converge at all, and I assume it has to do with a lack of identifiability because cut-points float and get in the way of the random effects.

I would appreciate any hint to make that model identifiable!

Thanks,

Tobi

data {

  int<lower=0> N;                                         // sample size

  int<lower=2> J;                                         // number of categories

  

  int<lower=0> n_age;                                     // number of age categories

  int<lower=0> n_race;                                    // number of race categories

  int<lower=0> n_edu;                                     // number of education catoegories

  int<lower=0> n_division;                                // number of Divisions

  int<lower=0> n_states;                                  // number of States

  int<lower=0> n_party_id;                                // number of party-categories

  int<lower=0> n_waves;                                   // number of Waves

 

  int<lower=1,upper=n_age> age[N];                        // age for each respondent

  int<lower=1,upper=n_race> race[N];                      // race for each respondent

  int<lower=1,upper=n_edu> edu[N];                        // education for each respondent

  vector<lower=0,upper=1>[N] female;                      // female for each respondent

  int<lower=1,upper=n_party_id> party_id[N];              // party id for each respondent

  int<lower=0,upper=n_states> state[N];                   // state for each respondent

  int<lower=0,upper=n_waves*2> wave_sex[N];               // wave_sex each respondent

  int<lower=0,upper=n_waves*n_age> wave_age[N];           // wave_age for each respondent

  int<lower=0,upper=n_waves*n_edu> wave_edu[N];           // wave_edu for each respondent

  int<lower=0,upper=n_waves*n_race> wave_race[N];         // wave_race for each respondent

  int<lower=0,upper=n_waves*n_party_id> wave_party_id[N]; // wave_party_id for each respondent

 //state level effects

  int<lower=1,upper=n_division> division[n_states];       // division for each state

  real<lower=0,upper=100> vote_prev_state[n_states];      // 2012 vote at the state-level

  int<lower=1,upper=J> y[N];                              // vote-choice

                                                            

}

parameters {

 real beta;                                               // common intercept

  real sex;                                               // female slope                    

  vector[n_age] a_age;                                    // age random effects

  vector[n_edu] a_edu;                                    // education random effects

  vector[n_race] a_race;                                  // race random effects

  vector[n_party_id] a_party_id;                          // party random effects

  vector[n_waves*2] a_wave_sex;                           // change for female                    

  vector[n_waves*n_age] a_wave_age;                       // change for sex

  vector[n_waves*n_edu] a_wave_edu;                        // change for edu

  vector[n_waves*n_race] a_wave_race;                      // change for race

  vector[n_waves*n_party_id] a_wave_party_id;              // change for party

                                                          //state level effects

  vector[n_division] a_division;                          // division random effects

  real gamma;                                             // previous vote- slope

  vector [n_states] a_state;                              // state effects

  ordered[J-1] tau;                                       // thresholds

  

  

                                                          // sigmas

  real<lower=0,upper=100> sigma;                          // sd for base-line parameters  

  real<lower=0,upper=100> sigma_change;                   // sd for dynamic parameters  

}

transformed parameters{                                   // state-level effects transformed

  real a_state_hat [n_states];

  for (s in 1:n_states)

a_state_hat[s]<-a_division[division[s]]+gamma*vote_prev_state[s];

}

model {

  vector[J] theta;

sigma ~ cauchy(0,2.5);                                    // half-cauchy priors, p. 54 for reference

sigma_change ~ cauchy(0,2.5);

a_age ~ normal (0, sigma);                                //base-line parameters

a_edu ~ normal (0, sigma);

a_race ~ normal (0, sigma);

a_party_id ~ normal (0, sigma);

beta ~ normal(0, sigma);

sex ~ normal(0, sigma);

gamma ~ normal(0, sigma);

a_division ~ normal (0, sigma);

a_wave_sex~ normal (0, sigma_change);                   //change parameters

a_wave_age~ normal (0, sigma_change);

a_wave_edu~ normal (0, sigma_change);

a_wave_race~ normal (0, sigma_change);

a_wave_party_id~ normal (0, sigma_change);

                                                       //#1st Level

a_state~ normal(a_state_hat,sigma);

for (i in 1:N) {                                      //Log Likelihood

real y_hat;

y_hat <- beta + sex * female[i]+a_state[state[i]]+

a_party_id[party_id[i]] + a_age[age[i]]+a_edu[edu[i]] + a_race[race[i]]+

+a_wave_sex[wave_sex[i]]+a_wave_age[wave_age[i]]+a_wave_edu[wave_edu[i]]+a_wave_race[wave_race[i]]+a_wave_party_id[wave_party_id[i]];

theta[1] <- 1 - Phi(y_hat - tau[1]);

for (k in 2:(J-1))

theta[k] <- Phi(y_hat - tau[k-1]) - Phi(y_hat - tau[k]);

theta[J] <- Phi(y_hat - tau[J-1]);

y[i] ~ categorical(theta);

} }

[Gregory]

I assume this has to do with the fact that you have an intercept term in there. To identify an ordinal model, you need to fix a cut point, or, alternatively, just drop the intercept.

[Tobias Konitzer]

I have tried leaving out the generic, non-random intercept with no avail. I have also tried fixing a cut-point at 0, with no avail. The mixing still goes completely haywire, i.e. does not move from the initial value.

Tobi

[Andrew Gelman]

Hi, to channel Andrew, let me quickly say:

No no no on real<lower=0,upper=100> sigma; and real<lower=0,upper=100> sigma_change;

It should just be real<lower=0> in each case.  Then you can use a normal(0,2.5) prior or whatever.

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

[Michael Betancourt]

Wow, you sound just like Andy!

[Bob Carpenter]

a_state_hat is essentially an intercept.

- Bob

[Tobias Konitzer]

Thanks for the comment! In my mind, all random effects (so essentially every parameter in my model) would get in the way of the cut-points. I am surprised that fixing one of the J-1 cut-points (and getting rid of that constraint on the variance parameters, thanks for the hind!) does not take care of this. Is there an intuitive solution to this? The combination of cut points and random effects, I think, confuses me (and the model). Again, many many thanks for maintaining this wonderful group!

  ordered[J-2] tau_raw;                                   // thresholds

  

  

                                                          // sigmas

  real<lower=0> sigma;                                    // sd for base-line parameters  

  real<lower=0> sigma_change;                             // sd for dynamic parameters  

}

transformed parameters{                                   // state-level effects transformed

  real a_state_hat [n_states];

    vector[J-1] tau;

  for (s in 1:n_states)

a_state_hat[s]<-a_division[division[s]]+gamma*vote_prev_state[s];

  tau[1]<-0;

  for (k in 2:(J-1)) 

    tau[k] <- tau_raw[k-1];

[Bob Carpenter]

You could also try doing the same thing with the other
varying intercepts ("random effects") as with sex, that
is use only (A - 1) of the A possibilities and set the other
one to 0. The models aren't identified other than through
priors otherwise.

Also, I'm not sure why you use the same hierarchical scales
sigma and sigma_change everywhere. Don't you expect the scale
of effects to vary?

Given that you're setting the first cutpoint to 0, you
probably want the rest of them to be positive,
so rather than

ordered[J-2] tau_raw;

you probably want

positive_ordered[J - 2] tau;

which ensures all tau[j] > 0.

I just defined tau directly. We won't need tau_raw, but can
just define the cumulative probabilty terms so that the calls
to Phi aren't repeated:

vector[J] phi_cumulative;
phi_cumulative[j] <- Phi(y_hat); // tau[1] = 0 by construction
for (j in 2:(J - 1))
phi_cumulative[j] <- Phi(y_hat - tau[j - 1]);
phi_cumulative[J] <- 1;

Then I'd use those to define theta:

theta[1] <- phi_cumulative[1];
for (j in 2:J)
theta[j] <- phi_cumulative[j] - phi_cumulative[j - 1];

It'll be much more efficient without the repeated calls to Phi
and there's some savings from not having to write out two versions
of tau. You could shave a little more off by defining theta[J] directly
and making phi_cumulative size J - 1 in stead of J, but that probably won't
be noticeable.

I'd also suggest writing the linear predictor as

y_hat
<- beta
+ sex * female[i]

+ a_state[state[i]]

+ a_party_id[party_id[i]]
+ a_age[age[i]]

+ a_edu[edu[i]]
+ a_race[race[i]]
+ a_wave_sex[wave_sex[i]]
+ a_wave_age[wave_age[i]]
+ a_wave_edu[wave_edu[i]]
+ a_wave_race[wave_race[i]]
+ a_wave_party_id[wave_party_id[i]];

I find that much easier to read than

> y_hat <- beta + sex * female[i]+a_state[state[i]]+
> a_party_id[party_id[i]] + a_age[age[i]]+a_edu[edu[i]] + a_race[race[i]]+
> +a_wave_sex[wave_sex[i]]+a_wave_age[wave_age[i]]+a_wave_edu[wave_edu[i]]+a_wave_race[wave_race[i]]+a_wave_party_id[wave_party_id[i]];

R and BUGS and JAGS make you write + at the end of a line because
they have parsers that treat end-of-line whitespace differently
than other whitespace, but start of line is traditional in code
and in math typesetting.

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

> <Ordered_dynamic.stan>

[Tobias Konitzer]

Hi All,

Thank you for the answers. Bob comments settles it for me, and the model converges fine now. However, the model produces too much smoothing (i.e. the wave_random_effect parameters capture no variation), so I would like to put more "natural" dynamics into the model. 

I thought I could simply put a random-walk-prior on each of the baseline coefficients including an evolution (variance) term, i.e. a_age ~N(a_age_(t-1), sigma_evolution), but now again the model does not converge (for now, I implemented only a version with a binary DV). Any idea what is going on here? Thanks again for these great answers, and fantastic maintenance of this site!

Tobi

data {

  int<lower=0> N;                                         // sample wave_size

  //int<lower=2> J;                                         // number of categories

  

  int<lower=0> n_age;                                     // number of age categories

  int<lower=0> n_race;                                    // number of race categories

  int<lower=0> n_edu;                                     // number of education catoegories

  int<lower=0> n_division;                                // number of Divisions

  int<lower=0> n_states;                                  // number of States

  int<lower=0> n_party_id;                                // number of party-categories

  int<lower=0> n_waves;                                   // number of Waves

  

  

  

  int<lower=0,upper=n_states> state[N];                   // state for each respondent

  int<lower=0,upper=1> gender[N];                         // gender each respondent

  int<lower=0,upper=n_age> age[N];                        // age for each respondent

  int<lower=0,upper=n_edu> edu[N];                        // edu for each respondent

  int<lower=0,upper=n_race> race[N];                      // race for each respondent

  int<lower=0,upper=n_party_id> party_id[N];              // party_id for each respondent

  //state level effects

  int<lower=1,upper=n_division> division[n_states];       // division for each state

  real<lower=0,upper=100> vote_prev_state[n_states];      // 2012 vote at the state-level

  

  int<lower=0,upper=1> y[N];                              // vote-choice

  int<lower=0,upper=N> in_wave[N];

  int<lower=0>        wave_size[n_waves];      

 int<lower=0> wave_lims[n_waves,2];

}

parameters {

  

  real beta [n_waves];                                  // coefficient for female; wave                    

  vector[n_age] a_age [n_waves];                        // coefficient for age; wave 

  vector[n_edu] a_edu [n_waves];                        // coefficient for edu; wave 

  vector[n_race] a_race [n_waves];                      // coefficient for race; wave 

  vector[n_party_id] a_party_id [n_waves];              // coefficient for party_id; wave 

  

  //state level effects

  vector[n_division] a_division;                        // division random effects

  real gamma;                                           // previous vote- slope

  vector [n_states] a_state;                            // state effects

  

  // sigmas

  real<lower=0> sigma;                                  // sd for base-line parameters  

  real<lower=0> sigma_evolution;                        // "evolution" term for dynamic parameters  

}

transformed parameters{                                   // state-level effects transformed

  real a_state_hat [n_states];

  for (s in 1:n_states)

    a_state_hat[s]<-a_division[division[s]]+gamma*vote_prev_state[s];

}

model {

  

  sigma_evolution ~ cauchy(0,2.5);                      // half-cauchy priors, p. 54 for reference

  sigma ~ cauchy(0,2.5);

  

  

  //#1st Level

    a_state~ normal(a_state_hat,sigma);

  //evolution

  for (k in 2:n_waves) {

    beta[k] ~ normal(beta[k-1], sigma_evolution[k]);

    a_age[k] ~ normal(a_age[k-1], sigma_evolution);

    a_edu[k] ~ normal(a_edu[k-1], sigma_evolution);

    a_race[k] ~ normal(a_race[k-1], sigma_evolution);

    a_party_id[k] ~ normal(a_party_id[k-1], sigma_evolution);

}

//Likelihood for each i in each wave

for(s in 1:n_waves){

      int obs_in_wave[wave_size[s]];                            

obs_in_wave<-segment(in_wave,wave_lims[s,1],wave_lims[s,2]) ;   // find observations in each wave

  for (i in 1:wave_size[s]) {

real y_hat;

    y_hat <- 

    beta[s]*gender[obs_in_wave[i]]

    +a_age[s,age[obs_in_wave[i]]]

    +a_edu[s,edu[obs_in_wave[i]]]

    +a_party_id[s,party_id[obs_in_wave[i]]]

    +a_race[s,race[obs_in_wave[i]]]

    +a_state[state[obs_in_wave[i]]];

    y[obs_in_wave[i]] ~ bernoulli_logit(y_hat);

  } 

  }

  }

[Gregory]

beta[k] ~ normal(beta[k-1], sigma_evolution[k]);

Not sure if this would be it, but there shouldn't be a [k] on that sigma_evolution I don't think.

[Tobias Konitzer]

Yes, sorry. This was a slip on my part, and the model I ran did not have the k subscript in the evolution parameter (yet does not converge, per R-Hat). Apologies! 

On Wednesday, March 23, 2016 at 7:26:21 AM UTC-7, Tobias Konitzer wrote:

[Bob Carpenter]

You might try putting priors on the k=1 case in your time
series, but that probably won't be enough. You might want
to also fit the different times separately and see if the fit
you gets indicate a time series. Or you might try adding one
effect at a time to see if you can get any of them to work.

- Bob

[Tobias Konitzer]

Hi All.

I have been playing around with the dynamic model for a long time. It turns out that using weakly informative priors as specified by Gelman (2008) to anchor the dynamic parameters helps; I also had to tweak the max_treedepth, but now the model converges beautifully, and has some of the desired properties (i.e. strong correlation between beta_t and beta_(t-1) that wears off slowly deeper down the temporal chain (model below). 

Two problems remain: 

First, my division random effects have a very wide posterior distribution for divisions in which I have no data (specifically, Division 3 and Division 6; I can post posterior results if desired). From my understanding, the posterior here should be pretty much the same as the prior, i.e. have 95% of the probability mass between -0.9 and 0.9, but the estimated standard deviation is >10. 

Second, I am trying to use the model results to get probability estimates for voting for specific demographic combinations, but I am unsure what metric to use to construct point estimates and credible intervals. My initial idea was to construct HPD intervals by doing something like  

 params_dem<-extract(fit)

best_dem<-which.max(params_dem$lp__)

and do the same for the 2.5% and 97.5% percentile. The problem here is that different model runs yield very different point estimates, which I guess is attributable to the built-in strong correlative structure of the model. I guess alternatively I can construct n.draws probability estimates for each demographic combination and then take the empirical percentiles, but even this yields substantively different point estimates. Having said that, the model converges perfectly in terms of the Rubin-Gelman diagnostic. I can't quite wrap my head around this, and any suggestions as to how to construct reliable point-estimates of probability and credible intervals based on the below model are more than welcome!

Many thanks again,

Tobi

The Model:

data {

  int<lower=0> N;                                         // sample wave_size

  int<lower=0> n_age;                                     // number of age categories

  int<lower=0> n_race;                                    // number of race categories

  int<lower=0> n_edu;                                     // number of education catoegories

  int<lower=0> n_division;                                // number of Divisions

  int<lower=0> n_states;                                  // number of States

  int<lower=0> n_party_id;                                // number of party-categories

  int<lower=0> n_waves;                                   // number of Waves

  int<lower=0,upper=n_division> division[N];              // Division for each respondent

  int<lower=0,upper=1> gender[N];                         // Gender each respondent

  int<lower=0,upper=n_age> age[N];                        // age for each respondent

  int<lower=0,upper=n_edu> edu[N];                        // edu for each respondent

  int<lower=0,upper=n_race> race[N];                      // race for each respondent

  int<lower=0,upper=n_party_id> party_id[N];              // party_id for each respondent

  real vote_prev_state[N];                                // 2012 state-vote at the individual-level

  int<lower=0,upper=1> y[N];                              // vote-choice

  int<lower=0,upper=N> in_wave[N];                        // wave of respondents

  int<lower=0>        wave_size[n_waves];                 // N individuals per Wave  

 int<lower=0> wave_lims[n_waves,2];                       // first individual in wave and individuals per wave                                                              (5*2 matrix)

}

parameters {

  

  real beta [n_waves];                                    // coefficient-array for female; wave               

  vector[n_age] a_age [n_waves];                          // coefficient-array for age; wave 

  vector[n_edu] a_edu [n_waves];                          // coefficient-array for edu; wave 

  vector[n_race] a_race [n_waves];                        // coefficient-array for race; wave 

  vector[n_party_id] a_party_id [n_waves];                // coefficient-array for party_id; wave 

  vector [n_division] a_division[n_waves];                // coefficient-array for division; wave 

  real gamma;                                             // previous vote- slope

}

model {

  

  gamma~cauchy(0,2.5);                                    // embed 95% probability mass between 0.01% and 99.9%                                                           // (Gelman 2008)

  a_division[1]~ cauchy(0,0.425);                         // embed 95% probability mass between 70% and 30%

  beta[1]~cauchy(0,2.5);                                  // embed 95% probability mass between 0.01% and 99.9%

  a_age[1]~cauchy(0,2.5);                                 // embed 95% probability mass between 0.01% and 99.9%

  a_edu[1]~cauchy(0,2.5);                                 // embed 95% probability mass between 0.01% and 99.9%

  a_race[1]~cauchy(0,2.5);                                // embed 95% probability mass between 0.01% and 99.9%

  a_party_id[1]~cauchy(0,2.5);                            // embed 95% probability mass between 0.01% and 99.9%

  for (k in 2:n_waves) {

    beta[k] ~ normal(beta[k-1], 0.05);                    //  assume 5% shift per demographic per wave covers                                                            // 95% of probability mass

    a_age[k] ~ normal(a_age[k-1], 0.05);                  // assume 5% shift per demographic per wave covers                                                             // 95% of probability mass

    a_edu[k] ~ normal(a_edu[k-1], 0.05);                  // assume 5% shift per demographic per wave covers 

                                                          // 95% of probability mass

    a_race[k] ~ normal(a_race[k-1], 0.05);                // assume 5% shift per demographic per wave covers                                                             // 95% of probability mass

    a_party_id[k] ~ normal(a_party_id[k-1], 0.05);        // assume 5% shift per demographic per wave covers                                                             // 95% of probability mass

    a_division[k] ~ normal(a_division[k-1], 0.05);

}

                                                          //Likelihood for each i in each wave

for(s in 1:n_waves){

      int obs_in_wave[wave_size[s]];                            

obs_in_wave<-segment(in_wave,wave_lims[s,1],wave_lims[s,2]) ;   // find observations in each wave

  for (i in 1:wave_size[s]) {

real y_hat;

    y_hat <- 

    beta[s]*gender[obs_in_wave[i]]

    +a_age[s,age[obs_in_wave[i]]]

    +a_edu[s,edu[obs_in_wave[i]]]

    +a_party_id[s,party_id[obs_in_wave[i]]]

    +a_race[s,race[obs_in_wave[i]]]+gamma*vote_prev_state[i]

    +a_division[s,division[obs_in_wave[i]]];

[Michael Betancourt]

Extreme quantiles, including min and max, are extremely noisy.  Stick

to at least 5/95, if not 10/90, to get more stable results.

[Andrew Gelman]

Yes, I agree.  And from more fundamental decision-analytic principles, I don’t think 2.5%/97.5% intervals make much sense.  We will need to change the default Stan output!  And also write a paper explaning why 2.5%/97.5% don’t make sense.

A

[Tobias Konitzer]

Thank you for the replies! I am looking forward to reading the paper... I assume the comments are uniquely meant for 

best_dem<-which.max(params_dem$lp__)

As far as I understand, assessing uncertainty by getting n.sim probability estimates in a logit model and then taking the  2.5%/97.5% is more valid? Or is this also too rigorous of an approach, and we would prefer 5/90 or something less restrictive even?

Many thanks again,

Tobi

[Bob Carpenter]

I'm not sure what Andrew's thinking, but part of the problem
is computational --- you need a lot of draws to estimate tail
quantiles (you can think of this as only 2.5% of the samples
coming into play to help estimate the lower bound of a 95%
interval).

Another issue is that tails can be more sensitive to priors
than central intervals or posterior means.

You can see all this through simulation by rerunning with different
seeds and different priors and looking at variance in the 2.5%
vs. 50% quantile estimates. More things to go in the book!

- Bob

[Andrew Gelman]

Hi, beyond the computational issue, my argument is that for decision making we want to know roughly where the parameter is but there’s no need to be requiring 97.5% certainty or anytihng close to that. Of course in prediction problems we might care about things like 100-year floods, but for basic parameter esitmation it seems to me that 50% intervals make more sense. Depends on context, of course, and I’ve struggled to figure out exactly how to write this up (even for blog, let alone in a paper or book).
A