---
title: "Expt 1 mice Pathology and Histology"
author: "Jasmine Clarkson"
date: "13/04/2021"
output: html_document
---

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


# load in the required packages 
```{r}
library(ggplot2)
library(MASS)
library(ordinal)
library(car)
library(ggeffects)
library(RColorBrewer)
library(dplyr)
library(emmeans)
library(readxl)
```


# load in the data set 
```{r}
d <- read_excel("Exp1HistoScores.xlsx")
View(d)
```

```{r}
# define factors in the dataset 
d$Strain <- as.factor(d$Strain)
d$CurveType <- as.factor(d$CurveType)
d$Rate <- as.factor(d$Rate)
d$Shape <- as.factor(d$Shape)
d$Area <- as.factor(d$Area)
d$Batch <- as.factor(d$Batch)
d$Cage <- as.factor(d$Cage)
```


```{r}
#descriptives to check correct n 
lapply(d[,c('CurveType', 'Strain', 'Rate', 'Shape', 'Score', 'Area', 'Cage')], table)

```

```{r}
#all results - for all areas of histology
ggplot(d, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)
  
```


```{r}
# want to subset it into each area of histology - i.e. right lung congestion etc 

RLC <- subset(d, Area=='RLC' )
RLH <- subset(d, Area=='RLH' )
RLA <- subset(d, Area=='RLA' )
LLC <- subset(d, Area=='LLC' )
LLH <- subset(d, Area=='LLH' )
LLA <- subset(d, Area=='LLA' )
MEC <- subset(d, Area=='MEC' )
MEH <- subset(d, Area=='MEH' )
MEI <- subset(d, Area=='MEI' )
TME <- subset(d, Area=='TME' )
MSC <- subset(d, Area=='MSC' )
MSH <- subset(d, Area=='MSH' )
RH <- subset(d, Area=='RH' )

```



Next want a plot for each area to see the potential differences - pooled across strain (n=8 in each)

```{r}

ggplot(RLC, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(RLH, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(RLA, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(LLC, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(LLH, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(LLA, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(MEC, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(MEH, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(MEI, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(TME, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(MSC, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)

ggplot(MSH, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)


ggplot(RH, aes(x = Shape, y = Score)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) + facet_wrap(~Rate)


```

### next we want to perform a cumulative link model - ordinal data  
need to have ordinal package loaded in for this. 

#RIGHT LUNG CONGESTION

```{r}
RLC$Score <- ordered(RLC$Score, levels = 0:5) # define Score as ordinal on a 6 point scale.

c <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=RLC, threshold = "equidistant")
summary(c)
```

```{r}
confint(c)
```

```{r}
emm1 <- emmeans(c, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c, pairwise ~ Shape|Strain, mode="mean.class") 

```

```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```


#RIGHT LUNG HAEMORRHAGE

```{r}
RLH$Score <- ordered(RLH$Score, levels = 0:5) 

c1 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain +(1|Cage), data=RLH, threshold = "equidistant")
summary(c1)
```


```{r}
confint(c1)
```

```{r}
emm1 <- emmeans(c1, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c1, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c1, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c1, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c1, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c1, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c1, pairwise ~ Shape|Strain, mode="mean.class") 
```

```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```


# Right lung atelectasis
```{r}
RLA$Score <- ordered(RLA$Score, levels = 0:5) 

c3 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=RLA, threshold = "equidistant")
summary(c3)
```


```{r}
confint(c3)
```

```{r}
emm1 <- emmeans(c3, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c3, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c3, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c3, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c3, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c3, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c3, pairwise ~ Shape|Strain, mode="mean.class") 
```

```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")

```


# left lung congestion 

```{r}
LLC$Score <- ordered(LLC$Score, levels = 0:5) 

c4 <- clmm(Score ~ Rate + Shape + Strain +  (1|Cage), data=LLC, threshold = "equidistant")
summary(c4)
```



```{r}
emm1 <- emmeans(c4, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c4, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c4, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c4, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c4, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c4, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c4, pairwise ~ Shape|Strain, mode="mean.class") 
```

```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```



#LEFT LUNG HAEMORRHAGE 

```{r}
LLH$Score <- ordered(LLH$Score, levels = 0:5) 

c5 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=LLH, threshold = "equidistant")
summary(c5)
```


```{r}
confint(c5)
```

```{r}
emm1 <- emmeans(c5, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c5, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c5, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c5, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c5, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c5, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c5, pairwise ~ Shape|Strain, mode="mean.class") 
```



```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```



#LEFT LUNG ATELEASIS

```{r}
LLA$Score <- ordered(LLA$Score, levels = 0:5) 

c6 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=LLA, threshold = "equidistant")
summary(c6)
```

```{r}
confint(c6)
```

```{r}
emm1 <- emmeans(c6, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c6, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c6, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c6, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c6, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c6, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c6, pairwise ~ Shape|Strain, mode="mean.class") 
```
```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```




#MIDDLE EAR CONGESTION 

```{r}
MEC$Score <- ordered(MEC$Score, levels = 0:5) 

c7 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=MEC, threshold = "equidistant")
summary(c7)
```

```{r}
confint(c7)
```

```{r}
emm1 <- emmeans(c7, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c7, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c7, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c7, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c7, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c7, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c7, pairwise ~ Shape|Strain, mode="mean.class") 
```




```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```

#MIDDLE EAR HAEMORRHAGE 

```{r}
MEH$Score <- ordered(MEH$Score, levels = 0:5) 

c8 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=MEH, threshold = "equidistant")
summary(c8)
```


```{r}
confint(c8)
```

```{r}
emm1 <- emmeans(c8, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c8, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c8, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c8, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c8, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c8, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c8, pairwise ~ Shape|Strain, mode="mean.class") 
```


```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```

#TYMPANIC MEM EOSINO

```{r}
TME$Score <- ordered(TME$Score, levels = 0:5) 

c9 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=TME, threshold = "equidistant")
summary(c9)
```

```{r}
confint(c9)
```

```{r}
emm1 <- emmeans(c9, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c9, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c9, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c9, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c9, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c9, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c9, pairwise ~ Shape|Strain, mode="mean.class") 
```


```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```

Gross pathological macroscopic scores 
# this does not include all scores only the ones where we see a digression of scores i.e. not all at 0 to see if we have any effect of decompression profile (rate and shape) on obtaining higher scores. 


#load in the data 
```{r}
d2 <- read_excel("Exp1PMScores.xlsx")
View(d2)

```


```{r}
# define factors in the dataset 
d2$Strain <- as.factor(d2$Strain)
d2$CurveType <- as.factor(d2$CurveType)
d2$Rate <- as.factor(d2$Rate)
d2$Shape <- as.factor(d2$Shape)
d2$Area <- as.factor(d2$Area)
d2$Cage <- as.factor(d2$Cage)
```

```{r}
Spleen <- subset(d2, Area=='Spleen' )
Stomach <- subset(d2, Area=='Stomach' )
smallIntes <- subset(d2, Area=='smallIntes' )
LeftLung <- subset(d2, Area=='LeftLung' )
RightLung <- subset(d2, Area=='RightLung' )
```


# INCREASE IN SPLEEN SIZE - not enough variation in scores to model. 
```{r}
Spleen$Score <- ordered(Spleen$Score, levels = 0:5) # define Score as ordinal on a 6 point scale.

c10 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=Spleen, link = "logit")
summary(c10)

```

```{r}
emm1 <- emmeans(c10, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c10, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c10, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c10, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c10, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c10, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c10, pairwise ~ Shape|Strain, mode="mean.class") 
```


```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```

# Small intestine data - not enough variation. 

```{r}
smallIntes$Score <- ordered(smallIntes$Score, levels = 0:5) # define Score as ordinal on a 6 point scale.

c11 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain +(1|Cage), data=smallIntes, link = "logit")
summary(c11)
```


# Right Lung discolouration - more variation (no model output for left lung due to limited variation of scores)

```{r}
RightLung$Score <- ordered(RightLung$Score, levels = 0:5) # define Score as ordinal on a 6 point scale.

c12 <- clmm(Score ~ Rate + Shape + Strain + Shape:Rate + Shape:Strain + Rate:Strain + (1|Cage), data=RightLung, threshold = "equidistant")
summary(c12)

```

```{r}
emm1 <- emmeans(c12, pairwise ~ Strain, mode="mean.class") 
emm2 <- emmeans(c12, pairwise ~ Rate, mode="mean.class") 
emm3 <- emmeans(c12, pairwise ~ Shape, mode="mean.class") 
emm4 <- emmeans(c12, pairwise ~ Rate|Shape, mode="mean.class") 
emm5 <- emmeans(c12, pairwise ~ Shape|Rate, mode="mean.class") 
emm6 <- emmeans(c12, pairwise ~ Rate|Strain, mode="mean.class") 
emm7 <- emmeans(c12, pairwise ~ Shape|Strain, mode="mean.class") 
```


```{r}
summary(emm1, type="mean.class")
summary(emm2, type="mean.class")
summary(emm3, type="mean.class")
summary(emm4, type="mean.class")
summary(emm5, type="mean.class")
summary(emm6, type="mean.class")
summary(emm7, type="mean.class")
```