https://github.com/charlyknz/HostParasite
Revision 9db9de4a3b6b4cbe79312bb291b5cd3de8b2d947 authored by charlyknz on 31 August 2021, 07:08:23 UTC, committed by GitHub on 31 August 2021, 07:08:23 UTC
1 parent 81ac900
Tip revision: 9db9de4a3b6b4cbe79312bb291b5cd3de8b2d947 authored by charlyknz on 31 August 2021, 07:08:23 UTC
Update README.md
Update README.md
Tip revision: 9db9de4
Host_BayesianAnalysis.R
#---
#R script to analyse parasite and temperature effects on host reproduction
#authors: Andrew Jackson & Charlotte Kunze
#---
#load packages
library(Hmisc)
library(cowplot)
library(tidyverse)
library(ggplot2)
library(scales)
library(rjags)
library(ggpubr)
#set the working directory in the same folder, our Rscript is saved
setwd(dirname(rstudioapi::getSourceEditorContext()$path))
#### import host data ####
host_repro_data <- read.csv2("sum_daphnia_noNA.csv", header = TRUE,
stringsAsFactors = FALSE) %>%
mutate(treat2 = treat)
#check import:
str(host_repro_data)
names(host_repro_data)
#first plot
ggplot(host_repro_data, aes(x = mean_temp, y = sum, group = treat2, color = treat2)) +
geom_point() +
geom_smooth()+
theme_bw()
#### Model formulation ####
#Source the beta function we used to model the TPCs which is held in an external `*.R` file as it is called from multiple notebooks.
source("betaFunction.R")
# Temperature T was replaced in the maths above with x when specifying the model.
model_beta_fun_poisson = '
model
{
# Likelihood
for (i in 1:n) {
y[i] ~ dpois(lambda[i])
# log10 link function
lambda[i] <- 10^(mu[i])
mu[i] <- c_m[treat[i]] *
( (x_max[treat[i]] - x[i]) / (x_max[treat[i]] - x_opt[treat[i]]) ) *
( (x[i] - x_min[treat[i]]) / (x_opt[treat[i]] - x_min[treat[i]]) ) ^
( (x_opt[treat[i]] - x_min[treat[i]]) / (x_max[treat[i]] - x_opt[treat[i]]) )
}
# Priors
for (j in 1:M) {
x_min[j] ~ dunif(0, 14)
x_opt[j] ~ dunif(x_min[j] + 0, x_min[j] + 25)
x_max[j] ~ dunif(x_opt[j] + 0, 40)
# these rather narrow priors worked during model testing
# x_min[j] ~ dunif(5, 14)
# x_opt[j] ~ dunif(15, 22)
# x_max[j] ~ dunif(23, 26)
c_m[j] ~ dunif(0, 10)
}
sigma ~ dunif(0, 100)
}
'
model_beta_cnx <- textConnection(model_beta_fun_poisson)
# prepare the new two-factor treatment column by combining usinig : notation
# to denote interactions between treat and inf
host_repro_data <- host_repro_data %>% mutate(treat2 = factor(treat):factor(inf))
model_data_repro <- list(y = host_repro_data$sum,
x = host_repro_data$mean_temp,
n = nrow(host_repro_data),
treat = host_repro_data$treat2,
treat_levels = levels(host_repro_data$treat2),
M = length(levels(host_repro_data$treat2)))
# The jags model is initialised using jags.model()
model_repro <- jags.model(model_beta_cnx, data = model_data_repro,
n.chains = 3, quiet = TRUE)
close(model_beta_cnx) # close the model connection as we dont need it any more
# We then use coda.samples to ask for posterior draws
# which we will use as a reflectin of the posterior.
output_repro <- coda.samples(model = model_repro,
variable.names = c("c_m", "x_max", "x_min", "x_opt"),
n.iter = 10000, thin = 5, quiet = TRUE)
gelman.diag(output_repro)
#summarise the model
post_smry_repro <- summary(output_repro)
print(post_smry_repro)
post_means_repro <- post_smry_repro$statistics[,1]
print(post_smry_repro)
#save model output in a text file
#sink("host_stats.csv")
print(summary(output_repro))
#sink()
#save model information in a dataframe
summary_model_repro<- summary(output_repro)
#sub data frame with only quantiles to add to the plot
repro_CI <- as.data.frame(summary_model_repro$quantiles) %>%
tibble::rownames_to_column() %>%
rename(dummy = rowname) %>%
mutate(treatment = ifelse(str_detect(dummy, '1'), paste('CS:I'),ifelse(str_detect(dummy, '2'), paste('CS:U'),ifelse(str_detect(dummy, '3'), paste('FLU:I'),ifelse(str_detect(dummy, '4'), paste('FLU:U'), ifelse(str_detect(dummy, '5'), paste('PULSE:I'),'PULSE:U'))) )))%>%
gather(key = 'Quantile', value = 'value', -dummy, -treatment) %>%
mutate(transform= ifelse(str_detect(dummy, 'c_m'), 10^value, value)) %>% #backtransform maximum spore data
select(-value) %>% #remove raw values to spread data
spread(key = Quantile, value = transform)
# specify a range of temps over which to evaluate our function
# add in the x_min values + a very small amount to ensure the
# lines are drawn at the lowest treshold. The resultant
# vector is sorted sequentially.
# find the indices in post_means_repro that match "x_min" and "x_max"
x_min_idx <- grep("x_min", names(post_means_repro))
x_max_idx <- grep("x_max", names(post_means_repro))
x_opt_idx <- grep("x_opt", names(post_means_repro))
c_m_idx <- grep("c_m" , names(post_means_repro))
# create the vector
xx <- seq(min(post_means_repro[x_min_idx]),
max(post_means_repro[x_max_idx]), length = 10^3)
# and the coresponding estimate based on the means of the posteriors
yy <- matrix(0, nrow = length(xx), ncol = model_data_repro$M)
for (i in 1:model_data_repro$M) {
yy[,i] <- b_fun(xx,
c_m = post_means_repro[c_m_idx[i]],
x_max = post_means_repro[x_max_idx[i]],
x_min = post_means_repro[x_min_idx[i]],
x_opt = post_means_repro[x_opt_idx[i]])
}
#transform the list into a data frame
host_data <- host_repro_data %>%
mutate(id = paste(treat2)) %>% #new column
separate(treat2, c('treat', 'inf'), ':') %>% #separate treat2 column into two new one
group_by(treat, inf, mean_temp) %>% #group our data frame after treatment, inf and mean temp
mutate(mean = mean(sum, na.rm = T), #calculate mean, sd, se for every treat, inf, mean_temp combination
sd = sd(sum, na.rm =T),
se = sd/sqrt(n())) %>%
distinct(treat, inf, mean_temp,id, mean, se ) #remove duplicates, keeps only unique entries for those columns
pred1 <- as.data.frame(cbind(yy, xx) ) #new dataframe containing predictions and temperature (named xx)
names(pred1) <-c('CS:I', 'CS:U', 'FLU:I', 'FLU:U', 'PULSE:I', 'PULSE:U','xx') #name the columns
pred <- pred1 %>% #dataframe consisting of our predictions
pivot_longer(-xx, names_to = "treat2", values_to = "sum") %>% #melt my prediction columns together
separate(treat2, c('treat', 'inf'), ':') %>% #separate column treat2 into two columns named treat and inf
mutate(trans_pred = 10^sum) #create new column named trans_pred with transformed predictions.
#### plots for MS ####
#### Figure 3 ####
label_treat <- c(CS = ' ', PULSE = ' ', FLU = ' ')
data1 <- filter(host_data, treat == 'CS')
pred_1 <- filter(pred, treat == 'CS')
baby1 <- ggplot(data1, aes(y=mean, x=mean_temp, col = treat))+
geom_line(data = pred_1, aes(x = xx, y = trans_pred, col = treat, linetype = inf), size = 0.8, alpha = 0.9)+
geom_point(aes(shape = id),size = 3.5)+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se), width =0.8)+
#t opt
geom_segment(x = repro_CI[19,3] , xend = repro_CI[19,7], y = 132, yend = 132, col = '#003C67FF', linetype = 6, size = 1.3)+ #constant
geom_segment(x = repro_CI[20,3] , xend = repro_CI[20,7], y = 129.5, yend = 129.5, col = '#4A6990FF', size = 1.3)+ #un constant
scale_y_continuous(limits = c(0,120), breaks = c(0,40, 80,120))+
scale_x_continuous(breaks = seq(10, 28,3), limits=c(8,30))+
scale_shape_manual(values=c( 16,1,17,2,15,0))+
scale_color_manual(values=c('#003C67FF'),labels = c('constant'))+
scale_linetype_manual(values=c(6,1), labels = c('exposed', 'unexposed'))+
coord_cartesian(ylim = c(0, 120), clip="off")+ #changes coord system to include the CI arrows in our plot
labs(x = '', y ='reproductive output', shape = ' ', linetype = 'Exposure',color = 'Treatment')+
facet_wrap(~treat, labeller = labeller(treat = label_treat))+
theme( panel.background = element_rect(fill = NA), #removes background
#panel.grid.major.y = element_line(color='grey', linetype = 'dashed', size=0.2),
panel.border= element_rect(colour = "black", fill=NA, size=0.5),
strip.background = element_rect(fill = NA),
legend.background = element_blank(),
legend.title = element_text(hjust=0), #moves text to the left
legend.key = element_blank(),
legend.position = 'none',
text = element_text(size=17),
plot.margin = unit(c(1,0.8,1.5,0), "cm"))
ggpar(baby1)
##ggplot of offspring in fluctuating treatment
data2 <- filter(host_data, treat == 'FLU') #- subset
pred_2 <- filter(pred, treat == 'FLU')
baby2 <- ggplot(data2, aes(y=mean, x=mean_temp, col = treat))+
geom_line(data = pred_2, aes(x = xx, y = trans_pred, col = treat, linetype = inf), size = 0.8, alpha = 0.9)+
geom_point(aes(shape = id),size = 3.5)+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se), width =0.8)+
#
#t opt
geom_segment(x = repro_CI[21,3] , xend = repro_CI[21,7], y = 132, yend = 132, col = '#EFC000FF', linetype = 6, size = 1.3)+ #flux
geom_segment(x = repro_CI[22,3] , xend = repro_CI[22,7], y = 129.5, yend = 129.5, col = '#c7b514', size = 1.3)+ #un flux
# scale adjusts
scale_y_continuous(breaks=NULL)+
scale_x_continuous(breaks = seq(10, 28,3), limits=c(8,30))+
scale_shape_manual(values=c( 17,2))+
scale_color_manual(values=c('#EFC000FF','#A73030FF'),labels = c( 'fluctuation', 'heat wave'))+
scale_linetype_manual(values=c(6,1), labels = c('exposed', 'unexposed'))+
coord_cartesian(ylim = c(0, 120), clip="off")+ #changes coord system to include the CI arrows in our plot
labs(x = 'mean Temperature (in °C)', y =' ', shape = ' ', linetype = 'Exposure',color = 'Treatment')+
facet_wrap(~treat, labeller = labeller(treat = label_treat))+
theme( panel.background = element_rect(fill = NA),
#panel.grid.major.y = element_line(color='grey', linetype = 'dashed', size=0.2),
panel.border= element_rect(colour = "black", fill=NA, size=0.5),
strip.background = element_rect(fill = NA),
legend.background = element_blank(),
legend.title = element_text(hjust=0),
legend.key = element_blank(),
legend.position = 'none',
text = element_text(size=17),
plot.margin = unit(c(1,0.8,1.5,0), "cm"))
ggpar(baby2)
#ggsave(plot = last_plot(), file = 'final_babies_per_treatment.tiff', width = 14, height = 10)
#pulse baby
data3 <- filter(host_data, treat == 'PULSE')
pred_3 <- filter(pred, treat == 'PULSE')
baby3 <- ggplot(data3, aes(y=mean, x=mean_temp, col = treat))+
geom_line(data = pred_3, aes(x = xx, y = trans_pred, col = treat, linetype = inf), size = 0.8, alpha = 0.9)+
geom_point(aes(shape = inf),size = 3.5)+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se), width =0.8)+
#C max
geom_segment(x = 31.7, xend = 31.7, y = repro_CI[1,3], yend = repro_CI[1,7], col = '#003C67FF' , size = 1.3)+#constant # inf
geom_segment(x = 32.5, xend = 32.5, y = repro_CI[3,3],yend =repro_CI[3,7], col ='#EFC000FF', size = 1.3)+#flux #inf
geom_segment(x = 33, xend = 33, y = repro_CI[5,3],yend = repro_CI[5,7],col = '#A73030FF', size = 1.3) +#pulse #inf
geom_segment(x = 31.9 , xend = 31.9, y = repro_CI[2,3], yend = repro_CI[2,7], col ='#4A6990FF', size = 1.3)+ #constant un
geom_segment(x = 32.7, xend = 32.7, y = repro_CI[4,3], yend = repro_CI[4,7], col ='#EFC000FF',size = 1.3)+ #flux un
geom_segment(x = 33 , xend = 33, y = repro_CI[6,3], yend = repro_CI[6,7], col = '#CD534CFF' , size = 1.3)+ #pulse un
#t opt
geom_segment(x = repro_CI[23,3] , xend = repro_CI[23,7], y = 129.5, yend = 129.5, col = '#A73030FF',linetype = 6, size = 1.3)+ #pulse
geom_segment(x = repro_CI[24,3] , xend = repro_CI[24,7], y = 129.5, yend = 129.5, col = '#CD534CFF', size = 1.3)+ #un pulse
# scale adjusts
scale_y_continuous(breaks=NULL)+
scale_x_continuous(breaks = seq(10, 28,3), limits=c(8,30))+
scale_shape_manual(values=c(15,0))+
scale_color_manual(values=c('#A73030FF'),labels = c( 'heat wave'))+
scale_linetype_manual(values=c(6,1), labels = c('exposed', 'unexposed'))+
coord_cartesian(ylim = c(0, 120), clip="off")+ #changes coord system to include the CI arrows in our plot
labs(x = ' ', y =' ', shape = ' ', linetype = 'Exposure',color = 'Treatment')+
facet_wrap(~treat, labeller = labeller(treat = label_treat))+
theme( panel.background = element_rect(fill = NA),
#panel.grid.major.y = element_line(color='grey', linetype = 'dashed', size=0.2),
panel.border= element_rect(colour = "black", fill=NA, size=0.5),
strip.background = element_rect(fill = NA),
legend.background = element_blank(),
legend.title = element_text(hjust=0),
legend.key = element_blank(),
legend.position = 'none',
text = element_text(size=17),
plot.margin = unit(c(1,0.8,1.5,0), "cm"))
ggpar(baby3)
library(cowplot)
allBAB<-plot_grid(baby1,baby2,baby3,nrow = 1)
ggsave(allBAB, file = 'allBABIES.tiff', width = 12, height = 7)
#### Figure 2c ####
#### only infected animals ###
inf_pred <-pred %>%
filter(inf == 'I')
inf_host_data <- host_data%>%
filter(inf == 'I')
plot <- ggplot(inf_host_data, aes(y=mean, x=mean_temp, col = treat, shape = treat))+
geom_line(data = inf_pred, aes(x = xx, y = trans_pred, col = treat), alpha = 0.9,linetype = 'longdash', size = 0.8)+
geom_point(size = 3.5)+ #change shape of my points
geom_errorbar(aes(ymin = mean-se, ymax = mean+se), width = .2)+
#C max
geom_segment(x = 31.8 , xend = 31.8, y = repro_CI[1,3], yend = repro_CI[1,7], col = '#003C67FF', size = 1.3)+ #constant
geom_segment(x = 32.2, xend = 32.2, y = repro_CI[3,3],yend =repro_CI[3,7], col ='#EFC000FF', size = 1.3)+ #flux
geom_segment(x = 32.6 , xend = 32.6, y = repro_CI[5,3],yend = repro_CI[5,7], col = '#A73030FF',size = 1.3)+ #pulse
#t opt
geom_segment(x = repro_CI[21,3] , xend = repro_CI[22,7], y = 135, yend = 135, col = '#003C67FF', size = 1.3)+ #constant
geom_segment(x = repro_CI[22,3] , xend = repro_CI[22,7], y = 133, yend = 133, col = '#EFC000FF', size = 1.3)+ #flux
geom_segment(x = repro_CI[23,3] , xend = repro_CI[23,7], y = 131, yend = 131, col = '#A73030FF',size = 1.3)+ #pulse
scale_x_continuous(breaks = seq(10, 28,3), limits=c(8,30))+
scale_y_continuous(limits = c(0,120), breaks = seq(0,120, 30))+
scale_shape_manual(values=c( 16,17,15), labels = c('constant', 'fluctuating', 'heat wave'))+
scale_color_manual(values=c('#003C67FF','#EFC000FF','#A73030FF'),labels = c('constant', 'fluctuating', 'heat wave'))+
labs(title = ' ',x = 'mean Temperature (in °C)', y = 'reproductive output', col = ' ', shape = ' ')+
coord_cartesian(ylim = c(0, 120), clip="off")+ #changes coord system to include the CI arrows in our plot
theme( panel.background = element_rect(fill = NA),
#panel.grid.major.y = element_line(color='grey', linetype = 'dashed', size=0.2),
panel.border= element_rect(colour = "black", fill=NA, size=0.5),
strip.background = element_rect(color = 'black', fill = 'grey95'),
legend.background = element_blank(),
legend.title = element_text(hjust=0),
legend.position =c(0.85,0.94),
legend.key = element_blank(),
text = element_text(size=21),
plot.margin = unit(c(1,2,1.5,1.5), "cm")) #t, r, b, l bzw. top, right, bottom, left
ggpar(plot)
Computing file changes ...