library(Cairo) #for cairo graphics to change font, save ggplots as pdfs etc
library(brms) #for fitting models in STAN
library(coda) #for diagnostics
library(bayesplot) #for diagnostics
library(rstan) #for interfacing with STAN
library(DHARMa) #for residual diagnostics
library(emmeans) #for marginal means etc
library(broom) #for tidying outputs
library(broom.mixed)
library(car) #for scatterplot matrix
library(tidybayes) #for more tidying outputs
library(ggeffects) #for partial plots
library(tidyverse) #for data wrangling etc
library(patchwork) #for combining graphs
library(ggplot2) #for plots
library(extrafont) #for extra text font families
library(grid) #for graphics
library(GGally) #for scatterplot matrices
library(knitr) #for kableGabazinedata = read_csv('01_Data/Gabazinedata.csv', trim_ws = TRUE) #trim white space
#remove senescent squid (experimental notes for each of these squid that they were close to death by senescence and thus behavioural results not reliable), rows 11, 13, 14 and 37 = squid 650, 652, 654 and 682
Gabazinedata_old <- Gabazinedata [c(-11,-13,-14,-37),]
save(Gabazinedata_old, file = '01_Data/Gabazinedata_old.RData')
#select just the variables interested in, melt data frame into long format (make one column called response which has my response variables of interest stacked on top of each other) and remove NAs
Gabazinedata_old_long_p = Gabazinedata_old %>%
select('Squid', 'CO2', 'Drug', 'Treatment', 'Acclimation_days', 'Tank_intro_day','System', 'Test_time_day_hours', 'Treatment_time_hours', 'Drug_test_no.', 'Behavioural_tank', 'ML_cm', "No._visits_A", "Soft_touch_no._no_0", "Aggressive_touch_no._no_0","Time_A_s_no_0", "Latency_soft_touch_s", "Latency_aggressive_touch_s", "Active_time_s", "Dist", "Vel", "Soft_touch_binomial", "Aggressive_touch_binomial") %>%
pivot_longer(cols = c("No._visits_A", "Soft_touch_no._no_0", "Aggressive_touch_no._no_0","Time_A_s_no_0", "Latency_soft_touch_s", "Latency_aggressive_touch_s", "Active_time_s", "Dist", "Vel", "Soft_touch_binomial", "Aggressive_touch_binomial"),
names_to = 'Response',
values_to = 'Response_measurement',
values_drop_na = TRUE) %>%
mutate(CO2 = factor(CO2, levels=c('Ambient', 'Elevated')),
Drug = factor(Drug, levels=c('Sham', 'Gabazine')),
Treatment = factor(Treatment, levels=c('Ambient_Sham', 'Ambient_Gabazine', 'Elevated_Sham', 'Elevated_Gabazine')),
fBehavioural_tank = factor(Behavioural_tank),
fDrug_test_no. = factor(Drug_test_no.),
fAcclimation_days = factor(Acclimation_days),
fSystem = factor(System),
fTank_intro_day = factor(Tank_intro_day),
Response = factor(Response))
str(Gabazinedata_old_long_p)
save(Gabazinedata_old_long_p, file = '01_Data/Gabazinedata_old_long_p.RData')
#Nest so one row for each response variable
Gabazinedata_old_response_p = Gabazinedata_old_long_p %>%
group_by(Response) %>%
nest()
#Update so datatype correct for each response var
datatype_function <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
data$data[[1]] %>%
mutate(Response_measurement = as.numeric(Response_measurement))
} else if (grepl('o.', response_var)) {
data$data[[1]] %>%
mutate(Response_measurement = as.integer(Response_measurement))
} else {
data$data[[1]] %>%
mutate(Response_measurement = as.numeric(Response_measurement))
}
}
Gabazinedata_old_response_p <- Gabazinedata_old_response_p %>%
mutate(data = map(Response, datatype_function))
save(Gabazinedata_old_response_p, file = '01_Data/Gabazinedata_old_response_p.RData')Picrotoxindata = read_csv('01_Data/Picrotoxindata.csv', trim_ws = TRUE) #trim white space
#Remove squid 516 (row 15) as has ML_cm missing - effects model selection
Picrotoxindata = Picrotoxindata [-15,]
save(Picrotoxindata, file = '01_Data/Picrotoxindata.RData')
##select just the variables interested in, melt data frame into long format (make one column called response which has my response variables of interest stacked on top of each other) and remove NAs
Picrotoxindata_long_p = Picrotoxindata %>%
select('Squid', 'CO2', 'Drug', 'Treatment', 'Acclimation_days', 'Tank_intro_day','System', 'Test_time_day_hours', 'Treatment_time_hours', 'Drug_test_no._dif_drugs', 'Behavioural_tank', 'ML_cm', "No._visits_A", "Soft_touch_no._no_0", "Aggressive_touch_no._no_0", "Time_A_s_no_0", "Latency_soft_touch_s", "Latency_aggressive_touch_s", "Active_time_s", "Dist", "Vel", "Soft_touch_binomial", "Aggressive_touch_binomial") %>%
pivot_longer(cols = c("No._visits_A", "Soft_touch_no._no_0", "Aggressive_touch_no._no_0", "Time_A_s_no_0", "Latency_soft_touch_s", "Latency_aggressive_touch_s", "Active_time_s", "Dist", "Vel", "Soft_touch_binomial", "Aggressive_touch_binomial"),
names_to = 'Response',
values_to = 'Response_measurement',
values_drop_na = TRUE) %>%
mutate(CO2 = factor(CO2, levels=c('Ambient', 'Elevated')),
Drug = factor(Drug, levels=c('Sham', 'Picrotoxin')),
Treatment = factor(Treatment, levels=c('Ambient_Sham', 'Ambient_Picrotoxin', 'Elevated_Sham', 'Elevated_Picrotoxin')),
fBehavioural_tank = factor(Behavioural_tank),
fDrug_test_no._dif_drugs = factor(Drug_test_no._dif_drugs),
fAcclimation_days = factor(Acclimation_days),
fSystem = factor(System),
fTank_intro_day = factor(Tank_intro_day),
Response = factor(Response))
str(Picrotoxindata_long_p)
save(Picrotoxindata_long_p, file = '01_Data/Picrotoxindata_long_p.RData')
#Nest so one row for each response variable
Picrotoxindata_response_p = Picrotoxindata_long_p %>%
group_by(Response) %>%
nest()
#Update so datatype correct for each response var
datatype_function <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
data$data[[1]] %>%
mutate(Response_measurement = as.numeric(Response_measurement))
} else if (grepl('o.', response_var)) {
data$data[[1]] %>%
mutate(Response_measurement = as.integer(Response_measurement))
} else {
data$data[[1]] %>%
mutate(Response_measurement = as.numeric(Response_measurement))
}
}
Picrotoxindata_response_p <- Picrotoxindata_response_p %>%
mutate(data = map(Response, datatype_function))
save(Picrotoxindata_response_p, file = '01_Data/Picrotoxindata_response_p.RData')load(file = '01_Data/Gabazinedata_old_response_p.RData')load(file = '01_Data/Picrotoxindata_response_p.RData')boxplots_function <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
ggplot(data = data$data[[1]], aes(x = Treatment, y = Response_measurement)) +
geom_boxplot() +
scale_y_continuous(response)
}
Gabazinedata_old_response_exploratory_p <- Gabazinedata_old_response_p %>%
mutate(boxplots = map(Response, boxplots_function))
pdf("02_Exploratoryanalysis_outputs/Gabazine_boxplots.pdf")
Gabazinedata_old_response_exploratory_p$boxplots
dev.off()stacked_function <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
data_plot <- data$data[[1]] %>%
mutate(Response_measurement = as.logical(Response_measurement))
plot_stacked <- ggplot(data = data_plot, aes(x = Treatment, fill = Response_measurement)) +
geom_bar(position = 'fill') +
scale_y_continuous(response_var)
} else {
}
}
Gabazinedata_old_response_exploratory_p <- Gabazinedata_old_response_exploratory_p %>%
mutate(plot_stacked = map(Response, stacked_function))
pdf("02_Exploratoryanalysis_outputs/Gabazine_plots_stacked_p.pdf")
Gabazinedata_old_response_exploratory_p$plot_stacked
dev.off()plot_scatter1_function <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
plot_scatter1 <- ggpairs(data$data[[1]], aes(colour = Treatment, alpha = 0.3),
columns = c("Response_measurement", "Treatment", "fBehavioural_tank", "fDrug_test_no.", "ML_cm", "Test_time_day_hours", "fAcclimation_days", "fTank_intro_day", "fSystem", "Treatment_time_hours"),
upper = list(continuous = "smooth", combo = "box", discrete = "facetbar"),
lower = list(continuous = "cor", combo = "dot", discrete = "facetbar"),
diag = list(continuous = 'densityDiag', discrete = 'barDiag'),
title = response)
}
Gabazinedata_old_response_exploratory_p <- Gabazinedata_old_response_exploratory_p %>%
mutate(plot_scatter1 = map(Response, plot_scatter1_function))
pdf("02_Exploratoryanalysis_outputs/Gabazine_plot_scatter1_p.pdf", width = 32, height = 24)
Gabazinedata_old_response_exploratory_p$plot_scatter1
dev.off()
plot_scatter2_function <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
plot_scatter2 <- ggpairs(data$data[[1]], aes(colour = Treatment, alpha = 0.3),
columns = c("Response_measurement", "Treatment", "Behavioural_tank", "Drug_test_no.", "ML_cm", "Test_time_day_hours", "Acclimation_days", "Tank_intro_day", "System", "Treatment_time_hours"),
upper = list(continuous = "smooth", combo = "box", discrete = "facetbar"),
lower = list(continuous = "cor", combo = "dot", discrete = "facetbar"),
diag = list(continuous = 'densityDiag', discrete = 'barDiag'),
title = response)
}
Gabazinedata_old_response_exploratory_p <- Gabazinedata_old_response_exploratory_p %>%
mutate(plot_scatter2 = map(Response, plot_scatter2_function))
pdf("02_Exploratoryanalysis_outputs/Gabazine_plot_scatter2_p.pdf", width = 32, height = 24)
Gabazinedata_old_response_exploratory_p$plot_scatter2
dev.off()
save(Gabazinedata_old_response_exploratory_p, file = '02_Exploratoryanalysis_outputs/Gabazinedata_old_response_exploratory_p.RData')model_co2drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_size <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_method <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no., family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no., family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no., family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_housing <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_system <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column in the Gabazinedata_old_response_p dataframe
Gabazinedata_old_response_models_01_binompoisgaus_p <- Gabazinedata_old_response_p %>%
mutate(model_co2drug = map(Response, model_co2drug),
model_size = map(Response, model_size),
model_method = map(Response, model_method),
model_drug = map(Response, model_drug),
model_housing = map(Response, model_housing),
model_system = map(Response, model_system))
save(Gabazinedata_old_response_models_01_binompoisgaus_p, file = '03_Modelfit_outputs/Gabazinedata_old_response_models_01_binompoisgaus_p.RData')
#check correct family for each response var
Gabazinedata_old_response_models_01_binompoisgaus_p$model_co2drugCalculate loo values for each model and each response variable
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug[[1]])
}
loo_size_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_size <- loo(data$model_size[[1]])
}
loo_method_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_method <- loo(data$model_method[[1]])
}
loo_drug_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug[[1]])
}
loo_housing_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing[[1]])
}
loo_system_function <- function(response) {
data <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_system <- loo(data$model_system[[1]])
}
#save each loo as a separate column
Gabazinedata_old_response_models_loo_01_binompoisgaus_p <- Gabazinedata_old_response_models_01_binompoisgaus_p %>%
mutate(
loo_co2drug_binompoisgaus = map(Response, loo_co2drug_function),
loo_size_binompoisgaus = map(Response, loo_size_function),
loo_method_binompoisgaus = map(Response, loo_method_function),
loo_drug_binompoisgaus = map(Response, loo_drug_function),
loo_housing_binompoisgaus = map(Response, loo_housing_function),
loo_system_binompoisgaus = map(Response, loo_system_function))
save(Gabazinedata_old_response_models_loo_01_binompoisgaus_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_01_binompoisgaus_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Gabazinedata_old_response_models_loo_01_binompoisgaus_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_binompoisgaus[[1]], data$loo_size_binompoisgaus[[1]], data$loo_method_binompoisgaus[[1]], data$loo_drug_binompoisgaus[[1]], data$loo_housing_binompoisgaus[[1]], data$loo_system_binompoisgaus[[1]])
}
Gabazinedata_old_response_models_loo_01_binompoisgaus_p <- Gabazinedata_old_response_models_loo_01_binompoisgaus_p %>%
mutate(loo_comp_binompoisgaus = map(Response, loo_comp_function))
save(Gabazinedata_old_response_models_loo_01_binompoisgaus_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_01_binompoisgaus_p.RData')Check loo values trustworthy
#check pareto k OK for each and that loo outputs are reliable
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_co2drug_binompoisgaus
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_size_binompoisgaus
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_method_binompoisgaus
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_drug_binompoisgaus
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_housing_binompoisgaus
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_system_binompoisgausCompare loo values to choose variables in model for each reponse variable
Gabazinedata_old_response_models_loo_01_binompoisgaus_p$loo_comp_binompoisgausPick chosen model
#Create new column of chosen model for each response variable
#No._visits_A = housing, Time_A_s_no_0 = co2drug, Active_time_s = co2drug, Dist = co2drug, Vel = co2drug, Soft_touch_binomial = co2drug, Aggressive_touch_binomial = co2drug, Soft_touch_no._no_0 = housing, Aggressive_touch_no._no_0 = method, Latency_soft_touch_s = co2drug, Latency_aggressive_touch_s = co2drug
Gabazinedata_old_response_chosen_model_01_binompoisgaus_p <- Gabazinedata_old_response_models_loo_01_binompoisgaus_p %>%
pivot_longer(c(model_co2drug, model_method, model_housing), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Gabazinedata_old_response_chosen_model_01_binompoisgaus_p <- Gabazinedata_old_response_chosen_model_01_binompoisgaus_p %>%
slice(c(3, 4, 7, 10, 13, 16, 19, 24, 26, 28, 31)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Gabazinedata_old_response_chosen_model_01_binompoisgaus_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_01_binompoisgaus_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Gabazinedata_old_response_chosen_model_01_binompoisgaus <- Gabazinedata_old_response_chosen_model_01_binompoisgaus_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag <- Gabazinedata_old_response_chosen_model_01_binompoisgaus %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_01_binompoisgaus.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag <- Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_01_binompoisgaus.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
data1 <- Gabazinedata_old_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_01_binompoisgaus_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Gabazinedata_old_response_chosen_model_01_binompoisgaus_diag$Response, diag_resid_function)#Count data poisson distribution no good, instead do negbinomial
model_co2drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_size <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_method <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no., family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_housing <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_system <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column
Gabazinedata_old_response_models_03_NB_p <- Gabazinedata_old_response_p %>%
mutate(model_co2drug_NB = map(Response, model_co2drug),
model_size_NB = map(Response, model_size),
model_method_NB = map(Response, model_method),
model_drug_NB = map(Response, model_drug),
model_housing_NB = map(Response, model_housing),
model_system_NB = map(Response, model_system))
save(Gabazinedata_old_response_models_03_NB_p, file = '03_Modelfit_outputs/Gabazinedata_old_response_models_03_NB_p.RData')
#check correct family for each response var
Gabazinedata_old_response_models_03_NB_p$model_co2drug_NBCalculate loo values for each model and each response variable
Gabazinedata_old_response_models_03_NB_p <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response %in% c('No._visits_A', 'Soft_touch_no._no_0', 'Aggressive_touch_no._no_0'))
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_NB[[1]])
}
loo_size_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_NB[[1]])
}
loo_method_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_NB[[1]])
}
loo_drug_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_NB[[1]])
}
loo_housing_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_NB[[1]])
}
loo_system_function <- function(response) {
data <- Gabazinedata_old_response_models_03_NB_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_NB[[1]])
}
#save each loo as a separate column
Gabazinedata_old_response_models_loo_03_NB_p <- Gabazinedata_old_response_models_03_NB_p %>%
mutate(
loo_co2drug_NB = map(Response, loo_co2drug_function),
loo_size_NB = map(Response, loo_size_function),
loo_method_NB = map(Response, loo_method_function),
loo_drug_NB = map(Response, loo_drug_function),
loo_housing_NB = map(Response, loo_housing_function),
loo_system_NB = map(Response, loo_system_function))
save(Gabazinedata_old_response_models_loo_03_NB_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_03_NB_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Gabazinedata_old_response_models_loo_03_NB_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_NB[[1]], data$loo_size_NB[[1]], data$loo_method_NB[[1]], data$loo_drug_NB[[1]], data$loo_housing_NB[[1]], data$loo_system_NB[[1]])
}
Gabazinedata_old_response_models_loo_03_NB_p <- Gabazinedata_old_response_models_loo_03_NB_p %>%
mutate(loo_comp_NB = map(Response, loo_comp_function))
save(Gabazinedata_old_response_models_loo_03_NB_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_03_NB_p.RData')Check loo values trustworthy
#check pareto k OK for each and that loo outputs are reliable
Gabazinedata_old_response_models_loo_03_NB_p$loo_co2drug_NB
Gabazinedata_old_response_models_loo_03_NB_p$loo_size_NB
Gabazinedata_old_response_models_loo_03_NB_p$loo_method_NB
Gabazinedata_old_response_models_loo_03_NB_p$loo_drug_NB
Gabazinedata_old_response_models_loo_03_NB_p$loo_housing_NB
Gabazinedata_old_response_models_loo_03_NB_p$loo_system_NBCompare loo values to choose variables in model for each reponse variable
Gabazinedata_old_response_models_loo_03_NB_p$loo_comp_NBPick chosen model
#Create new column of chosen model for each response variable
#No._visits_A = size, Soft_touch_no._no_0 = method, Aggressive_touch_no._no_0 = co2drug
Gabazinedata_old_response_chosen_model_03_NB_p <- Gabazinedata_old_response_models_loo_03_NB_p %>%
pivot_longer(c(model_co2drug_NB, model_method_NB, model_size_NB), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Gabazinedata_old_response_chosen_model_03_NB_p <- Gabazinedata_old_response_chosen_model_03_NB_p %>%
slice(c(3, 5, 7)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Gabazinedata_old_response_chosen_model_03_NB_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_03_NB_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Gabazinedata_old_response_chosen_model_03_NB <- Gabazinedata_old_response_chosen_model_03_NB_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_03_NB %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_diag_03_NB <- Gabazinedata_old_response_chosen_model_03_NB %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Gabazinedata_old_response_chosen_model_diag_03_NB, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_diag_03_NB.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_03_NB.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_diag_03_NB$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_diag_03_NB %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_diag_03_NB <- Gabazinedata_old_response_chosen_model_diag_03_NB %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Gabazinedata_old_response_chosen_model_diag_03_NB, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_diag_03_NB.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_03_NB.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_diag_03_NB$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data1 <- Gabazinedata_old_response_chosen_model_03_NB %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_03_NB_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Gabazinedata_old_response_chosen_model_diag_03_NB$Response, diag_resid_function)model_co2drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_size <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_method <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_drug <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no., family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_housing <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_system <- function(response) {
data <- Gabazinedata_old_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Dist|^Active_time_s', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column in the Gabazinedata_old_response_p dataframe
Gabazinedata_old_response_models_04_gammalog_p <- Gabazinedata_old_response_p %>%
mutate(model_co2drug_gammalog = map(Response, model_co2drug),
model_size_gammalog = map(Response, model_size),
model_method_gammalog = map(Response, model_method),
model_drug_gammalog = map(Response, model_drug),
model_housing_gammalog = map(Response, model_housing),
model_system_gammalog = map(Response, model_system))
save(Gabazinedata_old_response_models_04_gammalog_p, file = '03_Modelfit_outputs/Gabazinedata_old_response_models_04_gammalog_p.RData')
#check correct family for each response var
Gabazinedata_old_response_models_04_gammalog_p$model_co2drug_gammalogCalculate loo values for each model and each response variable
Gabazinedata_old_response_models_04_gammalog_p <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response %in% c('Time_A_s_no_0', 'Active_time_s', 'Dist', 'Latency_soft_touch_s', 'Latency_aggressive_touch_s'))
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_gammalog[[1]])
}
loo_size_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_gammalog[[1]])
}
loo_method_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_gammalog[[1]])
}
loo_drug_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_gammalog[[1]])
}
loo_housing_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_gammalog[[1]])
}
loo_system_function <- function(response) {
data <- Gabazinedata_old_response_models_04_gammalog_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_gammalog[[1]])
}
#save each loo as a separate column
Gabazinedata_old_response_models_loo_04_gammalog_p <- Gabazinedata_old_response_models_04_gammalog_p %>%
mutate(
loo_co2drug_gammalog = map(Response, loo_co2drug_function),
loo_size_gammalog = map(Response, loo_size_function),
loo_method_gammalog = map(Response, loo_method_function),
loo_drug_gammalog = map(Response, loo_drug_function),
loo_housing_gammalog = map(Response, loo_housing_function),
loo_system_gammalog = map(Response, loo_system_function))
save(Gabazinedata_old_response_models_loo_04_gammalog_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_04_gammalog_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Gabazinedata_old_response_models_loo_04_gammalog_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_gammalog[[1]], data$loo_size_gammalog[[1]], data$loo_method_gammalog[[1]], data$loo_drug_gammalog[[1]], data$loo_housing_gammalog[[1]], data$loo_system_gammalog[[1]])
}
Gabazinedata_old_response_models_loo_04_gammalog_p <- Gabazinedata_old_response_models_loo_04_gammalog_p %>%
mutate(loo_comp_gammalog = map(Response, loo_comp_function))
save(Gabazinedata_old_response_models_loo_04_gammalog_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_models_loo_04_gammalog_p.RData')Check loo values trustworthy
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_co2drug_gammalog
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_size_gammalog
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_method_gammalog
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_drug_gammalog
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_housing_gammalog
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_system_gammalogCompare loo values to choose variables in model for each reponse variable
Gabazinedata_old_response_models_loo_04_gammalog_p$loo_comp_gammalogPick chosen model
#Create new column of chosen model for each response variable
#Time_A_s_no_0 = co2drug, Active_time_s = co2drug, Dist = co2drug, Latency_soft_touch_s = co2drug, Latency_aggressive_touch_s = method
Gabazinedata_old_response_chosen_model_04_gammalog_p <- Gabazinedata_old_response_models_loo_04_gammalog_p %>%
pivot_longer(c(model_co2drug_gammalog, model_method_gammalog), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Gabazinedata_old_response_chosen_model_04_gammalog_p <- Gabazinedata_old_response_chosen_model_04_gammalog_p %>%
slice(c(1, 3, 5, 7, 10)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Gabazinedata_old_response_chosen_model_04_gammalog_p, file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_04_gammalog_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Gabazinedata_old_response_chosen_model_04_gammalog <- Gabazinedata_old_response_chosen_model_04_gammalog_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_04_gammalog %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_diag_04_gammalog <- Gabazinedata_old_response_chosen_model_04_gammalog %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Gabazinedata_old_response_chosen_model_diag_04_gammalog, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_diag_04_gammalog.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_04_gammalog.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_diag_04_gammalog$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Gabazinedata_old_response_chosen_model_diag_04_gammalog %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Gabazinedata_old_response_chosen_model_diag_04_gammalog <- Gabazinedata_old_response_chosen_model_diag_04_gammalog %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Gabazinedata_old_response_chosen_model_diag_04_gammalog, file = '05_Modelvalidation_outputs/Gabazinedata_old_response_chosen_model_diag_04_gammalog.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_04_gammalog.pdf", width = 32, height = 24)
Gabazinedata_old_response_chosen_model_diag_04_gammalog$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data1 <- Gabazinedata_old_response_chosen_model_04_gammalog %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_04_gammalog_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Gabazinedata_old_response_chosen_model_diag_04_gammalog$Response, diag_resid_function)#1st round: binomial for binomial data, poisson for count data, gaussian for continuous data
#Keep soft_touch_binomial, aggressive_touch_binomial
#Keep continuous Vel, Time_A_s_no_0
load(file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_01_binompoisgaus_p.RData')
#2nd round: negative binomial for count data
#Keep No._visits_A, soft_touch_no._no_0, aggressive_touch_no._no_0
load(file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_03_NB_p.RData')
#3rd round: Gamma, log for continuous data gaussian, identity is no good for
#Keep Active_time_s, Dist, latency_soft_touch, latency_aggressive_touch
load(file = '04_Variableselection_outputs/Gabazinedata_old_response_chosen_model_04_gammalog_p.RData')
#Tidy data
Gabazinedata_old_response_chosen_model_01_binompoisgaus_p <- Gabazinedata_old_response_chosen_model_01_binompoisgaus_p %>%
filter(Response %in% c('Soft_touch_binomial', 'Aggressive_touch_binomial', 'Vel', 'Time_A_s_no_0'))
Gabazinedata_old_response_chosen_model_03_NB_p
Gabazinedata_old_response_chosen_model_04_gammalog_p <- Gabazinedata_old_response_chosen_model_04_gammalog_p %>%
filter(Response != 'Time_A_s_no_0')
Gabazinedata_old_final_models <- bind_rows(Gabazinedata_old_response_chosen_model_01_binompoisgaus_p, Gabazinedata_old_response_chosen_model_03_NB_p, Gabazinedata_old_response_chosen_model_04_gammalog_p)
save(Gabazinedata_old_final_models, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models.RData')Calculate effect sizes for each comparison of interest, including prob of an effect, prob > 20% effect and 80% and 95% HPDI intervals Use exp() on models with log link = fold change Use exp() on models with logit link = odds ratio Linear models manually calculate fold change
#Make a dataset of effect sizes for each contrast = effect_data for each response variable and save in new column
#Models with log link use exp() to calculate fold change, models with gaussian, identity manually calculate fold change (Time_A_s_no_0 and Vel)
effect_data_function <- function(response) {
data <- Gabazinedata_old_final_models %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
co2_effect_draws <- emmeans(data$model_chosen[[1]], ~ CO2|Drug, type='link') %>%
gather_emmeans_draws() %>%
spread(key = CO2, value = .value) %>%
mutate(fold_change = Elevated/Ambient) %>%
select(-Ambient, -Elevated)
co2_effect_emmeans <- co2_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Elevated - Ambient') %>%
filter(Drug=='Sham') %>%
droplevels()
if (co2_effect_emmeans$fold_change[[1]] < 1) {
co2_effect_prob <- co2_effect_draws %>%
filter(Drug=='Sham') %>%
group_by(fold_change, Drug) %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
co2_effect_prob <- co2_effect_draws %>%
filter(Drug=='Sham') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
co2_effect_data <- co2_effect_prob %>% full_join(co2_effect_emmeans)
drug_effect_draws <- emmeans(data$model_chosen[[1]], ~ Drug|CO2, type='link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(fold_change = Gabazine/Sham) %>%
select(-Gabazine, -Sham)
drug_effect_emmeans_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Gabazine - Sham') %>%
droplevels()
if (drug_effect_emmeans_ambient$fold_change[[1]] < 1) {
drug_effect_prob_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_ambient <- drug_effect_prob_ambient %>% full_join(drug_effect_emmeans_ambient)
drug_effect_emmeans_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Gabazine - Sham') %>%
droplevels()
if (drug_effect_emmeans_elevated$fold_change[[1]] < 1) {
drug_effect_prob_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_elevated <- drug_effect_prob_elevated %>% full_join(drug_effect_emmeans_elevated)
drug_effect_data <- drug_effect_data_ambient %>% full_join(drug_effect_data_elevated)
interaction_effect_draws <- emmeans(data$model_chosen[[1]], ~Drug|CO2, type = 'link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(Gab_effect_fold_change = Gabazine / Sham) %>%
select(-Gabazine, -Sham) %>%
spread(key=CO2, value=Gab_effect_fold_change) %>%
mutate(Gab_effect_across_co2_fold_change = Elevated / Ambient) %>%
select(-Ambient, -Elevated)
interaction_effect_emmeans <- interaction_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Interaction')
if (interaction_effect_emmeans$Gab_effect_across_co2_fold_change[[1]] < 1) {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(Gab_effect_across_co2_fold_change<1)/n())*100,
Perc20 = (sum(Gab_effect_across_co2_fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
} else {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(Gab_effect_across_co2_fold_change>1)/n())*100,
Perc20 = (sum(Gab_effect_across_co2_fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
}
interaction_effect_data = interaction_effect_prob %>% full_join(interaction_effect_emmeans)
effect_data <- full_join(co2_effect_data, drug_effect_data) %>%
full_join(., interaction_effect_data)
}
else {
co2_effect_draws <- emmeans(data$model_chosen[[1]], revpairwise ~ CO2|Drug, type='link')$contrast %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))
co2_effect_emmeans <- co2_effect_draws %>%
group_by(contrast, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
dplyr::filter(Drug=='Sham') %>%
droplevels()
if (co2_effect_emmeans$exp.value[[1]] < 1) {
co2_effect_prob <- co2_effect_draws %>%
dplyr::filter(Drug=='Sham') %>%
group_by(contrast, Drug) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
co2_effect_prob <- co2_effect_draws %>%
dplyr::filter(Drug=='Sham') %>%
group_by(contrast, Drug) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
co2_effect_data <- co2_effect_prob %>% full_join(co2_effect_emmeans)
drug_effect_draws <- emmeans(data$model_chosen[[1]], revpairwise ~ Drug|CO2, type='link')$contrast %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))
drug_effect_emmeans_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
droplevels()
if (drug_effect_emmeans_ambient$exp.value[[1]] < 1) {
drug_effect_prob_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_ambient <- drug_effect_prob_ambient %>% full_join(drug_effect_emmeans_ambient)
drug_effect_emmeans_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
droplevels()
if (drug_effect_emmeans_elevated$exp.value[[1]] < 1) {
drug_effect_prob_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_elevated <- drug_effect_prob_elevated %>% full_join(drug_effect_emmeans_elevated)
drug_effect_data <- drug_effect_data_ambient %>% full_join(drug_effect_data_elevated)
interaction_effect_draws <- emmeans(data$model_chosen[[1]], ~Drug|CO2, type = 'link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(Gab_effect = Gabazine - Sham) %>%
dplyr::select(-Gabazine, -Sham) %>%
spread(key=CO2, value=Gab_effect) %>%
mutate(Gab_effect_across_co2 = Elevated - Ambient,
exp.Gab_effect_across_co2 = exp(Gab_effect_across_co2)) %>%
select(-Ambient, -Elevated)
interaction_effect_emmeans <- interaction_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Interaction')
if (interaction_effect_emmeans$exp.Gab_effect_across_co2[[1]] < 1) {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(exp.Gab_effect_across_co2<1)/n())*100,
Perc20 = (sum(exp.Gab_effect_across_co2<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
} else {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(exp.Gab_effect_across_co2>1)/n())*100,
Perc20 = (sum(exp.Gab_effect_across_co2>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
}
interaction_effect_data = interaction_effect_prob %>% full_join(interaction_effect_emmeans)
effect_data <- full_join(co2_effect_data, drug_effect_data) %>%
full_join(., interaction_effect_data)
}
}
Gabazinedata_old_final_models_effect_data <- Gabazinedata_old_final_models %>%
mutate(effect_data = map(Response, effect_data_function))
save(Gabazinedata_old_final_models_effect_data, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data.RData') #Partial plots, saved in a new column in data frame
#X axis limit defined by maximum upper CI * 1.2 and x axis labels defined for each response variable
#If used gaussian, identity (Time_A_s_no_0 and Vel) do NOT exponentiate value (no backtransformation needed), if used gamma, log or NB, log use exp for backtransforming, if used binomial, logit use plogis for backtransformation
plot_partial_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = .value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, (max(newdata$.upper)*1.2)))
} else if (grepl('binomial', response_var)) {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(plogis.value=plogis(.value)) %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(plogis.value=plogis(.value))%>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = plogis.value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, 1))
}
else {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value)) %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))%>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = exp.value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, (max(newdata$exp.value.upper)*1.2)))
}
if (response == 'Time_A_s_no_0') {
plot_partial <- plot_partial + scale_x_continuous('Time in Zone A (s) (no 0)', expand = c(0,0))
}
else if (response == 'Active_time_s') {
plot_partial <- plot_partial + scale_x_continuous('Active time (s)', expand = c(0,0))
}
else if (response == 'Vel') {
plot_partial <- plot_partial + scale_x_continuous('Average speed (cm/s)', expand = c(0,0))
}
else if (response == 'Soft_touch_binomial') {
plot_partial <- plot_partial + scale_x_continuous('Proportion of squid that touched mirror softly', expand = c(0,0))
}
else if (response == 'Aggressive_touch_binomial') {
plot_partial <- plot_partial + scale_x_continuous('Proportion of squid that touched mirror aggressively', expand = c(0,0))
}
else if (response == 'No._visits_A') {
plot_partial <- plot_partial + scale_x_continuous('No. of visits to Zone A', expand = c(0,0))
}
else if (response == 'Soft_touch_no._no_0') {
plot_partial <- plot_partial + scale_x_continuous('No. of soft mirror touches (no 0)', expand = c(0,0))
}
else if (response == 'Aggressive_touch_no._no_0') {
plot_partial <- plot_partial + scale_x_continuous('No. of aggressive mirror touches (no 0)', expand = c(0,0))
}
else if (response == 'Dist') {
plot_partial <- plot_partial + scale_x_continuous('Total distance moved (cm)', expand = c(0,0))
}
else if (response == 'Latency_soft_touch_s') {
plot_partial <- plot_partial + scale_x_continuous('Latency to first soft mirror touch (s)', expand = c(0,0))
}
else if (response == 'Latency_aggressive_touch_s') {
plot_partial <- plot_partial + scale_x_continuous('Latency to first aggressive mirror touch (s)', expand = c(0,0))
}
}
Gabazinedata_old_final_models_effect_data_plots <- Gabazinedata_old_final_models_effect_data_plots %>%
mutate(plot_partial = map(Response, plot_partial_function))
save(Gabazinedata_old_final_models_effect_data_plots, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots.RData')#Show effect sizes as fold change (by using exp() values from models with log link and manually calculating fold change for linear models(Time_A_s_no_0 and Vel)) - for binomial models effect size is an odds ratio
plot_caterpillar_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0', response_var)) {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=fold_change, x=CO2,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=fold_change, x=Drug,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=Gab_effect_across_co2_fold_change, x = contrast,
ymin = .lower, ymax = .upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=Gab_effect_across_co2_fold_change, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(y=fold_change, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=fold_change, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=Gab_effect_across_co2_fold_change, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', trans=scales::log2_trans()) +
coord_flip() +
theme_classic()
}
else if (grepl('^Vel', response_var)) {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=fold_change, x=CO2,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=fold_change, x=Drug,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=Gab_effect_across_co2_fold_change, x = contrast,
ymin = .lower, ymax = .upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=Gab_effect_across_co2_fold_change, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(y=fold_change, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=fold_change, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=Gab_effect_across_co2_fold_change, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', lim = c(0.499,2.05), breaks = c(0, 0.5, 1, 2), trans=scales::log2_trans(), labels = scales::number_format(accuracy = 0.1)) +
coord_flip() +
theme_classic()
}
else {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=exp.value, x=CO2,
ymin=exp.value.lower, ymax=exp.value.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.value, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=exp.value, x=Drug,
ymin=exp.value.lower, ymax=exp.value.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.value, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=exp.Gab_effect_across_co2, x = contrast,
ymin = exp.Gab_effect_across_co2.lower, ymax = exp.Gab_effect_across_co2.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.Gab_effect_across_co2, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(exp.value, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=exp.value, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=exp.Gab_effect_across_co2, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Gab effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', trans=scales::log2_trans()) +
coord_flip() +
theme_classic()
}
}
Gabazinedata_old_final_models_effect_data_plots <- Gabazinedata_old_final_models_effect_data_plots %>%
mutate(plot_caterpillar = map(Response, plot_caterpillar_function))
save(Gabazinedata_old_final_models_effect_data_plots, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots.RData') Look at partial and caterpillar plots created alongside quick overview partial and caterpillar plots made automatically to check if calculations were done correctly
#Partial plot and caterpillar plot alongside overview plots to check look like expected
#overview cat plot is a good comparison, but wont match exactly as CIs are with quantiles, and plotted on the link scale (rather than HPDI intervals and on the response scale)
check_plots_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == response)
plot_partial <- data$plot_partial[[1]]
plot_caterpillar <- data$plot_caterpillar[[1]]
overview_partial <- ggemmeans(data$model_chosen[[1]], terms = ~CO2 * Drug) %>% plot
overview_caterpillar <- mcmc_intervals(data$model_chosen[[1]], pars = c("b_CO2Elevated", "b_DrugGabazine", "b_CO2Elevated:DrugGabazine"), prob = 0.8, prob_outer = 0.95, point_est = "median")
check_plots <- (overview_partial + overview_caterpillar) / (plot_partial + plot_caterpillar)
ggsave(check_plots, file = paste("06_Modelinvestigation_outputs/check_plots_", response, ".pdf", sep = ""), width = 40, height = 20, unit = "cm", device = cairo_pdf)
}
map(Gabazinedata_old_final_models_effect_data_plots$Response, check_plots_function)#Sample size each treatment group, for each response variable
table_n_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == response)
data1 <- data$data
data1 %>%
as.data.frame() %>%
group_by(CO2, Drug) %>%
summarise(N = n()) %>%
kable(caption = response)
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_n = map(Response, table_n_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')
#Estimate +- 95% CI for each treatment group, across all response variables in a new column (table version of partial plot)
#On response scale
table_estimates_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB', model_name)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = prob,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('gamma', model_name)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = response,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = emmean,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('binomial', response_var)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = response,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_estimates = map(Response, table_estimates_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')
#Effect sizes +- 95% CI for each contrast, across all response variables in a new column (table version of caterpillar plot)
#exp for all models to show fold-change/odds ratio. Fold changes calculated manually for linear models
table_contrasts_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
options(knitr.kable.NA = '-')
effect_data <- data$effect_data[[1]]
if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
effect_data %>%
as.data.frame() %>%
filter(.width == 0.95) %>%
select(contrast, Drug, CO2, fold_change, .lower, .upper, Gab_effect_across_co2_fold_change, Perc, Perc20) %>%
mutate(Treatment = coalesce(Drug, CO2),
Estimate = coalesce(fold_change, Gab_effect_across_co2_fold_change),
'Lower HPDI' = .lower,
'Upper HPDI' = .upper) %>%
select(contrast, Treatment, Estimate, 'Lower HPDI', 'Upper HPDI', Perc, Perc20) %>%
rename(Contrast = contrast,
Prob = Perc,
'Prob(20%)' = Perc20) %>%
kable(caption = caption, digits = 2)
}
else {
effect_data %>%
as.data.frame() %>%
filter(.width == 0.95) %>%
select(contrast, Drug, CO2, exp.value, exp.value.lower, exp.value.upper, exp.Gab_effect_across_co2, exp.Gab_effect_across_co2.lower, exp.Gab_effect_across_co2.upper, Perc, Perc20) %>%
mutate(Treatment = coalesce(Drug, CO2),
Estimate = coalesce(exp.value, exp.Gab_effect_across_co2),
'Lower HPDI' = coalesce(exp.value.lower, exp.Gab_effect_across_co2.lower),
'Upper HPDI' = coalesce(exp.value.upper, exp.Gab_effect_across_co2.upper)) %>%
select(contrast, Treatment, Estimate, 'Lower HPDI', 'Upper HPDI', Perc, Perc20) %>%
rename(Contrast = contrast,
Prob = Perc,
'Prob(20%)' = Perc20) %>%
kable(caption = caption, digits = 2)
}
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_contrasts = map(Response, table_contrasts_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')
#Output for each table on link scale, including other explanatory variables included
table_summary_link_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
caption <- paste(response, " ~ CO2", " * Drug", caption1, caption2, sep = "")
model <- data$model_chosen[[1]]
table_summary_link <- model$fit %>%
tidyMCMC(estimate.method = 'median',
conf.int = TRUE, conf.method = 'HPDinterval') %>%
rename(Term = term,
Estimate = estimate,
SE = std.error,
'Lower HPDI' = conf.low,
'Upper HPDI' = conf.high) %>%
kable(caption = caption, digits = 2)
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_summary_link = map(Response, table_summary_link_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')
#Priors
table_priors_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB|gamma', model_name)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope', shape)
} else if (grepl('binomial', response_var)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope')
} else {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope', sigma)
}
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_priors = map(Response, table_priors_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')
#Model info (response var, explanatory vars, family and link) as well as priors used
table_modelinfo_priors_function <- function(response) {
data <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption_explanatory <- "CO2 * Drug"
} else if (grepl('size', model_name)) {
caption_explanatory <- "CO2 * Drug + Mantle length"
} else if (grepl('method', model_name)) {
caption_explanatory <- "CO2 * Drug + Behavioural tank + Time of test"
} else if (grepl('_drug', model_name)) {
caption_explanatory <- "CO2 * Drug + Drug test number"
} else if (grepl('housing', model_name)) {
caption_explanatory <- "CO2 * Drug + Number of acclimation days + Date introduced to treatment tank"
} else if (grepl('system', model_name)) {
caption_explanatory <- "CO2 * Drug + System"
}
if (grepl('NB', model_name)) {
caption_family_link <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption_family_link <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Vel', response_var)) {
caption_family_link <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption_family_link <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB|gamma', model_name)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept,
'shape prior' = shape) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior', 'shape prior')
} else if (grepl('binomial', response_var)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior')
} else {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept,
'sigma prior' = sigma) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior', 'sigma prior')
}
}
Gabazinedata_old_final_models_effect_data_plots_tables <- Gabazinedata_old_final_models_effect_data_plots_tables %>%
mutate(table_modelinfo_priors = map(Response, table_modelinfo_priors_function))
save(Gabazinedata_old_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Gabazinedata_old_final_models_effect_data_plots_tables.RData')boxplots_function <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
ggplot(data = data$data[[1]], aes(x = Treatment, y = Response_measurement)) +
geom_boxplot() +
scale_y_continuous(response)
}
Picrotoxindata_response_exploratory_p <- Picrotoxindata_response_p %>%
mutate(boxplots = map(Response, boxplots_function))
pdf("02_Exploratoryanalysis_outputs/Picrotoxin_boxplots_p.pdf")
Picrotoxindata_response_exploratory_p$boxplots
dev.off()stacked_function <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
data_plot <- data$data[[1]] %>%
mutate(Response_measurement = as.logical(Response_measurement))
plot_stacked <- ggplot(data = data_plot, aes(x = Treatment, fill = Response_measurement)) +
geom_bar(position = 'fill') +
scale_y_continuous(response_var)
} else {
}
}
Picrotoxindata_response_exploratory_p <- Picrotoxindata_response_exploratory_p %>%
mutate(plot_stacked = map(Response, stacked_function))
pdf("02_Exploratoryanalysis_outputs/Picrotoxin_plots_stacked_p.pdf")
Picrotoxindata_response_exploratory_p$plot_stacked
dev.off()
save(Picrotoxindata_response_exploratory_p, file = '02_Exploratoryanalysis_outputs/Picrotoxindata_response_exploratory_p.RData')plot_scatter1_function <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
plot_scatter1 <- ggpairs(data$data[[1]], aes(colour = Treatment, alpha = 0.3),
columns = c("Response_measurement", "Treatment", "fBehavioural_tank", "fDrug_test_no._dif_drugs", "ML_cm", "Test_time_day_hours", "fAcclimation_days", "fSystem", "Treatment_time_hours"),
upper = list(continuous = "smooth", combo = "box", discrete = "facetbar"),
lower = list(continuous = "cor", combo = "dot", discrete = "facetbar"),
diag = list(continuous = 'densityDiag', discrete = 'barDiag'),
title = response)
}
Picrotoxindata_response_exploratory_p <- Picrotoxindata_response_exploratory_p %>%
mutate(plot_scatter1 = map(Response, plot_scatter1_function))
pdf("02_Exploratoryanalysis_outputs/Picrotoxin_plot_scatter1_p.pdf", width = 32, height = 24)
Picrotoxindata_response_exploratory_p$plot_scatter1
dev.off()
plot_scatter2_function <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
plot_scatter2 <- ggpairs(data$data[[1]], aes(colour = Treatment, alpha = 0.3),
columns = c("Response_measurement", "Treatment", "Behavioural_tank", "Drug_test_no._dif_drugs", "ML_cm", "Test_time_day_hours", "Acclimation_days", "Tank_intro_day", "System", "Treatment_time_hours"),
upper = list(continuous = "smooth", combo = "box", discrete = "facetbar"),
lower = list(continuous = "cor", combo = "dot", discrete = "facetbar"),
diag = list(continuous = 'densityDiag', discrete = 'barDiag'),
title = response)
}
Picrotoxindata_response_exploratory_p <- Picrotoxindata_response_exploratory_p %>%
mutate(plot_scatter2 = map(Response, plot_scatter2_function))
pdf("02_Exploratoryanalysis_outputs/Picrotoxin_plot_scatter2_p.pdf", width = 32, height = 24)
Picrotoxindata_response_exploratory_p$plot_scatter2
dev.off()
save(Picrotoxindata_response_exploratory_p, file = '02_Exploratoryanalysis_outputs/Picrotoxindata_response_exploratory_p.RData')model_co2drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_size <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_method <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_housing <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
model_system <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('binomial', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = bernoulli('logit'))
} else if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = poisson(link = 'log'))
} else {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = gaussian('identity'))
}
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column in the Picrotoxindata_response_p dataframe
Picrotoxindata_response_models_01_binompoisgaus_p <- Picrotoxindata_response_p %>%
mutate(model_co2drug_binompoisgaus = map(Response, model_co2drug),
model_size_binompoisgaus = map(Response, model_size),
model_method_binompoisgaus = map(Response, model_method),
model_drug_binompoisgaus = map(Response, model_drug),
model_housing_binompoisgaus = map(Response, model_housing),
model_system_binompoisgaus = map(Response, model_system))
save(Picrotoxindata_response_models_01_binompoisgaus_p, file = '03_Modelfit_outputs/Picrotoxindata_response_models_01_binompoisgaus_p.RData')
#check correct family for each response var
Picrotoxindata_response_models_01_binompoisgaus_p$model_co2drug_binompoisgausCalculate loo values for each model and each response variable
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_binompoisgaus[[1]])
}
loo_size_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_binompoisgaus[[1]])
}
loo_method_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_binompoisgaus[[1]])
}
loo_drug_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_binompoisgaus[[1]])
}
loo_housing_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_binompoisgaus[[1]])
}
loo_system_function <- function(response) {
data <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_binompoisgaus[[1]])
}
#save each loo as a separate column
Picrotoxindata_response_models_loo_01_binompoisgaus_p <- Picrotoxindata_response_models_01_binompoisgaus_p %>%
mutate(
loo_co2drug_binompoisgaus = map(Response, loo_co2drug_function),
loo_size_binompoisgaus = map(Response, loo_size_function),
loo_method_binompoisgaus = map(Response, loo_method_function),
loo_drug_binompoisgaus = map(Response, loo_drug_function),
loo_housing_binompoisgaus = map(Response, loo_housing_function),
loo_system_binompoisgaus = map(Response, loo_system_function))
save(Picrotoxindata_response_models_loo_01_binompoisgaus_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_01_binompoisgaus_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Picrotoxindata_response_models_loo_01_binompoisgaus_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_binompoisgaus[[1]], data$loo_size_binompoisgaus[[1]], data$loo_method_binompoisgaus[[1]], data$loo_drug_binompoisgaus[[1]], data$loo_housing_binompoisgaus[[1]], data$loo_system_binompoisgaus[[1]])
}
Picrotoxindata_response_models_loo_01_binompoisgaus_p <- Picrotoxindata_response_models_loo_01_binompoisgaus_p %>%
mutate(loo_comp_binompoisgaus = map(Response, loo_comp_function))
save(Picrotoxindata_response_models_loo_01_binompoisgaus_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_01_binompoisgaus_p.RData')Check loo values trustworthy
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_co2drug_binompoisgaus
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_size_binompoisgaus
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_method_binompoisgaus
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_drug_binompoisgaus
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_housing_binompoisgaus
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_system_binompoisgausCompare loo values to choose variables in model for each reponse variable
Picrotoxindata_response_models_loo_01_binompoisgaus_p$loo_comp_binompoisgausPick chosen model
#Create new column of chosen model for each response variable
#No._visits_A = system, Soft_touch_no._no_0 = size, Time_A_s_no_0 = method, Latency_soft_touch_s = co2drug, Active_time_s = co2drug, Dist = co2drug, Vel = co2drug, Soft_touch_binomial = co2drug, Aggressive_touch_binomial = co2drug, Aggressive_touch_no._no_0 = housing, Latency_aggressive_touch_s = co2drug
Picrotoxindata_response_chosen_model_01_binompoisgaus_p <- Picrotoxindata_response_models_loo_01_binompoisgaus_p %>%
pivot_longer(c(model_co2drug_binompoisgaus, model_size_binompoisgaus, model_method_binompoisgaus, model_housing_binompoisgaus, model_system_binompoisgaus), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Picrotoxindata_response_chosen_model_01_binompoisgaus_p <- Picrotoxindata_response_chosen_model_01_binompoisgaus_p %>%
slice(c(5, 7, 13, 16, 21, 26, 31, 36, 41, 49, 51)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Picrotoxindata_response_chosen_model_01_binompoisgaus_p, file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_01_binompoisgaus_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Picrotoxindata_response_chosen_model_01_binompoisgaus < Picrotoxindata_response_chosen_model_01_binompoisgaus_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_01_binompoisgaus_diag <- Picrotoxindata_response_chosen_model_01_binompoisgaus %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Picrotoxindata_response_chosen_model_01_binompoisgaus_diag, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_01_binompoisgaus_diag.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_01_binompoisgaus.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_01_binompoisgaus_diag$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_01_binompoisgaus_diag %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_01_binompoisgaus_diag <- Picrotoxindata_response_chosen_model_01_binompoisgaus_diag %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Picrotoxindata_response_chosen_model_01_binompoisgaus_diag, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_01_binompoisgaus_diag.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_01_binompoisgaus.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_01_binompoisgaus_diag$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
data1 <- Picrotoxindata_response_chosen_model_01_binompoisgaus %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_01_binompoisgaus_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Picrotoxindata_response_chosen_model_01_binompoisgaus_diag$Response, diag_resid_function)model_co2drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_size <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_method <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_housing <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_system <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('No\\.', response_var, ignore.case = TRUE)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column
Picrotoxindata_response_models_03_NB_p <- Picrotoxindata_response_p %>%
mutate(model_co2drug_NB = map(Response, model_co2drug),
model_size_NB = map(Response, model_size),
model_method_NB = map(Response, model_method),
model_drug_NB = map(Response, model_drug),
model_housing_NB = map(Response, model_housing),
model_system_NB = map(Response, model_system))
save(Picrotoxindata_response_models_03_NB_p, file = '03_Modelfit_outputs/Picrotoxindata_response_models_03_NB_p.RData')
#check correct family for each response var
Picrotoxindata_response_models_03_NB_p$model_co2drug_NBCalculate loo values for each model and each response variable
Picrotoxindata_response_models_03_NB_p <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response %in% c('No._visits_A', 'Soft_touch_no._no_0', 'Aggressive_touch_no._no_0'))
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_NB[[1]])
}
loo_size_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_NB[[1]])
}
loo_method_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_NB[[1]])
}
loo_drug_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_NB[[1]])
}
loo_housing_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_NB[[1]])
}
loo_system_function <- function(response) {
data <- Picrotoxindata_response_models_03_NB_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_NB[[1]])
}
#save each loo as a separate column
Picrotoxindata_response_models_loo_03_NB_p <- Picrotoxindata_response_models_03_NB_p %>%
mutate(
loo_co2drug_NB = map(Response, loo_co2drug_function),
loo_size_NB = map(Response, loo_size_function),
loo_method_NB = map(Response, loo_method_function),
loo_drug_NB = map(Response, loo_drug_function),
loo_housing_NB = map(Response, loo_housing_function),
loo_system_NB = map(Response, loo_system_function))
save(Picrotoxindata_response_models_loo_03_NB_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_03_NB_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Picrotoxindata_response_models_loo_03_NB_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_NB[[1]], data$loo_size_NB[[1]], data$loo_method_NB[[1]], data$loo_drug_NB[[1]], data$loo_housing_NB[[1]], data$loo_system_NB[[1]])
}
Picrotoxindata_response_models_loo_03_NB_p <- Picrotoxindata_response_models_loo_03_NB_p %>%
mutate(loo_comp_NB = map(Response, loo_comp_function))
save(Picrotoxindata_response_models_loo_03_NB_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_03_NB_p.RData')Check loo values trustworthy
Picrotoxindata_response_models_loo_03_NB_p$loo_co2drug_NB
Picrotoxindata_response_models_loo_03_NB_p$loo_size_NB
Picrotoxindata_response_models_loo_03_NB_p$loo_method_NB
Picrotoxindata_response_models_loo_03_NB_p$loo_drug_NB
Picrotoxindata_response_models_loo_03_NB_p$loo_housing_NB
Picrotoxindata_response_models_loo_03_NB_p$loo_system_NBCompare loo values to choose variables in model for each reponse variable
Picrotoxindata_response_models_loo_03_NB_p$loo_comp_NBPick chosen model
#Create new column of chosen model for each response variable
#No._visist_A = system, Soft_touch_no._no_0 = size, Aggressive_touch_no._no_0 = housing
Picrotoxindata_response_chosen_model_03_NB_p <- Picrotoxindata_response_models_loo_03_NB_p %>%
pivot_longer(c(model_housing_NB, model_size_NB, model_system_NB), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Picrotoxindata_response_chosen_model_03_NB_p <- Picrotoxindata_response_chosen_model_03_NB_p %>%
slice(c(3, 5, 7)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Picrotoxindata_response_chosen_model_03_NB_p, file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_03_NB_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Picrotoxindata_response_chosen_model_03_NB <- Picrotoxindata_response_chosen_model_03_NB_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_03_NB %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_03_NB <- Picrotoxindata_response_chosen_model_03_NB %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Picrotoxindata_response_chosen_model_diag_03_NB, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_03_NB.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_03_NB.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_03_NB$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_diag_03_NB %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_03_NB <- Picrotoxindata_response_chosen_model_diag_03_NB %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Picrotoxindata_response_chosen_model_diag_03_NB, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_03_NB.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_03_NB.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_03_NB$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data1 <- Picrotoxindata_response_chosen_model_03_NB %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_03_NB_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Picrotoxindata_response_chosen_model_diag_03_NB$Response, diag_resid_function)model_co2drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_size <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_method <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_housing <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
model_system <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Time_A_s_no_0|^Latency|^Active_time_s|^Dist|^Vel', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 10000, warmup = 3000, thin = 5,
refresh = 0, seed = 814883)
} else {
}
}
#Fit each of the 6 models for each of the count response variables and output each model as a new column in the Picrotoxindata_response_p dataframe
Picrotoxindata_response_models_04_gammalog_p <- Picrotoxindata_response_p %>%
mutate(model_co2drug_gammalog = map(Response, model_co2drug),
model_size_gammalog = map(Response, model_size),
model_method_gammalog = map(Response, model_method),
model_drug_gammalog = map(Response, model_drug),
model_housing_gammalog = map(Response, model_housing),
model_system_gammalog = map(Response, model_system))
save(Picrotoxindata_response_models_04_gammalog_p, file = '03_Modelfit_outputs/Picrotoxindata_response_models_04_gammalog_p.RData')
#check correct family for each response var
Picrotoxindata_response_models_04_gammalog_p$model_co2drug_gammalogCalculate loo values for each model and each response variable
Picrotoxindata_response_models_04_gammalog_p <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response %in% c('Time_A_s_no_0', 'Latency_soft_touch_s', 'Latency_aggressive_touch_s'))
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_gammalog[[1]])
}
loo_size_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_gammalog[[1]])
}
loo_method_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_gammalog[[1]])
}
loo_drug_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_gammalog[[1]])
}
loo_housing_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_gammalog[[1]])
}
loo_system_function <- function(response) {
data <- Picrotoxindata_response_models_04_gammalog_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_gammalog[[1]])
}
#save each loo as a separate column
Picrotoxindata_response_models_loo_04_gammalog_p <- Picrotoxindata_response_models_04_gammalog_p %>%
mutate(
loo_co2drug_gammalog = map(Response, loo_co2drug_function),
loo_size_gammalog = map(Response, loo_size_function),
loo_method_gammalog = map(Response, loo_method_function),
loo_drug_gammalog = map(Response, loo_drug_function),
loo_housing_gammalog = map(Response, loo_housing_function),
loo_system_gammalog = map(Response, loo_system_function))
save(Picrotoxindata_response_models_loo_04_gammalog_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_04_gammalog_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Picrotoxindata_response_models_loo_04_gammalog_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_gammalog[[1]], data$loo_size_gammalog[[1]], data$loo_method_gammalog[[1]], data$loo_drug_gammalog[[1]], data$loo_housing_gammalog[[1]], data$loo_system_gammalog[[1]])
}
Picrotoxindata_response_models_loo_04_gammalog_p <- Picrotoxindata_response_models_loo_04_gammalog_p %>%
mutate(loo_comp_gammalog = map(Response, loo_comp_function))
save(Picrotoxindata_response_models_loo_04_gammalog_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_04_gammalog_p.RData')Check loo values trustworthy
Picrotoxindata_response_models_loo_04_gammalog_p$loo_co2drug_gammalog
Picrotoxindata_response_models_loo_04_gammalog_p$loo_size_gammalog
Picrotoxindata_response_models_loo_04_gammalog_p$loo_method_gammalog
Picrotoxindata_response_models_loo_04_gammalog_p$loo_drug_gammalog
Picrotoxindata_response_models_loo_04_gammalog_p$loo_housing_gammalog
Picrotoxindata_response_models_loo_04_gammalog_p$loo_system_gammalogCompare loo values to choose variables in model for each reponse variable
Picrotoxindata_response_models_loo_04_gammalog_p$loo_comp_gammalogPick chosen model
#Create new column of chosen model for each response variable
#Time_A_s_no_0 = co2drug, Latency_soft_touch_s = co2drug, Latency_aggresive_touch_s = system
Picrotoxindata_response_chosen_model_04_gammalog_p <- Picrotoxindata_response_models_loo_04_gammalog_p %>%
pivot_longer(c(model_co2drug_gammalog, model_system_gammalog), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Picrotoxindata_response_chosen_model_04_gammalog_p <- Picrotoxindata_response_chosen_model_04_gammalog_p %>%
slice(c(1, 3, 6)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Picrotoxindata_response_chosen_model_04_gammalog_p, file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_04_gammalog_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Picrotoxindata_response_chosen_model_04_gammalog <- Picrotoxindata_response_chosen_model_04_gammalog_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_04_gammalog %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_04_gammalog <- Picrotoxindata_response_chosen_model_04_gammalog %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Picrotoxindata_response_chosen_model_diag_04_gammalog, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_04_gammalog.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_04_gammalog.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_04_gammalog$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_diag_04_gammalog %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_04_gammalog <- Picrotoxindata_response_chosen_model_diag_04_gammalog %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Picrotoxindata_response_chosen_model_diag_04_gammalog, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_04_gammalog.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_04_gammalog.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_04_gammalog$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data1 <- Picrotoxindata_response_chosen_model_04_gammalog %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_04_gammalog_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Picrotoxindata_response_chosen_model_diag_04_gammalog$Response, diag_resid_function)data <- Picrotoxindata_response_p %>%
filter(Response == 'No._visits_A')
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = negbinomial('log'))
model_priors = c(
prior(normal(0, 8), class='Intercept'),
prior(normal(0,2.5), class='b'))
model_No._visits_A_system_priorsonly <- brm(formula, data = data$data[[1]],
prior = model_priors,
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99),
sample_prior = 'only')
prior_summary(model_No._visits_A_system_priorsonly)
save(model_No._visits_A_system_priorsonly, file = '03_Modelfit_outputs/model_No._visits_A_system_priorsonly.RData')
conditional_effects(model_No._visits_A_system_priorsonly) %>%
plot(points=TRUE)Plot shows how wide the priors allow the model to go (don’t want priors to influence data but don’t want to be too wide and sampling gets stuck out somewhere) Here, the priors specified allows No._visits_A to go from 0 to at least 5e+06 for each treatment group, including systems, so priors definitely not influencing the data
model_co2drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_size <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_method <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_housing <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = negbinomial('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_system <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('visits', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = negbinomial('log'))
model_priors = c(
prior(normal(0, 8), class='Intercept'),
prior(normal(0,2.5), class='b'))
brm(formula, data = data$data[[1]],
prior = model_priors,
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
#Fit each of the 6 models and save each model as a new column
Picrotoxindata_response_models_05_NBadjust_p <- Picrotoxindata_response_p %>%
mutate(model_co2drug_NBadjust = map(Response, model_co2drug),
model_size_NBadjust = map(Response, model_size),
model_method_NBadjust = map(Response, model_method),
model_drug_NBadjust = map(Response, model_drug),
model_housing_NBadjust = map(Response, model_housing),
model_system_NBadjust = map(Response, model_system))
save(Picrotoxindata_response_models_05_NBadjust_p, file = '03_Modelfit_outputs/Picrotoxindata_response_models_05_NBadjust_p.RData')
#check correct family for each response var
Picrotoxindata_response_models_05_NBadjust_p$model_co2drug_NBadjustCalculate loo values for each model and each response variable
Picrotoxindata_response_models_05_NBadjust_p <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response %in% c('No._visits_A'))
#Make a function to calculate loo for each model
loo_co2drug_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_co2drug <- loo(data$model_co2drug_NBadjust[[1]])
}
loo_size_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_size <- loo(data$model_size_NBadjust[[1]])
}
loo_method_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_method <- loo(data$model_method_NBadjust[[1]])
}
loo_drug_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_drug <- loo(data$model_drug_NBadjust[[1]])
}
loo_housing_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_housing <- loo(data$model_housing_NBadjust[[1]])
}
loo_system_function <- function(response) {
data <- Picrotoxindata_response_models_05_NBadjust_p %>%
filter(Response == response)
loo_system <- loo(data$model_system_NBadjust[[1]])
}
#save each loo as a separate column
Picrotoxindata_response_models_loo_05_NBadjust_p <- Picrotoxindata_response_models_05_NBadjust_p %>%
mutate(
loo_co2drug_NBadjust = map(Response, loo_co2drug_function),
loo_size_NBadjust = map(Response, loo_size_function),
loo_method_NBadjust = map(Response, loo_method_function),
loo_drug_NBadjust = map(Response, loo_drug_function),
loo_housing_NBadjust = map(Response, loo_housing_function),
loo_system_NBadjust = map(Response, loo_system_function))
save(Picrotoxindata_response_models_loo_05_NBadjust_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_05_NBadjust_p.RData')
#Now loop through and compare each loo with loo_compare saved as a new column
#(if just did one function to calculate loo each model would be called the same thing which isn't useful for loo_compare)
loo_comp_function <- function(response) {
data <- Picrotoxindata_response_models_loo_05_NBadjust_p %>%
filter(Response == response)
loo_compare(data$loo_co2drug_NBadjust[[1]], data$loo_size_NBadjust[[1]], data$loo_method_NBadjust[[1]], data$loo_drug_NBadjust[[1]], data$loo_housing_NBadjust[[1]], data$loo_system_NBadjust[[1]])
}
Picrotoxindata_response_models_loo_05_NBadjust_p <- Picrotoxindata_response_models_loo_05_NBadjust_p %>%
mutate(loo_comp_NBadjust = map(Response, loo_comp_function))
save(Picrotoxindata_response_models_loo_05_NBadjust_p, file = '04_Variableselection_outputs/Picrotoxindata_response_models_loo_05_NBadjust_p.RData')Check loo values trustworthy
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_co2drug_NBadjust
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_size_NBadjust
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_method_NBadjust
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_drug_NBadjust
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_housing_NBadjust
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_system_NBadjustCompare loo values to choose variables in model for each reponse variable
Picrotoxindata_response_models_loo_05_NBadjust_p$loo_comp_NBadjustPick chosen model
#Create new column of chosen model for each response variable
#No._visit_A = system
Picrotoxindata_response_chosen_model_05_NBadjust_p <- Picrotoxindata_response_models_loo_05_NBadjust_p %>%
pivot_longer(model_system_NBadjust, names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup()
Picrotoxindata_response_chosen_model_05_NBadjust_p <- Picrotoxindata_response_chosen_model_05_NBadjust_p %>%
slice(c(1)) %>%
select(Response, data, model_chosen_name, model_chosen)
save(Picrotoxindata_response_chosen_model_05_NBadjust_p, file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_05_NBadjust_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Picrotoxindata_response_chosen_model_05_NBadjust <- Picrotoxindata_response_chosen_model_05_NBadjust_p
#Chain Diagnostics
diag_chain_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_05_NBadjust %>%
filter(Response == response)
trace <- mcmc_plot(data$model_chosen[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen[[1]], type='neff_hist')
plot <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_05_NBadjust <- Picrotoxindata_response_chosen_model_05_NBadjust %>%
mutate(diag_chain = map(Response, diag_chain_function))
save(Picrotoxindata_response_chosen_model_diag_05_NBadjust, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_05_NBadjust.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_05_NBadjust.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_05_NBadjust$diag_chain
dev.off()
#Posterior probability checks
diag_pp_function <- function(response) {
data <- Picrotoxindata_response_chosen_model_diag_05_NBadjust %>%
filter(Response == response)
whole <- pp_check(data$model_chosen[[1]])
co2 <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = paste(response))
}
Picrotoxindata_response_chosen_model_diag_05_NBadjust <- Picrotoxindata_response_chosen_model_diag_05_NBadjust %>%
mutate(diag_pp = map(Response, diag_pp_function))
save(Picrotoxindata_response_chosen_model_diag_05_NBadjust, file = '05_Modelvalidation_outputs/Picrotoxindata_response_chosen_model_diag_05_NBadjust.RData')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_05_NBadjust.pdf", width = 32, height = 24)
Picrotoxindata_response_chosen_model_diag_05_NBadjust$diag_pp
dev.off()
#DHARMa residuals
diag_resid_function <- function(response) {
data1 <- Picrotoxindata_response_chosen_model_05_NBadjust %>%
filter(Response == response)
preds <- posterior_predict(data1$model_chosen[[1]], nsamples=250, summary=FALSE)
data2 <- data1$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf(paste("05_Modelvalidation_outputs/Gabazine_diag_resid_05_NBadjust_", response, ".pdf", sep = ""))
plot(resids)
dev.off()
}
map(Picrotoxindata_response_chosen_model_diag_05_NBadjust$Response, diag_resid_function)data <- Picrotoxindata_response_p %>%
filter(Response == 'Latency_aggressive_touch_s')
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = Gamma('log'))
model_priors = c(
prior(normal(0, 10), class='Intercept'),
prior(normal(0,2.5), class='b'))
model_latency_aggressive_priorsonly <- brm(formula, data = data$data[[1]],
prior = model_priors,
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99),
sample_prior = 'only')
prior_summary(model_latency_aggressive_priorsonly)
save(model_latency_aggressive_priorsonly, file = '03_Modelfit_outputs/model_latency_aggressive_priorsonly.RData')
conditional_effects(model_latency_aggressive_priorsonly) %>%
plot(points=TRUE)Plot shows how wide the priors allow the model to go (don’t want priors to influence data but don’t want to be too wide and sampling gets stuck out somewhere) Here, the priors specified allows latency_aggressive_touch_s to go from 0 to at least 6e+08 for each treatment group, including systems, so priors definitely not influencing the data
model_co2drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_size <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + ML_cm, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_method <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fBehavioural_tank + Test_time_day_hours, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_drug <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fDrug_test_no._dif_drugs, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_housing <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fAcclimation_days + Tank_intro_day, family = Gamma('log'))
brm(formula, data = data$data[[1]],
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
model_system <- function(response) {
data <- Picrotoxindata_response_p %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('Latency_aggressive', response_var)) {
formula <- bf(Response_measurement ~ CO2 * Drug + fSystem, family = Gamma('log'))
model_priors = c(
prior(normal(0, 10), class='Intercept'),
prior(normal(0,2.5), class='b'))
brm(formula, data = data$data[[1]],
prior = model_priors,
chains = 4, iter = 120000, warmup = 40000, thin = 80,
refresh = 0, seed = 814883,
control=list(adapt_delta=0.99))
} else {
}
}
#Fit each of the 6 models for Latency_aggressive_s and save each model as a new column
Picrotoxindata_response_models_06_gammalogadjust_p <- Picrotoxindata_response_p %>%
mutate(model_co2drug_gammalogadjust = map(Response, model_co2drug),
model_size_gammalogadjust = map(Response, model_size),
model_method_gammalogadjust = map(Response, model_method),
model_drug_gammalogadjust = map(Response, model_drug),
model_housing_gammalogadjust = map(Response, model_housing),
model_system_gammalogadjust = map(Response, model_system))
save(Picrotoxindata_response_models_06_gammalogadjust_p, file = '03_Modelfit_outputs/Picrotoxindata_response_models_06_gammalogadjust_p.RData')
#check correct family for each response var
Picrotoxindata_response_models_06_gammalogadjust_p$model_co2drug_gammalogadjustCalculate loo values for each model and each response variable, Check loo values trustworthy
Picrotoxindata_response_models_06_gammalogadjust_p <- Picrotoxindata_response_models_06_gammalogadjust_p %>%
filter(Response %in% c('Latency_aggressive_touch_s'))
(loo_co2drug <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_co2drug_gammalogadjust[[1]]))
(loo_size <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_size_gammalogadjust[[1]]))
(loo_method <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_method_gammalogadjust[[1]]))
(loo_drug <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_drug_gammalogadjust[[1]]))
(loo_housing <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_housing_gammalogadjust[[1]]))
(loo_system <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_system_gammalogadjust[[1]]))Compare loo values to choose variables in model for each reponse variable
loo_co2drug <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_co2drug_gammalogadjust[[1]])
loo_size <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_size_gammalogadjust[[1]])
loo_method <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_method_gammalogadjust[[1]])
loo_drug <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_drug_gammalogadjust[[1]])
loo_housing <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_housing_gammalogadjust[[1]])
loo_system <- loo(Picrotoxindata_response_models_06_gammalogadjust_p$model_system_gammalogadjust[[1]])
(loo_compare(loo_co2drug, loo_size, loo_method, loo_drug, loo_housing, loo_system))Pick chosen model
#Latency_aggressive_touch_s = system
Picrotoxindata_response_chosen_model_06_gammalogadjust_p <- Picrotoxindata_response_models_06_gammalogadjust_p %>%
pivot_longer(c(model_system_gammalogadjust), names_to = 'model_chosen_name', values_to = 'model_chosen') %>%
ungroup() %>%
select(Response, data, model_chosen_name, model_chosen)
save(Picrotoxindata_response_chosen_model_06_gammalogadjust_p, file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_06_gammalogadjust_p.RData')Chain diagnostics, posterior probability checks and DHARMa residuals
Picrotoxindata_response_chosen_model_06_gammalogadjust <- Picrotoxindata_response_chosen_model_06_gammalogadjust_p
#Chain Diagnostics
data <- Picrotoxindata_response_chosen_model_06_gammalogadjust %>%
filter(Response == 'Latency_aggressive_touch_s')
trace <- mcmc_plot(data$model_chosen_gammalogadjust[[1]], type='trace')
acf_bar <- mcmc_plot(data$model_chosen_gammalogadjust[[1]], type='acf_bar')
rhat_hist <- mcmc_plot(data$model_chosen_gammalogadjust[[1]], type='rhat_hist')
neff_hist <- mcmc_plot(data$model_chosen_gammalogadjust[[1]], type='neff_hist')
plot_diag <- (trace + acf_bar) / (rhat_hist + neff_hist) + plot_annotation(title = 'Latency_aggressive_touch_s')
pdf("05_Modelvalidation_outputs/Gabazine_diag_chain_06_gammalogadjust.pdf", width = 32, height = 24)
plot_diag
dev.off()
#Posterior probability checks
whole <- pp_check(data$model_chosen_gammalogadjust[[1]])
co2 <- pp_check(data$model_chosen_gammalogadjust[[1]], type = "violin_grouped", group = "CO2")
drug <- pp_check(data$model_chosen_gammalogadjust[[1]], type = "violin_grouped", group = "Drug")
plot <- whole + co2 + drug + plot_annotation(title = 'Latency_aggressive_touch_s')
pdf("05_Modelvalidation_outputs/Gabazine_diag_pp_06_gammalogadjust.pdf", width = 32, height = 24)
plot
dev.off()
#DHARMa residuals
preds <- posterior_predict(data$model_chosen_gammalogadjust[[1]], nsamples=250, summary=FALSE)
data2 <- data$data %>% as.data.frame()
resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = data2$Response_measurement,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE,
seed = 926330)
pdf("05_Modelvalidation_outputs/Gabazine_diag_resid_06_gammalogadjust.pdf")
plot(resids)
dev.off()#1st round: binomial for binomial data, poisson for count data, gaussian for continuous data
#Keep binomial Soft_touch_binomial, Aggressive_touch_binomial
#Keep continuous Active_time_s, Dist, Vel, Time_A_s_no_0
load(file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_01_binompoisgaus_p.RData')
#2nd round: negative binomial for count data
#Keep soft_touch_no._no_0, aggressive_touch_no._no_0
load(file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_03_NB_p.RData')
#3rd round: Gamma, log for continuous data which gaussian, identity is no good for
#Keep latency_soft_touch
load(file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_04_gammalog_p.RData')
#4th round: re-run negative binomial with parameters changed for those round 2 was no good for
#Keep No._visits_A
load(file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_05_NBadjust_p.RData')
#5th round: re-run gamma log with parameters changed for those round 3 was no good for
#Keep latency_aggressive_touch_s
load(file = '04_Variableselection_outputs/Picrotoxindata_response_chosen_model_06_gammalogadjust_p.RData')
#Tidy data
Picrotoxindata_response_chosen_model_01_binompoisgaus_p <- Picrotoxindata_response_chosen_model_01_binompoisgaus_p %>%
filter(Response %in% c('Soft_touch_binomial', 'Aggressive_touch_binomial', 'Active_time_s', 'Dist', 'Vel', 'Time_A_s_no_0'))
Picrotoxindata_response_chosen_model_03_NB_p <- Picrotoxindata_response_chosen_model_03_NB_p %>%
filter(Response %in% c('Soft_touch_no._no_0', 'Aggressive_touch_no._no_0'))
Picrotoxindata_response_chosen_model_04_gammalog_p <- Picrotoxindata_response_chosen_model_04_gammalog_p %>%
filter(Response == 'Latency_soft_touch_s')
Picrotoxindata_response_chosen_model_05_NBadjust_p
Picrotoxindata_response_chosen_model_06_gammalogadjust_p <- Picrotoxindata_response_chosen_model_06_gammalogadjust_p %>%
filter(Response == 'Latency_aggressive_touch_s')
Picrotoxindata_final_models <- bind_rows(Picrotoxindata_response_chosen_model_01_binompoisgaus_p, Picrotoxindata_response_chosen_model_03_NB_p, Picrotoxindata_response_chosen_model_04_gammalog_p, Picrotoxindata_response_chosen_model_05_NBadjust_p, Picrotoxindata_response_chosen_model_06_gammalogadjust_p)
save(Picrotoxindata_final_models, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models.RData')Calculate effect sizes for each comparison of interest, including prob of an effect, prob > 20% effect and 80% and 95% HPDI intervals Use exp() on models with log link = fold change Use exp() on models with logit link = odds ratio Linear models manually calculate fold change
#Make a dataset of effect sizes for each contrast = effect_data for each response variable and save in new column
#Models with log link use exp() to calculate fold change, models with gaussian, identity manually calculate fold change (Active_time_s, Dist, Vel, Time_A_s_no_0)
effect_data_function <- function(response) {
data <- Picrotoxindata_final_models %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Active_time|^Dist|^Vel|^Time_A_s_no_0', response_var)) {
co2_effect_draws <- emmeans(data$model_chosen[[1]], ~ CO2|Drug, type='link') %>%
gather_emmeans_draws() %>%
spread(key = CO2, value = .value) %>%
mutate(fold_change = Elevated/Ambient) %>%
select(-Ambient, -Elevated)
co2_effect_emmeans <- co2_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Elevated - Ambient') %>%
filter(Drug=='Sham') %>%
droplevels()
if (co2_effect_emmeans$fold_change[[1]] < 1) {
co2_effect_prob <- co2_effect_draws %>%
filter(Drug=='Sham') %>%
group_by(fold_change, Drug) %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
co2_effect_prob <- co2_effect_draws %>%
filter(Drug=='Sham') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
co2_effect_data <- co2_effect_prob %>% full_join(co2_effect_emmeans)
drug_effect_draws <- emmeans(data$model_chosen[[1]], ~ Drug|CO2, type='link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(fold_change = Picrotoxin/Sham) %>%
select(-Picrotoxin, -Sham)
drug_effect_emmeans_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Picrotoxin - Sham') %>%
droplevels()
if (drug_effect_emmeans_ambient$fold_change[[1]] < 1) {
drug_effect_prob_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_ambient <- drug_effect_draws %>%
filter(CO2=='Ambient') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_ambient <- drug_effect_prob_ambient %>% full_join(drug_effect_emmeans_ambient)
drug_effect_emmeans_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Picrotoxin - Sham') %>%
droplevels()
if (drug_effect_emmeans_elevated$fold_change[[1]] < 1) {
drug_effect_prob_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
summarize(Perc = (sum(fold_change<1)/n())*100,
Perc20 = (sum(fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_elevated <- drug_effect_draws %>%
filter(CO2=='Elevated') %>%
summarize(Perc = (sum(fold_change>1)/n())*100,
Perc20 = (sum(fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_elevated <- drug_effect_prob_elevated %>% full_join(drug_effect_emmeans_elevated)
drug_effect_data <- drug_effect_data_ambient %>% full_join(drug_effect_data_elevated)
interaction_effect_draws <- emmeans(data$model_chosen[[1]], ~Drug|CO2, type = 'link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(Picro_effect_fold_change = Picrotoxin / Sham) %>%
select(-Picrotoxin, -Sham) %>%
spread(key=CO2, value=Picro_effect_fold_change) %>%
mutate(Picro_effect_across_co2_fold_change = Elevated / Ambient) %>%
select(-Ambient, -Elevated)
interaction_effect_emmeans <- interaction_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Interaction')
if (interaction_effect_emmeans$Picro_effect_across_co2_fold_change[[1]] < 1) {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(Picro_effect_across_co2_fold_change<1)/n())*100,
Perc20 = (sum(Picro_effect_across_co2_fold_change<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
} else {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(Picro_effect_across_co2_fold_change>1)/n())*100,
Perc20 = (sum(Picro_effect_across_co2_fold_change>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
}
interaction_effect_data = interaction_effect_prob %>% full_join(interaction_effect_emmeans)
effect_data <- full_join(co2_effect_data, drug_effect_data) %>%
full_join(., interaction_effect_data)
}
else {
co2_effect_draws <- emmeans(data$model_chosen[[1]], revpairwise ~ CO2|Drug, type='link')$contrast %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))
co2_effect_emmeans <- co2_effect_draws %>%
group_by(contrast, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
dplyr::filter(Drug=='Sham') %>%
droplevels()
if (co2_effect_emmeans$exp.value[[1]] < 1) {
co2_effect_prob <- co2_effect_draws %>%
dplyr::filter(Drug=='Sham') %>%
group_by(contrast, Drug) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
co2_effect_prob <- co2_effect_draws %>%
dplyr::filter(Drug=='Sham') %>%
group_by(contrast, Drug) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
co2_effect_data <- co2_effect_prob %>% full_join(co2_effect_emmeans)
drug_effect_draws <- emmeans(data$model_chosen[[1]], revpairwise ~ Drug|CO2, type='link')$contrast %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))
drug_effect_emmeans_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
droplevels()
if (drug_effect_emmeans_ambient$exp.value[[1]] < 1) {
drug_effect_prob_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_ambient <- drug_effect_draws %>%
dplyr::filter(CO2=='Ambient') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_ambient <- drug_effect_prob_ambient %>% full_join(drug_effect_emmeans_ambient)
drug_effect_emmeans_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
droplevels()
if (drug_effect_emmeans_elevated$exp.value[[1]] < 1) {
drug_effect_prob_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value<1)/n())*100,
Perc20 = (sum(exp.value<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
} else {
drug_effect_prob_elevated <- drug_effect_draws %>%
dplyr::filter(CO2=='Elevated') %>%
group_by(contrast, CO2) %>%
summarize(Perc = (sum(exp.value>1)/n())*100,
Perc20 = (sum(exp.value>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
droplevels()
}
drug_effect_data_elevated <- drug_effect_prob_elevated %>% full_join(drug_effect_emmeans_elevated)
drug_effect_data <- drug_effect_data_ambient %>% full_join(drug_effect_data_elevated)
interaction_effect_draws <- emmeans(data$model_chosen[[1]], ~Drug|CO2, type = 'link') %>%
gather_emmeans_draws() %>%
spread(key = Drug, value = .value) %>%
mutate(Picro_effect = Picrotoxin - Sham) %>%
dplyr::select(-Picrotoxin, -Sham) %>%
spread(key=CO2, value=Picro_effect) %>%
mutate(Picro_effect_across_co2 = Elevated - Ambient,
exp.Picro_effect_across_co2 = exp(Picro_effect_across_co2)) %>%
select(-Ambient, -Elevated)
interaction_effect_emmeans <- interaction_effect_draws %>%
median_hdci (.width=c(0.8, 0.95)) %>%
mutate(contrast = 'Interaction')
if (interaction_effect_emmeans$exp.Picro_effect_across_co2[[1]] < 1) {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(exp.Picro_effect_across_co2<1)/n())*100,
Perc20 = (sum(exp.Picro_effect_across_co2<0.8)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
} else {
interaction_effect_prob <- interaction_effect_draws %>%
summarize(Perc = (sum(exp.Picro_effect_across_co2>1)/n())*100,
Perc20 = (sum(exp.Picro_effect_across_co2>1.2)/n())*100) %>%
mutate(Perc = round(Perc, 1),
Perc20 = round(Perc20, 1)) %>%
mutate(contrast = 'Interaction') %>%
droplevels()
}
interaction_effect_data = interaction_effect_prob %>% full_join(interaction_effect_emmeans)
effect_data <- full_join(co2_effect_data, drug_effect_data) %>%
full_join(., interaction_effect_data)
}
}
Picrotoxindata_final_models_effect_data <- Picrotoxindata_final_models %>%
mutate(effect_data = map(Response, effect_data_function))
save(Picrotoxindata_final_models_effect_data, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data.RData') #Partial plots, saved in a new column in data frame
#X axis limit defined by maximum upper CI * 1.2 and x axis labels defined for each response variable
#If used gaussian, identity (Active_time_s, Dist, Vel, Time_A_s_no_0) do NOT exponentiate value (no backtransformation needed), if used gamma, log or NB, log use exp for backtransforming, if used binomial, logit use plogis for backtransformation
plot_partial_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Active_time|^Dist|^Vel|^Time_A_s_no_0', response_var)) {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = .value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, (max(newdata$.upper)*1.2)))
} else if (grepl('binomial', response_var)) {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(plogis.value=plogis(.value)) %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(plogis.value=plogis(.value))%>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = plogis.value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, 1))
}
else {
newdata <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value)) %>%
group_by(CO2, Drug) %>%
median_hdci (.width=c(0.8, 0.95)) %>%
ungroup()
newdata2 <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type='link') %>%
gather_emmeans_draws() %>%
mutate(exp.value=exp(.value))%>%
group_by(CO2, Drug) %>%
ungroup()
CO2.labs <- c("Ambient CO2", "Elevated CO2")
names(CO2.labs) <- c("Ambient", "Elevated")
plot_partial <- ggplot() +
stat_halfeye(data = newdata2, aes(x = exp.value, y = fct_rev(Drug)),
point_interval = median_hdci, .width=c(0.8, 0.95),
slab_alpha = 0.4) +
facet_grid(CO2 ~ .,
switch = "y",
labeller = labeller(CO2 = CO2.labs)) +
scale_y_discrete("") +
theme_classic() +
coord_cartesian(xlim = c(0, (max(newdata$exp.value.upper)*1.2)))
}
if (response == 'Time_A_s_no_0') {
plot_partial <- plot_partial + scale_x_continuous('Time in Zone A (s) (no 0)', expand = c(0,0))
}
else if (response == 'Active_time_s') {
plot_partial <- plot_partial + scale_x_continuous('Active time (s)', expand = c(0,0))
}
else if (response == 'Vel') {
plot_partial <- plot_partial + scale_x_continuous('Average speed (cm/s)', expand = c(0,0))
}
else if (response == 'Soft_touch_binomial') {
plot_partial <- plot_partial + scale_x_continuous('Proportion of squid that touched mirror softly', expand = c(0,0))
}
else if (response == 'Aggressive_touch_binomial') {
plot_partial <- plot_partial + scale_x_continuous('Proportion of squid that touched mirror aggressively', expand = c(0,0))
}
else if (response == 'No._visits_A') {
plot_partial <- plot_partial + scale_x_continuous('No. of visits to Zone A', expand = c(0,0))
}
else if (response == 'Aggressive_touch_no._no_0') {
plot_partial <- plot_partial + scale_x_continuous('No. of aggressive mirror touches (no 0)', expand = c(0,0))
}
else if (response == 'Soft_touch_no._no_0') {
plot_partial <- plot_partial + scale_x_continuous('No. of soft mirror touches (no 0)', expand = c(0,0))
}
else if (response == 'Dist') {
plot_partial <- plot_partial + scale_x_continuous('Total distance moved (cm)', expand = c(0,0))
}
else if (response == 'Latency_soft_touch_s') {
plot_partial <- plot_partial + scale_x_continuous('Latency to first soft mirror touch (s)', expand = c(0,0))
}
else if (response == 'Latency_aggressive_touch_s') {
plot_partial <- plot_partial + scale_x_continuous('Latency to first aggressive mirror touch (s)', expand = c(0,0))
}
}
Picrotoxindata_final_models_effect_data_plots <- Picrotoxindata_final_models_effect_data_plots %>%
mutate(plot_partial = map(Response, plot_partial_function))
save(Picrotoxindata_final_models_effect_data_plots, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots.RData')#Show effect sizes as fold change (by using exp() values from models with log link and manually calculating fold change for linear models(Active_time_s, Dist, Vel, Time_A_s_no_0)) - for binomial models effect size is an odds ratio
plot_caterpillar_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == response)
response_var <- data$Response[[1]]
if (grepl('^Active_time|^Dist|^Time_A_s_no_0', response_var)) {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=fold_change, x=CO2,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=fold_change, x=Drug,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=Picro_effect_across_co2_fold_change, x = contrast,
ymin = .lower, ymax = .upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=Picro_effect_across_co2_fold_change, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(y=fold_change, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=fold_change, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=Picro_effect_across_co2_fold_change, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', trans=scales::log2_trans()) +
coord_flip() +
theme_classic()
}
else if (grepl('^Vel', response_var)) {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=fold_change, x=CO2,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=fold_change, x=Drug,
ymin=.lower, ymax=.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=fold_change, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=Picro_effect_across_co2_fold_change, x = contrast,
ymin = .lower, ymax = .upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=Picro_effect_across_co2_fold_change, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(y=fold_change, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=fold_change, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=Picro_effect_across_co2_fold_change, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', lim = c(0.499,2.05), breaks = c(0, 0.5, 1, 2), trans=scales::log2_trans(), labels = scales::number_format(accuracy = 0.1)) +
coord_flip() +
theme_classic()
}
else {
plot_caterpillar <- ggplot(data = data$effect_data[[1]]) +
geom_hline(yintercept=1, linetype='dashed', alpha = 0.5) +
geom_linerange(aes(y=exp.value, x=CO2,
ymin=exp.value.lower, ymax=exp.value.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.value, x=CO2),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=exp.value, x=Drug,
ymin=exp.value.lower, ymax=exp.value.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.value, x=Drug),
size = 3,
show.legend = FALSE) +
geom_linerange(aes(y=exp.Picro_effect_across_co2, x = contrast,
ymin = exp.Picro_effect_across_co2.lower, ymax = exp.Picro_effect_across_co2.upper,
size = factor(.width)),
show.legend = FALSE) +
geom_point(aes(y=exp.Picro_effect_across_co2, x=contrast),
size = 3,
show.legend = FALSE) +
geom_text(aes(exp.value, x = CO2, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=exp.value, x = Drug, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
geom_text(aes(y=exp.Picro_effect_across_co2, x = contrast, label = Perc, vjust = 0), nudge_x = 0.2, size = (5/14) * 8) +
scale_size_manual(values=c(1, 0.5)) +
scale_x_discrete('', labels = c('Sham' = expression(atop("",
atop(textstyle("CO"[2]~"effect"),
atop(textstyle(sham))))),
'Interaction' = expression(atop("",
atop(textstyle(""),
atop(textstyle(""),
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]),
atop(textstyle("vs"),
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))))))),
'Elevated' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("elevated CO"[2]))))),
'Ambient' = expression(atop("",
atop(textstyle("Picro effect"),
atop(textstyle("ambient CO"[2])))))),
limits = c('Interaction', 'Elevated', 'Ambient', 'Sham')) +
scale_y_continuous('Effect size', trans=scales::log2_trans()) +
coord_flip() +
theme_classic()
}
}
Picrotoxindata_final_models_effect_data_plots <- Picrotoxindata_final_models_effect_data_plots %>%
mutate(plot_caterpillar = map(Response, plot_caterpillar_function))
save(Picrotoxindata_final_models_effect_data_plots, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots.RData') Look at partial and caterpillar plots created alongside quick overview partial and caterpillar plots made automatically to check if calculations were done correctly
#Partial plot and caterpillar plot alongside overview plots to check look like expected
#overview cat plot is a good comparison, but wont match exactly as CIs are with quantiles, and plotted on the link scale (compared to HPDI intervals, and on response scale)
check_plots_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == response)
plot_partial <- data$plot_partial[[1]]
plot_caterpillar <- data$plot_caterpillar[[1]]
overview_partial <- ggemmeans(data$model_chosen[[1]], terms = ~CO2 * Drug) %>% plot
overview_caterpillar <- mcmc_intervals(data$model_chosen[[1]], pars = c("b_CO2Elevated", "b_DrugPicrotoxin", "b_CO2Elevated:DrugPicrotoxin"), prob = 0.8, prob_outer = 0.95, point_est = "median")
check_plots <- (overview_partial + overview_caterpillar) / (plot_partial + plot_caterpillar)
ggsave(check_plots, file = paste("06_Modelinvestigation_outputs/check_plots_", response, ".pdf", sep = ""), width = 40, height = 20, unit = "cm", device = cairo_pdf)
}
map(Picrotoxindata_final_models_effect_data_plots$Response, check_plots_function)#Sample size each treatment group, for each response variable
table_n_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == response)
data1 <- data$data
data1 %>%
as.data.frame() %>%
group_by(CO2, Drug) %>%
summarise(N = n()) %>%
kable(caption = response)
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots %>%
mutate(table_n = map(Response, table_n_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')
#Estimate +- 95% CI for each treatment group, across all response variables in a new column (table version of partial plot)
#On response scale
table_estimates_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Active_time_s|Dist|^Vel', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB', model_name)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = prob,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('gamma', model_name)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = response,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('^Time_A_s_no_0|Active_time_s|Dist|^Vel', response_var)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = emmean,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
else if (grepl('binomial', response_var)) {
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
table_estimates <- emmeans(data$model_chosen[[1]], ~CO2*Drug, type = 'response') %>%
as.data.frame() %>%
rename(Estimate = response,
'Lower HPDI' = lower.HPD,
'Upper HPDI' = upper.HPD) %>%
arrange(CO2) %>%
kable(caption = caption, digits = 2)
}
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots_tables %>%
mutate(table_estimates = map(Response, table_estimates_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')
#Effect sizes +- 95% CI for each contrast, across all response variables in a new column (table version of caterpillar plot)
#exp for all models to show fold-change/odds ratio. Fold changes calculated manually for linear models
table_contrasts_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Active_time|^Dist|^Vel|^Time_A_s_no_0', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
caption <- paste(caption_response, " ~ CO2 * Drug", caption1, caption2, sep = "")
options(knitr.kable.NA = '-')
effect_data <- data$effect_data[[1]]
if (grepl('^Active_time|^Dist|^Vel|^Time_A_s_no_0', response_var)) {
effect_data %>%
as.data.frame() %>%
filter(.width == 0.95) %>%
select(contrast, Drug, CO2, fold_change, .lower, .upper, Picro_effect_across_co2_fold_change, Perc, Perc20) %>%
mutate(Treatment = coalesce(Drug, CO2),
Estimate = coalesce(fold_change, Picro_effect_across_co2_fold_change),
'Lower HPDI' = .lower,
'Upper HPDI' = .upper) %>%
select(contrast, Treatment, Estimate, 'Lower HPDI', 'Upper HPDI', Perc, Perc20) %>%
rename(Contrast = contrast,
Prob = Perc,
'Prob(20%)' = Perc20) %>%
kable(caption = caption, digits = 2)
}
else {
effect_data %>%
as.data.frame() %>%
filter(.width == 0.95) %>%
select(contrast, Drug, CO2, exp.value, exp.value.lower, exp.value.upper, exp.Picro_effect_across_co2, exp.Picro_effect_across_co2.lower, exp.Picro_effect_across_co2.upper, Perc, Perc20) %>%
mutate(Treatment = coalesce(Drug, CO2),
Estimate = coalesce(exp.value, exp.Picro_effect_across_co2),
'Lower HPDI' = coalesce(exp.value.lower, exp.Picro_effect_across_co2.lower),
'Upper HPDI' = coalesce(exp.value.upper, exp.Picro_effect_across_co2.upper)) %>%
select(contrast, Treatment, Estimate, 'Lower HPDI', 'Upper HPDI', Perc, Perc20) %>%
rename(Contrast = contrast,
Prob = Perc,
'Prob(20%)' = Perc20) %>%
kable(caption = caption, digits = 2)
}
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots_tables %>%
mutate(table_contrasts = map(Response, table_contrasts_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')
#Output for each table on link scale, including other explanatory variables included
table_summary_link_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption1 <- ", "
} else if (grepl('size', model_name)) {
caption1 <- " + Mantle length, "
} else if (grepl('method', model_name)) {
caption1 <- " + Behavioural tank + Time of test, "
} else if (grepl('_drug', model_name)) {
caption1 <- " + Drug test number, "
} else if (grepl('housing', model_name)) {
caption1 <- " + Number of acclimation days + Date introduced to treatment tank, "
} else if (grepl('system', model_name)) {
caption1 <- " + System, "
}
if (grepl('NB', model_name)) {
caption2 <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption2 <- "Gamma (log)"
} else if (grepl('^Active_time|^Dist|^Vel|^Time_A_s_no_0', response_var)) {
caption2 <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption2 <- "Binomial (logit)"
}
caption <- paste(response, " ~ CO2", " * Drug", caption1, caption2, sep = "")
model <- data$model_chosen[[1]]
table_summary_link <- model$fit %>%
tidyMCMC(estimate.method = 'median',
conf.int = TRUE, conf.method = 'HPDinterval') %>%
rename(Term = term,
Estimate = estimate,
SE = std.error,
'Lower HPDI' = conf.low,
'Upper HPDI' = conf.high) %>%
kable(caption = caption, digits = 2)
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots_tables %>%
mutate(table_summary_link = map(Response, table_summary_link_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')
#Priors
table_priors_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB|gamma', model_name)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope', shape)
} else if (grepl('binomial', response_var)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope')
} else {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response) %>%
rename('Slope' = b) %>%
select('Response variable', Intercept, 'Slope', sigma)
}
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots_tables %>%
mutate(table_priors = map(Response, table_priors_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')
#Model info summary (response var, explanatory vars, family and link) as well as priors used
table_modelinfo_priors_function <- function(response) {
data <- Picrotoxindata_final_models_effect_data_plots_tables %>%
filter(Response == response)
model_name <- data$model_chosen_name[[1]]
response_var <- data$Response[[1]]
if (grepl('co2drug', model_name)) {
caption_explanatory <- "CO2 * Drug"
} else if (grepl('size', model_name)) {
caption_explanatory <- "CO2 * Drug + Mantle length"
} else if (grepl('method', model_name)) {
caption_explanatory <- "CO2 * Drug + Behavioural tank + Time of test"
} else if (grepl('_drug', model_name)) {
caption_explanatory <- "CO2 * Drug + Drug test number"
} else if (grepl('housing', model_name)) {
caption_explanatory <- "CO2 * Drug + Number of acclimation days + Date introduced to treatment tank"
} else if (grepl('system', model_name)) {
caption_explanatory <- "CO2 * Drug + System"
}
if (grepl('NB', model_name)) {
caption_family_link <- "Negative binomial (log)"
} else if (grepl('gamma', model_name)) {
caption_family_link <- "Gamma (log)"
} else if (grepl('^Time_A_s_no_0|^Active_time_s|^Vel', response_var)) {
caption_family_link <- "Gaussian (identity)"
} else if (grepl('binomial', response_var)) {
caption_family_link <- "Binomial (logit)"
}
if (grepl('Time_A_s_no_0', response_var)) {
caption_response <- 'Time in Zone A (s)'
} else if (grepl('No._visits_A$', response_var)) {
caption_response <- 'No. of visits to Zone A'
} else if (grepl('Soft_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror softly'
} else if (grepl('Latency_soft_touch_s', response_var)) {
caption_response <- 'Latency to first soft mirror touch (s)'
} else if (grepl('Soft_touch_no._no_0', response_var)) {
caption_response <- 'No. of soft mirror touches'
} else if (grepl('Aggressive_touch_binomial', response_var)) {
caption_response <- 'Proportion of squid that touched mirror aggressively'
} else if (grepl('Latency_aggressive_touch_s', response_var)) {
caption_response <- 'Latency to first aggressive mirror touch (s)'
} else if (grepl('Aggressive_touch_no._no_0', response_var)) {
caption_response <- 'No. of aggressive mirror touches'
} else if (grepl('Active_time_s', response_var)) {
caption_response <- 'Active time (s)'
} else if (grepl('Dist', response_var)) {
caption_response <- 'Total distance moved (cm)'
} else if (grepl('Vel', response_var)) {
caption_response <- 'Average speed (cm/s)'
} else {
caption_response <- response
}
if (grepl('NB|gamma', model_name)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept,
'shape prior' = shape) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior', 'shape prior')
} else if (grepl('binomial', response_var)) {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior')
} else {
df <- prior_summary(data$model_chosen[[1]]) %>%
as.data.frame()
df$prior <- sub("^$", "improper uniform", df$prior)
df %>%
filter(coef == '') %>%
select(prior, class) %>%
spread(class, prior) %>%
mutate('Response variable' = caption_response,
'Explanatory variables' = caption_explanatory,
'Family (link)' = caption_family_link) %>%
rename('Slope prior' = b,
'Intercept prior' = Intercept,
'sigma prior' = sigma) %>%
select('Response variable', 'Explanatory variables', 'Family (link)', 'Intercept prior', 'Slope prior', 'sigma prior')
}
}
Picrotoxindata_final_models_effect_data_plots_tables <- Picrotoxindata_final_models_effect_data_plots_tables %>%
mutate(table_modelinfo_priors = map(Response, table_modelinfo_priors_function))
save(Picrotoxindata_final_models_effect_data_plots_tables, file = '06_Modelinvestigation_outputs/Picrotoxindata_final_models_effect_data_plots_tables.RData')Put together figures for gabazine and picrotoxin experiment results, grouped by behaviour
theme.size = 8
axis.text.size = 7
geom.text.size = (5/14) * axis.text.size
My_Theme_Partial = theme(
axis.title = element_text(size = theme.size, family = "Arial", colour = "black", face = "bold"),
axis.text = element_text(size = axis.text.size, family = "Arial", colour = "black"),
strip.text = element_text(size = theme.size, family = "Arial", colour = "black", face = "bold"),
plot.subtitle = element_text(size = 12, family = "Arial", colour = "black", face = "bold"))
My_Theme_Caterpillar = theme(
axis.title.x = element_text(size = theme.size, family = "Arial", colour = "black", face = "bold"),
axis.text = element_text(size = axis.text.size, family = "Arial", colour = "black"),
plot.subtitle = element_text(size = 12, family = "Arial", colour = "black", face = "bold"))Go_data_Time_A_s_no_0 <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Time_A_s_no_0')
Go_text_Time_A_s_no_0 <- data.frame(
label = c("n=12", "n=12", "n=11", "n=16"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Time_A_s_no_0 <- Go_data_Time_A_s_no_0$plot_partial[[1]] +
labs(subtitle = "A") +
geom_text(data = Go_text_Time_A_s_no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('Time in Zone A (s)', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Time_A_s_no_0 <- Go_data_Time_A_s_no_0$plot_caterpillar[[1]] +
labs(subtitle = "B") +
scale_y_continuous('Effect size', lim = c(0.069,4.05), breaks = c(0.125, 0.25, 0.5, 1, 2, 4), labels = c('', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_No._visits_A <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'No._visits_A')
Go_text_No._visits_A <- data.frame(
label = c("n=12", "n=12", "n=13", "n=16"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_No._visits_A <- Go_data_No._visits_A$plot_partial[[1]] +
labs(subtitle = "C") +
geom_text(data = Go_text_No._visits_A, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
theme(axis.text.y = element_text(angle = 90, hjust = 0.5)) +
My_Theme_Partial
Go_plot_caterpillar_No._visits_A <- Go_data_No._visits_A$plot_caterpillar[[1]] +
labs(subtitle = "D") +
scale_y_continuous('Effect size', lim = c(0.059,16), breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16), labels = c('', '0.25', '', '1', '', '4', '', '16'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Time_A_s_no_0 <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Time_A_s_no_0')
P_text_Time_A_s_no_0 <- data.frame(
label = c("n=16", "n=18", "n=21", "n=21"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Time_A_s_no_0 <- P_data_Time_A_s_no_0$plot_partial[[1]] +
labs(subtitle = "E") +
geom_text(data = P_text_Time_A_s_no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90, hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('Time in Zone A (s)', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Time_A_s_no_0 <- P_data_Time_A_s_no_0$plot_caterpillar[[1]] +
labs(subtitle = "F") +
scale_y_continuous('Effect size', breaks = c(0.5, 1, 2), labels = c('0.5', '1', '2'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_No._visits_A <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'No._visits_A')
P_text_No._visits_A <- data.frame(
label = c("n=19", "n=18", "n=21", "n=22"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_No._visits_A <- P_data_No._visits_A$plot_partial[[1]] +
labs(subtitle = "G") +
geom_text(data = P_text_No._visits_A, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
theme(axis.text.y = element_text(angle = 90, hjust = 0.5)) +
My_Theme_Partial
P_plot_caterpillar_No._visits_A <- P_data_No._visits_A$plot_caterpillar[[1]] +
labs(subtitle = "H") +
scale_y_continuous('Effect size', breaks = c(0.25, 0.5, 1, 2), labels = c('0.25', '0.5', '1', '2'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
figure_space_use <- (wrap_elements(grid::textGrob(expression(bold("Gabazine experiment")~"(GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((Go_plot_partial_Time_A_s_no_0 +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(Go_plot_caterpillar_Time_A_s_no_0 +
theme(plot.margin = unit(c(0,30,5,0), "pt")))) /
((Go_plot_partial_No._visits_A +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(Go_plot_caterpillar_No._visits_A +
theme(plot.margin = unit(c(0,30,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2), 'cm'), widths = unit(9, "cm"))) |
(wrap_elements(grid::textGrob(expression(bold("Picrotoxin experiment")~"(non-specific GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((P_plot_partial_Time_A_s_no_0 +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(P_plot_caterpillar_Time_A_s_no_0 +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_No._visits_A +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(P_plot_caterpillar_No._visits_A +
theme(plot.margin = unit(c(0,10,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2), 'cm'), widths = unit(9, "cm")))
ggsave(figure_space_use, file = paste("07_Summary_figures_Gabazineold_Picrotoxin_outputs/figure_space_use2.pdf"), width = 21, height = 13, unit = "cm", device = cairo_pdf)Go_data_Soft_touch_binomial <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Soft_touch_binomial')
Go_text_Soft_touch_binomial <- data.frame(
label = c("n=22", "n=24", "n=22", "n=23"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Soft_touch_binomial <- Go_data_Soft_touch_binomial$plot_partial[[1]] +
labs(subtitle = "A") +
geom_text(data = Go_text_Soft_touch_binomial, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('Proportion of squid that \ntouched mirror softly', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Soft_touch_binomial <- Go_data_Soft_touch_binomial$plot_caterpillar[[1]] +
labs(subtitle = "B") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4), labels = c('', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Latency_soft_touch_s <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Latency_soft_touch_s')
Go_text_Latency_soft_touch_s <- data.frame(
label = c("n=13", "n=11", "n=12", "n=11"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Latency_soft_touch_s <- Go_data_Latency_soft_touch_s$plot_partial[[1]] +
labs(subtitle = "C") +
geom_text(data = Go_text_Latency_soft_touch_s, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('Latency to soft mirror touch', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Latency_soft_touch_s <- Go_data_Latency_soft_touch_s$plot_caterpillar[[1]] +
labs(subtitle = "D") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4), labels = c('', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Soft_touch_no._no_0 <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Soft_touch_no._no_0')
Go_text_Soft_touch_no._no_0 <- data.frame(
label = c("n=13", "n=11", "n=12", "n=11"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Soft_touch_no._no_0 <- Go_data_Soft_touch_no._no_0$plot_partial[[1]] +
labs(subtitle = "E") +
geom_text(data = Go_text_Soft_touch_no._no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('No. soft mirror touches', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Soft_touch_no._no_0 <- Go_data_Soft_touch_no._no_0$plot_caterpillar[[1]] +
labs(subtitle = "F") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8), labels = c('', '0.25', '', '1', '', '4', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Soft_touch_binomial <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Soft_touch_binomial')
P_text_Soft_touch_binomial <- data.frame(
label = c("n=27", "n=26", "n=26", "n=26"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Soft_touch_binomial <- P_data_Soft_touch_binomial$plot_partial[[1]] +
labs(subtitle = "G") +
geom_text(data = P_text_Soft_touch_binomial, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('Proportion of squid that \ntouched mirror softly', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Soft_touch_binomial <- P_data_Soft_touch_binomial$plot_caterpillar[[1]] +
labs(subtitle = "H") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8), labels = c('', '0.25', '', '1', '', '4', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Latency_soft_touch_s <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Latency_soft_touch_s')
P_text_Latency_soft_touch_s <- data.frame(
label = c("n=8", "n=14", "n=13", "n=18"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Latency_soft_touch_s <- P_data_Latency_soft_touch_s$plot_partial[[1]] +
labs(subtitle = "I") +
geom_text(data = P_text_Latency_soft_touch_s, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('Latency to soft mirror touch', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Latency_soft_touch_s <- P_data_Latency_soft_touch_s$plot_caterpillar[[1]] +
labs(subtitle = "J") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16), labels = c('', '0.25', '', '1', '', '4', '', '16'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Soft_touch_no._no_0 <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Soft_touch_no._no_0')
P_text_Soft_touch_no._no_0 <- data.frame(
label = c("n=8", "n=14", "n=13", "n=18"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Soft_touch_no._no_0 <- P_data_Soft_touch_no._no_0$plot_partial[[1]] +
labs(subtitle = "K") +
geom_text(data = P_text_Soft_touch_no._no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('No. soft mirror touches', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Soft_touch_no._no_0 <- P_data_Soft_touch_no._no_0$plot_caterpillar[[1]] +
labs(subtitle = "L") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4), labels = c('', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
figure_soft_mirror_touch <- (wrap_elements(grid::textGrob(expression(bold("Gabazine experiment")~"(GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((Go_plot_partial_Soft_touch_binomial +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(Go_plot_caterpillar_Soft_touch_binomial +
theme(plot.margin = unit(c(0,30,5,0), "pt")))) /
((Go_plot_partial_Latency_soft_touch_s +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(Go_plot_caterpillar_Latency_soft_touch_s +
theme(plot.margin = unit(c(0,30,5,0), "pt")))) /
((Go_plot_partial_Soft_touch_no._no_0 +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(Go_plot_caterpillar_Soft_touch_no._no_0 +
theme(plot.margin = unit(c(0,30,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2, 4.2), 'cm'), widths = unit(9, "cm"))) |
(wrap_elements(grid::textGrob(expression(bold("Picrotoxin experiment")~"(non-specific GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((P_plot_partial_Soft_touch_binomial +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(P_plot_caterpillar_Soft_touch_binomial +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Latency_soft_touch_s +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(P_plot_caterpillar_Latency_soft_touch_s +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Soft_touch_no._no_0 +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(P_plot_caterpillar_Soft_touch_no._no_0 +
theme(plot.margin = unit(c(0,10,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2, 4.2), 'cm'), widths = unit(9, "cm")))
ggsave(figure_soft_mirror_touch, file = paste("07_Summary_figures_Gabazineold_Picrotoxin_outputs/figure_soft_mirror_touch2.pdf"), width = 21, height = 18, unit = "cm", device = cairo_pdf)Go_data_Aggressive_touch_binomial <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Aggressive_touch_binomial')
Go_text_Aggressive_touch_binomial <- data.frame(
label = c("n=22", "n=24", "n=22", "n=23"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Aggressive_touch_binomial <- Go_data_Aggressive_touch_binomial$plot_partial[[1]] +
labs(subtitle = "A") +
geom_text(data = Go_text_Aggressive_touch_binomial, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('Proportion of squid that \ntouched mirror aggressively', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Aggressive_touch_binomial <- Go_data_Aggressive_touch_binomial$plot_caterpillar[[1]] +
labs(subtitle = "B") +
scale_y_continuous('Effect size', breaks = c(0.0625, 0.125, 0.25, 0.5, 1, 2, 4), labels = c('0.0625', '', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Latency_aggressive_touch_s <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Latency_aggressive_touch_s')
Go_text_Latency_aggressive_touch_s <- data.frame(
label = c("n=13", "n=11", "n=12", "n=11"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Latency_aggressive_touch_s <- Go_data_Latency_aggressive_touch_s$plot_partial[[1]] +
labs(subtitle = "C") +
geom_text(data = Go_text_Latency_aggressive_touch_s, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('Latency to aggressive \nmirror touch', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Latency_aggressive_touch_s <- Go_data_Latency_aggressive_touch_s$plot_caterpillar[[1]] +
labs(subtitle = "D") +
scale_y_continuous('Effect size', breaks = c(0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16), labels = c('0.03125', '', '', '0.25', '', '1', '', '4', '', '16'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Aggressive_touch_no._no_0 <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Aggressive_touch_no._no_0')
Go_text_Aggressive_touch_no._no_0 <- data.frame(
label = c("n=13", "n=11", "n=12", "n=11"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_plot_partial_Aggressive_touch_no._no_0 <- Go_data_Aggressive_touch_no._no_0$plot_partial[[1]] +
labs(subtitle = "E") +
geom_text(data = Go_text_Aggressive_touch_no._no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
scale_x_continuous('No. aggressive mirror touches', expand = c(0,0)) +
My_Theme_Partial
Go_plot_caterpillar_Aggressive_touch_no._no_0 <- Go_data_Aggressive_touch_no._no_0$plot_caterpillar[[1]] +
labs(subtitle = "F") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4), labels = c('', '0.25', '', '1', '', '4'), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Aggressive_touch_binomial <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Aggressive_touch_binomial')
P_text_Aggressive_touch_binomial <- data.frame(
label = c("n=27", "n=26", "n=26", "n=26"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Aggressive_touch_binomial <- P_data_Aggressive_touch_binomial$plot_partial[[1]] +
labs(subtitle = "G") +
geom_text(data = P_text_Aggressive_touch_binomial, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('Proportion of squid that \ntouched mirror aggressively', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Aggressive_touch_binomial <- P_data_Aggressive_touch_binomial$plot_caterpillar[[1]] +
labs(subtitle = "H") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8), labels = c('', '0.25', '', '1', '', '4', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Latency_aggressive_touch_s <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Latency_aggressive_touch_s')
P_text_Latency_aggressive_touch_s <- data.frame(
label = c("n=8", "n=14", "n=13", "n=18"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Latency_aggressive_touch_s <- P_data_Latency_aggressive_touch_s$plot_partial[[1]] +
labs(subtitle = "I") +
geom_text(data = P_text_Latency_aggressive_touch_s, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('Latency to aggressive \nmirror touch', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Latency_aggressive_touch_s <- P_data_Latency_aggressive_touch_s$plot_caterpillar[[1]] +
labs(subtitle = "J") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8), labels = c('', '0.25', '', '1', '', '4', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Aggressive_touch_no._no_0 <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Aggressive_touch_no._no_0')
P_text_Aggressive_touch_no._no_0 <- data.frame(
label = c("n=8", "n=14", "n=13", "n=18"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_plot_partial_Aggressive_touch_no._no_0 <- P_data_Aggressive_touch_no._no_0$plot_partial[[1]] +
labs(subtitle = "K") +
geom_text(data = P_text_Aggressive_touch_no._no_0, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
scale_x_continuous('No. aggressive mirror touches', expand = c(0,0)) +
My_Theme_Partial
P_plot_caterpillar_Aggressive_touch_no._no_0 <- P_data_Aggressive_touch_no._no_0$plot_caterpillar[[1]] +
labs(subtitle = "L") +
scale_y_continuous('Effect size', breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8), labels = c('', '0.25', '', '1', '', '4', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
figure_aggressive_mirror_touch <- (wrap_elements(grid::textGrob(expression(bold("Gabazine experiment")~"(GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((Go_plot_partial_Aggressive_touch_binomial +
theme(plot.margin = unit(c(0,15,5,0), "pt"))) |
(Go_plot_caterpillar_Aggressive_touch_binomial +
theme(plot.margin = unit(c(0,20,5,0), "pt")))) /
((Go_plot_partial_Latency_aggressive_touch_s +
theme(plot.margin = unit(c(0,15,5,0), "pt"))) |
(Go_plot_caterpillar_Latency_aggressive_touch_s +
theme(plot.margin = unit(c(0,20,5,0), "pt")))) /
((Go_plot_partial_Aggressive_touch_no._no_0 +
theme(plot.margin = unit(c(0,15,0,0), "pt"))) |
(Go_plot_caterpillar_Aggressive_touch_no._no_0 +
theme(plot.margin = unit(c(0,20,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.1, 4.1, 4.1), 'cm'), widths = unit(9, "cm"))) |
(wrap_elements(grid::textGrob(expression(bold("Picrotoxin experiment")~"(non-specific GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((P_plot_partial_Aggressive_touch_binomial +
theme(plot.margin = unit(c(0,15,5,0), "pt"))) |
(P_plot_caterpillar_Aggressive_touch_binomial +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Latency_aggressive_touch_s +
theme(plot.margin = unit(c(0,15,5,0), "pt"))) |
(P_plot_caterpillar_Latency_aggressive_touch_s +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Aggressive_touch_no._no_0 +
theme(plot.margin = unit(c(0,15,0,0), "pt"))) |
(P_plot_caterpillar_Aggressive_touch_no._no_0 +
theme(plot.margin = unit(c(0,10,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.1, 4.1, 4.1), 'cm'), widths = unit(9, "cm")))
ggsave(figure_aggressive_mirror_touch, file = paste("07_Summary_figures_Gabazineold_Picrotoxin_outputs/figure_aggressive_mirror_touch2.pdf"), width = 21, height = 18, unit = "cm", device = cairo_pdf)Go_text_Activity <- data.frame(
label = c("n=12", "n=12", "n=13", "n=16"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Gabazine", "Sham", "Gabazine")
)
Go_data_Active_time_s <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Active_time_s')
Go_plot_partial_Active_time_s <- Go_data_Active_time_s$plot_partial[[1]] +
labs(subtitle = "A") +
geom_text(data = Go_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
My_Theme_Partial
Go_plot_caterpillar_Active_time_s <- Go_data_Active_time_s$plot_caterpillar[[1]] +
labs(subtitle = "B") +
scale_y_continuous('Effect size', breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Dist <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Dist')
Go_plot_partial_Dist <- Go_data_Dist$plot_partial[[1]] +
labs(subtitle = "C") +
geom_text(data = Go_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
My_Theme_Partial
Go_plot_caterpillar_Dist <- Go_data_Dist$plot_caterpillar[[1]] +
labs(subtitle = "D") +
scale_y_continuous('Effect size', breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
Go_data_Vel <- Gabazinedata_old_final_models_effect_data_plots %>%
filter(Response == 'Vel')
Go_plot_partial_Vel <- Go_data_Vel$plot_partial[[1]] +
labs(subtitle = "E") +
geom_text(data = Go_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Gab", "Sham")) +
My_Theme_Partial
Go_plot_caterpillar_Vel <- Go_data_Vel$plot_caterpillar[[1]] +
labs(subtitle = "F") +
scale_y_continuous('Effect size', limits = c(0.45, 2.05), breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_text_Activity <- data.frame(
label = c("n=19", "n=18", "n=21", "n=22"),
CO2 = c("Ambient", "Ambient", "Elevated", "Elevated"),
x = c(0, 0, 0, 0),
y = c("Sham", "Picrotoxin", "Sham", "Picrotoxin")
)
P_data_Active_time_s <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Active_time_s')
P_plot_partial_Active_time_s <- P_data_Active_time_s$plot_partial[[1]] +
labs(subtitle = "G") +
geom_text(data = P_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
My_Theme_Partial
P_plot_caterpillar_Active_time_s <- P_data_Active_time_s$plot_caterpillar[[1]] +
labs(subtitle = "H") +
scale_y_continuous('Effect size', limits = c(0.45, 2.05), breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Dist <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Dist')
P_plot_partial_Dist <- P_data_Dist$plot_partial[[1]] +
labs(subtitle = "I") +
geom_text(data = P_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
My_Theme_Partial
P_plot_caterpillar_Dist <- P_data_Dist$plot_caterpillar[[1]] +
labs(subtitle = "J") +
scale_y_continuous('Effect size', breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
P_data_Vel <- Picrotoxindata_final_models_effect_data_plots %>%
filter(Response == 'Vel')
P_plot_partial_Vel <- P_data_Vel$plot_partial[[1]] +
labs(subtitle = "K") +
geom_text(data = P_text_Activity, mapping = aes(x = x, y = y, label = label),
size = geom.text.size,
family = "Arial",
colour = "black",
hjust = -0.5, vjust = 1.8) +
theme(axis.text.y = element_text(angle = 90,
hjust = 0.5)) +
scale_y_discrete("", labels=c ("Picro", "Sham")) +
My_Theme_Partial
P_plot_caterpillar_Vel <- P_data_Vel$plot_caterpillar[[1]] +
labs(subtitle = "L") +
scale_y_continuous('Effect size', limits = c(0.45, 2.05), breaks = c(0.25, 0.5, 1, 2, 4), labels = c('', '0.5', '1', '2', ''), trans=scales::log2_trans()) +
My_Theme_Caterpillar
figure_activity <- (wrap_elements(grid::textGrob(expression(bold("Gabazine experiment")~"(GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((Go_plot_partial_Active_time_s +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(Go_plot_caterpillar_Active_time_s +
theme(plot.margin = unit(c(0,20,5,0), "pt")))) /
((Go_plot_partial_Dist +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(Go_plot_caterpillar_Dist +
theme(plot.margin = unit(c(0,20,5,0), "pt")))) /
((Go_plot_partial_Vel +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(Go_plot_caterpillar_Vel +
theme(plot.margin = unit(c(0,20,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2, 4.2), 'cm'), widths = unit(9, "cm"))) |
(wrap_elements(grid::textGrob(expression(bold("Picrotoxin experiment")~"(non-specific GABA"[A]~"R antagonist)"),
gp = gpar(fontsize = 8,
fontfamily = "Arial",
col = "black",
fontface = "bold"))) /
((P_plot_partial_Active_time_s +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(P_plot_caterpillar_Active_time_s +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Dist +
theme(plot.margin = unit(c(0,10,5,0), "pt"))) |
(P_plot_caterpillar_Dist +
theme(plot.margin = unit(c(0,10,5,0), "pt")))) /
((P_plot_partial_Vel +
theme(plot.margin = unit(c(0,10,0,0), "pt"))) |
(P_plot_caterpillar_Vel +
theme(plot.margin = unit(c(0,10,0,0), "pt")))) +
plot_layout(heights = unit(c(0.3, 4.2, 4.2, 4.2), 'cm'), widths = unit(9, "cm")))
ggsave(figure_activity, file = paste("07_Summary_figures_Gabazineold_Picrotoxin_outputs/figure_activity2.pdf"), width = 21, height = 18, unit = "cm", device = cairo_pdf)