---
title: "NewAnalysis2"
author: "Jasmine Clarkson"
date: "10/02/2022"
output:
  html_document: default
  pdf_document: default
---

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



```{r}
library(readxl)
library(ggplot2)
library(plyr)
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(ordinal)
library (descr)
library(nnet)
library(VGAM)
library(RVAideMemoire)
```


## Import sorted dataset - only interested in complete survey responses. 219 participants total.

```{r}
d<- RodentSurveyFeb2022 <- read_excel("J:/BAHCM/LAPSProject/Survey Documents/Final Datasets/RodentSurveyFeb2022.xlsx", 
    sheet = "Completed Only - with consent")
View(RodentSurveyFeb2022)
```

```{r}
d$Sector <- as.factor(d$Sector)
d$Role_New <- as.factor(d$Role_New)
d$Experience_New <- as.factor(d$Experience_New)
d$Species <- as.factor(d$Species)
d$EuthanasiaUse <- as.factor(d$EuthanasiaUse)
d$Preference_Top <- as.factor(d$Preference_Top)
d$Preference_Multi_Top <- as.factor(d$Preference_Multi_Top)
d$Preference_BFT <-as.factor(d$Preference_BFT)
d$Preference_CO2 <-as.factor(d$Preference_CO2)
d$Preference_CD <-as.factor(d$Preference_CD)
d$Preference_ODA <-as.factor(d$Preference_ODA)
d$Confident_BFT <-as.factor(d$Confident_BFT)
d$Confident_CO2 <-as.factor(d$Confident_CO2)
d$Confident_CD <-as.factor(d$Confident_CD)
d$Confident_ODA <-as.factor(d$Confident_ODA)
d$Confident_ALL <-as.factor(d$Confident_ALL)
d$EuthanasiaMethods_ODA <-as.factor(d$EuthanasiaMethods_ODA)
d$EuthanasiaMethods_CD <-as.factor(d$EuthanasiaMethods_CD)
d$EuthanasiaMethods_CO2 <-as.factor(d$EuthanasiaMethods_CO2)
d$EuthanasiaMethods_BFT <-as.factor(d$EuthanasiaMethods_BFT)
d$EuthanasiaMethods_all <-as.factor(d$EuthanasiaMethods_all)
d$Preference_Top_humane <-as.factor(d$Preference_Top_humane)
d$Preference_Top_PrimaryReason <- as.factor(d$Preference_Top_PrimaryReason)

```

# ordinal logistic regression 
Here we want to look at the the average (mean) or mode score that each killing method was ranked (1-4)

## Cervical disloation across sector and role

```{r}
m1 <- clm(Preference_CD ~ Role_New, threshold = "equidistant", data=d)
summary(m1)
```

```{r}
emm1 <- emmeans(m1, pairwise ~ Role_New, mode="mean.class") 
summary(emm1, type="mean.class") 
# no real differences across role - all rank CD equally high 
#emm2<- emmeans(m1, pairwise ~ Sector, mode="mean.class") 
#summary(emm2, type="mean.class") 
# no huge difference across sector in ranking 
#emm3<- emmeans(m1, pairwise ~ EuthanasiaUse, mode="mean.class") 
#summary(emm3, type="mean.class")

#  lack of values - the data is not diverse enough for the model to detect an effect - likely because most people rank CD 1 or 2 - therefore removed sector and euthanasia use from the model and focused on role only. 
```
## co2 across sector and role 

```{r}
m2 <- clm(Preference_CO2 ~ Role_New + Sector, threshold = "equidistant", data=d)
summary(m2)
```

```{r}
emm1 <- emmeans(m2, pairwise ~ Role_New, mode="mean.class") 
summary(emm1, type="mean.class") 
# no  differences across role - all rank CO2 2 or 3
emm2<- emmeans(m2, pairwise ~ Sector, mode="mean.class") 
summary(emm2, type="mean.class") 
# no  difference across sector in ranking again ranked between 2 or 3
#emm3<- emmeans(m2, pairwise ~ EuthanasiaUse, mode="mean.class") 
#summary(emm3, type="mean.class")
```

## ODA across sector and role 

```{r}
m3 <- clm(Preference_ODA ~ Role_New + Sector, threshold = "equidistant", data=d)
summary(m3)
```

```{r}
emm1 <- emmeans(m3, pairwise ~ Role_New, mode="mean.class") 
summary(emm1, type="mean.class") 
# vets rank ODA higher in terms of preference - this is signif compared to management and technicians
emm2<- emmeans(m3, pairwise ~ Sector, mode="mean.class") 
summary(emm2, type="mean.class") 
# no  difference across sector 
#emm3<- emmeans(m3, pairwise ~ EuthanasiaUse, mode="mean.class") 
#summary(emm3, type="mean.class") # no effect so removed from model 

```

## blunt force trauma according to role and sector


```{r}
m4 <- clm(Preference_BFT ~ Role_New, threshold = "equidistant", data=d)
summary(m4)

# lack of convergence due to most people ranking it number 4 
```

```{r}
emm1 <- emmeans(m4, pairwise ~ Role_New, mode="mean.class") 
summary(emm1, type="mean.class") 
# no difference across role - ranked between 3 and 4 
#emm2<- emmeans(m4, pairwise ~ Sector, mode="mean.class") 
#summary(emm2, type="mean.class") 
# no  difference across sector 
#emm3<- emmeans(m4, pairwise ~ EuthanasiaUse, mode="mean.class") 
#summary(emm3, type="mean.class")

# no standard error as most people rate this a 4 and therefore there is not enough diversity in the data across role. SEM calculated manually in excel using filter variables. 

# management 0.86
# researcher 0.04
# regulatory 0
# technician 0.06
# vet 0.19

```


# what factors affect whether an indidivual is confident in using each given method

## blunt force trauma - probabilty of being confident in it according to preference

ive removed other factors like role, sector, use as they were non-sig
```{r}
m5 <- glm(Confident_BFT ~  Preference_BFT, family ="binomial", data=d)
summary(m5)

Anova(m5)
```

```{r}
emm1 <- emmeans(m5, pairwise ~ Preference_BFT, mode="response") 
summary(emm1, type="response") 

```
# other way around - is preference affected by confidence in the method? does confidence in a method make you more likely to have a preference for it? 

```{r}

m5a <- glm(Preference_BFT ~ Confident_BFT + Role_New + Sector + EuthanasiaUse , family ="binomial", data=d)
summary(m5a)

Anova(m5a)
```
```{r}
emm1 <- emmeans(m5a, pairwise ~ Confident_BFT, mode="response") 
summary(emm1, type="response") 

```

# CO2 confidence according to characteristics and preference 
```{r}
m6 <- glm(Confident_CO2 ~ Preference_CO2, family ="binomial", data=d)
summary(m6)

Anova(m6)
```
```{r}
emm1 <- emmeans(m6, pairwise ~ Preference_CO2, mode="response") 
summary(emm1, type="response") 
# looks like likelihood of stating confident in CO2 is related to preference 
```
# other way around - is preference affected by confidence in the method? does confidence in a method make you more likely to have a preference for it? 

```{r}
m6a <- glm(Preference_CO2 ~ Confident_CO2 +Role_New + Sector + EuthanasiaUse , family ="binomial", data=d)
summary(m6a)

Anova(m6a)
```
```{r}
emm1 <- emmeans(m6a, pairwise ~ Confident_CO2, mode="response") 
summary(emm1, type="response") 
# looks like likelihood of stating confident in CO2 is related to preference 
```




# ODA confidence with preference 

```{r}
m7 <- glm(Confident_ODA ~  Preference_ODA, family ="binomial", data=d)
summary(m7)

Anova(m7)
```
```{r}
emm1 <- emmeans(m7, pairwise ~ Preference_ODA, mode="response") 
summary(emm1, type="response") 
# positive relation between liklihood of being confident in the technique (ODA) and preference
#emm2 <- emmeans(m7,pairwise ~Role_New, mode="response")
#summary(emm2, type="response")
# no difference between any of the pairings - management and vets have a greater likelihood of being confident in ODA though not significantly so
```

# other way around - is preference affected by confidence in the method? does confidence in a method make you more likely to have a preference for it? 

```{r}
m7b <- glm(Preference_ODA ~ Confident_ODA + Role_New + Sector + EuthanasiaUse , family ="binomial", data=d)
summary(m7b)

Anova(m7b)
```

```{r}
emm1 <- emmeans(m7b, pairwise ~ Confident_ODA, mode="response") 
summary(emm1, type="response") 
# looks like likelihood of stating confident in CO2 is related to preference 
```


# CD with prefernce with confidence
```{r}
m7a <- glm(Confident_CD ~  Preference_CD, family ="binomial", data=d)
summary(m7a)

Anova(m7a)
```

```{r}
emm1 <- emmeans(m7a, pairwise ~ Preference_CD, mode="response") 
summary(emm1, type="response") 
# positive relation between liklihood of being confident in the technique (CD) and preference
#emm2 <- emmeans(m7a,pairwise ~Role_New, mode="response")
#summary(emm2, type="response")
# no difference between any of the pairings 
```


# other way around - is preference affected by confidence in the method? does confidence in a method make you more likely to have a preference for it? 
```{r}
m7c <- glm(Preference_CD~Confident_CD +Role_New + Sector + EuthanasiaUse  , family ="binomial", data=d)
summary(m7c)

Anova(m7c)
```

```{r}
emm1 <- emmeans(m7c, pairwise ~ Confident_CD, mode="response") 
summary(emm1, type="response") 
# looks like likelihood of stating confident in CO2 is related to preference 
```




## is confidence affected by availability of methods? 

## blunt force trauma
```{r}
m8 <- glm(Confident_BFT ~ Role_New + Sector + EuthanasiaUse + EuthanasiaMethods_BFT, family ="binomial", data=d)
summary(m8)

Anova(m8)
```

```{r}
emm1 <- emmeans(m8, pairwise ~ EuthanasiaMethods_BFT, mode="response") 
summary(emm1, type="response") 
# more likely to be confident if its available 

emm2 <- emmeans(m8,pairwise ~Role_New, mode="response")
summary(emm2, type="response")
```

# CO2

```{r}
m9 <- glm(Confident_CO2 ~ Role_New + Sector + EuthanasiaUse + EuthanasiaMethods_CO2, family ="binomial", data=d)
summary(m9)

Anova(m9)
```

```{r}
emm1 <- emmeans(m9, pairwise ~ EuthanasiaMethods_CO2, mode="response") 
summary(emm1, type="response")
# more likely be confident if available 
```
##ODA

```{r}
m10 <- glm(Confident_ODA ~ Role_New + Sector + EuthanasiaUse + EuthanasiaMethods_ODA, family ="binomial", data=d)
summary(m10)

Anova(m10)
```
``````{r}
emm1 <- emmeans(m10, pairwise ~ EuthanasiaMethods_ODA, mode="response") 
summary(emm1, type="response")
# more likely be confident if available 
```
## CD
```{r}
m11 <- glm(Confident_CD ~ Role_New + Sector + EuthanasiaUse + EuthanasiaMethods_CD, family ="binomial", data=d)
summary(m11)

Anova(m11)
```
```{r}
emm1 <- emmeans(m11, pairwise ~ EuthanasiaMethods_CD, mode="response") 
summary(emm1, type="response") # equally likely

emm2 <- emmeans(m11, pairwise ~ Role_New, mode="response") 
summary(emm2, type="response")
```

# how many people reported that all methods were available? 
```{r}
allmethods <- d%>%
  dplyr::group_by(EuthanasiaMethods_all) %>%
  dplyr::summarise(Frequency = n()) %>%
  dplyr::mutate(Percent = round(Frequency/sum(Frequency)*100, 1))

allmethods
```

# how many people rated confidence in all 4 methods? 

```{r}
coinallmethods <- d%>%
  dplyr::group_by(Confident_ALL) %>%
  dplyr::summarise(Frequency = n()) %>%
  dplyr::mutate(Percent = round(Frequency/sum(Frequency)*100, 1))

coinallmethods
```


## Reasons given for choice of preferred method  

how does ranking of each method affect reason chosen as humane
```{r}
m12 <- clm(Preference_Top_humane ~ Preference_BFT, threshold = "equidistant", data=d)
summary(m12)

```

```{r}
emm1 <- emmeans(m12, pairwise ~ Preference_BFT, mode="mean.class") 
summary(emm1, type="mean.class")

```

```{r}
m13 <- clm(Preference_Top_humane ~ Preference_CO2, threshold = "equidistant", data=d)
summary(m13)
```

```{r}
emm1 <- emmeans(m13, pairwise ~ Preference_CO2, mode="mean.class") 
summary(emm1, type="mean.class")
# this shows the mean score for humaneness according to CO2 ranking. INterestingly individuals that rank CO2 their preferred method dont state this is because they believe for it to be humane. 
```

```{r}
m14 <- clm(Preference_Top_humane ~ Preference_CD, threshold = "equidistant", data=d)
summary(m14)
```

```{r}
emm1 <- emmeans(m14, pairwise ~ Preference_CD, mode="mean.class") 
summary(emm1, type="mean.class")

# opposite to BFT and CO2 - people who chose CD as preferred method tend to score it due to humaneness ranking high
```

```{r}
m15 <- clm(Preference_Top_humane ~ Preference_ODA, threshold = "equidistant", data=d)
summary(m15)
```

```{r}
emm1 <- emmeans(m15, pairwise ~ Preference_ODA, mode="mean.class") 
summary(emm1, type="mean.class")
# this is the same as CD - people who rank ODA their preferred method tend to do so because they believe it is humane
```

# ASSESSING WHETHER METHOD (PREFERRED) AFFECTS THE REASON 

I.E IS THE RANKING ORDER FOR HUMANENESS DIFFERENT ACROSS THE METHODS
```{r}
r1  <- clm(Preference_Top_humane ~ Preference_Top, threshold = "equidistant", data=d)
summary(r1)
```
```{r}
emm1 <- emmeans(r1, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")


```

## Fast

```{r}
d$Preference_Top_fast <- as.factor(d$Preference_Top_fast)
r2  <- clm(Preference_Top_fast ~ Preference_Top, threshold = "equidistant", data=d)
summary(r2)
```

```{r}
emm1 <- emmeans(r2, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```
# easy

```{r}
d$Preference_Top_easy <- as.factor(d$Preference_Top_easy)
r3  <- clm(Preference_Top_easy ~ Preference_Top, threshold = "equidistant", data=d)
summary(r3)
```

```{r}
emm1 <- emmeans(r3, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```
# reliability

```{r}
d$Preference_Top_reliable <- as.factor(d$Preference_Top_reliable)
r4  <- clm(Preference_Top_reliable ~ Preference_Top, threshold = "equidistant", data=d)
summary(r4)
```

```{r}
emm1 <- emmeans(r4, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```

# Throughput

```{r}
d$Preference_Top_throughput <- as.factor(d$Preference_Top_throughput)
r5  <- clm(Preference_Top_throughput ~ Preference_Top, threshold = "equidistant", data=d)
summary(r5)
```

```{r}
emm1 <- emmeans(r5, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```


```{r}
d$Preference_Top_cost <- as.factor(d$Preference_Top_cost)
r6  <- clm(Preference_Top_cost ~ Preference_Top, threshold = "equidistant", data=d)
summary(r6)
```

```{r}
emm1 <- emmeans(r6, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```






```{r}
d$Preference_Top_safe <- as.factor(d$Preference_Top_safe)
r7  <- clm(Preference_Top_safe ~ Preference_Top, threshold = "equidistant", data=d)
summary(r7)
```

```{r}
emm1 <- emmeans(r7, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```




```{r}
d$Preference_Top_equipment <- as.factor(d$Preference_Top_equipment)
r8  <- clm(Preference_Top_equipment ~ Preference_Top, threshold = "equidistant", data=d)
summary(r8)
```

```{r}
emm1 <- emmeans(r8, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```

```{r}
d$Preference_Top_training <- as.factor(d$Preference_Top_training)
r9  <- clm(Preference_Top_training ~ Preference_Top, threshold = "equidistant", data=d)
summary(r9)
```

```{r}
emm1 <- emmeans(r9, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```

```{r}
d$Preference_Top_contact <- as.factor(d$Preference_Top_contact)
r10  <- clm(Preference_Top_contact ~ Preference_Top, threshold = "equidistant", data=d)
summary(r10)
```

```{r}
emm1 <- emmeans(r10, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```

```{r}
d$Preference_Top_operator <- as.factor(d$Preference_Top_operator)
r11  <- clm(Preference_Top_operator ~ Preference_Top, threshold = "equidistant", data=d)
summary(r11)
```

```{r}
emm1 <- emmeans(r11, pairwise ~ Preference_Top, mode="mean.class") 
summary(emm1, type="mean.class")

```


## binary logistic regression
looking at factors which might influence whether knew co2 fill method and flow rate 

```{r}
subset<- read_excel("J:/BAHCM/LAPSProject/Survey Documents/Final Datasets/SubsetCO2Data.xlsx")
View(subset)

subset$Role_New <- as.factor(subset$Role_New)
subset$Species <- as.factor(subset$Species)
subset$CO2filling_DontKnow <- as.numeric(subset$CO2filling_DontKnow)
subset$KnewFillRate <- as.numeric(subset$KnewFillRate)

```


```{r}

b1 <- glm(CO2filling_DontKnow ~ Role_New + Species, family = 'binomial', data=subset)
summary(b1)
```
```{r}
simulationOutputb1 <- simulateResiduals(fittedModel = b1, plot = T)

```



```{r}
Anova(b1)

# significant effect of role on liklihood of not knowing the fill method
```
```{r}
emm1 <- emmeans(b1, pairwise ~ Role_New, mode="response") 
summary(emm1, type="response")
```

# fill rate - whether people knew or not 

```{r}
b2 <- glm(KnewFillRate ~ Role_New + Species, family = 'binomial', data=subset)
summary(b2)
```

```{r}
simulationOutputb2 <- simulateResiduals(fittedModel = b2, plot = T)

```

```{r}
Anova(b2)
```


```{r}
emm1 <- emmeans(b2, pairwise ~ Role_New, mode="response") 
summary(emm1, type="response")

# highlights huge knowledge gap in researchers - less likely to know the fill rate 
```