---
title: "Resemblance of Legionella Untreated Controls"
author: "COST: Biofilm Regulatory Toolbox"
date: "8 Sept 2026"
output: word_document
---

Showing Resemblance of Legionella untreated control biofilms generated by an ISBR reactor at CBE.  These data were published by Bodle et al 2026: 

Bodle et al. Assessing a method for single- and multidomain biofilm growth in a novel biofilm reactor: shear stress, ruggedness, repeatability, and reproducibility. Biofilm 100372, 2026


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

# Get the Data

```{r}
d = read.csv("Parker_10-11_LegionellaControlData.csv",stringsAsFactors = TRUE)
summary(d)
```

# Plot the Data
```{r}
plot(d$Exp,log10(d$CFU.cm2),xlab="Experiment",ylab="log(CFU/cm2)")
```


# Mixed effects model

## Fit LMM to data
```{r}
m = lmer(log10(CFU.cm2) ~ (1|Exp),data=d)
summary(m)
confint(m)
```


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

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



## Summarize variance components
```{r}
VarCorr(m)
vc = as.data.frame(VarCorr(m))$vcov 

# Proportion of variance components
vc/sum(vc)
```

## Repeatability SD

Repeatablity SD = 0.5 log is considered acceptable for CFU assays.

```{r}
Sr = sqrt(sum(vc))
Sr
```


# Calculations by hand 
Mimicking what is done in Excel

## Get a CI for the mean
This CI will be different than what profile CI from lmer( ) 

```{r}
# Calculate the mean for each experiment 
mean.Exp = tapply(log10(d$CFU.cm2),d$Exp,mean)
mean.Exp

# Get confidence interval with a t-test
t.test(mean.Exp)
```

## Get the repeatability SD and variance components
```{r}
# Calculate the variance for each experiment 
Var.Exp=tapply(log10(d$CFU.cm2),d$Exp,var)

# Variance Component
Vcoupon = mean(Var.Exp)
Vcoupon
vc[2]   # Check against lmer()

Vexp = var(mean.Exp) - Vcoupon/length(mean.Exp)
Vexp
vc[1]   # Check against lmer()

# Repeatability SD
sqrt(Vexp + Vcoupon)
Sr      # Check against lmer()
```

# Assess Outlier Influence
Refit the LMM to the data without the outlier and compare results.


```{r}
m = lmer(log10(CFU.cm2) ~ (1|Exp),data=d[-3,])
summary(m)
confint(m)

# Check model fit - constant variance and normality of coupon LDs in each experiment
diagRES(m,res="p",num.panes=4,tests=TRUE)

# Check normality of experiment mean LDs
re=ranef(m)$Exp
qqnorm(re[,1])  
qqline(re[,1])

# Repeatability SD
sqrt(sum(as.data.frame(VarCorr(m))$vcov))
```








```{r}