knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
tic <- Sys.time()
library("tidyverse")
library("cowplot")
library("car")
#library("SchulzLabTools")
library("officer")

#path.to.save.pptx <- "S:/Data_Daniel/electrical_synapses_voltage_modification/docs/figures.pptx"
#path.to.save.figs <- "S:/Data_Daniel/electrical_synapses_voltage_modification/figs"
#path.to.pptx.template <- "S:/Data_Daniel/electrical_synapses_voltage_modification/docs/blank.pptx"
path.to.pptx.template <- "../../docs/blank.pptx"
path.to.save.pptx <- "../../docs/figures.pptx"
path.to.save.figs <- "../../figs"

Read Data -------------------------------------------------------------------

TODO: cowplot::ggsave() doesn't have the right parameters and isn't saving correctly.

tecc <- readRDS(file = paste0(getwd(), "/data/tecc.rds"))
tevc <- readRDS(file = paste0(getwd(), "/data/tevc.rds"))
a.peak <- readRDS(file = paste0(getwd(), "/data/a.peak.rds"))
htk.peak <- readRDS(file = paste0(getwd(), "/data/htk.peak.rds"))
htk.ls.peak <- readRDS(file = paste0(getwd(), "/data/htk.ls.peak.rds"))
a.end <- readRDS(file = paste0(getwd(), "/data/a.end.rds"))
htk.end <- readRDS(file = paste0(getwd(), "/data/htk.end.rds"))
htk.ls.end <- readRDS(file = paste0(getwd(), "/data/htk.ls.end.rds"))
my_pres<-read_pptx(path.to.pptx.template)
#knitr::kable(layout_summary(my_pres))

Figures ---------------------------------------------------------------------

Method Validation ==========================================================

How are Input Resistance and Membrane Resistance Related?

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = `Input Resistance`, y = r1))+
  geom_abline(intercept = 0, slope = 1, linetype = "dotted")+
  geom_point(aes(color = Treatment))+
  geom_smooth(method = lm, se = FALSE, linetype = "dashed")+
  labs(x = "Input Resistance",  y = "Membrane Resistance")+
  theme(legend.position = "bottom")+
  xlim(0,20)+
  ylim(0,20)+
  geom_text(x = 10, y = 1, label = lm_eqn(m = lm(r1 ~ `Input Resistance`, data = tecc[tecc$`Time Exposed` == 0, ])
                                          ), parse = TRUE)+
  labs(title = "Membrane Resistance vs Input Resistance: All Data")

cowplot::ggsave("val.r1vr11.all.png", plot = ggplot2::last_plot(), 
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.3, 
                height = 6.2,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Are R11 and R1 Correlated?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.3, height = 6.2)

How much overlap is there between the input resistances of each group?

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = `Input Resistance`, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Input Resistance By Group")

cowplot::ggsave("val.r11.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Systematic Differences in R11?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.5, height = 6.2)


ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = `Input Resistance`, fill = Condition))+
  geom_histogram(bins = 10)+
  labs(title = "Input Resistance is Unimodal")

cowplot::ggsave("val.r11.hist.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 4.85,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Histogram of R11") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 4.85)

Membrane resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0 & tecc$r1 > 0, ], aes(x = Condition, y = r1, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Membrane Resistance By Group")

cowplot::ggsave("val.r1.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Systematic Differences in R1?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.5, height = 6)

Transfer resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = r12, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Transfer Resistance By Group")

cowplot::ggsave("val.r12.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Systematic Differences in R12?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.5, height = 6)

Coupling coefficient across conditions

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = `Coupling Coefficient`, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Coupling Coefficient By Group")

cowplot::ggsave("val.cc.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Systematic Differences in CC?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.5, height = 6)

Coupling resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0 & tecc$rc>0, ], aes(x = Condition, y = rc, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Coupling Resistance By Group")

cowplot::ggsave("val.rc.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Systematic Differences in RC?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.5, height = 6)

Is input resistance of one LC predictive of the other's?

rect <- tecc[tecc$`Time Exposed` == 0,c("Experiment", "Inj_Cell", "Condition", "r11","r12","r1","rc","Coupling Coefficient")]
inj.lc5 <- rect[rect$Inj_Cell == "LC5",]
inj.lc4 <- rect[rect$Inj_Cell == "LC4",]
names(inj.lc5) <- c("Experiment","Inj_Cell","Condition",
                    "LC5.r11","LC5->LC4.r12","LC5.r1",
                    "LC5->LC4.rc","LC5->LC4.Coupling Coefficient")
names(inj.lc4) <- c("Experiment","Inj_Cell","Condition",
                    "LC4.r11","LC4->LC5.r12","LC4.r1",
                    "LC4->LC5.rc","LC4->LC5.Coupling Coefficient")
rect <- cbind(inj.lc4, inj.lc5)
ggplot(rect, aes(x = LC4.r11, y = LC5.r11))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm, linetype = 2)+
  ylim(0,8)+
  xlim(0,8)+
  geom_text(x = 4, y = 0, label = lm_eqn(m = lm(LC5.r11 ~ LC4.r11, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC5 Input Resistance by LC4 Input Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r11.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Does LC4 R11 Predict LC5 R11?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 5.86)

Is there directionality in membrane resistance?

ggplot(rect, aes(x = LC4.r1, y = LC5.r1))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm, linetype = 2)+
  ylim(0,15)+
  xlim(0,15)+
  geom_text(x = 7.5, y = 0, label = lm_eqn(m = lm(LC5.r1 ~ LC4.r1, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC5 Membrane Resistance by LC4 Membrane Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r1.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Does LC4 R1 Predict LC5 R1?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 5.86)

Is there directionality in transfer resistance?

ggplot(rect, aes(x = `LC4->LC5.r12`, y = `LC5->LC4.r12`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm, linetype = 2)+
  ylim(0,7)+
  xlim(0,7)+
  geom_text(x = 3.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.r12`~`LC4->LC5.r12`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Transfer Resistance by \nLC5->LC4 Transfer Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r12.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Does Transfer Resistance Differ by Direction?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 5.86)

Is there directionality in coupling resistance?

ggplot(rect, aes(x = `LC4->LC5.rc`, y = `LC5->LC4.rc`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm, linetype = 2)+
  ylim(0,11)+
  xlim(0,11)+
  geom_text(x = 5.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.rc`~`LC4->LC5.rc`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Coupling Resistance by \nLC5->LC4 Coupling Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.rc.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Does RC Differ by Direction?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 5.86)

Is there directionality in coupling coefficient?

ggplot(rect, aes(x = `LC4->LC5.Coupling Coefficient`, y = `LC5->LC4.Coupling Coefficient`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm, linetype = 2)+
  ylim(0,1)+
  xlim(0,1)+
  geom_text(x = 0.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.Coupling Coefficient`~`LC4->LC5.Coupling Coefficient`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Coupling Coefficient \nby LC5->LC4 Coupling Coefficient")+
  theme(legend.position = "")

cowplot::ggsave("val.cc.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Does Coupling Coefficent Differ by Direction?") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 6.44, height = 5.86)

How close a online and offline leak subtracted measures of HTK?

merge_htk_methods <- function(input.df = htk.peak,
                              input.df.ls = htk.ls.peak){
  #input.df = htk.peak
  #input.df.ls = htk.ls.peak

  input.df <- input.df[, !(names(input.df) %in% c("Recording"))]
  input.df.ls <- input.df.ls[, !(names(input.df.ls) %in% c("Recording"))]

  #input.df$online <- FALSE
  #input.df.ls$online <- TRUE

  input.df$Experiment <- as.character(input.df$Experiment)
  input.df$IN4 <- as.character(input.df$IN4)
  input.df$IN9 <- as.character(input.df$IN9)

  input.df$target.mV <- ((round((input.df$mV)/10))*10)
  input.df.ls$target.mV <- ((round((input.df.ls$mV)/10))*10)

  input.df <- input.df %>% rename(manual.nA = "nA")
  input.df <- input.df %>% rename(manual.mV = "mV")

  input.df.ls <- input.df.ls %>% rename(online.mV = "mV")
  input.df.ls <- input.df.ls %>% rename(online.nA = "nA")

  #dplyr::full_join(input.df, input.df.ls) %>% head()

  #full_join(input.df, input.df.ls, by = c("Experiment", "Time Exposed", "Cell", "target.mV")) %>% head()

  return(dplyr::full_join(input.df, input.df.ls))
}

htk.comparison <- merge_htk_methods(input.df = htk.peak,
                                      input.df.ls = htk.ls.peak)

All the data: peak

peak.htk.mv <- ggplot(htk.comparison, aes(x = manual.mV, y = online.mV))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  geom_text(x = -30, y = -80, label = lm_eqn(m = lm(`online.mV`~`manual.mV`, data = htk.comparison)), parse = TRUE)+
  theme(legend.position = "")+
  ylim(-80,20)+
  xlim(-80,20)+
  labs(title = "No Systematic Difference in Voltage Control")

cowplot::ggsave("val.man.pn.htk.mv.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

peak.htk.na <- ggplot(htk.comparison, aes(x = manual.nA, y = online.nA))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  geom_text(x = 225, y = 0, label = lm_eqn(m = lm(`online.nA`~`manual.nA`, data = htk.comparison)), parse = TRUE)+
  ylim(0,450)+
  xlim(0,450)+
  theme(legend.position = "")+
  labs(title = "Measures of HTK Transient Differ Based on \nLeak Subtraction Method")

cowplot::ggsave("val.man.pn.htk.na.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = "Differences Between Online and Offline Leak Subtraction") %>%
  ph_with_gg(type = "body", value = peak.htk.mv, width = 6.44, height = 5.86, index=1) %>% 
  ph_with_gg(type = "body", value = peak.htk.na, width = 6.44, height = 5.86, index=2)

Broken down by target voltage: peak

voltages <- c(-80, -70, -60, -50, -40, -30, -20, -10,  0, 10, 20)
compare.htk.methods <- list()
for (i in 1:length(voltages)){
  compare.htk.methods[[i]] <- 
    ggplot(htk.comparison[htk.comparison$target.mV == voltages[i], ], aes(x = manual.nA, y = online.nA))+
    geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
    geom_point(aes(color = interact))+
    geom_smooth(method = lm, color = "red")+
    labs(title = paste("At", as.character(voltages[i]), "mV"))+
    theme(legend.position = "")

}
cowplot::plot_grid(plotlist = compare.htk.methods)

cowplot::ggsave("val.man.pn.htk.multi.peak.mv.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.72, 
                height = 6.23,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Correspondence of nA Measurements from -80mV to +20mV") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.72, height = 6.23)

Looking at the differences: peak

htk.comparison[,"nA.diff"] <- (htk.comparison[,"online.nA"] - htk.comparison[,"manual.nA"])

lower.sd <- (mean(na.omit(htk.comparison$nA.diff)) - sd(na.omit(htk.comparison$nA.diff)))
upper.sd <- (mean(na.omit(htk.comparison$nA.diff)) + sd(na.omit(htk.comparison$nA.diff)))
q1 <- as.numeric(summary(htk.comparison$nA.diff)[2])
q2 <- as.numeric(summary(htk.comparison$nA.diff)[3])
q3 <- as.numeric(summary(htk.comparison$nA.diff)[5])

htk.comparison$target.mV <- as.factor(htk.comparison$target.mV)
ggplot(htk.comparison, aes(x = nA.diff, fill = target.mV))+
  geom_histogram(bins = 100)+
  geom_segment(aes(x = lower.sd,
                   y = -15,
                   xend = upper.sd, 
                   yend = -15), 
                   color = "red",
                   size = 1
               )+
  geom_segment(aes(x = mean(na.omit(htk.comparison$nA.diff)),
                   y = -10,
                   xend = mean(na.omit(htk.comparison$nA.diff)), 
                   yend = -20), 
                   color = "red",
                   size = 1
               )+
    geom_segment(aes(x = q1,
                   y = -5,
                   xend = q3, 
                   yend = -5), 
                   color = "blue",
                   size = 1
               )+
  geom_segment(aes(x = q2,
                   y = 0,
                   xend = q2, 
                   yend = -10), 
                   color = "blue",
                   size = 1)+
  labs(title = "Online less Offline nA: \nQuartiles Denoted in Blue \nStandard Deviation in Red")+
  theme(legend.position = "right")

cowplot::ggsave("val.man.pn.htk.hist.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Difference Between Measurements") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)

All the data: end

htk.comparison <- merge_htk_methods(input.df = htk.end,
                                      input.df.ls = htk.ls.end)
end.htk.mv <- ggplot(htk.comparison, aes(x = manual.mV, y = online.mV))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  theme(legend.position = "")+
  ylim(-80,20)+
  xlim(-80,20)+
  labs(title = "No Systematic Difference in Voltage Control")+
  geom_text(x = -30, y = -75, label = lm_eqn(m = lm(`online.mV`~`manual.mV`, data = htk.comparison)), parse = TRUE)

cowplot::ggsave("val.man.pn.htk.mv.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

end.htk.na <- ggplot(htk.comparison, aes(x = manual.nA, y = online.nA))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  ylim(0,450)+
  xlim(0,450)+
  theme(legend.position = "")+
  labs(title = "Measures of HTK Transient Differ Based on \nLeak Subtraction Method")+
  geom_text(x = 225, y = 0, label = lm_eqn(m = lm(`online.nA`~`manual.nA`, data = htk.comparison)), parse = TRUE)

cowplot::ggsave("val.man.pn.htk.na.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs ,
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = "Differences Between Online and Offline Leak Subtraction") %>%
  ph_with_gg(type = "body", value = end.htk.mv, width = 6.44, height = 5.86, index=1) %>% 
  ph_with_gg(type = "body", value = end.htk.na, width = 6.44, height = 5.86, index=2)

Broken down by target voltage: end

voltages <- c(-80, -70, -60, -50, -40, -30, -20, -10,  0, 10, 20)
compare.htk.methods <- list()
for (i in 1:length(voltages)){
  compare.htk.methods[[i]] <- 
    ggplot(htk.comparison[htk.comparison$target.mV == voltages[i], ], aes(x = manual.nA, y = online.nA))+
    geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
    geom_point(aes(color = interact))+
    geom_smooth(method = lm, color = "red")+
    labs(title = paste("At", as.character(voltages[i]), "mV"))+
    theme(legend.position = "")
}
cowplot::plot_grid(plotlist = compare.htk.methods)

cowplot::ggsave("val.man.pn.htk.multi.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.72, 
                height = 6.23,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Correspondence of nA Measurements from -80mV to +20mV") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1  ,width = 9.72, height = 6.23)

Looking at the differences: end

htk.comparison[,"nA.diff"] <- (htk.comparison[,"online.nA"] - htk.comparison[,"manual.nA"])

lower.sd <- (mean(na.omit(htk.comparison$nA.diff)) - sd(na.omit(htk.comparison$nA.diff)))
upper.sd <- (mean(na.omit(htk.comparison$nA.diff)) + sd(na.omit(htk.comparison$nA.diff)))
q1 <- as.numeric(summary(htk.comparison$nA.diff)[2])
q2 <- as.numeric(summary(htk.comparison$nA.diff)[3])
q3 <- as.numeric(summary(htk.comparison$nA.diff)[5])

htk.comparison$target.mV <- as.factor(htk.comparison$target.mV)
ggplot(htk.comparison, aes(x = nA.diff, fill = target.mV))+
  geom_histogram(bins = 100)+
  geom_segment(aes(x = lower.sd,
                   y = -15,
                   xend = upper.sd, 
                   yend = -15), 
                   color = "red",
                   size = 1
               )+
  geom_segment(aes(x = mean(na.omit(htk.comparison$nA.diff)),
                   y = -10,
                   xend = mean(na.omit(htk.comparison$nA.diff)), 
                   yend = -20), 
                   color = "red",
                   size = 1
               )+
    geom_segment(aes(x = q1,
                   y = -5,
                   xend = q3, 
                   yend = -5), 
                   color = "blue",
                   size = 1
               )+
  geom_segment(aes(x = q2,
                   y = 0,
                   xend = q2, 
                   yend = -10), 
                   color = "blue",
                   size = 1)+
  labs(title = "Online less Offline nA: \nQuartiles Denoted in Blue \nStandard Deviation in Red")+
  theme(legend.position = "right")

cowplot::ggsave("val.man.pn.htk.hist.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Difference Between Measurements") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)

Write Out Slide Deck --------------------------------------------------------

print(my_pres, target = paste0(path.to.save.pptx, "method.validation.pptx") ) 

Visualize Treatments =======================================================

my_pres<-read_pptx(path.to.pptx.template)
add_ionic_slides_default_size <- function(title.prefix = "Control Wave in TEA: HTK",
                                          return.plts.peak = return.plts.peak,
                                          return.plts.end = return.plts.end){
  my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
    ph_with_text(type = "title", str = paste(title.prefix, "")) %>%
    ph_with_gg(type = "body", value = return.plts.peak[[1]], index=1) %>% 
    ph_with_gg(type = "body", value = return.plts.end[[1]], index=2)

  my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
    ph_with_text(type = "title", str = paste(title.prefix,"")) %>%
    ph_with_gg(type = "body", value = return.plts.peak[[2]], index=1) %>% 
    ph_with_gg(type = "body", value = return.plts.end[[2]], index=2)

  my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
    ph_with_text(type = "title", str = paste(title.prefix,"")) %>%
    ph_with_gg(type = "body", value = return.plts.peak[[3]], index=1) %>% 
    ph_with_gg(type = "body", value = return.plts.end[[3]], index=2)  
}

1. Phase Angle

All groups

TREATMENT <- "Phase Shift"

return.plts <- list()
# Un-normalized ---------------------------------------------------------------
return.plts[[1]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = (rc),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "rc",
    title = paste(TREATMENT, "Coupling Resistance")
  ) +

  facet_grid(~Condition)

return.plts[[2]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "CC",
    title = paste(TREATMENT, "Coupling Coefficient")
  ) +

  facet_grid(~Condition)

return.plts[[3]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Input Resistance",
    title = paste(TREATMENT, "Input Resistance")
  ) +

  facet_grid(~Condition)

# Normalized ------------------------------------------------------------------
return.plts[[4]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = (rc / rc.0),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized rc",
    title = paste(TREATMENT, "Normalized Coupling Resistance")
  ) +

  facet_grid(~Condition)

return.plts[[5]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized CC",
    title = paste(TREATMENT, "Normalized Coupling Coefficient")
  ) +

  facet_grid(~Condition)

return.plts[[6]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance` / `Input Resistance.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized Input Resistance",
    title = paste(TREATMENT, "Normalized Input Resistance")
  ) +

  facet_grid(~Condition)

# Group interaction style plots -----------------------------------------------

return.plts[[7]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition, group = Condition
  )) +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  stat_summary(fun.y = mean, geom = "point", pch = 12, size = 5) +
  stat_summary(fun.y = mean, geom = "line", size = 1) +
  labs(color = "Phase \nAngle") +
  labs(title = "Interaction of Treatment:Time \non Coupling Coefficient (normalized)", x = "Time in Minutes", y = "Coupling Coefficient") +

  ylim(0, 2)

return.plts[[8]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition, group = interaction(Condition, `Time Exposed`)
  )) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  labs(title = "Coupling Coefficients over Time", x = "Time in Minutes", y = "Coupling Coefficient") +
  facet_grid(~Condition) +

  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "") +
  ylim(0, 1)

# -----------------------------------------------------------------------------
save_fig_list(input.list = return.plts,
                          title.prefix = "phaseshift",
                          file.type = ".png",
                          save.path = path.to.save.figs)


for (i in 1:length(return.plts)){
  print(i)
  my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
    ph_with_text(type = "title", str = "TITLE") %>%
    ph_with_gg(type = "body", value = return.plts[[i]], index=1)
}
ggplot(na.omit(tevc[tevc$Treatment == "Phase Shift" & tevc$`Time Exposed` < 90,]), aes(x = `Time Exposed`, y = slopes, color = Condition, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")+
  theme(legend.position = "")+
  facet_grid(~Condition)

cowplot::ggsave("phaseshift.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Phase Shift: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)

2. Small Phase Angles

#show_tecc_plots(TREATMENT = "Small Phase Angle")
TREATMENT <- "Small Phase Angle"

return.plts <- list()
# Un-normalized ---------------------------------------------------------------
return.plts[[1]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = (rc),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "rc",
    title = paste(TREATMENT, "Coupling Resistance")
  )

return.plts[[2]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "CC",
    title = paste(TREATMENT, "Coupling Coefficient")
  )

return.plts[[3]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Input Resistance",
    title = paste(TREATMENT, "Input Resistance")
  )

# Normalized ------------------------------------------------------------------
return.plts[[4]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = (rc / rc.0),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized rc",
    title = paste(TREATMENT, "Normalized Coupling Resistance")
  ) 

return.plts[[5]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized CC",
    title = paste(TREATMENT, "Normalized Coupling Coefficient")
  )

return.plts[[6]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance` / `Input Resistance.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized Input Resistance",
    title = paste(TREATMENT, "Normalized Input Resistance")
  ) 

# Group interaction style plots -----------------------------------------------

return.plts[[7]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = interact, group = Condition
  )) +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  stat_summary(fun.y = mean, geom = "point", pch = 12, size = 5) +
  stat_summary(fun.y = mean, geom = "line", size = 1) +
  labs(color = "Phase \nAngle") +
  labs(title = "Interaction of Treatment:Time \non Coupling Coefficient (normalized)", x = "Time in Minutes", y = "Coupling Coefficient") +

  ylim(0, 2)

return.plts[[8]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = Condition,
    y = `Coupling Coefficient`,
    color = interact, group = Condition
  )) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  labs(title = "Coupling Coefficients over Time", x = "Time in Minutes", y = "Coupling Coefficient") +
  facet_grid(~`Time Exposed`) +

  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "") +
  ylim(0, 1)
# -----------------------------------------------------------------------------
save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase",
                          file.type = ".png",
                          save.path = path.to.save.figs)

title.prefix <- "Small Phase Angle:"
my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix, "Coupling Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[1]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[4]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Coupling Coefficient")) %>%
  ph_with_gg(type = "body", value = return.plts[[2]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[5]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Input Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[3]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[6]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[7]], index=1)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[8]], index=1)
ggplot(tevc[tevc$Treatment == "Small Phase Angle",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("smallphase.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs,
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Small Phase Angle: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)
return.plts.peak <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "smallphase.a.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "smallphase.end.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Small Phase Angle: A",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

return.plts.peak <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "smallphase.htk.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "smallphase.htk.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Small Phase Angle: HTK",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

3. Inverted Wave

inv.plts <- show_tecc_plots(TREATMENT = "Inverted Wave")
save_fig_list(input.list = inv.plts,
                          title.prefix = "inv",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts <- inv.plts
title.prefix <- "Inverted Wave:"

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix, "Coupling Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[1]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[4]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Coupling Coefficient")) %>%
  ph_with_gg(type = "body", value = return.plts[[2]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[5]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Input Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[3]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[6]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[7]], index=1)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[8]], index=1)
ggplot(tevc[tevc$Treatment == "Inverted Wave",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("inv.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Inverted Wave: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)
return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "inv.a.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "inv.a.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Inverted Wave: A",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "inv.htk.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "inv.htk.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Inverted Wave: HTK",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

4. Aberrant Wave

ab.plts <- show_tecc_plots(TREATMENT = "Aberrant Wave")
save_fig_list(input.list = ab.plts,
                          title.prefix = "ab",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts <- ab.plts
title.prefix <- "Aberrant Wave:"

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix, "Coupling Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[1]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[4]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Coupling Coefficient")) %>%
  ph_with_gg(type = "body", value = return.plts[[2]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[5]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Input Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[3]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[6]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[7]], index=1)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[8]], index=1)
ggplot(tevc[tevc$Treatment == "Aberrant Wave",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("ab.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Aberrant Wave: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)
return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "ab.a.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "ab.a.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Aberrant Wave: A",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "ab.htk.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "ab.htk.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Aberrant Wave: HTK",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

5. Silent

sh.plts <- show_tecc_plots(TREATMENT = "Silent Control")
save_fig_list(input.list = sh.plts,
                          title.prefix = "sh",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts <- sh.plts
title.prefix <- "Silent Control:"

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix, "Coupling Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[1]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[4]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Coupling Coefficient")) %>%
  ph_with_gg(type = "body", value = return.plts[[2]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[5]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Input Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[3]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[6]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[7]], index=1)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[8]], index=1)
ggplot(tevc[tevc$Treatment == "Silent Control",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("sh.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Silent Control: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)
return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "sh.a.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "sh.a.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Silent Control: A",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "sh.htk.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "sh.htk.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Silent Control: HTK",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

6. Control in TEA

ctea.plts <- show_tecc_plots(TREATMENT = "Control Wave + TEA")
save_fig_list(input.list = ctea.plts,
                          title.prefix = "ctea",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts <- ctea.plts
title.prefix <- "Control Wave in TEA:"

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix, "Coupling Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[1]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[4]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Coupling Coefficient")) %>%
  ph_with_gg(type = "body", value = return.plts[[2]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[5]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
ph_with_text(type = "title", str = paste(title.prefix,"Input Resistance")) %>%
  ph_with_gg(type = "body", value = return.plts[[3]], index=1) %>% 
  ph_with_gg(type = "body", value = return.plts[[6]], index=2)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[7]], index=1)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = paste(title.prefix, "Coupling Coefficients Over Time")) %>%
  ph_with_gg(type = "body", value = return.plts[[8]], index=1)
ggplot(tevc[tevc$Treatment == "Control Wave + TEA",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("ctea.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = path.to.save.figs, 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
  ph_with_text(type = "title", str = "Control Wave in TEA: Voltage Clamp Current Fits") %>%
  ph_with_gg(type = "body", value = ggplot2::last_plot(), index=1)
return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "ctea.a.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "ctea.a.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Control Wave in TEA: A",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)

## ============================================================================

return.plts.peak <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts.peak,
                          title.prefix = "ctea.htk.peak",
                          file.type = ".png",
                          save.path = path.to.save.figs)

return.plts.end <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts.end,
                          title.prefix = "ctea.htk.end",
                          file.type = ".png",
                          save.path = path.to.save.figs)

add_ionic_slides_default_size(title.prefix = "Control Wave in TEA: HTK",
                              return.plts.peak = return.plts.peak,
                              return.plts.end = return.plts.end)
print(my_pres, target = paste0(path.to.save.pptx, "data.vis.pptx") ) 

Perform Statistical Tests ---------------------------------------------------

my_pres<-read_pptx(path.to.pptx.template)
tecc.r <- tecc[,c("Experiment", "Time Exposed", "Condition", "Treatment", "Inj_Cell", 
                  "r11", "r12", "r1", "rc", "Coupling Coefficient")]
tecc.r <- rename(tecc.r, cc = `Coupling Coefficient`)

tecc.r <- rename(tecc.r, Cell = Inj_Cell)
tevc.r <- tevc[, c("Experiment", "Time Exposed", "Condition", "Treatment", "Inj_Cell", "slopes")]
tevc.r <- rename(tevc.r, Cell = Inj_Cell)
prep_ioinc_for_bfdf <- function(input.df = a.peak,
                                closest.to = 0,
                                rename.mV = "a.peak.0.mV",
                                rename.nA = "a.peak.0.nA"){
  #a <- seq(from = -9.5, to = 9, by = 2)
  #closest.to = 5
  #min(abs(a - closest.to))
  #a[which(abs(a - closest.to) == min(abs(a - closest.to)))]
  #input.df = a.peak
  #closest.to = 0
  #rename.mV = "a.peak.0.mV"
  #rename.nA = "a.peak.0.nA"
  #INDIVIDUAL = 1
  #TIME = 1

  input.df$keep.me <- FALSE


  INDIVIDUALS <- unique(input.df$interact)
  for (INDIVIDUAL in 1:length(INDIVIDUALS)){
    TIMES <- unique(input.df[input.df$interact == INDIVIDUALS[INDIVIDUAL] ,"Time Exposed"])
    for (TIME in 1:length(TIMES)){
      print(paste(as.character(INDIVIDUAL), as.character(TIME)))

      #This is one of the most opaque lines I've written. It finds the row containing the voltage which is closest to the target voltage.
      input.df[input.df$interact == INDIVIDUALS[INDIVIDUAL] & input.df$`Time Exposed` == TIMES[TIME], ][which(abs(input.df[input.df$interact == INDIVIDUALS[INDIVIDUAL] & input.df$`Time Exposed` == TIMES[TIME] , "mV"] - closest.to) == min(abs(input.df[input.df$interact == INDIVIDUALS[INDIVIDUAL] & input.df$`Time Exposed` == TIMES[TIME] , "mV"] - closest.to))), "keep.me"] <- TRUE
    }
  }
  input.df <- input.df[input.df$keep.me == TRUE, ]
  input.df <- input.df[, c("Experiment", "mV", "nA", "Time Exposed", "Condition", "Treatment", "Cell")]
  names(input.df) <-  c("Experiment", rename.mV, rename.nA, "Time Exposed", "Condition", "Treatment", "Cell")
  return(input.df)
}

mk_pre_bfdf_factors_chr <- function(input.df = a.peak.r){
  input.df$Experiment <- as.character(input.df$Experiment)
  input.df$Condition <- as.character(input.df$Condition)
  input.df$Treatment <- as.character(input.df$Treatment)
  input.df$Cell <- as.character(input.df$Cell)
  return(input.df)
}

a.peak.r <- prep_ioinc_for_bfdf(input.df = a.peak,
                                closest.to = 0,
                                rename.mV = "a.peak.0.mV",
                                rename.nA = "a.peak.0.nA")
a.end.r <- prep_ioinc_for_bfdf(input.df = a.end,
                                closest.to = 0,
                                rename.mV = "a.end.0.mV",
                                rename.nA = "a.end.0.nA")
htk.peak.r <- prep_ioinc_for_bfdf(input.df = htk.peak,
                                closest.to = 0,
                                rename.mV = "htk.peak.0.mV",
                                rename.nA = "htk.peak.0.nA")
htk.end.r <- prep_ioinc_for_bfdf(input.df = htk.end,
                                closest.to = 0,
                                rename.mV = "htk.end.0.mV",
                                rename.nA = "htk.end.0.nA")

tecc.r <- mk_pre_bfdf_factors_chr(input.df = tecc.r)
tevc.r <- mk_pre_bfdf_factors_chr(input.df = tevc.r)
a.peak.r <- mk_pre_bfdf_factors_chr(input.df = a.peak.r)
a.end.r <- mk_pre_bfdf_factors_chr(input.df = a.end.r)
htk.peak.r <- mk_pre_bfdf_factors_chr(input.df = htk.peak.r)
htk.end.r <- mk_pre_bfdf_factors_chr(input.df = htk.end.r)

#make bfdf

bfdf <- full_join(tecc.r, tevc.r) %>% full_join(a.peak.r) %>% full_join(a.end.r) %>% full_join(htk.peak.r) %>% full_join(htk.end.r)

bfdf <- bfdf[!is.na(bfdf$Experiment),]
summary(bfdf)
#write.csv(bfdf, file = "C:/Users/dan/Desktop/bfdf.csv")
conditions <- unique(bfdf$Condition)
i = 10
temp <- bfdf[bfdf$Condition == conditions[i], c("Experiment", "Time Exposed", "cc")]

temp <- tecc.r[tecc.r$Condition == conditions[i], c("Experiment", "Cell", "Time Exposed", "cc")]
temp <- temp[!duplicated((temp)),]

#temp <- temp %>% spread(`Time Exposed`, cc)
#t.test(temp$`0`, temp$`60`, paired = TRUE)
temp$Individual <- interaction(temp$Experiment, temp$Cell)

# Testing Assumptions ---------------------------------------------------------
#independence
#univariate normality
car::qqPlot(temp[temp$`Time Exposed` == 0, "cc"])
shapiro.test(temp[temp$`Time Exposed` == 0, "cc"]) # if p<0.05, data is not normal
#homogeneity of variances (homoscedasticity)
bartlett.test(cc ~ `Time Exposed`, data = temp) # if p<0.05 variances are not homogeneous
car::leveneTest(cc ~ as.factor(`Time Exposed`), data = temp)
fligner.test(cc ~ as.factor(`Time Exposed`), data = temp) #non parametric version

fm <- aov(cc ~ `Time Exposed` + Error(Individual), data = temp)
summary(fm)
#NO IDEA WHY: this function if run will not read in the last curly brace. That is why it's here instead on in the code block below. when this chunk is run as a chunk it works.
mk_slides_for_stats <- function(input.df = bfdf[bfdf$Condition == "TEA_ControlWave", ],
                                dependent.variables = c(
                                  "r11", "r1", "rc", "cc",
                                  "slopes",
                                  "a.peak.0.nA", "a.end.0.nA",
                                  "htk.peak.0.nA", "htk.end.0.nA"
                                ),
                                condition.name = "TEA_ControlWave",
                                t1 = 0,
                                t2 = 60,
                                iterations = 1000,
                                use.seed = 2834959) {
  debugging <- FALSE
  if (debugging == TRUE) {
    input.df <- bfdf[bfdf$Condition == "TEA_ControlWave", ]
    dependent.variables <- c(
      "r11", "r1", "rc", "cc",
      "slopes",
      "a.peak.0.nA", "a.end.0.nA",
      "htk.peak.0.nA", "htk.end.0.nA"
    )
    condition.name <- "TEA_ControlWave"
    t1 <- 0
    t2 <- 60
    iterations <- 10
    use.seed <- 2834959
  }

  set.seed(use.seed)
  my_pres <- add_slide(x = my_pres, layout = "Section Header", master = "Office Theme") %>%
    ph_with_text(type = "title", str = as.character(condition.name))
  print("added section header")

  ## 1 select data ==============================================================
  temp <- input.df
  temp[, "Cell_Id"] <- interaction(temp[, "Experiment"], temp[, "Cell"])
  temp <- temp[, !(names(temp) %in% c("Experiment", "Condition", "Treatment", "Cell"))]
  temp <- temp[!duplicated(temp), ]

  Cells_at_t1 <- unique(temp[temp$`Time Exposed` == t1, "Cell_Id"])
  Cells_at_t2 <- unique(temp[temp$`Time Exposed` == t2, "Cell_Id"])
  keep_cells <- intersect(Cells_at_t1, Cells_at_t2)
  temp <- temp[(temp$`Time Exposed` == t1 | temp$`Time Exposed` == t2) &
    temp$Cell_Id %in% keep_cells, ]
  temp.backup <- temp
  for (i in 1:length(dependent.variables)) {
    print(i)
    temp <- temp.backup
    temp <- temp[, c("Cell_Id", "Time Exposed", dependent.variables[i])]
    temp <- spread(temp, `Time Exposed`, 3)
    temp <- temp[!(is.na(temp[2]) | is.na(temp[3])), ]
    ### outliers?
    temp$t1.has.outlier <- FALSE
    temp$t2.has.outlier <- FALSE
    temp[(temp[, 2] < (as.numeric(quantile(temp[, 2])[2]) - (1.5 * as.numeric(IQR(temp[, 2]))))), "t1.has.outlier"] <- TRUE
    temp[(temp[, 2] > (as.numeric(quantile(temp[, 2])[4]) + (1.5 * as.numeric(IQR(temp[, 2]))))), "t1.has.outlier"] <- TRUE
    temp[(temp[, 3] < (as.numeric(quantile(temp[, 3])[2]) - (1.5 * as.numeric(IQR(temp[, 3]))))), "t2.has.outlier"] <- TRUE
    temp[(temp[, 3] > (as.numeric(quantile(temp[, 3])[4]) + (1.5 * as.numeric(IQR(temp[, 3]))))), "t2.has.outlier"] <- TRUE

    if ((sum(temp$t1.has.outlier) > 0) | (sum(temp$t2.has.outlier) > 0)) {
      warning("Data contains outliers more extreme than q+1.5IQR")
      print(temp[temp$t1.has.outlier == TRUE | temp$t2.has.outlier == TRUE, ])
    }
    outliers <- temp[temp$t1.has.outlier == TRUE | temp$t2.has.outlier == TRUE, ]
    temp <- temp[temp$t1.has.outlier != TRUE & temp$t2.has.outlier != TRUE, 1:3]
    ## 2 test assumptions =========================================================
    ### univariate normality
    plt.qq <- ggplot(temp, aes(sample = temp[, 2])) + stat_qq()
    # car::qqPlot(temp[, 2])
    check.norm <- shapiro.test(temp[, 2]) # if p<0.05, data is not normal

    ## 3 run test =================================================================
    fm <- t.test(temp[, 2],
      temp[, 3],
      paired = TRUE
    )

    ## 4 get empirical p values ===================================================
    statistics <- array()
    statistics <- c(as.numeric(fm$statistic))

    all.values <- c(temp[, 2], temp[, 3])
    fake.temp <- temp
    for (iteration in 2:iterations) {
      mixed.values <- sample(all.values)
        fake.temp[, 3], paired = TRUE
      ))$statistic))
        fake.temp[, 2],
      statistics <- c(statistics, as.numeric((t.test(
      fake.temp[, 3] <- mixed.values[(nrow(fake.temp) + 1):length(mixed.values)]
      fake.temp[, 2] <- mixed.values[1:nrow(fake.temp)]
    }
    statistics <- as.data.frame(statistics)
    names(statistics) <- c("t.stat")
    statistics$more.extreme <- FALSE
    statistics[abs(statistics$t.stat) >= abs(statistics[1, "t.stat"]), "more.extreme"] <- TRUE
    plt.emp <- ggplot(statistics, aes(x = t.stat, fill = more.extreme)) + geom_histogram(bins = (iterations / 4)) + geom_vline(xintercept = as.numeric(statistics[1, "t.stat"]), linetype = 2) + labs(title = paste("Empirical P :", as.character(mean(statistics[, "more.extreme"]))))

    ## 5 add results to slides ====================================================
    #Added because this function is missing the last }
    my_pres <- add_slide(x = my_pres, layout = "Title and Content", master = "Office Theme") %>%
      ph_with_text(type = "title", str = paste(as.character(dependent.variables[i]), ":", "Omitted Outliers")) %>%
      ph_with_table(type = "body", value = outliers)

    my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
      ph_with_text(type = "title", str = paste(as.character(dependent.variables[i]), ":", "Checking for Normality")) %>%
      ph_with_text(type = "body", str = paste(
        "Shapiro-Wilk normality test: W=",
        as.character(check.norm$statistic),
        "p-value",
        as.character(check.norm$p.value)
      ), index = 1) %>%
      ph_with_gg(type = "body", value = plt.qq, index = 2)

    my_pres <- add_slide(x = my_pres, layout = "Two Content", master = "Office Theme") %>%
      ph_with_text(type = "title", str = paste(as.character(dependent.variables[i]), ":", "Paired t-test")) %>%
      ph_with_text(type = "body", str = paste(
        "Paired t-test: t=",
        as.character(fm$statistic),
        "df=",
        as.character(fm$parameter),
        "p-value",
        as.character(fm$p.value)
      ), index = 1) %>%
      ph_with_gg(type = "body", value = plt.emp, index = 2)

    print("slides added")
    }
  }
# [1] "Phase Offset 0 Degrees"   "Phase Offset 180 Degrees" "Phase Offset 90 Degrees" 
# [4] "Phase Offset 45 Degrees"  "Inverted Wave"            "TEA_ControlWave"         
# [7] "TEA_Silent"               "Aberrant Depol"           "hold_at_-53mV"           
#[10] "22.5_degrees"             "11.2_degrees"   
run.this.many.iterations <- 1000

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Phase Offset 0 Degrees", ],
                    dependent.variables = c(
                    "r11", "r1", "rc", "cc",
                    "slopes"
                    ),
                    condition.name = "Phase Offset 0 Degrees",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Phase Offset 45 Degrees", ],
                    dependent.variables = c(
                    "r11", "r1", "rc", "cc",
                    "slopes"
                    ),
                    condition.name = "Phase Offset 45 Degrees",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Phase Offset 90 Degrees", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes"),
                    condition.name = "Phase Offset 90 Degrees",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Phase Offset 180 Degrees", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes"),
                    condition.name = "Phase Offset 180 Degrees",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "22.5_degrees", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "22.5_degrees",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

if(FALSE == TRUE){
  #TODO this needs to accept TWO conditions (so we can use 80 minutes as baseline). 
  #Need to retool fcn to make this happen. 

  mod.bfdf <- bfdf[(bfdf$Condition == "22.5_degrees" & bfdf$`Time Exposed` == 80) | bfdf$Condition == "11.2_degrees",] #Mk disposable df
  mod.bfdf[mod.bfdf$`Time Exposed` == 81, "Time Exposed"] <- 80
  mod.bfdf <- mod.bfdf[!duplicated(mod.bfdf),]
  mod.bfdf$keep <- TRUE
  mod.bfdf[c(17,18),"keep"] <- FALSE
  mod.bfdf <- mod.bfdf[mod.bfdf$keep == TRUE,]
  mod.bfdf <- mod.bfdf[, !(names(mod.bfdf) %in% c("keep"))] #Without manual deduplication here the funciton is breaking because there are four rows with duplicated values. Error: Duplicate identifiers for rows (15, 17), (16, 18) 

mk_slides_for_stats(input.df = mod.bfdf,
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "11.2_degrees",
                    t1 = 80,
                    t2 = 140,
                    iterations = 10,
                    use.seed = 2834959)
  }

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Inverted Wave", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "Inverted Wave",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "Aberrant Depol", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "Aberrant Depol",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "hold_at_-53mV", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "hold_at_-53mV",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)

mk_slides_for_stats(input.df = bfdf[bfdf$Condition == "TEA_ControlWave", ],
                    dependent.variables = c("r11", "r1", "rc", "cc",
                                            "slopes", 
                                            "a.peak.0.nA", "a.end.0.nA", 
                                            "htk.peak.0.nA", "htk.end.0.nA"),
                    condition.name = "TEA_ControlWave",
                    t1 = 0,
                    t2 = 60,
                    iterations = run.this.many.iterations,
                    use.seed = 2834959)
## 1 select data ==============================================================
temp <- input.df[input.df$Condition == target.condition,]
temp[,"Cell_Id"] <- interaction(temp[,"Experiment"], temp[,"Cell"])
temp <- temp[, !(names(temp) %in% c("Experiment", "Condition", "Treatment", "Cell"))]
temp <- temp[!duplicated(temp),]

Cells_at_t1 <- unique(temp[temp$`Time Exposed` == t1 , "Cell_Id"])
Cells_at_t2 <- unique(temp[temp$`Time Exposed` == t2 , "Cell_Id"])
keep_cells <- intersect(Cells_at_t1, Cells_at_t2)
temp <- temp[(temp$`Time Exposed` == t1 | temp$`Time Exposed` == t2) & 
               temp$Cell_Id %in% keep_cells, ]
###outliers?

###make sure lengths are the same at each time point

## 2 test assumptions =========================================================
###independence
###univariate normality
car::qqPlot(temp[temp$`Time Exposed` == 0, dependent.variables[i]])
shapiro.test(temp[temp$`Time Exposed` == 0, dependent.variables[i]]) # if p<0.05, data is not normal
###homogeneity of variances (homoscedasticity)
bartlett.test(dependent.variables[i] ~ `Time Exposed`, data = temp) # if p<0.05 variances are not homogeneous
car::leveneTest(dependent.variables[i] ~ as.factor(`Time Exposed`), data = temp)
fligner.test(cc ~ as.factor(`Time Exposed`), data = temp) #non parametric version

## 3 run parametric or non parametric test ====================================

aov(dependent.variables[i] ~ `Time Exposed` + Error(Cell_Id), data = temp) %>% summary()


## 4 get empirical p values ===================================================
print(my_pres, target = paste0(path.to.save.pptx, "stat.tests.pptx") ) 


danielkick/esynvmod documentation built on May 17, 2019, 7:02 p.m.