#R-script for survey analysis survey <- read.delim("2013_survey.txt") View(survey) library(ggplot2) library (lme4) #graph actual trends using transect-level data, note that trends are not very obvious without including multivariate effects survey_trans <- read.delim("2013_survey_transectlevel.txt") View(survey_trans) Length<-ggplot(survey_trans, aes(x = Length, y = Prev2)) + ylim(c(0, 1))+ geom_point(aes(colour = Site))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("WD Prevalence")+xlab("Shoot Length (cm)") +theme(text = element_text(size=16)) +stat_smooth(method = 'lm', aes(colour = 'linear'), se = FALSE) Length Density<-ggplot(survey_trans, aes(x = Density, y = Prev2)) + ylim(c(0, 1))+ geom_point(aes(colour = Site))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("WD Prevalence")+xlab("Shoots per meter squared") +theme(text = element_text(size=16)) +stat_smooth(method = 'lm', aes(colour = 'linear'), se = FALSE) Density Epi<-ggplot(survey_trans, aes(x = Epi, y = Prev2)) + ylim(c(0, 1))+ geom_point(aes(colour = Site))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("WD Prevalence")+xlab("Epiphytes (mg per square cm of shoot)") +theme(text = element_text(size=16)) +stat_smooth(method = 'lm', aes(colour = 'linear'), se = FALSE) Epi require(gridExtra) grid.arrange(Length, Density, Epi) #Each factor alone DisMod.1<-glmer(Diseased~(Length)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.1) DisMod.2<-glmer(Diseased~(Density_decm)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.2) DisMod.3<-glmer(Diseased~(Epi.mg_cm2)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.3) #Combinations of two variables DisMod.4<-glmer(Diseased~(Length)+(Epi.mg_cm2)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.4) DisMod.5<-glmer(Diseased~(Length)+(Density_decm)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.5) DisMod.6<-glmer(Diseased~(Density_decm)+(Epi.mg_cm2)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.6) #Full model, Site as a random intercept DisMod.7<-glmer(Diseased~(Length)+(Density_decm)+(Epi.mg_cm2)+(1|Site), na.action = na.omit, data=survey, family=binomial) summary(DisMod.7) #Compare models, full is best AIC(DisMod.1, DisMod.2, DisMod.3, DisMod.4, DisMod.5, DisMod.6, DisMod.7) #Plot length effect, hold density at 200 shoots/m2 and Epiphytes at 1.8 mg/cm2 #dataset for predictions tmpdat <- survey[, c("Length", "Density_decm", "Epi.mg_cm2", "Site")] jvalues <- with(survey, seq(from = min(Length), to = max(Length), length.out = 100)) #predict results LengthOut <- lapply(jvalues, function(j) { tmpdat$Length<-j tmpdat$Density_decm<-rep.int(2, 1645) tmpdat$Epi.mg_cm2<-rep.int(1.8, 1645) predict(DisMod.7, newdata = tmpdat, type = "response") }) plotLengthOut<- t(sapply(LengthOut, function(x) { c(M = mean(x), quantile(x, c(0.25, 0.75))) })) # create data frame and rename plotLengthOut <- as.data.frame(cbind(plotLengthOut, jvalues)) colnames(plotLengthOut) <- c("PredictedProbability", "Lower", "Upper", "Length") ppLength<-ggplot(plotLengthOut, aes(x = Length, y = PredictedProbability)) + geom_linerange(aes(ymin = Lower, ymax = Upper)) + ylim(c(0, 1))+ geom_line(Density = 2)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("Probability of disease")+xlab("Shoot Length (cm)") +theme(text = element_text(size=16)) ppLength #Plot Density effect, hold Length at 70 mm and Epiphytes at 1.8 mg/cm2 #dataset for predictions tmpdat <- survey[, c("Length", "Density_decm", "Epi.mg_cm2", "Site")] jvalues <- with(survey, seq(from = min(Density_decm), to = max(Density_decm), length.out = 100)) #predict results DensOut <- lapply(jvalues, function(j) { tmpdat$Density_decm<-j tmpdat$Length<-rep.int(70, 1645) tmpdat$Epi.mg_cm2<-rep.int(1.8, 1645) predict(DisMod.7, newdata = tmpdat, type = "response") }) plotDensOut<- t(sapply(DensOut, function(x) { c(M = mean(x), quantile(x, c(0.25, 0.75))) })) # create data frame and rename jvalues<-jvalues*100 plotDensOut <- as.data.frame(cbind(plotDensOut, jvalues)) colnames(plotDensOut) <- c("PredictedProbability", "Lower", "Upper", "Density") ppDens<-ggplot(plotDensOut, aes(x = Density, y = PredictedProbability)) + geom_linerange(aes(ymin = Lower, ymax = Upper)) + ylim(c(0, 1))+ geom_line(Density = 2)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("Probability of disease")+xlab("Shoots per square meter") +theme(text = element_text(size=16)) ppDens #Plot Epiphyte effect, hold Length at 70 mm and Density at 200 shoots/m2 #dataset for predictions tmpdat <- survey[, c("Length", "Density_decm", "Epi.mg_cm2", "Site")] jvalues <- with(survey, seq(from = min(Epi.mg_cm2), to = max(Epi.mg_cm2), length.out = 100)) #predict results EpiOut <- lapply(jvalues, function(j) { tmpdat$Epi.mg_cm2<-j tmpdat$Length<-rep.int(70, 1645) tmpdat$Dens_decm<-rep.int(2.0, 1645) predict(DisMod.7, newdata = tmpdat, type = "response") }) plotEpiOut<- t(sapply(EpiOut, function(x) { c(M = mean(x), quantile(x, c(0.25, 0.75))) })) # create data frame and rename jvalues<-jvalues*100 plotEpiOut <- as.data.frame(cbind(plotEpiOut, jvalues)) colnames(plotEpiOut) <- c("PredictedProbability", "Lower", "Upper", "Epi") ppEpi<-ggplot(plotEpiOut, aes(x = Epi, y = PredictedProbability)) + geom_linerange(aes(ymin = Lower, ymax = Upper)) + ylim(c(0, 1))+ geom_line(Density = 2)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_rect(colour = "black", fill=NA, size=1))+ylab("Probability of disease")+xlab("Epiphyte load (mg per square cm of blade)") +theme(text = element_text(size=16)) ppEpi require(gridExtra) grid.arrange(ppLength, ppDens, ppEpi ) -------------------------------------------------------------------------------- #R-script for Phenol Analysis #Do disease, length, interaction or epiphyte predict phenols? phenol <- read.delim("Phenol_Survey_Final.txt") require(lme4) #Run several models and use AIC to pick the best fit #With disease*length interaction Phen1<-lmer(Phenols~diseased*Length+(1|Site), data=phenol) summary(Phen1) #no interaction Phen2<-lmer(Phenols~diseased+Length+(1|Site) , data=phenol) summary(Phen2) #no interaction, no length Phen3<-lmer(Phenols~diseased+(1|Site), data=phenol) summary(Phen3) #just length as a predictor Phen4<-lmer(Phenols~Length+(1|Site), data=phenol) summary(Phen5) AIC(Phen1, Phen2, Phen3, Phen4) #Phen3 is the best fit. Diagnostics look ok qqnorm(resid(Phen3)) plot(Phen3) plot(resid(Phen3) ~ fitted(Phen4), main="residual plot") abline(h=0) #re-run best-fitting model w/ lmertest to get the p-value require(lmerTest) Phen3.1<-lme(Phenols~diseased, data=phenol, random=~ 1|Site, na.action = na.omit) summary(Phen3.1)