1 Smelly armpit dataset

Smelly armpits are not caused by sweat, itself. The smell is caused by specific micro-organisms belonging to the group of Corynebacterium spp. that metabolise sweat. Another group of abundant bacteria are the Staphylococcus spp., these bacteria do not metabolise sweat in smelly compounds.

The CMET-group at Ghent University does research to on transplanting the armpit microbiome to save people with smelly armpits.

  • Proposed Therapy:

    1. Remove armpit-microbiome with antibiotics
    2. Influence armpit microbiome with microbial transplant, see this 2 minute talk on youtube
  • Experiment:

    • 20 students with smelly armpits are attributed to one of two treatment groups
    • placebo (only antibiotics)
    • transplant (antibiotica followed by microbial transplant).
    • The microbiome is sampled 6 weeks upon the treatment
    • The relative abundance of Staphylococcus spp. on Corynebacterium spp. + Staphylococcus spp. in the microbiome is measured via DGGE (Denaturing Gradient Gel Electrophoresis).

2 Goal

The overarching goal of this research was to assess if the relative abundance Staphylococcus spp. in the microbiome of the armpit is affected by transplanting the microbiome. To this end the researchers randomized patients to two treatment: A treatment with antibiotics only and a treatment with antibiotics and a microbial transplant.

In the tutorial on hypotheses testing we will use a formal statistical test to generalize the results from the sample to that of the population.

3 Import the dataset

# Load the libraries
library(tidyverse)

Import the data

ap <- read_csv("https://raw.githubusercontent.com/statOmics/PSLSData/main/armpit.csv")
## Rows: 20 Columns: 2
## ── Column specification ──────────────────────────────────────────────
## Delimiter: ","
## chr (1): trt
## dbl (1): rel
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(ap)
## Rows: 20
## Columns: 2
## $ trt <chr> "placebo", "placebo", "placebo", "placebo", "placebo", "…
## $ rel <dbl> 54.99208, 31.84466, 41.09948, 59.52064, 63.57341, 41.486…

4 Data Exploration

A crucial first step in a data analysis is to visualize and to explore the raw data.

ap %>%
  ggplot(aes(x = trt, y = rel, fill = trt)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  ylab("relative abundance (%)") +
  xlab("treatment group") +
  stat_summary(
    fun = mean, geom = "point",
    shape = 5, size = 3, color = "black",
  )

We clearly see that, on average, the subjects who had a microbial transplant have a higher relative abundance of Staphylococcus spp. But is this difference statistically significant so that we can generalized what we observe in the sample to the population?

We can test this with an unpaired, two-sample t-test, which falsifies the null hypothesis that there is on average no difference in relative abundance of Staphylococcus in the armpit microbiome between the transplant and the placebo group against the alternative hypothesis that there is a difference in average abundance of Staphyloccocus in the armpit microbiome between the transplant and placebo treatment.

But, before we can start the analysis, we must check if all assumptions to perform a t-test are met.

5 Analysis

5.1 Check the assumptions

  1. The observations are independent. This has to be guaranteed by the design.

  2. The data (rel) are normally distributed in each of the groups

  3. The variability within both groups is similar.

To check the normality assumption, we will use QQ plots.

ap %>%
  ggplot(aes(sample = rel)) +
  geom_qq() +
  geom_qq_line() +
  facet_grid(cols = vars(trt))

We can see that all of the data lies nicely around the quantile-quantile line (black line). As such, we may assume that our data are normally distributed.

For the third assumption, we must compare the within-group variability of both groups. We can do this visually:

ap %>%
ggplot(aes(x = trt, y = rel)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  ylab("relative abundance (%)") +
  xlab("treatment group") +
  stat_summary(
    fun = mean, geom = "point",
    shape = 5, size = 3, color = "black",
)

Here we can see that the interquartile range is approximately equal for groups.

As all three assumptions are met we may continue with performing the unpaired two-sample t-test.

5.2 Hypothesis test

output <- t.test(
  rel ~ trt,
  data = ap,
  conf.level = 0.95,
  var.equal = TRUE
)
output
## 
##  Two Sample t-test
## 
## data:  rel by trt
## t = -5.0334, df = 18, p-value = 8.638e-05
## alternative hypothesis: true difference in means between group placebo and group transplant is not equal to 0
## 95 percent confidence interval:
##  -31.53191 -12.96072
## sample estimates:
##    mean in group placebo mean in group transplant 
##                 44.15496                 66.40127

5.3 Conclusion

On average the relative abundance of Staphylococcus spp. in the microbiome of the armpit in the transplant group is extremely significantly different from that in the placebo group (\(p<<0.001\)). The relative abundance of Staphylococcus spp. is on average 22.2% larger in the transplant group than in the placebo group (95% CI [13.0,31.5]%).

6 Data Analysis with a linear model

6.1 Fit the model

We do not need to check the assumptions again because we have done this already in the data exploration when analysing the data using the t-test.

mod <- lm(rel ~ trt, data = ap)
mod
## 
## Call:
## lm(formula = rel ~ trt, data = ap)
## 
## Coefficients:
##   (Intercept)  trttransplant  
##         44.15          22.25

6.2 Inference

  1. hypothesis test
summary(mod)
## 
## Call:
## lm(formula = rel ~ trt, data = ap)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.714  -8.779  -1.848   6.981  19.419 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     44.155      3.125  14.128 3.49e-11 ***
## trttransplant   22.246      4.420   5.033 8.64e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.883 on 18 degrees of freedom
## Multiple R-squared:  0.5846, Adjusted R-squared:  0.5616 
## F-statistic: 25.33 on 1 and 18 DF,  p-value: 8.638e-05
  1. confidence interval
confint(mod)
##                  2.5 %   97.5 %
## (Intercept)   37.58905 50.72086
## trttransplant 12.96072 31.53191
  1. How do your results compare to the analysis with the t-test.

  2. Formulate your conclusion based on the output of the linear model.

LS0tCnRpdGxlOiAiRXhlcmNpc2UgNi4zOiBMaW5lYXIgcmVncmVzc2lvbiB2cyB0LXRlc3Qgb24gdGhlIGFybXBpdCBkYXRhc2V0IgphdXRob3I6ICJMaWV2ZW4gQ2xlbWVudCwgSmVyb2VuIEdpbGlzIGFuZCBNaWxhbiBNYWxmYWl0IgpkYXRlOiAic3RhdE9taWNzLCBHaGVudCBVbml2ZXJzaXR5IChodHRwczovL3N0YXRvbWljcy5naXRodWIuaW8pIgotLS0KCiMgU21lbGx5IGFybXBpdCBkYXRhc2V0CgpTbWVsbHkgYXJtcGl0cyBhcmUgbm90IGNhdXNlZCBieSBzd2VhdCwgaXRzZWxmLiBUaGUgc21lbGwgaXMgY2F1c2VkCmJ5IHNwZWNpZmljIG1pY3JvLW9yZ2FuaXNtcyBiZWxvbmdpbmcgdG8gdGhlIGdyb3VwIG9mCipDb3J5bmViYWN0ZXJpdW0gc3BwLiogdGhhdCBtZXRhYm9saXNlIHN3ZWF0LiBBbm90aGVyIGdyb3VwIG9mIGFidW5kYW50IGJhY3RlcmlhCmFyZSB0aGUgKlN0YXBoeWxvY29jY3VzIHNwcC4qLCB0aGVzZSBiYWN0ZXJpYSBkbyBub3QgbWV0YWJvbGlzZSBzd2VhdCBpbiBzbWVsbHkKY29tcG91bmRzLgoKVGhlIENNRVQtZ3JvdXAgYXQgR2hlbnQgVW5pdmVyc2l0eSBkb2VzIHJlc2VhcmNoIHRvIG9uIHRyYW5zcGxhbnRpbmcgdGhlIGFybXBpdAptaWNyb2Jpb21lIHRvIHNhdmUgcGVvcGxlIHdpdGggc21lbGx5IGFybXBpdHMuCgotIFByb3Bvc2VkIFRoZXJhcHk6CiAgCTEuIFJlbW92ZSBhcm1waXQtbWljcm9iaW9tZSB3aXRoIGFudGliaW90aWNzCiAgICAyLiBJbmZsdWVuY2UgYXJtcGl0IG1pY3JvYmlvbWUgd2l0aCBtaWNyb2JpYWwgdHJhbnNwbGFudCwgc2VlIHRoaXMgMiBtaW51dGUKICAgICAgIHRhbGsgb24gW3lvdXR1YmVdKGh0dHBzOi8veW91dHUuYmUvOVJJRnlxTFhkVncpCgotIEV4cGVyaW1lbnQ6CgogICAgLSAyMCBzdHVkZW50cyB3aXRoIHNtZWxseSBhcm1waXRzIGFyZSBhdHRyaWJ1dGVkIHRvIG9uZSBvZgogICAgICB0d28gdHJlYXRtZW50IGdyb3VwcwogICAgLSBwbGFjZWJvIChvbmx5IGFudGliaW90aWNzKQogICAgLSB0cmFuc3BsYW50IChhbnRpYmlvdGljYSBmb2xsb3dlZCBieSBtaWNyb2JpYWwgdHJhbnNwbGFudCkuCiAgICAtIFRoZSBtaWNyb2Jpb21lIGlzIHNhbXBsZWQgNiB3ZWVrcyB1cG9uIHRoZSB0cmVhdG1lbnQKICAgIC0gVGhlIHJlbGF0aXZlIGFidW5kYW5jZSBvZiAqU3RhcGh5bG9jb2NjdXMgc3BwLiogb24KICAgICAgKkNvcnluZWJhY3Rlcml1bSBzcHAuKiArICpTdGFwaHlsb2NvY2N1cyBzcHAuKiBpbiB0aGUKICAgICAgbWljcm9iaW9tZSBpcyBtZWFzdXJlZCB2aWEgREdHRSAoKkRlbmF0dXJpbmcgR3JhZGllbnQgR2VsCiAgICAgIEVsZWN0cm9waG9yZXNpcyopLgoKIyBHb2FsCgpUaGUgb3ZlcmFyY2hpbmcgZ29hbCBvZiB0aGlzIHJlc2VhcmNoIHdhcyB0byBhc3Nlc3MgaWYgdGhlIHJlbGF0aXZlIGFidW5kYW5jZQoqU3RhcGh5bG9jb2NjdXMgc3BwLioKaW4gdGhlIG1pY3JvYmlvbWUgb2YgdGhlIGFybXBpdCBpcyBhZmZlY3RlZCBieSB0cmFuc3BsYW50aW5nIHRoZSBtaWNyb2Jpb21lLgpUbyB0aGlzIGVuZCB0aGUgcmVzZWFyY2hlcnMgcmFuZG9taXplZCBwYXRpZW50cyB0byB0d28gdHJlYXRtZW50OgpBIHRyZWF0bWVudCB3aXRoIGFudGliaW90aWNzIG9ubHkgYW5kIGEgdHJlYXRtZW50IHdpdGgKYW50aWJpb3RpY3MgYW5kIGEgbWljcm9iaWFsIHRyYW5zcGxhbnQuCgpJbiB0aGUgdHV0b3JpYWwgb24gaHlwb3RoZXNlcyB0ZXN0aW5nIHdlIHdpbGwgdXNlIGEgZm9ybWFsIHN0YXRpc3RpY2FsIHRlc3QgdG8KZ2VuZXJhbGl6ZSB0aGUgcmVzdWx0cyBmcm9tIHRoZSBzYW1wbGUgdG8gdGhhdCBvZiB0aGUgcG9wdWxhdGlvbi4KCiMgSW1wb3J0IHRoZSBkYXRhc2V0CgpgYGB7ciwgbWVzc2FnZT1GQUxTRX0KIyBMb2FkIHRoZSBsaWJyYXJpZXMKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKSW1wb3J0IHRoZSBkYXRhCgpgYGB7cn0KYXAgPC0gcmVhZF9jc3YoImh0dHBzOi8vcmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbS9zdGF0T21pY3MvUFNMU0RhdGEvbWFpbi9hcm1waXQuY3N2IikKYGBgCgpgYGB7cn0KZ2xpbXBzZShhcCkKYGBgCgojIERhdGEgRXhwbG9yYXRpb24KCkEgY3J1Y2lhbCBmaXJzdCBzdGVwIGluIGEgZGF0YSBhbmFseXNpcyBpcyB0byB2aXN1YWxpemUgYW5kIHRvIGV4cGxvcmUgdGhlIHJhdwpkYXRhLgoKYGBge3J9CmFwICU+JQogIGdncGxvdChhZXMoeCA9IHRydCwgeSA9IHJlbCwgZmlsbCA9IHRydCkpICsKICBnZW9tX2JveHBsb3Qob3V0bGllci5zaGFwZSA9IE5BKSArCiAgZ2VvbV9wb2ludChwb3NpdGlvbiA9ICJqaXR0ZXIiKSArCiAgeWxhYigicmVsYXRpdmUgYWJ1bmRhbmNlICglKSIpICsKICB4bGFiKCJ0cmVhdG1lbnQgZ3JvdXAiKSArCiAgc3RhdF9zdW1tYXJ5KAogICAgZnVuID0gbWVhbiwgZ2VvbSA9ICJwb2ludCIsCiAgICBzaGFwZSA9IDUsIHNpemUgPSAzLCBjb2xvciA9ICJibGFjayIsCiAgKQpgYGAKCldlIGNsZWFybHkgc2VlIHRoYXQsIG9uIGF2ZXJhZ2UsIHRoZSBzdWJqZWN0cyB3aG8gaGFkIGEKbWljcm9iaWFsIHRyYW5zcGxhbnQgaGF2ZSBhIGhpZ2hlciByZWxhdGl2ZSBhYnVuZGFuY2Ugb2YKU3RhcGh5bG9jb2NjdXMgc3BwLiBCdXQgaXMgdGhpcyBkaWZmZXJlbmNlIHN0YXRpc3RpY2FsbHkKKnNpZ25pZmljYW50KiBzbyB0aGF0IHdlIGNhbiBnZW5lcmFsaXplZCB3aGF0IHdlIG9ic2VydmUKaW4gdGhlIHNhbXBsZSB0byB0aGUgcG9wdWxhdGlvbj8KCldlIGNhbiB0ZXN0IHRoaXMgd2l0aCBhbiB1bnBhaXJlZCwgdHdvLXNhbXBsZSB0LXRlc3QsIHdoaWNoIGZhbHNpZmllcyB0aGUgbnVsbApoeXBvdGhlc2lzIHRoYXQgdGhlcmUgaXMgb24gYXZlcmFnZSBubyBkaWZmZXJlbmNlIGluIHJlbGF0aXZlIGFidW5kYW5jZSBvZgoqU3RhcGh5bG9jb2NjdXMqIGluIHRoZSBhcm1waXQgbWljcm9iaW9tZSBiZXR3ZWVuIHRoZSB0cmFuc3BsYW50IGFuZCB0aGUKcGxhY2VibyBncm91cCBhZ2FpbnN0IHRoZSBhbHRlcm5hdGl2ZSBoeXBvdGhlc2lzIHRoYXQgdGhlcmUgaXMgYSBkaWZmZXJlbmNlCmluIGF2ZXJhZ2UgYWJ1bmRhbmNlIG9mICpTdGFwaHlsb2Njb2N1cyogaW4gdGhlIGFybXBpdCBtaWNyb2Jpb21lIGJldHdlZW4KdGhlIHRyYW5zcGxhbnQgYW5kIHBsYWNlYm8gdHJlYXRtZW50LgoKQnV0LCBiZWZvcmUgd2UgY2FuIHN0YXJ0IHRoZSBhbmFseXNpcywgd2UgbXVzdCBjaGVjayBpZiBhbGwgYXNzdW1wdGlvbnMgdG8KcGVyZm9ybSBhIHQtdGVzdCBhcmUgbWV0LgoKIyBBbmFseXNpcwoKIyMgQ2hlY2sgdGhlIGFzc3VtcHRpb25zCgoxLiBUaGUgb2JzZXJ2YXRpb25zIGFyZSBpbmRlcGVuZGVudC4gVGhpcyBoYXMgdG8gYmUKZ3VhcmFudGVlZCBieSB0aGUgZGVzaWduLgoKMi4gVGhlIGRhdGEgKHJlbCkgYXJlIG5vcm1hbGx5IGRpc3RyaWJ1dGVkIGluIGVhY2ggb2YgdGhlIGdyb3VwcwoKMy4gVGhlIHZhcmlhYmlsaXR5IHdpdGhpbiBib3RoIGdyb3VwcyBpcyBzaW1pbGFyLgoKVG8gY2hlY2sgdGhlIG5vcm1hbGl0eSBhc3N1bXB0aW9uLCB3ZSB3aWxsIHVzZSBRUSBwbG90cy4KCmBgYHtyfQphcCAlPiUKICBnZ3Bsb3QoYWVzKHNhbXBsZSA9IHJlbCkpICsKICBnZW9tX3FxKCkgKwogIGdlb21fcXFfbGluZSgpICsKICBmYWNldF9ncmlkKGNvbHMgPSB2YXJzKHRydCkpCmBgYAoKV2UgY2FuIHNlZSB0aGF0IGFsbCBvZiB0aGUgZGF0YSBsaWVzIG5pY2VseSBhcm91bmQgdGhlIHF1YW50aWxlLXF1YW50aWxlCmxpbmUgKGJsYWNrIGxpbmUpLiBBcyBzdWNoLCB3ZSBtYXkgYXNzdW1lIHRoYXQgb3VyIGRhdGEgYXJlIG5vcm1hbGx5IGRpc3RyaWJ1dGVkLgoKRm9yIHRoZSB0aGlyZCBhc3N1bXB0aW9uLCB3ZSBtdXN0IGNvbXBhcmUgdGhlIHdpdGhpbi1ncm91cAp2YXJpYWJpbGl0eSBvZiBib3RoIGdyb3Vwcy4gV2UgY2FuIGRvIHRoaXMgdmlzdWFsbHk6CgpgYGB7cn0KYXAgJT4lCmdncGxvdChhZXMoeCA9IHRydCwgeSA9IHJlbCkpICsKICBnZW9tX2JveHBsb3Qob3V0bGllci5zaGFwZSA9IE5BKSArCiAgZ2VvbV9wb2ludChwb3NpdGlvbiA9ICJqaXR0ZXIiKSArCiAgeWxhYigicmVsYXRpdmUgYWJ1bmRhbmNlICglKSIpICsKICB4bGFiKCJ0cmVhdG1lbnQgZ3JvdXAiKSArCiAgc3RhdF9zdW1tYXJ5KAogICAgZnVuID0gbWVhbiwgZ2VvbSA9ICJwb2ludCIsCiAgICBzaGFwZSA9IDUsIHNpemUgPSAzLCBjb2xvciA9ICJibGFjayIsCikKYGBgCgpIZXJlIHdlIGNhbiBzZWUgdGhhdCB0aGUgaW50ZXJxdWFydGlsZSByYW5nZSBpcyBhcHByb3hpbWF0ZWx5IGVxdWFsIGZvciBncm91cHMuCgpBcyBhbGwgdGhyZWUgYXNzdW1wdGlvbnMgYXJlIG1ldCB3ZSBtYXkgY29udGludWUgd2l0aApwZXJmb3JtaW5nIHRoZSB1bnBhaXJlZCB0d28tc2FtcGxlIHQtdGVzdC4KCiMjIEh5cG90aGVzaXMgdGVzdAoKYGBge3J9Cm91dHB1dCA8LSB0LnRlc3QoCiAgcmVsIH4gdHJ0LAogIGRhdGEgPSBhcCwKICBjb25mLmxldmVsID0gMC45NSwKICB2YXIuZXF1YWwgPSBUUlVFCikKb3V0cHV0CmBgYAoKIyMgQ29uY2x1c2lvbgoKT24gYXZlcmFnZSB0aGUgcmVsYXRpdmUgYWJ1bmRhbmNlIG9mICpTdGFwaHlsb2NvY2N1cyBzcHAuKiBpbiB0aGUgbWljcm9iaW9tZSBvZiB0aGUgYXJtcGl0IGluIHRoZSB0cmFuc3BsYW50IGdyb3VwIGlzIGV4dHJlbWVseSBzaWduaWZpY2FudGx5IGRpZmZlcmVudCBmcm9tIHRoYXQgaW4gdGhlIHBsYWNlYm8gZ3JvdXAgKCRwPDwwLjAwMSQpLiBUaGUgcmVsYXRpdmUgYWJ1bmRhbmNlIG9mICpTdGFwaHlsb2NvY2N1cyBzcHAuKiBpcyBvbiBhdmVyYWdlIGByIHJvdW5kKGRpZmYodC50ZXN0KHJlbH50cnQsZGF0YT1hcCx2YXIuZXF1YWw9VFJVRSkkZXN0aW1hdGUpLDEpYCUgbGFyZ2VyIGluIHRoZSB0cmFuc3BsYW50IGdyb3VwIHRoYW4gaW4gdGhlIHBsYWNlYm8gZ3JvdXAgKDk1XCUgQ0kgW2ByIHBhc3RlKGZvcm1hdCgtdC50ZXN0KHJlbH50cnQsZGF0YT1hcCx2YXIuZXF1YWw9VFJVRSkkY29uZi5pbnRbMjoxXSxkaWdpdHM9Mixuc21hbGw9MSksY29sbGFwc2U9IiwiKWBdJSkuCgojIERhdGEgQW5hbHlzaXMgd2l0aCBhIGxpbmVhciBtb2RlbAoKIyMgRml0IHRoZSBtb2RlbAoKV2UgZG8gbm90IG5lZWQgdG8gY2hlY2sgdGhlIGFzc3VtcHRpb25zIGFnYWluIGJlY2F1c2Ugd2UgaGF2ZSBkb25lIHRoaXMgYWxyZWFkeSBpbiB0aGUgZGF0YSBleHBsb3JhdGlvbiB3aGVuIGFuYWx5c2luZyB0aGUgZGF0YSB1c2luZyB0aGUgdC10ZXN0LgoKYGBge3J9Cm1vZCA8LSBsbShyZWwgfiB0cnQsIGRhdGEgPSBhcCkKbW9kCmBgYAoKIyMgSW5mZXJlbmNlCgoxLiBoeXBvdGhlc2lzIHRlc3QKCmBgYHtyfQpzdW1tYXJ5KG1vZCkKYGBgCgoyLiBjb25maWRlbmNlIGludGVydmFsCgpgYGB7cn0KY29uZmludChtb2QpCmBgYAoKMy4gSG93IGRvIHlvdXIgcmVzdWx0cyBjb21wYXJlIHRvIHRoZSBhbmFseXNpcyB3aXRoIHRoZSB0LXRlc3QuCgo0LiBGb3JtdWxhdGUgeW91ciBjb25jbHVzaW9uIGJhc2VkIG9uIHRoZSBvdXRwdXQgb2YgdGhlIGxpbmVhciBtb2RlbC4K