Fitting models of mixture distributions with continuous portions + discrete spikes

[Rachael Meager]

Hi Folks,

I'm trying to explicitly model some distributions that have very fat tails and a discrete spike at zero. The two tails can sometimes look very different to each other and I want to capture that. To give you some context: this is profit data at a household level from some developing countries. Some households in my data make positive business profits, some make negative, and some make zero or don't operate a business at all (which gets recorded as an exact zero). I need to keep all of these households in the sample because I want to study the distributional effects of a binary treatment variable on these people, which might include moving some people out of the spike into the tails or vice versa. I'm attaching some histograms as examples. The plots are zoomed in / truncated on both axes - the true spike is much higher and the true tails go much further out, but you get the idea.

After some discussions with econometrics and EECS folks around MIT, and using some prior information from economic theory that suggests income/profit/wealth variables typically follow power laws, I figured to start with a mixture model with a discrete spike at zero and 2 pareto distributions as the 2 tails (obviously the negative tail has to be a kind of "reverse" pareto). I want to do inference about the mixture assignment probabilities and the parameters that control the shapes of the tails. I know Stan can't do spike+slab models, so I wrote a Gibbs sampler. However, I don't like having to use conjugate priors on the pareto distribution parameters, and I eventually want to put all of these distributions into a big hierarchical model which might be a stretch of my Gibbs sampling abilities. The friendly person behind @mcmc_stan on twitter suggested there should be a workaround in stan, so I'm posting here to see what people think about this problem and/or if anyone has any suggestions.  I looked at the previous discussion topics for spike and slabs but they mostly seemed to be about discouraging people from using spike+slab priors -- and yes, I would absolutely never use spikes in my priors, but in this data there is a genuine discrete spike and I want to model it.

Thanks so much in advance for any input or advice you have! 

Cheers

Rachael 

[James Savage]

Hey Rachael - 

I have similar issues working with microfinance repayments data--lots of 0s, lots of due amounts, and different distributions on either side. 

Ben Goodrich suggested a work-around that I'm finding works very nicely: have a model for the type of person, and then a model for profits within each person type. So you could have a categorical logit model for whether the person falls into a loss-making/zero profit/profitable category (use k-1 parameterisation with informative priors), and then model the loss or profit with and exponential or lognormal distribution or whatnot, conditional on the category. This should be a lot easier than mixture modeling, and won't have all the weird label-switching problems. 

The disadvantage is that your parameters aren't as interpretable, so any treatment effects estimate will be the "derivative" of the posterior prediction with respect to the treatment. But that's not hard to simulate. 

Jim

[Michael Betancourt]

Spike-and-slab is a bad _prior_ for various reasons.  But mixtures of discrete and continuous

parameters is fine for the likelihood.  See the zero-inflation models in the manual for an example.

These mixtures also tend to be well-identified and fit well; I use zero-inflation all the time.

It’s just a small jump to the three-component model you want to build.

--
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.
<profit_histogram_plots_4.pdf><profit_histogram_plots_6.pdf><profit_histogram_plots_7.pdf>

[Andrew Gelman]

Regarding interpretable parameters:

I think it _is_ a good idea for parameters to be interpretable in the sense of having an understandable meaning within a generative model.

But I do _not_ think it necessary or often appropriate for parameters in a model to align with quantities of interest or relevant causal summaries.  For a simple example, consider a 2-part model which first has a logit probability of income exceeding zero, and then a log regression for income, conditional on it being zero (similar issues arise in survival-type models).  Then a total causal effect will be expressed as the difference in expected income, which tangles both parts of the model and is not a single coefficient.  And that’s fine:  people should abandon the idea that they have to set up their model so that a particular “beta” is the parameter they care about.

A

[Rachael Meager]

Hi All! Thanks so much for the replies! 

@Jim: You're always 1 step ahead of me! :) This multinomial logit sounds like a good plan. I was going to ask about how you model the spike without the dirac delta function, but then I realised it should be sufficient to simply sort people into the bin for "zero profit" using the logit, and then draw/assign them y <- 0 manually in the code. I think that *is* definitionally a dirac delta function, in fact. 

@Michael: I'll check out zero inflation models, thanks! 

@Andy: Thanks for the reality check on interpretable parameters! Sounds reasonable -- and after thinking about it, I would add that for more complex models it won't generally be possible to summarize the effect of interest in 1 parameter anyway!  Or perhaps not without some serious contortions. 

[Rachael Meager]

Update on this thread + a question!

So I got both bits of the model working well - the multinomial logit and the pareto fitting - but I am having issues when I try to (1) turn them hierarchical and (2) smush them together into a single model.

The main issue I am having today is that my hierarchical multinomial logit is very slow. I tried to write the model as close as I could to a hierarchical vectorised/ "efficient" version of the model on page 57 of the stan manual (Section 6.5). I'm attaching the R and stan files here. Likely due to my inept coding, stan takes about 8 hours to do 5000 draws which don't converge. Of course the non-vectorised version is worse... but when the code runs for 8 hours, my having "efficient" appended to the file name seems rather sad! 

Does anyone here spot the problem with my stan code that produces this super slow sampling? Jim, you said earlier that you would use k-1 parameterisation (which I assume isn't what the manual did nor what I did) -- but would that really affect the speed much here? 

I'm copying the stan code below just for ease of reference. The notation and parameters are the same as in the manual example on page 57, the only new thing is that I have K different "sites" forming the lowest level of my hierarchy.

data {

int M; 

int N;

int P;

int K;

int cat[N];

vector[P] x[N];

int site[N];

}

parameters {

matrix[M,P] beta;

real<lower=0> sigma[M,P];

matrix[M,P] beta_k[K];

}

model {

to_vector(beta) ~ normal(0,5);

for (m in 1:M){

  sigma[m] ~ cauchy(0,3);

  for (k in 1:K){

    beta_k[k,m] ~ normal(beta[m],sigma[m]);

  }}

for (n in 1:N)

cat[n] ~ categorical_logit(beta_k[site[n]] * x[n]);

}

[Andrew Gelman]

Is it ordered logit or unordered?  If it’s ordered, I suggest just treating the data as continuous and running a normal regression.

A

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

<multinomial-logit-hierarchical-efficient.stan><stan_for_multinomial_logit_hierarchical.R>

[Rachael Meager]

It is ordered in some sense, since this is a categorisation based on being negative, zero, or positive. But here we think the negative + positive guys have more in common with each other - and probably switch into each others' buckets more often - than with the zero guys, so I don't think we can really use the structure of ordered logits. Although I may be mistaken on this point, and would be happy to hear a counterargument!

(Also, am I right in thinking that a continuous normal regression would be the "linear probability model" (LPM) which could spit out answers that don't satisfy the axioms of probability? Just want to check my understanding.) 

I'd be more inclined to switch to hierarchical dirichlet models before an LPM unless anyone sees a big problem with that...

[Bob Carpenter]

> On Apr 14, 2016, at 4:37 PM, Rachael Meager <rachael...@gmail.com> wrote:
>
> It is ordered in some sense, since this is a categorisation based on being negative, zero, or positive. But here we think the negative + positive guys have more in common with each other - and probably switch into each others' buckets more often - than with the zero guys, so I don't think we can really use the structure of ordered logits. Although I may be mistaken on this point, and would be happy to hear a counterargument!


> (Also, am I right in thinking that a continuous normal regression would be the "linear probability model" (LPM) which could spit out answers that don't satisfy the axioms of probability? Just want to check my understanding.)

I don't know what "linear probabilty model" means, but we'd use an ordered logit,
which is well-behaved in terms of probability. See the manual.

> I'd be more inclined to switch to hierarchical dirichlet models before an LPM unless anyone sees a big problem with that...

The problem with Dirichlet is that it's not a very flexible
prior and it's hard to add predictors (compared to a logistic
regression with multi-lognormal prior, say).

- Bob

[Bob Carpenter]

The main problem is that you are only going to get weak identification
using a K-dimensional parameterization of a K-way multinomial output. Especially
if you include a hierarchical model with varying scale --- because the scale of
the predictors is not defined given that it's a K-dimensional parameterization.
You can move the K-1 dimensional parameterization, which makes the priors
asymmetric, but makes everything well behaved. It's the only thing you'd see
in a max likelihood setting.

You can also vectorize this

> for (m in 1:M){
> sigma[m] ~ cauchy(0,3);

by declaring sigma as a matrix and using to_vector(sigma) on the left.

- Bob

[Rachael Meager]

Thank you for that explanation, Bob! I will switch to the K-1 parameterisation -- I didn't realise that could affect the speed, but when you explained the situation, it makes a lot of sense. (Clearly I am not so familiar with these models...). I'll report back here with new code and let you know how it goes.

EDIT: whoops just realised we switched the indices here. It's an M-way multinomial logit, a K-way hierarchy. So we want to switch to M-1 not K-1. 

[Rachael Meager]

Bob, you are my hero! I changed to the suggested "#Categories-1" parameterisation (which, I note, would have been where I started had I just listened to Jim Savage... I have learned my lesson now.) So 5000 iterations of the new model takes 30 minutes, and it behaves extremely well (most Rhats are 1.02, though got a few Rhats at 1.09,  so I am re-running it now with 8000 just to check on it). 

the new code:

data {

int M;

int N;

int P;

int K;

int cat[N];

vector[P] x[N];

int site[N];

}

parameters {

matrix[M-1,P] beta;

matrix<lower=0>[M,P] sigma;

matrix[M,P] beta_k[K];

}

transformed parameters{

matrix[M,P] beta_full;

beta_full <- append_row(rep_row_vector(1, P), beta);

}

model {

to_vector(beta) ~ normal(0,5);

to_vector(sigma) ~ cauchy(0,3);

for (m in 1:M){

  for (k in 1:K){

    beta_k[k,m] ~ normal(beta_full[m],sigma[m]);

  }}

for (n in 1:N)

cat[n] ~ categorical_logit(beta_k[site[n]] * x[n]);

}

On Thursday, 14 April 2016 16:47:28 UTC-4, Bob Carpenter wrote:

[James Savage]

Good to hear Rachael. 

A few comments:

- You probably don't need 5k iterations--In a similar model I get autocorrelation of practically 0, perfect convergence, and typically run 4-500 iterations. How many independent draws from a distribution do you need to get an estimate of the mean? 5k draws might be useful for estimating tail quantiles, but I always take these with a grain of salt, dependent as they are on the choice of distributions. But I understand you're an economist and will probably have some seminar participant telling you to run half a million iterations with thinning of 1000 or something. 

- Your priors can probably be tighter, and this will speed things up a bit. If you have a coefficient of 10 (two standard deviations on your prior), that'll basically correspond to a probability of 1 for that case. So unless you have issues with incomplete overlap, I'd use something tighter--normal(0, 1 or 2 or something). Where I mightn't do this is if you have a theoretical reason why the incomplete overlap makes sense (your cat might be an absorbing state, like default or death). 

- This specification of the coefficient covariance matrix is diagonal. If you're only interested in what happens for the sites that you're measuring, then this isn't too bad (your posterior will have unmodelled covariance). In experiments I've found that explicitly modeling covariance in params doesn't make a big difference within observed levels. But if you want to generalise to new sites or model the group-level parameters directly, you'll want to specify a richer covariance structure. If you want to do this, look at the discussion from page 216 of the manual (don't use beta[m] ~ multi_normal). 

- Now you can model profits for cat 1 and 3 using the same beta covariance matrix, tying the state model and profits model together. I'd model the negative of the losses so that you can use positively-constrained distributions for this. 

Jim

[Rachael Meager]

Hi Jim! 

- Good point, I checked and you're right on this point. The previous bad convergence was probably also this scaling problem from the K parameters. If econ people want more iters, I can do that later.

- The tighter priors would be more alarming to my econ seminar participants, but I can give it a shot anyway, since as long as the sensitivity isn't too bad they won't mind.

- My current spec does a mediocre job inferring the beta parameters from test data. Nowhere near as good as I've seen stan do in all my other models. I wonder if the reparameterising trick on 216-18 would help with that too? I have only 5 sites in my fake test data, but all the sites actually have the same betas, which the model correctly detects -- but it doesn't do a great job of estimating any of them. This is with 1000 data points per site, which is certainly enough to expect good inference. I want to fix this before I build in the modelling for profits. I suspect this is a high-correlation-in-parameters sampling problem since it gets all of them kinda-wrong in almost exactly the same way, so I will try out the reparameterisation now. 

- At some point could I send you an email to discuss how to share information between the 2 bits of the models (the multinomial bit and the profit model bit)? I have mine set up and running independently now, but seems intuitive that I should be sharing info across them.

A note for anyone else reading this thread: the model I posted in my most recent comment before this comment has the first coefficient vector beta_k[1] normalised to 1, which then makes it hard to convert the betas to category probabilities. I have now changed it to normalize beta_k[K] to 0 which means I can use standard conversion formulas. 

[Rachael Meager]

I've fixed it now! I had a bug in my simulation of the test data. The model works great. I'm calling this thread completed now. Thanks again everyone!