---
title: "Experiment 1 LAPS mice preliminary analysis"
Author: "Jasmine Clarkson"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r include = FALSE}
#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)
```


```{r}
#load in the dataset
library(readr)
BehaviourExperiment1Mice <- read_csv("J:\\BAHCM\\LAPSProject\\Exp1MiceAnalysis\\BehaviourExperiment1Mice.csv",col_types = cols(Date = col_date(format = "%d/%m/%Y")))

# allows you to view the data
View(BehaviourExperiment1Mice)
```



```{r}
# Label the dataframe/factors
BehaviourExperiment1Mice$Date <- as.factor(BehaviourExperiment1Mice$Date)
BehaviourExperiment1Mice$Run <- as.factor(BehaviourExperiment1Mice$Run)
BehaviourExperiment1Mice$AnimalID <- as.factor(BehaviourExperiment1Mice$AnimalID)
BehaviourExperiment1Mice$Cage <- as.factor(BehaviourExperiment1Mice$Cage)
BehaviourExperiment1Mice$CurveType <- as.factor(BehaviourExperiment1Mice$CurveType)
BehaviourExperiment1Mice$Rate <- as.factor(BehaviourExperiment1Mice$Rate)
BehaviourExperiment1Mice$Shape <- as.factor(BehaviourExperiment1Mice$Shape)
BehaviourExperiment1Mice$Strain <- as.factor(BehaviourExperiment1Mice$Strain)
BehaviourExperiment1Mice$Batch <- as.factor(BehaviourExperiment1Mice$Batch)
BehaviourExperiment1Mice$Rate <- as.factor(BehaviourExperiment1Mice$Rate)
BehaviourExperiment1Mice$Shape <- as.factor(BehaviourExperiment1Mice$Shape)
```

## Latency to first gasp data 

This section is looking at the model that best fits the data so I start with fitting a linear mixed model. 

### Model 1 (Latency to First Gasp)

```{r}
# Fit a linear mixed model with weight as covariate, rate, shape and strain as fixed factors and cage nested within batch as random effects
m1 <-lmer(LatencyFirstGasp_s ~ Weight + Rate + Shape + Strain + Rate:Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(m1)
```
```{r}
# looking at the fit of the residuals model - not a great fit 
simulationOutputm1<- simulateResiduals(fittedModel = m1, plot = T)
plot(simulationOutputm1)
```

```{r}
# looking at the distibution 
hist(simulationOutputm1)
```

### Model 2 - (Latency to First Gasp)

```{r}
#  Negative binominal model - included weight as covariate but due to not having an effect and better fit of model when it is included as a random effect. Rate, shape and strain as fixed factors but breaks down the three way interaction into more meaningful comparisons. The model struggles for output if the three way interaction is included due to a lack of df. 
m2 <-glmmTMB(LatencyFirstGasp_s ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m2)
```
```{r}
# next we want to look at the fit of the model - this looks like a much better fit 
simulationOutputm2<- simulateResiduals(fittedModel = m2, plot = T)
plot(simulationOutputm2)

```


```{r}
# looking at the distibution 
hist(simulationOutputm2)
```
```{r}
# I need to also look at whether  the temperature and the humidity has any effect- put in as covariates in the model - in order for R to consider the variable a covariate you need to make sure it isn't labelled as a factor and also that it is a numerical variable. 
m2F <-lmer(LatencyFirstGasp_s ~ Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(m2F)
# this shows that when we include temperature and RH as covariates they have a non-sig effect so okay to not include these in the model. 
```
```{r}
# this shows that including these covariates makes the fit of the model much worse 
simulationOutputm2F<- simulateResiduals(fittedModel = m2F, plot = T)
plot(simulationOutputm2F)
```



## Model 2 output 
We see a main effect of rate, an interaction between rate and shape and rate and strain.

```{r}
# model output 
anovam2 <- Anova(m2)
anovam2
```
# Pairwise comparisons 
Enables us to see where the differences lie for our significant interactions.

```{r}
# looking at the effect of rate
# significantly longer latency to first gasp with 75m/s rate compared to 150m/s (irrespective of shape and strain)

emm1 <- emmeans(m2, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1, type = "response")
ss<- summary(emm1, type = "response")
ss
write.table(as.data.frame(ss),file="StatsSummFirstGaspRate.csv", quote=F,sep=",",row.names=F)
```
```{r}
# looking at the interaction between rate and shape with shape nested in rate - this shows that there is a difference in the latency to first gasp between the two rates only for the gradual shape but not for accelerated or linear. 
emm2 <- emmeans(m2, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm2, type = "response")
```
```{r}
# Looking at the effect of shape within rate 
emm <- emmeans(m2, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm, type = "response")

```

```{r}
# this looks at the interaction between rate and strain (significant in the main model) remember you can look at this interaction in both ways. 
# emm3 looks at the difference between the two rates across each strain independently. You can see the reason why we have a significant strain by rate interaction is because only c57s show a longer latency to first gasp with the slower 75m/s rate with 150m/s but the balbs do not. 
emm3 <- emmeans(m2, list(pairwise ~ Rate|Strain), adjust = "tukey")
summary(emm3, type ="response")

```
```{r}
# this looks at the interaction between rate and strain looking at the effect of strain across the two rates independently. i.e. is there a difference between our strains for 75m/s and 150m/s - can see here that strain doesn't effect responses at each rate independently. This comparison is most important for assessing our research question that within rate our two strains do not differ in their response. 
emm4 <- emmeans(m2, list(pairwise ~ Strain|Rate), adjust = "tukey")
summary(emm4, type = "response")
```

# Plots for latency to first gasp 
Here we graph the data to look at the trends we see. 

```{r}
library(ggplot2)
# subset the data for the latency to first gasp within the data frame i.e. 5th, 8th and 16th column
FirstGaspDataAll_strain_curvetype <- BehaviourExperiment1Mice[, c(5,6,8,16)]
FirstGaspDataAll_strain_curvetype$CurveType <- factor(FirstGaspDataAll_strain_curvetype$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
FirstGaspDataAll_strain_curvetype$Rate <- factor(FirstGaspDataAll_strain_curvetype$Rate, levels = c("75", "150"))
#plot the subsetted data                                                               
p3 <- ggplot(FirstGaspDataAll_strain_curvetype, aes(x=CurveType, y=LatencyFirstGasp_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to first gasp (s) ")

p3
library(viridis)
g3 <- p3 + scale_fill_viridis(option = "viridis", discrete = TRUE)
g3
```
```{r}
# produce a plot looking at the two strains differently at each curv to visualise the difference between strains.You can see that they show the same trends i.e. longer latency to first gasp with the slower rate.

p <- ggplot(FirstGaspDataAll_strain_curvetype, aes(x=CurveType, y=LatencyFirstGasp_s, fill=Strain)) +
  geom_boxplot(position=position_dodge(1)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
                 position=position_dodge(1),dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to first gasp (s) ")
p
```
## Latency to agonal gasp data

### Model 3 (Latency to Agonal Gasp)

Start with fitting a linear mixed model, included weight as a covraiate but due to non-significant effect instead included it as a random effect. 

```{r} 
# first included weight as covariate but non-significant effect but still want to account for it in the model so moved a random effect 
m3 <-lmer(LatencyAgonalGasp_s ~ Rate + Shape + Strain + Rate:Shape:Strain + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(m3)
AIC(m3)
```

```{r}
# looking at the fit of the residuals model - not a bad fit but not perfect
simulationOutputm3<- simulateResiduals(fittedModel = m3, plot = T)
plot(simulationOutputm3)

```
```{r}
# looking at the distibution 
hist(simulationOutputm3)
```
### Model 4 (LAtency to Agonal gasp) 
Instead try to fit a negative binomial model to see if this improves the fit of the model.

```{r}
m4 <-glmmTMB(LatencyAgonalGasp_s ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m4)
```
```{r}
# looking at the fit of the residuals model - seems like a better fit 
simulationOutputm4<- simulateResiduals(fittedModel = m4, plot = T)
plot(simulationOutputm4)
```
```{r}
# looking at the distibution = model 4 fits better in terms of qq plot and histo
hist(simulationOutputm4)
```
```{r}
# want to again check that temp and RH arent having an effect - include these as covariates and check that they dont explain some of the variance. 
m4F <-glmmTMB(LatencyAgonalGasp_s ~ Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m4F)
# this shows that no effect of temp and RH so again okay to drop these from the model. 
```




### Model 4 output 
```{r}
# model output - shows a main effect of rate but no other interactions. 
anovam4 <- Anova(m4)
anovam4
```

# Pairwise comparisons 

```{r}
# looking at the effect of rate
# significantly longer latency to agonal gasp with 75m/s rate compared to 150m/s (irrespective of shape and strain)

emm1 <- emmeans(m4, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1, type = "response")
```
```{r}
# looking at the interaction between rate and shape with rate nested in shape - this shows that there is a difference in the latency to agonal gasp between the two rates only for the gradual and linear shape (near signif) but not for accelerated. 
emm2 <- emmeans(m4, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm2, type = "response")
```
```{r}
# Looking at the effect of shape within rate 
emm <- emmeans(m4, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm, type = "response")
```




```{r}
# this looks at the interaction between rate and strain (non significant in the main model) remember you can look at this interaction in both ways. 
# emm3 looks at the difference between the two rates across each strain independently - shows that both strains have longer latencies overall between 75 and 150 and emphasises the effect of rate irrespective of strain
emm3 <- emmeans(m4, list(pairwise ~ Rate|Strain), adjust = "tukey")
summary(emm3, type ="response")
```
```{r}
#this looks at the interaction the opposite way, looks at the difference between the two strains for each rate independently. Demonstrates we have no strain effect when comparing both strains for each rate. 
emm4 <- emmeans(m4, list(pairwise ~ Strain|Rate), adjust = "tukey")
summary(emm4, type ="response")
```
```{r}
# out of interest wanted to look at the shape rate interaction the other way but this isn't informative as there is no effect of shape overall so comparing each shape with one another isnt important as you can see from the graphs below. 
emm5 <- emmeans(m4, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm5, type = "response")
```


# Plots for latency to agonal gasp

```{r}
#first subsetting the data for the appropriate columns 
LatencyAgonalDataAll_strain_curvetype <- BehaviourExperiment1Mice[, c(5,6,8,24)]
LatencyAgonalDataAll_strain_curvetype$CurveType <- factor(LatencyAgonalDataAll_strain_curvetype$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
LatencyAgonalDataAll_strain_curvetype$Rate <- factor(LatencyAgonalDataAll_strain_curvetype$Rate, levels = c("75", "150"))
#plot the graph
p3.5 <- ggplot(LatencyAgonalDataAll_strain_curvetype, aes(x=CurveType, y=LatencyAgonalGasp_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to agonal gasps (s) ")

p3.5
```


```{r}
p.5 <- ggplot(LatencyAgonalDataAll_strain_curvetype, aes(x=CurveType, y=LatencyAgonalGasp_s, fill=Strain)) +
  geom_boxplot(position=position_dodge(1)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center',
               position=position_dodge(1),dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to agonal gasps (s) ")
p.5
```


## Latency to motionless i.e. our measure of death from the behavioural observations

### Model 5 (Latency to motionless)

```{r}
# first included weight as covariate but significant effect therefore need to keep this in as a covariate; also tried this by including a three way interaction and didn't fit any better.  
m5 <-lmer(Motionless_s ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(m5)
AIC(m5)
```


```{r}
# this looks okay but not great 
simulationOutputm5<- simulateResiduals(fittedModel = m5, plot = T)
plot(simulationOutputm5)
```
```{r}
# looking at the distribution 
hist(simulationOutputm5)

```
## Model 6 (Latency to Motionless)
Again try and fit a negative binomial model and see if this fits any better 

```{r}
m6 <-glmmTMB(Motionless_s ~ Weight +  Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m6)
```

```{r}
#this looks like a better fit than model 5 in terms of overdispersion 
simulationOutputm6<- simulateResiduals(fittedModel = m6, plot = T)
plot(simulationOutputm6)
```
```{r}
# looking at the distribution - looks better than model 5 
hist(simulationOutputm6)
```
```{r}
# again just want to check that temp and RH dont have an effect when included as covariates in the model. Temperature looks to have a significant effect here but no effect of RH. 
m6F <-glmmTMB(Motionless_s ~ Temperature + RelativeHumidity + Weight +  Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m6F)

```
```{r}
# not a great fit when including these covariates - passable but not great
simulationOutputm6F<- simulateResiduals(fittedModel = m6F, plot = T)
plot(simulationOutputm6F)
```
## model 6F output
Temperature and RH included in the model 
```{r}
anovam6F <- Anova(m6F)
anovam6F

```

## Model 6 output 
(temperature and Rh not included in the model)

```{r}
anovam6 <- Anova(m6)
anovam6
```
# Pairwise comparisons 


```{r}
# looking at the effect on rate in latency to motionless ; 75m/s takes longer to motionless (irrespective of strain and shape)
emm1 <- emmeans(m6, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1, type = "response")
```
```{r}
# looking at the effect on rate in latency to motionless ; 75m/s takes longer to motionless (irrespective of strain and shape)
emm1F <- emmeans(m6F, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1F, type = "response")
```


```{r}
# looking at the differences between shapes - as marginal result from the model. No differences here but trend towards the gradual taking longer, followed by the linear and then the accelerated shapes. 

emm2 <- emmeans(m6, list(pairwise ~ Shape), adjust = "tukey")
summary(emm2, type = "response")
```
```{r}
# looking at the differences between shapes - as marginal result from the model. No differences here but trend towards the gradual taking longer, followed by the linear and then the accelerated shapes. 

emm2F <- emmeans(m6F, list(pairwise ~ Shape), adjust = "tukey")
summary(emm2F, type = "response")
```


```{r}
# in line with what we see with the other measures. Animals are faster to reach motionless state for the 150ms rate compared to 75ms for only the gradual and linear curves but not the accelerated shape
emm3 <- emmeans(m6, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm3, type = "response")
```
```{r}
# in line with what we see with the other measures. Animals are faster to reach motionless state for the 150ms rate compared to 75ms for only the gradual and linear curves but not the accelerated shape
emm3F <- emmeans(m6F, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm3F, type = "response")
```


```{r}
# Looking at the effect of shape within rate 
emm <- emmeans(m6, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm, type = "response")
```
```{r}
# Looking at the effect of shape within rate 
emmF <- emmeans(m6F, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emmF, type = "response")
```



```{r}
# just to confirm that our strains didn't behave differently at our two rates. (also checked shape and strain wasn't significant either)
emm4 <- emmeans(m6, list(pairwise ~ Strain|Rate), adjust = "tukey")
summary(emm4, type = "response")
```
```{r}
# just to confirm that our strains didn't behave differently at our two rates. (also checked shape and strain wasn't significant either)
emm4F <- emmeans(m6F, list(pairwise ~ Strain|Rate), adjust = "tukey")
summary(emm4F, type = "response")
```




```{r}
# just to confirm that our strains didn't behave differently according to our two rates when compared directly. both strains have longer latencies with the 75ms rate. 
emm4 <- emmeans(m6, list(pairwise ~ Rate|Strain), adjust = "tukey")
summary(emm4, type = "response")
```
```{r}
# just to confirm that our strains didn't behave differently according to our two rates when compared directly. both strains have longer latencies with the 75ms rate. 
emm4F <- emmeans(m6F, list(pairwise ~ Rate|Strain), adjust = "tukey")
summary(emm4F, type = "response")
```


# Latency to motionless plots

```{r}
MotionlessDataAll_strain_curvetype <- BehaviourExperiment1Mice[, c(5,6, 8,29)]
MotionlessDataAll_strain_curvetype$CurveType <- factor(MotionlessDataAll_strain_curvetype$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
MotionlessDataAll_strain_curvetype$Rate <- factor(MotionlessDataAll_strain_curvetype$Rate, levels = c("75", "150"))
# Coloured box and dots for strain by curve type
p.8 <- ggplot(MotionlessDataAll_strain_curvetype, aes(x=CurveType, y=Motionless_s, fill=Strain)) +
  geom_boxplot(position=position_dodge(1)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center',
               position=position_dodge(1),dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to motionless (s) ")
p.8
```
```{r}
p3.8 <- ggplot(MotionlessDataAll_strain_curvetype, aes(x=CurveType, y=Motionless_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to motionless (s) ")

p3.8
```
## Total Number of gasps 
This is not latency data but instead is count data therefore need to first fit a poisson model and see how it fits. We also need to consider the fact that we need to include an offset because we have different cycle lengths between our two rates and therefore need to adjust for that in the model. Total gasps excludes agonal gasps - but instead includes the sum of sporadic and rhythmic gasps.

### Model 7 (total gasps)
```{r}
m7 <- glmmTMB(TotalGasps ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = poisson) 
summary(m7)
```
```{r}
# this model does not fit well 
testDispersion(m7)
simulationOutputm7 <- simulateResiduals(fittedModel = m7, plot = T)
```

```{r}
hist(simulationOutputm7)
```



### model 8 (total gasps)
this model using nbinom 1 = quasi poisson distribution 

```{r}
m8 <- glmmTMB(TotalGasps ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom1) 
summary(m8)
```
```{r}
# better fit byt still not perfect 
testDispersion(m8)
simulationOutputm8 <- simulateResiduals(fittedModel = m8, plot = T)
```

```{r}
hist(simulationOutputm8)
```

## Model 9 (Total gasps)
This model uses the nbinom 2 - Negative Binomial distribution 

```{r}
m9 <- glmmTMB(TotalGasps ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m9)
```
```{r}
# this model fits much better 
testDispersion(m9)
simulationOutputm9 <- simulateResiduals(fittedModel = m9, plot = T)
```

```{r}
hist(simulationOutputm9)
```
```{r}
# included temp and RH in the model as covariates to check that these font have an effect - they dont so okay to drop them as above. 
m9F <- glmmTMB(TotalGasps ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(m9F)
```


## model output from model 9 (total gasps)

```{r}
anovam9 <- Anova(m9)
anovam9
```
## Pairwise Comparisons 


```{r}
# looking at the effect of rate - significantly more gasps for 150ms than 75ms despite the cycle taking longer for 75m/s - INTERESTING 

emm1 <- emmeans(m9, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1, type = "response")
```

```{r}
# this effect of rate is most prominent for the gradual curves 
emm2 <- emmeans(m9, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm2, type = "response")
```
```{r}
# Looking at the effect of shape within rate 
emm <- emmeans(m9, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm, type = "response")
```


```{r}
emm3 <- emmeans(m9, list(pairwise ~ Shape|Rate), adjust = "tukey")
summary(emm3, type = "response")
```
```{r}
# this shows that c57s gasp more for the 150m/s rate compared to the 75m/s rate
emm4 <- emmeans(m9, list(pairwise ~ Rate|Strain), adjust = "tukey")
summary(emm4, type = "response")
```
```{r}
# we  see no difference between the strains in total gasps when looking at each rate independently 
emm5 <- emmeans(m9, list(pairwise ~ Strain|Rate), adjust = "tukey")
summary(emm5, type = "response")
```

```{r}
# this shows that for each curve type we saw no effect between the strains in total gasps 
emm6 <- emmeans(m9, list(pairwise ~ Strain|Shape), adjust = "tukey")
summary(emm6, type = "response")
```
```{r}
# shows no effect here either 
emm7 <- emmeans(m9, list(pairwise ~ Shape|Strain), adjust = "tukey")
summary(emm7, type = "response")
```
## plots for Total Gasps 
An important thing to remember when looking at these graphs with the stats results (e.g we see a main effect of rate from the model) is the fact that we have included our offset for cycle length and therefore it can be misleading with them plotted i.e from the graphs below we see no difference in total gasps between rate or shape however, when considering our cycle length is longer for 75m/s the fact that we see the same amount of gasps for the faster rate with a shorter cycle length means that we do have an effect of time. i.e. the same number of gasps in a shorter period of time. 

```{r}
TotalGaspsDataAll_strain_curvetype <- BehaviourExperiment1Mice[, c(5,6, 8,30)]
TotalGaspsDataAll_strain_curvetype$CurveType <- factor(TotalGaspsDataAll_strain_curvetype$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
TotalGaspsDataAll_strain_curvetype$Rate <- factor(TotalGaspsDataAll_strain_curvetype$Rate, levels = c("75", "150"))
# Coloured box and dots for strain by curve type
p.9 <- ggplot(TotalGaspsDataAll_strain_curvetype, aes(x=CurveType, y=TotalGasps, fill=Strain)) +
  geom_boxplot(position=position_dodge(1)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center',
               position=position_dodge(1),dotsize = 0.5)+ labs(x="Decompression Curve",y="Total number of gasps ")
p.9
```
```{r}
# Just curve type boxplot with coloured dots for strain
p3.9 <- ggplot(TotalGaspsDataAll_strain_curvetype, aes(x=CurveType, y=TotalGasps, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Total number of gasps ")

p3.9

```
## Latency for cessation of rhythmic breathing 

### model 10 

```{r}
m10 <- lmer(CessRBreathing_s ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(m10)
AIC(m10) 
```
```{r}
simulationOutputm10 <- simulateResiduals(fittedModel = m10, plot = T)
```

```{r}
hist(simulationOutputm10)
```

### Model 11 
```{r}
# non significant effect of including weight as a covariate and if include as a random effect the model struggles for output. If remove from the model altogether we get better dispersion - I also tried to fit using nbinom1 but this gave worse fit than nbinom2.
m11 <- glmmTMB(CessRBreathing_s ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice, family= nbinom2) 
summary(m11)
AIC(m11) 
```

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

```{r}
hist(simulationOutputm11)
```
```{r}
# shows overdispersion still - need to try and solve this still 
testDispersion(m11)
```
```{r}
#Need to include temp and Rh as covariates - both non sig therefore okay to use the model 11 above - checked it doesn't make the fit any better by including - dispersion is still the same and sig quantile deviations
m11F <- glmmTMB(CessRBreathing_s ~ Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice, family= nbinom2) 
summary(m11F)
AIC(m11F) 

```
```{r}
# can see that includingCoV (temp and Rh) does not improve the fit of the model and are non-sig 
simulationOutputm11F <- simulateResiduals(fittedModel = m11F, plot = T)
hist(simulationOutputm11F)
testDispersion(m11F)
```



## Model output m11 (latency to cessation of rhythmic breathing)
This does not give the perfect model as we still have overdispersion - need to try and find an alternative 
```{r}
anovam11 <- Anova(m11)
anovam11
```

## Pairwise Comparisons 

```{r}
# this shows that latency to stop rhythmic breathing is longer with the 75m/s rate
emm1 <- emmeans(m11, list(pairwise ~ Rate), adjust = "tukey")
summary(emm1, type = "response")
```



```{r}
# balbs take longer for cessation of rhythmic breathing 
emm2 <- emmeans(m11, list(pairwise ~ Strain), adjust = "tukey")
summary(emm2, type="response")
```
```{r}
# longer latencies for cessation of breathing for gradual 75 compared to 150 for gradual curve only. 
emm3 <- emmeans(m11, list(pairwise ~ Rate|Shape), adjust = "tukey")
summary(emm3, type="response")
```



```{r}
# Looking at the effect of shape within rate this shows longer latencies for the grdaul curve compared to both accelerated and linear. 
emm4 <- emmeans(m11, list(pairwise ~ Shape|Rate), adjust ="tukey")
summary(emm4, type = "response")
```
## Plots for latency for cessation of rhythmic breathing 


```{r}
CessRBreathingDataAll_strain_curvetype <- BehaviourExperiment1Mice[, c(5,6,8,27)]
CessRBreathingDataAll_strain_curvetype$CurveType <- factor(CessRBreathingDataAll_strain_curvetype$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
CessRBreathingDataAll_strain_curvetype$Rate <- factor(CessRBreathingDataAll_strain_curvetype$Rate, levels = c("75", "150"))
# Coloured box and dots for strain by curve type
p.7 <- ggplot(CessRBreathingDataAll_strain_curvetype, aes(x=CurveType, y=CessRBreathing_s, fill=Strain)) +
  geom_boxplot(position=position_dodge(1)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center',
               position=position_dodge(1),dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to cessation of rhythmic breathing (s)")
p.7
```

```{r}
p3.7 <- ggplot(CessRBreathingDataAll_strain_curvetype, aes(x=CurveType, y=CessRBreathing_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to cessation of rhythmic breathing (s)")

p3.7
```

# Physiology - HR data  
here I took the raw HR values every 5 secs throughout the cycle length. I then calculate the average baseline per mouse and compute the change from baseline for each 5s interval. I then found the latency when the HR dropped below 50% of the baseline for each individual and then the point at which it was sustained below 50% of baseline. This was defined where we saw >50% decrease relative to baseline across at least 60s. 

```{r}
#load in the dataset
 library(readxl)
 HRLatencies <- read_excel("J:\\BAHCM\\LAPSProject\\Exp1MiceAnalysis\\HRLatencies.xlsx")


# allows you to view the data
 View(HRLatencies)
```


```{r}
# Label the dataframe/factors
HRLatencies$Strain <- as.factor(HRLatencies$Strain)
HRLatencies$Weight <- as.factor(HRLatencies$Weight)
HRLatencies$Rate <- as.factor(HRLatencies$Rate)
HRLatencies$Shape <- as.factor(HRLatencies$Shape)

```


## Model 1 - linear model 

```{r}
# Fit a linear mixed model with weight as covariate, rate, shape and strain as fixed factors and weight as a random efect (could consider getting rid of this)
m1 <-lmer(latencyFirstMinus50 ~  Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies) 
summary(m1)
```

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




##Model 2 - negative binomial model 
to see if this improves the fit of the data 
```{r}
m2 <- glmmTMB(latencyFirstMinus50 ~  Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies, family= nbinom2) 
summary(m2)
```


```{r}
# looking at the fit of the residuals model - a better fit 
simulationOutputm2<- simulateResiduals(fittedModel = m2, plot = T)
plot(simulationOutputm2)
```
```{r}
hist(simulationOutputm2)
```
```{r}
# include temp and Rh as covariates to check that they dont contribute some of the variability and check that it doesnt improve fit - nonsig therefore okay to drop as above model
m2F <- glmmTMB(latencyFirstMinus50 ~ Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies, family= nbinom2) 
summary(m2F)
```
```{r}
simulationOutputm2F<- simulateResiduals(fittedModel = m2F, plot = T)
plot(simulationOutputm2F)
hist(simulationOutputm2F)
```




```{r}
anovam2 <-Anova(m2)
anovam2
```

## pairwise comparisons

```{r}
#look at the differences between the strains - bakbs tend to have a higher latency to drop below the 50% threshold than C57s

emm1 <- emmeans(m2, list(pairwise ~ Strain), adjust ="tukey")
summary(emm1, type = "response")
```


```{r}
# look at the mean for rate - trend towards longer latency to bradycardia in the 75ms rate
emm2 <- emmeans(m2, list(pairwise ~ Rate), adjust ="tukey")
summary(emm2, type = "response")

```


```{r}
# pairwise for the shapes 
emm3 <- emmeans(m2, list(pairwise ~ Shape), adjust ="tukey")
summary(emm3, type = "response")

```

```{r}
# looked at this interaction both ways - no difference in latency between any of the rates and shapes, irrespective of which way we look at the interaction
emm4 <- emmeans(m2, list(pairwise ~ Shape|Rate), adjust ="tukey")
summary(emm4, type = "response")
```
```{r}
# this shows that the difference between the strains appears at the lower 75m/s rate. Balbs take longer to bradycardia than c57.
emm5 <- emmeans(m2, list(pairwise ~ Strain|Rate), adjust ="tukey")
summary(emm5, type = "response")

```

```{r}
HRLatencies$CurveType <- factor(HRLatencies$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
HRLatencies$Rate <- factor(HRLatencies$Rate, levels = c("75", "150"))
h1 <- ggplot(HRLatencies, aes(x=CurveType, y=latencyFirstMinus50, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + 
  labs(title = "First indication of bradycardia", x="Decompression Profile",y="Latency (s) to first 50% reduction in HR")

h1
```




## sustained bradycardia 
This was calculated where the latency to drop below 50% of baseline  HR was sustained for 60s or more. 

```{r}
# Fit a linear mixed model with rate, shape and strain as fixed factors and weight as a random efect (could consider getting rid of this)
m3 <-lmer(sustainedLatency ~ Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies) 
summary(m3)
```

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

```{r}
hist(simulationOutputm3)
```
```{r}
m4 <-glmmTMB(sustainedLatency ~  Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies, family=nbinom2) 
summary(m4)
```
```{r}
# this fits much better 
simulationOutputm4 <- simulationOutputm4<- simulateResiduals(fittedModel = m4, plot = T)
plot(simulationOutputm4)
```
```{r}
hist(simulationOutputm4)
```
```{r}
# check that temp and RH arent having an effect when included as covariates - no effect when included in the model so okay to drop them as above.
mF4 <-glmmTMB(sustainedLatency ~ Temperature + RelativeHumidity+ Rate + Shape + Strain + Rate:Shape + Rate: Strain + Shape: Strain + (1|Weight), data=HRLatencies, family=nbinom2) 
summary(m4F)

```
```{r}
simulationOutputm4F<- simulateResiduals(fittedModel = m4F, plot = T)
plot(simulationOutputm4F)
hist(simulationOutputm4F)
```




## MOdel output 

```{r}
# this shows that rate has a main effect - using model 4 
anovam4 <- Anova(m4)
anovam4
```
# pairwise comparisons


```{r}
# shows that 75ms rate takes longer to get sustained bradycardia compared to the 150ms rate 
emm1 <- emmeans(m4, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response")
```
```{r}
emm2 <- emmeans(m4, list(pairwise ~Strain), adjust ="tukey")
summary(emm2, type = "response")
```

```{r}
emm3 <- emmeans(m4, list(pairwise ~ Shape|Rate), adjust ="tukey")
summary(emm3, type = "response")
```

```{r}
emm4 <- emmeans(m4, list(pairwise ~ Strain|Rate), adjust ="tukey")
summary(emm4, type = "response")
```

```{r}
emm5 <- emmeans(m4, list(pairwise ~ Strain|Shape), adjust ="tukey")
summary(emm5, type = "response")
```

```{r}
HRLatencies$sustainedLatency <- as.numeric(HRLatencies$sustainedLatency)
HRLatencies$CurveType <- factor(HRLatencies$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
HRLatencies$Rate <- factor(HRLatencies$Rate, levels = c("75", "150"))
h2 <- ggplot(HRLatencies, aes(x=CurveType, y=sustainedLatency, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) +  
  labs(title = "Sustained bradycardia", x="Decompression Profile",y="Latency(s) to sustained 50% reduction in HR") # this was defined as <50% of baseline for 60s or more -not all mice showed this - 19 removed

h2
```
```{r}
# try and get the means for each curve type to see differences between the latencies 

Accelerated75 <- subset(HRLatencies, CurveType =='Accel 75')
Accelerated150 <- subset(HRLatencies, CurveType =='Accel 150')
Gradual75 <- subset(HRLatencies, CurveType =='Gradual 75')
Gradual150 <- subset(HRLatencies, CurveType =='Gradual 150')
Linear75 <- subset(HRLatencies, CurveType =='Linear 75')
Linear150 <- subset(HRLatencies, CurveType =='Linear 150')

na.omit(Accelerated75, Accelerated150, Gradual150, Gradual75, Linear150, Linear75)

attach(Accelerated75)
vars<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars)
describe(vars, skew = FALSE)

attach(Accelerated150)
vars1<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars1)
describe(vars1, skew = FALSE)

attach(Gradual75)
vars2<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars2)
describe(vars2, skew = FALSE)

attach(Gradual150)
vars3<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars3)
describe(vars3, skew = FALSE)

attach(Linear75)
vars4<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars4)
describe(vars4, skew = FALSE)

attach(Linear150)
vars5<-cbind(latencyFirstMinus50,sustainedLatency)
summary(vars5)
describe(vars5, skew = FALSE)

```



# Latency to Spontaneous Gasps

This looks at the latency to first spontaneous gasp - this should be pretty similar to the first gasp data as most gasps were defined spotaneous unless they were late in the cycle and were defined as agonal and therefore would not have a value here. 

### fit a linear mixed model and look at the fit.

```{r}
# Fit a linear mixed model with weight as covariate, rate, shape and strain as fixed factors and cage nested within batch as random effects
s1 <-lmer(LatencySporadicGasp_s ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(s1)
```

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


```{r}

hist(simulationOutputs1)
```

Not a great fit therefore want to look at fitting a different model to improve quantile deviations.

## Model 2 - negative binomial mixed model 

```{r}
s2 <- glmmTMB(LatencySporadicGasp_s ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice, family=nbinom2())
summary(s2)
```

```{r}
# much better fit in terms of residual error 
simulationOutputs2<- simulateResiduals(fittedModel = s2, plot = T)
plot(simulationOutputs2)
```

```{r}
hist(simulationOutputs2)
```
```{r}
# include temp and RH as covariates - no effect so okay to drop from the model as above. 
s2F <- glmmTMB(LatencySporadicGasp_s ~ Temperature + RelativeHumidity + Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice, family=nbinom2())
summary(s2F)

```



```{r}
# model outout - main effect of rate 
anovas2 <- Anova(s2)
anovas2
```
## Pairwise comparisons to look at differences between our factors 

```{r}
emm1 <- emmeans(s2, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```
```{r}
emm2 <- emmeans(s2, list(pairwise ~ Shape), adjust="tukey")
summary (emm2, type = "response")
```

```{r}
emm3 <- emmeans(s2, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s2, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```
```{r}
emm5 <- emmeans(s2, list(pairwise ~Shape|Rate), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s2, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```

## Plots of the data 
```{r}
LatencySporadicAlldata <- BehaviourExperiment1Mice[,c(5,6,8,18)]
LatencySporadicAlldata$CurveType <- factor(LatencySporadicAlldata$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p1 <- ggplot(LatencySporadicAlldata, aes(x=CurveType, y=LatencySporadicGasp_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to sponaneous gasp (s)")

p1



```
## Number of sporadic gasps data 
This includes the number of spradic gasps for the data - most of the initial gasps tended to be sporadic before becoming rhythmic then agonal. 

```{r}
# Fit a linear mixed model with poisson distribution with weight temp and humidity as covariates, rate, shape and strain as fixed factors and cage nested within batch as random effects- this was a terrible fit 
s4 <-glmmTMB(SporadicGasp ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = poisson) 
summary(s4)

```

```{r}
# this has a terrible fit 
simulationOutputs4 <- simulateResiduals(fittedModel = s4, plot = T)
plot(simulationOutputs4)
```
```{r}
hist(simulationOutputs4)
testDispersion(simulationOutputs4)
```
# next model - quasipoisson distribution - convergence problem 

```{r}
# tried removing covariates and doesnt help the convergence problem indicated. 
s5 <-glmmTMB(SporadicGasp ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom1()) 
summary(s5)
```

```{r}
simulationOutputs5 <- simulateResiduals(fittedModel = s5, plot = T ) # very overdispersed. 
testDispersion(s5)
```
## next model negative binomimal model 

```{r}
s6 <-glmmTMB(SporadicGasp ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom2()) 
summary(s6) # temp, RH and weight are non-sig - drop RH and temp from model 
```
```{r}
s7 <-glmmTMB(SporadicGasp ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage) + (1|Weight), data=BehaviourExperiment1Mice, family = nbinom2()) 
summary(s7)

```
```{r}
#nice fit - better fit than s6 with the covraiates included - therefore warrants removing  
simulationOutputs7 <- simulateResiduals(fittedModel = s7, plot = T ) 
hist(simulationOutputs7)
testDispersion(s7) # looks good - not overdispersed 
```


## Model Output 
```{r}
# model output - main effect of rate 
anovas7 <- Anova(s7)
anovas7
```
## Pairwise comparisons to look at differences between our factors 

```{r}
emm1 <- emmeans(s7, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```

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

```{r}
emm3 <- emmeans(s7, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s7, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```
```{r}
emm5 <- emmeans(s7, list(pairwise ~Shape|Rate), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s7, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```
```{r}
emm7 <- emmeans(s7, list(pairwise ~Strain|Shape), adjust='tukey')
summary(emm7, type="response")
```

## Plots of the data 
```{r}
LatencySporadicAlldata <- BehaviourExperiment1Mice[,c(5,6,8,19)]
LatencySporadicAlldata$CurveType <- factor(LatencySporadicAlldata$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p2 <- ggplot(LatencySporadicAlldata, aes(x=CurveType, y=SporadicGasp, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Curve",y="Total number of sporadic Gasps")

p2

```


## Latency to Rhythmic Gasp
This tended to come at the same time of cessation of rhythmic breathing where instead the animal would show rhythmic gasping instead before agonal gasping started. 


```{r} 
# linear mixed model - can drop temp and RH as non-sig but keep weight as covariate 
s8 <-lmer(LatencyRhythmicGasp_s ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(s8)
```

```{r}
# not a bad fit - quite nice 
simulationOutputs8 <-simulateResiduals(fittedModel = s8, plot = T)
hist(simulationOutputs8)
```

# try this with temp and rh removed as covariates 

```{r}
s9 <-lmer(LatencyRhythmicGasp_s ~ Weight + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(s9)
```
```{r}
# slightly better fit in terms of dispersion 
simulationOutputs9 <-simulateResiduals(fittedModel = s9, plot = T)
hist(simulationOutputs9)
```
# main output from model 

```{r}
# model output
anovas9 <- Anova(s9)
anovas9
```
## Pairwise comparisons to look at differences between our factors 

```{r}
emm1 <- emmeans(s9, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```

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

```{r}
emm3 <- emmeans(s9, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s9, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```
```{r}
emm5 <- emmeans(s9, list(pairwise ~Shape|Rate), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s9, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```
```{r}
emm7 <- emmeans(s9, list(pairwise ~Strain|Shape), adjust='tukey')
summary(emm7, type="response")

```

## Plots of the data 
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,21)]
x$CurveType <- factor(LatencySporadicAlldata$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p2 <- ggplot(x, aes(x=CurveType, y=LatencyRhythmicGasp_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Curve",y="Latency to rhythmic gasp (s)")

p2

```

## number of rhythmic gasps 
This again is taken into account in the total gasp data - might not be meaningful on its own. 

```{r}
# poisson distribution was a terrible fit. Therefore tried a quasi poisson distribution - much better - can consider removing temp and RH and weight as covariates 
s10 <-glmmTMB(RhythmicGasp ~ Weight + Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage), data=BehaviourExperiment1Mice, family = nbinom1)
summary(s10)

```

```{r}
# better fit than a poisson but still not ideal
simulationOutputs10 <- simulateResiduals(fittedModel=s10, plot=T)
```
```{r}
s11 <-glmmTMB(RhythmicGasp ~Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage) + (1|Weight), data=BehaviourExperiment1Mice, family = nbinom1)
summary(s11)
```
```{r}
# better fit with covariates removed 
simulationOutputs11 <- simulateResiduals(fittedModel=s11, plot=T)

```
# next try a negative binomial model to see if this improves the fit 

```{r}
s12 <-glmmTMB(RhythmicGasp ~Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain  + offset(log(CycleLength_s)) + (1|Batch/Cage) + (1|Weight), data=BehaviourExperiment1Mice, family = nbinom2)
summary(s12)
```
```{r}
# worse fit therefore model s11 provides the best fit 
simulationOutputs12 <- simulateResiduals(fittedModel=s12, plot=T)

```
```{r}
# model output
anovas11 <- Anova(s11)
anovas11
```
## Pairwise comparisons to look at differences between our factors 

```{r}
emm1 <- emmeans(s11, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```

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

```{r}
emm3 <- emmeans(s11, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s11, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```
```{r}
emm5 <- emmeans(s11, list(pairwise ~Shape|Rate), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s11, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```
```{r}
emm7 <- emmeans(s11, list(pairwise ~Strain|Shape), adjust='tukey')
summary(emm7, type="response")

```

## Plots of the data 
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,22)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=RhythmicGasp, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Curve",y="Total numbber of rhythmic gasps")

p3
```

## Urination
These next parameters have lots of missing data points i.e. not all mice urinated so wont have a latency - this is important in how we report the findings. first have to consider how many mice urinated and then when we report the findings from the analysis below we must say of those animals that did urinate ....


```{r}
sapply(BehaviourExperiment1Mice, function(x) sum(is.na(x))) # this returns the number of missing values in each column. so for Urination - we have 37 missing values.35 mice urinated etc... 
```

#First thing to do is to look at whether a certain curve profile is more likely to result in urination
# Binary Logistic Regression 


```{r}
logit_1 <- glm(Urinate ~ Rate + Shape + Strain, family=binomial(), data = BehaviourExperiment1Mice)
summary(logit_1)
```

```{r}
Anova(logit_1)
```
```{r}
emm1 <- emmeans(logit_1, list(pairwise ~ Strain), adjust ="tukey")
summary(emm1, type = "response")
```


```{r}
# mixed model - weight was non-sig so moved to random effect 
s13 <-lmer(Urination_s ~ Temperature + RelativeHumidity + Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage) + (1|Weight), data=BehaviourExperiment1Mice) 
summary(s13) 

```
```{r}
# not a great fit 
simulationOutputs13<- simulateResiduals(fittedModel = s13, plot = T)
plot(simulationOutputs13)
hist(simulationOutputs13)

```
```{r}
#struggling for output but better fit 
s14 <-glmmTMB(Urination_s ~  Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage) + (1|Weight)+(1|Temperature) + (1|RelativeHumidity), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(s14) 

```

```{r}
simulationOutputs14<- simulateResiduals(fittedModel = s14, plot = T)
plot(simulationOutputs14)
hist(simulationOutputs14)
```
```{r}
#taking the covariates out it now gives output 
s15 <-glmmTMB(Urination_s ~  Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage) + (1|Weight), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(s15) 
```
```{r}
simulationOutputs15<- simulateResiduals(fittedModel = s15, plot = T)
plot(simulationOutputs15)
hist(simulationOutputs15)
```
```{r}
Anova(s15) # best fit model 
```
## Pairwise comparisons to look at differences between our factors 

```{r}
emm1 <- emmeans(s15, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```

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

```{r}
emm3 <- emmeans(s15, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s15, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```


```{r}
emm5 <- emmeans(s15, list(pairwise ~Shape|Rate), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s15, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```
```{r}
emm7 <- emmeans(s15, list(pairwise ~Strain|Shape), adjust='tukey')
summary(emm7, type="response")

```

## Plots of the data 
As can see from the graph any strain effects reproted could be because of unequal representation across the two strains ... 
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,32)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=Urination_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to urination (s)")

p3
```
## Latency for Fore limb movement 
there are 44 NAs therefore only 28 mice showed this - again  all mice did not show this behaviour and again need to consider unequal representation across the two strains. 

```{r}
#Binary Logistic Regression
logit_2 <- glm(Fmovement ~ Rate + Shape + Strain, family=binomial(), data = BehaviourExperiment1Mice)
summary(logit_2)
```
```{r}
Anova(logit_2)
```
```{r}
emm1 <- emmeans(logit_2, list(pairwise ~Rate), adjust='tukey')
summary(emm1, type="response")

```

```{r}
# again struggling for output therefore include temp and Rh as random effects rather than covariates
s16 <-glmmTMB(ForeLimbMovement_s ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage) + (1|Weight) + (1|Temperature) + (1|RelativeHumidity), data=BehaviourExperiment1Mice, family = nbinom2) 
summary(s16) 
```
```{r}
# looks like an okay fit 
simulationOutputs16 <- simulateResiduals(fittedModel= s16, Plot = T)
plot(simulationOutputs16)
hist(simulationOutputs16)
```
```{r}
Anova(s16)
```

## Pairwise comparisons 

```{r}
emm1 <- emmeans(s16, list(pairwise ~ Rate), adjust ="tukey")
summary(emm1, type = "response") # this transfroms the data back into our original scale from the transfromed log scale 
```

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

```{r}
emm3 <- emmeans(s16, list(pairwise ~Strain), adjust="tukey")
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(s16, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm4, type="response")
```


```{r}
emm5 <- emmeans(s16, list(pairwise ~Rate|Shape), adjust='tukey')
summary(emm5, type="response")

```

```{r}
emm6 <- emmeans(s16, list(pairwise ~Strain|Rate), adjust='tukey')
summary(emm6, type="response")
```
```{r}
emm7 <- emmeans(s16, list(pairwise ~Strain|Shape), adjust='tukey')
summary(emm7, type="response")

```

## Plots of the data 
As can see from the graph any strain effects reported could be because of unequal representation across the two strains ... interesting here because could also look at whether a certain curve profile is associated with the presence of fore limb movement i.e. gradual 150 looks much more likely to have fore limb movement (could be a proxy for seizures)
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,41)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=ForeLimbMovement_s, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression Profile",y="Latency to fore limb movement (s)")

p3
```
## Looking at the rate of gasping within the first 4 mins
This includes the total number of spontaneous and rhythmic gasps in the first 4 mins and also the rate per min. 

This does not include agonal gasp data. This was obtained by a new observer analysis specifying intervals within the results section. 

First 5 min was too long an interval whereas - 4 minutes most animals have ceased rhythmic breathing and most have started agonal breathing therefore likely conservatively encompasses the anticipated conscious period.

###quick plots
### Number of gasps in the first 4 minutes 

```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,46)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=Gasp4min, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression profile",y="Total number of gasps within first 4 min")

p3

```
## rate of gasping per minute in these first four mins. 
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,47)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=RateGasp, fill=Strain)) +
  geom_boxplot(aes(middle = mean(RateGasp)),fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression profile",y="Rate of gasping per min for first 4 min")

p3
```
## Analysis of this data - Total number of gasps in the first 4 mins 

## fit a linear mixed model 
unlike previous total gasp we do not need to include an offset as this takes into account time i.e. all in a time frame of four mins. 
```{r}
# included weight, temp as RH as covariates and non-sig therefore removed 
g1 <-lmer(Gasp4min ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Weight)+(1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(g1)
```


```{r}
# not a great fit. 
simulationOutputg1 <- simulateResiduals(fittedModel= g1, Plot = T)
plot(simulationOutputg1)
hist(simulationOutputg1)
```


```{r}
g2 <-glmmTMB(Gasp4min ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain +(1|Batch/Cage), data=BehaviourExperiment1Mice, family=nbinom2()) 
summary(g2)
```
```{r}
#better fit 
simulationOutputg2 <- simulateResiduals(fittedModel= g2, Plot = T)
plot(simulationOutputg2)
hist(simulationOutputg2)
```
# main effects - main effect of rate. 
```{r}
Anova(g2)
```


## pairwise comparisons 

```{r}
# more gasps in the first 4 minutes for the faster rate. 
emm1 <- emmeans(g2, list(pairwise ~Rate), adjust='tukey')
summary(emm1, type="response")

```
```{r}
emm2 <- emmeans(g2, list(pairwise~Shape), adjust='tukey')
summary(emm2, type="response")
```
```{r}
emm2 <- emmeans(g2, list(pairwise~Strain), adjust='tukey')
summary(emm2, type="response")
```

```{r}
emm2 <- emmeans(g2, list(pairwise~Rate|Shape), adjust='tukey')
summary(emm2, type="response")
```
```{r}
emm2 <- emmeans(g2, list(pairwise~Shape|Rate), adjust='tukey')
summary(emm2, type="response")

```
## analysis for the rate per minute data - will reflect the exact same as the above. 

```{r}
# included weight, temp as RH as covariates and non-sig therefore removed 
g3 <-lmer(RateGasp ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(g3)
```

```{r}
# okay but not great
simulationOutputg3 <- simulateResiduals(fittedModel=g3, plot=T)
plot(simulationOutputg3)
hist(simulationOutputg3)
```

```{r}
g4 <-glmmTMB(RateGasp ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Batch/Cage), data=BehaviourExperiment1Mice, family=nbinom2()) 
summary(g4)
```
```{r}
#nice fit 
simulationOutputg4 <- simulateResiduals(fittedModel=g4, plot=T)
plot(simulationOutputg4)
hist(simulationOutputg4)
```

## main effects - main effect of rate
```{r}
Anova(g4)
```
## Pairwise comparisons 

```{r}
# more gasping per minute for the faster 150ms rate. 
emm1 <- emmeans(g4, list(pairwise~Rate), adjust='tukey')
summary(emm1, type="response")
```
```{r}
emm2 <- emmeans(g4, list(pairwise~Shape), adjust='tukey')
summary(emm2, type="response")
```

```{r}
emm3 <- emmeans(g4, list(pairwise~Strain), adjust='tukey')
summary(emm3, type="response")
```

```{r}
emm4 <- emmeans(g4, list(pairwise~Shape|Rate), adjust='tukey')
summary(emm4, type="response")
```
```{r}
emm5 <- emmeans(g4, list(pairwise~Rate|Shape), adjust='tukey')
summary(emm5, type="response")
```

## Gasping prior to cessation of rhythmic breathing 

this has been calculated for each individual mouse rather than using a blanket of 4 minutes 

# total number of gasps 
```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,48)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=GaspCessRB, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression profile",y="Total number of gasps prior to CRB")

p3

```

```{r}
x <- BehaviourExperiment1Mice[,c(5,6,8,49)]
x$CurveType <- factor(BehaviourExperiment1Mice$CurveType, levels = c("Accel 75", "Linear 75", "Gradual 75", "Accel 150", "Linear 150", "Gradual 150"))
p3 <- ggplot(x, aes(x=CurveType, y=RateCessRB_sec, fill=Strain)) +
  geom_boxplot(fill="white") +
  geom_dotplot(binaxis = 'y', stackdir = 'center',dotsize = 0.5) + labs(x="Decompression profile",y="Rate of gasping")

p3
```
## analysis in line with previous analysis - looking at frequency of gasps prior to cess r breathing

```{r}
g2 <-glmmTMB(GaspCessRB ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain +(1|Batch/Cage), data=BehaviourExperiment1Mice, family=nbinom2()) 
summary(g2)
```

```{r}
simulationOutputg2 <- simulateResiduals(fittedModel= g2, Plot = T)
plot(simulationOutputg2)
hist(simulationOutputg2)
```
```{r}
Anova(g2)

# more gasps in the first 4 minutes for the faster rate. however no more gasping prior to cessr overall between our two rates.  
emm1 <- emmeans(g2, list(pairwise ~Rate), adjust='tukey')
summary(emm1, type="response")
emm2 <- emmeans(g2, list(pairwise ~Rate|Strain), adjust='tukey') # more gasping prior ot cess RB in C57s only. 
summary(emm2, type="response")
```

# rate of gasping (gasps/duration Cess B)

```{r}
g3 <-lmer(RateCessRB_sec ~ Rate + Shape + Strain + Rate:Shape + Rate:Strain + Shape:Strain + (1|Weight) + (1|Batch/Cage), data=BehaviourExperiment1Mice) 
summary(g3)
```

```{r}
simulationOutputg3 <- simulateResiduals(fittedModel= g3, Plot = T)
plot(simulationOutputg3)
hist(simulationOutputg3)
```

```{r}

Anova(g3)

# more gasps in the first 4 minutes for the faster rate. however no more gasping prior to cess overall between our two rates.  
emm1 <- emmeans(g3, list(pairwise ~Rate), adjust='tukey')
summary(emm1, type="response")
emm2 <- emmeans(g3, list(pairwise ~Rate|Strain), adjust='tukey') # more gasping prior to cess RB in C57s only. 
summary(emm2, type="response")
```