Sample Size
Among the more important issues in Statistics are questions concerning the sample size required to achieve an answer of a certain quality. In this section we will study these issues.
Sample Size for Confidence Intervals
When we find an interval estimate of a parameter the width of the interval is an indicator of how well we are estimating the parameter. Specifically one consider 1/2 of the width, which is then called the error E.
Example (7.1.1)
Normal Mean
A \((1-\alpha)100\%\) confidence interval for the normal mean with unknown standard deviation is given by \(\bar{x}\pm t_{\alpha/2, n-1}s/\sqrt n\), so we have \(E=t_{\alpha/2, n-1}s/\sqrt n\).
Notice that E here is a random variable as it depends on the sample standard deviation s. It might therefore be better to consider \(E[E]=E[t_{\alpha/2, n-1}s/\sqrt n]=t_{\alpha/2, n-1}\sigma/\sqrt n\)
Example (7.1.2)
Binomial p
We previously found the the Clopper-Pearson intervals for a Binomial p to be
\[ \begin{aligned} &L(x) = qbeta(\alpha/2, x, n-x+1)\\ &U(x) = qbeta(1-\alpha/2, x, n-x+1) \end{aligned} \] so
\[E=\left(qbeta(1-\alpha/2, x, n-x+1)-qbeta(\alpha/2, x, n-x+1)\right)/2\]
Notice that one of the items in the calculation of the estimation error is the sample size n. So one can turn this around ask the question: what sample size is needed to achieve a specified error E?
Example (7.1.3)
Normal Mean
We have \(E=t_{\alpha/2, n-1}s/\sqrt n\), so we have \(n=\frac{t_{\alpha/2, n-1}^2s^2}{E^2}\). However, we can not find the right side because it includes n twice:
in \(t_{\alpha/2, n-1}\). Here one generally assumes that n will be sufficiently large so that \(t_{\alpha/2, n-1}\approx z_{\alpha/2}\), which means \(n>30\) or so.
in \(s^2\). Here we will need some idea of the population standard deviation. Generally one needs to have some information from similar experiments, or one has to do a pilot study.
What should E be? This depends on the experiment and what kind of error is acceptable to still have a useful.
**Numerical example*: say we know \(\sigma^2=12.7\), we want \(E=0.5\) and find a 90% confidence interval, the the required sample size is
round(qnorm(1-0.1/2)^2*12.7/0.5^2)## [1] 137Example (7.1.4)
Binomial p
We wish to check the parts in a shipment for faulty ones. Typically about 1 in 20 will be faulty. We want to find a 95% confidence interval with an error no more than 0.03.
E=function(n) {
x=round(n/20)
(qbeta(1-0.05/2, x, n-x+1)-qbeta(0.05/2, x, n-x+1))/2
}
n=20
repeat {
n=n+1
if(n%%10==0) cat(n," ",round(E(n), 4),"\n")
if(E(n)<0.03) break
}## 30 0.082
## 40 0.0627
## 50 0.0508
## 60 0.0524
## 70 0.0522
## 80 0.046
## 90 0.041
## 100 0.0414
## 110 0.0413
## 120 0.038
## 130 0.0352
## 140 0.0353
## 150 0.0352
## 160 0.0331
## 170 0.0312
## 180 0.0313n## [1] 189Sample Size when Testing
In this section we will consider the question of sample size if we plan on doing a hypothesis test. In this context the role of the estimation error is played by the power of the test. In other words we want a sufficiently large sample size so that the power of the test is reasonably large. This often means a power of at least 80%.
Of course the power of a test also depends on how wrong the null hypothesis is, that is the difference between the value of the parameter under the null hypothesis \(\theta_0\) and the true value \(\theta_1\). However, we do not know the true value \(\theta_1\)! In this case the expert in the subject matter needs to decide on an effect size, that is the smallest difference between \(\theta_0\) and \(\theta_1\) that is of practical importance. For example, say we test a new medication for a disease that using current treatments is cured in 10 day. If our new treatment cures the disease in 8 days, that seems quite important, but if it does it in 9.8 days, maybe not. So maybe the doctors tell us that anything less than a 1 day improvement is not enough.
Example (7.1.5)
Normal Mean
We previously found a test based on the test statistic \(T=\sqrt n (\bar{x}-\mu_0)/s\) and a rejection region \(\vert T\vert>t_{\alpha/2, n-1}\). So we know
\[ \begin{aligned} &P\left(\text{reject }H_0 \vert H_0\text{ is true}\right) = \\ &P\left(\vert \sqrt n (\bar{X}-\mu_0)/s\vert>t_{\alpha/2, n-1}\right\vert \mu=\mu_0)=1-\alpha \end{aligned} \]
now
\[ \begin{aligned} &\text{Power} = P\left(\text{reject }H_0 \vert H_0\text{ is false}\right) = \\ &P_{\mu_1}\left(\vert \sqrt n (\bar{X}-\mu_0)/s\vert>t_{\alpha/2, n-1}\right) = \\ &1-P_{\mu_1}\left(\vert \sqrt n (\bar{X}-\mu_0)/s\vert<t_{\alpha/2, n-1}\right) = \\ &1-P_{\mu_1}\left(-t_{\alpha/2, n-1}< \sqrt n (\bar{X}-\mu_0)/s<t_{\alpha/2, n-1}\right) = \\ &1-P_{\mu_1}\left(-t_{\alpha/2, n-1}< \sqrt n (\bar{X}-\mu_1+\mu_1-\mu_0)/s<t_{\alpha/2, n-1}\right) = \\ &1-P_{\mu_1}\left(-t_{\alpha/2, n-1}< \sqrt n (\bar{X}-\mu_1)/s+\sqrt n (\mu_1-\mu_0)/s<t_{\alpha/2, n-1}\right) =\\ &1-P_{\mu_1}\left(-t_{\alpha/2, n-1}-\sqrt n (\mu_1-\mu_0)/s< \sqrt n (\bar{X}-\mu_1)/s<t_{\alpha/2, n-1}-\sqrt n (\mu_1-\mu_0)/s\right) =\\ &1-\left(\Phi(t_{\alpha/2, n-1}-\sqrt n (\mu_1-\mu_0)/s)-\Phi(-t_{\alpha/2, n-1}-\sqrt n (\mu_1-\mu_0)/s)\right) \end{aligned} \]
pwr.normal=function(n, mu0=0,mu1=1,sigma=1,alpha=0.05) {
crit=qt(1-alpha/2, n-1)
tmp=sqrt(n)*(mu1-mu0)/sigma
100*(1-(pnorm(crit-tmp)-pnorm(-crit-tmp)))
}
mu1=seq(-1, 1, length=250)
plot(mu1, pwr.normal(20, mu1=mu1), type="l", lwd=2,col="blue", ylab="Power")
If we now fix mu1 we can also find n:
n=2
repeat {
n=n+1
if(pwr.normal(n, mu1=1)>80) break
}
cat(n, " ", round(pwr.normal(n, mu1=1)),"\n")## 10 82and so here a sample of size 10 is required for a power of 80%.
Example (7.1.6)
Binomial p
We saw before that a test can be based on \(Z=\frac{X-np_0}{\sqrt{np_0(1-p_0)}}\sim N(0,1)\) and we reject the null hypothesis if \(|Z|>z_{\alpha/2}\). Let’s use this to find the sample size required if \(p_0=0.25,p_1=0.3\) and \(\alpha=0.05\).
\[ \begin{aligned} &\text{Power} = P\left(\text{reject }H_0 \vert H_0\text{ is false}\right) = \\ &P\left(\vert \frac{X-np_0}{\sqrt{np_0(1-p_0)}}\vert>z_{\alpha/2}\vert p=p_1\right) = \\ &1-P\left(-z_{\alpha/2}< \frac{X-np_0}{\sqrt{np_0(1-p_0)}}<z_{\alpha/2}\vert p=p_1\right) \end{aligned} \]
Now
\[ \begin{aligned} &P\left(\frac{X-np_0}{\sqrt{np_0(1-p_0)}}<z\right) = \\ &P\left(\frac{X-np_0}{\sqrt{np_1(1-p_1))}}<z\sqrt{\frac{p_1(1-p_1))}{p_0(1-p_0))}}\right) = \\ &P\left(\frac{X-np_1}{\sqrt{np_1(1-p_1)}}<z\sqrt{\frac{p_1(1-p_1))}{p_0(1-p_0)}}+\frac{n(p_0-p_1)}{\sqrt{np_1(1-p_1)}}\right) =\\ &\Phi(z\sqrt{\frac{p_1(1-p_1))}{p_0(1-p_0)}}+\frac{n(p_0-p_1)}{\sqrt{np_1(1-p_1)}}) \end{aligned} \]
and so
pwr.binomial=function(n, p0=0.5,p1=0.6,alpha=0.05) {
tmp1=-qnorm(1-alpha/2)*sqrt(p1*(1-p0)/p0/(1-p1))+
n*(p0-p1)/sqrt(n*p1*(1-p1))
tmp2=qnorm(1-alpha/2)*sqrt(p1*(1-p0)/p0/(1-p1))+
n*(p0-p1)/sqrt(n*p1*(1-p1))
100*(1-(pnorm(tmp2)-pnorm(tmp1)))
}
p1=seq(0.1, 0.99, length=250)
plot(p1, pwr.binomial(20, p1=p1), type="l", lwd=2,col="blue", ylab="Power")
Also, if \(p1=0.6\) the sample size required is
n=10
repeat{
n=n+1
if(pwr.binomial(n, p1=0.6)>80) break
}
cat(n, " ", round(pwr.binomial(n, p1=0.6)),"\n")## 253 80Notice that this problem is not symmetric in \(p_0-p_1\):
n=10
repeat{
n=n+1
if(pwr.binomial(n, p1=0.4)>80) break
}
cat(n, " ", round(pwr.binomial(n, p1=0.4)),"\n")## 144 80