Introduction to statistical inference III: antibiotic resistance in carriage versus invasive disease

In the previous two posts I discussed an example study designed to compare the frequency of antibiotic resistance in invasive S. aureus and carried S. aureus (Table 1). The broad aim of the study was to determine whether carried and invasive bacteria differed genetically, with a focus on resistance to the antibiotic ciprofloxacin. In this post I will discuss confidence intervals, their relationship to hypothesis testing, and how to construct them in order to quantify the uncertainty in estimating parameters such as the frequencies of ciprofloxacin resistance in invasive and carried isolates.

	Susceptible	Resistant
Carried	89	11
Invasive	63	37

Table 1. Number of carried versus invasive S. aureus that are susceptible versus resistant to the antibiotic ciprofloxacin.

In the previous post I walked through a formal test of whether the frequency of antibiotic resistance differed in a meaningful way between invasive and carried isolates. Using the likelihood ratio test, it turned out that if there was no real difference in the frequency of antibiotic resistance between invasive and carried isolates, then the probability of observing differences in sampling frequency as large as or larger than those seen in Table 1 was \(p=1.1\times 10^{-5}\) - less than 0.05, the usual cut off for ruling out the null hypothesis.

In the first post, I described maximum likelihood estimation of the frequencies of antibiotic resistance in the two groups, which turned out to be simply the sampling frequencies. Here I will discuss how to quantify uncertainty in those estimates using confidence intervals, and show that there is a one-to-one relationship between hypothesis testing and confidence intervals.

Confidence intervals

When dealing with continuous parameters such as the population frequency of antibiotic resistance - as opposed to discrete variables such as a binary outcome - it is important to recognize that any point estimate is almost certainly wrong. That sounds like a strong statement but consider how many possible values a continuous parameter can take. The frequency of antibiotic resistance is constrained to lie between zero and one. So to the first decimal place, there are 11 possible parameter values (0.0, 0.1, ..., 1.0). To the second decimal place, there are 101 possible parameter values (0.00, 0.01, ..., 1.00). Since there is no limit to the number of decimal places a continuous number can have, it follows that there are in fact an infinite number of possible values between zero and one. The chance that your point estimate is correct to the first decimal place might be decent, to the second, reasonable, but by the time you get to the tenth let alone the hundredth, it is starting to look remote. In other words, the point estimate may be close, but it is unlikely to be spot on.

Quantifying the uncertainty in a point estimate is a good way of indicating how close to the truth it is likely to be. This is because under fairly general conditions, maximum likelihood parameters are consistent. So the smaller the uncertainty, the closer to the truth it will tend to be. In classical (aka "frequentist", as opposed to Bayesian) statistics, the standard way of representing uncertainty in a parameter is through a confidence interval. If you have ever sat through any statistics class, you will probably have heard of confidence intervals, and you will probably have been warned never to say something like this: "The probability that the 95% confidence interval contains the true parameter value is 95%." In fact this statement is valid, as long as you attribute the uncertainty to the confidence interval, and not to the parameter. This is because the confidence interval is a random variable in the same sense that the sample frequency and point estimates are random variables: even though there is no doubt about their values after the experiment, you did not know what they were going to be before the experiment, and if you repeated the experiment they could differ. This is in contrast to the parameter, which is always considered fixed (albeit unknown).

How do you construct a \(100 (1-\alpha)\%\) confidence interval so that it contains the truth \(100 (1-\alpha)\%\) of the time? Look at the likelihood curve for the frequency of antibiotic resistance in invasive isolates below (the following R code generated it):

Y = 37; n = 100; # the data
L = function(f) f^Y * (1-f)^(n-Y)
# Plot the likelihood of f
curve(L,0,1,xlab=expression(f),ylab="Likelihood of f",lwd=2)
# Indicate the position of the maximum likelihood estimate
f_MLE = Y/n
abline(v=f_MLE,lty=2)
# For every value of f, indicate whether a hypothesis test taking that value of f as the null hypothesis would be rejected. If it wouldn't, return the likelihood, otherwise return NA
rejectH0 = Vectorize(function(f) ifelse(pchisq(-2*log(L(f)/L(f_MLE)),1,lower.tail=FALSE)<=0.05,L(f),NA))
curve(rejectH0,add=TRUE,col=2,lwd=2,n=1001)
legend("topright",c("MLE","Rejection\nregion"), lwd=c(1,2),lty=c(2,1),col=c(1,2),bty="n")

Figure 1. Likelihood of the frequency of antibiotic resistance in invasive S. aureus, with the maximum likelihood estimate (MLE, vertical dashed line) and rejection region (red line) marked. The rejection region corresponds to those values of f that would be rejected as the true frequency of antibiotic resistance at a p-value threshold of \(\alpha=0.05\).

Remember that for the hypothesis test of \(f=f_0\) versus \(f\neq f_0\), the p-value threshold \(\alpha\), often taken to be 0.05, determines the frequency with which one would falsely reject the null hypothesis when it was true. In Figure 1, the values of f that would be rejected by this hypothesis test are indicated in red. So with probability 0.05 the red region contains the true value of f and with probability 0.95 the black region contains the true value of f.

This tells us that the answer to how you construct a \(100 (1-\alpha)\%\) confidence interval so that it contains the truth \(100 (1-\alpha)\%\) of the time, is to use the set of parameter values for which the null hypothesis would not be rejected in a hypothesis test with a p-value threshold of \(\alpha\) (black region in Figure 1.) 

The confidence interval can be found quickly by identifying the critical values of the hypothesis test, which occur at the boundary of the "acceptance region" (the null hypothesis is never really accepted, it is just not rejected) and the rejection region. When formally tested, the critical values produce p-values of exactly \(\alpha\). The chi-squared approximation to the likelihood ratio test tells us that this occurs when the difference in the log-likelihood between the alternative and null hypotheses is 1.92 \((\approx 2)\).

The following R code finds the critical values for \(\alpha=0.05\):

Y = 37; n = 100; # data
L = function(f) f^Y * (1-f)^(n-Y); # likelihood
l = function(f) log(L(f)); # log-likelihood
f_hat = Y/n; # MLE
qchisq(0.95,1)/2; # Required difference in log-likelihood
# Solve the equation l(f_hat)-l(f)=1.92, i.e.
# l(f_hat)-l(f)-1.92 = 0
# Lower bound of the 95% confidence interval
uniroot(function(f) l(f_hat)-l(f)-1.92,c(0,f_hat))$root# Upper bound of the 95% confidence interval
uniroot(function(f) l(f_hat)-l(f)-1.92,c(f_hat,1))$root

For the invasive isolates, this gives a 95% confidence interval of \((0.28,0.47)\) and for the carriage isolates, \((0.06,0.18)\). Recall that the MLEs were 0.37 and 0.11 respectively. Confidence intervals constructed like this from the likelihood ratio test will always include the MLE, for any \(\alpha\). This is because likelihood-based inference provides a coherent framework for obtaining parameter estimates, calculating confidence intervals and performing hypothesis tests.

Confidence intervals and the standard error

Sometimes confidence intervals are approximated via a quantity known as the standard error. The standard error is the standard deviation of a statistic - a statistic has a standard deviation because it is a random variable and can vary if the experiment is repeated. A statistic is just a function of the data, so all estimators, including MLEs, are statistics, meaning that (1) they can vary between repeated experiments and (2) a standard deviation could be calculated to express how much they vary and this quantity is called the standard error of the estimate. In other words, a standard error is a measure of uncertainty in an estimate, like a confidence interval.

The distribution of a statistic over repeated experiments is known as its sampling distribution and if the true parameter were known, the sampling distribution could be simulated directly. For example:

f = 0.4; # the "truth"
n = 100; # sample size
m = 10000; # number of simulations
X = replicate(m,rbinom(n,1,f)); # simulate data
Y = apply(X,2,sum); # number resistant isolates per dataset
f_hat = Y/n; # for each dataset, MLE of f
hist(f_hat,col=7,xlab=expression(hat(f)),prob=TRUE,main="Sampling distribution")
sd(f_hat); # standard error of the MLE

Figure 2. Sampling distribution of the maximum likelihood estimate of the frequency of antibiotic resistance.

The sampling distribution, if it could be known, could be used to calculate a confidence interval. For example, if one took the 2.5 and 97.5 percentiles of the sampling distribution of

Calculating confidence intervals by simulation

Often likelihood formulas are much more complex than the one presented here for a sample of independent Bernoulli-distributed random variables, to the extent that they cannot be written down as a simple formula, and they must be calculated numerically on a computer. When the likelihood can only be evaluated numerically, maximizing it involves an optimization algorithm, and finding the confidence interval requires a second step using a so-called root-finding algorithm. In some settings, optimization is easier than root-finding, so instead of finding the confidence intervals using an algorithm like uniroot() in the R code above, they are instead approximated by simulation.

The parametric bootstrap is a popular method of approximating confidence interval

In these situations, it might be time-consuming to evaluate the likelihood, even for a single combination of parameter values.  

To save time, the second step can be avoided by approximating the confidence interval using simulations.