---
title: "Reproducibility and Repeatability of the STM"
author: "COST: Biofilm Regulatory Toolbox"
date: "8 Sept 2026"
output: word_document
---

Showing Reproducibility and Repeatability of the STM to 6 different treatments against Pseudomonas biofilms in the CDC reactor.  Data generated by 6 laboratories and reported by Goeres et al 2019: 

Goeres, Walker, Buckingham-Meyer, Lorenz, Summers, Fritz, Goveia, Dickerman, Schultz, and Parker.  Development, standardization, and validation of a biofilm efficacy test: the single tube method.  Journal of Microbiological Methods.  165, 2019


# Housekeeping
```{r}
library(lme4)    # lmer()
source("diagRES.r")
library(lmerTest)
library(effects)
library(ggplot2)
```

# Get Data
```{r}
d = read.csv("Parker_10-11_STM_LRs_6Labs.csv",stringsAsFactors = TRUE)
d$Treatment = factor(d$Treatment,levels=c("Chlorine_LOW","Chlorine_HIGH", "Phenol_LOW", "Phenol_HIGH", "Quat-alcohol_LOW", "Quat-alcohol_HIGH"))
summary(d)
```

# Plot the data
```{r}
ggplot(d, aes(x = Chemical, y = LR, color = Chemical, group = Conc, shape = Conc)) +
  geom_point(
    position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6), 
    size = 3, 
    alpha = 1) +
stat_summary(
    fun = "mean", 
    geom = "crossbar", 
    position = position_dodge(width = 0.6), 
    size = 0.5)
```

# Get Reproducibility for each Chemical and Concentration Combination

## Fit the model
```{r}
m = lmer(LR ~ (1|Lab),data=d[d$Chemical=="Chlorine"&d$Conc=="HIGH",])
summary(m)
confint(m)  # for 2-sided 95% CI
confint(m, level=0.9)  # for upper 1-sided 95% CI 
```

## Reproducibility and Repeatability
```{r}
VarCorr(m)
vc = as.data.frame(VarCorr(m))$vcov 

# Proportion of variance components
vc/sum(vc)

# Reproducibility of the LRs
sqrt(sum(vc))

# Repeatability of the LRs (estimated by pooling over all labs)
sqrt(vc[2])
```

## Check model assumptions
Check constant variance and normality of coupon LDs in each experiment
```{r}
diagRES(m,res="p",num.panes=4,tests=FALSE)

# Check normality of lab mean LDs
re=ranef(m)$Lab
qqnorm(re[,1])  
qqline(re[,1])
```