######## data analysis for immigrant viability and fecundity in B. treatae, Linyi Zhang, Dec 2020 ###################### ####################### library(nlme) library(lme4) library(ggplot2) library(emmeans) library(plyr) library(cowplot) library(gridExtra) library(ggpubr) library(lattice) library(glmmADMB) library(glmmTMB) library(DHARMa) library(car) library(patchwork) ################### data.sum<-read.csv("sum.data2017new.csv") data.sum$type<-factor(data.sum$type,levels=c("Resident","Immigrant")) ####### data exploration ################ ##### 1. boxplot of group data ##### str(data.sum) bwplot(HR.OV ~ type| Plant * Treatment, data = data.sum, strip = strip.custom(bg = 'white'), cex = .5, xlab = "Type", ylab = "% host immune response", par.settings = list (box.rectangle = list(col = 1), box.umbrella = list(col = 1), plot.symbol = list(cex = .5, col = 1)), scales = list(x = list(relation = "same"), y = list(relation = "same"))) bwplot(gall.OV ~ type| Plant * Treatment, data = data.sum, strip = strip.custom(bg = 'white'), cex = .5, xlab = "Type", ylab = "% gall formation", par.settings = list (box.rectangle = list(col = 1), box.umbrella = list(col = 1), plot.symbol = list(cex = .5, col = 1)), scales = list(x = list(relation = "same"), y = list(relation = "same"))) ###### 2. Check the distribution of the data ######## histogram( ~ No.OV | cross*Treatment, type = "count", xlab = "No. OV scars", ylab = "Frequency", nint=30,layout=c(2,4), strip.left = strip.custom(bg = 'white'), strip = F, col.line = "black", col = "white", scales = list(x = list(relation = "same"), y = list(relation = "same"), draw = TRUE), data = data.sum) histogram( ~ No.HR | cross*Treatment, type = "count", xlab = "No. HR", ylab = "Frequency", nint=30,layout=c(2,4), strip.left = strip.custom(bg = 'white'), strip = F, col.line = "black", col = "white", scales = list(x = list(relation = "same"), y = list(relation = "same"), draw = TRUE), data = data.sum) histogram( ~ No.galls | cross*Treatment, type = "count", xlab = "No. HR", ylab = "Frequency", nint=30,layout=c(2,4), strip.left = strip.custom(bg = 'white'), strip = F, col.line = "black", col = "white", scales = list(x = list(relation = "same"), y = list(relation = "same"), draw = TRUE), data = data.sum) hist(data.sum$HR.OV) hist(data.sum$gall.OV) ### ####### data analysis ####### ###1. degree of host immune response ### ## analyze field and cage data per tree ## ### analysis field and cage data seperately ### HR.II.fieldm1<-glmmTMB(HR.OV~Plant*type+(1|Site),weights=No.OV,family=binomial, data=data.sum[data.sum$Treatment=="Field",]) simulationOutput1<-simulateResiduals(fittedModel=HR.II.fieldm1) testZeroInflation(simulationOutput1) summary(HR.II.fieldm1) emmeans(HR.II.fieldm1, list(pairwise~Plant*type),adjust="tukey") HR.II.cagem1<-glmmTMB(HR.OV~Plant*type+(1|Site),weights=No.OV,family=binomial, data=data.sum[data.sum$Treatment=="Cage",]) simulationOutput2<-simulateResiduals(fittedModel=HR.II.cagem1) testZeroInflation(simulationOutput2) summary(HR.II.cagem1) emmeans(HR.II.cagem1, list(pairwise~Plant*type),adjust="tukey") #### no differences in the trends, combine field and cage treatment #### HR.II.sm<-glmmTMB(HR.OV~Plant*type+Treatment+(1|Site),weights=No.OV,family=binomial,ziformula=~1, data=data.sum) summary(HR.II.sm) emmeans(HR.II.sm, list(pairwise~Plant*type),adjust="tukey") simulationOutput3<-simulateResiduals(fittedModel=HR.II.sm) testZeroInflation(simulationOutput3) Anova(HR.II.sm,type="III") ############################################################### ################################################################# #### pairwise comparison #### HR.II.cm.pair<-emmeans(HR.II.sm, list(pairwise~Plant*type),adjust="tukey") HR.II.cm.pair<-data.frame(HR.II.cm.pair$`pairwise differences of Plant, type`) HR.II.cm.pair$compare<-"immune.response" ###################################################################### ##### calculate the least square mean and standard error across trees #### HRII.sum<-lsmeans(HR.II.sm,~Plant*type) HRII.sum<-data.frame(summary(HRII.sum,type="response")) HRII.sum<-HRII.sum[order(HRII.sum$Plant,HRII.sum$type),] ## b. plot degree of host immune response ## hostplant<-c('Qv'="Quercus virginiana", 'Qg'="Quercus geminata") p_HROV.text<-data.frame(label=c("",""), Plant=c("Qv","Qg"), y=c(0.7,1.05),x=c(1.5,1.5)) HR.OV.II.plot<-ggplot(data.sum)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = HR.OV,fill = type),alpha = 0.80,outlier.shape = NA,show.legend = FALSE) + geom_point(aes(x=type,y=HR.OV,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,1.1),breaks=c(0,0.25,0.5,0.75,1),labels=c("0" = " 0", "0.25" = "25", "0.5" = "50" ,"0.75" = "75","1" = "100"))+ ylab("% host immune response")+ geom_text(data=p_HROV.text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none") HR.OV.II.plot ####### HR.OV.II plot without x-axis label # ########## HR.OV.II.plot1<-ggplot(data.sum)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = HR.OV,fill = type),alpha = 0.80,outlier.shape = NA,show.legend = FALSE) + geom_point(aes(x=type,y=HR.OV,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,1.1),breaks=c(0,0.25,0.5,0.75,1),labels=c("0" = "0", "0.25" = "25", "0.5" = "50" ,"0.75" = "75","1" = "100"))+ ylab("% host immune response")+ geom_text(data=p_HROV.text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), axis.text.x=element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust = 0.5,size=30))+ ggtitle("Host environment") HR.OV.II.plot1 ####### 2. proportion of gall formation ############# ###### a. analyze field and cage data seperately ########## ### field data analysis ### gall.II.fieldm<-glmmTMB(gall.OV~Plant*type+(1|Site),weights=No.OV,family=binomial,data=data.sum[data.sum$Treatment=="Field",]) summary(gall.II.fieldm) emmeans(gall.II.fieldm, list(pairwise~Plant*type),adjust="tukey") simulationOutput1<-simulateResiduals(fittedModel=gall.II.fieldm) testZeroInflation(simulationOutput1) gall.II.fieldm<-glmmTMB(gall.OV~Plant*type+(1|Site),weights=No.OV,family=binomial,ziformula=~1, data=data.sum[data.sum$Treatment=="Field",]) summary(gall.II.fieldm) emmeans(gall.II.fieldm, list(pairwise~Plant*type),adjust="tukey") simulationOutput1<-simulateResiduals(fittedModel=gall.II.fieldm) testZeroInflation(simulationOutput1) gall.IIfield.pair<-emmeans(gall.II.fieldm, list(pairwise~Plant*type),adjust="tukey") gall.IIfield.pair<-data.frame(gall.IIfield.pair$`pairwise differences of Plant, type`) gall.IIfield.pair$compare<-"gall.formation" ## field plot % gall formation ## hostplant<-c('Qv'="Quercus virginiana", 'Qg'="Quercus geminata") p_gallfield.text<-data.frame(label=c("","***"), Plant=c("Qv","Qg"), y=c(0.7,1.05),x=c(1.5,1.5)) gall.OV.field.plot<-ggplot(data.sum[data.sum$Treatment=="Field",])+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = gall.OV,fill = type),alpha = 0.80,outlier.shape = NA,show.legend = FALSE) + geom_point(aes(x=type,y=gall.OV,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,1.1),breaks=c(0,0.25,0.5,0.75,1),labels=c("0" = " 0", "0.25" = "25", "0.5" = "50" ,"0.75" = "75","1" = "100"))+ ylab("% gall formation")+ geom_text(data=p_gallfield.text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust=0.5,size=30))+ ggtitle("Host environment (field)") gall.OV.field.plot ###### cage % gall formation analysis #### gall.II.cagem<-glmmTMB(gall.OV~Plant*type+(1|Site),weights=No.OV,family=binomial,data=data.sum[data.sum$Treatment=="Cage",]) summary(gall.II.cagem) emmeans(gall.II.cagem, list(pairwise~Plant*type),adjust="tukey") simulationOutput2<-simulateResiduals(fittedModel=gall.II.cagem) testZeroInflation(simulationOutput2) gall.II.cagem<-glmmTMB(gall.OV~Plant*type+(1|Site),weights=No.OV,family=binomial,ziformula=~1, data=data.sum[data.sum$Treatment=="Cage",]) summary(gall.II.cagem) emmeans(gall.II.cagem, list(pairwise~Plant*type),adjust="tukey") simulationOutput2<-simulateResiduals(fittedModel=gall.II.cagem) testZeroInflation(simulationOutput2) gall.IIcage.pair<-emmeans(gall.II.cagem, list(pairwise~Plant*type),adjust="tukey") gall.IIcage.pair<-data.frame(gall.IIcage.pair$`pairwise differences of Plant, type`) gall.IIcage.pair$compare<-"gall.formation" ## cage plot % gall formation ## hostplant<-c('Qv'="Quercus virginiana", 'Qg'="Quercus geminata") p_gallcage.text<-data.frame(label=c("","**"), Plant=c("Qv","Qg"), y=c(0.7,1.05),x=c(1.5,1.5)) gall.OV.cage.plot<-ggplot(data.sum[data.sum$Treatment=="Cage",])+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = gall.OV,fill = type),alpha = 0.80,outlier.shape = NA,show.legend = FALSE) + geom_point(aes(x=type,y=gall.OV,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,1.1),breaks=c(0,0.25,0.5,0.75,1),labels=c("0" = " 0", "0.25" = "25", "0.5" = "50" ,"0.75" = "75","1" = "100"))+ ylab("")+ geom_text(data=p_gallcage.text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), axis.text.y = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust=0.5,size=30))+ ggtitle("Host environment (cage)") gall.OV.cage.plot gall.OV.field.plot|gall.OV.cage.plot ### b. combine data from field and cage to analynize ##### gall.II.sm<-glmmTMB(gall.OV~Plant*type+Treatment+(1|Site),weights=No.OV,family=binomial,ziformula=~1, data=data.sum) summary(gall.II.sm) Anova(gall.II.sm) gallII.Output0<-simulateResiduals(fittedModel=gall.II.sm) plot(gallII.Output0) testZeroInflation(gallII.Output0) emmeans(gall.II.sm, list(pairwise~Plant*type),adjust="tukey") gall.II.cm.pair<-emmeans(gall.II.sm, list(pairwise~Plant*type),adjust="tukey") gall.II.cm.pair<-data.frame(gall.II.cm.pair$`pairwise differences of Plant, type`) gall.II.cm.pair$compare<-"gall.formation" ##### calculate the least square mean and standard error across trees #### gallII.sum<-lsmeans(gall.II.sm,~Plant*type) gallII.sum<-data.frame(summary(gallII.sum,type="response")) gallII.sum<-gallII.sum[order(gallII.sum$Plant,gallII.sum$type),] gallII.sum$type<-factor(gallII.sum$type,levels=c("Resident","Immigrant")) ##################################################################################### ## b. plot proportion of gall formation ## p_gallOV.text<-data.frame(label=c("","***"), Plant=c("Qv","Qg"), y=c(0.7,1.05),x=c(1.5,1.5)) gall.II.plot<-ggplot(data.sum)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = gall.OV,fill = type),alpha = 0.80,outlier.shape = NA,show.legend = FALSE) + geom_point(aes(x=type,y=gall.OV,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,1.1),breaks=c(0,0.25,0.5,0.75,1),labels=c("0" = "0", "0.25" = "25", "0.5" = "50" ,"0.75" = "75","1" = "100"))+ ylab("% gall formation")+ geom_text(data=p_gallOV.text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=26), strip.text.x=element_blank(), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none") gall.II.plot plot_grid(HR.OV.II.plot1,gall.II.plot,nrow=2,labels=c("A","B"),align="v",label_size=38) #################################################################### ################################################################## #### calculation of the strength of immigrant inviability ### ######### do the calculation ####### sum.II<-ddply(data.sum,c("Site","Host","Plant","type"), summarize,No.OV=sum(No.OV),No.HR=sum(No.HR), No.galls=sum(No.galls), No.2mmgall=sum(No.2mmgall),No.2.62mmgall=sum(No.2.62mmgall),No.5.82mmgall=sum(No.5.82mmgall)) sum.II$gall.OV<-sum.II$No.galls/sum.II$No.OV sum.II.Qv<-sum.II[sum.II$Plant=="Qv",] sum.II.Qv.img<-sum.II.Qv[sum.II.Qv$type=="Immigrant",] sum.II.Qv.res<-sum.II.Qv[sum.II.Qv$type=="Resident",] n<-length(unique(sum.II.Qv.img$Site))*length(unique(sum.II.Qv.res$Site)) ### calculate RI between different host on host Qv n<-1 II.RI.Qv<-c() pair.Qv<-c() for(j in 1:3){ for (i in 1:3) { II.RI.Qv[n]<-1-2*sum.II.Qv.img$gall.OV[j]/(sum.II.Qv.res$gall.OV[i]+sum.II.Qv.img$gall.OV[j]) pair.Qv[n]<-paste(sum.II.Qv.res$Site[i],sum.II.Qv.img$Site[j],sep=" x ") n<-n+1 } } II.RI.Qv<-data.frame(Plant="Qv",pair=pair.Qv,II.RI=II.RI.Qv) ### calculate RI between same host on host Qv n<-1 II.shRI.Qv<-c() sh.pair.Qv<-c() for(j in 1:3){ for (i in 1:3) { II.shRI.Qv[n]<-1-2*sum.II.Qv.res$gall.OV[j]/(sum.II.Qv.res$gall.OV[i]+sum.II.Qv.res$gall.OV[j]) sh.pair.Qv[n]<-paste(sum.II.Qv.res$Site[i],sum.II.Qv.res$Site[j],sep=" x ") n<-n+1 } } II.shRI.Qv<-data.frame(Plant="Qv",pair=sh.pair.Qv,II.RI=II.shRI.Qv) ### calculate RI between different host on host Qg sum.II.Qg<-sum.II[sum.II$Plant=="Qg",] sum.II.Qg.img<-sum.II.Qg[sum.II.Qg$type=="Immigrant",] sum.II.Qg.res<-sum.II.Qg[sum.II.Qg$type=="Resident",] n<-length(unique(sum.II.Qg.img$Site))*length(unique(sum.II.Qg.res$Site)) n<-1 II.RI.Qg<-c() pair.Qg<-c() for(j in 1:length(unique(sum.II.Qg.img$Site))){ for (i in 1:length(unique(sum.II.Qg.res$Site))) { II.RI.Qg[n]<-1-2*sum.II.Qg.img$gall.OV[j]/(sum.II.Qg.res$gall.OV[i]+sum.II.Qg.img$gall.OV[j]) pair.Qg[n]<-paste(sum.II.Qg.res$Site[i],sum.II.Qg.img$Site[j],sep=" x ") n<-n+1 } } II.RI.Qg<-data.frame(Plant="Qg",pair=pair.Qg,II.RI=II.RI.Qg) ### calculate RI between same host on host Qg n<-1 II.shRI.Qg<-c() sh.pair.Qg<-c() for(j in 1:length(unique(sum.II.Qg.res$Site))){ for (i in 1:length(unique(sum.II.Qg.res$Site))) { II.shRI.Qg[n]<-1-2*sum.II.Qg.res$gall.OV[j]/(sum.II.Qg.res$gall.OV[i]+sum.II.Qg.res$gall.OV[j]) sh.pair.Qg[n]<-paste(sum.II.Qg.res$Site[i],sum.II.Qg.res$Site[j],sep=" x ") n<-n+1 } } II.shRI.Qg<-data.frame(Plant="Qg",pair=sh.pair.Qg,II.RI=II.shRI.Qg) II.RI<-rbind(II.RI.Qv,II.RI.Qg) II.RI$Wasp[II.RI$Plant=="Qv"]<-"Wasp from host Qg to Qv" II.RI$Wasp[II.RI$Plant=="Qg"]<-"Wasp from host Qv to Qg" II.RI$Plant<-factor(II.RI$Plant,levels=c("Qg","Qv")) II.RI.sum<-ddply(II.RI,c("Plant","Wasp"),summarize,II.RI.mean=mean(II.RI), II.RI.se=sd(II.RI)/sqrt(length(II.RI))) II.shRI<-rbind(II.shRI.Qv,II.shRI.Qg) II.shRI<-II.shRI[II.shRI$II.RI>0,] ### plot the strength of immigrant inviability ### II.RI.plot<-ggplot()+ geom_bar(data=II.RI.sum,aes(x=Plant,y=II.RI.mean,fill=Plant),colour="black",stat="identity", width = 0.5)+ geom_errorbar(data=II.RI.sum,aes(x=Plant,ymin=II.RI.mean-II.RI.se, ymax=II.RI.mean+II.RI.se),width=0.2)+ geom_point(data=II.RI,aes(x=Plant,y=II.RI), size = 2, alpha=0.8, show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(-0.65,0.85),breaks=c(-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8))+ scale_x_discrete(labels=c("Qv" = "Q. virginiana", "Qg" = "Q. geminata"))+ geom_hline(yintercept=0,linetype="solid",color="black")+ ylab("Reprodutive isolation")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), axis.text.x=element_text(face="italic",size=30), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(size=34,hjust = 0.5))+ ggtitle("Immigrant viability") II.RI.plot ######################################################################### #### immigrant fecundity #### tibia<-read.csv("tibia length.csv") tibia$Host<-"Qv" tibia$Host[tibia$site %in% c("ABS","DCK","LL")]<-"Qg" tibia1<-tibia[tibia$treatment %in% c("t","c"),] tibia1$type<-"Resident" tibia1$type[tibia1$treatment=="t"]<-"Immigrant" tibia1$type<-as.factor(tibia1$type) tibia1$type<-factor(tibia1$type,levels=c("Resident","Immigrant")) tibia1$Plant<-tibia1$Host tibia1$Plant[tibia1$type=="Immigrant" & tibia1$Host=="Qv"]<-"Qg" tibia1$Plant[tibia1$type=="Immigrant" & tibia1$Host=="Qg"]<-"Qv" tibia1$Plant<-as.factor(tibia1$Plant) tibia1$Plant<-factor(tibia1$Plant,levels=c("Qg","Qv")) tibia1$cross<-paste(tibia1$Plant,tibia1$type,sep="_") tibia1$No.egg<--225.4+423.9*tibia1$tibia.length.mm ### average number of eggs in each treatment ##### sum.egg<-ddply(tibia1,c("Plant","type"),summarize,No.egg=mean(No.egg)) ########################################## ##### read the tube number information #### tube.number<-read.csv("Raw data from GRH.csv") head(tube.number) tube.number<-tube.number[,-13] head(tibia1) tube1<-tube.number[,c("Tube.NO","Plant","Tree.No","Site","Treatment")] colnames(tube1)[1]<-"tube" tibia1[is.na(tibia1$tube),"tube"]<-344 ##### combine the tibia measurement with tube number information ##### tibia.new<-merge(tibia1[,c("tube","treatment","tibia.length.mm","in.vitro.ID","Enter.by","No.egg")],tube1, by="tube",all.x = TRUE) tibia.new<-tibia.new[!is.na(tibia.new$Treatment),] colnames(tibia.new)[2]<-"code" tibia.new$count<-1 tibia.new$Host<-"Qv" tibia.new$Host[tibia.new$Site %in% c("LL","DCK","ABS")]<-"Qg" tibia.new$type<-"Resident" tibia.new$type[tibia.new$Plant!=tibia.new$Host]<-"Immigrant" tibia.new<-tibia.new[tibia.new$Treatment %in% c("Cage","Field"),] tibia.new$cross<-paste(tibia.new$Plant,tibia.new$type,sep="_") tibia.new$Host<-factor(tibia.new$Host,levels=c("Qg","Qv")) tibia.new$type<-factor(tibia.new$type,levels=c("Resident","Immigrant")) str(tibia.new) individual.sum<-ddply(tibia.new,c("Treatment","Plant","type","Site"),summarize,number=sum(count)) #### don't run write.csv(individual.sum,"Table S3. Number of individuals per site per treatment.csv") #### ### data analysis ### ## data exploration ## bwplot(No.egg ~ type| Plant * Treatment, data = tibia.new, strip = strip.custom(bg = 'white'), cex = .5, xlab = "Type", ylab = "No.eggs", par.settings = list (box.rectangle = list(col = 1), box.umbrella = list(col = 1), plot.symbol = list(cex = .5, col = 1)), scales = list(x = list(relation = "same"), y = list(relation = "same"))) histogram( ~ No.egg | cross*Treatment, type = "count", xlab = "No.eggs", ylab = "Frequency", nint=30,layout=c(2,4), strip.left = strip.custom(bg = 'white'), strip = F, col.line = "black", col = "white", scales = list(x = list(relation = "same"), y = list(relation = "same"), draw = TRUE), data = tibia.new) ##### library(nlme) ##### egg.fieldm<-lme(No.egg~Plant*type,random=~1|Site,data=tibia.new[tibia.new$Treatment=="Field",]) summary(egg.fieldm) emmeans(egg.fieldm, list(pairwise~Plant*type),adjust="tukey") eggfield.sum<-lsmeans(egg.fieldm,~Plant*type) eggfield.sum<-data.frame(summary(eggfield.sum,type="response")) eggfield.sum<-eggfield.sum[order(eggfield.sum$Plant,eggfield.sum$type),] Egg.field.plot<-ggplot(data=eggfield.sum)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_linerange(aes( x=type,ymin= lower.CL,ymax=upper.CL,fill = type),show.legend = FALSE) + geom_point(aes(x=type,y=lsmean,fill = type), size = 5, shape = 21, show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(100,300),breaks=seq(100,300,50))+ ylab("Predicted number of eggs")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 26,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust=0.5,size=30))+ ggtitle("Host environment (field)") Egg.field.plot egg.cagem<-lme(No.egg~Plant*type,random=~1|Site,data=tibia.new[tibia.new$Treatment=="Cage",]) summary(egg.cagem) emmeans(egg.cagem, list(pairwise~Plant*type),adjust="tukey") eggcage.sum<-lsmeans(egg.cagem,~Plant*type) eggcage.sum<-data.frame(summary(eggcage.sum,type="response")) eggcagesum<-eggcage.sum[order(eggcage.sum$Plant,eggcage.sum$type),] Egg.cage.plot<-ggplot(data=eggcage.sum)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_linerange(aes( x=type,ymin= lower.CL,ymax=upper.CL,fill = type),show.legend = FALSE) + geom_point(aes(x=type,y=lsmean,fill = type), size = 5, shape = 21, show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(100,300),breaks=seq(100,300,50))+ ylab("")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 26,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), axis.text.y = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust=0.5,size=30))+ ggtitle("Host environment (cage)") Egg.cage.plot Egg.field.plot|Egg.cage.plot ############################################################### egg.m1<-lme(No.egg~Plant*type+Treatment,random=~1|Site,data=tibia.new) summary(egg.m1) Anova(egg.m1,type="III") emmeans(egg.m1, list(pairwise~Plant*type),adjust="tukey") egg.field1<-lme(No.egg~Plant*type,random=~1|Site,data=tibia.new[tibia.new$Treatment=="Field",]) summary(egg.field1) emmeans(egg.field1, list(pairwise~Plant*type),adjust="tukey") egg.cage1<-lme(No.egg~Plant*type,random=~1|Site,data=tibia.new[tibia.new$Treatment=="Cage",]) summary(egg.cage1) emmeans(egg.cage1, list(pairwise~Plant*type),adjust="tukey") tibia.m1<-lme(tibia.length.mm~Plant*type+Treatment,random=~1|Site,data=tibia.new) summary(tibia.m1) emmeans(tibia.m1, list(pairwise~Plant*type),adjust="tukey") egg.m1.pair<-emmeans(egg.m1, list(pairwise~Plant*type),adjust="tukey") egg.m1.pair<-data.frame(egg.m1.pair$`pairwise differences of Plant, type`) egg.m1.pair$compare<-"egg.number" ##### get all the pairwise comparison result out ###### pairwise.result<-rbind(HR.II.cm.pair,gall.II.cm.pair,egg.m1.pair) write.csv(pairwise.result,"pairwise.result.csv") ##### calculate the least square mean and standard error of egg numbers and body size #### eggIF.sum<-lsmeans(egg.m1,~Plant*type) eggIF.sum<-data.frame(summary(eggIF.sum,type="response")) eggIF.sum<-eggIF.sum[order(eggIF.sum$Plant,eggIF.sum$type),] tibiaIF.sum<-lsmeans(tibia.m1,~Plant*type) tibiaIF.sum<-data.frame(summary(tibiaIF.sum,type="response")) tibiaIF.sum<-tibiaIF.sum[order(tibiaIF.sum$Plant,tibiaIF.sum$type),] #### plot the number of eggs #### egg.p_text<-data.frame(label=c("*",""), Plant=c("Qv","Qg"), y=c(380,1.05),x=c(1.5,1.5)) Egg.NO.plot<-ggplot(tibia1)+ facet_grid(~Plant,labeller=as_labeller(hostplant))+ geom_boxplot(aes( x=type,y = No.egg,fill = type),alpha = 0.80,show.legend = FALSE) + geom_point(aes(x=type,y=No.egg,fill = type), size = 5, shape = 21, position = position_jitterdodge(),show.legend = FALSE) + scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(0,400),breaks=seq(0,400,100))+ ylab("Predicted number of eggs")+ geom_text(data= egg.p_text,aes(x=x,y=y,label=label),size=14)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 26,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(hjust=0.5,size=30))+ ggtitle("Host environment") Egg.NO.plot ###### calculate the strength of immigrant fecundity ######## egg.sum<-ddply(tibia1,c("site","Plant","type"),summarize, No.egg=mean(No.egg)) egg.sum.Qv<-egg.sum[egg.sum$Plant=="Qv",] egg.sum.Qg<-egg.sum[egg.sum$Plant=="Qg",] egg.sum.Qv.img<-egg.sum.Qv[egg.sum.Qv$type=='Immigrant',] egg.sum.Qv.res<-egg.sum.Qv[egg.sum.Qv$type=='Resident',] egg.sum.Qg.img<-egg.sum.Qg[egg.sum.Qg$type=='Immigrant',] egg.sum.Qg.res<-egg.sum.Qg[egg.sum.Qg$type=='Resident',] ### calculate RI between different host comparison on host Qv ### n<-1 IF.RI.Qv<-c() pair.Qv<-c() for(j in 1:3){ for (i in 1:3) { IF.RI.Qv[n]<-1-2*egg.sum.Qv.img$No.egg[j]/(egg.sum.Qv.res$No.egg[i]+egg.sum.Qv.img$No.egg[j]) pair.Qv[n]<-paste(egg.sum.Qv.res$site[i],egg.sum.Qv.img$site[j],sep=" x ") n<-n+1 } } IF.RI.Qv<-data.frame(Plant="Qv",pair=pair.Qv,IF.RI=IF.RI.Qv) ### calculate RI between same host comparison on host Qv ### n<-1 IF.shRI.Qv<-c() sh.pair.Qv<-c() for(j in 1:3){ for (i in 1:3) { IF.shRI.Qv[n]<-1-2*egg.sum.Qv.res$No.egg[j]/(egg.sum.Qv.res$No.egg[i]+egg.sum.Qv.res$No.egg[j]) sh.pair.Qv[n]<-paste(egg.sum.Qv.res$site[i],egg.sum.Qv.res$site[j],sep=" x ") n<-n+1 } } IF.shRI.Qv<-data.frame(Plant="Qv",pair=sh.pair.Qv,RI=IF.shRI.Qv) ### calculate RI between different host comparison on host Qg ### n<-1 IF.RI.Qg<-c() pair.Qg<-c() for(j in 1:3){ for (i in 1:3) { IF.RI.Qg[n]<-1-2*egg.sum.Qg.img$No.egg[j]/(egg.sum.Qg.res$No.egg[i]+egg.sum.Qg.img$No.egg[j]) pair.Qg[n]<-paste(egg.sum.Qg.res$site[i],egg.sum.Qg.img$site[j],sep=" x ") n<-n+1 } } IF.RI.Qg<-data.frame(Plant="Qg",pair=pair.Qg,IF.RI=IF.RI.Qg) ### calculate RI between same host comparison on host Qg ### n<-1 IF.shRI.Qg<-c() sh.pair.Qg<-c() for(j in 1:3){ for (i in 1:3) { IF.shRI.Qg[n]<-1-2*egg.sum.Qg.res$No.egg[j]/(egg.sum.Qg.res$No.egg[i]+egg.sum.Qg.res$No.egg[j]) sh.pair.Qg[n]<-paste(egg.sum.Qg.res$site[i],egg.sum.Qg.res$site[j],sep=" x ") n<-n+1 } } IF.shRI.Qg<-data.frame(Plant="Qg",pair=sh.pair.Qg,RI=IF.shRI.Qg) IF.RI<-rbind(IF.RI.Qv,IF.RI.Qg) IF.RI$Wasp[IF.RI$Plant=="Qv"]<-"Wasp from host Qg to Qv" IF.RI$Wasp[IF.RI$Plant=="Qg"]<-"Wasp from host Qv to Qg" IF.RI.sum<-ddply(IF.RI,c("Plant","Wasp"),summarize,IF.RI.mean=mean(IF.RI), IF.RI.se=sd(IF.RI)/sqrt(length(IF.RI))) IF.RI.sum$Plant<-factor(IF.RI.sum$Plant,levels=c("Qg","Qv")) IF.shRI<-rbind(IF.shRI.Qv,IF.shRI.Qg) IF.shRI<-IF.shRI[IF.shRI$RI>0,] ###### plot the strength of immigrant fecundity #### IF.RI.plot<-ggplot()+ geom_bar(data=IF.RI.sum,aes(x=Plant,y=IF.RI.mean,fill=Plant),colour="black",stat="identity", width = 0.5)+ geom_errorbar(data=IF.RI.sum,aes(x=Plant,ymin=IF.RI.mean-IF.RI.se, ymax=IF.RI.mean+IF.RI.se,fill=Plant),width=0.2)+ geom_point(data=IF.RI,aes(x=Plant,y=IF.RI), size = 2,alpha=0.8,show.legend = FALSE)+ scale_fill_manual(values=c("grey","white"))+ scale_y_continuous(limits=c(-0.65,0.85),breaks=c(-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8))+ scale_x_discrete(labels=c("Qv" = "Q. virginiana", "Qg" = "Q. geminata"))+ geom_hline(yintercept=0,linetype="solid",color="black")+ ylab("")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background=element_blank(), axis.line=element_line(colour="black"), text = element_text(size = 28,colour = "black"), axis.text=element_text(colour="black"), strip.text = element_text(face = "italic",size=28), axis.text.x=element_text(face="italic",size=30), axis.title.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), legend.position = "none", plot.title = element_text(size=34,hjust = 0.5))+ ggtitle("Immigrant fecundity") IF.RI.plot ############################################################ #### plot the immigrant inviability and immigrant fecundity together #### plot_grid(II.RI.plot,IF.RI.plot,nrow=1,labels=c("A","B"),align="v",label_size=38)