traits<-read.csv((file.choose()), header=T) #choose plant_functional_traits.csv head(traits) ###load packages library(ggplot2) library(cowplot) library(lme4) library(arm) library(blmeco) library(MuMIn) ### translation paper -> script #Mean cover -> diameter #Height -> height #Bladelength -> leaflength #Leaf thickness -> leafthick #SLA -> SLA #LDMC -> LDMC ###subsets for Betula nana and Vaccinium uliginosum betnan<-traits[traits$species=="Betula nana",] vaculi<-traits[traits$species=="Vaccinium uliginosum",] alltraits<-as.data.frame(rbind(betnan, vaculi)) head(alltraits) ###subsets for burnt/control fire<-subset(alltraits, fire=="fire") nofire<-subset(alltraits, fire=="nofire") ###NA's? subset(fire, is.na(height)) subset(fire, is.na(diametermean)) subset(fire, is.na(leaflength)) subset(fire, is.na(leafthick)) subset(fire, is.na(LDMC)) subset(fire, is.na(SLA)) ###data preparation for height, leaf length and leaf thickness model rateheight<-(fire$height-nofire$height) rateleaflength<-(fire$leaflength-nofire$leaflength) rateleafthick<-(fire$leafthick-nofire$leafthick) ratedat<-fire ratedat$rateheight<-rateheight ratedat$rateleaflength<-rateleaflength ratedat$rateleafthick<-rateleafthick height<-as.data.frame(cbind(trait=rep("height", 60), rate=as.character(as.numeric(rateheight)), year=as.character(ratedat$year), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) leaflength<-as.data.frame(cbind(trait=rep("leaflength", 60), rate=as.character(as.numeric(rateleaflength)), year=as.character(ratedat$year), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) leafthick<-as.data.frame(cbind(trait=rep("leafthick", 60), rate=as.character(as.numeric(rateleafthick)), year=as.character(ratedat$year), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) ###data preparation for diamater model traits$diametermean #delete NA's traitsdiameter<-traits[-c(145, 146), ] traitsdiameter betnand<-traitsdiameter[traitsdiameter$species=="Betula nana",] vaculid<-traitsdiameter[traitsdiameter$species=="Vaccinium uliginosum",] alltraitsd<-as.data.frame(rbind(betnand, vaculid)) fired<-subset(alltraitsd, fire=="fire") nofired<-subset(alltraitsd, fire=="nofire") ratedatd<-fired ratediameter<-(fired$diametermean-nofired$diametermean) ratedatd$ratediameter<-ratediameter length(ratedatd$ratediameter) diameter<-as.data.frame(cbind(trait=rep("diameter", 59), rate=as.character(as.numeric(ratediameter)), year=as.character(ratedatd$year), plot=as.character(ratedatd$plot), species=as.character(ratedatd$species))) ###data preparation for SLA LDMC model ###delete rows, where NA (in SLA & LDMC) and rows of related pair traits1<-traits[-c(7,8,15,16,19,20,31,32,151,152), ] ###subsets betnan<-traits1[traits1$species=="Betula nana",] dim(betnan) vaculi<-traits1[traits1$species=="Vaccinium uliginosum",] dim(vaculi) alltraits1<-as.data.frame(rbind(betnan, vaculi)) head(alltraits1) dim(alltraits1) #claculate rates fire<-subset(alltraits1, fire=="fire") dim(fire) nofire<-subset(alltraits1, fire=="nofire") dim(nofire) #check subset(fire, is.na(LDMC)) subset(fire, is.na(SLA)) rateSLA<-(fire$SLA-nofire$SLA) rateLDMC<-(fire$LDMC-nofire$LDMC) ratedat1<-fire ratedat1$rateSLA<-rateSLA ratedat1$rateLDMC<-rateLDMC SLA<-as.data.frame(cbind(trait=rep("SLA", 55), rate=as.character(as.numeric(rateSLA)), year=as.character(ratedat1$year), plot=as.character(ratedat1$plot), species=as.character(ratedat1$species))) LDMC<-as.data.frame(cbind(trait=rep("LDMC", 55), rate=as.character(as.numeric(rateLDMC)), year=as.character(ratedat1$year), plot=as.character(ratedat1$plot), species=as.character(ratedat1$species))) ###models/results/plots############################################################### #####diameter#################################################################### diameter$rate<-as.numeric(as.character(diameter$rate)) diameter$species<-as.factor(diameter$species) diameter$year<-as.factor(diameter$year) mod1<-lmer(rate~species+year+species:year+(1|plot), diameter, REML=T) summary(mod1) r.squaredGLMM(mod1) mod1 print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###check model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), diameter, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(diameter$species)), year=factor(c("12","28","44"), levels=levels(diameter$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_diameter.csv") ###results head(bsim@fixef) sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim diameter_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) diameter_out ###plot diameter$rate<-as.numeric(diameter$rate) legend_title<-"years after fire" diameter$species <- factor(diameter$species, levels=c( "Vaccinium uliginosum", "Betula nana")) plot1<-ggplot(data = diameter, aes(x =species, y =rate))+ labs(x="")+ labs(y="Mean cover (cm)")+ scale_y_continuous(limits = c(-50, 160))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(axis.text.x = element_blank())+ theme(legend.position="n")+ #theme(legend.position=c(0.8, 0.2))+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot1 plot1 dev.off() #####height#################################################################### height$rate<-as.numeric(as.character(height$rate)) height$rate height$species<-as.factor(height$species) height$year<-as.factor(height$year) mod1<-lmer(rate~species+year+species:year+(1|plot), height, REML=T) summary(mod1) r.squaredGLMM(mod1) mod1 print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###check model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), height, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(height$species)), year=factor(c("12","28","44"), levels=levels(height$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_height.csv") ###results bsim@fixef sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) height_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) height_out ###plot height$rate<-as.numeric(height$rate) legend_title<-"years after fire" height$species <- factor(height$species, levels=c("Vaccinium uliginosum", "Betula nana")) plot2<-ggplot(data = height, aes(x =species, y =rate))+ labs(x="")+ labs(y="Height (cm)")+ scale_y_continuous(limits = c(-55, 75))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(axis.text.x = element_blank())+ theme(legend.position="n")+ #theme(legend.position=c(0.8, 0.2))+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot2 dev.off() #####leaflength#################################################################### leaflength$rate<-as.numeric(as.character(leaflength$rate)) leaflength$rate leaflength$species<-as.factor(leaflength$species) leaflength$year<-as.factor(leaflength$year) mod1<-lmer(rate~species+year+species:year+(1|plot), leaflength, REML=T) summary(mod1) r.squaredGLMM(mod1) mod1 print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), leaflength, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Vaccinium uliginosum","Betula nana"), levels=levels(leaflength$species)), year=factor(c("12","28","44"), levels=levels(leaflength$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_leaflength.csv") bsim@fixef sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) leaflength_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) leaflength_out leaflength$rate<-as.numeric(leaflength$rate) legend_title<-"years after fire" leaflength$species <- factor(leaflength$species, levels=c("Vaccinium uliginosum", "Betula nana")) plot3<-ggplot(data = leaflength, aes(x =species, y =rate))+ labs(x="")+ labs(y="Blade length (cm)")+ scale_y_continuous(limits = c(-1, 1.25))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(axis.text.x = element_blank())+ theme(legend.position="n")+ #theme(legend.position=c(0.8, 0.2))+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot3 dev.off() #####leafthick#################################################################### leafthick$rate<-as.numeric(as.character(leafthick$rate)) leafthick$rate leafthick$species<-as.factor(leafthick$species) leafthick$year<-as.factor(leafthick$year) mod1<-lmer(rate~species+year+species:year+(1|plot), leafthick, REML=T) r.squaredGLMM(mod1) summary(mod1) mod1 print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###check model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), leafthick, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(leafthick$species)), year=factor(c("12","28","44"), levels=levels(leafthick$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_leafthick.csv") ###results bsim@fixef sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) leafthick_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) leafthick_out leafthick$rate<-as.numeric(leafthick$rate) legend_title<-"years after fire" ###plot leafthick$species <- factor(leafthick$species, levels=c("Vaccinium uliginosum", "Betula nana")) plot4<-ggplot(data = leafthick, aes(x =species, y =rate))+ labs(x="")+ labs(y="Leaf thickness (mm)")+ scale_y_continuous(limits = c(-0.11, 0.15))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(axis.text.x = element_blank())+ #theme(legend.position=c(0.8, 0.2))+ theme(legend.position="n")+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot4 dev.off() #####SLA#################################################################### SLA$rate<-as.numeric(as.character(SLA$rate)) SLA$rate SLA$species<-as.factor(SLA$species) SLA$year<-as.factor(SLA$year) mod1<-lmer(rate~species+year+species:year+(1|plot), SLA, REML=T) summary(mod1) mod1<-lm(rate~species+year+species:year, SLA) print(mod1, correlation=TRUE) r.squaredGLMM(mod1) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), SLA, REML=T) summary(mod1) r.squaredGLMM(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Betula nana","Vaccinium uliginosum"), levels=levels(SLA$species)), year=factor(c("12","28","44"), levels=levels(SLA$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_SLA.csv") ###results bsim@fixef sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) SLA_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) SLA_out ###plot SLA$rate<-as.numeric(SLA$rate) legend_title<-"years after fire" SLA$species <- factor(SLA$species, levels=c("Vaccinium uliginosum", "Betula nana")) plot5<-ggplot(data = SLA, aes(x =species, y =rate))+ labs(x="")+ labs(y = expression(paste("SLA (cm"^{2}, "/g)")))+ scale_y_continuous(limits = c(-100, 100))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(legend.position="n")+ theme(axis.text.x = element_text(face = "italic"))+ #theme(legend.position=c(0.8, 0.2))+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot5 dev.off() #####LDMC#################################################################### LDMC$rate<-as.numeric(as.character(LDMC$rate)) LDMC$rate LDMC$species<-as.factor(LDMC$species) LDMC$year<-as.factor(LDMC$year) mod1<-lmer(rate~species+year+species:year+(1|plot), LDMC, REML=T) summary(mod1) mod1 print(mod1, correlation=TRUE) r.squaredGLMM(mod1) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rate~species+year+species:year+(1|plot), LDMC, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) bsim round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(species=factor(c("Betula nana","Vaccinium uliginosum"), levels=levels(LDMC$species)), year=factor(c("12","28","44"), levels=levels(LDMC$year))) newdat Xmat<-model.matrix(~species+year+species:year, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) newdat write.table(newdat,"newdat_LDMC.csv") ###results bsim@fixef sum((bsim@fixef[,1]+bsim@fixef[,2]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,3]+bsim@fixef[,5]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,2]+bsim@fixef[,4]+bsim@fixef[,6]>0)/nsim) sum((bsim@fixef[,1]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,3]>0)/nsim) sum((bsim@fixef[,1]+bsim@fixef[,4]>0)/nsim) #estimates sum((bsim@fixef[,1])>0)/nsim sum((bsim@fixef[,2])>0)/nsim sum((bsim@fixef[,3])>0)/nsim sum((bsim@fixef[,4])>0)/nsim sum((bsim@fixef[,5])>0)/nsim sum((bsim@fixef[,6])>0)/nsim round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) LDMC_out<-as.data.frame(round(apply(bsim@fixef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) LDMC_out ###plot LDMC$rate<-as.numeric(LDMC$rate) legend_title<-"years after fire:" LDMC$species <- factor(LDMC$species, levels=c("Vaccinium uliginosum", "Betula nana")) plot6<-ggplot(data = LDMC, aes(x =species, y =rate))+ labs(x="")+ labs(y = expression(paste("LDMC (mg/g)")))+ scale_y_continuous(limits = c(-500, 500))+ geom_hline(yintercept=0, col="grey40")+ #scale_color_manual( legend_title,values=c("grey40", "yellow4", "yellowgreen", "green4", "darkgreen", "turquoise4"))+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), alpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(year)), position=position_dodge(width=0.5), alpha=0.5) + geom_point(data=newdat, aes(x = species, y = fit, group=year) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = species, y = fit, group=year, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ theme(axis.text.x = element_text(face = "italic"))+ theme(legend.position=c(0.5, 0.07),legend.background = element_rect(fill=alpha(0.1)),legend.direction = "horizontal",legend.title=element_text(size=7), legend.text=element_text(size=7))+ #theme(legend.position="n")+ theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ theme(plot.margin = unit(c(0.1,0.1,0.1, 0.1), "cm")) plot6 dev.off() #############################fig paper############################################################ ###fig.3 plot_grid(plot1, plot2, plot3 ,plot4, plot5, plot6, ncol = 2,align="v", rel_widths = 1, rel_heights = c(1,1,1.08), labels = c('A', 'B', 'C', 'D', 'E', 'F'))+ theme(plot.margin = unit(c(0,0.3,0.3, 0.3), "cm")) fig3<-plot_grid(plot1, plot2, plot3 ,plot4, plot5, plot6, ncol = 2,align="v", rel_widths = 1, rel_heights = c(1,1,1.08), labels = c('A', 'B', 'C', 'D', 'E', 'F'))+ cairo_ps("Figure 3.eps", height=8, width=6, fallback_resolution = 800) ############################table estimates###################################################### estimates_traits<-rbind(diameter_out, height_out, leaflength_out, leafthick_out, SLA_out, LDMC_out) write.table(estimates_traits, "estimates_traits.csv", sep=",") ################################################################################ ##############################numeric analysis################################# ratedat$yearN<-as.numeric(as.character(ratedat$year)) ratedatd$yearN<-as.numeric(as.character(ratedatd$year)) ratedat1$yearN<-as.numeric(as.character(ratedat1$year)) ###data preparation for height, leaf length and leaf thickness model height<-as.data.frame(cbind(trait=rep("height", 60), rate=as.numeric(rateheight), yearN=as.numeric(ratedat$yearN), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) height$yearN<-as.numeric(as.character(height$year)) height$rate<-as.numeric(as.character(height$rate)) height$species<-as.factor(height$species) height$species<-factor(height$species, levels=c("Vaccinium uliginosum", "Betula nana")) leaflength<-as.data.frame(cbind(trait=rep("leaflength", 60), rate=as.numeric(rateleaflength), yearN=as.numeric(ratedat$yearN), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) leaflength$yearN<-as.numeric(as.character(leaflength$year)) leaflength$rate<-as.numeric(as.character(leaflength$rate)) leaflength$species<-as.factor(leaflength$species) leaflength$species<-factor(leaflength$species, levels=c("Vaccinium uliginosum", "Betula nana")) leafthick<-as.data.frame(cbind(trait=rep("leafthick", 60), rate=as.numeric(rateleafthick), yearN=as.numeric(ratedat$yearN), plot=as.character(ratedat$plot), species=as.character(ratedat$species))) leafthick$yearN<-as.numeric(as.character(leafthick$year)) leafthick$rate<-as.numeric(as.character(leafthick$rate)) leafthick$species<-as.factor(leafthick$species) leafthick$species<-factor(leafthick$species, levels=c("Vaccinium uliginosum", "Betula nana")) ###data preparation for diamater model diameter<-as.data.frame(cbind(trait=rep("diameter", 59), rate=as.numeric(ratediameter), yearN=as.numeric(as.character(ratedatd$yearN)), plot=as.character(ratedatd$plot), species=as.character(ratedatd$species))) diameter$yearN<-as.numeric(as.character(diameter$year)) diameter$rate<-as.numeric(as.character(diameter$rate)) diameter$species<-as.factor(diameter$species) diameter$species<-factor(diameter$species, levels=c("Vaccinium uliginosum", "Betula nana")) ###data preparation for SLA LDMC model SLA<-as.data.frame(cbind(trait=rep("SLA", 55), rate=as.numeric(rateSLA), yearN=as.numeric(ratedat1$yearN), plot=as.character(ratedat1$plot), species=as.character(ratedat1$species))) SLA$yearN<-as.numeric(as.character(SLA$year)) SLA$rate<-as.numeric(as.character(SLA$rate)) SLA$species<-as.factor(SLA$species) SLA$species<-factor(SLA$species, levels=c("Vaccinium uliginosum", "Betula nana")) LDMC<-as.data.frame(cbind(trait=rep("LDMC", 55), rate=as.numeric(rateLDMC), yearN=as.numeric(ratedat1$yearN), plot=as.character(ratedat1$plot), species=as.character(ratedat1$species))) LDMC$yearN<-as.numeric(as.character(LDMC$year)) LDMC$rate<-as.numeric(as.character(LDMC$rate)) LDMC$species<-as.factor(LDMC$species) LDMC$species<-factor(LDMC$species, levels=c("Vaccinium uliginosum", "Betula nana")) ###activate for 2nd scenario #height$yearN[height$yearN==44]<-100 #leaflength$yearN[leaflength$yearN==44]<-100 #leafthick$yearN[leafthick$yearN==44]<-100 #diameter$yearN[diameter$yearN==44]<-100 #SLA$yearN$yearN[SLA$yearN$yearN==44]<-100 #LDMC$yearN[LDMC$yearN==44]<-100 #########################diameter############################################ ###model mod1<-lmer(ratediameter~yearN+species+yearN:species+(1|plot), diameter, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(ratediameter~yearN+species+yearN:species+(1|plot), diameter, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) diameter$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(diameter$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_diameter<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim) #########################height############################################ ###model mod1<-lmer(rateheight~yearN+species+yearN:species+(1|plot), height, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rateheight~yearN+species+yearN:species+(1|plot), height, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) height$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(height$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_height<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim) #########################leaflength############################################ ###model mod1<-lmer(rateleaflength~yearN+species+yearN:species+(1|plot), leaflength, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rateleaflength~yearN+species+yearN:species+(1|plot), leaflength, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) leaflength$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(leaflength$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_leaflength<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim) #########################leafthick############################################ ###model mod1<-lmer(rateleafthick~yearN+species+yearN:species+(1|plot), leafthick, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rateleafthick~yearN+species+yearN:species+(1|plot), leafthick, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) leafthick$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(leafthick$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_leafthick<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim) #########################SLA############################################ ###model mod1<-lmer(rateSLA~yearN+species+yearN:species+(1|plot), SLA, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rateSLA~yearN+species+yearN:species+(1|plot), SLA, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) SLA$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(SLA$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_SLA<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim) #########################LDMC############################################ ###model mod1<-lmer(rateLDMC~yearN+species+yearN:species+(1|plot), LDMC, REML=T) summary(mod1) r.squaredGLMM(mod1) summary(mod1) print(mod1, correlation=TRUE) summary(mod1)$sigma round(fixef(mod1),3) ranef(mod1) ###model assumptions par(mfrow=c(2,2)) scatter.smooth(fitted(mod1), resid(mod1));abline(h=0, lty=2) title("Tuckey-Anscombe Plot") qqnorm(resid(mod1, main="normalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) qqnorm(ranef(mod1)$plot[,1], main="normap QQ random effects") qqline(ranef(mod1)$plot[,1]) ###simulation mod1<-lmer(rateLDMC~yearN+species+yearN:species+(1|plot), LDMC, REML=T) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@fixef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) LDMC$species newdat<-expand.grid(species=factor(c("Vaccinium uliginosum", "Betula nana"), levels=levels(LDMC$species)), yearN=seq(12,44,length=100)) newdat Xmat<-model.matrix(~yearN+species+yearN:species, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@fixef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% fixef(mod1) write.table(newdat,"newdat.csv") newdat_LDMC<-newdat ###results fixef(mod1) summary(mod1) bsim@fixef[,2] ###slope Vaccinium sum((bsim@fixef[,2]>0)/nsim) ###slope Betula sum((bsim@fixef[,2]+bsim@fixef[,4]>0)/nsim)