# Effect of uncertainty in a parameter of a probability model on decisions

I am exposing several experimental organisms to a toxin in the food. If the dosage is high enough, they will die. I am examining the number that are dead after a day. For practical reasons, the dosage they get is not directly under my control. I have decided to model the number that are dead after one day as a binomial distribution, Bin(n,p).

In order to figure out what p is for my model, I did an experiment with 300 organisms and counted the number that had died after a day.

In the experiment of interest, I have exposed them to a potential enhancer of the the toxic effects. I would like to know if substantially more of the organisms have died. For a particular result m, I can easily compute the probability of m or greater. If this is less than 1% then I accept that the putative enhancer works otherwise I reject the hypothesis that the putative enhancer is an actual enhancer of the toxic effects.

I am wondering if my analysis of the putative enhancer needs to take into account the uncertainty in the value of p. If so, how would I do it? (I am good with R and able to program if necessary.)

Update: I just want to emphasize that I only care if my enhancer works i.e. more organisms die than I'm expecting. Any answer should leave out cases where fewer than expected organisms die.

• Is the exposure dose going to be the same across your conditions? Otherwise, you will have level of dosage as a confound. – Michelle Jan 30 '12 at 22:58
• @Michelle They are both fed from the same batch of food. I can't control how much of the food is consumed, but if two organisms consume the same amount, they will get the same dosage. – Henry B. Jan 31 '12 at 2:28
• If they are social animals like rats, you may wish to consider separate caging. If they are bacteria, meh. – Michelle Jan 31 '12 at 5:49

A very well described question. In mathematical terms we have two binomial populations:

$$x|npI\sim bin(n,p)$$

$$y|mqI\sim bin(m,q)$$

Where $x$ is the number of deaths in the control group with $n$ observations, $y$ is the number of deaths in the group of size $m$ which received the enhancer, and $I$ is the prior information - which includes the binomial assumption. The hypothesis is that $q>p$. The alternative which is most obvious is $q\leq p$ (i.e. we are not quesitoning the binomial model). I presume your "naive" approach is to consider $\hat{p}=\frac{x}{n}$ as fixed and then test if $q>\hat{p}$ using the cummulative binomial distribution.

To do the test and account for the uncertainty in $p$, the best approach is to integrate them out using a prior distribution which describes what is known about them (independently of what the data tells you about them). To do this we require the joint posterior for $p,q$.

$$p(pq|xynmI)\propto p(pq|I)p(xy|pqnmI)$$

Now I think it is reasonable to suppose the likehood factors, as you have assumed a binomial distribution, which is conditionally independent. So we have:

$$p(xy|nmpqI)\propto p^{x}(1-p)^{n-x}q^{y}(1-q)^{m-y}$$

Now because you have a large $n$ then the particular prior used won't matter unless it is quite strong. Assuming this is the case, then we can use a uniform prior $p(p|I)=1$. However we also require the conditonal prior $p(q|pI)$. This is likely to be non-uniform, but as before unless it is on the order of $m$ (i.e an $m$ degree polynomial) it won't be able to greatly affect the result. So i will use a flat prior also. This means the posterior is proportional to the likelihood which is is the product of two beta distributions. The probability we are after is given by:

$$Pr(q>p|xynmI)=\frac{ \int_0^1\int_p^1 q^{y}(1-q)^{m-y}dq p^{x}(1-p)^{n-x} dp }{ \int_0^1\int_0^1 q^{y}(1-q)^{m-y} dq p^{x}(1-p)^{n-x} dp }$$

$$=\frac{(n+1)!}{x!(n-x)!} \int_0^1\left[\frac{(m+1)!}{y!(m-y)!}\int_0^{1-p} (1-q)^{y}q^{m-y}dq\right] p^{x}(1-p)^{n-x} dp$$

The term inside the brackets can be converted to a sum via repeated use of integration by parts (see e.g. wikipedia), so we can write:

$$Pr(q>p|xynmI)=\frac{(n+1)!}{x!(n-x)!} \int_0^1\left[\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!}p^{k}(1-p)^{m+1-k}\right] p^{x}(1-p)^{n-x} dp$$ $$=\frac{(n+1)!}{x!(n-x)!}\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!} \int_0^1p^{x+k}(1-p)^{n-x+m+1-k} dp$$ $$=\frac{(n+1)!}{x!(n-x)!}\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!}\frac{(x+k)!(n+m+1-x-k)!}{(n+m+2)!}$$

We can arrange this into 3 combinatorical functions (bit like hypergeometric pdf), and show how the uncertainty affects the answer as follows:

$$Pr(q>p|xynmI)=\frac{1}{{n+m+2 \choose n+1}}\sum_{k=0}^{y}{x+k \choose k}{n-x+m+1-k \choose m+1-k}$$

I think there is an intuitive way that one could reason this out as a good test from a frequentist point of view, but I can't think of it myself. For the "sanity check" note that it is a strictly increasing function of $y$ - any additional death in the enhanced group results in an increased probability for $q>p$. Additionally, because it is a Bayesian test, you can use it sequentially, and continue testing the enhanced cases until the probability is sufficiently high or low. Further, for large $a$ and fixed $b$ we have ${a+b\choose b}\approx \frac{a^b}{b!}$, and so as $n,x$ become large such that $\hat{p}=\frac{x}{n}\approx p$ we get:

$$Pr(q>p|xynmI)\approx\frac{(m+1)!}{n^{m+1}}\sum_{k=0}^{y}\frac{x^k}{k!}\frac{(n-x)^{m+1}}{(m+1-k)!}=\sum_{k=0}^{y}{m+1\choose k}\hat{p}^k(1-\hat{p})^{m+1-k}$$

Which is basically the "naive" approach just by plugging in the estimate (and another nice equivalence between bayes and frequentist methods). For your case where $n=300$ using open office on my pc had no problem calculating the choose functions with $m=300$ also (even though it can't calculate above $170!$). So there isn't any need to use the approximation, apart from gaining an intuitive understanding the result. They are both simple sums.

When $n>>m$ and $x>>y$, then the binomial plug-in method will work fine. When $m$ is also large, you are probably better off using the normal approximation to the binomial sum above, and then you get the standard test for difference in normal means from two independent populations - just compare the difference in proportions to the standard deviation of the difference.

Note that if you do wish to use a non-uniform prior for the test, a computationally easy way is to assign two independent beta distributions, so we have $p(pq|I)\propto p^{a_p-1}(1-p)^{b_p-1}q^{a_q-1}(1-q)^{b_q-1}$ (with all parameters integers) and the test becomes:

$$\frac{1}{{n+m+a_p+b_p+a_q+b_q-2 \choose n+a_p+b_p-1}}\sum_{k=0}^{y+a_q-1}{x+k+a_p-1 \choose k}{n-x+m-k+b_p+a_q+b_q-2 \choose m-k+a_q+b_q-1}$$

UPDATE

In response to the comments, it is quite simple to write an R function to calculate the above probability as follows:

prob.q.gt.p <- function(x,y,n,m,a_p=1,b_p=1,a_q=1,b_q=1){
prob <- 0
sum.to <- y+a_q-1
for(k in 0:sum.to){
prob <- prob + choose(x+k+a_p-1,k) * choose(n+m+b_p+a_q+b_q-2-x-k,m+a_q+b_q-1-k)
}
prob <- prob / choose(n+m+a_p+b_p+a_q+b_q-2,n+a_p+b_p-1)
return(prob)
}
#example run - should return value of 0.0336
prob.q.gt.p(20,0,300,50)


Now I don't know what observed values you intend to use, except that $n=300$. Now suppose $x=20$ ($20$ deaths) for death rate of $6.7\text{%}$, and also that $m=50$ (expose $50$ new organisms to increased toxin). The table below shows the probability that $q>p$ for various values of $y$

$$\begin{array}{c|c}y & \frac{y}{m}\times 100\text{%}& Pr(q>p|y,x=20,n=300,m=50,I)\\ 0 & 0\text{%} & 0.0336 \\ 1 & 2\text{%} & 0.1424 \\ 2 & 4\text{%} & 0.3238\\ 3 & 6\text{%} & 0.5308 \\ 4 & 8\text{%} & 0.7126 \\ 5 & 10\text{%} & 0.8433 \\ 6 & 12\text{%} & 0.9231 \\ 7 & 14\text{%} & 0.9658 \\ 10 & 20\text{%} & 0.9982 \\ 15 & 30\text{%} & >0.9999 \end{array}$$

• This sounds like what I want, but is a lot more complicated than I'm able to understand. Is there a computational way to get at this? – Henry B. Feb 3 '12 at 14:44
• The easiest way is to compute the expression $Pr(q>p|xynmI)$ for various values of $x,y,n,m$. I'll add a table if you like. – probabilityislogic Feb 3 '12 at 21:36
• Thanks for this answer I asked another related question here, stats.stackexchange.com/questions/22288/… – Henry B. Feb 5 '12 at 6:10

Yes, the analysis needs to take into account the uncertainty in the proportion p. You can do that by arranging the data as a 2×2 contingency table and applying Pearson's chi-squared test or Fisher's exact test.

• +1 Other tests can be found by Googling "test of proportions". – whuber Jan 30 '12 at 21:40
• My understanding of the Fisher exact test is it computes the probability of what I observed then it adds up all cases as likely or less likely than my observations. The problem is I only care if my enhancer works i.e. more organisms die than I'm expecting. I would assume that standard Fisher exact test adds up all unlikely cases including the cases where much fewer organisms die than expected. If fewer die then I don't care. – Henry B. Jan 31 '12 at 2:21

You might want to also consider estimating the variance under a normal distribution approximation. You have high n's in your samples and if your proportions aren't too extreme that will allow you to generate standard deviations, and even confidence intervals. The variance in that case would be N*p*(1-p). Googling "binomial distribution confidence interval" will get you some good links.