hierarchical model with binary outcome

[Joshua Wiley]

Dear List,

I am working on a hierarchical model of child mortality (0/1).  Observations (N = 408,706) are nested within mothers (n = 141,999), cities (n = 3,806), districts (n = 443), and states (n = 26).  The outcome is quite unequal, with about 91.4% 0s and 8.6% 1s.

I would love any feedback on how tractable working with data of this size is in Stan.  Right now, I am trying a simple intercept only model, which I tried parameterizing two different ways.

To Stan's credit, both seem to be running.  The time it takes for a single iteration is a bit too long though --- a few minutes for 1 iteration.  Even if the effective sample size was exactly the number of iterations, it would take a couple days to get to 1,000, even in this simple model.  Suggestions would be welcome, it also may just be unrealistic to run multiple levels with this much data quickly.

Cheers,

Josh (code below):

First a simple way, about how I would think about the model/equations:

# simple way
data {
  int<lower=0> n4;
  int<lower=0> n3;
  int<lower=0> n2;
  int<lower=0> n1;
  int<lower=0> N;


  int<lower=1, upper=n4> ID4[N];
  int<lower=1, upper=n3> ID3[N];
  int<lower=1, upper=n2> ID2[N];
  int<lower=1, upper=n1> ID1[N];


  int<lower=0, upper=1> y[N];
}
parameters {
  real<lower=0> sigma_1;
  real<lower=0> sigma_2;
  real<lower=0> sigma_3;
  real<lower=0> sigma_4;


  vector[n4] int_4;
  vector[n3] int_3;
  vector[n2] int_2;
  vector[n1] int_1;
  real intercept;
}
transformed parameters {
  vector[N] p;


  for (i in 1:N) {
    p[i] <- inv_logit(intercept + int_4[ID4[i]] + int_3[ID3[i]] + int_2[ID2[i]] + int_1[ID1[i]]);
  }
}
model {
  sigma_4 ~ cauchy(0, 2);
  sigma_3 ~ cauchy(0, 2);
  sigma_2 ~ cauchy(0, 2);
  sigma_1 ~ cauchy(0, 2);


  int_4 ~ normal(0, sigma_4);
  int_3 ~ normal(0, sigma_3);
  int_2 ~ normal(0, sigma_2);
  int_1 ~ normal(0, sigma_1);
  intercept ~ normal(0, 3);


  y ~ bernoulli(p);
}

And here trying to reparameterize to scale and shift the parameters so the priors come from N(0, 1).

# Shift and Scale Approach
data {

  int<lower=0> n4;
  int<lower=0> n3;
  int<lower=0> n2;
  int<lower=0> n1;
  int<lower=0> N;


  int<lower=1, upper=n4> ID4[N];
  int<lower=1, upper=n3> ID3[N];
  int<lower=1, upper=n2> ID2[N];
  int<lower=1, upper=n1> ID1[N];


  int<lower=0, upper=1> y[N];
}
parameters {
  real<lower=0> sigma_1;
  real<lower=0> sigma_2;
  real<lower=0> sigma_3;
  real<lower=0> sigma_4;


  vector[n4] int_4raw;
  vector[n3] int_3raw;
  vector[n2] int_2raw;
  vector[n1] int_1raw;
  real intercept;
}
transformed parameters {
  vector[N] p;
  vector[n4] int_4;
  vector[n3] int_3;
  vector[n2] int_2;
  vector[n1] int_1;


  int_4 <- sigma_4 * int_4raw;
  int_3 <- sigma_3 * int_3raw;
  int_2 <- sigma_2 * int_2raw;
  int_1 <- sigma_1 * int_1raw;


  for (i in 1:N) {
    p[i] <- intercept + int_4[ID4[i]] + int_3[ID3[i]] + int_2[ID2[i]] + int_1[ID1[i]];
  }
}
model {
  sigma_4 ~ cauchy(0, 2);
  sigma_3 ~ cauchy(0, 2);
  sigma_2 ~ cauchy(0, 2);
  sigma_1 ~ cauchy(0, 2);


  int_4raw ~ normal(0, 1);
  int_3raw ~ normal(0, 1);
  int_2raw ~ normal(0, 1);
  int_1raw ~ normal(0, 1);
  intercept ~ normal(0, 3);


  y ~ bernoulli_logit(p);
}

[Christian Roy]

My first question would be: Do you include compilation in your running time? 

Some of my model can easily take up to >20 minutes to compile and run for a few iterations (25-50 samples). However, once the model is compiled the numbers of iterations by hours is much more reasonable. A model with a few thousands iteration is bound to take some time to compile... 

Chris

----------------------

 

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

[Joshua Wiley]

Hi Chris,

I did not include compile time (and it only takes a minute or two).

Cheers,

Josh

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

--
Joshua F. Wiley
Ph.D. Student, UCLA Department of Psychology
http://joshuawiley.com/
Senior Analyst, Elkhart Group Ltd.
http://elkhartgroup.com
Office: 260.673.5518

[Bob Carpenter]

I'd recommend removing the transformed parameter block, and replacing
the sampling statement in the model block with:

for (i in 1:N)

p[i] ~ bernoulli_logit(intercept + int_4[ID4[i]] + int_3[ID3[i]] + int_2[ID2[i]] + int_1[ID1[i]]);

Using vector arguments doesn't speed up bernoulli(_logit) as there are
no shared parameters.

Your model's only identified through priors, which can cause
substantial computational slowdown if the priors aren't particularly
tight. It's OK to center them all around zero and push the means into
the shared intercept.

If you've fit the data using some other GLM software, you might want
to see what the values are for coefficients to make sure the prior for
the intercept is reasonable.

Even if this doesn't give you a big speedup, you may not need 1000 iterations.
You only need to run until your MCMC std error is below where you want it.
And you can run in parallel if you have enough memory, though you have a very
large number of parameters (150K just for the mothers!), which could be prohibitive.

Run time will be much slower during warmup, then go faster per iteration
after that. You might also not need much warmup for this model, but I'm
guessing the mother parameters are going to be relatively hard to identify
given the low counts per mother.

And if worse comes to worst, 1000 minutes is still less than 24 hours :-)

- Bob

[Christian Roy]

Have you tried with a really basic model? Like only the intercept first, then with only the mothers, then mothers nested in cities... etc.?

Just to get an idea of how the simplest model perform? 

It could also give you a goo idea of much time is needed to deal with the size of the data.

Chris

-----------------

[Bob Carpenter]

Compile time only depends on the size of the model specification,
and the complexity of operations involved (inverse wishart is more
compile intense than addition), but does not depend on the size of
the data.

Reading data into the model will depend on the size of the data, but
that's only done once.

In addition, the 9% vs. 91% balance isn't enough to worry about. Should
be no problem for logistic regression.

- Bob

[Bob Carpenter]

That's a good idea. You can also see how well it fits
with just some of the coefficients. You'll need the overall
intercept the way your model works in all cases. I have to imagine
the mothers are dominating the compute time.

- Bob

[Joshua Wiley]

Thanks!  These are very helpful suggestions, and good to know about bernoulli(_logit) and vectorization.

I have run this model in MCMCglmm in R before, so have some other numbers for reference.

I am trying these suggestions out and will post back how it goes.

Thanks! :)

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

[Luc Coffeng]

Hi Joshua,

From your code, I understand (and correct me if I'm wrong) that you are not using individual level predictors, right? If you are sticking to this intercept-only model, you could crank down significantly on computation by collapsing you data from (0/1) bernoulli responses to (k/N_kids) binomial responses, where N_kids is the number of children a mother had, and k is the number of kids that died. Then just model the binomial likelihood instead of the bernoulli:

for (i in 1:N_mothers)
    k[i] ~ binomial_logit(N_kids[i], intercept + int_4[ID4[i]] + int_3[ID3[i]] + int_2[ID2[i]] + int_1[ID1[i]]);

If on average mothers in your data have 2 kids, this means that you can compress your data (and number of likelihood evaluations) by a factor two. If at some point, you want to add individual-level predictors and these predictors are categorical (e.g. gender), you might to still apply aggregation, but now at the level of kid-gender | mother. Of course, if the number of kids per mother is low, this won't really have a significant impact.

Good luck,

Luc

[Joshua Wiley]

Luc,

That is a great idea --- I will eventually add predictors, but they would be at the level of mother or above.  The average kids per mother was a bit over 2, so that would actually reduce the data a lot.

Thank you! :)

Josh

--

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

[Bob Carpenter]

The binomial should be quite a bit faster --- it saves doing multiple
additions and inverse-logit calculations, which is where the
time's going to be spent in this model. It does add some multiplications,
though. So it'll be interesting to me to see what the difference is if
the average N_kids is 2.

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