#Importing Data
Raw_Data <- read.table("DaMN_Datafile.csv", header = T, sep = ',')

#Importing Libraries 
library(ggplot2) 
library(car)
library(ggpubr)

###############################As factor########################################

# Change the structure of the data from "chr" to factors.
# Change the names too.
Raw_Data$NP_treatment <- factor(Raw_Data$NP_treatment, 
                            levels = c("ZERO", "LOW", "HIGH"),
                            labels = c("Zero", "Low", "High"))

#### Line 18-20 changes the labels to numeric values####
#Data$NP_treatment <- factor(Data$NP_treatment, 
#                            levels = c("ZERO", "LOW", "HIGH"),
#                            labels = c("0", "5", "20"))

Raw_Data$Treatment <- factor(Raw_Data$Treatment,
                         levels = c("A_ZERO_METS" , "A_ZERO____" , "A_LOW_METS" ,
                                    "A_LOW____" , "A_HIGH_METS" , "A_HIGH____" ),
                         labels = c("Zero_Mets", "Zero", "Low_Mets", 
                                    "Low", "High_Mets", "High"))

# Infection treatment = Inf_treatment. 
#Individuals are divided according to their infection status: Control vs Infected
Raw_Data$Inf_treatment <- factor(Raw_Data$Inf_treatment,
                             levels = c("NO_METS","METS"),
                             labels = c("Control", "Infected"))



Raw_Data$ID.No. <- as.factor(Raw_Data$ID.No.)
Raw_Data$Rep. <- as.factor(Raw_Data$Rep.)
Raw_Data$Inf_treatment <- as.factor(Raw_Data$Inf_treatment)
Raw_Data$H_errors <- as.factor(Raw_Data$H_errors)

#################################Data Manipulation##############################

#Remove Handling errors
Data_2 <- subset(Raw_Data, H_errors == 0 )

#Remove Background mortality. Main data!!!
Data_1 <- subset(Data_2, Age_death > 3) 

#Subsets of the main Data

# Only the exposed to the parasite treatment.
Mets <- subset(Data_1, Inf_treatment == "Infected")

# Only the parasite-free treatments.
Nomets <- subset(Data_1, Inf_treatment == "Control") 

#The ones that died before day 8 post exposure
Died_bf_day8 <- subset(Mets, Retrieved == 0) 

# Only the exposed and infected from the parasite individuals.
InfH <- subset(Mets, Spore_yield > 0)

# Combine the the three datasets. 
#Individuals that died after day-8 and not infected were excluded + Control. 
d <- rbind(Died_bf_day8, InfH, Nomets)

#Only the ones that survived day-8 and were successfully infected + Control
Successfully_inf <- rbind(InfH, Nomets)

#From the infection treatment only the ones that survived after day 8.
Inf_after_day8 <- subset(Mets, Retrieved == 1) 

str(Data_1)

################################1.Host Lifespan################################

errorbar_HL <- stat_summary(data = d,
                         fun.data = "mean_se",
                         geom = "errorbar",
                         width = 0.2,
                         position = position_dodge(width = 0.6),
                         colour = "black", alpha = 0.5)

lines_HL <- stat_summary(data = d, fun.data = "mean_se",
                      geom = "line",
                      position = position_dodge(width=0.6),
                      colour = "black",
                      alpha = 0.3)

mean_HL <- stat_summary(data = d, fun.data = "mean_se",
                     geom = "point",
                     position = position_dodge(width=0.6),
                     pch = 21,
                     size = 1.5)

#PLot - Host Lifespan
A1 <- ggplot(d, aes(NP_treatment, Age_death, fill = NP_treatment)) +
  facet_wrap(~ Inf_treatment) +
  geom_jitter(width = 0.2, alpha = 0.25, pch = 20, size = 0.7) +
  scale_y_continuous(breaks=seq(0, 30, 5)) +
  scale_fill_grey(start = 1, end = 0) +
  lines_HL + 
  errorbar_HL +
  mean_HL + 
  labs( x = "\nNanoplastics concentration",
        y = "\nAge at death (days)",
        title = "\nHost lifespan" ) + 
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "none")
A1

# Save plot.
ggsave("plot_DAMN_Host_Lifespan.tiff", 
       units = "in", 
       width = 4.5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')

model_HL <- aov(Age_death ~ NP_treatment * Inf_treatment, data = d)

Anova(model_HL, type = 2)

##  % Variance explained
412.2 + 3769.7 + 59.6 + 3052.2
(412.2 / 7293.7)*100 #NPs
(3769.7 / 7293.7)*100 #Inf
(59.6 / 7293.7)*100 #Interaction
(3052.2 / 7293.7)*100 #Res

TukeyHSD(model_HL)

hist(resid(model_HL))
plot(model_HL)



##############################2.Host Fecundity##################################

#Only the ones that reproduced at least once
Reproduced_at_least1 <- subset(d, Total_juv > 0) 

errorbar_HF <- stat_summary(data = Reproduced_at_least1,
                           fun.data = "mean_se",
                           geom = "errorbar",
                           width = 0.2,
                           position = position_dodge(width = 0.6),
                           colour = "black", alpha = 0.5)

lines_HF <- stat_summary(data = Reproduced_at_least1, fun.data = "mean_se",
                        geom = "line",
                        position = position_dodge(width=0.6),
                        colour = "black",
                        alpha = 0.3)

mean_HF <- stat_summary(data = Reproduced_at_least1, fun.data = "mean_se",
                       geom = "point",
                       position = position_dodge(width=0.6),
                       pch = 21,
                       size = 1.5)

C1 <- ggplot(Reproduced_at_least1, aes(NP_treatment, Total_juv, fill = NP_treatment)) +
  facet_wrap(~ Inf_treatment) +
  geom_jitter(width = 0.1, alpha = 0.25, pch = 20, size = 0.7) +
  lines_HF +
  errorbar_HF +
  mean_HF + 
  coord_cartesian(ylim=c(0, 18)) +
  scale_fill_grey(start = 1, end = 0) +
  scale_y_continuous(breaks=seq(0, 20, 4)) +
  labs( x = "\nNanoplastics concentration", 
        y = "\nTotal number of offspring (per host)", 
        title = "\nHost fecundity") +
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "bottom")

C1 # Host fecundity

ggsave("plot_Host_fecundity.tiff", 
       units = "in", 
       width = 4.5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')


model_HF <- aov(Total_juv ~ NP_treatment * Inf_treatment, data = Reproduced_at_least1)

summary(model_HF)

Anova(model_HF, type = 2)
#  % Variance explained
104.24 + 1343.83 + 71.03 + 683.16
(104.24 / 2202.26)*100 #NP
(1343.83 / 2202.26)*100 #Inf
(71.03 / 2202.26)*100 #Interaction
(683.16 / 2202.26)*100 #Res

TukeyHSD(model_HF)

hist(resid(model_HF))
plot(model_HF)

####################3.Fecundity until the 3rd clutch############################
#
#ggplot(Reproduced_at_least1,aes(NP_treatment, Juv_3_clutch,
#             fill = NP_treatment, group = NP_treatment)) +
#  facet_wrap(~ Inf_treatment) +
#  coord_cartesian(ylim=c(0, 10)) +
#  scale_y_continuous(breaks=seq(0, 10, 2)) +
#  geom_jitter(width = 0.1, alpha = 0.25, pch = 20, size = 0.7) +
#  scale_fill_grey(start = 1, end = 0) +
#  lines_1 +
#  errorbar_1 +
#  mean_1 +
#  labs(x = "\nNanoplastics concentration (mg/L)", 
#       y = "\nTotal number of offspring (per host)",
#       title = "\nFecundity until the 3rd clutch") +
#  theme_bw() +
#  theme(panel.grid.major = element_line(),
#        panel.grid.minor = element_blank(),
#        plot.title = element_text(hjust = 0.5, size = 15),
#        legend.title = element_blank(), 
#        legend.position = "none")
#
#ggsave("plot_Fecundity_at_first_three_clutches.tiff", 
#       units = "in", 
#       width = 4.5, 
#       height = 4, 
#       dpi = 600, 
#       compression = 'lzw')
#
#
#model_F3C <- aov(Juv_3_clutch ~ NP_treatment * Inf_treatment, data = Reproduced_at_least1)
#
#Anova(model_F3C, type = 2)
#
#TukeyHSD(model_F3C)
#plot(model_F3C)
#hist(f$Juv_3_clutch)
#hist(resid(model_F3C))
#
######################4.Host that Reached Maturity##############################


# Create new column; 0 | 1 did not reached maturity | reached maturity.
# maturity = The individuals that reached maturity
d$maturity <-  (d$Age_at_Mat > 0)
d$maturity <- as.numeric(d$maturity,
                          levels = c("TRUE","FALSE"),
                          labels = c("1", "0"))


str(d)

#HRM stands for Host Reached Maturity#

binom_HRM <- subset(d,  maturity == "0" | maturity == "1")


bn_HRM <- glm(formula = binom_HRM$maturity ~ Inf_treatment * NP_treatment, 
             family = binomial(link = "logit"), 
             data = binom_HRM )
summary(bn_HRM)


Anova(bn_HRM, type="II")

mean_HRM <- stat_summary(data = binom_HRM, fun.data = "mean_se",
                        geom = "bar", 
                        size = 0.6,
                        width = 0.2,
                        colour = "black",
                        position = "dodge")

errorbar_HRM <- stat_summary(data = binom_HRM,
                            fun.data = "mean_se",
                            geom = "errorbar",
                            width = 0.05,
                            colour = "darkgrey",
                            position = position_dodge(width = 0.2))

# Plot
B1 <- ggplot(binom_HRM, aes( Inf_treatment, maturity, fill = NP_treatment)) + 
  coord_cartesian(ylim=c(0, 1)) +
  scale_y_continuous(breaks=seq(0, 1, 0.25)) +
  scale_fill_grey(start = 1, end = 0) +
  mean_HRM +
  errorbar_HRM +
  labs( x = "\nNanoplastics concentration",
        y = "\nProportion of Daphnia ",
        title = "\nHost that reached maturity" ) +
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(),
        legend.position = "none")
B1
# Save plot
ggsave("plot_Host_R_mat.tiff", 
       units = "in", 
       width = 6, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')

####Same plot but axis values are flipped. ####


#B1_1 <- ggplot(binom_HRM, aes( NP_treatment, maturity, color = maturity, fill = NP_treatment)) + 
#  facet_wrap(~Inf_treatment) +
#  coord_cartesian(ylim=c(0, 1)) +
#  scale_y_continuous(breaks=seq(0, 1, 0.25)) +
#  scale_fill_grey(start = 1, end = 0) +
#  mean_11 +
#  errorbar_11 +
#  labs( x = "\nNanoplastics concentration",
#        y = "\nHost reaching maturity",
#        title = "\nHost that reached maturity" ) +
#  theme_bw() +
#  theme(panel.grid.major = element_line(),
#        panel.grid.minor = element_blank(),
#        plot.title = element_text(hjust = 0.5, size = 15),
#        legend.title = element_blank(), 
#        legend.position = "none")
#
# 
# Save plot
#ggsave("plot_Host_R_mat_1_1.tiff", 
#       units = "in", 
#       width = 6, 
#       height = 4, 
#       dpi = 600, 
#       compression = 'lzw')

### Combine graphs for Host life traits ####
figure_2  <- ggarrange(A1, B1, C1,
                       labels = c("A", "B", "C"),
                       ncol = 1 , nrow = 3,
                       common.legend = TRUE,
                       legend="bottom")

figure_2

ggsave("plot_figure_2.tiff", 
       units = "in", 
       width = 4, 
       height = 11, 
       dpi = 600, 
       compression = 'lzw')
###################5.Proportion of Host Viability###########################
#HV stands for Host Viability#

binom_HV <- subset(Mets,  Retrieved == "0" | Retrieved == "1")
bn_HV <- glm(formula = Retrieved ~ NP_treatment,
             family = binomial(link = "logit"),
             data = binom_HV)

summary(bn_HV)
Anova(bn_HV, type = 2)

mean_HV <- stat_summary(data = binom_HV, fun.data = "mean_se",
                       geom = "bar",
                       size = 0.5,
                       width = 0.1,
                       colour = "black")

errorbar_HV <- stat_summary(data = binom_HV,
                           fun.data = "mean_se",
                           geom = "errorbar",
                           width = 0.05,
                           position = position_dodge(width = 0.1),
                           colour = "darkgrey", alpha = 1)

# Plot
HV <- ggplot(binom_HV, aes( NP_treatment, Retrieved, fill = NP_treatment, group = NP_treatment)) + 
  coord_cartesian(ylim=c(0, 1)) +
  scale_y_continuous(breaks=seq(0, 1, 0.25)) +
  scale_fill_grey(start = 1, end = 0) +
  mean_HV +
  errorbar_HV +
  labs( x = "\nNanoplastics concentration",
        y = "\nHost survival probability \n(until day 9 post- inoculation)",
        title = "\nHost viability" ) + 
  theme_classic() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "none")
HV
# Save plot
ggsave("plot_Host_Viability_Individuals.tiff", 
       units = "in", 
       width = 5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')

#####################6.Parasite Reproduction####################################
# PR stand for Parasite Reproduction#


errorbar_PR <- stat_summary(data = InfH,  # Dataframe InfH: only the ones that were infected
                           fun.data = "mean_se",
                           geom = "errorbar",
                           width = 0.1,
                           position = position_dodge(width = 0.6),
                           colour = "black", alpha = 0.7)

lines_PR <- stat_summary(data = InfH, fun.data = "mean_se",
                        geom = "line",
                        position = position_dodge(width=0.6),
                        colour = "black",
                        alpha = 0.3)
mean_PR <- stat_summary(data = InfH, fun.data = "mean_se",
                       geom = "point",
                       position = position_dodge(width=0.6),
                       pch = 21,
                       size = 1.5)

PR <- ggplot(InfH, aes(NP_treatment, Spore_yield, fill = NP_treatment, group = NP_treatment)) +
  scale_y_continuous(breaks=seq(0, 40000, 5000)) +
  geom_jitter(width = 0.1, alpha = 0.5, pch = 20, size = 0.7) +
  scale_fill_grey(start = 1, end = 0) +
  errorbar_PR +
  lines_PR +
  mean_PR +
  labs(x = "\nNanoplastics concentration",
       y = "\nSpore yield per infected host",
       title = "\nParasite reproduction") +
  #Significant difference#
  #annotate("text", x = 1, y= 20500,  label = "a      ", size = 3) + 
  #annotate("text", x = 2, y= 25500, label = "a      ", size = 3) +  
  #annotate("text", x = 3, y= 8300, label = "b      ", size = 3) +   
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "none")

PR
ggsave("plot_Parasite_reproduction.tiff", 
       units = "in", 
       width = 4, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')


model_PR <- aov(Spore_yield ~ NP_treatment, data = InfH)

summary(model_PR)
plot(model_PR)
hist(resid(model_PR))

Anova(model_PR,  type = 2)
TukeyHSD(model_PR)

#% Variance expained 
1840926463 + 2881274756
(1840926463 / 4722201219)*100
(2881274756 / 4722201219)*100



######################7.Proportion of Infection#################################
#PI stands for Parasite infectivity#

binom_PI <- subset(Inf_after_day8, Inf_mets == "1" | Inf_mets == "0")

bn_PI <- glm(formula = Inf_mets ~ NP_treatment,
           family = binomial(link = "logit"), 
           data = binom_PI)

summary(bn_PI)
Anova(bn_PI, type = 2)

mean_PI <- stat_summary(data = binom_PI, fun.data = "mean_se",
                       geom = "bar",
                       size = 0.5,
                       width = 0.1,
                       colour = "black")

errorbar_PI <- stat_summary(data = binom_PI,
                           fun.data = "mean_se",
                           geom = "errorbar",
                           width = 0.05,
                           position = position_dodge(width = 0.3),
                           colour = "darkgrey", alpha = 1)



PI <- ggplot(binom_PI, aes( x = NP_treatment, y = Inf_mets, fill = NP_treatment)) +
  coord_cartesian(ylim=c(0, 1)) +
  scale_y_continuous(breaks=seq(0, 1, 0.25)) +
  scale_fill_grey(start = 1, end = 0) +
  mean_PI +
  errorbar_PI +
  labs(x = "\nNanoplastics concentration", 
       y = "\nProportion of infected host",
       title = "\nParasite infectivity") +
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "none")
PI
ggsave("plot_Parasite_infectivity.tiff", 
       units = "in", 
       width = 5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')


HV # Host vialability plot
PI # Parasite infectivity plot
PR # Parasite reproduction plot


### Combine the graphs for parasite traits####
figure_1  <- ggarrange(HV, PI, PR,
                       labels = c("A", "B", "C"),
                       ncol = 3 , nrow = 1)
figure_1

ggsave("plot_figure_1.tiff", 
       units = "in", 
       width = 14, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')









###Supplement material####

################################Sup_Host Lifespan################################

Sup_errorbar_HL <- stat_summary(data = Successfully_inf,
                            fun.data = "mean_se",
                            geom = "errorbar",
                            width = 0.2,
                            position = position_dodge(width = 0.6),
                            colour = "black", alpha = 0.5)

Sup_lines_HL <- stat_summary(data = Successfully_inf, fun.data = "mean_se",
                         geom = "line",
                         position = position_dodge(width=0.6),
                         colour = "black",
                         alpha = 0.3)

Sup_mean_HL <- stat_summary(data = Successfully_inf, fun.data = "mean_se",
                        geom = "point",
                        position = position_dodge(width=0.6),
                        pch = 21,
                        size = 1.5)

#PLot - Host Lifespan
Sup_A1 <- ggplot(Successfully_inf, aes(NP_treatment, Age_death, fill = NP_treatment)) +
  facet_wrap(~ Inf_treatment) +
  geom_jitter(width = 0.2, alpha = 0.25, pch = 20, size = 0.7) +
  scale_y_continuous(breaks=seq(0, 30, 5)) +
  scale_fill_grey(start = 1, end = 0) +
  Sup_lines_HL + 
  Sup_errorbar_HL +
  Sup_mean_HL + 
  labs( x = "\nNanoplastics concentration",
        y = "\nAge at death (days)",
        title = "\nHost lifespan" ) + 
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "none")
Sup_A1

# Save plot.
ggsave("plot_DAMN_Sup_Host_Lifespan.tiff", 
       units = "in", 
       width = 4.5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')


##############################Sup_Host Fecundity##################################

#Only the ones that reproduced at least once
Sup_Reproduced_at_least1 <- subset(Successfully_inf, Total_juv > 0) 

Sup_errorbar_HF <- stat_summary(data = Sup_Reproduced_at_least1,
                            fun.data = "mean_se",
                            geom = "errorbar",
                            width = 0.2,
                            position = position_dodge(width = 0.6),
                            colour = "black", alpha = 0.5)

Sup_lines_HF <- stat_summary(data = Sup_Reproduced_at_least1, fun.data = "mean_se",
                         geom = "line",
                         position = position_dodge(width=0.6),
                         colour = "black",
                         alpha = 0.3)

Sup_mean_HF <- stat_summary(data = Sup_Reproduced_at_least1, fun.data = "mean_se",
                        geom = "point",
                        position = position_dodge(width=0.6),
                        pch = 21,
                        size = 1.5)

Sup_C1 <- ggplot(Sup_Reproduced_at_least1, aes(NP_treatment, Total_juv, fill = NP_treatment)) +
  facet_wrap(~ Inf_treatment) +
  geom_jitter(width = 0.1, alpha = 0.25, pch = 20, size = 0.7) +
  Sup_lines_HF +
  Sup_errorbar_HF +
  Sup_mean_HF + 
  coord_cartesian(ylim=c(0, 18)) +
  scale_fill_grey(start = 1, end = 0) +
  scale_y_continuous(breaks=seq(0, 20, 4)) +
  labs( x = "\nNanoplastics concentration", 
        y = "\nTotal number of offspring (per host)", 
        title = "\nHost fecundity") +
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(), 
        legend.position = "bottom")

Sup_C1 # Sup_Host fecundity

ggsave("plot_Sup_Host_fecundity.tiff", 
       units = "in", 
       width = 4.5, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')



######################Sup_Host that Reached Maturity##############################


# Create new column; 0 | 1 did not reached maturity | reached maturity.
# maturity = The individuals that reached maturity
Successfully_inf$maturity <-  (Successfully_inf$Age_at_Mat > 0)
Successfully_inf$maturity <- as.numeric(Successfully_inf$maturity,
                         levels = c("TRUE","FALSE"),
                         labels = c("1", "0"))


str(Successfully_inf)

#HRM stands for Host Reached Maturity#

Sup_binom_HRM <- subset(Successfully_inf,  maturity == "0" | maturity == "1")


Sup_bn_HRM <- glm(formula = Sup_binom_HRM$maturity ~ Inf_treatment * NP_treatment, 
              family = binomial(link = "logit"), 
              data = Sup_binom_HRM )
summary(Sup_bn_HRM)


Sup_mean_HRM <- stat_summary(data = Sup_binom_HRM, fun.data = "mean_se",
                         geom = "bar", 
                         size = 0.6,
                         width = 0.2,
                         colour = "black",
                         position = "dodge")

Sup_errorbar_HRM <- stat_summary(data = Sup_binom_HRM,
                             fun.data = "mean_se",
                             geom = "errorbar",
                             width = 0.05,
                             colour = "darkgrey",
                             position = position_dodge(width = 0.2))

# Plot
Sup_B1 <- ggplot(Sup_binom_HRM, aes( Inf_treatment, maturity, fill = NP_treatment)) + 
  coord_cartesian(ylim=c(0, 1)) +
  scale_y_continuous(breaks=seq(0, 1, 0.25)) +
  scale_fill_grey(start = 1, end = 0) +
  Sup_mean_HRM +
  Sup_errorbar_HRM +
  labs( x = "\nNanoplastics concentration",
        y = "\nProportion of Daphnia ",
        title = "\nHost that reached maturity" ) +
  theme_bw() +
  theme(panel.grid.major = element_line(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 15),
        legend.title = element_blank(),
        legend.position = "none")
Sup_B1
# Save plot
ggsave("plot_Sup_Host_R_mat.tiff", 
       units = "in", 
       width = 6, 
       height = 4, 
       dpi = 600, 
       compression = 'lzw')


Sup_figure_3  <- ggarrange(Sup_A1, Sup_B1, Sup_C1,
                       labels = c("A", "B", "C"),
                       ncol = 1 , nrow = 3,
                       common.legend = TRUE,
                       legend="bottom")

Sup_figure_3

ggsave("plot_Sup_figure_3.tiff", 
       units = "in", 
       width = 4, 
       height = 11, 
       dpi = 600, 
       compression = 'lzw')






##################8.Spore yield / Age post infection############################

#InfH$cal <-  (InfH$Spore_yield / InfH$Age_post_inf)
#str(InfH)
#
#errorbar_cal <- stat_summary(data = InfH,  # Dataframe InfH: only the ones that were infected
#                             fun.data = "mean_se",
#                             geom = "errorbar",
#                             width = 0.1,
#                             position = position_dodge(width = 0.6),
#                             colour = "black", alpha = 0.7)
#
#lines_cal <- stat_summary(data = InfH, fun.data = "mean_se",
#                          geom = "line",
#                          position = position_dodge(width=0.6),
#                          colour = "black",
#                          alpha = 0.3)
#mean_cal <- stat_summary(data = InfH, fun.data = "mean_se",
#                         geom = "point",
#                         position = position_dodge(width=0.6),
#                         pch = 21,
#                         size = 1.5)
#
#ggplot(InfH, aes(NP_treatment, cal, fill = NP_treatment, group = NP_treatment)) +
#  scale_y_continuous(breaks=seq(0, 4000, 300)) +
#  geom_jitter(width = 0.1, alpha = 0.5, pch = 20, size = 0.7) +
#  scale_fill_grey(start = 1, end = 0) +
#  errorbar_cal +
#  lines_cal +
#  mean_cal +
#  labs(x = "\nNanoplastics concentration",
#       y = "\nSpore yield/ \nAge post infection",
#       title = "\nParasite reproduction") +
#  theme_bw() +
#  theme(panel.grid.major = element_line(),
#        panel.grid.minor = element_blank(),
#        plot.title = element_text(hjust = 0.5, size = 15),
#        legend.title = element_blank(), 
#        legend.position = "none")
#
#ggsave("plot_Parasite_Age_post_inf.tiff", 
#       units = "in", 
#       width = 6, 
#       height = 4, 
#       dpi = 600, 
#       compression = 'lzw')

#model_cal <- aov(cal ~ NP_treatment, data = InfH)
#
#summary(model_cal)
#plot(model_cal)
#hist(resid(model_cal))
#Anova(model_cal)
#TukeyHSD(model_cal)