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

Pathogenic bacteria such as Staphylococcus aureus, Streptococcus pneumoniae and Neisseria meningitidis are major causes of communicable disease, but more often than not, they are carried by people as part of their body's natural microbiota with no ill effect, to the point that they wouldn't even know they were carriers without being screened. Comparisons of asymptomatically carried bacteria to those that have caused invasive disease could reveal the underlying reasons why some people suffer illness, while others do not. Invasive bacteria may differ genetically from those carried asymptomatically - for example they may contain virulence genes or genes that confer resistance to antimicrobial medicines. How do we go about testing for these differences, and if we detect a difference, how do we quantify whether it is meaningful?

To test for a difference in the frequency of a characteristic, such as antibiotic resistance, between different groups, first it is necessary to understand how frequency is estimated - in a rigorous statistical sense - in a single population. Here I will use this example to talk about formal approaches to statistical inference in a simple setting.

Formal approach to inference

Suppose you wanted to estimate the frequency of ciprofloxacin resistance in invasive S. aureus and compare it to the frequency in carried S. aureus. So you isolated bacteria from the blood of 100 people with septicaemia and from the noses of 100 people known to be asymptomatic carriers, you grew the bacteria in the presence and absence of the antibiotic, and tabulated the results (Table 1).

	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.

Estimating the frequency of drug resistance in the invasive bacteria seems totally straightforward - you just read it off the table: 37/(63+37) = 0.37 (i.e. 37%). But all this tells you is the frequency of drug resistance in your sample. Why is that a reasonable estimate of the frequency in the wider population of invasive S. aureus? And how good an estimate is it really?

The first step is to realize that implicit to any survey of this sort is an abstraction of the world and the way it works, an abstraction which could be represented as a set of simplifying assumptions. These assumptions motivated you to carry out the study in the first place, and they represent the ways in which you think your findings will apply to the world more generally.

For example, you probably think that

• The frequency of antibiotic resistance in invasive S. aureus is likely to be steady from day to day - although the exact time-scale over which it is stable you do not know, and
• The frequency of antibiotic resistance in invasive S. aureus is similar from place to place - although the geographic area across which the frequency is similar you do not know.

This abstract place and time from which your sample was drawn, and to which it could be used for prediction, is usually called the population.

Further, you probably think that the world behaves in a reasonably predictable way so that

• The susceptibility of S. aureus to drugs in the body is same as in the laboratory,
• The act of sampling S. aureus from the blood or nose does not affect susceptibility, and
• By sampling S. aureus from the blood of one individual and measuring its susceptibility, you are not influencing the susceptibility of the S. aureus that you will sample from other individuals.

This working knowledge about how the world works actually corresponds to a set of implicit assumptions, which can be represented precisely using a formal mathematical model.

One of the aims of science is to explain the complex world we observe around us in simple terms. This approach is formalised by the parameters of a mathematical model, which can be used to explain or predict complex observations.

In the current example, the formal approach is to estimate a single quantity - the population frequency of antibiotic resistance in invasive S. aureus - in order to explain hundreds of individual observations of antibiotic resistance or susceptibility. This is the parameter of the model. Let's call it f, and figure out how to estimate it.

The Bernoulli model

If we dissect the survey behind Table 1, we can consider the process that led to a single observation. Each observation had one of only two values - susceptible or resistant to the antibiotic. In a formal sense, the process of testing each bacterial isolate for antibiotic resistance was an experiment, and these were the possible outcomes of the experiment.

An experiment with two possible outcomes, the results of which were unknown before the test was carried out, can be modelled using a probabilistic model (also called a stochastic model, or a statistical model). The key idea is to treat the outcome as if it were random at the time it was measured. Let's call the random outcome X, and refer to it as a random variable.

The meaning of random is open to interpretation. Obviously if resistance were genetically determined, e.g. by the presence of a specific allele, then the outcome was never really random, it's just that the person doing the experiment did not know what it would be. The idea of randomness can represent a lack of knowledge, as well as something that is "truly" random (if such a thing exists - e.g. is the outcome of a coin toss truly random?)

The probabilistic model in question is called the Bernoulli model and has a single parameter, the probability of one outcome instead of the other. Suppose we use the Bernoulli model to model the outcome of the test for antibiotic resistance. What is the probability of antibiotic resistance?

Assuming the isolate was representative of the wider population of invasive S. aureus, it would be natural to assume that the probability of antibiotic resistance was equal to the population frequency of antibiotic resistance, f.

This model can be represented using formal notation as follows. Let X be the random variable representing the outcome of the test for resistance, and let X = 1 if the isolate is resistant and X = 0 if it is susceptible. We can write down the probability of a specific observation of X under the Bernoulli model as follows

$$ \begin{array}{ccc} Pr(X=x) & = & \left\{ \begin{array}{cc} f & \textrm{if $x=1$} \\ 1-f & \textrm{if $x=0$} \end{array} \right. \end{array} $$

This, the Bernoulli model, is the simplest probabilistic model of a random variable that has a parameter. Try simulating from it in R:

n = 100; # this is the number of simulated experiments
f = 0.5; # f is the probability of resistance, from 0-1
X = rbinom(n,1,f); # n simulated outcomes (1=resistant)
barplot(table(X),xlab="Outcome",ylab="Number of experiments")

Figure 1. Results of Bernoulli simulations in R using n = 100 independently simulated experiments with probability of antibiotic resistance f = 0.5 (left) and f = 0.9 (right).

Obviously a single observation is not much use in parameter estimation. Before we can estimate the parameter, we have to extend the approach to model all the data as random variables. Let's number the invasive isolates 1 to n and index them by i so that \(X_i\) refers to the observed outcome in the ith invasive isolate. If we are prepared to assume the observations are independent, then we can use the laws of probability to write down the joint probability of all the X's as follows:

$$ \begin{array}{ccc}Pr(X_1,X_2,\dots,X_n|f) & = & \prod_{i=1}^n Pr(X_i|f) \\
& = & \prod_{i=1}^n \left\{ \begin{array}{cc} f & \textrm{if $X_i=1$} \\ 1-f & \textrm{if $X_i=0$} \end{array} \right.\end{array} $$

In words, the top line of this equation states that the joint probability of \(X_1\) to \(X_n\) conditional on f equals the product of the marginal probabilities of \(X_1\) to \(X_n\) conditional on f. Notice that I have made the conditioning on f explicit by using the vertical bar notation, although it was implicit in the definition of the Bernoulli model above. This change of notation is to give added emphasis to the role that the parameter plays in the probability formula. When dealing with probabilities, the idea of conditioning on a variable means considering it fixed. In other words, the above equation states the joint probability of the X's, which are considered random, for a fixed value of the parameter f.

We can simplify the formula if we define Y to be the total number of isolates, out of n, that were resistant, i.e. \( Y=\sum_{i=1}^n X_i\). Substituting Y into the equation gives

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

As an aside, note that we call Y a statistic because it is a summary of the observed data \(X_1\) to \(X_n\). Any summary of data is called a statistic, and since the data are considered random variables, the statistic is also considered a random variable, even though it is a deterministic function of the data. The treatment of data and statistics as random variables contrasts to that of parameters, which are traditionally treated as fixed, even though their values are unknown. It can be potentially confusing to think of things that are known (the data) as random and things that are unknown (the parameters) as non-random.

Classical approaches to parameter estimation

The common thread linking various approaches to parameter estimation is to find parameter values that are consistent with the observed data, in the sense that a model with those parameter values would produce data that are like the ones observed.

I will only consider principled methods of parameter estimation that are based on explicit probabilistic models. This is because it is easy to evaluate how good the performance of such methods are with respect to a specific biological model, and there has been extensive theoretical research into the performance of these methods in general. This definition excludes optimization methods that are not based on explicit probabilistic models, such as parsimony and support vector machines, and ad hoc algorithms such as distance-based clustering and tree-building methods.

For the time being, I am only going to talk about classical methods of inference, postponing a discussion of Bayesian inference, which is the other major approach to probabilistic inference. Classical methods, which include the method of moments and maximum likelihood, are also widely referred to as frequentist methods.

Method of moments

To explain this method I first need to explain the difference between a mean and an expected value, two terms that are often used interchangeably. The sample mean is of course familiar, and in the example it equals

$$ \bar{X} = \frac{1}{n} \sum_{i=1}^n X_i = \frac{Y}{N} $$

The sample mean \(\bar{X}\) is, like Y, a statistic and therefore, like the data itself, considered to be a random variable. This implies that if you sampled two sets of 100 invasive isolates and calculated the sample mean \(\bar{X}\) in each, you could very well get different values. But if you continued to sample sets of 100 invasive isolates and you calculated the mean of the sample means, that quantity would begin to converge on a value known as the expected value.

Expected values, which can be thought of as long-run means, the mean of an infinite sample, or as population means, can be calculated from the probability distribution, and they depend only on the parameters. In the method of moments, you equate a statistic such as the sample mean to its expected value, and rearrange the equation to solve for the parameters. Often the statistics used are a class of statistic called moments which include the mean, variance and skewness of a distribution.

In other words, the method of moments finds parameter values that are consistent with the observed data in the specific sense that if a model with those parameter values were used to simulate data many times, then the long-run mean of the simulated data would match the data actually observed.

In the Bernoulli model, the expected value of X is defined as

$$ \begin{array}{ccl} E(X) & = & \sum_{x=0}^1 x Pr(X=x) \\
& = & 0 \times (1-f) + 1 \times f \\
& = & f \end{array} $$

To show that f is the long-run sample mean of \(\bar{X}\), try this R code:

n = 100
f = 0.75
Xbar = mean(rbinom(n,1,f))
for(i in 2:30) Xbar = c(Xbar,mean(rbinom(n,1,f)))
plot(1:30,Xbar,xlab="Sample number",ylab=expression(bar(X)),ylim=c(0,1))
lines(1:30,cumsum(Xbar)/(1:30),col=2,type="o")
abline(h=0.75,lty=2)
legend("bottom",c("Sample mean","Running mean of sample means","Expected value"),col=c(1,2,1),bty="n",pch=c(21,21,NA),lty=c(0,0,2))

Figure 2. Relationship between sample mean and expected value.

Suppose we use the sample mean as our statistic to estimate the parameter f of the Bernoulli distribution, then by the method of moments we have

$$ \bar{X} = E(X) = f $$
which, if we rearrange for f gives \(\hat{f}=\bar{X}\), proving that for the Bernoulli distribution, the method of moments estimate of f is the sample mean. That was a lot of work for a very intuitive result, but it does give a justification of taking the sample mean to estimate the frequency of antibiotic resistance in invasive S. aureus.

However, there are various problems with the method of moments. When there are multiple parameters, it may not be possible to obtain analytical results (i.e. explicit formulas) for the parameter estimates in terms of the statistics, so numerical techniques are required. There is ambiguity because any statistic can be used, and different parameter estimates can be obtained when different statistics are used. There is no objective way to choose the statistics, which becomes harder the more parameters there are. And sometimes there is no combination of parameter values that can yield expected values equal to the observed statistics. The bottom line is that there are better methods available.

Method of maximum likelihood

Maximum likelihood is widely considered to be the best method for statistical inference in a classical setting. The likelihood is a function that defines the probability of observing the data you actually observed given particular values of the parameters. Likelihood is therefore just another way of looking at probability: some people would say that probability is a function of the data, whereas the likelihood is a function of the parameters. (In truth they are the same thing, and a function of both). Anyway, the likelihood function for n independently sampled Bernoulli random variables is just

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

where Y was defined earlier as the total number of resistant invasive isolates. In Table 1, there were 37 resistant and 63 susceptible invasive isolates, i.e. \(Y=37\) and \(n=100\). For this data, we can look at how the likelihood changes as a function of the parameter f using the following R code:

Y = 37; n = 100; # the data
L = function(f) f^Y * (1-f)^(n-Y)
curve(L,0,1,xlab=expression(f),ylab="Likelihood of f")

Figure 3. Likelihood of f, the frequency of antibiotic resistance in invasive S. aureus, given the data in Table 1.

The method of maximum likelihood estimates the parameters by maximizing the likelihood function with respect to those parameters. Mathematically, this is written as

$$ \hat{f} = \textrm{arg max}_f\,L(f) $$

For a single parameter you can plot the likelihood and actually see where the maximum value lies, and while there are algorithms that will try to find the maximum of a function automatically, these algorithms are slow compared to finding an analytic solution (i.e. a formula). When there are multiple parameters, the algorithms for numerically maximizing the likelihood are prone to getting stuck in local maxima. In general, the more you can do analytically, the better.

To find the maximum likelihood estimate (MLE) analytically, we need to use high school maths. At the maximum of a function, the first derivative of will be zero. So by differentiating the likelihood once, it is possible to find one or more stationary points that correspond to the local and global minima and maxima of the likelihood. In practice, it is usually easier to differentiate the logarithm of the likelihood, \(l(f)=log(L(f))\), which works because the maximum log-likelihood occurs at the same value of the parameters as the maximum likelihood (the technical reason for this convenient fact being that the logarithm is a strictly monotonic transformation).

For the Bernoulli distribution,

$$ \begin{array}{ccl} \frac{d\,l(f)}{df} & = & \frac{d}{df} \left\{ Y \log{f} + (n-Y) \log{(1-f)} \right\} \\
& = & \frac{Y}{f} - \frac{n-Y}{1-f} \end{array} $$

Setting this equal to zero and solving for f we get \(\hat{f}=Y/N\), which can be obtained using the Maple code

solve(diff(Y*log(f)+(n-Y)*log(1-f),f),f);

So in the example, the maximum likelihood estimate is equal to the method of moments estimate, both of which justify the intuitive use of the sample frequency to estimate the population frequency of ciprofloxacin resistance in invasive S. aureus, which was 0.37. It follows that the maximum likelihood estimate for the frequency of ciprofloxacin resistance in carried bacteria is 0.11. How do we assess how good these estimates are, in terms of the uncertainty surrounding them, and how do we decide whether 0.11 is meaningfully different from 0.37? These questions will be the focus of the next tutorial.