I am curious about performing poisson test where I have uncertainty about my count. For example, I expect to see 15 bunny rabbits per hike. On a given hike, I positively identify 19 bunny rabbits, and I think I see 4 more bunny rabbits, but I can't be sure. How do I include my maybe bunny rabbits into my poisson test? Is there a more difficult to pronounce test I need to use in this case?

Edit 1

I'm not sure this is the right approach, and I'm not sure what to do with the resulting standard error.

> install.packages(Amelia);
> library(Amelia);
> q <- data.frame("bunny-count"=c(20,19,23,22,19,22))
> se <- data.frame("se"=c(0.0625,0.0625,0.0625,0.0625,0.0625,0.0625))
> mi.meld(q, se)
[1,]    20.83333

[1,] 1.861456

One solution is to assign each sighting a probability. If these probabilities are independent, then you simply sample the cases - each according to a Bernoulli distribution with that probability - to repeatedly create multiple datasets with 19 to 23 bunnies. If they are not independent, e.g. "because if number 20 was a bunny then so was 21, but I just don't remember whether something with ears like that is a hare instead", you may need some fancier sampling scheme.

Then you analyze each of these multiple datasets in the same way you would analyze a dataset without any uncertainty in it and combine the results according to Rubin's rule. You can then do any inference you need on the combined results.

| cite | improve this answer | |
  • $\begingroup$ I found an R package named Amelia with a function called mi.meld; I think that seems to be what your'e talking about. It accepts a data frame of observations, and another of standard errors. If I count my maybe bunny rabbits as a single Bernouli trials, with p=q=0.5, then the standard error for a single maybe buny is $\sqrt{p,q}=0.5$, and four of them is $0.5^{4}=0.0625$. Then I'm thinking the mi.meld input would be a data frame with several counts, and a data frame of the same dimensions, filled with 0.0625. Does that sound right? $\endgroup$ – cjohnson318 Feb 20 '19 at 22:19
  • $\begingroup$ mi.meld may well implement Rubin's rule (Amelia = a package for MI, so should have functions for that). However, you want estimates of what you are interested in based on your certain and uncertain data. E.g. if bunnies per hour ~ Poisson(mean rate * time) and you want to know the mean rate, you would fit a Poisson model repeatedly on your (partially) imputed data (observed counts + $Y_{\text{imputed},i} \sim \text{Bernoulli}(n=4, \pi=0.5)$). You then put the estimated log mean rates with SEs into Rubin's rule and back-transform the result to a mean rate/time unit. $\endgroup$ – Björn Feb 21 '19 at 4:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.