#contains basic steps from the raw csv file


adu <- read.csv(file="C:/Users/host database.csv",header=TRUE, sep=",")

adu$Growth <- adu$end_length - adu$start_length
adu$inf <- as.factor(adu$inf)
adu$tempi <-adu$temp
adu$temp <- as.factor(adu$temp)
#the duration of the experiment
adu$date <- as.Date(adu$date, format = "%d/%m/%Y")
adu$end_date <- as.Date(adu$end_date, format = "%d/%m/%Y")
adu$days<-as.numeric(difftime(adu$end_date, adu$date, units = "days"))

#survival and cox regression:
library(survival)
adu$inf <- factor(adu$inf, levels = c("0", "1", "2"))

# Create a survival object
surv_obj <- Surv(time = adu$days, event = adu$died)
#the models:
cox_model <- coxph(surv_obj ~ inf + tempi, data = adu)
cox_modeli <- coxph(surv_obj ~ inf * tempi, data = adu)

# Set options to display coefficients without scientific notation
options(scipen = 999)

summary(cox_model)
summary(cox_modeli)
AIC(cox_model)
AIC(cox_modeli)

# Compute residuals for the best model
martingale_res <- residuals(cox_model, type = "martingale")
deviance_res <- residuals(cox_model, type = "deviance")

# Plot Martingale residuals
plot_martingale <- function(residuals) {
  plot(residuals, pch = 16, xlab = "Observation", ylab = "Martingale Residuals", 
       main = "Martingale Residuals")
  abline(h = 0, col = "red", lty = 2)
}
plot_martingale(martingale_res)

# Plot Deviance residuals
plot_deviance <- function(residuals) {
  plot(residuals, pch = 16, xlab = "Observation", ylab = "Deviance Residuals", 
       main = "Deviance Residuals")
  abline(h = 0, col = "red", lty = 2)
}
plot_deviance(deviance_res)

#just 26:
c26<-subset(adu, temp==26)
surv_obj26 <- Surv(time = c26$days, event = c26$died)
# Fit Cox proportional hazards model
cox_model26 <- coxph(surv_obj26 ~ inf , data = c26)
summary(cox_model26)




#removing mussels that died before dissection and mussels that were not representative of their infection category for growth, condition and reproduction models:
ass<- subset(adu, success==1)
as<- subset(ass, died==0)

as$growth<-as$Growth / 10
asg<-subset(as, Growth>=0)
asr<- subset(as, CI2_fixed<10)
asr$CI1<- as.numeric(asr$CI1)

#colours so that control is green and infected mussels blue and purple based on parasite:
inf_colors3 <- c("0" = "green", "1" = "blue", "2" = "purple")

#growth model selection
#was data deficient: g0 <- lmer(growth ~  parasite_number * days *tempi *inf+ (1 | tank), data = as)
g0 <- lmer(growth ~  parasite_number + days *tempi *inf+ (1 | tank), data = as)
summary(g0)
g11<- lm(growth ~  parasite_number * days *tempi *inf, data = as)
summary(g11)
g1<- lm(growth ~  parasite_number + days *tempi *inf, data = as)
summary(g1)
g2<- lm(growth ~  parasite_number + days +tempi *inf, data = as)
summary(g2)
g3<- lm(growth ~  parasite_number + days +tempi +inf, data = as)
summary(g3)
g4<- lm(growth ~  days +tempi +parasite_number *inf, data = as)
summary(g4)
g51<- lm(growth ~  days *tempi *inf, data = as)
summary(g51)
g52<- lm(growth ~  days *tempi +inf, data = as)
summary(g52)
g53<- lm(growth ~  days +tempi *inf, data = as)
summary(g53)
g5<- lm(growth ~  days +tempi +inf, data = as)
summary(g5)
g6 <- lmer(growth ~  days +tempi +inf+ (1 | tank), data = as)
summary(g6)
anova(g1, g2, g3, g4, g5)

AIC(g0)
AIC(g11)
AIC(g1)
AIC(g2)
AIC(g3)
AIC(g4)
AIC(g51)
AIC(g52)
AIC(g53)
AIC(g5)
AIC(g6)

# Checking the fit of the best model
residuals <- residuals(g5)

# Create diagnostic plots
par(mfrow=c(2,2))
plot(g5)

#three influential points detected. First trying out how much they affect the results: 
# Create a vector of numbers to exclude
exclude_numbers <- c("E48", "E114", "E122")

# Subset the data to exclude the influential points
asg_subset <- subset(asg, !(number %in% exclude_numbers))

g5wo_outliers<- lm(growth ~  days +tempi +inf, data = asg_subset)
summary(g5wo_outliers)

#reproduction model selection:
# Fit a negative binomial mixed-effects model

# Sample mean and variance excluding NAs
mean_val <- mean(as$repro, na.rm = TRUE)
var_val <- var(as$repro, na.rm = TRUE)

mean_val
var_val

# Method of Moments estimation of theta
theta_mom <- mean_val^2 / (var_val - mean_val)

#first with the random effect:
g1 <- glmer(repro ~  days *tempi *inf+ (1 | tank), family = negative.binomial(theta = 3.0), data = as)
summary(g1)
g2 <- glmer(repro ~  days +tempi *inf+ (1 | tank), family = negative.binomial(theta = 3.0), data = as)
summary(g2)
g2 <- glmer.nb(repro ~ days + tempi * inf + (1 | tank), data = as)
g3 <- glmer(repro ~  days +tempi +inf+ (1 | tank), family = negative.binomial(theta = 3.0), data = as)
summary(g3)
#all with tank result in singularity issues!

# Fit the model without any random effects using glm
r0<- glm(repro ~ days + parasite_number + tempi * inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r01<- glm(repro ~ days * parasite_number * tempi * inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r02<- glm(repro ~ days  + tempi * inf * parasite_number, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r03<- glm(repro ~ days + parasite_number + tempi + inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r1<- glm(repro ~ days * tempi * inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r12<- glm(repro ~ days + tempi * inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
r13<- glm(repro ~ days * tempi + inf, 
          family = negative.binomial(theta = 3.0), 
          data = as)
r2<- glm(repro ~ days + tempi + inf, 
         family = negative.binomial(theta = 3.0), 
         data = as)
summary(r0)
summary(r01)
summary(r02)
summary(r03)
summary(r1)
summary(r12)
summary(r13)
summary(r2)

#condition:
g1<- lm(CI2_fixed ~  parasite_number + days *tempi *inf, data = asr)
summary(g1)
g11<- lm(CI2_fixed ~  parasite_number * days *tempi *inf, data = asr)
summary(g11)
g2<- lm(CI2_fixed ~  parasite_number + days +tempi *inf, data = asr)
summary(g2)
g3<- lm(CI2_fixed ~  parasite_number + days +tempi +inf, data = asr)
summary(g3)
g4<- lm(CI2_fixed~  days +tempi +parasite_number *inf, data = asr)
summary(g4)
g5<- lm(CI2_fixed ~  days +tempi +inf, data = asr)
summary(g5)
g6 <- lmer(CI2_fixed ~  days +tempi +inf+ (1 | tank), data = asr)
summary(g6)
g7<- lm(CI2_fixed ~  days +tempi *inf, data = asr)
summary(g7)
g71<- lm(CI2_fixed ~  days *tempi +inf, data = asr)
summary(g71)
g8 <- lmer(CI2_fixed ~  days +tempi *inf+ (1 | tank), data = asr)
summary(g8)
anova(g1, g2, g3, g4, g5)

AIC(g1)
AIC(g11)
AIC(g2)
AIC(g3)
AIC(g4)
AIC(g5)
AIC(g6)
AIC(g7)
AIC(g71)
AIC(g8)

# Checking the fit of the best model
residuals <- residuals(g7)

# Create diagnostic plots
par(mfrow=c(2,2))
plot(g7)



#Figures

#kaplan-meier:
# Create a survival object
surv_obj <- Surv(time = adu$days, event = adu$died)

# Fit Kaplan-Meier survival curves by 'inf'
km_fit <- survfit(surv_obj ~ inf, data = adu)

# Plot Kaplan-Meier survival curves
plot(km_fit, col = c("blue", "red", "green"), lty = c(1, 2, 3),
     xlab = "Time (Days)", ylab = "Survival Probability",
     main = "Kaplan-Meier Survival for all temperatures")

# Fit Kaplan-Meier survival curves by 'temp'
km_fit <- survfit(surv_obj ~ temp, data = adu)

# Plot Kaplan-Meier survival curves
plot(km_fit, col = c("blue", "red", "green", "purple", "grey"), lty = c(1, 2, 3, 4, 5),
     xlab = "Time (Days)", ylab = "Survival Probability",
     main = "Kaplan-Meier Survival for all temperatures")

# Fit Kaplan-Meier survival curves by 'temp' and 'inf'
km_fit <- survfit(surv_obj ~ temp+inf, data = adu)

# Plot Kaplan-Meier survival curves
plot(km_fit, col = c("blue", "red", "green", "purple", "grey"), lty = c(1, 2, 3, 4, 5),
     xlab = "Time (Days)", ylab = "Survival Probability",
     main = "Kaplan-Meier Survival for all temperatures")


# Subset the data
c26 <- subset(adu, temp == 26)
c18 <- subset(adu, temp == 18)
c14 <- subset(adu, temp == 14)

# Function to create the plots
create_plots <- function() {
  # Set up a 1x3 layout for the plots with adjusted margins
  par(mfrow = c(1, 3), mar = c(4, 4, 2, 1) + 0.1)  

  # Custom parameters for thicker lines and increased font size
  par(cex = 1.2, cex.lab = 1.2, cex.axis = 1.0, cex.main = 1.5, lwd = 2)

  # Plot for 14 C
  surv_obj14 <- Surv(time = c14$days, event = c14$died)
  km_fit14 <- survfit(surv_obj14 ~ inf, data = c14)
  plot(km_fit14, col = c("green", "blue", "purple"), 
       xlab = "", ylab = "Survival Probability",
       main = "14 C", lwd = 5, xlim = c(0, 140))  
  grid()

  # Plot for 18 C
  surv_obj18 <- Surv(time = c18$days, event = c18$died)
  km_fit18 <- survfit(surv_obj18 ~ inf, data = c18)
  plot(km_fit18, col = c("green", "blue", "purple"), 
       xlab = "Time (days)", ylab = "",  
       main = "18 C", lwd = 5, xlim = c(0, 140), yaxt = "n")  
  grid()

  # Add legend with increased font size
  legend("bottom", legend = c(expression("Control"), expression(italic("M. orientalis")), expression(italic("M. intestinalis"))),
         col = c("green", "blue", "purple"), title = "Exposure", cex = 1.2, lwd = 5,
         text.font = 2)  

  # Plot for 26 C
  surv_obj26 <- Surv(time = c26$days, event = c26$died)
  km_fit26 <- survfit(surv_obj26 ~ inf, data = c26)
  plot(km_fit26, col = c("green", "blue", "purple"), 
       xlab = "", ylab = "",  
       main = "26 C", lwd = 5, xlim = c(0, 140), yaxt = "n")  
  grid()
}

# Save as EPS
setEPS()
postscript("output.eps", width = 14, height = 7)
create_plots()
dev.off()

# Save as TIFF
tiff("output.tiff", width = 14, height =7, units = "in", res = 300)
create_plots()
dev.off()


#growth figure

#mussel growth over time without excluding negative values:
t10<- subset(as, temp==10)
t14<- subset(as, temp==14)
t18<- subset(as, temp==18)
t22<- subset(as, temp==22)
t26<- subset(as, temp==26)

#colours so that control is green and infected mussels blue and purple based on parasite:
inf_colors3 <- c("0" = "green", "1" = "blue", "2" = "purple")

new_legend_labels <- c(expression("Control"), expression(italic("M. orientalis")), expression(italic("M. intestinalis")))

p10g <- ggplot(t10, aes(x = days, y = growth, color = inf)) +
  geom_point(size = 3) +
  labs(title = "     10 C",
       x = "Days",
       y = "") +
  scale_y_continuous(breaks = seq(0, 2, by = 1),
                     limits = c(-0.5, 2)) +
  scale_color_manual(values = inf_colors3, name = "", labels = new_legend_labels) +  # Add labels here
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

# Ensure the legend is positioned at the bottom
p10g <- p10g + theme(legend.position = "bottom")

p14g <- ggplot(t14, aes(x = days, y = growth, color = inf)) +
  geom_point(size = 3) +
  labs(title = "     14 C",
       x = "",
       y = "") +
    scale_y_continuous(breaks = seq(0, 2, by = 1),
                     limits = c(-0.5, 2)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p18g <- ggplot(t18, aes(x = days, y = growth, color = inf)) +
  geom_point(size = 3) +
  labs(title = "     18 C",
       x = "",
       y = "Growth (mm)") +
  scale_y_continuous(breaks = seq(0, 2, by = 1),
                     limits = c(-0.5, 2)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p22g <- ggplot(t22, aes(x = days, y = growth, color = inf)) +
  geom_point(size = 3) +
  labs(title = "     22 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(0, 2, by = 1),
                     limits = c(-0.5, 2)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p26g <- ggplot(t26, aes(x = days, y = growth, color = inf)) +
  geom_point(size = 3) +
  labs(title = "b.   26 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(0, 2, by = 1),
                     limits = c(-0.5, 2)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

#version where the left side of Figure 3 has growth over time:
as$gtime<-as$growth / as$days

#growth flipped:
growthall2 <- ggplot(as, aes(x = temp, y = gtime, fill = factor(inf))) +
  geom_boxplot() +
  coord_flip() +  # To flip the x and y axes
scale_y_continuous(breaks = seq(-0.01, 0.025, by = 0.01),
                   limits = c(-0.01, 0.025)) +
  ggtitle("a.") +
  xlab("Temperature ( C)") +  # Add x-axis label
  ylab("Mussel growth (mm/day)") +   # Add y-axis label
  scale_fill_manual(values = inf_colors3)

# Increase font size 
custom_theme <- theme(
  axis.title = element_text(size = 14),     
  axis.text = element_text(size = 12),       
  legend.text = element_text(size = 12),     
  legend.title = element_text(size = 12),    
  plot.title = element_text(size = 12),      
  strip.text = element_text(size = 14)       
)

# Apply the custom theme to each plot
growthall2 <- growthall2 + theme(legend.position = "none") + custom_theme
p10g <- p10g + theme(legend.position = "bottom") + custom_theme
p14g <- p14g + theme(legend.position = "none") + custom_theme
p18g <- p18g + theme(legend.position = "none") + custom_theme
p22g <- p22g + theme(legend.position = "none") + custom_theme
p26g <- p26g + theme(legend.position = "none") + custom_theme

# Define the layout for reproall and other plots
growthall_layout <- ggarrange(growthall2, ncol = 1)

other_plots <- ggarrange(p26g, p22g, p18g, p14g, p10g, 
                         ncol = 1, nrow = 5, heights = c(1, 1, 1, 1, 1.35))

# Combine reproall and other plots
grouped_plot <- ggarrange(growthall_layout, other_plots,
                          widths = c(1, 1), 
                          common.legend = TRUE, legend = "bottom")

# Display the grouped plot
print(grouped_plot)

ggsave("growth0706final.tiff", grouped_plot, width = 11, height = 8, units = "in", dpi = 300)


#condition index figure :
#mussel CI2 over time:
t10<- subset(asr, temp==10)
t14<- subset(asr, temp==14)
t18<- subset(asr, temp==18)
t22<- subset(asr, temp==22)
t26<- subset(asr, temp==26)

#colours so that control is green and infected mussels blue and purple based on parasite:
inf_colors3 <- c("0" = "green", "1" = "blue", "2" = "purple")
#condition flipped:
CI2all <- ggplot(asr, aes(x = temp, y = CI2_fixed, fill = factor(inf))) +
  geom_boxplot() +
  coord_flip() +  # To flip the x and y axes
scale_y_continuous(breaks = seq(1, 8.5, by = 1),
                     limits = c(0.5, 8.5)) +
  ggtitle("a.") +
  xlab("Temperature ( C)") +  # Add x-axis label
  ylab("Mussel condition") +   # Add y-axis label
  scale_fill_manual(values = inf_colors3)

CI2all

new_legend_labels <- c(expression("Control"), expression(italic("M. orientalis")), expression(italic("M. intestinalis")))

p10g <- ggplot(t10, aes(x = days, y = CI2_fixed, color = inf)) +
  geom_point(size = 3) +
 # geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     10 C",
       x = "Days",
       y = "") +
  scale_y_continuous(breaks = seq(1, 8.5, by = 3),
                     limits = c(0.5, 8.5)) +
  scale_color_manual(values = inf_colors3, name = "", labels = new_legend_labels) +  # Add labels here
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

# Ensure the legend is positioned at the bottom
p10g <- p10g + theme(legend.position = "bottom")

p14g <- ggplot(t14, aes(x = days, y = CI2_fixed, color = inf)) +
  geom_point(size = 3) +
 # geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     14 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(1, 8.5, by = 3),
                     limits = c(0.5, 8.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p18g <- ggplot(t18, aes(x = days, y = CI2_fixed, color = inf)) +
  geom_point(size = 3) +
  #geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     18 C",
       x = "",
       y = " Condition ") +
  scale_y_continuous(breaks = seq(1, 8.5, by = 3),
                     limits = c(0.5, 8.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p22g <- ggplot(t22, aes(x = days, y = CI2_fixed, color = inf)) +
  geom_point(size = 3) +
  #geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     22 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(1, 8.5, by = 3),
                     limits = c(0.5, 8.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

p26g <- ggplot(t26, aes(x = days, y = CI2_fixed, color = inf)) +
  geom_point(size = 4) +
 # geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "b.   26 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(1, 8.5, by = 3),
                     limits = c(0.5, 8.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
scale_x_continuous(limits = c(55, 145))

# Increase font size 
custom_theme <- theme(
  axis.title = element_text(size = 14),     
  axis.text = element_text(size = 12),       
  legend.text = element_text(size = 12),     
  legend.title = element_text(size = 12),    
  plot.title = element_text(size = 12),      
  strip.text = element_text(size = 14)       
)

# Apply the custom theme to each plot
CI2all <- CI2all + theme(legend.position = "none") + custom_theme
p10g <- p10g + theme(legend.position = "bottom") + custom_theme
p14g <- p14g + theme(legend.position = "none") + custom_theme
p18g <- p18g + theme(legend.position = "none") + custom_theme
p22g <- p22g + theme(legend.position = "none") + custom_theme
p26g <- p26g + theme(legend.position = "none") + custom_theme

# Define the layout for reproall and other plots
CI2all_layout <- ggarrange(CI2all, ncol = 1)

other_plots <- ggarrange(p26g, p22g, p18g, p14g, p10g, 
                         ncol = 1, nrow = 5, heights = c(1, 1, 1, 1, 1.35))

# Combine reproall and other plots
grouped_plot <- ggarrange(CI2all_layout, other_plots,
                          widths = c(1, 1), 
                          common.legend = TRUE, legend = "bottom")

# Display the grouped plot
print(grouped_plot)
ggsave("condition0706final.tiff", grouped_plot, width = 11, height = 8, units = "in", dpi = 300)

#reproduction figure:
t10<- subset(as, temp==10)
t14<- subset(as, temp==14)
t18<- subset(as, temp==18)
t22<- subset(as, temp==22)
t26<- subset(as, temp==26)

# Repro flipped:
reproall <- ggplot(as, aes(x = temp, y = repro, fill = factor(inf))) +
  geom_boxplot() +
  coord_flip() +  # To flip the x and y axes
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(0, 3)) +
  ggtitle("a.") +
  xlab("Temperature ( C)") +  # Add x-axis label
  ylab("Reproduction index") +   # Add y-axis label
  scale_fill_manual(values = inf_colors3)

reproall

#alternative as a dot plot:

library(dplyr)

# count observations
as_summary <- as %>%
  group_by(temp, repro, inf) %>%
  summarise(count = n()) %>%
  ungroup()

reproall <- ggplot(as_summary, aes(x = temp, y = repro, size = count, fill = factor(inf))) +
  geom_point(shape = 21, color = "black", alpha = 0.7, 
             position = position_dodge(width = 0.5)) +  
  coord_flip() +  # To flip the x and y axes
  scale_y_continuous(breaks = seq(0, 3, by = 1), limits = c(0, 3)) +
  ggtitle("a.") +
  xlab("Temperature ( C)") +  
  ylab("Reproduction index") +   
  scale_fill_manual(values = inf_colors3) +  
  scale_size_continuous(range = c(2, 10)) +  
  theme_minimal()  

# Display the plot
print(reproall)

new_legend_labels <- c(expression("Control"), expression(italic("M. orientalis")), expression(italic("M. intestinalis")))

# Create p10g plot with new legend labels
p10g <- ggplot(t10, aes(x = days, y = repro, color = inf)) +
  geom_point(position = position_jitter(width = 2, height = 0.15), size = 3) +
  labs(title = "     10 C",
       x = "Days",
       y = "") +
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(-0.5, 3.5)) +
  scale_color_manual(values = inf_colors3, name = "", labels = new_legend_labels) +  # Add labels here
  theme_minimal() +
  theme(
    plot.title = element_text(size = 10)
  ) + 
  scale_x_continuous(limits = c(55, 145))

# Ensure the legend is positioned at the bottom
p10g <- p10g + theme(legend.position = "bottom")

p14g <- ggplot(t14, aes(x = days, y = repro, color = inf)) +
  geom_point(position = position_jitter(width = 2, height = 0.15), size = 3) +  # Adding jitter  #geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     14 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(-0.5, 3.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
  scale_x_continuous(limits = c(55, 145))

p18g <- ggplot(t18, aes(x = days, y = repro, color = inf)) +
  geom_point(position = position_jitter(width = 2, height = 0.15), size = 3) +  # Adding jitter # geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     18 C",
       x = "",
       y = "Reproduction") +
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(-0.5, 3.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
  scale_x_continuous(limits = c(55, 145))

p22g <- ggplot(t22, aes(x = days, y = repro, color = inf)) +
  geom_point(position = position_jitter(width = 2, height = 0.15), size = 3) +  # Adding jitter  #geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "     22 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(-0.5, 3.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
  scale_x_continuous(limits = c(55, 145))

p26g <- ggplot(t26, aes(x = days, y = repro, color = inf)) +
  geom_point(position = position_jitter(width = 2, height = 0.15), size = 3) +  # Adding jitter
  # geom_smooth(method = "lm", se = FALSE, aes(group = inf)) +
  labs(title = "b.   26 C",
       x = "",
       y = "") +
  scale_y_continuous(breaks = seq(0, 3, by = 1),
                     limits = c(-0.5, 3.5)) +
  scale_color_manual(values = inf_colors3) +
  theme_minimal()+
  theme(
    plot.title = element_text(size = 10)   ) + 
  scale_x_continuous(limits = c(55, 145))

# Increase font size 
custom_theme <- theme(
  axis.title = element_text(size = 14),     
  axis.text = element_text(size = 12),       
  legend.text = element_text(size = 12),     
  legend.title = element_text(size = 12),    
  plot.title = element_text(size = 12),      
  strip.text = element_text(size = 14)       
)

reproall <- reproall + theme(legend.position = "none")+ custom_theme
p10g <- p10g + theme(legend.position = "bottom")+ custom_theme 
p14g <- p14g + theme(legend.position = "none")+ custom_theme
p18g <- p18g + theme(legend.position = "none")+ custom_theme
p22g <- p22g + theme(legend.position = "none")+ custom_theme
p26g <- p26g + theme(legend.position = "none")+ custom_theme

# Define the layout for reproall and other plots
reproall_layout <- ggarrange(reproall, ncol = 1)

other_plots <- ggarrange(p26g, p22g, p18g, p14g,  p10g, 
                         ncol = 1, nrow = 5, heights = c(1, 1, 1, 1, 1.4))

# Combine reproall and other plots
grouped_plot <- ggarrange(reproall_layout, other_plots,
                          widths = c(1, 1), 
                          common.legend = TRUE, legend = "bottom")

# Display the grouped plot
print(grouped_plot)
ggsave("repro706final.tiff", grouped_plot, width = 11, height = 8, units = "in", dpi = 300)

#supplementary figure 1 parasite number
inf_colors <- c("2" = "blue", "3" = "purple")
as$inf<-as.numeric(as$inf)
as_filtered<-subset(as, inf>1)

inf_colors3 <- c("0" = "green", "1" = "blue", "2" = "purple")

as_filtered$inf<-as.factor(as_filtered$inf)
inf_colors <- c("2" = "blue", "3" = "purple")


# Create the boxplot with counts
palladults<-ggplot(as_filtered, aes(x = temp, y = as.numeric(parasite_number), fill = inf)) + 
  geom_boxplot() +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  ggtitle("Parasite load by temperature")+
  xlab("Temperature ( C)") +  # Add x-axis label
  ylab("Parasite number")     +   # Add y-axis label
  scale_fill_manual(values = inf_colors) 

palladults
#separate plots for different temps over time:
t10<- subset(as_filtered, temp==10)
t14<- subset(as_filtered, temp==14)
t18<- subset(as_filtered, temp==18)
t22<- subset(as_filtered, temp==22)
t26<- subset(as_filtered, temp==26)

#juveniles included:
p10<- ggplot(t10, aes(x = days, y = parasite_number, color = inf)) +
  geom_point(size = 5) +  # Set the size to your preferred value
  labs(title = "Parasite load at 10 C",
       x = "Days",
       y = "Parasite number") +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  scale_color_manual(values = inf_colors) +
  theme_minimal()
p14<- ggplot(t14, aes(x = days, y = parasite_number, color = inf)) +
  geom_point(size = 5) +  # Set the size to your preferred value
  labs(title = "Parasite load at 14 C",
       x = "Days",
       y = "Parasite number") +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  scale_color_manual(values = inf_colors) +
  theme_minimal()
p18<- ggplot(t18, aes(x = days, y = parasite_number, color = inf)) +
  geom_point(size = 5) +  # Set the size to your preferred value
  labs(title = "Parasite load at 18 C",
       x = "Days",
       y = "Parasite number") +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  scale_color_manual(values = inf_colors) +
  theme_minimal()
p22<- ggplot(t22, aes(x = days, y = parasite_number, color = inf)) +
  geom_point(size = 5) +  # Set the size to your preferred value
  labs(title = "Parasite load at 22 C",
       x = "Days",
       y = "Parasite number") +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  scale_color_manual(values = inf_colors) +
  theme_minimal()
p26<- ggplot(t26, aes(x = days, y = parasite_number, color = inf)) +
  geom_point(size = 5) +  # Set the size to your preferred value
  labs(title = "Parasite load at 26 C",
       x = "Days",
       y = "Parasite number") +
  scale_y_continuous(breaks = seq(0, 25, by = 5),
                     limits = c(0, 25)) +
  scale_color_manual(values = inf_colors) +
  theme_minimal()

# Set the new legend labels
#new_legend_labels <- c("Mytilicola orientalis", "Mytilicola intestinalis")
new_legend_labels <- c(expression(italic("M. orientalis")), expression(italic("M. intestinalis")))

palladults <- palladults + scale_fill_manual(values = inf_colors, name = "Infection", labels = new_legend_labels)
ggarrange(palladults, p10, p14, p18, p22, p26,
         ncol = 2, nrow = 3, common.legend = TRUE, legend = "bottom")


#Means and 95% confidence intervals for growth:

as$gtime<-as$growth / as$days

gor<-subset(as, inf==2)
gint<-subset(as, inf==3)
gcon<-subset(as, inf==1)

summary(gor$growth)
summary(gcon$growth)
summary(gint$growth)

summary(gor$gtime)
summary(gcon$gtime)
summary(gint$gtime)

gtime <- as$gtime

#all:
mean_gtime <- mean(as$gtime)
ci <- t.test(as$gtime)$conf.int
print(paste("Mean:", mean_gtime))
print(paste("95% Confidence Interval:", ci))

#controls:
mean_gtime <- mean(gcon$gtime)
ci <- t.test(gcon$gtime)$conf.int
print(paste("Mean:", mean_gtime))
print(paste("95% Confidence Interval:", ci))

# ori:
mean_gtime <- mean(gor$gtime)
ci <- t.test(gor$gtime)$conf.int
print(paste("Mean:", mean_gtime))
print(paste("95% Confidence Interval:", ci))

#int:
mean_gtime <- mean(gint$gtime)
ci <- t.test(gint$gtime)$conf.int
print(paste("Mean:", mean_gtime))
print(paste("95% Confidence Interval:", ci))



#checking for possible bias in the condition values by mussel size:
plot(asr$end_length, asr$CI2_fixed)
cor_test_result <- cor.test(asr$end_length, asr$CI2_fixed)
print(cor_test_result)