goodness of fit for mixed model

[Linas Mockus]

Hi,

I am trying to utilize generated quantities for goodness of fit of a mixed model.  The random effect is the intercept for a lot (batch) - variable alpha.  There are several fixed effects but the one related with intercept is called coeff[1].  When I check goodness of fit using Perason's SSQ for existing lots the fit is perfect and p-value is very close to 0.5.  However, when I check goodness of fit for new lots the the fit is very bad.  Attached is stan code, R code for fit, and the fits for existing and new lots.  Can anybody explain what is happening.

Thanks for any insight,
Linas

[Luc Coffeng]

Just to check if I'm understanding your code: your model predicts counts using a Poisson, with a log-linear predictor that involves both a random intercept (alpha) and slope (gamma*diameter), where you've set the prior mean for the random slope to zero (why?)).

In the generated quantities you simulate both posterior predictive values (i.e. "counts_rep": random counts given the fitted mu for each data point) and mixed prior-posterior predictive values (i.e. "counts_rep1": random counts given independently re-sampled random effects and given the fitted fixed effects).

My first thought is there may be considerable correlation between the two random effects alpha and gamma, which you are not capturing in your posterior predictive checks. Solution would be to explicitly model alpha and and gamma from a bivariate normal distribution, and using the estimated hyperparameters of this bivariate distribution to generate alpha and gamma in the generated quantities block.

I'm not sure how to interpret the plots that you've produced. What does a data point represent? The average of the posterior? I would somehow expect "error-bars" for each data point to illustrate the variation of the discrepancy measure over all the samples of the posterior. A good starting point for me is usually to plot the replicated data (y-axis) against the observed data (x-axis), and then go on to discrepancy measures. It might also be helpful to plot the results of different lots with different colors or symbols.

Best,

Luc

[Bob Carpenter]

I have no idea about the modeling, but as far as the Stan
code goes, rather than writing this:

log_lik [i] <- poisson_log (counts [i], exp (mu [i]));

you could write:

log_lik[i] <- poisson_log_log(counts[i], mu[i]);

This takes a parameter on the log scale and does the exp transform
internally, which is less likely to cause arithmetic problems.

(I find the spaces make the code very hard to read.)

There's no reason to define the mu as a transformed parameter rather
than a local variable unless you need it in the output.

We don't have a poisson_log_rng, but we probably should.
Basically, anywhere you see an exp() in a program you should
be nervous about arithmetic behavior. Ditto an inverse() for
matrices.

- 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.
> For more options, visit https://groups.google.com/d/optout.
> <ExistingLots.jpg><NewLots.jpg><poisson.cauchy.stan><fit.r>

[Linas Mockus]

Luc,

Thanks for the insight.  Actually gamma is to model overdispersion.  I took this model from WinBUGS and it works better that negative binomial, log normal etc.

Technically total intercept is coeff[1] (fixed effect) and sigma_alpha*alpha[lot[i]] (random effect). 

Your second paragraph is correct.

I am not sure that alpha and gamma are physically correlated since alpha is between lot variation and gamma is within lot variation.  I guess I should try bivariate normal.  Can you suggest an example I could use?

The plots I derived from Bayesian modeling using WinBUGS.

I am not sure I understand the sentence "A good starting point for me is usually to plot the replicated data (y-axis) against the observed data (x-axis)".  For each draw counts_rep is different.

Thank you,
Linas

--
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/mMSr_EyK-hI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Linas Mockus]

Bob,

Thanks a lot. I will replace wit poisson_log_log.

I need mu to output for goodness of fit.

Thank you,
Linas
765-494-0344

[Luc Coffeng]

Thanks for the clarification. I know see that gamma is vector[N], not vector[NL]. Not so sure whether you need to pursue the bivariate modeling of alpha and gamma. You might want to check first whether the posterior estimates of each alpha and gamma are correlated between lots. If so, you may want to model alpha and gamma bivariately, specifying gamma conditional on alpha:

Copied and pasted from: http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_conditional_expectation

Here, X1 is gamma[i], X2 is the lot-specific alpha[lot[i]], sigma1 and sigma2 are the marginal standard deviations of gamma and alpha (both 1 in your model), mu2 is the mean of prior for alpha (zero), and rho is a correlated parameter that needs to be estimated. The implementation in your model is straightforward if you let go the Matt trick (the multiplication of a N(0,1) random value with some scale parameter). If you want to hold on to the Matt trick because your dataset is small, you would have to figure out some smart coding, which should be doable.

What I meant by "A good starting point..." was to plot the posterior of the replicated data (as mean or median + error bars for 95% credible interval) on the y-axis against the observed data on the x-axis. Telling me that you derived the plots from "Bayesian Modeling with WinBUGS" doesn't really help me understand them. Can you be more specific about what you plotted?

Best,

Luc

[Luc Coffeng]

Ignore my suggestion about the bivariate modeling. If alpha and gamma represent between and within lot variation, you certainly don't want to tangle them up. Bivariate modeling of alpha and gamma would make the model less identifiable, as fordifferent values of the rho parameter, you can arrive at exactly the same posterior density (given a uniform prior on rho) by rescaling the alpha and gamma parameters. So I guess that only leaves my question about what the plots represent!

[Linas Mockus]

Luc,

Thanks a lot.  When in R I did cor(sigma_alpha,sigma) it was 0.001.  I guess alpha and gamma are expressed using Matt's trick: alpha~N(0,1),gamma~N(0,1).  Does it make sense to go any further?

I am sorry, the correct reference for goodness of fit I am using is copied from 8.4.2 of  http://www.mbr-pwrc.usgs.gov/software/kerybook/R_WB_code.txt.

Thank you,
Linas

--

[Luc Coffeng]

I think you're fine modeling alpha and gamma as is. Forget my comments about Matt tricks and bivariate distributions. This was started by me misunderstanding what gamma represents.

The link you provided is too much text for me to go through and try to interpret. I can't tell from the code what the plotted variables represent. Can you please paraphrase and explain what you did? Also, providing a plot of the posterior (mean and 95% CI) against observed data, preferably color-coded by lot, would be very helpful to try and help understand the issues of your model.

[Luc Coffeng]

Sorry, just to clarify my last comment: I meant plotting the posterior of the replicated data against the observed data!

[Christian Roy]

If I understand correctly what Linas is trying to calculate is a Bayesian P-value. 

Basically in the generated block in stan it would look like this:

generated quantities{

vector[N] Pearson_resid;

vector[N] Deviance;

real fit;

vector[N] pred_counts;

vector[N] pred_resid;

vector[N] pred_Deviance;

real new_fit;

Pearson_resid <-  (counts-mu)/sqrt(mu);     //Pearson residuals

Deviance <- pow(Pearson_resid,2);

pred_counts <- poisson_rng (exp (mu));

pred_resid <- (pred_counts-mu)/sqrt(mu);   //Pearson residuals

pred_Deviance <- pow(pred_resid,2);

fit <- sum(Deviance);

new_fit <- sum(new_Deviance);

}

and he his plotting new_fit in function of fit.

Hope this help... but I've tried to do this from memory so there's probably a mistake somewhere.....

Chris

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

On Tue, Jul 29, 2014 at 8:23 PM, Luc Coffeng <lucco...@gmail.com> wrote:

I think you're fine modeling alpha and gamma as is. Forget my comments about Matt tricks and bivariate distributions. This was started by me misunderstanding what gamma represents.

The link you provided is too much text for me to go through and try to interpret. I can't tell from the code what the plotted variables represent. Can you please paraphrase and explain what you did? Also, providing a plot of the posterior (mean and 95% CI) against observed data, preferably color-coded by lot, would be very helpful to try and help understand the issues of your model.

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

[Luc Coffeng]

Again looking at the second plot of NewLots.jpg, I start wondering what went into the calculation of the discrepancy metric (squared residuals?) for each axis. I'm guessing that for the replicate data, you did (data_rep - mu_rep)^2, where mu_rep includes the resampled random effects (i.e. mixed prior-posterior predictive check);  and that for the observed data, you calculated (data - mu_rep)^2.

If this is so, one should expect that discrepancy metric for the observed data is always going to be much and much larger than for the replicated data. So I'm not sure (and really don't know) whether you are calculating your discrepancy metrics in the right way. I often get terribly confused when I try to think about posterior predictive checks, and it's happening again now :). Really, we want someone like Andrew Gelman to chime in at this point!

What I think is the correct way to calculated discrepancy metric T for your mixed prior-posterior predictive check is:

T_rep = (data_rep - mu_fit)^2, where data_rep includes resampled random effects
T_obs = (data_obs - mu_fit)^2

T = T_obs - T_rep

Then, if T for any data-point is often > 0, it means that the observed data are rather far away from the model. You could also sum your T_rep and T_obs per cluster and then compare whether some clusters lie within a plausible distance from the model.

[Luc Coffeng]

Thanks Chris,

This seems to confirm my point: the discrepancy metric should include the fitted value of mu for both the replicate and observed data (i.e. not mu_rep in case of mixed prior-posterior predictive checks).

Thanks,

Luc

[Linas Mockus]

Luc, Chris,

Thanks a lot.  I used Perason's SSQ since the distribution is poisson, not normal.

Anyway, attached is plot of replicated vs observed (red - 0.5 quantile, green - 0.25, and 0.975 quantiles for existing batches).  Looks like awesome fit.

The fit for new batches is terrible!!!

I suppose there is no reason to draw each lot - the clusters are easily separable.

Thank you,
Linas
 

On 7/29/2014 8:48 PM, Luc Coffeng wrote:

Thanks Chris,

This seems confirms my point: the discrepancy metric should include the fitted value of mu for both the replicate and observed data (i.e. not mu_rep in case of mixed prior-posterior predictive checks).

Thanks,

Luc

On Tuesday, 29 July 2014 17:41:24 UTC-7, Chris Roy wrote:

--
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/mMSr_EyK-hI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Linas Mockus]

fit.data$counts is observed data on x axis.  q is quantile.  The R code is:

la <- extract (fit.poisson.mix [[1]], permuted = TRUE)
q <- array (NA, dim = c (length (fit.data$counts), 3))
for (i in 1 : length (fit.data$counts)) {
    q [i, ] <- quantile (la$counts_rep [, i], probs = c (0.25, 0.5, 0.975))
}
plot (fit.data$counts, q [, 2], col = "red")
lines (fit.data$counts, q [, 1], col = "green")
lines (fit.data$counts, q [, 3], col = "green")
for (i in 1 : length (fit.data$counts)) {
    q [i, ] <- quantile (la$counts_rep1 [, i], probs = c (0.25, 0.5, 0.975))
}
plot (fit.data$counts, q [, 2], col = "red")
points (fit.data$counts, q [, 1], col = "green")
points (fit.data$counts, q [, 3], col = "green")

Thank you,
Linas
 

[Luc Coffeng]

So a common issue with posterior predictive checks for hierarchical models is that the the classic approach gives results that are far too optimistic (e.g. compared to cross-validation). When you only replicate data from the fitted mu's, you are essentially overestimating the model fit in your PPC. That's exactly what you see in the first plot: perfect correlation. No single model is really that good! The second model looks much more like something that you would expect from a cross-validation exercise, although it's hard to see which sets of green and red data-points actually represent a single observations. Error-bars would be great :), as in dx.doi.org/10.1371/journal.pntd.0002168.g004, just to blatantly quote myself. I don't have the code around though, otherwise I'd happily share it with you.

Now about that discrepancy test that you used for NewLots.jpg in your original post: did you calculate Pearson residuals relative to mu_rep or mu_fit? I.e., what "mu" are you feeding into the function in "fit.r" from your original post? You should be feeding it fitted values of mu and replicated data based on replicated values of mu (i.e. mu including new random effects).

Hope this clarifies things.

Luc

[Luc Coffeng]

As an example: here's a link to a paper that compares cross-validation to classic and mixed posterior predictive checks: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2675184/

[Linas Mockus]

Luc,

Thanks for the insight.  I think that I need to look for a better model since the full predictive distribution is too optimistic.  My observation is that for a wide range of observed values prediction is pretty much the same.  Are there any techniques to make prediction more sensitive to observed values?  My first thought was that all measurements in the lot have the same diameter and length.  What about adding random effects to diameter and length?

The plot with error bars with each lot in different color is attached.

Thank you,
Linas

On 7/29/2014 9:38 PM, Luc Coffeng wrote:

As an example: here's a link to a paper that compares cross-validation to classic and mixed posterior predictive checks: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2675184/

[Linas Mockus]

Luc,

For existing lots I used mu to calculate Pearson residuals.  For new lots I used mu_rep to calculate Pearson residuals.  Hopefully I got at least this right.

Thank you,
Linas
765-494-0090 

--

[Luc Coffeng]

Hi Linus,

Like I said, you should always be using the fitted mu to calculate your discrepancy test, also for your "new lots". So for mixed prior/posterior predictive checks, the Pearson residuals are: (data_obs - mu_fit) / mu_fit and (data_rep - mu_fit) / mu_fit, where data_rep is based on mu_rep.

"I think that I need to look for a better model since the full predictive distribution is too optimistic."

The way you calculate the residuals for your "existing lots" is too optimistic. That in itself does not mean that your model is wrong. You need to make a better comparison, which is why you should be (and are) doing a mixed prior/posterior predictive check (i.e. as above in bold). .

"My observation is that for a wide range of observed values prediction is pretty much the same.  Are there any techniques to make prediction more sensitive to observed values?"

Keep in mind that have a Poisson model. So when you simulate replicate data from mu 1.6 and 1.8, the generated values will look pretty much the same because of Poisson sampling error. Given that your model also included overdispersion, the variation in your simulated replicate data should be large! You may want to check whether the overdispersion term is really necessary in your case.

That plot looks great. For interpretability, you may want to log-transform both axes and add an identity line (y=x). Judging from the plot, it seems that the random intercepts explains most of the variation in the data (the data seem very neatly grouped by color), and within lots, there seem to be relatively little variation in predicted values. Do the length and diameter covariates really explain anything?

[Linas Mockus]

Luc,

Thanks a lot for the insight you have provided.  Attached is error bar plot with random effects for diameter and length. 

SSQ of Pearson residuals or just SSQ for new batches don't provide decent  p-value (not equal 0 or 1).  But thanks for reminding that mu_fit should be used both for existing and new batches.

I am novice and knowing how to check for overdispersion will be very helpful.

Thank you,
Linas
765-494-0090 

On 7/30/2014 11:39 AM, Luc Coffeng wrote:

Hi Linus,

Like I said, you should always the fitted mu to calculate your discrepancy test, also for your "new lots". So for mixed prior/posterior predictive checks, the Pearson residuals are: (data_obs - mu_fit) / mu_fit and (data_rep - mu_fit) / mu_fit, where data_rep is based on mu_rep.

[Luc Coffeng]

I'm not sure what you mean with "random effects for diameter and length". Previously, you had fixed effects for diameter and length. Did you change your model? Just in case: in your posterior predictive checks, you shouldn't resample the coefficients for diameter and length. If you do (and also resample the sigma's), you step into the arena of prior predictive checks (which I've never used before, but there's plenty of literature on it if you want to find out what it's for).

What do you mean by decent in "SSQ of Pearson residuals or just SSQ for new batches don't provide decent  p-value (not equal 0 or 1)"? You want your p-value to be far away from 0 and 1, right? I.e. you don't want large differences between the observed and replicated data, relative to the fitted model.

Regarding checking for overdispersion, the first thing would be to use common sense: check the scale of the overdispersion parameter. Is it close to zero? Next, you could use an information criterion to compare models with and without the overdispersion parameter (the WAIC has been suggested before elsewhere on this list). You might even be able to do some sort posterior predictive check with a clever discrepancy test to see whether the overdispersion makes sense (although I would have to think longer about that, have to run now).

Luc

[Linas Mockus]

Luc,

Tanks for providing insight.  It was very helpful.

Would you please let me know if the model fit is reasonable?  I actually added batch random effects for diameter and length (stan code is attached).

NewLots show predicted vs. observed for mixed predictive: each lot in different color, black circles - 0.5 quantile, error bar - 0.025 and 0.975 quantile (variable pred_m_...)

ExistingLots shows full posterior predictive check - too optimistic (variable pred_...)

The Bayesian p-value by usinq SSQ of Pearson residuals is also shown. p-value is 0.49 for existing lots (variable pred_SSQ) and 0.93 for new lots (variable pred_m_SSQ).

I also looked how well mean and standard deviation is predicted for a given diameter and length.  Figures are also attached.  For existing batches mean and sd are predicted perfectly.  For new batches mean is predicted quite well while standard deviation is under estimated.

Let me know your thoughts.

Thank you,
Linas
765-494-0090 

[Luc Coffeng]

Hi Linus,

Nice clean code, easy to read :).

First of all, you might want to check how reasonable the assumption is that the three random effects are not correlated. In the model section, your code says:

intercept ~ normal (0, 1);
	d ~ normal (0, 1);
	l ~ normal (0, 1);

If in fact these random effects are correlated (check by scattering the draws from the posterior against each other) - which they probably are - you get into trouble when you generate replicate data based on resampled random effects. In other words, in your mixed prior/posterior predictive checks, your resampled random effects are assumed to be independent, whereas they probably are not. Therefore, I recommend that you model the three normalized random effects through a multivariate normal. See elsewhere on this list for advice about efficiently sampling from multivariate normal distributions.

Second, I strongly recommend that you plot data and replicate data on the log-scale (the first two figures; not necessary for residuals). It's hard to see what's happening near the origin. Also, add an identity line (y = x) to each of the scatter plots. Otherwise it's hard to see what's happening, especially in the graphs *SSQ.jpg (you might want to make the scales and ranges the same for the y and x axes).

I'm not sure that "ExistingLots.jpg" is showing what it's supposed to. What do the axes represent? If it's replicated data (pred_counts) versus observed data, it should look similar in shape to "NewLots.jpg", but with the means approximately lined up along the line y=x. I would also expect to see at least some variation (error bars) due to the random Poisson samples.

I presume that in the *SSQ.jpg figures, each open circle represents a draw from the posterior? To interpret these figures, you need the line y=x and equal scale and ranges for the axes. Judging from the p-values (0.93 for new lots), something funny is going on. This could be due to the independent sampling of random effects, so you should first fix that (see above). After having done so, (and regardless of whether it solves the issue,) you might want to produce sums of squared residuals per lot and calculate p-values for those to see whether any of lots are deviating from the general pattern.

Regarding checking the means and sd's, I'm not sure I understand what you are trying to do. What did you calculate means and sd's for? Are the histograms the posterior? Of what? Is the red line the mean or sd of the data?

Best,

Luc

[Linas Mockus]

Luc,

Thanks.  I have update the code as you have suggested but really didn't see any visual effect.  WAIC is remained similar.  However I learned how to sample from a multivariate distribution and caught a typo in the stan user's manual.  I think both learnings are very valuable.

Attached is the plot in log scale with identity line for new lots (generated from pred_m_...).
The R code to generate this plot:
for (i in 1 : N) q [i, ] <- quantile (la$pred_m_counts [, i], probs = c (0.025, 0.5, 0.975))
plot (log (x), log (q [, 2]), col = "black", 
    xlab = "log (observed)", xlim = lim,
    ylab = "log (predicted)", ylim = lim)
up <- q [, 3]
low <-  q [, 1]
segments (log (x), log (low), log (x), log (up), col = cols [lot.i []])
abline (0, 1, col = "black", lwd = 2)


Attached is the plot in log scale for existing lots (generated from pred_...).
The R code to generate this plot:
for (i in 1 : N) q [i, ] <- quantile (la$pred_counts [, i], probs = c (0.025, 0.5, 0.975))
plot (log (x), log (q [, 2]), col = "black", 
    xlab = "log (observed)", xlim = lim,
    ylab = "log (predicted)", ylim = lim, type = "l")
lines (log (x), log (q [, 1]), col = "green")
lines (log (x), log (q [, 3]), col = "green")

Maybe I messed up somewhere but New* and Existing* look quite different.

For existing lots the p-value was 0.5: mean (la$pred_SSQ > la$SSQ)

For new lots p-value was 0.8: mean (la$pred_m_SSQ > la$SSQ)

I saw somewhere that folks check if observed mean and sd (in red) lies within the distribution of predicted means and sd.

Let me know your thoughts (if you still have patience - I apologize, but I am novice in this area).

Thank you,
Linas
765-494-0090 

[Luc Coffeng]

Hi Linas,

It looks like you did things right. In the full posterior predictive checks, it seems that poisson sampling error is relatively small compared to the marginal variance of the data, hence the "smooth" line in "ExistingLots.jpg". If you'd zoom in enough, you'd probably see the Poisson sampling error. But like I said, using full posterior predictive checks for a hierarchical model overestimates the model's predictive validity.

I seems that modeling the correlation of random effects did improve the predictive validity of the model, as the Bayesian p-value for the mixed PPC changed from 0.93 to 0.80, which would make me feel much more comfortable making predictions. Just curious: what were the correlations of the random effects?

To even better understand the model and the data, you may still want to evaluate the SSQ within each lot. Probably, one or a few lots are responsible for the results you see in NewLotsSSQ.jpg. My guess would be that the lot in the left bottom corner of NewLots.jpg is contributing considerably to the total SSQ.

Good luck finishing up!

Luc

[Linas Mockus]

Luc,

Looks like I am having convergence issues I didn't have when all random effects were considered independent.  Is there any way to improve convergence?

Thank you,
Linas
765-494-0090 

[Luc Coffeng]

Can you be more specific about what the issue is?