---
title: "Ruggedness of STM"
author: "COST: Biofilm Regulatory Toolbox"
date: "8 Sept 2026"
output: word_document
---

Showing Ruggedness of the STM to different treatments against Pseudomonas biofilms in the CDC reactor.  Data generated by CBE 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("STM_Ruggedness.csv",stringsAsFactors = TRUE)
d$Power_Time=factor(d$Power_Time,levels=c("20_25", "20_35", "80_30", "200_25", "200_35"))   # Get the order correct for ggplot()
summary(d)
```

# Plot the LDs vs Covariates 

View the data a few different ways.

```{r}
ggplot(d, aes(x = Degas, y = LD, group = Power_Time, color = as.factor(Power), shape = as.factor(Time))) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6), 
    size = 3) 

ggplot(d, aes(x = as.factor(Power), y = LD, group = Time_Degas, color = as.factor(Time), shape = Degas)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6), 
    size = 3) 

ggplot(d, aes(x = as.factor(Time), y = LD, group = Power_Degas, color = as.factor(Power), shape = Degas)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6), 
    size = 3)
```


# Assess Ruggedness

Centering at SOP values is crucial for interpretation of the interactions!


## Fit Mixed Effects Multiple Regression w All interactions
```{r}
m = lmer(LD ~ Degas*I(Power-80)*I(Time - 30) + (1|Exp),data=d)
summary(m)
confint(m)
```


## Investigate 3-way interaction
Model each Degas level separately, have to drop the RE due to Exp bc each degas was on a separate Exp:
```{r}
m.DegasN = lm(LD ~ I(Power-80)*I(Time - 30),data=d[d$Degas=="N",])
summary(m.DegasN)
confint(m.DegasN)

m.DegasY = lm(LD ~ I(Power-80)*I(Time - 30),data=d[d$Degas=="Y",])
summary(m.DegasY)
confint(m.DegasY)
```


## Plot Regression Predictions
```{r}
plot(allEffects(m.DegasN))
plot(allEffects(m.DegasY))
```


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