# required packages ----------- library(rcompanion) library(bestNormalize) library(car) library(emmeans) library(bestglm) library(LEAP) library(multcomp) # loading data ----------- data_RGR_PierisTW <- read.csv("data.RGR.Taiwan.csv") data_RGR_PierisJP <- read.csv("data.RGR.Japan.csv") data_RGR_HzAm <- read.csv("data.RGR.Hz.America.csv") data_RGR_PierisAm <- read.csv("data.RGR.Pr.America.csv") data_FAW <- read.csv("data.RGR.Sf.America.csv") data_RGR_TniAm <- read.csv("data.RGR.Tn.America.csv") # Fig. 3a Pieris rapae (Taiwan) ------------ data_RGR_PierisTW <- data_RGR_PierisTW[,c(1,5)] colnames(data_RGR_PierisTW) <- c("trt","response") data_RGR_PierisTW <- data_RGR_PierisTW[complete.cases(data_RGR_PierisTW),] # testing model assumptions qqnorm(data_RGR_PierisTW$response) qqline(data_RGR_PierisTW$response) plotNormalHistogram(data_RGR_PierisTW$response) # transformation trdata_RGR_PierisTW <- data_RGR_PierisTW bestNormalize(trdata_RGR_PierisTW$response) orderNorm_obj <- orderNorm(trdata_RGR_PierisTW$response) trdata_RGR_PierisTW$response<-predict(orderNorm_obj) # ANOVA lm_RGR_PierisTW = lm(response ~ trt, data=trdata_RGR_PierisTW) Anova(lm_RGR_PierisTW, type="II") # multiple comparisons lq_RGR_PierisTW = emmeans(lm_RGR_PierisTW,~trt) cld(lq_RGR_PierisTW, by=NULL, alpha=.05, Letters=letters, reverse = FALSE) # Fig. 3b Pieris rapae (Japan) ------------ data_RGR_PierisJP <- data_RGR_PierisJP[,c(1,5)] colnames(data_RGR_PierisJP) <- c("trt","response") data_RGR_PierisJP <- data_RGR_PierisJP[complete.cases(data_RGR_PierisJP),] # testing model assumptions qqnorm(data_RGR_PierisJP$response) qqline(data_RGR_PierisJP$response) plotNormalHistogram(data_RGR_PierisJP$response) # transformation trdata_RGR_PierisJP <- data_RGR_PierisJP bestNormalize(trdata_RGR_PierisJP$response) orderNorm_obj <- orderNorm(trdata_RGR_PierisJP$response) trdata_RGR_PierisJP$response<-predict(orderNorm_obj) # ANOVA lm_RGR_PierisJP = lm(response ~ trt, data=trdata_RGR_PierisJP) Anova(lm_RGR_PierisJP, type="II") # multiple comparisons lq_RGR_PierisJP = emmeans(lm_RGR_PierisJP,~trt) cld(lq_RGR_PierisJP, by=NULL, alpha=.05, Letters=letters, reverse = TRUE) # Fig. 3c Pieris rapae (America) ------------ data_RGR_PierisAm <- data_RGR_PierisAm[,c(1,2,7)] colnames(data_RGR_PierisAm) <- c("wtr","location","response") data_RGR_PierisAm <- data_RGR_PierisAm[complete.cases(data_RGR_PierisAm),] data_RGR_PierisAm$wtr <- as.factor(data_RGR_PierisAm$wtr) data_RGR_PierisAm$location <- as.factor(data_RGR_PierisAm$location) data_RGR_PierisAm$response <- as.numeric(data_RGR_PierisAm$response) # testing model assumptions qqnorm(data_RGR_PierisAm$response) qqline(data_RGR_PierisAm$response) plotNormalHistogram(data_RGR_PierisAm$response) # transformation trdata_RGR_PierisAm <- data_RGR_PierisAm bestNormalize(trdata_RGR_PierisAm$response) orderNorm_obj <- orderNorm(trdata_RGR_PierisAm$response) trdata_RGR_PierisAm$response<-predict(orderNorm_obj) # ANOVA lm_RGR_PierisAm = lm(response ~ wtr+location, data=trdata_RGR_PierisAm) Anova(lm_RGR_PierisAm, type="II") # multiple comparisons lq_RGR_PierisAm = emmeans(lm_RGR_PierisAm,~wtr+location) cld(lq_RGR_PierisAm, by=NULL, alpha=.05, Letters=letters, reverse = TRUE) # Fig. 3d Trichoplusia ni (America) ------------ data_RGR_TniAm <- data_RGR_TniAm[,c(1,2,7)] colnames(data_RGR_TniAm) <- c("wtr","location","response") data_RGR_TniAm <- data_RGR_TniAm[complete.cases(data_RGR_TniAm),] data_RGR_TniAm$wtr <- as.factor(data_RGR_TniAm$wtr) data_RGR_TniAm$location <- as.factor(data_RGR_TniAm$location) data_RGR_TniAm$response <- as.numeric(data_RGR_TniAm$response) #test model assumptions qqnorm(data_RGR_TniAm$response) qqline(data_RGR_TniAm$response) plotNormalHistogram(data_RGR_TniAm$response) trdata_RGR_TniAm <- data_RGR_TniAm # ANOVA lm_RGR_TniAm = lm(response ~ wtr*location, data=trdata_RGR_TniAm) Anova(lm_RGR_TniAm, type="II") # multiple comparisons lq_RGR_TniAm = emmeans(lm_RGR_TniAm,~wtr*location) cld(lq_RGR_TniAm, by=NULL, alpha=.05, Letters=letters, reverse = TRUE) # Fig. 3e Spodoptera frugiperda (America) ------------ data_FAW<-data_FAW[c(1,2,3,4,8)] colnames(data_FAW) <- c("wtr","Location","str_w", "str_i","response") data_FAW<-as.data.frame(data_FAW) data_FAW$wtr <- as.factor(data_FAW$wtr) data_FAW$Location <- as.factor(data_FAW$Location) data_FAW$response <- as.factor(data_FAW$response) # choosing model res.bestglm<-bestglm(data_FAW, family = binomial, IC= "AIC") res.bestglm$BestModels lbr0<-glm(formula = response ~ 1, data =data_FAW, family = binomial(link ="logit")) summary(lbr0) lbr1<-glm(formula = response ~ wtr, data =data_FAW, family = binomial(link ="logit")) summary(lbr1) lbr2<-glm(formula = response ~ wtr+Location, data =data_FAW, family = binomial(link ="logit")) A<-summary(lbr2) lbr3<-glm(formula = response ~ wtr*Location*str_w, data =data_FAW, family = binomial(link ="logit")) summary(lbr3) lbr4<-glm(formula = response ~ wtr+Location+str_w, data =data_FAW, family = binomial(link ="logit")) summary(lbr4) Anova_FAW_lbr<-anova(lbr0,lbr1,lbr2,lbr4, test = "LRT") Anova_FAW_lbr