# multiplication of probabilities of independent events with confidence intervals in archaeology

I am an archaeologist and I am devising rules for attributing objects to a certain region based a combination of region-specific traits. Particularly I would like to estimate probability that objects displaying a combination of two traits originate from the certain region. I consider (arbitrarily, I admit) the traits I observe to occur independently and I establish their occurrence rates from the limited number of securely attributable objects, presently known (hence, I cannot just increase the sample size). Say, there are 34 objects of known origin displaying trait A and 29 of these are from region X. And there are 84 objects of known origin displaying trait B and 75 of these are from region X. It should be noted that very few objects appear in both samples (i. e., traits A and B but rarely co-occur, whence my arbitrary judgement on their independence). Using the Wilson's test with Yates' continuity correction (the prop.test function in R) I estimate that the probability that trait A occurs outside region X is between 0.06 and 0.32, and the probability that trait B occurs outside region X is between 0.05 and 0.20 at 95% confidence level.

prop.test(34-29, 34, conf.level = 0.95, correct=TRUE) # for trait A
prop.test(84-75, 84, conf.level = 0.95, correct=TRUE) # for trait B


Now I want to calculate the probability that a combination of traits A and B occurs outside region X. Would it be correct to say that it may be defined by the products of the lower and upper confidence limits of the independent traits? Would it be correct to say that P(A ∩ B) lies within the confidence limits 0.06*0.05=0.00 and 0.32*0.20=0.06 and that thus an object displaying both traits A and B stems from region X with probability P(A ∩ B) > 1-0.06=94% at 95% confidence level (with the assumption that the two traits are distributed independently)? Or do I miss something? Does the confidence level change when I multiply the probabilities?

UPDATE In a comment below it is suggested that I'd better apply another method instead. However, I reckon that if there is no mistake in my calculations, they could explain in very simple terms why the combination of two traits on a single monument provides sufficient evidence and occurrence of only one trait does not.

UPDATE 2: From what I observe applying the method suggested by @jwimberley to my data, the confidence limits of the product seem be narrower than the products of the original lower and upper confidence limits. Would it be safe to say that this should always be the case? (Or in other words that the confidence level corresponding to the product of the original confidence limits would be => 95%)?

• I bet you would be happier with a comprehensive, multivariate model to predict region. This would potentially save you work, especially if you add variables, and concisely account for major dependencies that exist between predictors. – rolando2 Jan 17 '17 at 23:46
• Thank you for your suggestion. Indeed, I could as well use a binomial generalized linear model and get prediction probabilities and standard errors for all possible combinations of traits. However, I wonder if there is an answer to my question. – greenb Jan 18 '17 at 10:31
• For basically I need to explain in very simple terms how my method works and why the combination of two traits on a single monument provides sufficient evidence and occurrence of only one trait does not. – greenb Jan 18 '17 at 11:51
• Have you verified that these traits are indeed independent? – jwimberley Jan 19 '17 at 14:49
• Possible duplicate: math.stackexchange.com/questions/765887/… Although the answer there perhaps does not completely address this question – jwimberley Jan 19 '17 at 15:09

Leaving aside some of the questions addressed in the comments (is there a better procedure? have you verified that the traits are indeed independent), here is a direct answer: No, the confidence interval of the product is not equal to the product of the separate confidence intervals. The first confidence interval essentially states that $$\int_{p_A^L}^{p_A^U} f_A(p_A) \, dp_A = 0.95$$ and the second that $$\int_{p_B^L}^{p_B^U} f_B(p_B) \, dp_B = 0.95$$ Note that these are integrals of probability densities of probabilities!. Now, what's the density of the product probability $\rho = p_A p_B$? From a standard result, this is $$f_{P}(\rho) = \int_0^1 \frac{f_A(p) f_B(\rho/p)}{p} \, dp$$ The CDF of this distribution is $$F_{P}(\rho) = \int_0^1 \frac{f_A(p)}{p} \left( \int_0^\rho f_B(p'/p) \, dp' \right)\, dp \\ = \int_0^1 f_A(p) F_B(\rho/p) \, dp$$ It does not follow that $$\int_{p_A^L p_B^L}^{p_B^U p_B^U} f_P(\rho) \, d\rho = 0.95$$ which is what you estimate. What you must solve is $$\int_{p_P^L}^{p_P^U} f_P(\rho) \, d\rho = 0.95$$ for the new confidence interval bounds $p_P^L$ and $p_P^U$

I'm not sure if there is a model-independent way of estimating this. Here's one idea: some probability confidence intervals are based on methods modeling the densities $f_A$ and $f_B$ as normal distributions or beta distributions. I understand (I believe) that Wilson's test is based on a frequentist method with a normal approximation; an alternative is a Bayesian method using Beta distributions. If you obtain the parameters of these normal distributions or beta distributions, you can analytically (unlikely) or numerically (more likely) compute the above integrals and find your confidence interval.

For example, the binom package in R performs Bayesian estimation of binomial proportions with binom.bayes.

> print(A <- binom.bayes(34-29, 34, conf.level = 0.95))
method x  n shape1 shape2      mean      lower     upper  sig
1  bayes 5 34    5.5   29.5 0.1571429 0.04866337 0.2771654 0.05
> print(B <- binom.bayes(84-75, 84, conf.level = 0.95))
method x  n shape1 shape2      mean      lower     upper  sig
1  bayes 9 84    9.5   75.5 0.1117647 0.04949479 0.1793428 0.05


We care most about the shape parameters; with these we have $f_A$ and $f_B$ and can numerically compute $f_P$. This is easiest via MC integration:

dbetaprod <- Vectorize(function(x,alpha1,beta1,alpha2,beta2) {
p <- rbeta(1000,alpha1,beta1)
mean(dbeta(x/p,alpha2,beta2)/p)
},"x")

pbetaprod <- Vectorize(function(x,alpha1,beta1,alpha2,beta2) {
p <- rbeta(1000,alpha1,beta1)
mean(pbeta(x/p,alpha2,beta2))
},"x")


You can then plot the density and CDF to find the confidence interval:

It looks like the 95% confidence window upper bound is about 0.04 (4%). You can modify this procedure to suit your specific needs.

• Thank you so much for your extensive answer. From what I observe applying the method you suggested to my data, the confidence limits of the product seem be narrower than the products of the original lower and upper confidence limits. Would it be safe to say that this should always be the case? (Or in other words that the confidence level corresponding to the product of the original confidence limits would be => 95%)? – greenb Jan 20 '17 at 9:12
• I don't think it would be safe to say this, no. – jwimberley Jan 20 '17 at 12:45