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

Showing Responsiveness of the MBEC to different concentrations of an oxidizer against Pseudomonas biofilms.  Data generated by a single lab from an 8-laboratory study and reported by Parker et al 2019: 

Parker, mbec ...


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

# Get Data
```{r}
d = read.csv("MBEC_LRforR_Lab1.csv",stringsAsFactors = TRUE)
summary(d)

# Pick a Single Chemical
levels(d$Chemical)
dt = d[d$Chemical=="Oxidizer",]
```

# Plot the LRs for all Chemicals

## With a Linear Fit
```{r}
ggplot(d,aes(y=LR,x=log2.Conc.,color=Chemical,shape=Chemical, linetype=Chemical))+geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6),size=2)+geom_smooth(method="lm",se=FALSE)
```

## With Smoother Fit
```{r}
ggplot(d,aes(y=LR,x=log2.Conc.,color=Chemical,shape=Chemical, linetype=Chemical))+geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6),size=2)+geom_smooth(method="loess",se=FALSE)
```


# How to Model A Single Chemical Over Multiple Experiments

Check for Responsiveness to Concentration Increase  (a quantitative covariate).

## Plot the LRs by Concentration with Linear Fit
```{r}
dt$Day = as.factor(dt$Day)
ggplot(dt,aes(y=LR,x=log2.Conc.,color=Day,shape=Day, linetype=Day))+geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6),size=4)+geom_smooth(method="lm",se=FALSE)
```

## Plot the LRs by Concentration with Smoother Fit
```{r}
dt = d[d$Chemical=="Oxidizer",]
dt$Day = as.factor(dt$Day)
ggplot(dt,aes(y=LR,x=log2.Conc.,color=Day,shape=Day, linetype=Day))+geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6),size=4)+geom_smooth(method="loess",se=FALSE)
```


# LMM Estimate Slope Pooled Overall Experiments

## Fit the LMM
```{r}
mean(dt$log2.Conc.)
m = lmer(LR ~ I(log2.Conc.-mean(log2.Conc.)) + (1|Day),data=dt)
anova(m)
summary(m)
confint(m)
```

## Plot LMM Predictions with Data

```{r}
plot(dt$log2.Conc,dt$LR,lwd=2, xlab="Concentration",ylab="LR")
p = predict(m,re.form=~0)
lines(unique(dt$log2.Conc),unique(p),col="red",lwd=2)

plot(allEffects(m))
```



## Check Model Assumptions 
Check constant variance and normality of LRs
```{r}
diagRES(m,res="p",num.panes=4,tests=FALSE)

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


# GAM Smoother Regression

Compare the Chemicals (a categorical factor)


## Fit the GAM
```{r}
g = gamm(LR ~ ti(log2.Conc.),random=list(Day=~1),data=dt)
summary(g$gam)
confint(g$gam)
```

## Plot GAM predictions with data
```{r}
plot(dt$log2.Conc,dt$LR,lwd=2, xlab="Concentration",ylab="LR")
pg = predict(g$gam,re.form=~0)
lines(unique(dt$log2.Conc),unique(pg),col="red",lty=1,lwd=2)
```

## Plot comparison of GAM and LMM fits

```{r}
plot(dt$log2.Conc,dt$LR,lwd=2, xlab="Concentration",ylab="LR")
pg = predict(g$gam,re.form=~0)
lines(unique(dt$log2.Conc),unique(p),col="red",lwd=2)
lines(unique(dt$log2.Conc),unique(pg),col="green",lty=2,lwd=2)
legend("topleft",c("LMM","GAM"),col=c("red","green"),lty=1:2)
```

## Check Model Assumptions 
Check constant variance and normality of LRs
```{r}
diagRES(g$lme,res="p",num.panes=4,tests=FALSE)
gam.check(g$gam)
plot.gam(g$gam)
```