modeling normal measurement noise for counts in stan

[Linas Mockus]

Hi,

Does anybody have experience with modeling measurement noise for count data in stan.  Assuming that true count data from observation to observation are distributed according to poisson (perhaps with the addition of overdispersion parameter as in Gelman & Hill) measuring instrument adds some normal additive/multiplicative/log normal noise.  It could be expressed as:
true_counts[i]~poisson(lambda[i])
observed_counts[i]~normal(true_counts[i],sigma)

Problem with this is that Stan doesn't allow integer parameters true_counts.  observed_counts are real since the instrument returns real values (current).

I appreciate any comments on this.

Thanks,
Linas

[Luc Coffeng]

Hi Linas,

To clarify: you say that the observed_counts data are reals, rather than integers. Does this mean that the data are counts per some unit (volume, time, etc)? Similarly, does it perhaps mean that your instrument is actually counting something but returns an expected count as a result of some regression it is performing internally (e.g. as in qPCR)? Either way, I would start with a log-normal distribution for the observed_counts (which might actually be rates). It would be great of course if your instrument would actually return the data that it's basing its "result" on in the first place. That way, you could incorporate all that information in your Stan model.

What do you need the true_counts for? Is it ok to just produce them in the generated quantities block or do you need them somewhere else in your model? In the latter case, it sounds like you need to marginalize over all possible outcomes of true_counts, which doesn't sound pretty.

Best,

Luc

[Bob Carpenter]

On Jul 3, 2014, at 1:08 PM, Linas Mockus <linasm...@gmail.com> wrote:

> Hi,
>
> Does anybody have experience with modeling measurement noise for count data in stan. Assuming that true count data from observation to observation are distributed according to poisson (perhaps with the addition of overdispersion parameter as in Gelman & Hill) measuring instrument adds some normal additive/multiplicative/log normal noise. It could be expressed as:
> true_counts[i]~poisson(lambda[i])
> observed_counts[i]~normal(true_counts[i],sigma)

I'm confused about the data here --- is the observed_counts
variable continuous?

> Problem with this is that Stan doesn't allow integer parameters true_counts. observed_counts are real since the instrument returns real values (current).

Right --- there's no way to code this up directly. And it's
probably going to be hard to marginalize out.

> I appreciate any comments on this.

Pardon my questions here --- I'm not an expert on these models.

If the observed counts are actually continuous, then why use
a count based model at all?

If the observed counts are discrete, then why not just use
an overdispersed model like negative_binomial --- that also gives
you one parameter for the amount of extra dispersion.

An alternative would be to roll your own overdispersed model
by doing something like this:

true_expected_counts[i] ~ lognormal(lambda[i], sigma);

observed_counts[i] ~ poisson(true_expected_counts[i]);

- Bob



> Thanks,
> Linas
>
> --
> 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.

[Linas Mockus]

Hi Luc,

Observed counts are particle counts in liquid per volume.  I think the idea is to pass the light and measure the current.  Then based on the current reading and instrument calibration (assuming that particles are spherical) infer particle counts.  I always thought that count data is best represented by poisson (with extension to overdispersion) and not log normal.  In fact poisson with overdesispersion works very well (Bayesian p-value is close to 0.5, Pearson residuals look random, histogram of mean and standard deviation is replicated - actual mean and sd are in the middle of replicated mean and sd histograms - basically I did lot of model checking suggested by Gelman).  However, when I start predictions I see large uncertainty. Therefore I thought about improving the model by adding additional layer for instrument.

Technically true counts are not needed, I just wanted to separate count variability between samples and measurement error.

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

[Luc Coffeng]

Hi Linas,

Thanks for the clarification. Based on what you are saying, I would argue that your data are not counts but concentrations (i.e. rates), and therefore Poisson, negative binomial or whatever other discrete distribution you can think of are not formally appropriate. In my experience, concentrations/rates are often treated as following a log-normal distribution, with sigma representing your tool's measurement error on the log-scale. If you think there is overdispersion, I suggest you fit a log-T-distribution or a hierarchical model, depending on the setup of your experiment. The large uncertainty you are getting in your estimates probably comes from your Poisson distribution (sigma^2 = mu). BTW, as you may know, the measurement error of lab instruments in general is usually not constant, but relatively higher (beyond log-normal) for very low (and sometimes high) values of observations. Also, an important cause of so-called "overdispersion" or "zero-inflation" is actually that most instruments have a minimum detection threshold. You can take account of this by modeling observations under the detection limit ("zeroes") as right-censored observations (using the CDF given expectation, measurement error, and the detection limit (assuming you know this)).

Best,

Luc

[Linas Mockus]

Luc,

Thanks.  I kind of agree that concentration might be a good alternative - at least worth trying.  The data is longitudinal - sampled every few months - thus I am using varying slopes varying intercepts model (data comes from several batches).  I hope that lognormal will reduce the prediction uncertainty.  

Thank you,
Linas

On 7/3/2014 5:03 PM, Luc Coffeng wrote:

Hi Linas,

Thanks for the clarification. Based on what you are saying, I would argue that your data are not counts but concentrations (i.e. rates), and therefore Poisson, negative binomial or whatever other discrete distribution you can think of are not formally appropriate. In my experience, concentrations/rates are often treated as following a log-normal distribution, with sigma representing your tool's measurement error on the log-scale. If you think there is overdispersion, I suggest you fit a log-T-distribution or a hierarchical model, depending on the setup of your experiment. The large uncertainty you are getting in your estimates probably comes from your Poisson distribution (sigma^2 = mu). BTW, as you may know, the measurement error of lab instruments in general is usually not constant, but relatively higher for very low (and sometimes high) values of observations. Also, an important cause of s0-called "overdispersion" or "zero-inflation" is actually that most instruments have a minimum detection threshold. You can take account of this by modeling observations under the detection limit ("zeroes") as right-censored observations (using the CDF given expectation and measurement error).

Best,

Luc

[Luc Coffeng]

I think the issue here is appropriately quantifying the uncertainty in predictions rather than "reducing" it. If there is a lot of uncertainty due to measurement error from your instrument, this should be propagated into predictions (assuming that by predicting you mean that you are trying to predict new data points). You could use posterior predictive checks (or better, cross-validation, if possible) to see how well the model actually predicts the data. Then you can also better decide whether you like the poisson or log-normal better (I'm suspecting the log-normal does a better job).

[Bob Carpenter]

I'm sort of channeling for Andrew here, but I just
love all this discussion of model checking!

And I'm a big fan of adding measurement error models,
especially if you can tailor them to known properties of
the instrument.

And a big +1 on appropriately quantifying the uncertainty.

- Bob

On Jul 3, 2014, at 8:01 PM, Luc Coffeng <lucco...@gmail.com> wrote:

> I think the issue here is rather appropriately quantifying the uncertainty in predictions than "reducing" it. If there is a lot of uncertainty due to measurement error from your instrument, this should be propagated into predictions (assuming that by predicting you mean that you are trying to predict new data points). You could use posterior predictive checks (or better, cross-validation, if possible) to see how well the model actually predicts the data. Then you can also better decide whether you like the poisson or log-normal better (I'm suspecting the log-normal does a better job).
>
>
> --

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

[Linas Mockus]

Luc/Bob, many thanks.

I did quick check and log-normal predictions are somewhat similar to poisson (in case of poisson I had to round the data).  In fact Pearson residuals are larger as well as waic than in poisson (in poisson I have used overdispersion parameter as in Gelman & Hill).  Bayesian p-value remains similar although nice round shape in graphical posterior predictive check for sum of squared Pearson residuals became oval for log-normal.  Mean and sd prediction histograms remained very similar.  I wonder if log t distribution could do a better job.  What df should I start with? 

I agree with Bob - measurement error model may help.  It looks that now measurement errors are incorporated in randomness of counts.  I wonder what would happen if we separate them. Probably errors are caused by particles not being oval (while instruments thinks they are oval), sample preparation error etc.

Thank you,
Linas

On 7/3/2014 8:01 PM, Luc Coffeng wrote:

I think the issue here is rather appropriately quantifying the uncertainty in predictions than "reducing" it. If there is a lot of uncertainty due to measurement error from your instrument, this should be propagated into predictions (assuming that by predicting you mean that you are trying to predict new data points). You could use posterior predictive checks (or better, cross-validation, if possible) to see how well the model actually predicts the data. Then you can also better decide whether you like the poisson or log-normal better (I'm suspecting the log-normal does a better job).