|
I posed the question:
A client has taken a random sample of 108 (n) from a population of 150
(N). Out of this 108 there are only three which are positive. Under
normal circumstances you would use the normal approximation to the
binomial (correcting for sampling from a finite population) to get a
confidence interval, but with (np<5) in this case it is no longer
valid. All comments or suggestions would be welcome.
Many thanks to all who took time to reply - I enclose all comments
below.
Dear Alan,
I suggest the use of a bayesian approach. The prior distribution is
binomial (as the positive-negative outcome can be regarded as a
bernoulli
trial), and thus the posterior distribution will be beta, from which an
appropriate confidence interval can be calculated. This method will
never
result in negative probabilities.
with regards
Karoly Acsai
[log in to unmask]
The software StatXact does exact confidence interval calculations
for binomial data. There are some useful charts for binomial
confidence intervals in Elementary Statistics Tables by H R Neave
From the chart I estimate the 95% confidence interval to be
approximately (0.005 , 0.08)
Dick Brown
[log in to unmask]
Dear Alan,
I think it is possible to look at previous messages that where sent to
allstat. A similar query as yours has already been submitted to the
list.
You should have a look at Volker Knecht's messages on 29th May and 31st
May
2002 and Luigi Santoro's summary on 12th February 2002.
Regards,
Eve
[log in to unmask]
Hi Alan,
I tend to use the following reference.
Leslie Daly (1992). Simple SAS macros for the calculation of exact
binomial
and poisson confidence limits. Comput.Biol.Med. Vol22, No 5, pp351-361
The confidence interval for p=0 is [0, 1-(a/2)^(1/n)] = [0, 0.03358] i.e
your break even point is 4 positive samples. Infact providing n >= 8 the
break even point is 4 which tends to make sense for the np=5 rule of
thumb.
r = 0 lower = 0
upper = 1-(a/2)^(1/n)
r=n lower = (a/2)^(1/n)
upper = 1
else lower = 1- Beta( 1-a/2 ; n-r-1 ; r )
upper = Beta (1-a/2 ; r+1 ; n-r )
where Beta is the Beta distribution.
Dave.
[log in to unmask]
For x positives, observed in 108 draws
from a population of 150 with p positives;
you can carry out an exact calculation of
Pr[X=x|p], the probability of
drawing x positives.
Using combination theory,
Pr[X=x|p] = xCp*(108-x)C(150-p)/108C150
where aCb is the number of combinations
for choosing a from b; aCb = b!/a!*(b!-a!).
These probabilities
can be used to determine an interval in several ways.
1) Classical approach:
Carry out a test to reject or accept each p in turn.
Where Pr[X<=x|p] <0.025 (or Pr[X>=x|p] <0.025),
p is in the lower (or upper) rejection region, and
can be excluded from the Confidence Interval (CI).
2) Bayesian approach:
Provide weak prior probabilities of p
taking any of a range of values
up to 150 (eg Pr[P=p] = 1/151 for all
0<=p<=150). Use Bayes' Theorem to estimate
posterior probabilities Pr[P=p|X=3]. Establish
the Credibility Interval (CrI)
using the same approach as (1).
3) Likelihood ratio:
Take p_max, the point estimate of p as
the value with the maximum likelihood Pr[X=3|p].
So, Pr[X=3|p_max] >= Pr[X=3|p] for all p.
Work out the likelihood rato for all p of interest.
LR = Pr[X=3|p] / Pr[X=3|p_max]
Accept p as within the likelihood interval if
LR is "large enough". As I remember,
(and without checking), exp(-2) is often considered large
enough, as with Normal errors, it gives an interval of +/- 2SD.
[log in to unmask]
Use a Bayesian approach - which can give a Credible Interval for the
parameter p, and also (presumably what is wanted ?) a 'predictive'
interval
for the number of positives in the 42 as yet unobserved cases - hence an
interval for the total number of positives. BUGS is suitable software.
http://www.mrc-bsu.cam.ac.uk/bugs/
If you only want a CI for p, then see robert newcombe's website (or any
number of sites that will do it 'while you wait' (at own risk!))
http://archive.uwcm.ac.uk/uwcm/ms/Robert2.html
[log in to unmask]
Hi Allan,
It is a misorientation to be thinking in terms of normal
approximations in this kind of context anyway. This kind
of sampling (I infer from your description that it is
without replacement) is best approached by exact computation
using the hypergeometric distribution.
Notation: sampling n = 108 from a population of size N = 150
containing M = ?? positives, and getting m = 3 positives in
the sample.
First, the maximum likelihood estimate of M is M = 4:
P(m=3 | M=3) = 0.37
P(m=3 | M=4) = 0.42
P(m=3 | M=5) = 0.30 ...
Second, necessarily M >= 3.
Third, a value of M as low as M=3 cannot be rejected:
see P(m=3 | M=3) = 0.37 above. Hence for any reasonable
confidence level the lower confidence limit is M=3.
Fourth, the smallest value of M such that P(m <= 3 | M)
is <= a small probability alpha can be read from:
M = 10: P(m <= 3 | M) = 0.005
M = 9: P(m <= 3 | M) = 0.015
M = 8: P(m <= 3 | M) = 0.039
M = 7: P(m <= 3 | M) = 0.096
so for a 99% upper confidence limit you have M=10,
for 97.5% you have M=9, and for 95% you have M=8.
Hope this helps!
Ted.
[log in to unmask]
You're right to note that this is a case in which the ordinary Wald
interval is inappropriate. I'm attaching an Excel spreadsheet that
calculates CIs for proportions and their differences using good
methods. Its use is self-explanatory.
Hope this helps.
Robert Newcombe.
[log in to unmask]
Hi Alan,
An interesting example of an uncommon problem when you want a
confidence interval for a parameter which can take on discrete values
only.
I would define the ci by the usual trick of inverting a significance
test.
If we denote by k the number of positives in the full population then x
has a hypergeometric distribution:
x ~ hyper(x, m, k, n)
where m is the size of the full population (i.e. your N =150).
The usual trick here is to use two one-sided hypothesis tests with
significance level alpha/2 to give the lower and upper confidence
limits.
Look first at the upper limit.
The question to ask is, for what values of K would we reject a null
hypothesis that k=K compared to an alternative that k<K.
The critical region here consists of small values of x, so we need the
probability of seeing X=3 or some smaller (more extreme) value.
Thus we would exclude any k for which
P(X<x) = cdf('Hyper', x, m, k, n) < alpha/2
For the lower limit, we want to look exclude values of k for which
P(X>x) = 1 - cdf('Hyper', x-1, m, k, n) < alpha/2
Note that as x is discrete we need to use x-1 in this formula.
SAS has a function for the cdf of the hypergeometric distribution:
data cdf;
m = 150;
n = 108;
x = 3;
do k = 3 to 10;
cdf = cdf('hyper', x, m, k, n);
cdfr = 1 - cdf('hyper', x-1, m, k, n);
output;
end;
run;
proc print data=cdf;
run;
which produces output:
Obs m n x k cdf cdfr
1 150 108 3 3 1.00000 0.37032
2 150 108 3 4 0.73549 0.68773
3 150 108 3 5 0.43112 0.86601
4 150 108 3 6 0.21596 0.94797
5 150 108 3 7 0.09643 0.98127
6 150 108 3 8 0.03938 0.99366
7 150 108 3 9 0.01495 0.99796
8 150 108 3 10 0.00533 0.99937
Therefore, for a 95% confidence interval with equal tail area
probabilities we would exclude k>= 9, and values k<3 are impossible
given the data, so a conservative 95% ci for k is (3, 4, 5, 6, 7, 8).
The actual confidence level for this is 1-(0.01495 + 0) = 0.985
You could also argue for the range (3, 4, 5, 6, 7) as this has a
confidence level 1 - (0.03938 + 0) = 0.96 which is less conservative
but does not have symmetric tail probabilities.
However, for this problem, a Bayesian analysis with a uniform prior for
k is very easy to do and interpret.
The following little program does the calculations:
data Bayes;
m = 150;
n = 108;
x = 3;
do k = 0 to m;
likelihood = pdf('hyper', x, m, k, n);
prior = 1/(m+1);
output;
end;
run;
* calculate posterior distribution for k;
proc sql;
create table Bayes1 as
select bayes.*, prior*likelihood/sum(prior*likelihood) as posterior
from Bayes
order by k
;
quit;
* calculate cdf of posterior distribution;
data Bayes2;
set bayes1;
cdf + posterior;
run;
proc print data=Bayes2;
where k <= 10;
run;
and produces this output;
Obs m n x k likelihood prior posterior
cdf
1 150 108 3 0 0.00000 .006622517 0.00000
0.00000
2 150 108 3 1 0.00000 .006622517 0.00000
0.00000
3 150 108 3 2 0.00000 .006622517 0.00000
0.00000
4 150 108 3 3 0.37032 .006622517 0.26732
0.26732
5 150 108 3 4 0.42322 .006622517 0.30550
0.57282
6 150 108 3 5 0.29712 .006622517 0.21448
0.78730
7 150 108 3 6 0.16393 .006622517 0.11833
0.90563
8 150 108 3 7 0.07770 .006622517 0.05609
0.96172
9 150 108 3 8 0.03303 .006622517 0.02385
0.98556
10 150 108 3 9 0.01291 .006622517 0.00932
0.99488
11 150 108 3 10 0.00471 .006622517 0.00340
0.99828
Thus a conservative 95% highest posterior density credible interval for
k is (3, 4, 5, 6, 7) with posterior probability that k is in this range
= 0.96172.
Here, as is often the case, the Bayesian and frequentist analyses leads
to broadly similar conclusions.
Best wishes
Tim Auton
[log in to unmask]
Alan
What a strange thing to do! If you only have 150 why not either use a
smaller sample, or else do the lot?
The problem of exact CIs for a binomial parameter is standard (if not
always used) when the population is not finite. It is in books such as
'Armitage and Berry' (Statistical methods in Medical research- now in
its 4th edition with an extra author) or 'Statistics with Confidence'
(published by the British medical Journal); they actually use different
formulae.
But the fact that you have a finite population complicates things. The
effect is to reduce the variance of estimates by a factor (N-n)/N; in
this case 42/150 so in the absence of some complex combinatorics I would
use this to correct the normal approximation.
However someone may have a more rigorous approach!
Regards
Phil
[log in to unmask]
I assume it was sampling without replacement? You could report 0.02 -
0.30
as a 100% CI (!), but I guess you want something other than that. Check
out
hypergeometric sampling.
Quentin Burrell
[log in to unmask]
Take a look at:
http://members.aol.com/rcknodt/pubpage.htm
Modstat does over 300 statistical tests and routines. It even has a
routine
to help the user select the proper test to use with the data. It will
also
suggest alternative tests in the event that the conditions of the first
test
are not met or are in doubt.
One of the routines calculates the exact asymmetric confidence interal
for
times when the estimated proportions are very close to 0 or 1. You can
also
write a routine in SPSS to do the work. I believe that you might also
use
Minitab vut I'm not sure if the results are accurate using that
software.
A free two-month trial is available. Simply e-mail me your name and
complete
postal mailing address.
The one time registration fee, should you decide to continue using the
software after the trial period is only $22.00 for individuals and
$17.00 for
students.
Dr. Robert C. Knodt
[log in to unmask]
Looks like a place where the bootstrap for finite population samples
would be useful. If available to you have a look at Lunneborg, CE "Data
analysis by resampling."
Clifford E. Lunneborg
[log in to unmask]
This advice may be worth what it cost you (0, pounds or $, whatever) :)
A test is either positice or negative (not positive). Hence binomial
distribution.
Why not go directly to the binomial distribution & calculations.
Sampling & sampling theory I believe does that.
Cheers,
Jay
[log in to unmask]
Once again, many thanks to all who replied,
Alan.
|
