rm(list=ls())
#setwd("")
require(dabestr)
require(dplyr)
require(MASS)
require(car)
require(Cairo)
require(tidyr)
require(ggplot2)


######################
# figure of thermal ramping rates
######################

# N.B. Thermal ramping data are missing for cods 24-27 and the beginning of the trials for cods 21-23 (the thermometer was erroneously turned off at the start of the trials).

#-#-#-#-#-#-#-#-#-#-#-#-#  import data #-#-#-#-#-#-#-#-#-#-#-#-#

ramp <- read.csv("Thermal ramping data.csv", header=T)
	dim(ramp)
	head(ramp)
	str(ramp)
	summary(ramp)

#switch table format from wide to long for plotting
ramp_long <- gather(ramp, codID, temperature, cod2:cod39, factor_key=TRUE)
	str(ramp_long)
	summary(ramp_long)

#-#-#-#-#-#-#-#-#-#-#-#-#  create figure #-#-#-#-#-#-#-#-#-#-#-#-#

P_ramp <- ggplot(ramp_long, aes(time, temperature))
P_ramp + geom_point(aes(colour = factor(codID)), alpha = 1/5) + 
	labs(x = "Time (s)", y = "Temperature (°C)") + theme(legend.position="bottom", legend.title = element_blank()) +
	guides(col = guide_legend(nrow = 3))




######################
# estimation statistics and graphs
######################


#-#-#-#-#-#-#-#-#-#-#-#-#  import data #-#-#-#-#-#-#-#-#-#-#-#-#

#load dataset with duplicate treatments for plotting purposes (see https://www.estimationstats.com/#/)
#duplicates are necessary in this case because "dabestr" doesn't allow all contrasts from a one-way anova with three factor levels
#this is a workaround to avoid having to produce two disctinct plots: a multi two-group (control vs ambient, control vs cooled) and an additional single two-group (ambient vs cooled)
#I opened an issue on the respr Github repo and Joses Ho mentioned that this might be addressed in a later version of the package

cool <- read.csv("Brain_cooling_data_for_fig.csv", header=T)
	dim(cool)
	head(cool)
	str(cool)
	summary(cool)

#select only the three groups with "_1" to remove duplicates and look at diagnostic plots
cool2 <- cool %>%
	filter(treatment == "control_1" | treatment =="ambient_1" | treatment == "cooled_1")
  hist(cool2$CTmax)
  qqp(cool2$CTmax, "norm")

# remove two statistical outliers from dataset with duplicates (to re-run stats and examine their influence)
cool3 <- cool[cool$CTmax > 24.7,]
	dim(cool3)
	summary(cool3)

# examine mean ans SD by group for Ctmax, mass and TL
cool2 %>% 
	group_by(treatment) %>% 
	summarise(CTmax_mean = mean(CTmax), CTmax_sd = sd(CTmax), TL_mean = mean(TL), 
	TL_sd = sd(TL), mass_mean = mean(mass), mass_sd = sd(mass))%>% 
	as.data.frame()

# examine mean ans SD by group for Ctmax without outliers
cool3 %>% 
	group_by(treatment) %>% 
	summarise(CTmax_mean = mean(CTmax), CTmax_sd = sd(CTmax)) %>% 
	as.data.frame()



#-#-#-#-#-#-#-#-#-#-#-#-# estimation stats for CTmax #-#-#-#-#-#-#-#-#-#-#-#-#

multi.two.group.CTmax <- 					#notice that subscripts _1 and _2 are used for this function to produced the desired comparisons; without subscripts the user gets an error message that factor levels are duplicated
  dabest(cool, treatment, CTmax, 
         idx = list(c("control_1", "ambient_1"),
         c("control_2", "cooled_1"),
	  c("ambient_2", "cooled_2")),
	  paired = FALSE, id.column = NULL, ci = 95,reps = 5000, func = mean, seed = 12345
         )
multi.two.group.CTmax #to look at the stats
plot(multi.two.group.CTmax) # estimation plot


#-#-#-#-#-#-#-#-#-#-#-#-# estimation stats for CTmax WITHOUT OUTLIERS #-#-#-#-#-#-#-#-#-#-#-#-#

multi.two.group.CTmax.NO <- 					#notice that subscripts _1 and _2 are used for this function to produced the desired comparisons; without subscripts the user gets an error message that factor levels are duplicated
  dabest(cool3, treatment, CTmax, 
         idx = list(c("control_1", "ambient_1"),
         c("control_2", "cooled_1"),
	  c("ambient_2", "cooled_2")),
	  paired = FALSE, id.column = NULL, ci = 95,reps = 5000, func = mean, seed = 12345
         )
multi.two.group.CTmax.NO #to look at the stats
plot(multi.two.group.CTmax.WO) # estimation plot


#-#-#-#-#-#-#-#-#-#-#-#-# estimation stats for TL #-#-#-#-#-#-#-#-#-#-#-#-#

multi.two.group.TL <- 					#notice that subscripts _1 and _2 are used for this function to produced the desired comparisons; without subscripts the user gets an error message that factor levels are duplicated
  dabest(cool, treatment, TL, 
         idx = list(c("control_1", "ambient_1"),
         c("control_2", "cooled_1"),
	  c("ambient_2", "cooled_2")),
	  paired = FALSE, id.column = NULL, ci = 95,reps = 5000, func = mean, seed = 12345
         )
multi.two.group.TL #to look at the stats
plot(multi.two.group.TL) # estimation plot


#-#-#-#-#-#-#-#-#-#-#-#-# estimation stats for mass #-#-#-#-#-#-#-#-#-#-#-#-#

multi.two.group.mass <- 					#notice that subscripts _1 and _2 are used for this function to produced the desired comparisons; without subscripts the user gets an error message that factor levels are duplicated
  dabest(cool, treatment, mass, 
         idx = list(c("control_1", "ambient_1"),
         c("control_2", "cooled_1"),
	  c("ambient_2", "cooled_2")),
	  paired = FALSE, id.column = NULL, ci = 95,reps = 5000, func = mean, seed = 12345
         )
multi.two.group.mass #to look at the stats
plot(multi.two.group.mass) # estimation plot