### Worker microcolonies (male larvae) of Bombus impatiens reaered with exposure to Nosema bombi or Chlorothalonil. Questions: How does single or co-exposure affect development time? How does single or co-exposure affect adult size? How does single or co-exposure affect adult protein? How does co-exposure to chlorothalonil affect infection and spore loads of Nosema? How does single or co-exposure affect survival? ### Clear workspace rm(list = ls(all = TRUE)) ### Set for correct SS (change for survival) options(contrasts = c("contr.sum", "contr.poly")) ### Load packages (if not present then need to install.packages()). library(MASS) library(lme4) library(ggplot2) library(car) library(fitdistrplus) library(DHARMa) library(gridExtra) library(emmeans) library(bbmle) library(GGally) library(lattice) library(survival) library(effects) library(glmmTMB) library(coxme) ###Read in data dat37<-read.csv("Chlorothalonil_Nosema_individual_data.csv",header=T) ###Check data and adjust as necessary summary(dat37) head(dat37) dim(dat37) str(dat37) ### Specify correct data types dat37$source_colony<-as.factor(dat37$source_colony) dat37$microcolony<-as.factor(dat37$microcolony) dat37$nosema<-as.factor(dat37$nosema) dat37$chlorothalonil<-as.factor(dat37$chlorothalonil) dat37$development<-as.numeric(as.character(dat37$development)) dat37$survival<-as.numeric(as.character(dat37$survival)) dat37$total_inf<-as.numeric(as.character(dat37$total_inf)) dat37$extra_spores<-as.numeric(as.character(dat37$extra_spores)) dat37$intra_bi<-as.factor(dat37$intra_bi) dat37$extra_bi<-as.factor(dat37$extra_bi) dat37$infection<-as.factor(dat37$infection) ###Recheck data summary(dat37) head(dat37) dim(dat37) str(dat37) ### Create subsets for different purposes. dat37surv<-droplevels(subset(dat37,purpose=="survival")) summary(dat37surv) dim(dat37surv) dat37infAll<-droplevels(subset(dat37,purpose=="biochem"& nosema=="1"|purpose=="qPCR" )) summary(dat37infAll) dim(dat37infAll) dat37bio<-droplevels(subset(dat37,purpose=="biochem"&!is.na(protein_ug_ml))) summary(dat37bio) dim(dat37bio) ### DEVELOPMENT (using all data) ### (a) Across all treatment combinations ggplot(dat37, aes(x = development)) + geom_density() ggplot(dat37, aes(x = development)) + geom_density() + facet_wrap(nosema~chlorothalonil) ggplot(dat37, aes(y = development, x = nosema, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,38) ### Test distributions ###Normal qqp(dat37$development, "norm") ### log normal qqp(dat37$development, "lnorm") ### Gamma gamma <- fitdistr(dat37$development, "gamma") qqp(dat37$development, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]]) plotdist(dat37$development, histo = TRUE, demp = TRUE) descdist(dat37$development, discrete=F, boot=500) ### Gamma and normal both have potential issues, but fit much better than lnorm. # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models Dev_mod1<-lmer(development~nosema*chlorothalonil+(1|microcolony:source_colony), data=dat37, REML=T) ###Graphing the residuals of model. plot(fitted(Dev_mod1), residuals(Dev_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(Dev_mod1), residuals(Dev_mod1))) Dev_mod1_simres <- simulateResiduals(Dev_mod1) plot(Dev_mod1_simres) #### Residuals of normal good summary(Dev_mod1) drop1(Dev_mod1,test="Chisq") Dev_mod2<-update(Dev_mod1,.~. -nosema:chlorothalonil) summary(Dev_mod2) drop1(Dev_mod2,test="Chisq") Dev_mod3<-update(Dev_mod2,.~. -chlorothalonil) summary(Dev_mod3) drop1(Dev_mod3,test="Chisq") Dev_mod4<-update(Dev_mod3,.~. -nosema) summary(Dev_mod4) AIC(Dev_mod1, Dev_mod2, Dev_mod3, Dev_mod4) AICtab(Dev_mod1, Dev_mod2, Dev_mod3, Dev_mod4) ### Confirm stepwise LRT anova(Dev_mod1, Dev_mod2, Dev_mod3, Dev_mod4) ### Best model on AIC is with no fixed effects ### Check on random effects. If effect needed, get with anova of nested models removing ranef Dev_ran<-ranef(Dev_mod1,condVar=TRUE) a <- dotplot(Dev_ran)[["microcolony:source_colony"]] grid.arrange(a, nrow=1) ### emmeans to report (or plot) based on nosema and chlorothalonil combination Dev_emmeans<-as.data.frame(emmeans(Dev_mod1, ~nosema:chlorothalonil)) ### Sample sizes Dev_ss<-as.data.frame(xtabs(~nosema+chlorothalonil, model.frame(Dev_mod1))) Dev_emmeans<-cbind(Dev_emmeans,Dev_ss[3]) Dev_emmeans ### Plot trt<-c("Control","Exposed") chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") dev_plot<-ggplot(Dev_emmeans, aes(x = nosema, y = emmean, color = chlorothalonil)) + geom_jitter(data=dat37, aes(x = nosema, y = development, col = chlorothalonil),position = position_jitterdodge(jitter.width = 0.4, jitter.height = 0, dodge.width = 0.6, seed = NA), alpha=0.1)+ geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=0.75, position = position_dodge(width=0.6))+ ylab("Development time (days \u00b1 95% C.I.)") + xlab (expression(paste(italic("Nosema bombi")," treatment")))+ ylim(0,40) + scale_color_manual(labels = c("Control","Exposed"), values = c("#C0C0C0","#2F4F4F"), name="Chlorothalonil:") + scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") dev_plot ### (b) Within Nosema exposed infected v non-infected based on infection evidence (qPCR or spore) dat37nos<-droplevels(subset(dat37, nosema==1)) summary(dat37nos) dim(dat37nos) ggplot(dat37nos, aes(x = development)) + geom_density() + facet_wrap(infection~chlorothalonil) ggplot(dat37nos, aes(y = development, x = infection, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,38) ### Test distributions ###Normal qqp(dat37nos$development, "norm") ### log normal qqp(dat37nos$development, "lnorm") ### Gamma gamma <- fitdistr(dat37nos$development, "gamma") qqp(dat37nos$development, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]]) plotdist(dat37nos$development, histo = TRUE, demp = TRUE) descdist(dat37nos$development, discrete=F, boot=500) ### Gamma and normal both have issues, but fit much better than lnorm. # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models NosDev_mod1<-lmer(development~infection*chlorothalonil+(1|microcolony:source_colony), data=dat37nos, REML=T) ### Check diagnostics plot(fitted(NosDev_mod1), residuals(NosDev_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(NosDev_mod1), residuals(NosDev_mod1))) NosDev_mod1_simres <- simulateResiduals(NosDev_mod1) plot(NosDev_mod1_simres) summary(NosDev_mod1) drop1(NosDev_mod1,test="Chisq") NosDev_mod2<-update(NosDev_mod1,.~. -infection:chlorothalonil) summary(NosDev_mod2) drop1(NosDev_mod2,test="Chisq") NosDev_mod3<-update(NosDev_mod2,.~. -chlorothalonil) summary(NosDev_mod3) drop1(NosDev_mod3,test="Chisq") NosDev_mod4<-update(NosDev_mod3,.~. -infection) summary(NosDev_mod4) AIC(NosDev_mod1, NosDev_mod2, NosDev_mod3, NosDev_mod4) AICtab(NosDev_mod1, NosDev_mod2, NosDev_mod3, NosDev_mod4) ### Confirm stepwise LRT anova(NosDev_mod1, NosDev_mod2, NosDev_mod3, NosDev_mod4) ### Support for best model with infection ### Check on random effects NosDev_ran<-ranef(NosDev_mod1,condVar=TRUE) a <- dotplot(NosDev_ran)[["microcolony:source_colony"]] grid.arrange(a, nrow=1) ### emmeans to report (or plot) based on nosema and chlorothalonil combination NosDev_emmeans<-as.data.frame(emmeans(NosDev_mod1, ~infection:chlorothalonil)) ### Sample sizes NosDev_ss<-as.data.frame(xtabs(~infection+chlorothalonil, model.frame(NosDev_mod1))) NosDev_emmeans<-cbind(NosDev_emmeans,NosDev_ss[3]) NosDev_emmeans ### Sample sizes NosDev_ss<-as.data.frame(xtabs(~infection+chlorothalonil, model.frame(NosDev_mod1))) NosDev_emmeans<-cbind(NosDev_emmeans,NosDev_ss[3]) NosDev_emmeans ### emmeans to report (or plot) based on infection NosDevInf_emmeans<-emmeans(NosDev_mod3, pairwise~infection) plot(NosDevInf_emmeans, comparisons=T) NosDevInf_emmeans<-as.data.frame(emmeans(NosDev_mod3, ~infection)) ### Sample sizes NosDevInf_ss<-as.data.frame(xtabs(~infection, model.frame(NosDev_mod3))) NosDevInf_emmeans<-cbind(NosDevInf_emmeans,NosDevInf_ss[2]) NosDevInf_emmeans ### Plot trt<-c("Exposed only","Infected") chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") nos_dev_plot<-ggplot(NosDevInf_emmeans, aes(x = infection, y = emmean)) + geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=0.75, position = position_dodge(width=0.6))+ ylab("Development time (days \u00b1 95% C.I.)") + xlab (expression(paste(italic("Nosema bombi")," status")))+ ylim(0,38) + scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="none") nos_dev_plot #################################################################################################################################################### ### Size - wing size as a surrogate - (using all data) ggplot(dat37, aes(x = body_size)) + geom_density() ggplot(dat37, aes(x = body_size)) + geom_density() + facet_wrap(nosema~chlorothalonil) ggplot(dat37, aes(y = body_size, x = nosema, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,3.5) ### Test distributions # qqp requires estimates of the parameters of the negative binomial, Poisson # and gamma distributions. You can generate estimates using the fitdistr # function. Save the output and extract the estimates of each parameter as below. ###Normal qqp(dat37$body_size, "norm") plotdist(dat37$body_size, histo = TRUE, demp = TRUE) descdist(dat37$body_size, discrete=F, boot=500) ### Normal looks good # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models Size_mod1<-lmer(body_size~nosema*chlorothalonil+ (1|microcolony:source_colony), data=dat37, REML=T) ###Graphing the residuals of model. plot(fitted(Size_mod1), residuals(Size_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(Size_mod1), residuals(Size_mod1))) Size_mod1_simres <- simulateResiduals(Size_mod1) plot(Size_mod1_simres) summary(Size_mod1) drop1(Size_mod1,test="Chisq") Size_mod2<-update(Size_mod1,.~. -nosema:chlorothalonil) summary(Size_mod2) drop1(Size_mod2,test="Chisq") Size_mod3<-update(Size_mod2,.~. -chlorothalonil) summary(Size_mod3) drop1(Size_mod3,test="Chisq") Size_mod4<-update(Size_mod3,.~. -nosema) summary(Size_mod4) AIC(Size_mod1, Size_mod2, Size_mod3, Size_mod4) AICtab(Size_mod1, Size_mod2, Size_mod3, Size_mod4) anova(Size_mod1, Size_mod2, Size_mod3, Size_mod4) ### Best model on AIC is with no fixed effects ### Check random effects Size_ran<-ranef(Size_mod1,condVar=TRUE) a <- dotplot(Size_ran)[["microcolony:source_colony"]] grid.arrange(a, nrow=1) Size_emmeans<-as.data.frame(emmeans(Size_mod1, ~nosema:chlorothalonil)) ### Sample sizes Size_ss<-as.data.frame(xtabs(~nosema+chlorothalonil, model.frame(Size_mod1))) Size_emmeans<-cbind(Size_emmeans,Size_ss[3]) Size_emmeans ### Plot trt<-c("Control","Exposed") chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") size_plot<-ggplot(Size_emmeans, aes(x = nosema, y = emmean, color = chlorothalonil)) + geom_jitter(data=dat37, aes(x = nosema, y = body_size, col = chlorothalonil),position = position_jitterdodge(jitter.width = 0.4, jitter.height = 0, dodge.width = 0.6, seed = NA), alpha=0.1)+ geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=0.75, position = position_dodge(width=0.6))+ ylab("Size based on radial cell length\n (mm \u00b1 95% C.I.)") + xlab (expression(paste(italic("Nosema bombi")," treatment")))+ ylim(0,3.5) + ### To break one of my pet peaves or not? ### ylim(2,3.5) + scale_color_manual(labels = c("Control","Exposed"), values = c("#C0C0C0","#2F4F4F"), name="Chlorothalonil:") + scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") size_plot ### (b) Within Nosema exposed infected v non-infected based on evidence of infection ggplot(dat37nos, aes(x = body_size)) + geom_density() + facet_wrap(infection~chlorothalonil) ggplot(dat37nos, aes(y = body_size, x = infection, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,3.5) ### Test distributions ###Normal qqp(dat37nos$body_size, "norm") plotdist(dat37nos$body_size, histo = TRUE, demp = TRUE) descdist(dat37nos$body_size, discrete=F, boot=500) ### Normal good # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models NosSize_mod1<-lmer(body_size~infection*chlorothalonil+(1|microcolony:source_colony), data=dat37nos, REML=T) ### Check diagnostics plot(fitted(NosSize_mod1), residuals(NosSize_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(NosSize_mod1), residuals(NosSize_mod1))) NosSize_mod1_simres <- simulateResiduals(NosSize_mod1) plot(NosSize_mod1_simres) summary(NosSize_mod1) drop1(NosSize_mod1,test="Chisq") NosSize_mod2<-update(NosSize_mod1,.~. -infection:chlorothalonil) summary(NosSize_mod2) drop1(NosSize_mod2,test="Chisq") NosSize_mod3<-update(NosSize_mod2,.~. -chlorothalonil) summary(NosSize_mod3) drop1(NosSize_mod3,test="Chisq") NosSize_mod4<-update(NosSize_mod3,.~. -infection) summary(NosSize_mod4) AIC(NosSize_mod1, NosSize_mod2, NosSize_mod3, NosSize_mod4) AICtab(NosSize_mod1, NosSize_mod2, NosSize_mod3, NosSize_mod4) ### Confirm stepwise LRT anova(NosSize_mod1, NosSize_mod2, NosSize_mod3, NosSize_mod4) ### Check on random effects NosSize_ran<-ranef(NosSize_mod1,condVar=TRUE) a <- dotplot(NosSize_ran)[["microcolony:source_colony"]] a ### emmeans to report (or plot) based on nosema and chlorothalonil combination NosSize_emmeans<-as.data.frame(emmeans(NosSize_mod1, ~infection:chlorothalonil)) ### Sample sizes NosSize_ss<-as.data.frame(xtabs(~infection+chlorothalonil, model.frame(NosSize_mod1))) NosSize_emmeans<-cbind(NosSize_emmeans,NosSize_ss[3]) NosSize_emmeans ### emmeans to report (or plot) based on infection NosSizeInf_emmeans<-emmeans(NosSize_mod3, pairwise~infection) plot(NosSizeInf_emmeans, comparisons=T) NosSizeInf_emmeans<-as.data.frame(emmeans(NosSize_mod3, ~infection)) ### Sample sizes NosSizeInf_ss<-as.data.frame(xtabs(~infection, model.frame(NosSize_mod3))) NosSizeInf_emmeans<-cbind(NosSizeInf_emmeans,NosSizeInf_ss[2]) NosSizeInf_emmeans ### Plot trt<-c("Exposed only","Infected") # chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") nos_size_plot<-ggplot(NosSizeInf_emmeans, aes(x = infection, y = emmean)) + geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=0.75, position = position_dodge(width=0.6))+ ylab("Size based on radial cell length\n (mm \u00b1 95% C.I.)") + xlab (expression(paste(italic("Nosema bombi")," status")))+ # ylim(0,3.5) + ylim(2.7,2.9)+ scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="none") nos_size_plot #################################################################################################################################################### ### Infection - (using InfAll data set) ### Plot to investigate data (histograms, density plot, boxplot) ggplot(dat37infAll, aes(x = total_inf, col=chlorothalonil)) + geom_histogram()+ facet_wrap(~chlorothalonil) ggplot(dat37infAll, aes(x = total_inf, col=chlorothalonil)) + geom_density()+ facet_wrap(~chlorothalonil) ggplot(dat37infAll, aes(y = total_inf, x=chlorothalonil)) + geom_boxplot() ggplot(subset(dat37infAll,total_inf>0), aes(y = total_inf, x=chlorothalonil)) + geom_boxplot() ggplot(subset(dat37infAll,total_inf>0), aes(x = total_inf, col=chlorothalonil)) + geom_histogram()+ facet_wrap(~chlorothalonil) ggplot(subset(dat37infAll,total_inf>0), aes(x = extra_spores, col=chlorothalonil)) + geom_density()+ facet_wrap(~chlorothalonil) ggplot(subset(dat37infAll,extra_spores>0), aes(y = prop_spores, x=chlorothalonil)) + geom_boxplot() ### Get descriptive statistics for infection and relationship between total (intra) infection and extra spores ### Total proportion of infections by qPCR mean(as.numeric(as.character(dat37infAll$intra_bi)), na.rm=T) ### 0.59 ### Relationship between infection by qPCR and extra spores to report R^2, F value, p, etc. log transform both for fit. mod_inf<-lm(log(extra_spores)~log(total_inf),data=subset(dat37infAll,extra_spores>0)) ### Check model fit par(mfrow=c(2,2)) plot(mod_inf) par(mfrow=c(1,1)) plot(fitted(mod_inf), residuals(mod_inf), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(mod_inf), residuals(mod_inf))) mod_inf_simres <- simulateResiduals(mod_inf) plot(mod_inf_simres) ### Get relevant output summary(mod_inf) Anova(mod_inf, test="F", type="III") inf_spore_plot<-ggplot(subset(dat37infAll, extra_spores>0), aes(x=log(total_inf), y=log(extra_spores)))+ geom_point()+ geom_smooth(method=lm, col="darkgreen")+ ylab(expression(paste("Log ", italic("Nosema bombi")," spores"))) + xlab (expression(paste("Log ", italic("Nosema bombi")," total infection")))+ theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="none") inf_spore_plot #### Infection 0/1 ggplot(dat37infAll, aes(y = as.numeric(infection)-1, x = chlorothalonil)) + geom_bar(stat = "summary", fun.y = "mean")+ ylim(0,1) inf_bi_mod1<-glmer(infection~ body_size*chlorothalonil+ (1|microcolony:source_colony), data=dat37infAll, family="binomial") inf_bi_mod1_simres <- simulateResiduals(inf_bi_mod1) plot(inf_bi_mod1_simres) inf_bi_ran<-ranef(inf_bi_mod1,condVar=TRUE) dotplot(inf_bi_ran)[["microcolony:source_colony"]] summary(inf_bi_mod1) drop1(inf_bi_mod1,test="Chisq") inf_bi_mod2<-update(inf_bi_mod1,.~. -body_size:chlorothalonil) summary(inf_bi_mod2) drop1(inf_bi_mod2,test="Chisq") inf_bi_mod3<-update(inf_bi_mod2,.~. -chlorothalonil) summary(inf_bi_mod3) drop1(inf_bi_mod3,test="Chisq") inf_bi_mod4<-update(inf_bi_mod3,.~. -body_size) summary(inf_bi_mod4) AIC(inf_bi_mod1, inf_bi_mod2, inf_bi_mod3, inf_bi_mod4) AICtab(inf_bi_mod1, inf_bi_mod2, inf_bi_mod3, inf_bi_mod4) ### Confirm stepwise LRT anova(inf_bi_mod1, inf_bi_mod2, inf_bi_mod3, inf_bi_mod4) ### Get effect of size size_fit <- as.data.frame(effect('body_size', inf_bi_mod3, xlevels = 100)) ### Plot Size_inf_plot<-ggplot() + geom_point(aes(body_size, as.numeric(as.character(infection))), dat37infAll, size=2, alpha=0.4)+ geom_line(aes(body_size, fit), size_fit, size = 1, col = 'red') + geom_ribbon(aes(body_size, ymin = lower, ymax = upper), size_fit, alpha = 0.1) + ylab("Estimated proportion with infection (\u00b1 95% C.I.)") + xlab ("Size based on radial cell length (mm)")+ scale_size_area()+ theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15)) Size_inf_plot ### emmeans to report (or plot) based on chlorothalonil Chl_inf_emmeans<-as.data.frame(emmeans(inf_bi_mod1, ~chlorothalonil, type="response")) ### Sample sizes Chl_inf_ss<-as.data.frame(xtabs(~chlorothalonil, model.frame(inf_bi_mod1))) Chl_inf_emmeans<-cbind(Chl_inf_emmeans,Chl_inf_ss[2]) Chl_inf_emmeans #### Total infection inf<-subset(dat37infAll,!is.na(total_inf)) plotdist(inf$total_inf, histo = TRUE, demp = TRUE) descdist(inf$total_inf, discrete=T, boot=500) dat37inf<-droplevels(subset(dat37infAll,infection==1)) dim(dat37inf) summary(dat37inf) head(dat37inf) str(dat37inf) plotdist(dat37inf$total_inf, histo = TRUE, demp = TRUE) descdist(dat37inf$total_inf, discrete=T, boot=500) Infnb = glmmTMB(total_inf~ body_size*chlorothalonil+ (1|microcolony:source_colony), data=dat37inf, family="nbinom1") ### Negative binomial fits (checked with nbimon2, poisson, log linked, log transformed...all others have fit issues) Infnb_mod1_simres <- simulateResiduals(Infnb) plot(Infnb_mod1_simres) summary(Infnb) drop1(Infnb, test="Chisq") Infnb2<-update(Infnb,.~. -body_size:chlorothalonil) summary(Infnb2) drop1(Infnb2, test="Chisq") Infnb3<-update(Infnb2,.~. -body_size) summary(Infnb3) #### Check with different optimizer due to warning of false convergence https://cran.r-project.org/web/packages/glmmTMB/vignettes/troubleshooting.html Infnb_check = glmmTMB(total_inf~ chlorothalonil+ (1|microcolony:source_colony), data=dat37inf, family="nbinom1",control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))) summary(Infnb_check) drop1(Infnb3, test="Chisq") Infnb4<-update(Infnb3,.~. -chlorothalonil) AIC(Infnb, Infnb2, Infnb3, Infnb4) AICtab(Infnb, Infnb2, Infnb3, Infnb4) anova(Infnb, Infnb2, Infnb3, Infnb4) ### emmeans to report (or plot) based on chlorothalonil Chl_tot_inf_emmeans<-as.data.frame(emmeans(Infnb3, ~chlorothalonil, type="response")) ### Sample sizes Chl_tot_inf_ss<-as.data.frame(xtabs(~chlorothalonil, model.frame(Infnb3))) Chl_tot_inf_emmeans<-cbind(Chl_tot_inf_emmeans,Chl_tot_inf_ss[2]) Chl_tot_inf_emmeans Chl_tot_inf_plot<-ggplot(Chl_tot_inf_emmeans, aes(x = chlorothalonil, y = response)) + geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=.75, position = position_dodge(width=0.6))+ ylab("Total infection intensity \u00b1 95% C.I.") + xlab ("Chlorothalonil Treatment")+ ylim(0,100000000)+ scale_x_discrete(labels = c("Control", "Exposed")) + #scale_color_manual(labels = c("Control","Chlorothalonil"), values = c("#C0C0C0","#2F4F4F"), name="Treatment:") + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") Chl_tot_inf_plot ####### Spores in those infected individuals ggplot(dat37inf, aes(y = extra_spores, x=chlorothalonil)) + geom_boxplot() dat37inf_sp<-droplevels(subset(dat37inf,extra_spores>0)) dim(dat37inf_sp) ggplot(dat37inf, aes(y = as.numeric(extra_bi)-1, x = chlorothalonil)) + geom_bar(stat = "summary", fun.y = "mean")+ ylim(0,1) ggplot(dat37inf_sp, aes(y = extra_spores, x=chlorothalonil)) + geom_boxplot() #### Spores 0/1 ### Total proportion infected with spores mean(as.numeric(as.character(dat37inf$extra_bi)), na.rm=T) ### 0.69 sp_bi_mod1<-glmer(extra_bi~ body_size*chlorothalonil+ (1|microcolony:source_colony), data=dat37inf, family="binomial") sp_bi_mod1_simres <- simulateResiduals(sp_bi_mod1) plot(sp_bi_mod1_simres) sp_bi_ran<-ranef(sp_bi_mod1,condVar=TRUE) dotplot(sp_bi_ran)[["microcolony:source_colony"]] summary(sp_bi_mod1) drop1(sp_bi_mod1,test="Chisq") sp_bi_mod2<-update(sp_bi_mod1,.~. -body_size:chlorothalonil) summary(sp_bi_mod2) drop1(sp_bi_mod2,test="Chisq") sp_bi_mod3<-update(sp_bi_mod2,.~. -body_size) summary(sp_bi_mod3) drop1(sp_bi_mod3,test="Chisq") sp_bi_mod4<-update(sp_bi_mod3,.~. -chlorothalonil) summary(sp_bi_mod4) AIC(sp_bi_mod1, sp_bi_mod2, sp_bi_mod3, sp_bi_mod4) AICtab(sp_bi_mod1, sp_bi_mod2, sp_bi_mod3, sp_bi_mod4) ### Confirm stepwise LRT anova(sp_bi_mod1, sp_bi_mod2, sp_bi_mod3, sp_bi_mod4) ### emmeans to report (or plot) based on chlorothalonil Chl_sp_emmeans<-as.data.frame(emmeans(sp_bi_mod3, ~chlorothalonil, type="response")) ### Sample sizes Chl_sp_ss<-as.data.frame(xtabs(~chlorothalonil, model.frame(sp_bi_mod3))) Chl_sp_emmeans<-cbind(Chl_sp_emmeans,Chl_sp_ss[2]) Chl_sp_emmeans Chl_sp_emmeans_plot<-ggplot(Chl_sp_emmeans, aes(x = chlorothalonil, y = prob)) + geom_pointrange(aes(ymin = asymp.LCL, ymax = asymp.UCL),size=.75, position = position_dodge(width=0.6))+ ylab(expression(paste(italic("Nosema bombi")," infection probability \u00b1 95% C.I."))) + xlab ("Chlorothalonil Treatment")+ ylim(0,1)+ scale_x_discrete(labels = c("Control","Exposed")) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") Chl_sp_emmeans_plot ###Number of spores plotdist(dat37inf_sp$extra_spores, histo = TRUE, demp = TRUE) descdist(dat37inf_sp$extra_spores, discrete=T, boot=500) #### Note: trying to include infection intensity as a covariate in the following does not work as models fail Sp_nb = glmmTMB(extra_spores~ body_size*chlorothalonil+ (1|microcolony:source_colony), data=dat37inf_sp, family="nbinom1") ### Negative binomial fits (checked with nbimon2, poisson, log linked, log transformed...all others have fit issues or not as good as nbinom1) plot(fitted(Sp_nb), residuals(Sp_nb), xlab = "Fitted Values", ylab = "Residuals", main ="Negative binomial") abline(h = 0, lty = 2) Sp_nb_mod1_simres <- simulateResiduals(Sp_nb) plot(Sp_nb_mod1_simres) summary(Sp_nb) drop1(Sp_nb, test="Chisq") Sp_nb2<-update(Sp_nb,.~. -body_size:chlorothalonil) summary(Sp_nb2) drop1(Sp_nb2, test="Chisq") Sp_nb3<-update(Sp_nb2,.~. -body_size) summary(Sp_nb3) drop1(Sp_nb3, test="Chisq") Sp_nb4<-update(Sp_nb3,.~. -chlorothalonil) summary(Sp_nb4) AIC(Sp_nb, Sp_nb2, Sp_nb3, Sp_nb4) AICtab(Sp_nb, Sp_nb2, Sp_nb3, Sp_nb4) anova(Sp_nb, Sp_nb2, Sp_nb3, Sp_nb4) ### emmeans to report (or plot) based on chlorothalonil Chl_tot_sp_emmeans<-as.data.frame(emmeans(Sp_nb3, ~chlorothalonil, type="response")) plot(emmeans(Sp_nb3, pairwise~chlorothalonil, type="response"), comparisons=T) ### Sample sizes Chl_tot_sp_ss<-as.data.frame(xtabs(~chlorothalonil, model.frame(Sp_nb3))) Chl_tot_sp_emmeans<-cbind(Chl_tot_sp_emmeans,Chl_tot_sp_ss[2]) Chl_tot_sp_emmeans Chl_tot_sp_emmeans_plot<-ggplot(Chl_tot_sp_emmeans, aes(x = chlorothalonil, y = response)) + geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=.75, position = position_dodge(width=0.6))+ ylab(expression(paste(italic("Nosema bombi")," spores per bee \u00b1 95% C.I."))) + xlab ("Chlorothalonil Treatment")+ ylim(0,5000000)+ scale_x_discrete(labels = c("Control","Exposed")) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") Chl_tot_sp_emmeans_plot #### Look at proportion spores ggplot(dat37inf_sp, aes(y = prop_spores, x=chlorothalonil)) + geom_boxplot() ggplot(dat37inf_sp, aes(x = prop_spores, fill=chlorothalonil)) + geom_histogram() plotdist(dat37inf_sp$prop_spores, histo = TRUE, demp = TRUE) descdist(dat37inf_sp$prop_spores, discrete=F, boot=500) prop_mod<-lmer(prop_spores~ body_size*chlorothalonil+ (1|microcolony:source_colony), data=dat37inf_sp, REML=T) plot(fitted(prop_mod), residuals(prop_mod), xlab = "Fitted Values", ylab = "Residuals", main ="Negative binomial") abline(h = 0, lty = 2) prop_mod_simres <- simulateResiduals(prop_mod) plot(prop_mod_simres) prop_mod_ran<-ranef(prop_mod,condVar=TRUE) dotplot(prop_mod_ran)[["microcolony:source_colony"]] summary(prop_mod) drop1(prop_mod, test="Chisq") prop_mod2<-update(prop_mod,.~. -body_size:chlorothalonil) summary(prop_mod2) drop1(prop_mod2, test="Chisq") prop_mod3<-update(prop_mod2,.~. -body_size) summary(prop_mod3) drop1(prop_mod3, test="Chisq") prop_mod4<-update(prop_mod3,.~. -chlorothalonil) summary(prop_mod4) AIC(prop_mod, prop_mod2, prop_mod3) AICtab(prop_mod, prop_mod2, prop_mod3) anova(prop_mod, prop_mod2, prop_mod3, prop_mod4) ### emmeans to report (or plot) based on chlorothalonil Chl_prop_sp_emmeans<-as.data.frame(emmeans(prop_mod3, ~chlorothalonil, type="response")) plot(emmeans(prop_mod3, pairwise~chlorothalonil, type="response"), comparisons=T) ### Sample sizes Chl_prop_sp_ss<-as.data.frame(xtabs(~chlorothalonil, model.frame(prop_mod3))) Chl_prop_sp_emmeans<-cbind(Chl_prop_sp_emmeans,Chl_tot_sp_ss[2]) Chl_prop_sp_emmeans Chl_prop_sp_emmeans_plot<-ggplot(Chl_prop_sp_emmeans, aes(x = chlorothalonil, y = emmean)) + geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),size=.75, position = position_dodge(width=0.6))+ ylab("Proportion of spore production \u00b1 95% C.I.") + xlab ("Chlorothalonil treatment")+ ylim(0,0.06)+ scale_x_discrete(labels = c("Control", "Exposed")) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") Chl_prop_sp_emmeans_plot ##### Total protein ggplot(dat37bio, aes(x = protein_ug_ml)) + geom_density() ggplot(dat37bio, aes(x = protein_ug_ml)) + geom_density() + facet_wrap(nosema~chlorothalonil) ggplot(dat37bio, aes(y = protein_ug_ml, x = nosema, fill=chlorothalonil)) + geom_boxplot() qqp(dat37bio$protein_ug_ml, "norm") plotdist(dat37bio$protein_ug_ml, histo = TRUE, demp = TRUE) descdist(dat37bio$protein_ug_ml, discrete=F, boot=500) ### Normal looks ok, but check # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models Prot_mod1<-lmer(protein_ug_ml~body_size*nosema*chlorothalonil+ (1|microcolony:source_colony), data=dat37bio, REML=T) ###Graphing the residuals of model. plot(fitted(Prot_mod1), residuals(Prot_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(Prot_mod1), residuals(Prot_mod1))) Prot_mod1_simres <- simulateResiduals(Prot_mod1) plot(Prot_mod1_simres) ### Check random effects Prot_ran<-ranef(Prot_mod1,condVar=TRUE) dotplot(Prot_ran)[["microcolony:source_colony"]] summary(Prot_mod1) drop1(Prot_mod1,test="Chisq") Prot_mod2<-update(Prot_mod1,.~. -body_size:nosema:chlorothalonil) summary(Prot_mod2) drop1(Prot_mod2,test="Chisq") Prot_mod3<-update(Prot_mod2,.~. -nosema:chlorothalonil) summary(Prot_mod3) drop1(Prot_mod3,test="Chisq") Prot_mod4<-update(Prot_mod3,.~. -body_size:chlorothalonil) summary(Prot_mod4) drop1(Prot_mod4,test="Chisq") Prot_mod5<-update(Prot_mod4,.~. -chlorothalonil) summary(Prot_mod5) drop1(Prot_mod5,test="Chisq") Prot_mod6<-update(Prot_mod5,.~. -body_size:nosema) summary(Prot_mod6) drop1(Prot_mod6,test="Chisq") Prot_mod7<-update(Prot_mod6,.~. -nosema) summary(Prot_mod7) drop1(Prot_mod7,test="Chisq") Prot_mod8<-update(Prot_mod7,.~. -body_size) summary(Prot_mod8) AIC(Prot_mod1, Prot_mod2, Prot_mod3, Prot_mod4,Prot_mod5,Prot_mod6,Prot_mod7,Prot_mod8) AICtab(Prot_mod1, Prot_mod2, Prot_mod3, Prot_mod4,Prot_mod5,Prot_mod6,Prot_mod7,Prot_mod8) anova(Prot_mod1, Prot_mod2, Prot_mod3, Prot_mod4,Prot_mod5,Prot_mod6,Prot_mod7,Prot_mod8) ###Get emmeans of Nosema/Chloro treatments Prot_emmeans<-as.data.frame(emmeans(Prot_mod2, ~nosema:chlorothalonil)) ### Sample sizes Prot_ss<-as.data.frame(xtabs(~nosema+chlorothalonil, model.frame(Prot_mod2))) Prot_emmeans<-cbind(Prot_emmeans,Prot_ss[3]) Prot_emmeans #### Infected v not protein ggplot(dat37nos, aes(y = protein_ug_ml, x = infection, fill=chlorothalonil)) + geom_boxplot() ### Test distributions ###Normal qqp(dat37nos$protein_ug_ml, "norm") Nos_prot_na<-dat37nos$protein_ug_ml[!is.na(dat37nos$protein_ug_ml)] plotdist(Nos_prot_na, histo = TRUE, demp = TRUE) descdist(Nos_prot_na, discrete=F, boot=500) ### Normal looks ~ good # Fit linear mixed effects model with all effects. Use drop1 to perform LRT, anova, AIC and AICtab to compare models NosProt_mod1<-lmer(protein_ug_ml~body_size*infection*chlorothalonil+ (1|microcolony:source_colony), data=dat37nos, REML=T) ###Graphing the residuals of model. plot(fitted(NosProt_mod1), residuals(NosProt_mod1), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(NosProt_mod1), residuals(NosProt_mod1))) NosProt_mod1_simres <- simulateResiduals(NosProt_mod1) plot(NosProt_mod1_simres) ### Check random effects NosProt_ran<-ranef(NosProt_mod1,condVar=TRUE) dotplot(NosProt_ran)[["microcolony:source_colony"]] summary(NosProt_mod1) drop1(NosProt_mod1,test="Chisq") NosProt_mod2<-update(NosProt_mod1,.~. -body_size:infection:chlorothalonil) summary(NosProt_mod2) drop1(NosProt_mod2,test="Chisq") NosProt_mod3<-update(NosProt_mod2,.~. -infection:chlorothalonil) summary(NosProt_mod3) drop1(NosProt_mod3,test="Chisq") NosProt_mod4<-update(NosProt_mod3,.~. -body_size:infection) summary(NosProt_mod4) drop1(NosProt_mod4,test="Chisq") NosProt_mod5<-update(NosProt_mod4,.~. -body_size:chlorothalonil) summary(NosProt_mod5) drop1(NosProt_mod5,test="Chisq") NosProt_mod6<-update(NosProt_mod5,.~. -chlorothalonil) summary(NosProt_mod6) drop1(NosProt_mod6,test="Chisq") NosProt_mod7<-update(NosProt_mod6,.~. -infection) summary(NosProt_mod7) drop1(NosProt_mod7,test="Chisq") NosProt_mod8<-update(NosProt_mod7,.~. -body_size) summary(NosProt_mod8) AIC(NosProt_mod1, NosProt_mod2, NosProt_mod3, NosProt_mod4,NosProt_mod5,NosProt_mod6,NosProt_mod7,NosProt_mod8) AICtab(NosProt_mod1, NosProt_mod2, NosProt_mod3, NosProt_mod4,NosProt_mod5,NosProt_mod6,NosProt_mod7,NosProt_mod8) anova(NosProt_mod1, NosProt_mod2, NosProt_mod3, NosProt_mod4,NosProt_mod5,NosProt_mod6,NosProt_mod7,NosProt_mod8) ### Get emmeans of infection/chloro combination NosProt_emmeans<-as.data.frame(emmeans(NosProt_mod2, ~infection:chlorothalonil)) ### Sample sizes NosProt_ss<-as.data.frame(xtabs(~infection+chlorothalonil, model.frame(NosProt_mod2))) NosProt_emmeans<-cbind(NosProt_emmeans,NosProt_ss[3]) NosProt_emmeans ################################################################################################################################## ### Survival - (using Survival data) ### Set contrasts to give hazard ratios relative to reference - return for other cases options(contrasts = c("contr.treatment", "contr.poly")) ### Nosema exposed v not ggplot(dat37surv, aes(x = survival)) + geom_density() + facet_wrap(nosema~chlorothalonil) ggplot(dat37surv, aes(x = survival)) + geom_histogram() + facet_wrap(nosema~chlorothalonil) ggplot(dat37surv, aes(y = survival, x = nosema, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,38) ### All individuals in survival Sur_mod<-coxme(Surv(survival,death_y_n)~body_size*nosema*chlorothalonil+(1|source_colony/microcolony),data=dat37surv) summary(Sur_mod) Anova(Sur_mod, test="LR") Sur_mod2<-update(Sur_mod,.~. -body_size:nosema:chlorothalonil) summary(Sur_mod2) Anova(Sur_mod2, test="LR") Sur_mod3<-update(Sur_mod2,.~. -body_size:chlorothalonil) summary(Sur_mod3) Anova(Sur_mod3, test="LR") Sur_mod4<-update(Sur_mod3,.~. -nosema:chlorothalonil) summary(Sur_mod4) Anova(Sur_mod4, test="LR") Sur_mod5<-update(Sur_mod4,.~. -body_size:nosema) summary(Sur_mod5) Anova(Sur_mod5, test="LR") Sur_mod6<-update(Sur_mod5,.~. -chlorothalonil) summary(Sur_mod6) Anova(Sur_mod6, test="LR") Sur_mod7<-update(Sur_mod6,.~. -body_size) summary(Sur_mod7) Anova(Sur_mod7, test="LR") Sur_mod8<-update(Sur_mod7,.~. -nosema) summary(Sur_mod8) AIC(Sur_mod,Sur_mod2,Sur_mod3,Sur_mod4,Sur_mod5,Sur_mod6,Sur_mod7,Sur_mod8) AICtab(Sur_mod,Sur_mod2,Sur_mod3,Sur_mod4,Sur_mod5,Sur_mod6,Sur_mod7,Sur_mod8) ### Get estimated survival hazards for nosema/chlorothalonil combo to report or plot Surv_emmeans<-emmeans(Sur_mod3,~nosema:chlorothalonil, type="response") plot(Surv_emmeans, comparisons=T) Surv_emmeans<-as.data.frame(Surv_emmeans) Surv_ss<-as.data.frame(xtabs(~nosema+chlorothalonil, model.frame(Sur_mod3))) Surv_ss Surv_emmeans<-cbind(Surv_emmeans,Surv_ss[3]) Surv_emmeans trt<-c("Control","Exposed") chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") Surv_emmeans_plot<-ggplot(Surv_emmeans, aes(x = nosema, y = response, color = chlorothalonil)) + geom_pointrange(aes(ymin = asymp.LCL, ymax = asymp.UCL),size=0.75, position = position_dodge(width=0.6))+ ylab("Estimated survival hazard \u00b1 95% C.I.") + xlab (expression(paste(italic("Nosema bombi")," treatment")))+ ylim(0,1.5) + scale_color_manual(labels = c("Control","Exposed"), values = c("#C0C0C0","#2F4F4F"), name="Chlorothalonil:") + scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") Surv_emmeans_plot ##### Survival within Nosema exposed infected v not dat37surv_nos<-droplevels(subset(dat37surv, nosema==1)) summary(dat37surv_nos) dim(dat37surv_nos) ggplot(dat37surv_nos, aes(x = survival)) + geom_density() + facet_wrap(infection~chlorothalonil) ggplot(dat37surv_nos, aes(x = survival)) + geom_histogram() + facet_wrap(infection~chlorothalonil) ggplot(dat37surv_nos, aes(y = survival, x = infection, fill=chlorothalonil)) + geom_boxplot()+ ylim(0,38) NosSur_mod<-coxme(Surv(survival,death_y_n)~body_size*infection*chlorothalonil+(1|source_colony/microcolony),data=dat37surv_nos) summary(NosSur_mod) Anova(NosSur_mod, test="LR") NosSur_mod2<-update(NosSur_mod,.~. -body_size:infection:chlorothalonil) summary(NosSur_mod2) Anova(NosSur_mod2, test="LR") NosSur_mod3<-update(NosSur_mod2,.~. -body_size:chlorothalonil) summary(NosSur_mod3) Anova(NosSur_mod3, test="LR") NosSur_mod4<-update(NosSur_mod3,.~. -body_size:infection) summary(NosSur_mod4) Anova(NosSur_mod4, test="LR") NosSur_mod5<-update(NosSur_mod4,.~. -infection:chlorothalonil) summary(NosSur_mod5) Anova(NosSur_mod5, test="LR") NosSur_mod6<-update(NosSur_mod5,.~. -chlorothalonil) summary(NosSur_mod6) Anova(NosSur_mod6, test="LR") NosSur_mod7<-update(NosSur_mod6,.~. -infection) summary(NosSur_mod7) Anova(NosSur_mod7, test="LR") NosSur_mod8<-update(NosSur_mod7,.~. -body_size) summary(NosSur_mod8) AIC(NosSur_mod,NosSur_mod2,NosSur_mod3,NosSur_mod4,NosSur_mod5,NosSur_mod6,NosSur_mod7,NosSur_mod8) AICtab(NosSur_mod,NosSur_mod2,NosSur_mod3,NosSur_mod4,NosSur_mod5,NosSur_mod6,NosSur_mod7,NosSur_mod8) ### Get estimated survival hazards for nosema/chlorothalonil combo to report or plot NosSurv_emmeans<-emmeans(NosSur_mod4,~infection:chlorothalonil, type="response") plot(NosSurv_emmeans, comparisons=T) NosSurv_emmeans<-as.data.frame(NosSurv_emmeans) NosSurv_ss<-as.data.frame(xtabs(~infection+chlorothalonil, model.frame(NosSur_mod4))) NosSurv_ss NosSurv_emmeans<-cbind(NosSurv_emmeans,NosSurv_ss[3]) NosSurv_emmeans trt<-c("Uninfected","Infected") chl_col<-c("#C0C0C0","#2F4F4F","#C0C0C0","#2F4F4F") NosSurv_emmeans_plot<-ggplot(NosSurv_emmeans, aes(x = infection, y = response, color = chlorothalonil)) + geom_pointrange(aes(ymin = asymp.LCL, ymax = asymp.UCL),size=0.75, position = position_dodge(width=0.6))+ ylab("Estimated survival hazard \u00b1 95% C.I.") + xlab (expression(paste(italic("Nosema bombi")," status")))+ ylim(0,2) + scale_color_manual(labels = c("Control","Exposed"), values = c("#C0C0C0","#2F4F4F"), name="Chlorothalonil:") + scale_x_discrete(labels = trt) + theme_classic(base_size=20)+ theme(axis.text=element_text(size=15),axis.title=element_text(size=15), axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),title=element_text(size=15))+ theme(legend.position="top") NosSurv_emmeans_plot