---
title: "Experiment 2 stats"
author: "Jasmine Clarkson"
date: "15/06/2022"
output: html_document
---
This excludes the one outlier - re-runs the analyses without the long decompression run for writing up purposes. 
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
#load in the packages
library(lme4)
library(lmerTest)
library(emmeans)
library(pbkrtest)
library(multcomp)
library(car)
library(tidyverse)
library(glmmTMB)
library(ggpubr)
library(DHARMa)
library(psych)
library(viridis)
library(fishualize)
library(ggplot2)
library(ordinal)
library (descr)
library(readxl)
```

## load in the data

```{r}
#load in the dataset
library(readxl)
Exp2 <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\Exp2_fullbehaviour.xlsx")
# allows you to view the data
View(Exp2)

```


```{r}
Exp2$Treatment <- as.factor(Exp2$Treatment)
Exp2$Intervention <- as.factor(Exp2$Intervention)
Exp2$Batch <- as.factor(Exp2$Batch)
Exp2$Cage <- as.factor(Exp2$Cage)
Exp2$CycleLengthS <- as.numeric(Exp2$CycleLengthS)
```

# Proportion data - for key behaviours 

## Proportion of mice ear scratching

```{r}
sum_data <- Exp2 %>% group_by(Treatment, Intervention) %>% summarise(xfrequency = sum(EarScratchBinary=='1'))
sum_data
```



```{r}
logit_1 <- glm(EarScratchBinary ~ Intervention + Treatment + Intervention:Treatment, family=binomial(), data = Exp2)
summary(logit_1)
```


```{r}
simulationOutput <- simulateResiduals(fittedModel = logit_1, plot = T)
```

```{r}
Anova(logit_1)
```


```{r}
emm1 <- emmeans(logit_1, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(logit_1, list(pairwise ~ Intervention), adjust ="tukey")
emm3 <- emmeans(logit_1, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# proportion of mice facial grimacing - changes in ear posture 

```{r}
logit_2 <- glm(TotalGrimScore ~ RH + Temp + Intervention + Treatment + Intervention:Treatment, family=binomial(), data = Exp2)
summary(logit_2)
```

```{r}
simulationOutput <- simulateResiduals(fittedModel = logit_2, plot = T)
```


```{r}
Anova(logit_2)
emm1 <- emmeans(logit_2, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(logit_2, list(pairwise ~ Intervention), adjust ="tukey")
emm3 <- emmeans(logit_2, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# Proportion of mice head flicking

```{r}
logit_3 <- glm(HeadFlickBin ~ RH + Temp + Intervention + Treatment + Intervention:Treatment, family=binomial(), data = Exp2)
summary(logit_3)
```

```{r}
simulationOutput <- simulateResiduals(fittedModel = logit_3, plot = T)

```

```{r}
Anova(logit_3)
emm1 <- emmeans(logit_3, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(logit_3, list(pairwise ~ Intervention), adjust ="tukey")
emm3 <- emmeans(logit_3, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```



# Latency data  

# latency to ear scratch 

```{r}
# Fit a linear mixed model with weight as covariate (no effect so moved to random effect), RH and temp no effect so removed to improve fit 
m1 <-lmer(EarScratchLat ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m1)
```

```{r}
# looking at the fit of the residuals model - not a great fit 
simulationOutputm1<- simulateResiduals(fittedModel = m1, plot = T)
plot(simulationOutputm1)
hist(simulationOutputm1)
```

```{r}
# model output 
anovam1 <- Anova(m1)
anovam1
```

```{r}
emm1 <- emmeans(m1, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m1, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m1, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# latency to head flick

```{r}
m2 <-glmmTMB(HeadFlick ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage),  data=Exp2, family=nbinom2()) 
summary(m2)
```

```{r}
simulationOutputm2<- simulateResiduals(fittedModel = m2, plot = T)
plot(simulationOutputm2)
hist(simulationOutputm2)
```

```{r}
Anova(m2)
emm1 <- emmeans(m2, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m2, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m2, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# Latency to gasping 

```{r}
m3 <-lmer(GaspingLatency ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m3)
```

```{r}
simulationOutputm3<- simulateResiduals(fittedModel = m3, plot = T)
plot(simulationOutputm3)
hist(simulationOutputm3)
```

```{r}
Anova(m3)
emm1 <- emmeans(m3, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m3, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m3, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")

```

# latency to facial grimace - should do proprtions of mice showing a facial grimace as power for this model is small. 

```{r}
m4.1 <-lmer(FacialGrimaceLat ~  Treatment + Intervention + Treatment:Intervention  + (1|Batch/Cage), data=Exp2) 
summary(m4.1)
```

```{r}
simulationOutputm4.1<- simulateResiduals(fittedModel = m4.1, plot = T)
plot(simulationOutputm4.1)
hist(simulationOutputm4.1)
```


```{r}
Anova(m4.1)
emm1 <- emmeans(m4.1, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m4.1, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m4.1, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```



# latency to LOB 

```{r}
m4 <-lmer(LOBLatency ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m4)
```

```{r}
simulationOutputm4<- simulateResiduals(fittedModel = m4, plot = T)
plot(simulationOutputm4)
hist(simulationOutputm4)
```

```{r}
Anova(m4)
emm1 <- emmeans(m4, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m4, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m4, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## latency to LOP 

```{r}
m5 <-lmer(LOPLatency ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m5)
```

```{r}
simulationOutputm5<- simulateResiduals(fittedModel = m5, plot = T)
plot(simulationOutputm5)
hist(simulationOutputm5)
```

```{r}
Anova(m5)
emm1 <- emmeans(m5, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m5, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m5, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

## motionless latency 

```{r}
m6 <-lmer(MotionlessLatency ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m6)
```


```{r}
simulationOutputm6<- simulateResiduals(fittedModel = m6, plot = T)
plot(simulationOutputm6)
hist(simulationOutputm6)
```


```{r}
Anova(m6)
emm1 <- emmeans(m6, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m6, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m6, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# Latency to Agonal gasping 


```{r}
m7 <-lmer(AgonalLatency ~  Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m7)
```

```{r}
simulationOutputm7<- simulateResiduals(fittedModel = m7, plot = T)
plot(simulationOutputm7)
hist(simulationOutputm7)
```

```{r}
Anova(m7)
emm1 <- emmeans(m7, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m7, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m7, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## convulsion latency

```{r}
m8 <-lmer(ConvlLatency ~ RH + Temp + Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m8)
```

```{r}
simulationOutputm8<- simulateResiduals(fittedModel = m8, plot = T)
plot(simulationOutputm8)
hist(simulationOutputm8)
```

```{r}
Anova(m8)
emm2<- emmeans(m8, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")
```

## latency to head up

```{r}
m9 <-glmmTMB(HeadUpLat ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m9)
```

```{r}
simulationOutputm9<- simulateResiduals(fittedModel = m9, plot = T)
plot(simulationOutputm9)
hist(simulationOutputm9)
```

```{r}
Anova(m9)
emm1 <- emmeans(m9, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m9, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m9, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

## recumbent latency 

```{r}
m10 <-glmmTMB(RecumbentLat ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m10)
```



```{r}
simulationOutputm10<- simulateResiduals(fittedModel = m10, plot = T)
plot(simulationOutputm10)
hist(simulationOutputm10)
```

```{r}
Anova(m10)
emm1 <- emmeans(m10, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m10, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m10, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

## Duration data 

# ear scratching duration 

```{r}
m11 <-glmmTMB(EarScratchDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family = nbinom2()) 
summary(m11)
```

```{r}
simulationOutputm11<- simulateResiduals(fittedModel = m11, plot = T)
plot(simulationOutputm11)
hist(simulationOutputm11)
```

```{r}
Anova(m11)
emm1 <- emmeans(m11, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m11, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m11, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# Grooming bout duration 

```{r}
m12 <-glmmTMB(GroomingBoutDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m12)
```


```{r}
simulationOutputm12<- simulateResiduals(fittedModel = m12, plot = T)
plot(simulationOutputm12)
hist(simulationOutputm12)
```


```{r}
Anova(m12)
emm1 <- emmeans(m12, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m12, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m12, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# last agonal gasp - latency 

```{r}
m13 <-glmmTMB(LastAgonal ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m13)
```


```{r}
simulationOutputm13<- simulateResiduals(fittedModel = m13, plot = T)
plot(simulationOutputm13)
hist(simulationOutputm13)
```


```{r}
Anova(m13)
emm1 <- emmeans(m13, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m13, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m13, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

## Count/frequency data 

## head flick freq

for this we need to decide if we account for differences in cycle length therefore need to include cycle length as an offset in the models. 

```{r}
m14 <- lmer(HeadFlickFrq ~ RH + Temp +Treatment + Intervention + Treatment:Intervention  + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m14)
 #offset(log(CycleLengthS))
```

```{r}
simulationOutputm14<- simulateResiduals(fittedModel = m14, plot = T)
plot(simulationOutputm14)
hist(simulationOutputm14)
```

```{r}
Anova(m14)
emm1 <- emmeans(m14, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m14, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m14, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# the above analysis incorporates all individuals including those that didnt do it -
therefore want to perform a subset analysis focused on the mean frequency of those individsuals that performed the behaviour

```{r}
# filter out the zeros 
s1d <- subset(Exp2, HeadFlickBin=='1')

s1 <- lmer(HeadFlickFrq ~ RH + Temp +Treatment + Intervention + Treatment:Intervention  + (1|Weight) + (1|Batch/Cage), data=s1d) 
summary(s1)
```


```{r}
simulationOutputs1<- simulateResiduals(fittedModel = s1, plot = T)
plot(simulationOutputs1)
hist(simulationOutputs1)
```


```{r}
Anova(s1)
emm1 <- emmeans(s1, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(s1, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(s1, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# grooming bout frequency
number of 'normal' full grooming events

```{r}
m15 <- lmer(GroomingBFrq ~ RH + Temp + Treatment + Intervention + Treatment:Intervention  + (1|Batch/Cage), data=Exp2) 
summary(m15)
```

```{r}
simulationOutputm15<- simulateResiduals(fittedModel = m15, plot = T)
plot(simulationOutputm15)
hist(simulationOutputm15)
```



```{r}
Anova(m15)
emm1 <- emmeans(m15, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m15, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m15, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## same as previous - filter out those animals that did not perfrom this behaviour 
 GROOMING BOUTS/EVENTS SUMMED TOGETHER

```{r}
# filter out the zeros 
s2d <- subset(Exp2, TotalGroomComp > 0)

s2 <- lmer(TotalGroomComp ~ RH + Temp +Treatment + Intervention + Treatment:Intervention  + (1|Weight) + (1|Batch/Cage), data=s2d) 
summary(s2)
```



```{r}
simulationOutputs2<- simulateResiduals(fittedModel = s2, plot = T)
plot(simulationOutputs2)
hist(simulationOutputs2)
```


```{r}
Anova(s2)
emm1 <- emmeans(s2, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(s2, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(s2, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```



## number of gasps 

```{r}
m16 <- lmer(Gasp ~ RH + Temp + Treatment + Intervention + Treatment:Intervention  + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m16)
```


```{r}
simulationOutputm16<- simulateResiduals(fittedModel = m16, plot = T)
plot(simulationOutputm16)
hist(simulationOutputm16)
```


```{r}
Anova(m16)
emm1 <- emmeans(m16, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m16, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m16, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# latency ataxia 

```{r}
m17 <-lmer(AtaxiaLatency ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m17)
```


```{r}
simulationOutputm17<- simulateResiduals(fittedModel = m17, plot = T)
plot(simulationOutputm17)
hist(simulationOutputm17)
```


```{r}
Anova(m17)
emm1 <- emmeans(m17, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m17, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m17, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```



## number of ear scratches (frequency) - this inlcudes animals including those that never performed this behaviour. 

```{r}
m18 <-glmmTMB(EarScratchFrq ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m18)
```


```{r}
simulationOutputm18<- simulateResiduals(fittedModel = m18, plot = T)
plot(simulationOutputm18)
hist(simulationOutputm18)
```


```{r}
Anova(m18)
emm1 <- emmeans(m18, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m18, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m18, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## filtering out the animals that did not perform ear scratching 

```{r}
# filter out the zeros 
s3d <- subset(Exp2, EarScratchBinary ==1)

s3 <- glmmTMB(EarScratchFrq ~ Treatment + Intervention + Treatment:Intervention  + (1|Weight) + (1|Batch/Cage), family=nbinom2, data=s3d) 
summary(s3)
```

```{r}
simulationOutputs3<- simulateResiduals(fittedModel = s3, plot = T)
plot(simulationOutputs3)
hist(simulationOutputs3)
```


```{r}
Anova(s3)
emm1 <- emmeans(s3, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(s3, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(s3, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## number of rears 

```{r}
m19 <-glmmTMB(RearFrq ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m19)
```

```{r}
simulationOutputm19<- simulateResiduals(fittedModel = m19, plot = T)
plot(simulationOutputm19)
hist(simulationOutputm19)
```


```{r}
Anova(m19)
emm1 <- emmeans(m19, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m19, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m19, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## number of digging bouts 


```{r}
m20 <-glmmTMB(DigFrq ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m20)
```

```{r}
simulationOutputm20<- simulateResiduals(fittedModel = m20, plot = T)
plot(simulationOutputm20)
hist(simulationOutputm20)
```

```{r}
Anova(m20)
emm1 <- emmeans(m20, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m20, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m20, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## number of head up postures

```{r}
m21 <-glmmTMB(HeadUpFrq ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m21)
```

```{r}
simulationOutputm21<- simulateResiduals(fittedModel = m21, plot = T)
plot(simulationOutputm21)
hist(simulationOutputm21)
```


```{r}
Anova(m21)
emm1 <- emmeans(m21, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m21, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m21, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## defecation 

```{r}
m22 <-glmmTMB(DefecateFrq ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m22)
```
```{r}
simulationOutputm22<- simulateResiduals(fittedModel = m22, plot = T)
plot(simulationOutputm22)
hist(simulationOutputm22)
```


```{r}
Anova(m22)
emm1 <- emmeans(m22, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m22, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m22, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


# number of agonal gasps 

```{r}
# wouldnt converge with random factors included
m23 <-glmmTMB(AgonalFrq ~ Treatment + Intervention + Treatment:Intervention, data=Exp2, family=nbinom2()) 
summary(m23)
```
```{r}
simulationOutputm23<- simulateResiduals(fittedModel = m23, plot = T)
plot(simulationOutputm23)
hist(simulationOutputm23)
```

```{r}
Anova(m23)
emm1 <- emmeans(m23, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m23, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m23, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## latency to apnoea 
```{r}
m24 <-lmer(ApnoeaLatency ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m24)
```

```{r}
simulationOutputm24<- simulateResiduals(fittedModel = m24, plot = T)
plot(simulationOutputm24)
hist(simulationOutputm24)
```


```{r}
Anova(m24)
emm1 <- emmeans(m24, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m24, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m24, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```



## exploration of the peridphery - total duration (as a measure of activity?) 

```{r}
m25 <-glmmTMB(ExplPerDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family = nbinom2()) 
summary(m25)
```


```{r}
simulationOutputm25<- simulateResiduals(fittedModel = m25, plot = T)
plot(simulationOutputm25)
hist(simulationOutputm25)
```


```{r}
Anova(m25)
emm1 <- emmeans(m25, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m25, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m25, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# head up duration 

```{r}
m26 <-glmmTMB(HeadUpDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m26)
```

```{r}
simulationOutputm26<- simulateResiduals(fittedModel = m26, plot = T)
plot(simulationOutputm26)
hist(simulationOutputm26)
```


```{r}
Anova(m26)
emm1 <- emmeans(m26, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m26, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m26, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## digging duration 

```{r}
m27 <-glmmTMB(DigDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2, family=nbinom2()) 
summary(m27)
```
```{r}
simulationOutputm27<- simulateResiduals(fittedModel = m27, plot = T)
plot(simulationOutputm27)
hist(simulationOutputm27)
testDispersion(simulationOutputm27)
```


```{r}
Anova(m27)
emm1 <- emmeans(m27, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m27, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m27, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

## convulsion duration (decompression only)


```{r}
m28 <-lmer(ConvulsionDur ~ Intervention  + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m28)
```


```{r}
simulationOutputm28<- simulateResiduals(fittedModel = m28, plot = T)
plot(simulationOutputm28)
hist(simulationOutputm28)
testDispersion(simulationOutputm28)
```


```{r}
Anova(m28)

emm2<- emmeans(m28, list(pairwise ~ Intervention), adjust = "tukey")

summary(emm2, type = "response")

```

## latency for abnormal breathing (taken shortest latency and pooled irrespective of fast, slow or irregular)

```{r}
m29 <-lmer(AbnormalCombined ~ RH + Temp + Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(m29)
```

```{r}
simulationOutputm29<- simulateResiduals(fittedModel = m29, plot = T)
plot(simulationOutputm29)
hist(simulationOutputm29)
testDispersion(simulationOutputm29)
```


```{r}
Anova(m29)
emm1 <- emmeans(m29, list(pairwise ~ Treatment), adjust = "tukey")
emm2<- emmeans(m29, list(pairwise ~ Intervention), adjust = "tukey")
emm3<- emmeans(m29, list(pairwise ~ Intervention|Treatment), adjust = "tukey")
              
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


## plots of the data 

# full behavioural overview from emmeans from the models. 

```{r}
#load in the dataset
library(readxl)
Emms <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\AllLatenciesEMMs.xlsx")
# allows you to view the data
View(Emms)

Behaviour<- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\FullBehav_Long.xlsx")
# allows you to view the data
View(Behaviour)
```

```{r}
Emms$Behaviour <- factor(Emms$Behaviour, levels = c("Head flick", "HeadUp", "Gasp", "Ataxia", "Ear posture", "Ear scratch", "LOB", "Recumbent", "Apnoea", "LOP", "Agonal", "Convulsion", "Last agonal", "Motionless"))
Emms$Treatment <- factor(Emms$Treatment, levels = c("CO2", "Decom", "Sham"))
Emms$Intervention <- factor(Emms$Intervention, levels = c("Analgesic", "Anxio", "Saline"))

Emms$facets = factor(Emms$Treatment, labels = c("CO[2]","Decompression", 'SHAM'))
```


```{r}
g1 <- ggplot(Emms,aes(x=Mean, y=Behaviour, group = Intervention, colour=Intervention)) + geom_point(size = 1.5, position = position_dodge(width = 0.8)) + geom_errorbar(mapping = aes(x=Mean, y=Behaviour, xmin= Mean -SE, xmax= Mean+SE), size=0.5, width=0.8, position = position_dodge(width = 1)) + theme(text = element_text(size = 14)) + labs(x="Latency (s)") + xlim(0,400)
g1

g2 <- g1 +facet_wrap(~facets, labeller = label_parsed) + theme(text = element_text(size = 16), axis.text.x = element_text(angle=360, hjust = 1), panel.background = element_rect(), panel.grid.major = element_line(linetype = 'solid',colour = "white"), axis.line = element_line(colour = "black"))
g2


mean_LOP <- data.frame(Treatment = c("CO[2]","Decompression", 'SHAM'), z = c(80.2, 273.7, NA))# derived from emmeans for mean LOP
mean_LOP$facets = factor(mean_LOP$Treatment, labels = c("CO[2]","Decompression", 'SHAM'))

g3 <- g2 +  geom_vline(data= mean_LOP, aes (xintercept = z), color="red")# can add more lines for cycle phases too just need to know ther mean latencies.  
g3

ggsave("EMMsBehaviouralOverview.tiff", units="mm", width=300, height=200, dpi=600, compression = 'lzw')

```

## pool across treatments but not intervention

Makes for a clearer graph for the sequence of events 
need to then annotate the graph with symbols if we find intervention effects 


```{r}

# load in new data from EMMs focused on mean Latencies for treatment only
Emms2 <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsAllLatenciesTreatmentOnly.xlsx")
# allows you to view the data
View(Emms2)
Emms2$facets = factor(Emms2$Treatment, labels = c("CO[2]","Decompression", 'Sham'))
Emms2$Behaviour <- factor(Emms2$Behaviour, levels = c("Head flick", "HeadUp", "Ear scratch", "Abnormal breathing","Gasp",  "Ataxia", "Ear posture","LOB", "Recumbent", "Apnoea", "LOP", "Agonal", "Convulsion", "Last agonal", "Motionless"), labels = c("Head flick", "Head Up", "Ear scratch", "Abnormal breathing","Gasp",  "Ataxia", "Ear posture","Loss of Balance", "Recumbent", "Apnoea", "Loss of Posture", "Agonal gasp", "Convulsion", "Last agonal gasp", "Motionless") )


g1 <- ggplot(Emms2,aes(x=Mean, y=Behaviour, group = Treatment, colour=Treatment)) + geom_point(size = 1.5, position = position_dodge(width = 0.8)) + geom_text(aes(label=n, family = n), nudge_x=55) + geom_errorbar(mapping = aes(x=Mean, y=Behaviour, xmin= Mean -SE, xmax= Mean+SE), size=0.5, width=0.8, position = position_dodge(width = 1)) + theme(text = element_text(size = 14)) + labs(x="Latency (s)") + xlim(0,420)+ scale_color_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))
g1

g2 <- g1 +facet_wrap(~facets, labeller = label_parsed) + theme( axis.text.x = element_text(angle=360, hjust = 1), panel.background = element_rect(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
g2


mean_LOP <- data.frame(Treatment = c("CO[2]","Decompression", 'Sham'), z = c(80.2, 273.7, NA))# derived from emmeans for mean ☻LOP
mean_LOP$facets = factor(mean_LOP$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

mean_cycle <-data.frame(Treatment = c("CO[2]","Decompression", 'Sham'), z = c(114, 303, 360))
mean_cycle$facets = factor(mean_cycle$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

g3 <- g2 +  geom_vline(data= mean_LOP, aes (xintercept = z), color="red")+ theme(panel.background = element_blank(), legend.position = "none") # can add more lines for cycle phases too just need to know ther mean latencies.  
g4 <-g3 +  geom_vline(data= mean_cycle, aes (xintercept = z), color="grey", linetype = "longdash")+ theme(panel.background = element_blank(), legend.position = "none")
g4

ggsave("EMMsBehaviouralOverview_treatmentNEW.tiff", units="mm", width=190, height=150, dpi=600, compression = 'lzw')
```



#latency bar charts for behaviours - respiratory distress 

```{r}
Gasp <- subset(Emms2, Behaviour=='Gasp')

r1 <- ggplot(Gasp[!is.na(Gasp$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Gasp")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,300) + theme_bw() 
r2 <- r1 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
r2

Agonal <- subset(Emms2, Behaviour=='Agonal gasp')

r3 <- ggplot(Agonal[!is.na(Agonal$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Agonal Gasp")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,400) + theme_bw() 
r4 <- r3 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
r4

Apnoea <- subset(Emms2, Behaviour=='Apnoea')

r5 <- ggplot(Apnoea[!is.na(Apnoea$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Apnoea")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,300) + theme_bw() 
r6 <- r5 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
r6


LastAg <- subset(Emms2, Behaviour=='Last agonal gasp')

r7 <- ggplot(LastAg[!is.na(LastAg$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to last agonal gasp")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,400) + theme_bw() 
r8 <- r7 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
r8
```

```{r}
# grid them together 
# grid these figures together - where all three treatment groups are present

figure <- ggarrange(r2, r6, r4, r8,  ncol = 2, nrow = 2, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D"))
figure

ggsave("Resp Latencies Talk.tiff", units="mm", width=250, height=200, dpi=600, compression = 'lzw')
```

## bar charts for latencies for death behaviours

```{r}
Atax <- subset(Emms2, Behaviour=='Ataxia')

h1 <- ggplot(Atax[!is.na(Atax$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Ataxia")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,300) + theme_bw() 
h2 <- h1 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h2

LOB <- subset(Emms2, Behaviour=='Loss of Balance')

h3 <- ggplot(LOB[!is.na(LOB$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Loss of Balance")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,300) + theme_bw() 
h4 <- h3 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h4

LOP <- subset(Emms2, Behaviour=='Loss of Posture')

h5 <- ggplot(LOP[!is.na(LOP$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Loss of Posture")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,300) + theme_bw() 
h6 <- h5 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h6

Conv <- subset(Emms2, Behaviour=='Convulsion')

h7 <- ggplot(Conv[!is.na(LOP$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Convulsion")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,400) + theme_bw() 
h8 <- h7 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h8

Recum <- subset(Emms2, Behaviour=='Recumbent')

h9 <- ggplot(Recum[!is.na(Recum$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Recumbent")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,400) + theme_bw() 
h10 <- h9 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h10

Mot <- subset(Emms2, Behaviour=='Motionless')

h11 <- ggplot(Mot[!is.na(Mot$Mean), ], aes(x=Treatment, y=Mean, fill = (Treatment), ymin =Mean-SE, ymax = Mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2,
                position = position_dodge(0.9))+ ggtitle("Latency to Motionless")+ ylab("Latency (s)") + labs(x="Terminal treatment",y=" Latency (s)") + ylim(0,400) + theme_bw() 
h12 <- h11 + theme(legend.position="none", axis.title.x = element_blank()) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
h12

```

```{r}
# grid them together 
# grid these figures together - where all three treatment groups are present

figure <- ggarrange(h2, h4, h6, h10, h8, h12, ncol = 3, nrow = 2, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D", "E", "F"))
figure

ggsave("Hypoxia Latencies Talk.tiff", units="mm", width=250, height=200, dpi=600, compression = 'lzw')
```


## mean pressure versus the behaviours seen in decompression 

```{r}
library(readxl)
Pressure<- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\DecomPressure.xlsx")
# allows you to view the data
View(Pressure)
```


```{r}
#plotting each profile 

p <- ggplot(data=Pressure, aes(x=Secs, y=Pressure)) +
 
  geom_point()+
  labs(title = "Decompression profiles", x="Time (s)",y="Pressure (mBar)") + ggeasy::easy_center_title()
p

# time shifted mouse 8.2 batch 2 where the valve wasn't open correctly until approx 90s in to the cycle
```



```{r}
library(plyr)
summary <- ddply(Pressure, c('Secs'),summarise,N = length(Pressure),
               mean = mean(Pressure),
               sd   = sd(Pressure),
               se   = sd / sqrt(N))
                
summary
write.table(summary, file = "J:/BAHCM/LAPSProject/Experiment 2 Mice/Analysis/MeanPressure.csv", sep=",", row.names=FALSE)

# can also look at this for each batch separately and plot accordingly - shows that they were veyr consistent. 
```




```{r}
p <- ggplot(data=summary, aes(x=Secs, y=mean)) +
  geom_line()+ geom_ribbon(aes(x = Secs, ymin = mean-sd, ymax = mean+sd),alpha=0.2, color='#29AF7FFF') +
  xlim(0,300) + (ylim(200, 1000))+
  labs(title = " Gradual decompression", x="Time (s)",y="Pressure (mBar)") + ggeasy::easy_center_title()+ 
  scale_y_continuous("Mean pressure (mBar +/- 95% CI)", sec.axis = sec_axis(~ (. -63) * 0.0224119530, name = expression("Atmos O"[2]* " equivalent (%)")))+ geom_vline(aes (xintercept = 60))+ geom_vline(aes (xintercept = 273.7, color="red"))

p1 <- p + geom_vline(xintercept = 76.9, colour="blue", linetype = "longdash") # adding mean ear scratching latency
p1
  
# still need to time sort the one profile above which should narrow the confint - done
```



## correlation between ear scratching and head flick - same individuals?

```{r}
ggscatter(Exp2, x = "EarScratchLat", y = "HeadFlick",col=("Treatment"),shape = ("Intervention"),
          xlab = 'Latency to ear scratch (s)', ylab = "Latency to head flick (s)") + geom_abline(intercept=0, slope=1)

# 8 sham animals did both and 13 decom animals 
```

# looking at the CO2 data profiles 

```{r}
library(readxl)
Cdata<- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\CO2Data.xlsx")
# allows you to view the data
View(Cdata)
```



```{r}
p1 <- ggplot(data=Cdata, aes(x=Secs, y=CO2)) +
 
  geom_point()+
  labs(title = "Carbon dioxide runs", x="Time (s)",y="% CO2") + ggeasy::easy_center_title()
p1
```


```{r}
library(plyr)
summary3 <- ddply(Cdata, c('Secs'),summarise,N = length(CO2),
               mean = mean(CO2),
               sd   = sd(CO2),
               se   = sd / sqrt(N))
                
summary3
write.table(summary3, file = "J:/BAHCM/LAPSProject/Experiment 2 Mice/Analysis/MeanCO2.csv", sep=",", row.names=FALSE)

# n=47 for starting time points because no data recorded for one mouse in batch 1 - recording started after 2 mins (120s).
```


```{r}
p2 <- ggplot(data=summary3, aes(x=Secs, y=mean)) +
  geom_line()+ geom_ribbon(aes(x = Secs, ymin = mean-sd, ymax = mean+sd),alpha=0.2, color='steelblue2') +
  xlim(0,300) + (ylim(200, 1000))+
  labs(title = "Carbon Dioxide Runs", x="Time (s)",y="% CO2") + ggeasy::easy_center_title()+ 
  scale_y_continuous("Mean CO2 (% +/- 95% CI)") + geom_vline(aes (xintercept = 60))+ geom_vline(aes (xintercept = 80.2, color="red"))

p2
```

## oxygen concentration 

```{r}
summary2 <- ddply(Cdata, c('Secs'),summarise,N = length(O2),
               mean = mean(O2),
               sd   = sd(O2),
               se   = sd / sqrt(N))
                
summary2
write.table(summary2, file = "J:/BAHCM/LAPSProject/Experiment 2 Mice/Analysis/MeanO2.csv", sep=",", row.names=FALSE)

p3 <- ggplot(data=summary2, aes(x=Secs, y=mean)) +
  geom_line()+ geom_ribbon(aes(x = Secs, ymin = mean-sd, ymax = mean+sd),alpha=0.2, color='steelblue2') +
  xlim(0,300) + (ylim(200, 1000))+
  labs(title = "Carbon Dioxide Runs", x="Time (s)",y="% O2") + ggeasy::easy_center_title()+ 
  scale_y_continuous("Mean O2 (% +/- 95% CI)") + geom_vline(aes (xintercept = 60))+ geom_vline(aes (xintercept = 80.2, color="red"))

p3

```

##for carbon dioxide - put CO2 and O2 on the same graph 

## Methods Figure of terminal treatment overview 

```{r}
coeff <- 4.4

p5 <- ggplot(data=summary3, aes(x=Secs, y=mean)) +
  geom_line()+ geom_ribbon(aes(x = Secs, ymin = mean-sd, ymax = mean+sd),alpha=0.2, color='#481567FF') +
  xlim(0,300) + (ylim(200, 1000))+
  labs(title = "Carbon Dioxide", x="Time (s)",y="% CO2") + ggeasy::easy_center_title()+ 
  scale_y_continuous(expression(" Mean CO"[2]* " (% +/- 95% CI)") , sec.axis = sec_axis(~ (.) / 4.4, name = expression(" O"[2]* " concentration (% +/- 95% CI)"))) + geom_vline(aes (xintercept = 60))+ geom_vline(aes (xintercept = 80.2, color="red")) + geom_line(data = summary2, aes(x=Secs, y=mean*coeff))+ geom_ribbon(data = summary2, aes(x = Secs, ymin = (mean-sd)*coeff, ymax = (mean+sd)*coeff),alpha=0.2, color='steelblue2')+ geom_label(label= "CO2" , x =28, y=30, colour = "white", fill='#481567FF') + geom_label(label= "O2" , x =30, y=65, colour = "white", fill= 'steelblue2')

p5
```

  


# map this onto same graph as decompression profiles? - see the similarities of decom equiv o2 conc with co2 reduction in o2


```{r}
coeff <- 0.021
p4 <- ggplot(data=summary, aes(x=Secs, y=mean)) +
  geom_line()+ geom_ribbon(aes(x = Secs, ymin = mean-sd, ymax = mean+sd),alpha=0.2, color='steelblue2') +
  xlim(0,300) + (ylim(200, 1000))+
  labs(title = "Decompression vs CO2 comparison", x="Time (s)",y="Pressure (mBar)") + ggeasy::easy_center_title()+ 
  scale_y_continuous("Mean pressure (mBar +/- 95% CI)", sec.axis = sec_axis(~ (. -63) * 0.0224119530, name = expression("Atmos O"[2]* " equivalent (%)"))) + 
geom_line(data = summary2, aes(x=Secs, y=mean/coeff), color='Red')+ geom_vline(data= mean_LOP, aes (xintercept = z)) + geom_text(aes(x=80.2, label="CO2mean LOP\n", y=350), colour="black", angle=90, text=element_text(size=16)) +
  geom_text(aes(x=273.7, label="Decom mean LOP\n", y=600), colour="black", angle=90, text=element_text(size=16))
p4

ggsave("Terminal Treatment profiles O2 conc.tiff", units="mm", width=250, height=150, dpi=600, compression = 'lzw')
```


```{r}
# grid these figures together 

figure <- ggarrange(p, p3,p2, ncol = 3, nrow = 2, common.legend = T, legend="bottom", labels = c( "A.", "B.", "C."))
figure
ggsave("Terminal Treatment profiles Manuscript.tiff", units="mm", width=190, height=250, dpi=600, compression = 'lzw')
```


```{r}
figure <- ggarrange(p, p5, ncol = 2, nrow = 1, common.legend = T, legend="none", labels = c( "A.", "B."))
figure
ggsave("Terminal Treatments info.tiff", units="mm", width=190, height=100, dpi=600, compression = 'lzw')
```





```{r}
# looking at the pressure profiles of the mice that ear scratched 
ears <- subset(Pressure, MouseID== c("4.1","7.3",'9.2','14.3', '3.2', '6.1', '20.1', Batch='1')) 
ears2 <- subset(Pressure, MouseID==c("7.2","12.2",'9.2','19.1', '2.1', '4.2', '7.1', '12.1', '16.1', '20.3', '2.3', '20.1', Batch ='2'))

ears$MouseID <- as.factor(ears$MouseID)
ears2$MouseID <- as.factor(ears2$MouseID)

ggplot(data=ears, aes(x=Secs, y=Pressure, color=MouseID)) +
geom_point()+
  labs(title = "Decompression profiles for mice ear scratching in Batch 1", x="Time (s)",y="Pressure (mBar)") + ggeasy::easy_center_title()

ggplot(data=ears2, aes(x=Secs, y=Pressure, color=MouseID)) +
geom_point()+
  labs(title = "Decompression profiles for mice ear scratching in Batch 2", x="Time (s)",y="Pressure (mBar)") + ggeasy::easy_center_title()
```


## additional preliminary models 
to check that time to reach 200mBar and time to increase CO2 did not affect outcomes 
This requires the data to be subsetted into decompression and carbon dioxide only.


```{r}
decom <- subset(Exp2, Treatment %in% c("Decom"))
Co2 <- subset(Exp2, Treatment %in% c("CO2"))
```

run the models for each behaviour - include the time to 200 and co2 time increase as covariate to determine whether significant 

#proportions models
```{r}
logit_1 <- glm(EarScratchBinary ~ Timeto200 + Intervention, family=binomial(), data = decom)
summary(logit_1)
Anova(logit_1)

logit_2 <- glm(HeadFlickBin ~ Timeto200 + Intervention, family=binomial(), data = decom)
summary(logit_2)
Anova(logit_2)

logit_3 <- glm(TotalGrimScore ~ Timeto200 + Intervention, family=binomial(), data = decom)
summary(logit_3)
Anova(logit_3)
```

## latency models - decom 

```{r}
p1 <-lmer(EarScratchLat ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p1)
Anova(p1)

p2 <-lmer(LOBLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p2)
Anova(p2)

p3 <-lmer(GaspingLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p3)
Anova(p3)

p4 <-lmer(AgonalLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p4)
Anova(p4)

p5 <-lmer(ApnoeaLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p5)
Anova(p5)

p6 <-lmer(ConvlLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p6)
Anova(p6)

p7 <-lmer(FacialGrimaceLat ~ Timeto200 + Intervention + (1|Batch), data=decom) 
summary(p7)
Anova(p7)

p8 <-lmer(AtaxiaLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p8)
Anova(p8)

p9 <-lmer(LOPLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p9)
Anova(p9)

p10 <-lmer(RecumbentLat ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p10)
Anova(p10)

p11 <-lmer(MotionlessLatency ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p11)
Anova(p11)

p12 <-lmer(AbnormalCombined ~ Timeto200 + Intervention + (1|Weight) + (1|Batch/Cage), data=decom) 
summary(p12)
Anova(p12)
```

# look at the pairwise effects between interventions 

```{r}

emm2<- emmeans(p2, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")

emm2<- emmeans(p4, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")

emm2<- emmeans(p5, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")

emm2<- emmeans(p9, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")

emm2<- emmeans(p4, list(pairwise ~ Intervention), adjust = "tukey")
summary(emm2, type = "response")
```

## Carbon Dioxide - models to look at CoV of time of gas increase


```{r}


p2 <-lmer(LOBLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p2)
Anova(p2)

p3 <-lmer(GaspingLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p3)
Anova(p3)

p4 <-lmer(AgonalLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p4)
Anova(p4)

p5 <-lmer(ApnoeaLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p5)
Anova(p5)



p7 <-lmer(FacialGrimaceLat ~ CO2IncS + Intervention + (1|Batch), data=Co2) 
summary(p7)
Anova(p7)

p8 <-lmer(AtaxiaLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p8)
Anova(p8)

p9 <-lmer(LOPLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p9)
Anova(p9)


p11 <-lmer(MotionlessLatency ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p11)
Anova(p11)

p12 <-lmer(AbnormalCombined ~ CO2IncS + Intervention + (1|Weight) + (1|Batch/Cage), data=Co2) 
summary(p12)
Anova(p12)
```


## Figures for data summary


### proportions data - use probability data or actual numbers?
# ear scratch 

```{r}
EarPro <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsEarScratchProb.xlsx")
EarPro$facets = factor(EarPro$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

z1 <- ggplot(EarPro, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) + 
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))+ geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9))+ ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Probability") + ylim(0,1) + facet_wrap(~facets, labeller = label_parsed) +  ggtitle("Ear Scratching") + theme_bw() + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3)
z2 <- z1 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
z2
```

##  Head Flick 

```{r}
HFPro <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsHeadFlickProb.xlsx")
HFPro$facets = factor(HFPro$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

z3 <- ggplot(HFPro, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9))+  ggtitle("Head Flick") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Probability") + ylim(0,1.2) + facet_wrap(~facets, labeller = label_parsed) + theme_bw() + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3)

z4 <- z3 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
z4
```



## facial grimace - changes in ear posture 
```{r}
GrimPro <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsGrimProb.xlsx")
GrimPro$facets = factor(GrimPro$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

z5 <- ggplot(GrimPro, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9))+  ggtitle("Changes in Ear Posture") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Probability") +  ylim(-0.1,1) + facet_wrap(~facets, labeller = label_parsed) + theme_bw() + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3)

z6 <- z5 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 12))
z6
```


```{r}
# grid these figures together - where all three treatment groups are present

figure <- ggarrange(z2, z4, z6,  ncol = 2, nrow = 2, common.legend = T, legend="none", labels = c( "A.", "B.", "C."))
figure

ggsave("ProportionData MS.tiff", units="mm", width=250, height=200, dpi=600, compression = 'lzw')
```


## Duration data - Figures

### Ear Scratching Duration 

```{r}
EarDur <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsEarScratchDur.xlsx")
EarDur$facets = factor(EarDur$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

f1 <- ggplot(EarDur, aes(x=Intervention, y=mean, fill = (Treatment))) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))  +
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + 
 ylab("Latency (s)") + ggtitle("Total Duration ear scratching")+labs(x="Pharmacological Intervention",y=" duration spent ear scratching (s)") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() + geom_text(aes(label = Freq), vjust = 0) 
f2 <- f1 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
f2
```


## ear scratch as a proprotion of conscious time - deleted this sheet- need to create with new emmeans output
```{r}
EarDurProp <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsEarScartchDurAsProportion.xlsx")
EarDurProp$facets = factor(EarDurProp$Treatment, labels = c("CO[2]","Decompression", 'Sham'))

e1 <- ggplot(EarDurProp, aes(x=Intervention, y=mean, fill = (Treatment))) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9))  +
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + 
 ylab("Latency (s)") + ggtitle("Proportion of conscious time ear scratching")+ labs(x="Pharmacological Intervention",y=" proportion of conscious time") + facet_wrap(~facets, labeller = label_parsed) + theme_bw()
e2 <- e1 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
e2
```


## Grooming bout duration

```{r}
GroomB <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsGroomBDur.xlsx")
GroomB$facets = factor(GroomB$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f3 <- ggplot(GroomB, aes(x=Intervention, y=mean, fill = (Treatment),  ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Duration spent grooming") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Duration (s)") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f4 <- f3 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f4
```


## exploration of the periphery 

```{r}
Expl <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsExplDur.xlsx")
Expl$facets = factor(Expl$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f5 <- ggplot(Expl, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Duration of Exploration") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Duration (s)") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f6<- f5 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f6

```

# duration of head up posture 


```{r}
HeadUp <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsHeadUplDur.xlsx")
HeadUp$facets = factor(HeadUp$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f7 <- ggplot(HeadUp, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,position = position_dodge(0.9)) +
  ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Duration (s)") + ggtitle("Duration of head up") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f8<- f7 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 11))
f8
```

# frequency of head up 

```{r}
HeadUp <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsHeadUplFreq.xlsx")
HeadUp$facets = factor(HeadUp$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))


h1 <- ggplot(HeadUp, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of head up postures") + ylab("Number of events") + labs(x="Pharmacological Intervention",y="Number of events") +facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
h2<- h1 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 11))
h2
```

```{r}
# grid these figures together - head up figure manuscript

figure <- ggarrange(h2, f8,  ncol = 2, nrow = 1, common.legend = T, legend="none", labels = c( "A.", "B."))
figure

ggsave("Head up MS.tiff", units="mm", width=250, height=110, dpi=600, compression = 'lzw')
```


## duration of digging 

```{r}
Dig <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsDiggingDur.xlsx")
Dig$facets = factor(Dig$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))


f9 <- ggplot(Dig, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Duration of digging") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Duration (s)") +facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f10<- f9 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f10


```




## convulsion duration (decompression only)


```{r}
Convl <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsConvlDur.xlsx")
Convl$facets = factor(Convl$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f11 <- ggplot(Convl, aes(x=Intervention, y=mean, fill = (Treatment))) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Duration of convulsions (s)") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f12<- f11 + theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))+ theme(text = element_text(size = 14))
f12
```

## grid them together to make a multipanel 


```{r}
# grid these figures together - where all three treatment groups are present

figure <- ggarrange(f2, f8, f12, f4, f6, f10,  ncol = 2, nrow = 3, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D.", "E.", "F."))
figure

ggsave("Duration Data MS.tiff", units="mm", width=200, height=250, dpi=600, compression = 'lzw')
```


## count/frequency data 

# number of ear scratches

only plotted for animals that performed the behaviour - which is what the annotated numbers refer to

```{r}
EarFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsEarscratchFreq.xlsx")
EarFrq$facets = factor(EarFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f13 <- ggplot(EarFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of ear scratches") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Number of events") + ylim(0,8)+facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f14<- f13 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 14))
f14
```

## number of head flicks
this was not modelled on the emms from the analysis looking at only those mice that performed the behaviour instead this is modelled on mean out of all mice - this was due to lack of ability to calculate means. 

```{r}
HFFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsHFFreq_allmice.xlsx")
HFFrq$facets = factor(HFFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f15 <- ggplot(HFFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of Head Flicks") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Number of events") + ylim(0,6) + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f16<- f15 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 14))
f16
```

## number of head flicks 
this is calculated from the means from model focused only on the animals that performed the behaviour - less susceptible to zero inflation

```{r}
HFFrq2 <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsHFFreq_subset.xlsx")
HFFrq2$facets = factor(HFFrq2$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f17 <- ggplot(HFFrq2, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of head flicks")+ ylab("Latency (s)") + labs(x="Pharmacological Intervention",y=" Number of events") + ylim(0,6) + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f18<- f17 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 14))
f18
```


## grooming bout frequency 
like the above two behaviours - plotted from emms from the model excluding the animals that didnt perform the behaviour


```{r}
GBFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsGroomFrq.xlsx")
GBFrq$facets = factor(GBFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f19 <- ggplot(GBFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of grooming bouts") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Number of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f20<- f19 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f20
```

## number of digging events 

```{r}
DigFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsDigFrq.xlsx")
DigFrq$facets = factor(DigFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f27 <- ggplot(DigFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of digging events") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Number of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f28<- f27 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f28
```


```{r}
nrow(s2d[s2d$Treatment == "CO2" & s2d$Intervention == "Saline",])
nrow(s2d[s2d$Treatment == "Decom" & s2d$Intervention == "Saline",])
nrow(s2d[s2d$Treatment == "Sham" & s2d$Intervention == "Saline",])
nrow(s2d[s2d$Treatment == "CO2" & s2d$Intervention == "Analgesic",])
nrow(s2d[s2d$Treatment == "Decom" & s2d$Intervention == "Analgesic",])
nrow(s2d[s2d$Treatment == "Sham" & s2d$Intervention == "Analgesic",])
nrow(s2d[s2d$Treatment == "CO2" & s2d$Intervention == "Anxio",])
nrow(s2d[s2d$Treatment == "Decom" & s2d$Intervention == "Anxio",])
nrow(s2d[s2d$Treatment == "Sham" & s2d$Intervention == "Anxio",])

```

## number of gasps 

```{r}
GaspFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsGaspFrq.xlsx")
GaspFrq$facets = factor(GaspFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f21 <- ggplot(GaspFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean, ymax=mean+SE), width=.2, 
                position = position_dodge(0.9)) + 
 ggtitle("Number of gasps") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Number of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f22<- f21 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 14)) +ylim(0,120)
f22
```

# number of agonal gasps 

```{r}
AgGaspFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsAgonalGaspFrq.xlsx")
AgGaspFrq$facets = factor(AgGaspFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f23 <- ggplot(AgGaspFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Numbeer of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f24<- f23 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 14))
f24
```


## rearing 

```{r}
RearFrq <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\EMMsRearFrq.xlsx")
RearFrq$facets = factor(RearFrq$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

f25 <- ggplot(RearFrq, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of rears") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Number of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
f26<- f25 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + geom_text(aes (y = stat(ymax), label = n), vjust = -0.5, size =3) + theme(text = element_text(size = 12))
f26
```

# grid together 

```{r}
# grid these figures together - where all three treatment groups are present

figure <- ggarrange(f14, f18, f20, f22, f24, f26,  ncol = 2, nrow = 3, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D.", "E.", "F."))
figure

ggsave("Frequency Data MS.tiff", units="mm", width=250, height=300, dpi=600, compression = 'lzw')
```


## figures for talk - not able to include all data 
```{r}
# grid these figures together - for talk - abnormal behaviours 

figure <- ggarrange(f14, f18,  ncol = 2, nrow = 1, common.legend = T, legend="none", labels = c( "A.", "B.", "C."))
figure

ggsave("Frequency Data Talk Abnormal.tiff", units="mm", width=250, height=100, dpi=600, compression = 'lzw')
```

```{r}
# grid these figures together - for talk - normal behaviours 

figure <- ggarrange(f26, f20, f10, f6, f4,  ncol = 3, nrow = 2, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D.", "E."))
figure

ggsave("Frequency & Dur Data Talk Normal.tiff", units="mm", width=300, height=200, dpi=600, compression = 'lzw')

```

```{r}
# grid these figures together - for manuscript - illustrating normal behaviours 

figure <- ggarrange(f26, f20, f28, f10, f4, f6,  ncol = 2, nrow = 3, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D.", "E.", "F."))
figure

ggsave("Manuscript normal behaviour.tiff", units="mm", width=250, height=330, dpi=600, compression = 'lzw')

```
```{r}
# grid these figures together - for poster 

figure <- ggarrange(f20, f4, f10, f6,  ncol = 2, nrow = 2, common.legend = T, legend="none")
figure

ggsave("Poster2.tiff", units="mm", width=300, height=200, dpi=600, compression = 'lzw')

```

## plot escape attempts 


```{r}
escape <-Exp2[, c(4,5,25)] # selects only those columns to plot 
escape$facets = factor(escape$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

g7 <- ggplot(escape, aes(x=Intervention, y=EscapeAtt, fill = (Treatment))) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 1) + ggtitle("Number of escape attempts") + labs(x="Pharmacological Intervention",y=" Frequency") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 

g8<- g7 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF")) + theme(text = element_text(size = 14))
g8

escape2 <- read_excel("J:\\BAHCM\\LAPSProject\\Experiment 2 Mice\\Analysis\\MeansEscapeAttempt.xlsx")
escape2$facets = factor(escape2$Treatment, labels = c("CO[2]" , "Decompression", 'Sham'))

q2 <- ggplot(escape2, aes(x=Intervention, y=mean, fill = (Treatment), ymin =mean-SE, ymax = mean+SE)) +
   geom_bar(stat ="identity", colour = "black", position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width=.2,
                position = position_dodge(0.9)) + ggtitle("Number of escape attempts") + ylab("Latency (s)") + labs(x="Pharmacological Intervention",y="Number of events") + facet_wrap(~facets, labeller = label_parsed) + theme_bw() 
q3<- q2 + theme(legend.position="none", axis.title.x = element_blank(),  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + geom_text(aes (label = n), vjust = -0.5, size =3)+ scale_fill_manual(values=c("#481567FF","#29AF7FFF", "#DCE319FF"))  + theme(text = element_text(size = 14)) + ylim(0,4)
q3


```

```{r}
# grid escape and gasps together for talk 

figure <- ggarrange(f22, g8,  ncol = 2, nrow = 1, common.legend = T, legend="none", labels = c( "A.", "B."))
figure

ggsave("Gasp and Escape Talk.tiff", units="mm", width=250, height=100, dpi=600, compression = 'lzw')

```


```{r}
# grid together for manuscript

figure <- ggarrange(f14, f18, f22, q3,  ncol = 2, nrow = 2, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D."))
figure

ggsave("Frequency data.tiff", units="mm", width=250, height=200, dpi=600, compression = 'lzw')
```

```{r}
# grid together for poster

figure <- ggarrange(f22, f14,  ncol = 1, nrow = 2, common.legend = T, legend="none")
figure

ggsave("Poster1.tiff", units="mm", width=150, height=250, dpi=600, compression = 'lzw')
```

--------
none of this is included in manuscript - just additional exploration of the data 
## ear scartch duration - absolute vs proportion of conscious time - not included /correct at present

```{r}

figure <- ggarrange(f2, e2,  ncol = 2, nrow = 1, common.legend = T, legend="none", labels = c( "A.", "B.", "C.", "D."))
figure

ggsave("Ear Scratch Dur data Manuscript.tiff", units="mm", width=250, height=100, dpi=600, compression = 'lzw')
```





## duration data - calculate the proportion of CONSCIOUS time (to LOP) spent doing these behaviours 
this should then allow us to make directly comparable 
 
 
# exploration of the periphery - as a proportion of conscious time 
```{r}
a1 <-lmer(ExplProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(a1)

simulationOutput <- simulateResiduals(fittedModel = a1, plot = T)
```

```{r}
Anova(a1)
emm1 <- emmeans(a1, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a1, list(pairwise ~ Intervention), adjust ="tukey")
emm3 <- emmeans(a1, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```
 
# dig as a proportion of LOP time 

```{r}
a2 <-lmer(DigProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(a2)

simulationOutput <- simulateResiduals(fittedModel = a2, plot = T)

```

```{r}
Anova(a2)
emm1 <- emmeans(a2, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a2, list(pairwise ~ Intervention), adjust ="tukey")
emm3 <- emmeans(a2, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# ears scratch as a proportion of LOP time 
this is important finding as - absolute duration are higher with decom compared to sham but not when you account for the available time to do so!! 

```{r}
a3 <-lmer(EarSProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2)
summary(a3)

simulationOutput <- simulateResiduals(fittedModel = a3, plot = T)
```

```{r}
Anova(a3)
emm1 <- emmeans(a3, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a3, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(a3, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# head up proportion of LOP time 

```{r}
a4 <-lmer(HeadUpProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(a4)

simulationOutput <- simulateResiduals(fittedModel = a4, plot = T)
```


```{r}
Anova(a4)
emm1 <- emmeans(a4, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a4, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(a4, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# grooming as a proportion of LOP time 

```{r}
a5 <-lmer(GroomProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(a5)

simulationOutput <- simulateResiduals(fittedModel = a5, plot = T)
```

```{r}
Anova(a5)
emm1 <- emmeans(a5, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a5, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(a5, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

# Inactive as a proportion of LOP time 

```{r}
a6 <-lmer(InactiveProportion ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage), data=Exp2) 
summary(a6)

simulationOutput <- simulateResiduals(fittedModel = a6, plot = T)

```

```{r}
Anova(a6)
emm1 <- emmeans(a6, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(a6, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(a6, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```

```{r}
summary <- Exp2 %>% count(Treatment, Intervention, GaspingLatency, sort = TRUE)
```


# ataxic duration 
this was calculated manually by subtracting latency to ataxic from latency to LOP 


```{r}
m0.1 <-lmer(AtaxicDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage),  data=Exp2) 
summary(m0.1)
```

```{r}
simulationOutputm0.1<- simulateResiduals(fittedModel = m0.1, plot = T)
plot(simulationOutputm0.1)
```


```{r}
Anova(m0.1)
emm1 <- emmeans(m0.1, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(m0.1, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(m0.1, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```


#abnormal breathing duration
```{r}
m0.2 <-glmmTMB(AbnormalDur ~ Treatment + Intervention + Treatment:Intervention + (1|Weight) + (1|Batch/Cage),  data=Exp2, family=nbinom2()) 
summary(m0.2)
```





```{r}
simulationOutputm0.2<- simulateResiduals(fittedModel = m0.2, plot = T)
plot(simulationOutputm0.2)
```


```{r}
Anova(m0.2)
emm1 <- emmeans(m0.2, list(pairwise ~ Treatment), adjust ="tukey")
emm2 <- emmeans(m0.2, list(pairwise ~ Treatment), adjust ="tukey")
emm3 <- emmeans(m0.2, list(pairwise ~ Intervention|Treatment), adjust ="tukey")
summary(emm1, type = "response")
summary(emm2, type = "response")
summary(emm3, type = "response")
```