#load package---------------------
install_package_cr <- function(package_name) {
  for (pk in package_name){
    if(!requireNamespace(pk, quietly = TRUE)){install.packages(pk)}
  }
  invisible(lapply(package_name, require, character.only = TRUE))
}
packages <- c("readr", "ggplot2", "statmod", "multcomp", "nnet", "car", "plyr")

install_package_cr(packages)

#load data---------------
setwd("set directory to folder where the data is stored")

#Data from Muscas et. al 2019
EO_original <- read_csv("EO_original.csv")
#Data from Muscas et. al 2019 + 100mg/kg Lovastatin data from Osterweil et. al.
EO_original_100 <- read_csv("EO_original_100.csv")
#Data from Muscas et. al 2019 + all groups from Osterweil et. al.
EO_original_all <- read_csv("EO_original_all.csv")


#setting baseline factors as WT Veh, setting the order of factors-------------------

EO_original$incidence <- as.factor(EO_original$incidence)
EO_original$treatment <- factor(EO_original$treatment, level = c("veh", "lova", "simvlow", "simvhigh"))
EO_original$genotype <- factor(EO_original$genotype, level = c("WT", "KO"))
EO_original$severity <- factor(EO_original$severity, level = c("no", "run", "clonic", "tonic"))

EO_original_100$genotype <- factor(EO_original_100$genotype, levels = c("WT", "KO"))
EO_original_100$treatment <- factor(EO_original_100$treatment, levels = c("veh", "lova","lov100", "simvhigh", "simvlow"))
EO_original_100$incidence <- factor(EO_original_100$incidence, levels = c("no", "yes"))
EO_original_100$severity <- factor(EO_original_100$severity, level = c("no", "run", "clonic", "tonic"))

EO_original_all$genotype <- factor(EO_original_all$genotype, levels = c("WT", "KO"))
EO_original_all$treatment <- factor(EO_original_all$treatment, levels = c("veh", "lova","lov100", "lov100c", "lov10c", "lov30c", "simvhigh", "simvlow"))
EO_original_all$incidence <- factor(EO_original_all$incidence, levels = c("no", "yes"))
EO_original_all$severity <- factor(EO_original_all$severity, level = c("no", "run", "clonic", "tonic"))

#creating a version of Muscas et. al 2019 + 100mg/kg Lovastatin data from Osterweil et. al. with one lovastatin group------

EO_original_100_m <- EO_original_100
EO_original_100_m$treatment <- revalue(EO_original_100_m$treatment, c("lov100" = "lova"))


#Logistics regression with data from Muscas et. al 2019-----------------

#full model
ags_ALL <- glm(incidence ~ genotype*treatment, data=EO_original, family="binomial")
#model without interaction
ags_GT <- glm(incidence ~ genotype + treatment, data=EO_original, family="binomial")

#effect testing
Anova(ags_ALL,type ="II")

#posthoc testing
tmp <- expand.grid(treatment = c("veh", "lova", "simvlow", "simvhigh"),
                   genotype = c("WT", "KO"))
X <- model.matrix( ~ genotype * treatment, data = tmp)
Tukey <- contrMat(table(EO_original$treatment), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow=nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(EO_original$genotype)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow(Tukey), ncol=ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(EO_original$genotype)[2], rownames(K2), sep = ":")
K <- rbind(K1,K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

#results
posthoc_tukey <- glht(ags_ALL, linfct = K%*%X)
summary(posthoc_tukey)


#Multinomial regression with data from Muscas et. al 2019-----------------

#full model
mult_ALL <- multinom(severity ~ genotype + treatment + genotype:treatment, data = EO_original)
#model without interaction
mult_GT <- multinom(severity ~  genotype+treatment, data = EO_original)

#effect testing
Anova(mult_ALL,type ="II")
z <- summary(mult_GT)$coefficients/summary(mult_GT)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2


#Logistics regression with data from Muscas et. al 2019 + 100mg/kg data from Osterweil et al-----------------

#full model
ags_ALL_100 <- glm(incidence ~ genotype*treatment, data=EO_original_100, family="binomial")
#model without interaction
ags_GT_100 <- glm(incidence ~ genotype + treatment, data=EO_original_100, family="binomial")

#effect testing
Anova(ags_ALL_100,type ="II")

tmp <- expand.grid(treatment = c("veh", "lova","lov100", "simvhigh", "simvlow"),
                   genotype = c("WT", "KO"))
X <- model.matrix( ~ genotype * treatment, data = tmp)
Tukey <- contrMat(table(EO_original_100$treatment), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow=nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(EO_original_100$genotype)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow(Tukey), ncol=ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(EO_original_100$genotype)[2], rownames(K2), sep = ":")
K <- rbind(K1,K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

#results
posthoc_tukey_ALL_100 <- glht(ags_ALL_100, linfct = K%*%X)
summary(posthoc_tukey_ALL_100)


#Logistics regression with data from Muscas et. al 2019 + 100mg/kg data from Osterweil et al (one lovastatin group)-----------------
ags_ALL_100_m <- glm(incidence ~ genotype*treatment, data=EO_original_100_m, family="binomial")

#effect testing
Anova(ags_ALL_100_m,type ="II")

tmp <- expand.grid(treatment = c("veh", "lova", "simvlow", "simvhigh"),
                   genotype = c("WT", "KO"))
X <- model.matrix( ~ genotype * treatment, data = tmp)
Tukey <- contrMat(table(EO_original_100_m$treatment), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow=nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(EO_original_100_m$genotype)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow(Tukey), ncol=ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(EO_original_100_m$genotype)[2], rownames(K2), sep = ":")
K <- rbind(K1,K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

#results
posthoc_tukey_100_m <- glht(ags_ALL_100_m, linfct = K%*%X)
summary(posthoc_tukey_100_m)

#Logistics regression with data from Muscas et. al 2019 + all groups from Osterweil et al-----------------

#full model
ags_ALL_all <- glm(incidence ~ genotype*treatment, data=EO_original_all, family="binomial")
#model without interaction
ags_GT_all <- glm(incidence ~ genotype + treatment, data=EO_original_all, family="binomial")

#effect testing
Anova(ags_ALL_all,type ="II")

tmp <- expand.grid(treatment = c("veh", "lova","lov100", "lov100c", "lov10c", "lov30c", "simvhigh", "simvlow"),
                   genotype = c("WT", "KO"))
X <- model.matrix( ~ genotype * treatment, data = tmp)
Tukey <- contrMat(table(EO_original_all$treatment), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow=nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(EO_original_all$genotype)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow(Tukey), ncol=ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(EO_original_all$genotype)[2], rownames(K2), sep = ":")
K <- rbind(K1,K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

#results
posthoc_tukey_ALL_all <- glht(ags_ALL_all, linfct = K%*%X)
summary(posthoc_tukey_ALL_all)