library(MuMIn) library(ggplot2) temperature<-read.csv((file.choose()), header=T) #tempterature_permafrost_paper.csv ###data preparation #calculate mean for plots meantemp12<-aggregate(temp12 ~ plot, data=temperature, mean) meanal<-aggregate(active.layer ~ plot, data=temperature, mean) meanfire<-aggregate(fire ~ plot, data=temperature, mean) meanyear<-aggregate(year ~ plot, data=temperature, mean) temperature$pairs<-as.numeric(gsub("P", "", temperature$pairs)) meanpairs<-aggregate(pairs ~ plot, data=temperature, mean) meantemp<-as.data.frame(cbind(temp12=meantemp12$temp12, meanal,fire=meanfire$fire, year=meanyear$year, pairs=meanpairs$pairs )) #calculate rates fire<-subset(meantemp, fire==1) nofire<-subset(meantemp, fire==0) ratetemp12<-(fire$temp12-nofire$temp12) rateal<-(fire$active.layer-nofire$active.layer) ratedata<-fire ratedata$ratetemp12<-ratetemp12 ratedata$rateal<-rateal ratedata$year_sf<-ratedata$year[ratedata$year == 1] <- 12 ratedata$year_sf<-ratedata$year[ratedata$year == 0.5] <- 28 ratedata$year_sf<-ratedata$year[ratedata$year == 0] <- 44 temp12<-as.data.frame(cbind(type=rep("temp12", 30), rate=as.character(as.numeric(ratetemp12)), year=as.character(ratedata$year), plot=as.character(ratedata$plot))) tempal<-as.data.frame(cbind(type=rep("al", 30), rate=as.character(as.numeric(rateal)), year=as.character(ratedata$year), plot=as.character(ratedata$plot))) ###models ###temperature ###################################################################### fyear<-as.factor(temp12$year) ratetemp12<-as.numeric(as.character(temp12$rate)) mod1<-lm(ratetemp12~fyear, temp12) summary(mod1) r.squaredGLMM(mod1) ###check model assumptions par(mfrow=c(2,2)) plot(mod1) ###simulation mod1<-lm(ratetemp12~fyear, temp12) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@coef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(fyear=factor(c("12","28","44"), levels=levels(allrates$fyear))) newdat Xmat<-model.matrix(~fyear, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@coef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% coef(mod1) newdat #results bsim@coef sum(bsim@coef[,1]>0)/nsim sum(bsim@coef[,2]+bsim@coef[,1]>0)/nsim sum(bsim@coef[,3]+bsim@coef[,1]>0)/nsim restemp12<-round(apply(bsim@coef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) #estimates sum((bsim@coef[,1])>0)/nsim sum((bsim@coef[,2])>0)/nsim sum((bsim@coef[,3])>0)/nsim round(apply(bsim@coef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) st_out<-as.data.frame(round(apply(bsim@coef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) st_out ###plot fyear<-as.factor(temp12$year) ratetemp12<-as.numeric(as.character(temp12$rate)) legend_title<-"years after fire" plot2<-ggplot(data = temp12, aes(x =fyear, y =ratetemp12))+ labs(x="")+ labs(y="Soil temperature (°C)")+ scale_y_continuous(limits = c(-5, 12))+ geom_hline(yintercept=0, col="grey40")+ 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(fyear)), position=position_dodge(width=0.5), alpha=0.5) + scale_x_discrete(labels=c("12" = "12", "28" = "28","44" = ">44"))+ geom_point(data=newdat, aes(x = fyear, y = fit) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = fyear, y = fit, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ # theme(legend.position="n")+ theme(legend.position=c(0.5, 0.1),legend.background = element_rect(fill=alpha("white",1)),legend.direction = "horizontal",legend.title=element_text(size=7), legend.text=element_text(size=7))+ #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() ###permafrost thaw depth################################################################ fyear<-as.factor(tempal$year) raten<-as.numeric(as.character(tempal$rate)) raten mod1<-lm(raten~fyear, tempal) summary(mod1) mod1 print(mod1, correlation=TRUE) summary(mod1)$sigma round(coef(mod1),3) ###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="normtempalQQ")) qqline(resid(mod1)) scatter.smooth(fitted(mod1),sqrt(abs(resid(mod1)))) ###simulation mod1<-lm(raten~fyear, tempal) summary(mod1) nsim<-2000 bsim<-sim(mod1, n.sim=nsim) str(bsim) round(apply(bsim@coef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) newdat<-expand.grid(fyear=factor(c("12","28","44"), levels=levels(allrates$fyear))) newdat Xmat<-model.matrix(~fyear, data=newdat) fitmat<-matrix(ncol=nsim, nrow=nrow(newdat)) for(i in 1:nsim) fitmat[,i]<-Xmat %*% bsim@coef[i,] newdat$lower<-apply(fitmat, 1, quantile, prob=0.025) newdat$upper<-apply(fitmat, 1, quantile, prob=0.975) newdat$fit <- Xmat %*% coef(mod1) newdat ###plot fyear<-as.factor(tempal$year) raten<-as.numeric(as.character(tempal$rate)) legend_title<-"years after fire" plot3<-ggplot(data = tempal, aes(x =fyear, y =raten))+ labs(x="")+ labs(y="Permafrost thaw depth (cm)")+ scale_y_continuous(limits = c(-60, 58))+ geom_hline(yintercept=0, col="grey40")+ scale_color_manual( legend_title,values=c("red4", "lightpink2", "azure4"), labels=c("12", "28", ">44"))+ #geom_point(aes(color=factor(ftype), tempalpha=factor(fyear)), position=position_dodge(width=0.5)) + geom_point(aes(color=factor(fyear)), position=position_dodge(width=0.5), alpha=0.5) + scale_x_discrete(labels=c("12" = "12", "28" = "28","44" = ">44"))+ geom_point(data=newdat, aes(x = fyear, y = fit) , position=position_dodge(width=0.5), col="black", shape=16, size=1.5)+ geom_linerange(data=newdat, aes(x = fyear, y = fit, ymin=lower, ymax=upper), position=position_dodge(width=0.5), col="black")+ theme_bw()+ 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() #results bsim@coef sum(bsim@coef[,1]>0)/nsim sum(bsim@coef[,2]+bsim@coef[,1]>0)/nsim sum(bsim@coef[,3]+bsim@coef[,1]>0)/nsim resal<-round(apply(bsim@coef, 2, quantile, prob=c(0.025, 0.5, 0.975)),3) resal #estimates sum((bsim@coef[,1])>0)/nsim sum((bsim@coef[,2])>0)/nsim sum((bsim@coef[,3])>0)/nsim round(apply(bsim@coef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3) ptd_out<-as.data.frame(round(apply(bsim@coef, 2 ,quantile, prob=c(0.025, 0.5, 0.975)),3)) ptd_out ######################################################figure1############ plot_grid(plot2, plot3 , ncol = 2, rel_widths = 1, rel_heights = 1, labels = c('A', 'B'))+ theme(plot.margin = unit(c(0,0.3,0.3, 0.3), "cm")) fig1<-plot_grid(plot2, plot3 , ncol = 2, rel_widths = 1, rel_heights = 1, labels = c('A', 'B'))+ cairo_ps("Figure 1.eps", height=3, width=6, fallback_resolution = 800) dev.off() ##################extract estimates soil########################### estimates_soil<-rbind(st_out, ptd_out) write.table(estimates_soil, "estimates_soil.csv", sep=",")