### To accompany "Crowding reduces per-capita parasite infection risk in a butterfly host"

# Experimental data analysis

# Clear workspace
rm(list=ls()) 
graphics.off()

# Load necessary libraries 
library(tidyverse)
library(lubridate)
library(survival)
library(survminer)
library(lme4)
library(lmerTest)
library(devtools)
library(sjPlot)
library(coxme)
library(multcomp)
library(EnvStats)
library(ggthemes)
library(multcompView)
library(modelr)

##################################################################################################

###### EXPERIMENT 1 ###########

##################################################################################################


exp1<-read.csv("Experiment1.csv")

dfi<-exp1 %>% 
  mutate(Infected = InfectedLS) %>% 
  filter(Treatment=="infected")

dfii<-dfi %>% 
  filter(lSL5>0)

df1<-exp1 %>% 
  mutate(Infected = InfectedLS) %>% 
  mutate(Success = ifelse(Treatment == "uninfected", "success", "fail")) %>% 
  mutate(Success = ifelse(Infected == 1, "success", Success)) %>% 
  mutate(Dose = factor(Dose)) %>% 
  mutate(Treatment = as.factor(Treatment))


dfnew<-df1 %>% 
  filter(Treatment == "infected" & Success == "fail") 

#removed failed inoculations
df1<-df1 %>% 
  anti_join(dfnew) %>% 
  mutate(NewTreatment = Treatment)

dfnew<-dfnew %>% 
  mutate(NewTreatment = "uninfected", Dose = "0")

df2 <- df1 %>% 
  filter(Treatment == "infected") %>% 
  group_by(TubeID, Density, Dose) %>% 
  summarize(totalInf = sum(Infected,na.rm = T))

df3<-df1 %>% 
  filter(Treatment == "infected") %>% 
  group_by(TubeID, Density, Dose) %>% 
  summarize(totalInfSpores = sum(lSL5, na.rm = T))


#added failed inoculations back 
df1<-df1 %>% bind_rows(dfnew) %>% 
  filter(!is.na(Adult_death_date))

df1$NewTreatment<-factor(df1$NewTreatment, levels = c("uninfected","infected"), labels = c( "control","exposed"))

df1 <- df1 %>% 
  mutate(Infected = InfectedLS)%>% 
  filter(!is.na(Treatment)) 


# Model fitting infection status #
summary(model1<-lmer(Infected ~ Dose * Density + (1|TubeID), data = dfi))
dfi$interaction<-interaction(dfi$Density,dfi$Dose)

# Tukey post-hoc tests
m1<-glmer(Infected~interaction+(1|TubeID),data=dfi, family = "binomial")
summary(glht(m1, mcp(interaction="Tukey")))
cld(glht(m1, mcp(interaction="Tukey")))

# Model fitting parasite load #
summary(model1.55<-lmer(lSL5 ~ Dose * Density + (1|TubeID), data = dfii))
dfii$interaction<-interaction(dfii$Density,dfii$Dose)

# Tukey post-hoc tests
m1.55<-glmer(lSL5~interaction+(1|TubeID),data=dfii)
g1<-summary(glht(m1.55, mcp(interaction="Tukey")))
summary(glht(m1.55, mcp(interaction="Tukey")))
cld(glht(m1.55, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(model1,  model1.55, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"),   digits = 2, digits.p = 2, file = "Exp1Table1.doc")
cld(glht(m1.55, linfct=mcp(interaction="Tukey")))

# Model fitting total number of infected #
summary(modelX<-lm(totalInf ~ Dose * Density, data = df2))
df2$interaction<-interaction(df2$Density,df2$Dose)

# Tukey post-hoc tests
modelX<-lm(totalInf~interaction,data=df2)
summary(glht(modelX, mcp(interaction="Tukey")))
cld(glht(modelX, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(modelX, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp1Table2.doc")

# Model fitting cumulative spore numbers #
summary(modelX2<-lm(totalInfSpores ~ Dose * Density, data = df3))
df3$interaction<-interaction(df3$Density,df3$Dose)

# Summary 
sjPlot::tab_model(modelX2, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"),  digits = 2, digits.p = 2, file = "Exp1Table3.doc")

# Tukey post-hoc tests
modelX2<-lm(totalInfSpores~interaction,data=df3)
summary(glht(modelX2, mcp(interaction="Tukey")))
cld(glht(modelX2, mcp(interaction="Tukey")))

sjPlot::tab_model(modelX2, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp1Table4.doc")

# Model fitting adult lifespan
summary(survmod<-lmer(AdultLongevity ~ Dose * Density + (1|TubeID), data = df1))
df1$interaction<-interaction(df1$Density,df1$Dose)
msurv<-lmer(AdultLongevity~interaction+(1|TubeID),data=df1)

# Tukey post-hoc tests
summary(glht(msurv, mcp(interaction="Tukey")))
cld(glht(msurv, mcp(interaction="Tukey")))
sjPlot::tab_model(survmod, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"),digits = 2, digits.p = 2, file = "Exp1Table5.doc")

 


##################################################################################################

###### EXPERIMENT 2 ###########

##################################################################################################

exp2<-read.csv("Experiment2.csv")

exp2$Dose <- factor(exp2$Dose, levels = c("control", "low", "high"))
exp2$Density <- factor(exp2$Density, levels = c("single", "double", "ten"))

exp2i<-exp2 %>% 
  filter(Treatment == 1) 

df2e <- exp2 %>% 
  filter(Treatment == 1) %>% 
  group_by(TubeID, Density, Dose) %>% 
  summarize(totalInf = sum(InfectedLS,na.rm = T))

df3e<-exp2 %>% 
  filter(Treatment == 1) %>% 
  group_by(TubeID, Density, Dose) %>% 
  summarize(totalInfSpores = sum(lSL5,na.rm = T))


# Model fitting infection status
summary(model1<-lmer(InfectedLS ~ Dose * Density +(1|TubeID), data =exp2i))
exp2i$interaction<-interaction(exp2i$Density,exp2i$Dose)

# Tukey post-hoc tests
m4<-lm(InfectedLS~interaction,data=exp2i, family = "binomial")
summary(glht(m4, mcp(interaction="Tukey")))
cld(glht(m4, mcp(interaction="Tukey")))

# Model fitting parasite load
exp2ii<-exp2i %>% 
  filter(lSL5>0)

summary(model1.2<-lmer(lSL5 ~ Dose*Density+ (1|TubeID), data =exp2ii))

# Tukey post-hoc tests
exp2ii$interaction<-interaction(exp2ii$Density,exp2ii$Dose)
m3<-lmer(lSL5~interaction+(1|TubeID),data=exp2ii)
summary(glht(m3, mcp(interaction="Tukey")))
cld(glht(m3, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(model1,  model1.2, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"),   digits = 2, digits.p = 2,  file = "Exp2Table1.doc")

dfexp <-exp2 %>%  
  mutate(Infected = InfectedLS) %>% 
  mutate(Success = ifelse(Treatment == 0, "success", "fail")) %>% 
  mutate(Success = ifelse(Infected == 1, "success", Success)) 

dfnew2<-dfexp %>% 
  filter(Treatment == 1 & Success == "fail") 

dfexp<-dfexp %>% 
  anti_join(dfnew2)

dfexp<-dfexp %>% 
  mutate(NewTreatment = Treatment)

dfnew2<-dfnew2%>% 
  mutate(NewTreatment = 0, Dose = "control")

dfexp<-dfexp %>% bind_rows(dfnew2) %>% 
  filter(!is.na(Adult_death_date))

dfexp$NewTreatment<-factor(dfexp$NewTreatment, levels = c(0,1), labels = c("control", "exposed"))

dfiii<-dfexp %>% 
  filter(lSL5>0)


# Model fitting total number of infected ####
summary(modelXe<-lm(totalInf ~ Dose * Density, data = df2e))
df2e$interaction<-interaction(df2e$Density,df2e$Dose)

# Tukey post-hoc tests
modelXe<-lm(totalInf~interaction,data=df2e)
summary(glht(modelXe, mcp(interaction="Tukey")))
cld(glht(modelXe, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(modelX2, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp2Table2.doc")


# Model fitting cumulative spore numbers #
summary(modelX2e<-lm(totalInfSpores ~ Dose * Density, data = df3e))
df3e$interaction<-interaction(df3e$Density,df3e$Dose)

# Summary 
sjPlot::tab_model(modelX2e, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp2Table3.doc")

# Tukey post-hoc tests
modelX2e<-lm(totalInfSpores~interaction,data=df3e)
summary(glht(modelX2e, mcp(interaction="Tukey")))
cld(glht(modelX2e, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(modelX2e, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp2Table4.doc")

# Model fitting adult lifespan
summary(model9<-lmer(AdultLongevity~  Dose*Density + (1|TubeID), data = dfexp))
dfexp$interaction<-interaction(dfexp$Density,dfexp$Dose)

dfexpa<-dfexp %>% 
  filter(!is.na(AdultLongevity))

# Tukey post-hoc tests
m3.1<-lmer(AdultLongevity ~ interaction + (1|TubeID),data=dfexpa)
summary(glht(m3.1, mcp(interaction="Tukey")))
cld(glht(m3.1, mcp(interaction="Tukey")))

# Summary 
sjPlot::tab_model(model9, CSS = css_theme("cells"),show.est = TRUE,show.se = TRUE, show.ci = FALSE, show.stat = TRUE, show.df = TRUE,   col.order = c("est", "se", "stat", "df.error","p"), digits = 2, digits.p = 2, file = "Exp2Table5.doc")


####

df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)



df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)




df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)




df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)




df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)




df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)




df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores per tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)



df3n2<-df3 %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores\nper tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,25)) +
  scale_fill_manual(values = c("#f4f4f4", "#c4c2c2"))+  facet_grid(~Density)


df3n2e<-df3e %>% 
  group_by(Density, Dose) %>% 
  summarize(totalInfSporesmean = mean(totalInfSpores, na.rm = T), se = sd(totalInfSpores)/sqrt(n()), n=n())

ggplot(df3n2e, aes(x = as.factor(Dose), y = totalInfSporesmean, fill = Dose)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8, color = "black") +
  geom_errorbar(aes(ymin = totalInfSporesmean - se, ymax = totalInfSporesmean + se),
                position = position_dodge(width = 0.9), width = 0.25) +
  labs(x = "Caterpillar density", y = "Cumulative log10 parasite\nspore load per tube") +
  theme_classic()+
  labs(x = "", title = "", y = "Cumulative spores\nper tube", fill = "")+
  theme(strip.text.x = element_text(size = 15, face = "bold.italic"),
        axis.text=element_text(size=15),
        title =element_text(size=16, face='bold'),
        axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=15),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=15),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = "none",
        panel.spacing.x = unit(0, "pt"), strip.background = element_rect(color="black", size=0.5, linetype="solid"), plot.title = element_text(hjust = 0.5))+
  scale_y_continuous(expand = c(0,0),limits = c(0,60)) +
  scale_fill_manual(values = c("#8e8e8e", "#595858"))+  facet_grid(~Density)