Speeding up Dynamic Regression Model with many time-periods and big sample size

[Tobias Konitzer]

HI All,

I have written a dynamic model predicting Democratic vote (attached). The problem: With a sample size of about 600,000 and 46 time-periods the model takes forever to run (about a week per 50 iterations). Is there any obvious ways to speed up the model?

Thanks,

Tobi

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_dma;                                     // number of dmas

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

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

  int<lower=0,upper=1> female[N];                         // female 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

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

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

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

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

  int<lower=0>        day_size[n_day];                 // N individuals per Wave  

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

}

parameters {

  

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

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

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

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

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

  vector [n_dma] a_dma[n_day];                // coefficient-array for division; wave 

  real gamma;                                             // previous vote- slope

  real alpha;                                             //intercept

   real <lower=0> sigma_dma;

  real <lower=0> sigma_age;

  real <lower=0> sigma_edu;

  real <lower=0> sigma_race;

  real <lower=0> sigma_party_id;

  real <lower=0> sigma_beta_evo;

  real <lower=0> sigma_dma_evo;

  real <lower=0> sigma_age_evo;

  real <lower=0> sigma_edu_evo;

  real <lower=0> sigma_race_evo;

  real <lower=0> sigma_party_id_evo;

}

model {

  sigma_dma ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_age ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_edu ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_race ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_party_id ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_beta_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

     sigma_dma_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_age_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_edu_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_race_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    sigma_party_id_evo ~normal(0,2.5);                          // within_division variance can move from 0.1% and 99%

    

  gamma~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%   

    alpha~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%                                                           // (Gelman 2008)

// (Gelman 2008)

  beta[1]~normal(0,5);                                  // embed 95% probability mass between 0.1% and 99%

  a_dma[1]~ normal(0,sigma_dma);                         // embed 95% probability mass between 70% and 30%

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

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

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

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

  

  for (k in 2:n_day) {

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

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

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

    // 95% of probability mass

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

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

    a_dma[k] ~ cauchy(a_dma[k-1], sigma_dma_evo);

    

  }

  //Likelihood for each i in each wave

  for(s in 1:n_day){

    int obs_in_day[day_size[s]];                            

    

    obs_in_day<-segment(in_day,day_lims[s,1],day_lims[s,2]) ;   // find observations in each wave

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

      real y_hat;

      y_hat <- alpha

        +beta[s]*female[obs_in_day[i]]

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

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

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

      +a_race[s,race[obs_in_day[i]]]+gamma*vote_prev_dma[i]

      +a_dma[s,dma[obs_in_day[i]]];

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

    } 

  }

}

[Bob Carpenter]

Moving to non-centered parameterizations will probably help for
group-level parameters where you don't have many elements.

So instead of

gamma ~ foo(loc, scale);

use

gamma_raw ~ foo(0, 1);
gamma <- loc + scale * gamma_raw;

where gamma_raw is the parameter and gamma is a local variable.

Also, these:

> for (k in 2:n_day) {
> beta[k] ~ cauchy(beta[k-1], sigma_beta_evo);

are much more efficient as:

beta[2:n_day] ~ cauchy(beta[1:(n_day - 1)], sigma_beta_evo);

This:

> for (i in 1:day_size[s]) {
> real y_hat;
> y_hat <- alpha
> +beta[s]*female[obs_in_day[i]]
> +a_age[s,age[obs_in_day[i]]]
> +a_edu[s,edu[obs_in_day[i]]]
> +a_party_id[s,party_id[obs_in_day[i]]]
> +a_race[s,race[obs_in_day[i]]]+gamma*vote_prev_dma[i]
> +a_dma[s,dma[obs_in_day[i]]];
> y[obs_in_day[i]] ~ bernoulli_logit(y_hat);
> }

can be vectorized as:

y[obs_in_day] ~ bernoulli_logit(alpha
+ beta[s] * female[obs_in_day]
+ ...);

If you make female a vector.

If female[...] is an indicator, then it's acutally more efficient to
build up that linear predictor using conditionals to avoid the unnecessary
multiplication by zero and useless derivative propagations.

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

Hi!

Thanks, these make a lot of sense, especially with regard to vectorization. however, I get an error 

No matches for:

  vector[] ~ cauchy(vector[], real)

Although one of the listed argument signatures for cauchy is indeed vector ~ cauchy(vector, real)

?

Thanks,

Tobi

[Bob Carpenter]

vector[] is an array of vectors. So you need either one
more level of looping or you need to apply to_vector if
they have the same dims:

vector[5] y[3];
vector[5] mu[3];

y ~ cauchy(mu, 1.0)

gives the error you see.

But if you do this:

matrix[3,5] y;
matrix[3,5] mu;

to_vector(y) ~ cauchy(to_vector(mu), 1.0);

it'll work.

- Bob

[Tobias Konitzer]

Hi!

Thanks again. It took me a while to figure out how to arrange the random-effect matrices etc., but here is the model I came up with. I am not sure about my implementation of the non-centered parametrization and would love a quick comment! Apart from that, this settles it for me. I will keep you in the loop wrt run-time (650,000 N sample, 48 waves)!

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_dma;                                     // number of dmas

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

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

  vector <lower=0,upper=1> [N]female;                       // female 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

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

  vector[N] vote_prev_dma;                                // 2012 state-vote at the individual-level

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

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

  int<lower=0>        day_size[n_day];                    // N individuals per Wave  

  int<lower=0> day_lims[n_day,2];                         // first individual in wave and individuals per wave                                   

}

parameters {

  

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

  matrix[n_day,n_age] a_age_raw;                          // coefficient-matrix for age; wave 

  matrix[n_day,n_edu] a_edu_raw;                          // coefficient-matrix for edu; wave 

  matrix[n_day,n_race] a_race_raw;                        // coefficient-matrix for race; wave 

  matrix[n_day,n_party_id] a_party_id_raw ;               // coefficient-matrix for party_id; wave 

  matrix[n_day,n_dma] a_dma;                              // coefficient-matrix for DMA; wave 

  real gamma;                                             // previous vote- slope

  real alpha;                                             //intercept

   real <lower=0> sigma_dma;

  real <lower=0> sigma_age;

  real <lower=0> sigma_edu;

  real <lower=0> sigma_race;

  real <lower=0> sigma_party_id;

  real <lower=0> sigma_beta_evo;

  real <lower=0> sigma_dma_evo;

  real <lower=0> sigma_age_evo;

  real <lower=0> sigma_edu_evo;

  real <lower=0> sigma_race_evo;

  real <lower=0> sigma_party_id_evo;

}

transformed parameters{

  matrix[n_day,n_age] a_age;                            //efficient non-centered paramterization

  matrix[n_day,n_edu] a_edu;

  matrix[n_day,n_race] a_race;

  matrix[n_day,n_party_id] a_party_id;

  a_age[1]<-a_age_raw[1]*sigma_age;

  a_edu[1]<-a_edu_raw[1]*sigma_edu;

  a_race[1]<-a_race_raw[1]*sigma_race;

  a_party_id[1]<-a_party_id_raw[1]*sigma_party_id;

for (i in 2:n_day){

  a_age[i]<-a_age_raw[i];

  a_edu[i]<-a_edu_raw[i];

  a_race[i]<-a_race_raw[i];

  a_party_id[i]<-a_party_id_raw[i];

}

}

model {

  sigma_dma ~normal(0,2.5);                             // within_DMA variance can move from 0.1% and 99%

    sigma_age ~normal(0,2.5);                           // within_Age variance can move from 0.1% and 99%

    sigma_edu ~normal(0,2.5);                           // within_edu variance can move from 0.1% and 99%

    sigma_race ~normal(0,2.5);                          // within_race variance can move from 0.1% and 99%

    sigma_party_id ~normal(0,2.5);                      // within_party variance can move from 0.1% and 99%

    sigma_beta_evo ~normal(0,2.5);                      // Evolution Variance for Female

    sigma_dma_evo ~normal(0,2.5);                       // Evolution Variance for DMA

    sigma_age_evo ~normal(0,2.5);                       // Evolution Variance for Age

    sigma_edu_evo ~normal(0,2.5);                       // Evolution Variance for Edu

    sigma_race_evo ~normal(0,2.5);                      // Evolution Variance for race

    sigma_party_id_evo ~normal(0,2.5);                  // Evolution Variance for Party

    

  gamma~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%   

  alpha~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%                                                           /

  beta[1]~normal(0,5);                                  // embed 95% probability mass between 0.1% and 99%

  a_dma[1]~ normal(0,1);                                // embed 95% probability mass between 70% and 30%

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

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

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

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

                                                        // Evolution Sampling:

 

 beta[2:n_day] ~ cauchy( beta[1:(n_day - 1)], sigma_beta_evo); 

 to_vector(a_age[2:n_day]) ~  cauchy( to_vector(a_age[1:(n_day - 1)]), sigma_age_evo);     

 to_vector(a_edu[2:n_day]) ~  cauchy(to_vector( a_edu[1:(n_day - 1)]), sigma_edu_evo);     

 to_vector(a_race[2:n_day]) ~ cauchy(to_vector( a_race[1:(n_day - 1)]), sigma_race_evo);     

 to_vector(a_party_id[2:n_day]) ~ cauchy( to_vector(a_party_id[1:(n_day - 1)]), sigma_party_id_evo);     

to_vector(a_dma[2:n_day]) ~ cauchy( to_vector(a_dma[1:(n_day - 1)]), sigma_dma_evo);     

                                                    

  for(s in 1:n_day){                                    // Likelihood for each  each wave vectorized

   int obs_in_day[day_size[s]];                            

                                                        // find observations in each wave:

    obs_in_day<-segment(in_day,day_lims[s,1],day_lims[s,2]);

                                                        // Likelihood:

    y[obs_in_day] ~ bernoulli_logit(alpha

        +beta[s]*female[obs_in_day]

      +to_vector(a_age[s,age[obs_in_day]])

      +to_vector(a_edu[s,edu[obs_in_day]])

      +to_vector(a_party_id[s,party_id[obs_in_day]])

      +to_vector(a_race[s,race[obs_in_day]])

      +gamma*vote_prev_dma[obs_in_day]

      +to_vector(a_dma[s,dma[obs_in_day]]));

  }}

[Bob Carpenter]

These pieces are still centered:

to_vector(a_age[2:n_day]) ~ cauchy( to_vector(a_age[1:(n_day - 1)]), sigma_age_evo);
to_vector(a_edu[2:n_day]) ~ cauchy(to_vector( a_edu[1:(n_day - 1)]), sigma_edu_evo);
to_vector(a_race[2:n_day]) ~ cauchy(to_vector( a_race[1:(n_day - 1)]), sigma_race_evo);
to_vector(a_party_id[2:n_day]) ~ cauchy( to_vector(a_party_id[1:(n_day - 1)]), sigma_party_id_evo);
to_vector(a_dma[2:n_day]) ~ cauchy( to_vector(a_dma[1:(n_day - 1)]), sigma_dma_evo);

You should also check for these that the noise really does have crazy long
tails. You'll find that if you do something like posterior predictive
checks, you'll go way way way out into the tails with a Cauchy.

Student t with d.o.f. = 4 is a compromise between Cauchy and normal.

- Bob

> On Apr 26, 2016, at 9:08 PM, Tobias Konitzer <tobias....@gmail.com> wrote:
>
> Hi!
>

> Thanks again. It took me a while to figure out how to arrange the random-effect matrices etc., but here is the model I came up with. I am not sure about my implementation of the non-centered parametrization and would love a quick comment! Apart from that, this settles it for me. Many thanks again!

[Tobias Konitzer]

Hi All, 

I have implemented the changes, and have settled on a t-distribution in the dynamic part as a good compromise indeed, thanks for that tip! However, the model still takes forever (literally <10 iterations per day). I have tested the same model with smaller N and less time-periods, and it converges well, albeit is running slowly. Could it be that N>600,000 and 46 time-periods push stan to the limits? Model attached, many thanks again for all the excellent suggestions!

Tobi

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_dma;                                     // number of dmas

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

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

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

  vector <lower=0,upper=1> [N]female;                     // female 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

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

  vector[N] vote_prev_dma;                                // 2012 state-vote at the individual-level

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

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

  int<lower=0>        day_size[n_day];                    // N individuals per Wave  

  int<lower=0> day_lims[n_day,2];                         // first individual in wave and individuals per wave                                   

}

parameters {

  

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

  matrix[n_day,n_age] a_age_raw;                          // coefficient-matrix for age; wave 

  matrix[n_day,n_edu] a_edu_raw;                          // coefficient-matrix for edu; wave 

  matrix[n_day,n_race] a_race_raw;                        // coefficient-matrix for race; wave 

  matrix[n_day,n_party_id] a_party_id_raw ;               // coefficient-matrix for party_id; wave 

  matrix[n_day,n_dma] a_dma_raw;                          // coefficient-matrix for DMA; wave 

  real gamma;                                             // previous vote- slope

  real alpha;                                             //intercept

   real <lower=0> sigma_dma;

  real <lower=0> sigma_age;

  real <lower=0> sigma_edu;

  real <lower=0> sigma_race;

  real <lower=0> sigma_party_id;

  real <lower=0> sigma_beta_evo;

  real <lower=0> sigma_dma_evo;

  real <lower=0> sigma_age_evo;

  real <lower=0> sigma_edu_evo;

  real <lower=0> sigma_race_evo;

  real <lower=0> sigma_party_id_evo;

}

transformed parameters{

  matrix[n_day,n_age] a_age;                            //efficient non-centered paramterization

  matrix[n_day,n_edu] a_edu;

  matrix[n_day,n_race] a_race;

  matrix[n_day,n_party_id] a_party_id;

  matrix[n_day,n_dma] a_dma;

  a_age[1]<-a_age_raw[1]*sigma_age;

  a_edu[1]<-a_edu_raw[1]*sigma_edu;

  a_race[1]<-a_race_raw[1]*sigma_race;

  a_party_id[1]<-a_party_id_raw[1]*sigma_party_id;

  a_dma[1]<-a_dma_raw[1]*sigma_dma;

for (i in 2:n_day){

  a_age[i]<-a_age_raw[i]*sigma_age_evo;

  a_edu[i]<-a_edu_raw[i]*sigma_edu_evo;

  a_race[i]<-a_race_raw[i]*sigma_race_evo;

  a_party_id[i]<-a_party_id_raw[i]*sigma_party_id_evo;

  a_dma[i]<-a_dma_raw[i]*sigma_dma_evo;

}

}

model {

  sigma_dma ~normal(0,2.5);                             // within_DMA variance can move from 0.1% and 99%

    sigma_age ~normal(0,2.5);                           // within_Age variance can move from 0.1% and 99%

    sigma_edu ~normal(0,2.5);                           // within_edu variance can move from 0.1% and 99%

    sigma_race ~normal(0,2.5);                          // within_race variance can move from 0.1% and 99%

    sigma_party_id ~normal(0,2.5);                      // within_party variance can move from 0.1% and 99%

    sigma_beta_evo ~normal(0,2.5);                      // Evolution Variance for Female

    sigma_dma_evo ~normal(0,2.5);                       // Evolution Variance for DMA

    sigma_age_evo ~normal(0,2.5);                       // Evolution Variance for Age

    sigma_edu_evo ~normal(0,2.5);                       // Evolution Variance for Edu

    sigma_race_evo ~normal(0,2.5);                      // Evolution Variance for race

    sigma_party_id_evo ~normal(0,2.5);                  // Evolution Variance for Party

    

  gamma~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%   

  alpha~normal(0,5);                                    // embed 95% probability mass between 0.01% and 99.9%                                                           /

  beta[1]~normal(0,5);                                  // embed 95% probability mass between 0.1% and 99%

  a_dma[1]~ normal(0,1);                                // embed 95% probability mass between 70% and 30%

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

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

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

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

                                                        // Evolution Sampling:

 

 beta[2:n_day] ~ student_t(nu, beta[1:(n_day - 1)], 1); 

 to_vector(a_age[2:n_day]) ~  student_t(nu, to_vector(a_age[1:(n_day - 1)]), 1);     

 to_vector(a_edu[2:n_day]) ~  student_t(nu,to_vector( a_edu[1:(n_day - 1)]), 1);     

 to_vector(a_race[2:n_day]) ~ student_t(nu,to_vector( a_race[1:(n_day - 1)]), 1);     

 to_vector(a_party_id[2:n_day]) ~ student_t(nu, to_vector(a_party_id[1:(n_day - 1)]), 1);     

to_vector(a_dma[2:n_day]) ~ student_t(nu, to_vector(a_dma[1:(n_day - 1)]), 1);     

[Andrew Gelman]

10 iterations per day . . . that does sound slow!  This reminds me of various discussions we’ve had about data-splitting algorithms such as here:

http://www.stat.columbia.edu/~gelman/research/unpublished/ep.pdf

A

[Bob Carpenter]

That would be pushing Stan to the limit even if the model were
more cleanly identified. As is, this linear predictor

> y[obs_in_day] ~ bernoulli_logit(alpha
> +beta[s]*female[obs_in_day]
> +to_vector(a_age[s,age[obs_in_day]])
> +to_vector(a_edu[s,edu[obs_in_day]])
> +to_vector(a_party_id[s,party_id[obs_in_day]])
> +to_vector(a_race[s,race[obs_in_day]])
> +gamma*vote_prev_dma[obs_in_day]
> +to_vector(a_dma[s,dma[obs_in_day]]));
> }}

is only weakly identified because of all the varying
intercepts.

Have you thought about only including intercepts for some
of the alternatives? That is, if there are 5 age cohorts,
only include predictors for four of them and roll the last
into the intercept --- it can help with identification. But
even with this kind of recoding, this scale of data's probably
more than Stan can handle.

Otherwise, I don't see any alternative other than subsampling
the data to fit, fitting pieces of the model separately, or
moving to an approximate algorithm like ADVI (which we don't
understand well in terms of how fits relate to the full Bayesian
posterior).

- Bob

[Tobias Konitzer]

Thanks for the comments! Could you explain what you mean by "fitting pieces of the model separately"? The way the model is written (implicitly Markov Chain), I should be able to run the model separately for each time period. But I have no idea how I would go about it...

Many thanks again,

Tobi

[Bob Carpenter]

I meant fitting different models to different time periods.
It's just another way to subsample the data. Or what
Andrew calls "the secret weapon".

http://andrewgelman.com/2009/05/24/handy_statistic/

- Bob

[Julian King]

1) For the dynamic part of the model, why not try a normal distribution instead of a Cauchy/t distribution?  If this gives you much faster results, it tells you something about the evolution process.

2) Is there any way you can aggregate the individuals into groups (eg # of people in a certain age bin, # of people in a certain party ID etc) and then use a binomial likelihood instead of a bernoulli likelihood?  This way you will have a much smaller number of data points to loop over.   This effectively means you are doing binomial regression.  Assuming you have 10 age bins, 5 edu bins, 3 race bins, 2 party ID bins, 2 DMA bins (I don't know what DMA is), and 2 gender bins, this is a maximum of 1200 unique bins.  This means that your likelihood has 500x less terms in the calculation, so you should expect a substantial speedup.  Explanation for this mathematically is here: https://en.wikipedia.org/wiki/Logistic_regression#In_terms_of_binomial_data

3) To go further down this path, if you do binomial regression and your sample size in each bin is large enough, you can approximate the binomial likelihood with a Gaussian, and then apply the Kalman filter to solve your problem in one go (as a Dynamic Linear Model).  Note that Stan now appears to have a function for this directly - gaussian_dlm_obs, so the effort in coding the DLM should be fairly low. 

4) As noted by others you are going to run into identifiability problems here because of multicollinearity in the predictors.  I ran into this problem recently - see here.  You need to drop one out of each of the predictors so that your predictors are linearly independent.  Although the Bayesian priors in theory solve the identifiability issue, you'll find that your posteriors have extremely high correlations on these parameters (probably >0.99), so that the standard NUTS sampling will be very slow.  A dense mass matrix will solve this, but you'll never be able to get enough posterior samples to estimate the mass matrix reliability (ie chicken and egg problem), so in practice dropping one out of each of the predictors and folding it into the intercept as suggested is a much better option. 



The advantage of using the Kalman filter for the evolving parameters is that it integrates out the random walk of the parameters directly, giving an analytical formula for the posterior. 

[Julian King]

Also, note that by doing this problem as suggested via binomial regression, rather than logistic regression, the number of terms in the likelihood becomes independent of sample size.

[Tobias Konitzer]

Hi All,

Sorry it took me a while to respond/implement changes.

Unfortunately, I have over 100,000 bins, so I think Julian's suggestion does not give me much!

I have tried removing one condition for all random effect variables and have this be swallowed by the non-random intercept. Another idea that came to mind is to run the model piecewise, i.e. separately for the first ten time periods and then for the second ten time periods, and then "set" the distributions in the second run for t=1 to Normal with the observed mean and standard deviation of the previous run. Is this a legitimate approach?

On the other hand, the following model for N=80,000 and t=8 still runs on about 10 iterations a day, which seems crazy slow:

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

  real <lower=0> sigma_division;

  real <lower=0> sigma_age;

  real <lower=0> sigma_edu;

  real <lower=0> sigma_race;

  real <lower=0> sigma_party_id;

  real <lower=0> sigma_evolution_beta;

  real <lower=0> sigma_evolution_age;

  real <lower=0> sigma_evolution_edu;

  real <lower=0> sigma_evolution_race;

  real <lower=0> sigma_evolution_party_id;

  real <lower=0> sigma_evolution_division;

}

model {

  sigma_division ~normal(0,1.25);                          // within_division variance can move from 0.1% and 99%

  sigma_age ~normal(0,1.25);                          // within_division variance can move from 0.1% and 99%

  sigma_edu ~normal(0,1.25);                          // within_division variance can move from 0.1% and 99%

  sigma_race ~normal(0,1.25);                          // within_division variance can move from 0.1% and 99%

  sigma_party_id ~normal(0,1.25);                          // within_division variance can move from 0.1% and 99%

  sigma_evolution_beta~normal(0,0.1);

  sigma_evolution_age~normal(0,0.1);

  sigma_evolution_edu~normal(0,0.1);

  sigma_evolution_race~normal(0,0.1);

  sigma_evolution_party_id~normal(0,0.1);

  sigma_evolution_division~normal(0,0.1);

  beta[1]~normal(0,2.5);                                  // embed 95% probability mass between 0.1% and 99%

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

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

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

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

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

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

  for (k in 2:n_waves) {

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

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

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

                                                          // 95% of probability mass

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

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

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

}

                                                          //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]]];

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

  } 

  }

  }

Thanks,

Tobi

On Monday, April 25, 2016 at 10:44:45 AM UTC-7, Tobias Konitzer wrote:

[Andrew Gelman]

Hi, I have 3 thoughts. First, yes, subset the data for sure.  We have ideas of how to do this more formally (see here: http://www.stat.columbia.edu/~gelman/research/unpublished/ep.pdf), but even a simple partitioning will often make sense.  Second, strip down your model.  Remove most of these varying intercepts. Start with a model with only one of these sigmas.  That should run fast.  Then add them, one at a time, and see how that goes.

Third, once you start to add in these factors, start by giving them fixed values of sigma.  For example, sigma=0.1 or 0.2.  Remember, excluding a factor is just sigma=0.  In lots of these examples, these sigmas aren't so well identified anyway.  Run the models with fixed sigmas, that should be fast.  In fact, with fixed sigmas this is just a logistic regression model with a prior, you could even do point estimation, perhaps.  Then replaces some of the fixed sigmas with priors, etc.  Build up to it.  Ideally the entire process should just take a few days.

10 iterations a day is ridiculous:  there's no way you should trust your model enough to sit around and wait like that.  I'd have more confidence in a model with externally-supplied sigmas that could run fast so you could try out different possibilities.

Also this example could be a good candidate for GMO, perhaps.

A

[Bob Carpenter]

> On May 16, 2016, at 9:07 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>
> ...

> Also this example could be a good candidate for GMO, perhaps.

Or ADVI. You'd have to make sure the approximate answers
you get are close to the real ones.

- Bob