22.2.1 Capturing the population parameter

A plausible range of values for the population parameter is called a confidence interval. Using only a point estimate is like fishing in a murky lake with a spear, and using a confidence interval is like fishing with a net. We can throw a spear where we saw a fish, but we will probably miss. On the other hand, if we toss a net in that area, we have a good chance of catching the fish.

If we report a point estimate, we probably will not hit the exact population parameter. On the other hand, if we report a range of plausible values – a confidence interval – we have a good shot at capturing the parameter.

Exercise: If we want to be very certain we capture the population parameter, should we use a wider interval or a smaller interval?⁸⁶

22.2.2 Constructing a confidence interval

A point estimate is our best guess for the value of the parameter, so it makes sense to build the confidence interval around that value. The standard error, which is a measure of the uncertainty associated with the point estimate, provides a guide for how large we should make the confidence interval.

Generally, what you should know about building confidence intervals is laid out in the following steps:

1. Identify the parameter you would like to estimate (for example, \(\mu\)).

2. Identify a good estimate for that parameter (sample mean, \(\bar{X}\)).

3. Determine the distribution of your estimate or a function of your estimate.

4. Use this distribution to obtain a range of feasible values (confidence interval) for the parameter. (For example if \(\mu\) is the parameter of interest and we are using the CLT, then \(\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\sim \textsf{Norm}(0,1)\). We can solve the equation for \(\mu\) to find a reasonable range of feasible values.)

Let’s do an example to solidify these ideas.

Constructing a 95% confidence interval for the mean
When the sampling distribution of a point estimate can reasonably be modeled as normal, the point estimate we observe will be within 1.96 standard errors of the true value of interest about 95% of the time. Thus, a 95% confidence interval for such a point estimate can be constructed:

\[ \hat{\theta} \pm\ 1.96 \times SE_{\hat{\theta}}\] Where \(\hat{\theta}\) is our estimate of the parameter and \(SE_{\hat{\theta}}\) is the standard error of that estimate.

We can be 95% confident this interval captures the true value. The 1.96 can be found using the qnorm() function. If we want .95 in the middle, that leaves 0.025 in each tail. Thus we use .975 in the qnorm() function.

qnorm(.975)

## [1] 1.959964

Exercise:
Compute the area between -1.96 and 1.96 for a normal distribution with mean 0 and standard deviation 1.

pnorm(1.96)-pnorm(-1.96)

## [1] 0.9500042

In mathematical terms, the derivation of this confidence is as follows:

Let \(X_1,X_2,...,X_n\) be an iid sequence of random variables, each with mean \(\mu\) and standard deviation \(\sigma\). The central limit theorem tells us that \[ \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\overset{approx}{\sim}\textsf{Norm}(0,1) \]

If the significance level is \(0\leq \alpha \leq 1\), then the confidence level is \(1-\alpha\). Yes \(\alpha\) is the same as the significance level in hypothesis testing. Thus \[ \mbox{P}\left(-z_{\alpha/2}\leq {\bar{X}-\mu\over \sigma/\sqrt{n}} \leq z_{\alpha/2}\right)=1-\alpha \]

where \(z_{\alpha/2}\) is such that \(\mbox{P}(Z\geq z_{\alpha/2})=\alpha/2\), where \(Z\sim \textsf{Norm}(0,1)\), see Figure 22.1.

Figure 22.1: The pdf of a standard normal distribution showing idea of how to develop a confidence interval.

So, we know that \((1-\alpha)*100\%\) of the time, \({\bar{X}-\mu\over \sigma/\sqrt{n}}\) will be between \(-z_{\alpha/2}\) and \(z_{\alpha/2}\).

By rearranging the expression above and solving for \(\mu\), we get: \[ \mbox{P}\left(\bar{X}-z_{\alpha/2}{\sigma\over\sqrt{n}}\leq \mu \leq \bar{X}+z_{\alpha/2}{\sigma\over\sqrt{n}}\right)=1-\alpha \]

Be careful with the interpretation of this expression. As a reminder \(\bar{X}\) is the random variable here. The population mean, \(\mu\), is NOT a variable. It is an unknown parameter. Thus, the above expression is NOT a probabilistic statement about \(\mu\), but rather about the random variable \(\bar{X}\).

Nonetheless, the above expression gives us a nice interval for “reasonable” values of \(\mu\) given a particular sample.

A \((1-\alpha)*100\%\) confidence interval for the mean is given by: \[ \mu\in\left(\bar{x}\pm z_{\alpha/2}{\sigma\over\sqrt{n}}\right) \]

Notice in this equation we are using the lower case \(\bar{x}\), the sample mean, and thus nothing is random in the interval. Thus we will not use probabilistic statements about confidence intervals when we calculate numerical values from data for the upper and/or lower limits.

In most applications, the most common value of \(\alpha\) is 0.05. In that case, to construct a 95% confidence interval, we would need to find \(z_{0.025}\) which can be found quickly with qnorm():

qnorm(1-0.05/2)

## [1] 1.959964

qnorm(.975)

## [1] 1.959964

22.2.3 Body Temperature Example

Example:
Find a 95% confidence interval for the body temperature data from last lesson.

We need the mean, standard deviation, and sample size from this data. The following R code calculates the confidence interval, make sure you can follow the code.

temperature %>%
  favstats(~temperature,data=.) %>%
  select(mean,sd,n) %>%
  summarise(lower_bound=mean-qt(0.975,129)*sd/sqrt(n),
            upper_bound=mean+qt(0.975,129)*sd/sqrt(n))

##   lower_bound upper_bound
## 1      98.122    98.37646

The 95% confidence interval for \(\mu\) is \((98.12,98.38)\). We are 95% confident that \(\mu\), the average human body temperature, is in this interval. Alternatively and equally relevant, we could say that 95% of similarly constructed intervals will contain the true mean, \(\mu\). It is important to understand the use of the word confident and not the word probability.

There is a link between hypothesis testing and confidence intervals. Remember when we used this data in a hypothesis test, the null hypothesis was \(H_0\): The average body temperature is 98.6 \(\mu = 98.6\). This null hypothesized value is not in the interval, so we could reject the null hypothesis with this confidence interval.

We could also use R to find the confidence interval and conduct the hypothesis test. Read about the function t_test() in the help menu to determine why we used the mu option.

t_test(~temperature,data=temperature,mu=98.6)

## 
##  One Sample t-test
## 
## data:  temperature
## t = -5.4548, df = 129, p-value = 2.411e-07
## alternative hypothesis: true mean is not equal to 98.6
## 95 percent confidence interval:
##  98.12200 98.37646
## sample estimates:
## mean of x 
##  98.24923

Or if you just want the interval:

confint(t_test(~temperature,data=temperature,mu=98.6))

##   mean of x  lower    upper level
## 1  98.24923 98.122 98.37646  0.95

In reviewing the hypothesis test for a single mean, you can see how this confidence interval was formed by inverting the test statistic. As a reminder, the following equation inverts the test statistic.

\[ \mbox{P}\left(\bar{X}-z_{\alpha/2}{\sigma\over\sqrt{n}}\leq \mu \leq \bar{X}+z_{\alpha/2}{\sigma\over\sqrt{n}}\right)=1-\alpha \]

22.2.4 One-sided Intervals

If you remember the hypothesis test for temperature in the central limit theorem lesson, you may be crying foul. That was a one-sided hypothesis test and we just conducted a two-sided test. So far, we have discussed only “two-sided” intervals. These intervals have an upper and lower bound. Typically, \(\alpha\) is apportioned equally between the two tails. (Thus, we look for \(z_{\alpha/2}\).)

In “one-sided” intervals, we only bound the interval on one side. We construct one-sided intervals when we are concerned with whether a parameter exceeds or stays below some threshold. Building a one-sided interval is similar to building two-sided intervals, except rather than dividing \(\alpha\) into two, you simply apportion all of \(\alpha\) to the relevant side. The difficult part is to determine if we need an upper bound or lower bound.

For the body temperature study, the alternative hypothesis was that the mean was less than 98.6. In our confidence interval, we want to find the largest value the mean could be and thus we want the upper bound. We are trying to reject the hypothesis by showing an alternative that is smaller than the null hypothesized value. Finding the lower limit does not help us since the confidence interval indicates an interval that starts at the lower value and is unbounded above. Let’s just make up some numbers; suppose the lower confidence bound is 97.5. All we know is the true average temperature is this value or greater. This is not helpful. However, if we find an upper confidence bound and the value is 98.1, we know the true average temperature is most likely no larger than this value. This is much more helpful.

Repeating the analysis with this in mind.

temperature %>%
  favstats(~temperature,data=.) %>%
  select(mean,sd,n) %>%
  summarise(upper_bound=mean+qt(0.95,129)*sd/sqrt(n))

##   upper_bound
## 1    98.35577

confint(t_test(~temperature,data=temperature,alternative="less"))

##   mean of x lower    upper level
## 1  98.24923  -Inf 98.35577  0.95

Notice the upper bound in the one-sided interval is smaller than the upper bound in the two-sided interval since all 0.05 is going into the upper tail.

22.3.1 Permutation test for two proportions

We start with the null hypothesis which is two-sided since we don’t know if the treatment is harmful or beneficial.

\(H_0\): The treatment and outcome are independent. \(p_{control} - p_{trmt} = 0\) or \(p_{control} = p_{trmt}\).
\(H_A\): The treatment and outcome are dependent \(p_{control} \neq p_{trmt}\).

We will use \(\alpha = 0.05\).

The test statistic is the difference in proportions of patients with stroke in the control and treatment groups.

obs<-diffprop(outcome30~group,data=stent)
obs

##    diffprop 
## -0.09005271

To calculate the p-value, we will shuffle the treatment and control labels because under the null hypothesis, there is no difference.

set.seed(2027)
results <- do(10000)*diffprop(outcome30~shuffle(group),data=stent)

Figure 22.2 a visual summary of the distribution of the test statistics generated under the null hypothesis, the sampling distribution.

results %>%
  gf_dhistogram(~diffprop,fill="cyan",color="black") %>%
  gf_vline(xintercept =obs ) %>%
  gf_theme(theme_classic()) %>%
  gf_labs(title="Sampling distribution of randomization test",
          x="Difference in proportions",y="")

Figure 22.2: Sampling distribution of the difference in proportions.

We next calculate the p-value. We will calculate it as if it were a one-sided test and then double the result to account for the fact that we would reject with a similar value in the opposite tail. Note that the prop1() includes the observed value in the calculation of the p-value.

2*prop1(~(diffprop<=obs),data=results)

##  prop_TRUE 
## 0.00259974

Based on the data, if there were no difference between the treatment and control groups, the probability of the observed differences in proportion of strokes being - 0.09 or more extreme is 0.0026. This is too unlikely, so we reject that there is no difference between control and stroke groups.

22.3.2 Hypothesis test for two proportions using normal model

We must check two conditions before applying the normal model to a generic test of \(\hat{p}_1 - \hat{p}_2\). First, the sampling distribution for each sample proportion must be nearly normal, and secondly, the samples must be independent. Under these two conditions, the sampling distribution of \(\hat{p}_1 - \hat{p}_2\) may be well approximated using the normal model.

The hypotheses are the same as above.

22.3.2.2 Standard error

For our research question the conditions must be verified. Because each group is a simple random sample from less than 10% of the population, the observations are independent, both within the samples and between the samples. The success-failure condition also holds for each sample, at least 10 in each cell is the easiest way to think about it. Because all conditions are met, the normal model can be used for the point estimate of the difference in proportion of strokes

\[p_{control} - p_{trmt} = 0.05726872 - 0.14732143 = - 0.090\] The standard error of the difference in sample proportions is \[ SE_{\hat{p}_1 - \hat{p}_2} = \sqrt{SE_{\hat{p}_1}^2 + SE_{\hat{p}_2}^2}\] \[ = \sqrt{\frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2}}\] where \(p_1\) and \(p_2\) represent the population proportions, and \(n_1\) and \(n_2\) represent the sample sizes.

The calculation of the standard error for our problem must be done carefully. Remember in hypothesis testing, we assume the null hypothesis is true; this means the proportions of strokes must be the same.

\[SE = \sqrt{\frac{p(1-p)}{n_{control}} + \frac{p(1-p)}{n_{trmt}}}\] We don’t know the exposure rate, \(p\), but we can obtain a good estimate of it by pooling the results of both samples: \[\hat{p} = \frac{\text{# of successes}}{\text{# of cases}} = \frac{13 + 33}{451} = 0.102\] This is called the pooled estimate of the sample proportion, and we use it to compute the standard error when the null hypothesis is that \(p_{control} = p_{trmt}\).

\[SE \approx \sqrt{\frac{\hat{p}(1-\hat{p})}{n_{control}} + \frac{\hat{p}(1-\hat{p})}{n_{trmt}}}\]

\[SE \approx \sqrt{\frac{0.102(1-0.102)}{227} + \frac{0.102(1-0.102)}{224}} = 0.0285\]

The test statistic is \[Z = \frac{\text{point estimate} - \text{null value}}{SE} = \frac{-.09 - 0}{0.0285} = - 3.16 \]

The p-value is

2*pnorm(-3.16)

## [1] 0.001577691

Which is close to what we got with permutation test. This should not surprise us as the sampling distribution under the permutation test looked normal.

Figure 22.3 plots the empirical sampling distribution from the permutation test again with a normal density curve overlayed.

results %>%
  gf_dhistogram(~diffprop,fill="cyan",color="black") %>%
  gf_vline(xintercept =obs ) %>%
  gf_dist("norm",sd=0.0285,color="red") %>%
  gf_theme(theme_classic()) %>%
  gf_labs(title="Sampling distribution of randomization test",
          subtitle="Reference normal distribution in red",
          x="Difference in proportions")

Figure 22.3: The sampling distribution of the randomization test with a normal distribution plotted in red.

22.3.3 Confidence interval for two proportions using normal model

The conditions for applying the normal model have already been verified, so we can proceed to the construction of the confidence interval. Remember the form of the confidence interval is

\[\text{point estimate} \ \pm\ z^{\star}SE\]

Our point estimate is -0.09. The standard error is different since we can’t assume the proportion of strokes are equal. We will estimate the standard error from

\[SE = \sqrt{\frac{p_{control}(1-p_{control})}{n_{control}} + \frac{p_{trmt}(1-p_{trmt})}{n_{trmt}}}\]

\[SE \approx \sqrt{\frac{0.057(1-0.057)}{227} + \frac{0.15(1-0.15)}{224}} = 0.0284\]

It is close to the pooled value because of the nearly equal sample sizes.

The critical value is found from the normal quantile.

qnorm(.975)

## [1] 1.959964

The 95% confidence interval is

\[ - 0.09 \ \pm\ 1.96 \times 0.0284 \quad \to \quad (-0.146,- 0.034)\] We are 95% confident that the difference in proportions of strokes in the control and treatment groups is between -0.146 and -0.034. Since this does not include zero, we are confident they are different. This supports the hypothesis tests. This confidence interval is not an accurate method for smaller samples sizes. This is because the actual coverage rate, the percentage of intervals that contain the true population parameter, will not be the nomial coverage rate. This means it is not true that 95% of similarly constructed 95% confidence intervals will contain the true parameter. This because the pooled estimate of the standard error is not accurate for small sample sizes. For the example above, the sample sizes are large and the performance of the method should be adequate.

Of course, R has a built in function to calculate the hypothesis test and confidence interval for two proportions.

prop_test(outcome30~group,data=stent)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  tally(outcome30 ~ group)
## X-squared = 9.0233, df = 1, p-value = 0.002666
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.03022922 0.14987619
## sample estimates:
##    prop 1    prop 2 
## 0.9427313 0.8526786

The p-value is a little different from the one we calculated and closer to the randomization test, which is an approximation of the exact permutation test, because a correction factor was applied. Read online about this correction to learn more. We run the code below with the correction factor off and get the same p-value as we calculated above. The confidence interval is a little different because the function used no stroke as its success event, but since zero is not in the interval, we get the same conclusion.

prop_test(outcome30~group,data=stent,correct=FALSE)

## 
##  2-sample test for equality of proportions without continuity
##  correction
## 
## data:  tally(outcome30 ~ group)
## X-squared = 9.9823, df = 1, p-value = 0.001581
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.03466401 0.14544140
## sample estimates:
##    prop 1    prop 2 
## 0.9427313 0.8526786

Essentially, confidence intervals and hypothesis tests serve similar purposes, but answer slightly different questions. A confidence interval gives you a range of feasible values of a parameter given a particular sample. A hypothesis test tells you whether a specific value is feasible given a sample. Sometimes you can informally conduct a hypothesis test simply by building an interval and observing whether the hypothesized value is contained in the interval. The disadvantage to this approach is that it does not yield a specific \(p\)-value. The disadvantage of the hypothesis test is that it does not give a range of values for the test statistic.

As with hypothesis tests, confidence intervals are imperfect. About 1-in-20 properly constructed 95% confidence intervals will fail to capture the parameter of interest. This is a similar idea to our Type 1 error.