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

In the previous post 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 cover how to formally compare the frequencies to test for a difference between 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.


As a preliminary to testing for differences between invasive and carried bacteria, in the last post I talked about casting the problem in the form of a probabilistic model with explicit parameters, the frequencies in the two groups. I touched on the relationship between the concrete sample and the more abstract concept of the population which it represents, and how the model parameters could be formally estimated.

Maximum likelihood (ML) is the de facto standard method for parameter estimation in "classical" (i.e. non-Bayesian) statistics. One important result is that under fairly general conditions, ML estimators are consistent. This means that as you accumulate more and more independent observations, the estimator is guaranteed to converge to the true parameter value. However, this result might not apply when the true value is near the edge of a parameter range, it does not guarantee good results for small samples, and it is not true that ML estimates are necessarily unbiased.

Hypothesis testing

Testing theories or ideas is an important part of the scientific method, and the textbook approach, taught to children in high school, in which a proposed hypothesis is compared to a simpler "null" hypothesis is intimately related to the theory of statistical hypothesis testing developed by figures including Fisher, Neyman and Pearson.

In the previous post, a Bernoulli model was constructed to describe the process of assaying antibiotic resistance in invasive isolates. The outcome of each assay was encoded as \(X=0\) for susceptible or as \(X=1\) for resistant, and the n isolates were indexed by i so that \(X_i\) refered to the observed outcome in the ith invasive isolate. Using the laws of probability, the joint probability of the n observations was written as

$$ \begin{array}{ccc} Pr(X_1,X_2,\dots,X_n|f) & = & f^Y (1-f)^{n-Y} \end{array} $$

where Y was defined to be the total number of isolates, out of n, that were resistant, i.e. \(Y=\sum_{i=1}^n X_i\), and f was the parameter representing the true frequency of resistance in the population of invasive isolates.

The likelihood of the the parameter, f, given the data is equal to the probability of the data given the parameter. Using differentiation, it was possible to find the value of f that maximized the log-likelihood (equivalent to maximizing the likelihood), which gave a maximum likelihood estimate, or MLE, of \(Y/N\). Applied separately to the invasive and carried isolates, the estimated frequency of antibiotic resistance was 0.37 and 0.11 respectively.

The biological question now is whether 0.37 and 0.11 are meaningfully different? To address this question rigorously, we have to formulate a hypothesis test. That means coming up with a null model and, possibly, an alternative model.

Significance testing

Significance testing is ubiquitous in classical statistics, and it forms the basis of the hypothesis test (although it is more general) so I will present it first. In significance testing, advocated by Fisher, one wishes to stress-test an explicit probabilistic model to see if it is consistent with empirical observations. Only a single model - which you might call the null model or reference model - is made explicit. Significance testing does not need to involve parameter estimation. For example, suppose I wanted to test my hypothesis that the frequency of ciprofloxacin resistance in invasive S. aureus was 0.5. That might seem arbitrary but suppose I'd read in the literature that 50% of Canadian S. aureus were ciprofloxacin resistant, and I wanted to test it over here.

The central tenet of significance testing is to define a test statistic or discrepancy statistic that is sensitive to departures from the null model, and calculate the probability of observing a value of the discrepancy statistic as or more extreme than what was actually observed. This is known as the p-value. In a sense, the p-value measures how "surprising" the observations were according to the null model.

To make things more concrete, let's imagine taking Y, the number of resistant isolates, as the discrepancy statistic. The observed value of Y was 37 for the invasive isolates. This is less than the expected value, \(nf=50\), so we define "as or more extreme" to be "equal to or less than" \(Y_{obs}=37\).

This p-value can be calculated by simulation, using the following R code:

m = 10000; # number of independent simulations
n = 100; # number of independent observations
f = 0.5; # frequency of resistant isolates under the null hypothesis
Y_obs = 37; # observed value of the discrepancy statistic
Y = sum(rbinom(n,1,f))
for(i in 2:m) Y = c(Y,sum(rbinom(n,1,f)))
barplot(table(factor(Y,levels=0:n))/m,xlab="Number of resistant isolates",ylab="Probability",col="orange3",main="Null distribution")
abline(v=Y_obs,lty=2,col="skyblue3")
(p = sum(Y<=Y_obs)/m); # "raw" p-value
(p = (1+sum(Y<=Y_obs))/(1+m)); # "conventional" p-value

Figure 1. Null distribution of the discrepancy statistic, the number of resistant isolates (orange barplot) with the observed value indicated (blue vertical dashed line).




In 10,000 simulations, I simulated discrepancy statistics as or more extreme as that observed 51 times, giving a "raw" p-value of \(51/10000=0.0051\). However, it is conventional when using simulation to treat the observed value as if one of the simulated values, which ensures the estimated p-value can never be exactly zero (which if interpreted literally would imply the observed discrepancy statistic was impossible, rather than exceedingly unlikely). This correction gives a p-value of 0.0052.

Currently, the p-value is one-tailed, meaning that it defines extreme values of the discrepancy statistic in one direction only - smaller than the expected value. Before collecting the data, we would not have known whether resistance was more or less frequent than 0.5, and so we would have been prepared to test values of the discrepancy statistic larger than the expected value as well as smaller. To obtain a two-tailed p-value, we need to define what an equivalently large deviation in the other direction would be. If we define it as a deviation that had an equal one-tailed p-value (but in the other direction), then to get the two-tailed p-value we simply double the one-tailed p-value. The test we have constructed is actually the binomial test, and is implemented in R as binom.test.

So the p-value is the probability, under the null hypothesis, of observing a discrepancy statistic as or more extreme as that observed, where extremity is defined in one or both directions. How small must a p-value be to consider the observed data inconsistent with the model? Typically a threshold of \(\alpha=0.05\) is used. This means that if the null hypothesis were in fact true, we would be wrong one in twenty times.

In our example, the null hypothesis is rejected at the \(\alpha=0.05\) level in the two-tailed test \((p=0.01)\). If we applied to test to the carriage isolates, the p-value would be even smaller.

Can we use significance testing to test whether the carriage and invasive isolates differ in the frequency of antibiotic resistance? The answer is yes, and I will show how, but it is not the best way, as I will explain afterwards.

There is some ambiguity about how to set up the significance test. For example, I could test whether the frequency of resistance in carried isolates is consistent with that estimated for invasive isolates (0.37), or vice versa (0.11). Starting with the former, I can use the above R code to estimate a two-tailed p-value of 0.0002, overwhelmingly rejecting the null hypothesis. For the latter, I have to be careful to reverse the definition of extreme in the one-tailed test, before converting it to obtain a two-tailed p-value of 0.0002. So in this example, I get the same p-value, given the accuracy of the simulation approach. But in general, it is possible to get different results depending which way round you do the test.

Likelihood ratio tests

There are various criticisms of significance testing, some of the fiercest involving disagreements between Bayesian and classical statisticians. But one thing both sides agree on is that explicitly defining the alternative hypothesis leads to more powerful statistical tests, meaning that the risk of incorrectly favouring the null hypothesis is reduced.

Now the null and alternative hypotheses have to model all the observations. We can accommodate this by making the notation more complicated. Let \(X_i\) represent the outcome of the antibiotic resistance assay of an invasive isolate, where i indexes the isolates from 1 to n, and let \(W_i\) represent the outcome of the antibiotic resistance assay of a carried isolate, where i indexes the isolates from 1 to m, which may be different to n. Let \(Y=\sum_{i=1}^n X_i\) be the number of resistant invasive isolates and \(Z=\sum_{i=1}^m W_i\) be the number of resistant carriage isolates. Let \(f_X\) and \(f_W\) be the population frequency of antibiotic resistance in invasive and carriage isolates respectively. Now we can write a likelihood for the full data as:

$$ \begin{array}{rcl} Pr(X_1,\dots,X_n,W_1,\dots,W_m|f_X,f_W) & = & Pr(X_1,\dots,X_n|f_X) Pr(W_1,\dots,W_m|f_W) \\&=& f_X^Y (1-f_X)^{n-Y} f_W^Z (1-f_W)^{m-Z} \end{array} $$

In the null hypothesis, \(f_X=f_W\) and in the alternative hypothesis, \(f_X\neq f_W\). Maximum likelihood estimates of the parameters are used, which gives \(\hat{f}=(11+37)/(100+100)=0.24\) under the null hypothesis and \(\hat{f}_X=0.37, \hat{f}_W=0.11\) under the alternative hypothesis.

The likelihood ratio test is a significance test advocated by Neyman and Pearson, where the discrepancy statistic is defined to be the ratio of likelihoods:

$$ R = \frac{L(\hat{f}|H_0)}{L(\hat{f}_X,\hat{f}_W|H_1)} = \frac{\hat{f}^Y \left(1-\hat{f}\right)^{n-Y} \hat{f}^Z \left(1-\hat{f}\right)^{m-Z}}{\hat{f}_X^Y \left(1-\hat{f}_X\right)^{n-Y} \hat{f}_W^Z \left(1-\hat{f}_W\right)^{m-Z}} $$

This looks complicated but we can calculate it for our data using R:

fhat = function(n,Y) Y/n; # MLE given number of isolates in total and resistant
L = function(n,Y) fhat(n,Y)^Y*(1-fhat(n,Y))^(n-Y); # maximum likelihood given number of isolates in total and resistant
n = 100; # observed number of invasive isolates
Y = 37; # observed number of resistant invasive isolates
m = 100; # observed number of carried isolates
Z = 11; # observed number of resistant carried isolates
(R = L(n+m,Y+Z)/(L(n,Y)*L(m,Z))); # likelihood ratio

This gives \(R=6.3\times 10^{-5}\), which is fine, but what is the null distribution of the test statistic? Under similar conditions used to justify maximum likelihood, the null distribution of \(-2\log{(R)}\) can be approximated by the chi-squared distribution with degrees of freedom equal to the difference in the number of parameters between the models. I won't explain this distribution now, except to say first that it only works when the two models are nested, meaning that the null hypothesis is a special case of the alternative, and that when the conditions are met, this test is the uniformly most powerful test, meaning that there is no other choice of test statistic that can beat it.

In our case the null hypothesis has one parameter (\(f\)) and the alternative has two (\(f_X\) and \(f_W\)). R can be used to calculate the p-value, which is one-tailed because the likelihood of the null hypothesis will always be smaller than the alternative because it is more constrained, so \(R\) will always be less than one, with smaller values indicating a larger departure from the null hypothesis.

pchisq(-2*log(R),df=1,lower.tail=FALSE)

This gives a p-value of \(1.1\times 10^{-5}\), smaller than before, although that is partly because with 10,000 simulations the smallest p-value calculable is \(1\times 10^{-4}\). In fact, this likelihood ratio test is equivalent to a binary logistic regression, which can be coded in R as

resp = rep(c(1,0,1,0),times=c(37,100-37,11,100-11)); # response variable: resistant (1) or susceptible (0)
expl = rep(1:0,each=100); # explanatory variable: invasive (1) or carried (0)
table(expl,resp); # compare to Table 1
pchisq(glm(resp~1,family=binomial)$deviance-glm(resp~expl,family=binomial)$deviance,df=1,lower.tail=FALSE)

Even though the assumptions underlying the chi-squared approximation may not always hold, particularly when parameter estimates are near the edge of the parameter range, or when sample sizes are small, the null distribution for the likelihood ratio test statistic can always be simulated directly under the null hypothesis to obtain a valid test whose discrepancy statistic is well motivated.