5  Statistical Analysis

Although the result seems obvious, we might think of using a t-test to determine whether or not the difference we see in Figure 4.2 is statistically significant at some threshold level, or a linear model to estimate the size of the difference, with confidence intervals. In either case, it is good practice to check whether the distribution of the data meets the necessary assumptions for the analysis.

# Read experiment data
dfm <- readr::read_tsv(file.path("assets", "data", "farm_scenario_1.tsv"),
                       col_types="ffdd") |>
  dplyr::mutate(ratio_diff = round(CFU_NA_AMP / CFU_NA, 4))

DT::datatable(dfm)

Table 5.1: Simulated experimental data for the BM110 microbiology laboratory experiment. There are 25 samples each from farms A and B, each plated on nutrient agar (NA) and nutrient agar plus Ampicillin (NA/Amp) plates to calculate CFU/mL values in the presence and absence of antibiotic selection. Columns indicate sample number, farm identifier, the CFU/mL obtained from the NA plate (CFU_NA), and the CFU/mL obtained from the NA/Amp plate (CFU_NA_AMP).

5.1 Checking normality

Both t-tests and linear models make similar assumptions about data normality (i.e. that the data fit a normal distribution). For linear models the assumption is that the residuals are distributed normally, and for a pairwise (i.e. t-test) the equivalent assumption is that the observed data are distributed normally around the mean. Since we have a small amount of data, we can use the Shapiro-Wilk normality test to check this, as in the code below.

# Shapiro-Wilk test for farm A NA/Amp:NA ratio data
shapiro.test(dfm |> dplyr::filter(farm == "A") |> dplyr::pull(ratio_diff))

    Shapiro-Wilk normality test

data:  dplyr::pull(dplyr::filter(dfm, farm == "A"), ratio_diff)
W = 0.98593, p-value = 0.9725

# Shapiro-Wilk test for farm B NA/Amp:NA ratio data
shapiro.test(dfm |> dplyr::filter(farm == "B") |> dplyr::pull(ratio_diff))

    Shapiro-Wilk normality test

data:  dplyr::pull(dplyr::filter(dfm, farm == "B"), ratio_diff)
W = 0.97526, p-value = 0.7781

The Shapiro-Wilk test is a hypothesis test that, like a t-test, returns a p-value that we can use to check against a desired threshold (e.g. p < 0.05) or, more robustly, report directly. For the Shapiro-Wilk test, a low (e.g. <0.05) p-value indicates that the data likely do not conform to a normal distribution. Here however, p values are high for both farms and there is no suggestion that the data are not normal. We can proceed to use a t-test or linear model.

5.2 Checking for unequal variance (heteroscedasticity)

An assumption of standard t-tests and linear models is that the groups of data have the same variance - the same degree of dispersion in the data. If this assumption is not met, then variants of the statistical test can usually be employed to accommodate the deviation.

As we have two groups (farms A and B) we can use the F-test to check for equal variance. In the F-test, we are checking to see if the ratio of variances (the value F) is plausibly equal to unity (i.e. 1).

# Test for equal variance in the ratio values from farms A and B
var.test(ratio_diff ~ farm, data=dfm)

    F test to compare two variances

data:  ratio_diff by farm
F = 1.2092, num df = 24, denom df = 24, p-value = 0.6455
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.5328429 2.7439383
sample estimates:
ratio of variances 
          1.209168 

The result here has a 95% confidence interval that includes the value 1, and a p-value of 0.6455 (which is greater than 0.05) so we can conclude that there is no evidence of unequal variance. We can use a standard t-test or linear model.

5.3 A t-test

Our hypothesis was stated as

We hypothesise that, given bacterial communities sampled from farms A (no antibiotic use) and B (prophylactic antimicrobial use), the proportion of antimicrobial-resistant bacteria will be greater in the communities from farm B.

Or, formally:

• \(H_0\) (null hypothesis): There is either (i) no difference between the proportions of antimicrobial-resistant bacteria in samples from farms A (no antibiotic use) and B (antimicrobial use), or (ii) the samples from farm A have a higher proportion of antimicrobial-resistant bacteria.
• \(H_1\) (null hypothesis): The samples from farm B have a higher proportion of antimicrobial-resistant bacteria.

This means that we are explicitly interested in directionality: we want to know if there is evidence to support the proportion of antimicrobial-resistant bacteria from farm B being specifically higher than the proportion from farm A - not just different to the proprtion farm A.

To test this outcome, we must use a one-tailed t-test. This test rejects the null hypothesis when the proportion of antimicrobial-resistant bacteria in samples from farm A is less than that in samples from farm B.

# One-tailed test of ratio values from farms A and B
# This test rejects the null hypothesis when the proportion of
# antimicrobial-resistant bacteria in samples from farm A is less than that
# in samples from farm B
t.test(ratio_diff ~ farm, data=dfm, alternative="less")

    Welch Two Sample t-test

data:  ratio_diff by farm
t = -82.169, df = 47.574, p-value < 2.2e-16
alternative hypothesis: true difference in means between group A and group B is less than 0
95 percent confidence interval:
       -Inf -0.2445435
sample estimates:
mean in group A mean in group B 
        0.04850         0.29814 

The result from this t-test is very strongly in favour of rejecting the null hypothesis and concluding that the samples from farm B have an increased proportion of antimicrobial-resistant bacteria.

5.4 A linear model

We don’t have to use a t-test to understand whether our data support our hypothesis. t-tests are special cases of a more general class of model - linear models - that are much more powerful and widely applicable. We can use a linear model in almost exactly the same way as the t-test above.

A key difference between linear models and t-tests is the amount of information that model provides us with. Rather than simply testing a hypothesis of difference, we can directly estimate effect sizes.

# Linear model of ratio values from farms A and B
model <- lm(ratio_diff ~ farm, data=dfm)  # fit the model
emmeans::emmeans(model, ~farm)   # means and confidence intervals for farms

 farm emmean      SE df lower.CL upper.CL
 A    0.0485 0.00215 48   0.0442   0.0528
 B    0.2981 0.00215 48   0.2938   0.3025

Confidence level used: 0.95 

emmeans::contrast(emmeans::emmeans(model, ~farm),
                  method="pairwise")

 contrast estimate      SE df t.ratio p.value
 A - B       -0.25 0.00304 48 -82.169 <0.0001

confint(emmeans::contrast(emmeans::emmeans(model, ~farm),
                  method="pairwise"))  # estimate difference between farms

 contrast estimate      SE df lower.CL upper.CL
 A - B       -0.25 0.00304 48   -0.256   -0.244

Confidence level used: 0.95 

In this case, the coefficients returned tell us that:

1. the mean ratio for farm A (Intercept) is 0.0485, with a 95% confidence interval [0.0442,0.0528]
2. the mean ratio for farm B is 0.2981, with a 95% confidence interval [0.2938,0.3025]

and the contrast tells us that the estimated difference between ratios for farms A and B is 0.25, with a 95% confidence interval of [-0.256,-0.244].

5.5 Interpreting the statistical results

The t-test gives us the following information:

• mean of farm A: 0.04850
• mean of farm B: 0.29814
• p-value < 2.2e-1 (i.e. we should reject the null hypothesis of the one-tailed test)

The linear model gives us the following information:

• mean of farm A: 0.0485 (95% CI: [0.0442,0.0528])
• mean of farm B: 0.2981 (95% CI: [0.2938,0.3025])
• difference between ratio for farm A and farm B (i.e. “effect size”): 0.25 (95% CI [-0.256,-0.244])
• p-value <0.0001 (i.e. we should reject the null hypothesis that the means are equal)