4
$\begingroup$

Let's say I want to built a confidence interval for a proportion ($p$). I have seen that many times this is made from the standard error:

$p\pm 1.96\times SE(p)$

taking SE from

$SE(p)=\sqrt{(p(1-p)/n)}$

But I guess it is not reliable if I have a small $n$ or $p=1$. So, if had 5 positive animals from 5 sampled animals ($p=1$) taken from a large population (e.g. 20000), this procedure does not make much sense for me. Shouldn't be better to calculate the Clopper-Pearson interval or Wilson Score Interval?

$\endgroup$
  • 1
    $\begingroup$ I think you've answered your own question.... $\endgroup$ – Nick Cox Nov 25 '14 at 18:25
  • 3
    $\begingroup$ I would suggest moving to a Bayesian solution, which can handle this sort of extreme event. $\endgroup$ – Xi'an Nov 25 '14 at 20:18
  • 3
    $\begingroup$ Some of the references in the wikipedia article on binomial proportion confidence intervals are worth reading. The article itself gives a reasonable overview of the choices. If your $n$ was a little larger I'd suggest the rule of three as something reasonably simple. $\endgroup$ – Glen_b -Reinstate Monica Nov 25 '14 at 21:56
  • 3
    $\begingroup$ I disagree with the first part of @heropup's advice. Clopper-Pearson doesn't work so well as Jeffreys or Wilson methods and even the Jeffreys method has a frequentist interpretation. See e.g. projecteuclid.org/download/pdf_1/euclid.ss/1009213286 $\endgroup$ – Nick Cox Nov 26 '14 at 17:08
  • 3
    $\begingroup$ You might want to check out Agresti's work: Agresti, A. and Caffo, B. (2000). Simple and effective confidence intervals for proportions and differences of proportions result from adding two successes and two failures. The American Statistician, 54(4):280–288. —and— Agresti, A. and Coull, B. A. (1998). Approximate is better than “exact” for interval estimation of binomial proportions. The American Statistician, 52(2):119–126. $\endgroup$ – Alexis Nov 26 '14 at 17:40
3
$\begingroup$

I'm going to presume to rephrase the question in two ways. If you don't like the new question please let me know in the comments and I can adjust.

Given:

  • An estimated total population of 20,000 which exist in two states appropriate to a binomial model (or even a SIR model)
  • That five samples taken at (presumably random) from this population were "positive"

Determine:

  • What is the likely rate of occurrence in the population? Include CI if possible.
  • If rate cannot be determined, what bounds can be placed on it?

Requirements:

  • If a programming language like "R" is used, provide code.
  • Describe both the method and some justification for it.

(Working)

Thoughts on approach:

  • I don't know the rate, so I'm going to use resampling to find the rate where 50% the time a random draw will yield 5/5 positive samples.
  • After I do that, I will use that rate and resampling to determine the variation of results. I will lock in the rate, redraw 1000 times, and look at the distribution of results.
  • I'm assuming that the sick creatures are tagged, and you aren't just putting sick fish back into the river. This means sampling is done without replacement.
  • This is about disease in a natural population so an SIR model is much more appropriate to the "physics" of the phenomena. I need to think about this and maybe do something eventually.
  • I don't believe in "pristine" results. It allows only reproduction of exact results, but is i.m.o. dishonest about the path so it does not allow application of the sometimes inelegant actual path and denies access to innovation in two ways. The broken road gets to all destinations that along roads that are similar to the one described. When the student is faced with a real and "non-pristine" problem they can become dis-spirited and abandon the correct but inelegant approach. I hope to capture more of the road without too much of the rocks.

Code:

#housekeeping
rm(list = ls())
plot.new()
dev.off()

#make candidate rates

#how many steps in rate
myratecount <- 100+1; 
myrates <- seq(from=0, to=1, length.out=myratecount);
myrates

#how many repeats to get values
myloops <- 1000

#total population
mypop <- 20000

#predeclare before loop
mu <- matrix(nrow=myratecount,ncol=myloops)

for (i in 1:length(myrates)) {

  print(i)

  #make distribution samples
  y <- rbinom(n=mypop*30,prob=myrates[i],size=5)/5

  for (j in 1:myloops){

    #sample from the distribution
    #draw samples
    y2 <- sample(y, size=5, replace = F, prob = NULL)

    #find mean
    mu[i,j] <- mean(y2)

  }

}

#compute ensemble mean and sd
M <- rowMeans(mu,dims=1,na.rm=F)
S <- apply(mu,1,sd)

#linear fit of central tendency
fit <- lm(M~myrates)
M1 <- predict.lm(fit)

#1-sigma UCL but not non-physical
# S1 <- M1+2.132*2*S
S1 <- M1+2.132*2*S
S1[which(S1>1)] <-1

#1-sigma LCL but not non-physical
S2 <- M1-2.132*2*S
S2[which(S2<0)] <- 0


#Make Plot
plot(myrates, M,
     main="Population rate vs. Sample Rate",
     xlab="Population rate",
     ylab="Mean Sample rate")
lines(myrates,predict.lm(fit), col="Green")
lines(M1,S1,col="Red")
lines(M1,S2,col="Red")


#annotate plot
grid()

leg.txt <- c("Data for column means", 
             "Line fit of column means",
             "2.132*2-Sigma UCL",
             "2.132*2-Sigma LCL")

legend(x=0,y=1,leg.txt,
       col=c(1,3,2,2),
       pch=c(1,NA,NA,NA),
       lty=c(NA,1,1,1))

ind <- which(S1==1,arr.ind=T)
ind2 <- which(ind==min(ind))
ind2
S1[ind[ind2]]
M1[ind[ind2]]
points(M1[ind[ind2]],S1[ind[ind2]],
       pch=20,col="Red")
points(1,1,
       pch=20,col="Red")

points(M1[ind[ind2]],S1[ind[ind2]],
       pch=20,col="Red")
points(mean(c(M1[ind[ind2]],1)),mean(c(S1[ind[ind2]],1)),
       pch=20,col="Blue")

print(mean(c(M1[ind[ind2]],1)))

The result of this is the following plot. Remember that the red lines are LOWER than textbook 3-sigma or 3.5-sigma Upper control limits (UCL) or Lower control limits (LCL).

enter image description here

If this was a perfect model then you could have only a 96% rate in your population and still get 5 samples that were positive. The (HUGE) problem with that is that we need to have a scale of sigma appropriate to our population - we are trying to draw a box around where the data can live.

If I look at a t-table for 5 samples (nu=4) and find the 95% CI scale then I get a value of 2.132. This means that we can expect the majority of variation in data to occur with a window of +/- 2.132*Sample_standard_deviations around the central tendency. If we scale the sample_standard_deviation by this value then it yields a truer 1-sigma boundary. Lets make the truer 2-sigma boundary by multiplying 2.132 by 2. If I modify the code (not shown) so that the UCL and LCL are mean +/- 2.132*2*sigma then the following results.

enter image description here

The really surprising (for me) consequence here is that the scaled 2-sigma UCL is at 100% below a population mean of 60%. A few numeric lookups gives the value as 0.57546, or about 57%. It is not entirely out of the realm of consideration that though you got 5/5 positive, the actual population rate is around 57%. If you have to have a highest plausible rate then pick 100%. If you have to pick a lowest then pick 57%. The midpoint between 100 and 57%, one that is equidistant to the extremes of error is about 78%. If you have to make a guess then 78% is what I would put as the least likely erroneous rate, all else being equal.

Thought:

It is worthy to contemplate "what if I had one more sample". If you have 6 samples (aka nu=5) then your t-scale on variation goes from 2.132 to 2.015. Your minimum plausible threshold goes from 57% up to 66%. That is about a 9% change in estimate for a single sample. So now you can relate an increase in sample size by one to a value in the estimate.

I have to ask about selection bias. If you have a buddy who is a below-average hunter, and he is tranquilizing coyotes and taking blood samples, and this is where your data comes from then there might be a selection bias issue. He is going to have an easier time going after the weakest/slowest. If that were the case, and you could get estimates of the pack-sizes then you might be able to do some sort of extreme-value-distribution related analysis. It would change each positive from being 1/1 to 1/pack-size. The 5 samples would, hopefully, come from different packs.

...

SIR thoughts to come later .. possibly never. It is too complex, and I don't know anything about the nature of the creature or the nature of the "positive". You could be looking at hormone levels in fish, counting parasites on shrimp, looking at bark-beetles in southwest pine trees, or Ebola patients in Mali. A general model to blindly handle that is ... too big for me to work on right now.

$\endgroup$
  • $\begingroup$ It's a very nice answer. Supositions you've done are right. In my work, these samples often are collected from hunters, so a type of bias could be acting. But the response is the presence/abscence of antibodies, so I've presumed no bias with the response (it may be different if I'd check exposure factors). Regarding SIR, I guess it should be the next step, once you've detected a condition worths to model and it'll also require more data. Overall, I feel that confidence intervals, in the way that I've often seen them reported in literature, are not always very reliable. $\endgroup$ – BonScott Nov 27 '14 at 15:26
  • $\begingroup$ Analytic forms are always right, unless they are not. The trick then is to get a library of methods that are good at testing when the analytic forms apply and when they do not. $\endgroup$ – EngrStudent - Reinstate Monica Nov 27 '14 at 18:54

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.