library(car) library(ggplot2) library(multcomp) library(fBasics) library(e1071) library(lme4) library(lmerTest) ################################################################################################ ################################################################################################ ################################################################################################ #################### TRICHOCARPA DATA ######################## ################################################################################################ ################################################################################################ ################################################################################################ tritrait<-read.csv("Henning_PeerJ_trait_data.csv") tritrait names(tritrait) tritrait$Treatment<-factor(tritrait$Treatment, levels=c("C", "GM30", "GM41", "BT03")) boxplot(tritrait$Initial.Root.SA) # not normal normalTest(tritrait$Initial.Root.SA, method=c( "da" )) # only need to look at omnibus value. logINrootSA<-(log(tritrait$Initial.Root.SA+1)) normalTest(logINrootSA, method=c( "da" )) boxdelt_LeafSAplot(logINrootSA) # logINrootSA NORMAL boxplot(tritrait$Init_Leaf_SA) # normal normalTest(tritrait$Init_Leaf_SA, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$shoot) # normal normalTest(tritrait$shoot, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$root.biomass) # normal normalTest(tritrait$root.biomass, method=c( "da" )) # only need to look at omnibus value. logrootmass<-(log(tritrait$root.biomass+1)) normalTest(logrootmass, method=c( "da" )) boxplot(tritrait$Whole.plant) # normal normalTest(tritrait$Whole.plant, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$percentroot) # not normal normalTest(tritrait$percentroot, method=c( "da" )) logperroot<-(log(tritrait$percentroot+1)) boxplot(logperroot) normalTest(zperroot, method=c( "da" )) ### add logperroot to data table boxplot(tritrait$percentshoot) # not normal normalTest(tritrait$percentshoot, method=c( "da" )) logpershoot<-(log(tritrait$percentshoot+1)) boxplot(logpershoot) normalTest(logpershoot, method=c( "da" )) ### add logpershoot to data table boxplot(tritrait$leaf.weight) # normal normalTest(tritrait$leaf.weight, method=c( "da" )) boxplot(tritrait$percent.leaf) # not normal normalTest(tritrait$percent.leaf, method=c( "da" )) logperleaf<-(log(tritrait$percent.leaf+1)) boxplot(logperleaf) normalTest(logperleaf, method=c( "da" )) ### add logpershoot to data table boxplot(tritrait$rootleaf) # not normal normalTest(tritrait$rootleaf, method=c( "da" )) logrootleaf<-(log(tritrait$rootleaf+1)) boxplot(logrootleaf) normalTest(logrootleaf, method=c( "da" )) ### add logrootleaf to data table boxplot(tritrait$rootshoot) # normal normalTest(tritrait$rootshoot, method=c( "da" )) # only need to look at omnibus value. logrs<-((tritrait$rootshoot)^(1/5)) normalTest(logrs, method=c( "da" )) boxplot(logrs) boxplot(tritrait$SPAD) # normal normalTest(tritrait$SPAD, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$CH_leafN) # normal normalTest(tritrait$CH_leafN, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$delt_RootSA) # normal normalTest(tritrait$delt_RootSA, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$delt_RootL) # normal normalTest(tritrait$delt_RootL, method=c( "da" )) # only need to look at omnibus value. boxplot(tritrait$Delt_stem) # normal normalTest(tritrait$Delt_stem, method=c( "da" )) # only need to look at omnibus value. sqdeltstem<-((tritrait$Delt_stem)^(1/2)) normalTest(sqdeltstem, method=c( "da" )) boxplot(sqdeltstem) boxplot(tritrait$delt_LeafSA) # normal normalTest(tritrait$delt_LeafSA, method=c( "da" )) # only need to look at omnibus value. logdeltleaf<-(log(tritrait$delt_LeafSA+1)) normalTest(logdeltleaf, method=c( "da" )) boxplot(logdeltleaf) boxplot(tritrait$SRL) # normal normalTest(tritrait$SRL, method=c( "da" )) # only need to look at omnibus value. logSRL<-(log(tritrait$SRL+1)) normalTest(logSRL, method=c( "da" )) boxplot(tritrait$SLA) # normal normalTest(tritrait$SLA, method=c( "da" )) # only need to look at omnibus value. logSLA<-(log(tritrait$SLA+1)) normalTest(logSLA, method=c( "da" )) boxplot(logSLA) tritrait["logINrootSA"]<-NA tritrait$logINrootSA<-logINrootSA tritrait["logrootmass"]<-NA tritrait$logrootmass<-logrootmass tritrait["logrs"]<-NA tritrait$logrs<-logrs tritrait["logSLA"]<-NA tritrait$logSLA<-logSLA tritrait["logSRL"]<-NA tritrait$logSRL<-logSRL tritrait["sqdeltstem"]<-NA tritrait$sqdeltstem<-sqdeltstem tritrait["logdeltleaf"]<-NA tritrait$logdeltleaf<-logdeltleaf tritrait["logperroot"]<-NA tritrait$logperroot<-logperroot tritrait["logpershoot"]<-NA tritrait$logpershoot<-logpershoot tritrait["logperleaf"]<-NA tritrait$logperleaf<-logperleaf tritrait["logrootleaf"]<-NA tritrait$logrootleaf<-logrootleaf tritrait$Treatment<-factor(tritrait$Treatment, levels=c("C", "GM30", "GM41", "BT03")) names(tritrait) tritrait SRA<-(tritrait$Fin_Sum.of.Root.SA/tritrait$root.biomass) SRA tritrait["SRA"]<-NA tritrait$SRA<-SRA normalTest(SRA, method=c( "da" )) logSRA<-(SRA^(1/4)) normalTest(logSRA, method=c( "da" )) boxplot(logSRA) tritrait["logSRA"]<-NA tritrait$logSRA<-logSRA #### SRA mod<-lmer(SRA~Treatment+(1|EP.run), data=tritrait) # full model summary(mod) mod1<-lmer(logSRA~Treatment +(1|EP.run), REML=FALSE, data=tritrait) # full model mod2<-lmer(logSRA~1 +(1|EP.run), REML=FALSE,data=tritrait) # full model anova(mod1,mod2) summary(mod2) summary(mod1) mod<-lm(SRL~Treatment, data=tritrait) # full model summary(mod) # Shoot biomass - shoot ########################################################################################################################### tritmod<-lm(tritrait$shoot~tritrait$Treatment+tritrait$Init_Leaf_SA) Anova(tritmod) ##################### shoot biomass mixed model ############## mod<-lmer(shoot~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) newmod<-lmer(shoot~Treatment + Init_Leaf_SA +(1|EP.run), data=tritrait) summary(newmod) # testing significance of adding treatment of my mod mod1<-lmer(shoot~Treatment +Init_Leaf_SA + (1|EP.run), data=tritrait, REML=FALSE) summary(mod1) modb<-lmer(shoot~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(shoot~1 +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(shoot~Init_Leaf_SA+(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) anova(moda, modb) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # plot(difflsmeans(mod)) plot(difflsmeans(newmod)) par(mar=c(7,6,6,3)) boxplot(shoot~Treatment, ylim=c(0,200),ylab="Shoot biomass (mg)", lwd=5, xaxt="n",col=c("darkgreen"), data=tritrait, cex.lab=3.5, cex.axis=2.5) axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp."),cex.axis=2.1) text(1, 130, "B",cex=3) text(2, 181, "A",cex=3) text(3, 164, "AB",cex=3) text(4, 197, "AB",cex=3) ########################### root biomass -########################################################################################################################### ######### ROOT BIOMASS MIXED MODEL ###### mod<-lmer(logrootmass~Treatment + (1|EP.run), data=tritrait) # full model summary(mod) print(mod) newmod<-lmer(logrootmass~Treatment +Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(newmod) # testing significance of adding treatment of my mod mod1<-lmer(logrootmass~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logrootmass~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logrootmass~Initial.Root.SA+(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logrootmass~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) summary(mod1) library(lmerTest) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # boxplot(root.biomass~Treatment, data=tritrait, ylim=c(0,90), ylab="root biomass (mg)", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") difflsmeans(newmod) par(mar=c(7,6,6,3)) boxplot(root.biomass~Treatment, data=tritrait,xaxt="n", ylim=c(0,90), lwd=5, ylab="Root biomass (mg)", col=c("darkgreen"), cex.axis=2.5, cex.lab=3.5) axis(1, at=1:4, c("Control","Pseudomonas GM30", "Psudomonas GM41", "Burkholdaria BT03" ),cex.axis=2) ################################### Total biomass - ########################################################################################################################### ######### Whole BIOMASS MIXED MODEL ###### mod<-lmer(Whole.plant~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) newmod<-lmer(Whole.plant~Treatment +Initial.Root.SA+(1|EP.run), data=tritrait) # full model summary(newmod) # testing significance of adding treatment of my mod mod1<-lmer(Whole.plant~Treatment +Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(Whole.plant~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(Whole.plant~1 +Initial.Root.SA+(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(Whole.plant~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) summary(modc) summary(modb) library(lmerTest) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # par(mar=c(6,5,5,3)) boxplot(Whole.plant~Treatment, data=tritrait, xaxt="n",ylim=c(0,250),lwd=4, ylab="Total biomass (mg)", col=c("darkgreen"), cex.axis=2, cex.lab=2) axis(1, at=1:4, c("Non-inoclated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) ############### percent root model ################################################################################################################################## ######### percent ROOT MIXED MODEL ###### mod<-lmer(logperroot~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(logperroot~Treatment +tritrait$Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logperroot~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logperroot~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logperroot~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logperroot~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) # NO MODELS WERE SIGNIFICANTLY DIFFERENT HERE boxplot(percentroot~Treatment, data=tritrait, xaxt="n", ylab="percent root", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") axis(1, at=1:4, c("Control", "GM30", "GM41", "BT03" )) leveneTest(percentroot~Treatment, data=tritrait) leveneTest(logperroot~Treatment, data=tritrait) ############### percent shoot model ################################################################################################################################## tritmod<-lm(tritrait$logpershoot~tritrait$Treatment+tritrait$Initial.Root.SA) Anova(tritmod) ######### percent SHOOT MIXED MODEL ###### mod<-lmer(logpershoot~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(logpershoot~Treatment +tritrait$Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logpershoot~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logpershoot~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logpershoot~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logpershoot~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) boxplot(percentshoot~Treatment, data=tritrait, ylim=c(0,100), xaxt="n", ylab="percent shoot", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") axis(1, at=1:4, c("Control", "GM30", "GM41", "BT03" )) leveneTest(percentshoot~Treatment, data=tritrait) leveneTest(logpershoot~Treatment, data=tritrait) ############### percent leaf model ################################################################################################################################## ######### percent leaf MIXED MODEL ###### mod<-lmer(logperleaf~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(logperleaf~Treatment +tritrait$Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logperleaf~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logperleaf~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logperleaf~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logperleaf~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) boxplot(percent.leaf~Treatment, data=tritrait, ylim=c(0,100), xaxt="n", ylab="percent leaf", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") axis(1, at=1:4, c("Control", "GM30", "GM41", "BT03" )) leveneTest(percent.leaf~Treatment, data=tritrait) leveneTest(logperleaf~Treatment, data=tritrait) ############### root:leaf ratio model ################################################################################################################################## ######### root: leaf MIXED MODEL ###### mod<-lmer(logrootleaf~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(logrootleaf~Treatment +tritrait$Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logrootleaf~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logrootleaf~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logrootleaf~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logrootleaf~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) boxplot(rootleaf~Treatment, data=tritrait, xaxt="n", ylab="root to leaf ratio", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") axis(1, at=1:4, c("Control", "GM30", "GM41", "BT03" )) leveneTest(rootleaf~Treatment, data=tritrait) leveneTest(logrootleaf~Treatment, data=tritrait) ############### leaf weight ratio model ################################################################################################################################## ######### leaf.weight MIXED MODEL ###### mod<-lmer(leaf.weight~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(leaf.weight~Treatment +tritrait$Initial.Root.SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(leaf.weight~Treatment + Initial.Root.SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(leaf.weight~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(leaf.weight~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(leaf.weight~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) # leaf weight treatments significant ########### boxplot(leaf.weight~Treatment, data=tritrait, xaxt="n", ylab="leaf weight (mg)", col=c("gray100", "gray20", "gray40", "gray80"),main="P. trichocarpa") axis(1, at=1:4, c("Control", "GM30", "GM41", "BT03" )) ######### ROOT SHOOT MIXED MODEL ###### mod<-lmer(logrs~Treatment +(1|EP.run), data=tritrait) # full model newmod<-lmer(logrs~Treatment +Init_Leaf_SA + (1|EP.run), data=tritrait) # full model summary(mod) summary(newmod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logrs~Treatment + Init_Leaf_SA + (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logrs~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logrs~Init_Leaf_SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logrs~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) summary(mod1) library(lmerTest) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # boxplot(rootshoot~Treatment, data=tritrait) leveneTest(rootshoot~Treatment, data=tritrait) ######### SPAD MIXED MODEL ###### mod<-lmer(SPAD~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(SPAD~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(SPAD~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # boxplot(SPAD~Treatment, lwd=5, data=tritrait, ylim=c(15,40), ylab="Chlorophyll content (SPAD)", col=c("darkgreen"), cex.axis=2.5, cex.lab=3.5, xaxt="n") axis(1, at=1:4, c("Non-inoclated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) # ch_leaf number "CH_leafN" ######### Leaf number MIXED MODEL ###### mod<-lmer(CH_leafN~Treatment + (1|EP.run), data=tritrait) # full model summary(mod) newmod<-lmer(CH_leafN~Treatment + Init_Leaf_num+ (1|EP.run), data=tritrait) # full model summary(newmod) print(mod) boxplot(CH_leafN~Treatment, data=tritrait) # testing significance of adding treatment of my mod mod1<-lmer(CH_leafN~Treatment +Init_Leaf_num+(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(CH_leafN~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(CH_leafN~Init_Leaf_num +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(CH_leafN~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) summary(mod1) difflsmeans(newmod) # Post-hoc test to determine differences among my fixed effects # plot(difflsmeans(newmod)) # delt_rootSA "delt_RootSA" tritmod<-lm(tritrait$delt_RootSA~tritrait$Treatment) Anova(tritmod) ######### Root growth MIXED MODEL ###### mod<-lmer(delt_RootSA~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(delt_RootSA~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(delt_RootSA~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # plot(difflsmeans(mod)) boxplot(delt_RootSA~Treatment, xaxt="n", lwd=5, data=tritrait, ylim=c(0,.56), cex.lab=3.5, cex.axis=2.5, ylab="Root growth (cm^2/day)", col=c("darkgreen")) axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2.1) text(1,0.22, "B", cex=3) text(2,0.5, "A", cex=3) text(3,0.56, "B", cex=3) text(4,0.31, "B",cex=3) #delt_rootL "delt_RootL" ######### Root growth MIXED MODEL ###### mod<-lmer(delt_RootL~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) newmod<-lmer(delt_RootL~Treatment + Initial.Root.SA+(1|EP.run), data=tritrait) # full model summary(newmod) names(tritrait) # testing significance of adding treatment of my mod mod1<-lmer(delt_RootL~Treatment + Initial.Root.SA+(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(delt_RootL~Treatment +(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(delt_RootL~Initial.Root.SA +(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(delt_RootL~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # #delt_stem "logdeltstem" ######### stem growth MIXED MODEL ###### mod<-lmer(sqdeltstem~Treatment + (1|EP.run), data=tritrait) # full model summary(mod) newmod<-lmer(sqdeltstem~Treatment + Init_Leaf_SA+ (1|EP.run), data=tritrait) # full model summary(newmod) # testing significance of adding treatment of my mod mod1<-lmer(sqdeltstem~Treatment + Init_Leaf_SA+ (1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(sqdeltstem~Treatment + (1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(sqdeltstem~Init_Leaf_SA+ (1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(sqdeltstem~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # boxplot(Delt_stem~Treatment, data=tritrait, ylim=c(0,.50), ylab="stem growth (cm/day)", col=c("gray100", "gray20", "gray40", "gray80")) # SRL "logSRL" ######### SRL MIXED MODEL ###### mod<-lmer(logSRL~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logSRL~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logSRL~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # # SLA"logSLA" ######### SLA MIXED MODEL ###### mod<-lmer(logSLA~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(logSLA~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logSLA~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # plot(difflsmeans(mod)) # LEAF GROWTH - "logdeltleaf" ######### leaf growth MIXED MODEL ###### mod<-lmer(logdeltleaf~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) leveneTest(logdeltleaf~Treatment, data=tritrait) newmod<-lmer(logdeltleaf~Treatment + Init_Leaf_SA+(1|EP.run), data=tritrait) # full model summary(newmod) # testing significance of adding treatment of my mod mod1<-lmer(logdeltleaf~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(logdeltleaf~Treatment +Init_Leaf_SA+(1|EP.run), data=tritrait, REML=FALSE) modc<-lmer(logdeltleaf~Init_Leaf_SA+(1|EP.run), data=tritrait, REML=FALSE) modb<-lmer(logdeltleaf~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda, modb, modc) summary(mod1) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # boxplot(delt_LeafSA~Treatment, xaxt="n", data=tritrait, ylim=c(0,1.75), ylab="Leaf growth (cm^2/day)", col=c("darkgreen"), lwd=5, cex.axis=2.5, cex.lab=3.5) axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) text(1, .65, "B", cex=3) text(2, .82, "A", cex=3) text(3, 1.62, "A", cex=3) text(4, 1.35, "AB", cex=3) # change in average leaf size names(tritrait) boxplot(tritrait$CH_avg_leaf) normalTest(tritrait$CH_avg_leaf, method=c( "da" )) log_chg_leaf_size<-(log(tritrait$CH_avg_leaf+1)) normalTest(log_chg_leaf_size, method=c( "da" )) boxplot(log_chg_leaf_size) mod<-lmer(log_chg_leaf_size~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(log_chg_leaf_size~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(log_chg_leaf_size~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) # Post-hoc test to determine differences among my fixed effects # # final average leaf size boxplot(tritrait$Final_avg_leaf) normalTest(tritrait$Final_avg_leaf, method=c( "da" )) mod<-lmer(Final_avg_leaf~Treatment +(1|EP.run), data=tritrait) # full model summary(mod) print(mod) # testing significance of adding treatment of my mod mod1<-lmer(Final_avg_leaf~Treatment +(1|EP.run), data=tritrait, REML=FALSE) moda<-lmer(Final_avg_leaf~1 +(1|EP.run), data=tritrait, REML=FALSE) anova(mod1,moda) difflsmeans(mod) boxplot(Final_avg_leaf~Treatment, data=tritrait) ################################################################################################ ################## Figs for the trait paper ################### ################################################################################################ ################# Biomass Panel ####################### par(mar=c(5,6,1,1), srt=0) #tiff(filename="Frontiers_Henning_2016_Figure3", width=4, height=9, res=300) #par(mfrow=c(4,1)) par(mfrow=c(3,1)) par(mar=c(0,6,0.5,1)) boxplot(Whole.plant~Treatment, data=tritrait, xaxt="n",ylim=c(0,255),lwd=4, ylab="Total dry weight (mg)", col=c("gray100", "gray25","gray25", "gray25"), cex.axis=1.5, cex.lab=1.6) text(1, 177, "A",cex=3) text(2, 224, "A",cex=3) text(3, 230, "A",cex=3) text(4, 250, "A",cex=3) text(0.5, 242, "A",cex=3.5) #axis(1, at=1:4, c("Non-inoclated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) par(mar=c(0,6,0,1)) boxplot(leaf.weight~Treatment, data=tritrait, ylim=c(0,210), xaxt="n",lwd=4, ylab="Dry leaf weight (mg)", col=c("gray100", "gray25","gray25", "gray25"),cex.lab=1.5, cex.axis=1.6) #boxplot(shoot~Treatment, ylim=c(0,225),ylab="Shoot biomass (mg)", lwd=3, xaxt="n",col=c("gray100", "gray20", "gray40", "gray80"), data=tritrait, cex.lab=2, cex.axis=2) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp."),cex.axis=2) text(1, 124, "B",cex=3) text(2, 169, "A",cex=3) text(3, 155, "AB",cex=3) text(4, 182, "AB",cex=3) text(0.5, 195, "B",cex=3.5) par(mar=c(2,6,0,1)) boxplot(root.biomass~Treatment, data=tritrait,xaxt="n", ylim=c(0,90), lwd=4, ylab="Dry root weight (mg)",col=c("gray100", "gray25","gray25", "gray25"), cex.axis=1.6, cex.lab=1.5) text(1, 45, "A",cex=3) text(2, 58, "A",cex=3) text(3,67, "A",cex=3) text(4, 68, "A",cex=3) text(0.5, 82, "C",cex=3.5) axis(1, at=1:4,label=FALSE) labels<- c("Control","GM30", "GM41", "BT03" ) text(1:4,par("usr")[3] - 6, cex=1.6, labels = labels, xpd = TRUE) ################# architecture panel ######################## par(mfrow=c(4,1)) par(mar=c(0,6,0.5,1)) boxplot(CH_leafN~Treatment, xaxt="n", data=tritrait, ylab="Number of new leaves", ylim=c(-6,11), col=c("gray100", "gray25","gray25", "gray25"), lwd=4, cex.axis=1.6, cex.lab=1.5) text(1, 7.6, "B", cex=3) text(2, 8.6, "A", cex=3) text(3, 8.6, "AB", cex=3) text(4,9.6, "AB", cex=3) text(0.5, 10, "A", cex=3.5) par(mar=c(0,6,0,1)) boxplot(delt_LeafSA~Treatment, xaxt="n", data=tritrait, ylim=c(0,1.75), ylab=expression(paste("Leaf growth (cm"^2, " d"^-1,")")), col=c("gray100", "gray25","gray25", "gray25"), lwd=4, cex.axis=1.6, cex.lab=1.5) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Psudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) text(1, .74, "B", cex=3) text(2, .91, "A", cex=3) text(3, 1.69, "A", cex=3) text(4, 1.41, "AB", cex=3) text(0.5, 1.63, "B", cex=3.5) par(mar=c(0,6,0,1)) boxplot(Delt_stem~Treatment, data=tritrait, ylim=c(0,.53), xaxt="n", ylab=expression(paste("Stem growth (cm"^1, " d"^-1,")")), lwd=4, cex.axis=1.6, cex.lab=1.5,col=c("gray100", "gray25","gray25", "gray25")) text(1, 0.32, "A", cex=3) text(2, 0.50, "A", cex=3) text(3, 0.31, "A", cex=3) text(4,0.13, "A", cex=3) text(0.5, 0.45, "C", cex=3.5) par(mar=c(2,6,0,1)) boxplot(delt_RootSA~Treatment, xaxt="n", lwd=3, data=tritrait, ylim=c(0,.65), cex.lab=1.5, cex.axis=1.6, ylab=expression(paste("Root growth (cm"^2, " d"^-1,")")), col=c("gray100", "gray25","gray25", "gray25")) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Pseudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) text(1,0.27, "B", cex=3) text(2,0.55, "A", cex=3) text(3,0.605, "B", cex=3) text(4,0.365, "B",cex=3) text(0.5, 0.58, "D", cex=3.5) axis(1, at=1:4,label=FALSE) labels<- c("Control","GM30", "GM41", "BT03" ) text(1:4,par("usr")[3] - .07, cex=1.6, labels = labels, xpd = TRUE) ################################################################################################################################################## ############################################# PHYSIOLOGY DATA ########################################################################### ################################################################################################################################################## physd<-read.csv("Henning_PeerJ_Phys.csv") physd boxplot(physd$carboxylase) logcarb<-(sqrt(physd$carboxylase)) #### NORMAL### boxplot(logcarb) boxplot(physd$Asat) boxplot(physd$Amax) mod<- lm(logcarb~Treatment, data=physd) # non sig Anova(mod) physd["logcarb"]<-NA physd$logcarb<-logcarb mod<- lm(Asat~Treatment, data=physd) # non sig Anova(mod) mod<- lm(Amax~Treatment, data=physd) # no sig Anova(mod) ##################################################### ############### PHYS TRAIT PANEL ################### physd physd$Treatment<-factor(physd$Treatment, levels=c("Control", "GM30_alone", "GM41_alone", "BT03_alone")) physd par(mfrow=c(3,1)) par(mar=c(0,6,0.5,1)) boxplot(SPAD~Treatment, lwd=3, data=tritrait, ylim=c(15,40), ylab="Chlorophyll content (SPAD)", col=c("gray100", "gray25","gray25", "gray25"), cex.axis=1.5, cex.lab=1.6, xaxt="n") text(1, 32.5, "A",cex=3) text(2, 36, "A",cex=3) text(3, 39, "A",cex=3) text(4, 36.25, "A",cex=3) text(0.5, 39, "A",cex=3.5) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Pseudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) par(mar=c(0,6,0,1)) boxplot(physd$logcarb~physd$Treatment, xaxt="n", lwd=4,ylim=c(0.05,0.32), cex.axis=1.5, cex.lab=1.6, col=c("gray100", "gray25","gray25", "gray25"), ylab= expression(paste("ɸCO"[2]," (Quantum yield of photo.)"))) text(4, 0.215, "A",cex=3) text(1, 0.211, "A",cex=3) text(2, 0.247, "A",cex=3) text(3, 0.288, "A",cex=3) text(0.5, 0.303, "B",cex=3.5) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Pseudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) par(mar=c(2,6,0,1)) boxplot(physd$Amax~physd$Treatment, xaxt="n", lwd=4, ylim=c(3,17), col=c("gray100", "gray25","gray25", "gray25"), cex.axis=1.5, cex.lab=1.6, ylab=expression(paste("A"[max],"(µmol CO"[2]," m"^-2," s"^-1,")") )) #axis(1, at=1:4, c("Non-inoculated Control","Pseudomonas sp.1", "Pseudomonas sp.2", "Burkholdaria sp." ),cex.axis=2) ylab=expression(paste("Leaf growth (cm"^2, " d"^-1,")")), c text(4, 12.2, "A",cex=3) text(1, 13.7, "A",cex=3) text(2, 14.8, "A",cex=3) text(3, 15.3, "A",cex=3) text(0.5, 15.8, "C",cex=3.5) axis(1, at=1:4,label=FALSE) labels<- c("Control","GM30", "GM41", "BT03" ) text(1:4,par("usr")[3] - 1, cex=1.8, labels = labels, xpd = TRUE) ################################################################################################ ################################################################################################ # Microbial Abundance Data ################################################################################################ ################################################################################################ cfu<-read.csv("Henning_PeerJ_microbial_abundance_data.csv") cfu names(cfu) boxplot(CFU.g.tissue~Treatments*tissue, data=cfu) logCFU<-log(cfu$CFU.g.tissue+1) boxplot(logCFU~tissue*Treatments, data=cfu) cfu["logCFU"]<-NA cfu$logCFU<-logCFU cfu$Treatments<-factor(cfu$Treatments, levels=c("Control", "GM30", "GM41", "Bt03")) boxplot(CFU.g.tissue~tissue*Treatments, data=cfu, las=2) mod<-lm(CFU.g.tissue~Treatments*dilution, data=cfu) Anova(mod) summary(mod) first<-cfu[which(cfu$dilution=="1"),] second<-cfu[which(cfu$dilution=="2"),] third<-cfu[which(cfu$dilution=="3"),] mod<-lm(CFU.g.tissue~Treatments*tissue, data=first) Anova(mod) summary(mod) mod2<-lm(CFU.g.tissue~Treatments*tissue, data=second) Anova(mod2) summary(mod2) mod3<-lm(CFU.g.tissue~Treatments*tissue, data=third) Anova(mod3) summary(mod3)