require(ggplot2);require(dismo);require(gbm);require(effects);require(data.table);require(plotmo);require(plyr) library(pROC); require(dplyr);require(zoo);require(strucchange);require(car);require(cowplot);require(vegan) ## datafile with annual recruitment rate at each site, annual climate in years following fire, fire severity, distance to seed source etc. dat<-read.csv("Davis_et_al_recruitment_plot_data.csv") ## time series of climate data at all sites from 1981-2015 dat.c<-read.csv("Davis_et_al_climate_timeseries.csv") # create column for presence or absence of any recuitment in each year dat$PA<-ifelse(dat$density>0,1,0) ########## calculate z-scores for deficit and vpd zscore<-function(x){(x-mean(x,na.rm=T))/sd(x,na.rm=T)} # vpd and def z-score in climate dataset dat.c<- data.frame(dat.c %>% group_by(site) %>% mutate_at(c("vpd_mean678","def"), zscore)) # add climate data to tree dataset dat<-merge(dat, dat.c, by=c("site","year")) head(dat) ## subset data so that just have trees that established following fire and before 2016 (we didn't have 2016 climate data) # note that CO sites climate and tree data only go up to 2011 when sites were sampled dat<-dat[dat$year<2016,] dat<-dat[dat$time_since_fire>0,] ##### PIPO models ### datpi<-dat[dat$species=="PIPO",] # select just sites with PIPO datpi$site<-factor(datpi$site) # this removes PSME only sites from the levels in the dataframe datpi$fire<-factor(datpi$fire) # this removes PSME only fires from the levels in the dataframe datpi$cv_group<-factor(datpi$cv_group) # this removes PSME only cv.groups from the levels in the dataframe datpi<-transform(datpi, cv.num=match(cv_group, unique(cv_group))) #give each fire/cv_group a number to use in gbm.step to assign folds ## calculate and add annual recruitment rate threshold for each region qc<-function(x){quantile(x,c(0.25,0.5))} tapply(datpi[datpi$density>0,]$density,datpi[datpi$density>0,]$region,qc) datpi$q25<-ifelse(datpi$region=="CO",0.0016,ifelse(datpi$region=="NR",0.0042,ifelse(datpi$region=="CA",0.0017,0.0013))) # create binomial response based on regional annual recruitment rate thresholds datpi$PA25<-ifelse(datpi$density>=datpi$q25,1,0) mean(datpi$PA25) tapply(datpi$PA25,list(datpi$PA25,datpi$region),length) names(datpi) #vars<-c(7,19,28,32,33,34,35,38)# all variables vars<-c(7,19,28,33,35) # simplified final model ## binomial gbm model m.pi<-gbm.step(data=datpi, gbm.x = vars, gbm.y = 41, family = "bernoulli", tree.complexity = 3, learning.rate = 0.002, bag.fraction = 0.75, fold.vector=datpi[,39],n.folds=length(unique(datpi[,39]))) # note that it is important to include the fire as the folds for cv, without it the model is completely overfit and the # cv AUC is artificially high. ## look at partial dependency plots and relative importance gbm.plot(m.pi, write.title = FALSE) ## model AUC datpi$pr<-predict.gbm(m.pi, datpi, n.trees=m.pi$gbm.call$best.trees, type="response") roccurve <- roc(datpi$PA25 ~ datpi$pr) plot(roccurve) auc(roccurve) ############################################################## ##### Validate model by holding out each site and calculating accuracy & AUC library(hmeasure) # H measure is a coherent AUC alternative # create a function to hold out each site, fit the model, predict to the holdpout site, and caluclate accuracy and AUC holdout1.site<-function(sites,dat,metrics){ for (i in sites){ test.sites<- i metrics[i,5]<-test.sites dat.test<-dat[dat$site==test.sites,] dat.train<-dat[dat$site!=test.sites,] metrics[i,6]<-paste(dat.test[1,"region"]) gbm1 <- gbm(PA25~time_since_fire+dist_PIPO_seed_source+dnbr+vwc_mean_driestM+vpd_mean678, data=dat.train, distribution="bernoulli", n.trees=trees, shrinkage=lr,interaction.depth=id, bag.fraction = bf) actual <- dat.test$PA25 preds.g <- predict.gbm(gbm1, dat.test, n.trees=trees, type="response") if(sum(dat.test$PA25)>0 & sum(dat.test$PA25)/length(dat.test$PA25)!=1){ roccurve <- roc(dat.test$PA25 ~ preds.g) thresh=coords(roccurve, "b", ret="t", best.method="youden") thresh=thresh[1] metrics[i,4]<-thresh results <- HMeasure(actual,preds.g) metrics[i,2]<-results$metrics['H'] metrics[i,3]<-auc(roccurve) } else { metrics[i,2:3]<-NA thresh=0.25 metrics[i,4]<-thresh } pr <- as.numeric(predict.gbm(gbm1, dat.test, n.trees=trees, type="response")>thresh) #actual;pr;preds.g;dat.test$proportion cm = as.matrix(table(Actual = actual, Predicted = pr)) n = sum(cm) # number of instances diag = diag(cm) # number of correctly classified instances per class if (sum(actual)==0 & sum(pr)==length(pr) | sum(actual)==length(actual) & sum(pr)==0){ accuracy=0 k=NA } else{ accuracy = sum(diag) / n dat.test$pr<-pr } metrics[i,1]<-accuracy } assign("metrics.sum", metrics, envir = .GlobalEnv) } # set settings from above gbm.step model bf=0.75 lr=0.002 trees=1550 id=3 # get names of PIPO sites sites<-levels(datpi$site) # make data frame for output metrics metrics<-setNames(data.frame(matrix(ncol =6, nrow = 0)), c("acc","H","AUC","thresh","site","region")) # run function holdout1.site(sites,datpi,metrics) metrics.sum colMeans(metrics.sum[c(1:4)],na.rm=T) # get means by region and test for differences # accuracy metrics.sum$region<-as.factor(metrics.sum$region) tapply(metrics.sum$acc,metrics.sum$region,mean) an<-aov(metrics.sum$acc~metrics.sum$region) summary(an) TukeyHSD(an) #AUC tapply(metrics.sum$AUC,metrics.sum$region,mean,na.rm=T) an<-aov(metrics.sum$AUC~metrics.sum$region) summary(an) TukeyHSD(an) # H measure tapply(metrics.sum$H,metrics.sum$region,mean,na.rm=T) an<-aov(metrics.sum$H~metrics.sum$region) summary(an) TukeyHSD(an) ############################################## ### predict recruitment probability for PIPO with time since fire = 1 ## add region to dat.c from dat.pi dataframe, then get rid of sites without region # which is a round about way of selecting only sites used for PIPO model in the dat.c dataframe to then used for hindcasting dat.c.pi<-dat.c nm <- c("region") dat.c.pi[nm] <- lapply(nm, function(x) datpi[[x]][match(dat.c.pi$site, datpi$site)]) dat.c.pi<-dat.c.pi[is.na(dat.c.pi$region)==F,] #take out sites not used in PIPO analysis dat.c.pi$site<-factor(dat.c.pi$site) # make regional means for dat.c.pi for plotting purposes dat.c.pi.reg<- dat.c.pi %>% group_by(region,year) %>% summarize_at(c("tskin_max","vpd_mean678","vwc_mean345","vwc_mean_driestM"), mean) dat.c.pi.reg<-data.frame(dat.c.pi.reg) ### function to predict probability of pipo regen. ("pr") pred.brt1<-function(sites,dat,nd.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PIPO_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean_driestM<-rep(dat.test$vwc_mean_driestM) nd$vpd_mean678<-rep(dat.test$vpd_mean678) # current vpd nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.pi, nd, n.trees=m.pi$gbm.call$best.trees, type="response") nd$Site<-i nd.sum<-rbind(nd.sum,nd) } assign("nd.fin.cur", nd.sum, envir = .GlobalEnv) } nd.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_PIPO_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.pi$site) pred.brt1(sites,dat.c.pi,nd.sum) head(nd.fin.cur) ### summarize by region ### nd.reg.cur<- nd.fin.cur %>% group_by(region,year) %>% summarize(pr=mean(pr)) nd.reg.cur<-data.frame(nd.reg.cur) head(nd.reg.cur) # set colors cols <- c("Both"="#0570b0","PIPO"="#74a9cf","PSME"="#bdc9e1","CA"="#CC79A7","CO"="#000000","NR"="#009E73","SW"="#E69F00")#new ones # plot recruitment probability ggplot(nd.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability\n")+xlab("Year")+ #geom_point(data=nd.reg.cur,aes(year,pulse,color=Region))+ theme_classic()+ theme(legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5)) ## get sd of annual recruitment probability by region nd.reg.cur<-nd.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.reg.cur<-data.frame(nd.reg.cur) ggplot(nd.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_point(shape=1)+ geom_smooth(aes(fill=region),se=F,size=0.4)+ scale_colour_manual(values = cols)+ xlab("Year")+ ylab("SD of recruitment probability")+ theme(#legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) ######################################################## ### predict to get cumulative recruitment probability 5-10 yrs post-fire pred.brt10<-function(sites,dat,nd.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=seq(1,10,1),dist_PIPO_seed_source=50, dnbr=412,year=dat.test$year) nd$vwc_mean_driestM<-rep(dat.test$vwc_mean_driestM,each=10) nd$vpd_mean678<-rep(dat.test$vpd_mean678,each=10) # current vpd nd$region<-rep(dat.test$region,each=10) nd$pr<-predict.gbm(m.pi, nd, n.trees=m.pi$gbm.call$best.trees, type="response") nd$site<-i nd.sum<-rbind(nd.sum,nd) } assign("nd.510.pipo", nd.sum, envir = .GlobalEnv) } nd.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_PIPO_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.pi$site) pred.brt10(sites,dat.c.pi,nd.sum) head(nd.510.pipo) # now calculate pr5 and pr10 (cumulative probability in first 5 or 10 years postfire) for each site nd.510.pipo2<- nd.510.pipo %>% group_by(region,year,site) %>% summarize(pr=sum(pr)) nd.510.pipo2<-data.frame(nd.510.pipo2) head(nd.510.pipo2) head(nd.510.pipo);head(nd.510.pipo2) nd.510.pipo2$pr5<-NA nd.510.pipo2$pr10<-NA nd.510.pipo2$FireYear<-NA for (i in 1:nrow(nd.510.pipo2)){ date1 <- nd.510.pipo2[i, "year"];region1 <- nd.510.pipo2[i, "region"];site<-nd.510.pipo2[i,"site"] nd.510.pipo2[i,"FireYear"]<-date1-1 if (date1<2012){ getrow1<-which(nd.510.pipo$year==date1 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==1 & nd.510.pipo$site==site) getrow2<-which(nd.510.pipo$year==date1+1 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==2 & nd.510.pipo$site==site) getrow3<-which(nd.510.pipo$year==date1+2 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==3 & nd.510.pipo$site==site) getrow4<-which(nd.510.pipo$year==date1+3 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==4 & nd.510.pipo$site==site) getrow5<-which(nd.510.pipo$year==date1+4 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==5 & nd.510.pipo$site==site) nd.510.pipo2[i,"pr5"]<-nd.510.pipo[getrow1,"pr"]+nd.510.pipo[getrow2,"pr"]+nd.510.pipo[getrow3,"pr"]+nd.510.pipo[getrow4,"pr"]+nd.510.pipo[getrow5,"pr"] } if (date1<2007){ getrow6<-which(nd.510.pipo$year==date1+5 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==6 & nd.510.pipo$site==site) getrow7<-which(nd.510.pipo$year==date1+6 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==7 & nd.510.pipo$site==site) getrow8<-which(nd.510.pipo$year==date1+7 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==8 & nd.510.pipo$site==site) getrow9<-which(nd.510.pipo$year==date1+8 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==9 & nd.510.pipo$site==site) getrow10<-which(nd.510.pipo$year==date1+9 & nd.510.pipo$region==region1 & nd.510.pipo$time_since_fire==10 & nd.510.pipo$site==site) nd.510.pipo2[i,"pr10"]<-nd.510.pipo[getrow1,"pr"]+nd.510.pipo[getrow2,"pr"]+nd.510.pipo[getrow3,"pr"]+nd.510.pipo[getrow4,"pr"]+nd.510.pipo[getrow5,"pr"]+ nd.510.pipo[getrow6,"pr"]+nd.510.pipo[getrow7,"pr"]+nd.510.pipo[getrow8,"pr"]+nd.510.pipo[getrow9,"pr"]+nd.510.pipo[getrow10,"pr"] } } head(nd.510.pipo2);tail(nd.510.pipo2) # get mean by region nd.510.pipo3<- nd.510.pipo2 %>% group_by(region,FireYear) %>% summarize(pr5=mean(pr5),pr10=mean(pr10)) nd.510.pipo3<-data.frame(nd.510.pipo3) ## plot pr5 & pr10 head(nd.510.pipo3) cols <- c("Both"="#0570b0","PIPO"="#74a9cf","PSME"="#bdc9e1","CA"="#CC79A7","CO"="#000000","NR"="#009E73","SW"="#E69F00") pipo.rp5<-ggplot(nd.510.pipo3,aes(FireYear,pr5,color=region))+geom_line()+ ylab("Cumulative recruitment probability\nyears 1-5 post-fire")+xlab("Fire year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position=c(0.97,1),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2010,5),limits=c(1980,2010)) pipo.rp5 pipo.rp10<-ggplot(nd.510.pipo3,aes(FireYear,pr10,color=region))+geom_line()+ ylab("Cumulative recruitment probability\nyears 1-10 post-fire")+xlab("Fire year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position="none", #legend.position=c(0.97,1.15),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2010,5),limits=c(1980,2010)) pipo.rp10 ########################################################################### ## use change point analysis to detect shifts in time series of pr, pr5, and pr10 ## results in Table S5 ################################################## ## test for breakpoints in pr 1 yr postfire head(nd.reg.cur) nd.reg.cur$FireYear<-nd.reg.cur$year-1 # extract time series for just CA d.ca<-ts (nd.reg.cur[nd.reg.cur$region=="CA",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) # extract time series for just CO d.co<-ts (nd.reg.cur[nd.reg.cur$region=="CO",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.co) fs.vpd.co <- Fstats(d.co ~ 1) plot(fs.vpd.co) sctest(fs.vpd.co) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.co) lines(breakpoints(fs.vpd.co)) breakpoints(fs.vpd.co) # extract time series for just NR d.nr<-ts (nd.reg.cur[nd.reg.cur$region=="NR",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) # extract time series for just SW d.sw<-ts (nd.reg.cur[nd.reg.cur$region=="SW",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ##################################### ## test for breakpoints in pr5 head(nd.510.pipo3) ## make timeseries for CA d.ca<-ts (nd.510.pipo3[nd.510.pipo3$region=="CA",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) plot(fs.vpd.ca,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) ## make timeseries for CO d.co<-ts (nd.510.pipo3[nd.510.pipo3$region=="CO",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.co) fs.vpd.co <- Fstats(d.co ~ 1) plot(fs.vpd.co) plot(fs.vpd.co,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.co) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.co) lines(breakpoints(fs.vpd.co)) breakpoints(fs.vpd.co) ## make timeseries for NR d.nr<-ts (nd.510.pipo3[nd.510.pipo3$region=="NR",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) plot(fs.vpd.nr,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) ## make timeseries for SW d.sw<-ts (nd.510.pipo3[nd.510.pipo3$region=="SW",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) plot(fs.vpd.sw,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ########################################## ## test for breakpoints in pr10 #### make time series for CA d.ca<-ts (nd.510.pipo3[nd.510.pipo3$region=="CA",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) plot(fs.vpd.ca,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) #### make time series for CO d.co<-ts (nd.510.pipo3[nd.510.pipo3$region=="CO",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.co) fs.vpd.co <- Fstats(d.co ~ 1) plot(fs.vpd.co) plot(fs.vpd.co,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.co) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.co) lines(breakpoints(fs.vpd.co)) breakpoints(fs.vpd.co) #### make time series for NR d.nr<-ts (nd.510.pipo3[nd.510.pipo3$region=="NR",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) plot(fs.vpd.nr,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) #### make time series for SW d.sw<-ts (nd.510.pipo3[nd.510.pipo3$region=="SW",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) plot(fs.vpd.sw,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ################################################################################################## ########################################## ###### Repeat all for PSME ############################# ################################################## datps<-dat[dat$species=="PSME" & dat$region!="CO",] datps$site<-factor(datps$site) datps$fire<-factor(datps$fire) datps$region<-factor(datps$region) datps$cv_group<-factor(datps$cv_group) datps<-transform(datps, cv.num=match(cv_group, unique(cv_group))) ## calculate and add annual recruitment rate threshold for each region qc<-function(x){quantile(x,c(0.25,0.5))} tapply(datps[datps$density>0,]$density,datps[datps$density>0,]$region,qc) datps$q25<-ifelse(datps$region=="CA",0.0017,ifelse(datps$region=="NR",0.0056,0.0012)) # create binomial response based on regional annual recruitment rate thresholds datps$PA25<-ifelse(datps$density>=datps$q25,1,0) names(datps) #vars<-c(7,20,28,32,33,34,35,38)# all variables vars<-c(7,20,28,34,32) # simplified model final ## binomial m.ps<-gbm.step(data=datps, gbm.x = vars, gbm.y = 41, # response is "PA25" family = "bernoulli", tree.complexity = 3, learning.rate = 0.001, bag.fraction = 0.75, fold.vector=datps[,39],n.folds=length(unique(datps[,39]))) #fold vector is "cv.num" ## plot partial dependency plots & see relative influence gbm.plot(m.ps, write.title = FALSE) ## Model AUC datps$pr<-predict.gbm(m.ps, datps, n.trees=m.ps$gbm.call$best.trees, type="response") roccurve <- roc(datps$PA25 ~ datps$pr) coords(roccurve, "b", ret="t", best.method="youden") plot(roccurve) auc(roccurve) ############################################################## ##### Validate model by holding out each site and calculating accuracy & AUC # create a function to hold out each site, fit the model, predict to the holdpout site, and caluclate accuracy and AUC holdout1.site<-function(sites,dat,metrics){ for (i in sites){ test.sites<- i metrics[i,5]<-test.sites dat.test<-dat[dat$site==test.sites,] dat.train<-dat[dat$site!=test.sites,] metrics[i,6]<-paste(dat.test[1,"region"]) gbm1 <- gbm(PA25~time_since_fire+dist_PSME_seed_source+dnbr+vwc_mean345+tskin_max, data=dat.train, distribution="bernoulli", n.trees=trees, shrinkage=lr,interaction.depth=id, bag.fraction = bf) actual <- dat.test$PA25 preds.g <- predict.gbm(gbm1, dat.test, n.trees=trees, type="response") if(sum(dat.test$PA25)>0 & sum(dat.test$PA25)/length(dat.test$PA25)!=1){ roccurve <- roc(dat.test$PA25 ~ preds.g) thresh=coords(roccurve, "b", ret="t", best.method="youden") thresh=thresh[1] metrics[i,4]<-thresh results <- HMeasure(actual,preds.g) metrics[i,2]<-results$metrics['H'] metrics[i,3]<-auc(roccurve) } else { metrics[i,2:3]<-NA thresh=0.23 metrics[i,4]<-thresh } pr <- as.numeric(predict.gbm(gbm1, dat.test, n.trees=trees, type="response")>thresh) #actual;pr;preds.g;dat.test$proportion cm = as.matrix(table(Actual = actual, Predicted = pr)) n = sum(cm) # number of instances diag = diag(cm) # number of correctly classified instances per class if (sum(actual)==0 & sum(pr)==length(pr) | sum(actual)==length(actual) & sum(pr)==0){ accuracy=0 k=NA } else{ accuracy = sum(diag) / n dat.test$pr<-pr } metrics[i,1]<-accuracy } assign("metrics.sum.ps", metrics, envir = .GlobalEnv) } # set settings from above gbm.step model bf=0.75 lr=0.001 trees=1400 id=3 # get names of PSME sites sites<-levels(datps$site) # make data frame for output metrics metrics<-setNames(data.frame(matrix(ncol =6, nrow = 0)), c("acc","H","AUC","thresh","site","region")) # run function holdout1.site(sites,datps,metrics) metrics.sum.ps colMeans(metrics.sum.ps[c(1:4)],na.rm=T) # get means by region and test for differences # accuracy metrics.sum.ps$region<-as.factor(metrics.sum.ps$region) tapply(metrics.sum.ps$acc,metrics.sum.ps$region,mean) an<-aov(metrics.sum.ps$acc~metrics.sum.ps$region) summary(an) TukeyHSD(an) #AUC tapply(metrics.sum.ps$AUC,metrics.sum.ps$region,mean,na.rm=T) an<-aov(metrics.sum.ps$AUC~metrics.sum.ps$region) summary(an) TukeyHSD(an) # H measure tapply(metrics.sum.ps$H,metrics.sum.ps$region,mean,na.rm=T) an<-aov(metrics.sum.ps$H~metrics.sum.ps$region) summary(an) TukeyHSD(an) ################################################### ####### predict recruitment probability first year after fire ## add region to dat.c from dat.ps dataframe, then get rid of sites without region # which is a round about way of selecting only sites used for PSME model in the dat.c dataframe to then used for hindcasting dat.c.ps<-dat.c nm <- c("region") dat.c.ps[nm] <- lapply(nm, function(x) datps[[x]][match(dat.c.ps$site, datps$site)]) dat.c.ps<-dat.c.ps[is.na(dat.c.ps$region)==F,] #take out sites not used in psme analysis dat.c.ps$site<-factor(dat.c.ps$site) # make regional means for dat.c.ps for plotting purposes dat.c.ps.reg<- dat.c.ps %>% group_by(region,year) %>% summarize_at(c("tskin_max","vpd_mean678","vwc_mean345","vwc_mean_driestM"), mean) dat.c.ps.reg<-data.frame(dat.c.ps.reg) ## make function to predict all sites based on actual annual climate with fixed severity and DSS ### pred.brt<-function(sites,dat,nd.ps.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PSME_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean345<-rep(dat.test$vwc_mean345) nd$tskin_max<-rep(dat.test$tskin_max) nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.ps, nd, n.trees=m.ps$gbm.call$best.trees, type="response") nd$site<-i nd.ps.sum<-rbind(nd.ps.sum,nd) } assign("nd.ps.fin.cur", nd.ps.sum, envir = .GlobalEnv) } nd.ps.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_psme_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.ps$site) pred.brt(sites,dat.c.ps,nd.ps.sum) head(nd.ps.fin.cur) ### summarize by region ### nd.ps.reg.cur<- ddply(nd.ps.fin.cur,.(region,year),summarize,pr=mean(pr)) # plot recruitment probability ggplot(nd.ps.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability\n")+xlab("Year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5)) ## get sd of annual recruitment probability by region nd.ps.reg.cur<-nd.ps.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.ps.reg.cur<-data.frame(nd.ps.reg.cur) ggplot(nd.ps.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_point(shape=1)+ geom_smooth(aes(fill=region),se=F,size=0.4)+ xlab("Year")+ ylab("SD of recruitment probability")+ scale_colour_manual(values = cols)+ theme(#legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) ######################################################## ### predict to get cumulative recruitment probability 5-10 yrs post-fire pred.brt<-function(sites,dat,nd.ps.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=seq(1,10,1),dist_PSME_seed_source=50, dnbr=412,year=dat.test$year) nd$vwc_mean345<-rep(dat.test$vwc_mean345,each=10) nd$tskin_max<-rep(dat.test$tskin_max,each=10) nd$region<-rep(dat.test$region,each=10) nd$pr<-predict.gbm(m.ps, nd, n.trees=m.ps$gbm.call$best.trees, type="response") nd$site<-i nd.ps.sum<-rbind(nd.ps.sum,nd) } assign("nd.510.psme", nd.ps.sum, envir = .GlobalEnv) } nd.ps.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_psme_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.ps$site) pred.brt(sites,dat.c.ps,nd.ps.sum) head(nd.510.psme) # now calculate cumulative probability of recruitment 5 (pr5) and 10 (pr10) years postfire nd.510.psme2<- nd.510.psme %>% group_by(region,year,site) %>% summarize(pr=sum(pr)) nd.510.psme2<-data.frame(nd.510.psme2) head(nd.510.psme2) head(nd.510.psme);head(nd.510.psme2) nd.510.psme2$pr5<-NA nd.510.psme2$pr10<-NA nd.510.psme2$FireYear<-NA for (i in 1:nrow(nd.510.psme2)){ date1 <- nd.510.psme2[i, "year"];region1 <- nd.510.psme2[i, "region"];site<-nd.510.psme2[i,"site"] nd.510.psme2[i,"FireYear"]<-date1-1 if (date1<2012){ # if just want pr 5 date<2012, if want pr10 need date<2007 getrow1<-which(nd.510.psme$year==date1 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==1 & nd.510.psme$site==site) getrow2<-which(nd.510.psme$year==date1+1 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==2 & nd.510.psme$site==site) getrow3<-which(nd.510.psme$year==date1+2 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==3 & nd.510.psme$site==site) getrow4<-which(nd.510.psme$year==date1+3 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==4 & nd.510.psme$site==site) getrow5<-which(nd.510.psme$year==date1+4 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==5 & nd.510.psme$site==site) nd.510.psme2[i,"pr5"]<-nd.510.psme[getrow1,"pr"]+nd.510.psme[getrow2,"pr"]+nd.510.psme[getrow3,"pr"]+nd.510.psme[getrow4,"pr"]+nd.510.psme[getrow5,"pr"] } if (date1<2007){ getrow6<-which(nd.510.psme$year==date1+5 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==6 & nd.510.psme$site==site) getrow7<-which(nd.510.psme$year==date1+6 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==7 & nd.510.psme$site==site) getrow8<-which(nd.510.psme$year==date1+7 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==8 & nd.510.psme$site==site) getrow9<-which(nd.510.psme$year==date1+8 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==9 & nd.510.psme$site==site) getrow10<-which(nd.510.psme$year==date1+9 & nd.510.psme$region==region1 & nd.510.psme$time_since_fire==10 & nd.510.psme$site==site) nd.510.psme2[i,"pr10"]<-nd.510.psme[getrow1,"pr"]+nd.510.psme[getrow2,"pr"]+nd.510.psme[getrow3,"pr"]+nd.510.psme[getrow4,"pr"]+nd.510.psme[getrow5,"pr"]+ nd.510.psme[getrow6,"pr"]+nd.510.psme[getrow7,"pr"]+nd.510.psme[getrow8,"pr"]+nd.510.psme[getrow9,"pr"]+nd.510.psme[getrow10,"pr"] } } head(nd.510.psme2);tail(nd.510.psme2) # get regional means nd.510.psme3<- nd.510.psme2 %>% group_by(region,FireYear) %>% summarize(pr5=mean(pr5),pr10=mean(pr10)) nd.510.psme3<-data.frame(nd.510.psme3) ## plot pr5 & pr10 head(nd.510.psme3) cols <- c("Both"="#0570b0","PIPO"="#74a9cf","PSME"="#bdc9e1","CA"="#CC79A7","CO"="#000000","NR"="#009E73","SW"="#E69F00")#new ones psme.rp5<-ggplot(nd.510.psme3,aes(FireYear,pr5,color=region))+geom_line()+ ylab("Cumulative recruitment probability\nyears 1-5 post-fire")+xlab("Fire year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position=c(0.97,1),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2010,5),limits=c(1980,2010)) psme.rp5 psme.rp10<-ggplot(nd.510.psme3,aes(FireYear,pr10,color=region))+geom_line()+ ylab("Cumulative recruitment probability\nyears 1-10 post-fire")+xlab("Fire year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position="none", #legend.position=c(0.97,1.15),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2010,5),limits=c(1980,2010)) psme.rp10 ########################################################################### ## use change point analysis to detect shifts in time series of pr, pr5, and pr10 ## results in Table S5 ## test for breakpoints in pr 1 yr postfire head(nd.ps.reg.cur) nd.ps.reg.cur$FireYear<-nd.ps.reg.cur$year-1 #make time series for CA d.ca<-ts (nd.ps.reg.cur[nd.ps.reg.cur$region=="CA",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) #make time series for NR d.nr<-ts (nd.ps.reg.cur[nd.ps.reg.cur$region=="NR",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) plot(fs.vpd.nr,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) #make time series for SW d.sw<-ts (nd.ps.reg.cur[nd.ps.reg.cur$region=="SW",]$pr, start=c(1980), end=c(2014), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) plot(fs.vpd.sw,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ############################################# ## test for breakpoints in pr5 head(nd.510.psme3) #create CA timeseries d.ca<-ts (nd.510.psme3[nd.510.psme3$region=="CA",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) plot(fs.vpd.ca,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) #create NR timeseries d.nr<-ts (nd.510.psme3[nd.510.psme3$region=="NR",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) plot(fs.vpd.nr,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) #create SW timeseries d.sw<-ts (nd.510.psme3[nd.510.psme3$region=="SW",]$pr5, start=c(1980), end=c(2010), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) plot(fs.vpd.sw,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ############################ ## test for breakpoints in pr10 # make CA timeseries d.ca<-ts (nd.510.psme3[nd.510.psme3$region=="CA",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.ca) fs.vpd.ca <- Fstats(d.ca ~ 1) plot(fs.vpd.ca) plot(fs.vpd.ca,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.ca) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.ca) lines(breakpoints(fs.vpd.ca)) breakpoints(fs.vpd.ca) # make NR timeseries d.nr<-ts (nd.510.psme3[nd.510.psme3$region=="NR",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.nr) fs.vpd.nr <- Fstats(d.nr ~ 1) plot(fs.vpd.nr) plot(fs.vpd.nr,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.nr) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.nr) lines(breakpoints(fs.vpd.nr)) breakpoints(fs.vpd.nr) # make SW timeseries d.sw<-ts (nd.510.psme3[nd.510.psme3$region=="SW",]$pr10, start=c(1980), end=c(2005), frequency=1) plot(d.sw) fs.vpd.sw <- Fstats(d.sw ~ 1) plot(fs.vpd.sw) plot(fs.vpd.sw,pval=T,ylim=c(0,0.05)) sctest(fs.vpd.sw) ## visualize the breakpoint implied by the argmax of the F statistics plot(d.sw) lines(breakpoints(fs.vpd.sw)) breakpoints(fs.vpd.sw) ######################################################################## ########### Eveness of observed recruitment over time ### #################################################### even<-function(prop,tsf){H<-diversity(prop) H/log(max(tsf))} ########################### ## make evenness scatterplots dat<-read.csv("Davis_et_al_recruitment_plot_data.csv") dat<-dat[dat$time_since_fire>0,] dat<-dat[dat$year<2016,] fn <- function(x) x/sum(x) dat<-ddply(dat,.(site,species),mutate,Prop=fn(count)) dt<-ddply(dat[dat$species=="PIPO" & dat$PIPO_taken>9 | dat$species=="PSME" & dat$PSME_taken>9,],.(site,species),mutate,even=even(Prop,time_since_fire)) head(dt) dat.site<-ddply(dt,.(site,region,species,PIPO_taken,PSME_taken),summarize,shrub=mean(mean_shrub_cover), ev=mean(even),def30=mean(def30),dnbr=mean(dnbr)) dat.site<-dat.site[dat.site$region!="CO",] head(dat.site) dat.vars<-dat.site[,c(3,6,7,8,9)] ev1a<-lm(ev~.*.,data=dat.vars) Anova(ev1a) #def and species*shrub ev1<-lm(ev~shrub*species+def30,data=dat.site) summary(ev1) ## even vs deficit nd<-expand.grid(species=c("PIPO","PSME"),shrub=c(20,80),def30=seq(300,1000,1)) ndd <- cbind(nd, predict(ev1,nd,se=TRUE)) ndd <- within(ndd, { LL <- fit - (1.96 * se.fit) UL <- fit + (1.96 * se.fit) }) labels <- c(PIPO = "Ponderosa pine", PSME = "Douglas-fir") evenD<-ggplot()+geom_point(data=dat.site,aes(x=def30,y=ev))+ facet_wrap(~species,labeller=labeller(species = labels))+ geom_line(data=ndd,aes(def30,fit,linetype=as.factor(shrub)))+ geom_ribbon(data=ndd,aes(x=def30,ymin = LL,ymax = UL,group=as.factor(shrub)),alpha = .25)+ theme_bw()+ scale_linetype_discrete("Shrub\ncover (%)")+ ylab("Evenness of recruitment (J)")+ xlab("30-year mean deficit (mm)")+ theme(panel.grid.minor=element_blank(),panel.grid.major=element_blank(), axis.text.x = element_text(size=10,color=1), axis.text.y = element_text(size=10,color=1), axis.title = element_text(size=11), legend.text= element_text(size=10), legend.title= element_text(size=10), strip.text.x = element_text(size=10), strip.background = element_rect(fill="white"), legend.key.height=unit(0.47,"cm"), legend.justification=c(1,1),legend.position=c(1,1), legend.background = element_rect(fill = "transparent", colour = "transparent")) evenD ######################################################################### ## plot partial dependency plots ######################### ## add degree function add_degree <- function(x) { paste(x, "(ÂșC)",sep = "") } ## 25% threshold model (main model in text) vars<-c(7,19,28,33,35) # simplified final model ## binomial m.pi<-gbm.step(data=datpi, gbm.x = vars, gbm.y = 41, family = "bernoulli", tree.complexity = 3, learning.rate = 0.002, bag.fraction = 0.75, fold.vector=datpi[,39],n.folds=length(unique(datpi[,39]))) gbm.plot(m.pi) plotdata.tsf <- plot.gbm(m.pi, i.var = c(1), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") tsf<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.4, label="22%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.13,0.45))+ geom_rug(data=datpi,aes(time_since_fire,PA25),sides="b") plotdata.dss <- plot.gbm(m.pi, i.var = c(2), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dss<-ggplot(plotdata.dss,aes(dist_PIPO_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(88,0.4, label="32%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.13,0.45))+ geom_rug(data=datpi,aes(dist_PIPO_seed_source,PA25),sides="b") plotdata.vwc <- plot.gbm(m.pi, i.var = c(4), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vwc<-ggplot(plotdata.vwc,aes(vwc_mean_driestM,y))+geom_line()+ ylab("Probability")+ xlab("Soil moist. driest month (vwc)")+ geom_text(aes(0.04,0.4, label="13%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.13,0.45))+ geom_rug(data=datpi,aes(vwc_mean_driestM,PA25),sides="b") vwc plotdata.vpd<- plot.gbm(m.pi, i.var = c(5), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vpd<-ggplot(plotdata.vpd,aes(vpd_mean678,y))+geom_line()+ ylab("Probability")+ xlab("Mean summer VPD (z)")+ geom_text(aes(1.8,0.4, label="11%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.13,0.45))+ geom_rug(data=datpi,aes(vpd_mean678,PA25),sides="b") vpd plotdata.dnbr<- plot.gbm(m.pi, i.var = c(3), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(380,0.4, label="22%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.13,0.45))+ geom_rug(data=datpi,aes(dnbr,PA25),sides="b") dnbr pipo.pd<-plot_grid(dss,tsf,dnbr,vwc,vpd,nrow=2,vjust = 1.2) pipo.pd ## psme names(datps) vars<-c(7,20,28,34,32) # simplified model final ## binomial m.ps<-gbm.step(data=datps, gbm.x = vars, gbm.y = 41, family = "bernoulli", tree.complexity = 3, learning.rate = 0.001, bag.fraction = 0.75, fold.vector=datps[,39],n.folds=length(unique(datps[,39]))) summary(m.ps) gbm.plot(m.ps, write.title = FALSE) plotdata.tsf <- plot.gbm(m.ps, i.var = c(1), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tsf.ps<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.32, label="35%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.16,0.34))+ geom_rug(data=datps,aes(time_since_fire,PA25),sides="b") plotdata.dss <- plot.gbm(m.ps, i.var = c(2), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dss.ps<-ggplot(plotdata.dss,aes(dist_PSME_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(85,0.32, label="21%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.16,0.34))+ geom_rug(data=datps,aes(dist_PSME_seed_source,PA25),sides="b") plotdata.vwc <- plot.gbm(m.ps, i.var = c(4), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") vwc.ps<-ggplot(plotdata.vwc,aes(vwc_mean345,y))+geom_line()+ ylab("Probability")+ xlab("Spring soil moist. (vwc)")+ geom_text(aes(0.33,0.32, label="20%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.16,0.34))+ geom_rug(data=datps,aes(vwc_mean345,PA25),sides="b") plotdata.tskin<- plot.gbm(m.ps, i.var = c(5), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tskin.ps<-ggplot(plotdata.tskin,aes(tskin_max,y))+geom_line()+ ylab("Probability")+ xlab(add_degree("Max. surface temp. "))+ geom_text(aes(58,0.32, label="14%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.16,0.34))+ geom_rug(data=datps,aes(tskin_max,PA25),sides="b") tskin.ps plotdata.dnbr<- plot.gbm(m.ps, i.var = c(3), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr.ps<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(420,0.32, label="10%"),size=3.4,check_overlap = TRUE)+ #8% for dnbr.norm theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.16,0.34))+ geom_rug(data=datps,aes(dnbr,PA25),sides="b") psme.pd<-plot_grid(tsf.ps,dss.ps,vwc.ps,tskin.ps,dnbr.ps,nrow=2,vjust = 1.2) psme.pd ## plot partial dependency plots both species ### plot_grid(pipo.pd,psme.pd,nrow=2,labels="AUTO",vjust=1.1) ################################## ## Presence absence only models ## pipo names(datpi) vars<-c(7,19,28,33,35) # simplified final model ## binomial m.pi<-gbm.step(data=datpi, gbm.x = vars, gbm.y = 31, family = "bernoulli", tree.complexity = 3, learning.rate = 0.002, bag.fraction = 0.75, fold.vector=datpi[,39],n.folds=length(unique(datpi[,39]))) gbm.plot(m.pi) ## plot partial dependency plots plotdata.tsf <- plot.gbm(m.pi, i.var = c(1), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") tsf<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.49, label="34%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.21,0.52))+ geom_rug(data=datpi,aes(time_since_fire,PA),sides="b") plotdata.dss <- plot.gbm(m.pi, i.var = c(2), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dss<-ggplot(plotdata.dss,aes(dist_PIPO_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(88,0.49, label="30%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.21,0.52))+ geom_rug(data=datpi,aes(dist_PIPO_seed_source,PA),sides="b") plotdata.vwc <- plot.gbm(m.pi, i.var = c(4), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vwc<-ggplot(plotdata.vwc,aes(vwc_mean_driestM,y))+geom_line()+ ylab("Probability")+ xlab("Soil moist. driest month (vwc)")+ geom_text(aes(0.085,0.47, label="9%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.21,0.52))+ #scale for PAT5 geom_rug(data=datpi,aes(vwc_mean_driestM,PA),sides="b") plotdata.vpd<- plot.gbm(m.pi, i.var = c(5), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vpd<-ggplot(plotdata.vpd,aes(vpd_mean678,y))+geom_line()+ ylab("Probability")+ xlab("Mean summer VPD (z)")+ geom_text(aes(1.8,0.49, label="5%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.21,0.52))+ geom_rug(data=datpi,aes(vpd_mean678,PA),sides="b") vpd plotdata.dnbr<- plot.gbm(m.pi, i.var = c(3), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(800,0.49, label="22%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.21,0.52))+ geom_rug(data=datpi,aes(dnbr,PA),sides="b") pipo.pd<-plot_grid(tsf,dss,dnbr,vwc,vpd,nrow=2,vjust = 1.2) pipo.pd ## psme names(datps) vars<-c(7,20,28,34,32) # simplified model final ## binomial m.ps<-gbm.step(data=datps, gbm.x = vars, gbm.y = 31, family = "bernoulli", tree.complexity = 3, learning.rate = 0.001, bag.fraction = 0.75, fold.vector=datps[,39],n.folds=length(unique(datps[,39]))) gbm.plot(m.ps, write.title = FALSE) summary(m.ps) plotdata.tsf <- plot.gbm(m.ps, i.var = c(1), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tsf.ps<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.4, label="30%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.20,0.42))+ geom_rug(data=datps,aes(time_since_fire,PA),sides="b") plotdata.dss <- plot.gbm(m.ps, i.var = c(2), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dss.ps<-ggplot(plotdata.dss,aes(dist_PSME_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(85,0.4, label="19%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.20,0.42))+ geom_rug(data=datps,aes(dist_PSME_seed_source,PA),sides="b") plotdata.vwc <- plot.gbm(m.ps, i.var = c(4), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") vwc.ps<-ggplot(plotdata.vwc,aes(vwc_mean345,y))+geom_line()+ ylab("Probability")+ xlab("Spring soil moist. (vwc)")+ geom_text(aes(0.33,0.4, label="23%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.20,0.42))+ geom_rug(data=datps,aes(vwc_mean345,PA),sides="b") plotdata.tskin<- plot.gbm(m.ps, i.var = c(5), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tskin.ps<-ggplot(plotdata.tskin,aes(tskin_max,y))+geom_line()+ ylab("Probability")+ xlab(add_degree("Max. skin temp. "))+ geom_text(aes(58,0.4, label="13%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.20,0.42))+ geom_rug(data=datps,aes(tskin_max,PA),sides="b") tskin.ps plotdata.dnbr<- plot.gbm(m.ps, i.var = c(3), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr.ps<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(800,0.4, label="14%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.20,0.42))+ geom_rug(data=datps,aes(dnbr,PA),sides="b") psme.pd<-plot_grid(tsf.ps,vwc.ps,dss.ps,dnbr.ps,tskin.ps,nrow=2,vjust = 1.2) psme.pd ## plot partial dependency plots both species ### plot_grid(pipo.pd,psme.pd,nrow=2,labels="AUTO",vjust=1.1) ########################################################### ## predict for 1 yr postfire and constant dss, severity ####################### pred.brt<-function(sites,dat,nd.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PIPO_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean_driestM<-rep(dat.test$vwc_mean_driestM) nd$vpd_mean678<-rep(dat.test$vpd_mean678) nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.pi, nd, n.trees=m.pi$gbm.call$best.trees, type="response") nd$site<-i nd.sum<-rbind(nd.sum,nd) } assign("nd.fin.cur", nd.sum, envir = .GlobalEnv) } nd.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_PIPO_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.pi$site) pred.brt(sites,dat.c.pi,nd.sum) ### summarize by region ### nd.reg.cur<- nd.fin.cur %>% group_by(region,year) %>% summarize(pr=mean(pr)) nd.reg.cur<-data.frame(nd.reg.cur) head(nd.reg.cur) ## get sd of annual recruitment probability by region nd.reg.cur<-nd.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.reg.cur<-data.frame(nd.reg.cur) #psme pred.brt<-function(sites,dat,nd.ps.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PSME_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean345<-rep(dat.test$vwc_mean345) nd$tskin_max<-rep(dat.test$tskin_max) nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.ps, nd, n.trees=m.ps$gbm.call$best.trees, type="response") nd$Site<-i nd.ps.sum<-rbind(nd.ps.sum,nd) } assign("nd.ps.fin.cur", nd.ps.sum, envir = .GlobalEnv) } nd.ps.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_psme_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.ps$site) pred.brt(sites,dat.c.ps,nd.ps.sum) ### summarize by region ### nd.ps.reg.cur<- ddply(nd.ps.fin.cur,.(region,year),summarize,pr=mean(pr)) ## get sd of annual recruitment probability by region nd.ps.reg.cur<-nd.ps.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.ps.reg.cur<-data.frame(nd.ps.reg.cur) ### RP plots cols <- c("Both"="#0570b0","PIPO"="#74a9cf","PSME"="#bdc9e1","CA"="#CC79A7","CO"="#000000","NR"="#009E73","SW"="#E69F00")#new ones pipo.rp<-ggplot(nd.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability")+xlab("Year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position=c(0.97,1.1),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5)) pipo.rp psme.rp<-ggplot(nd.ps.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability")+xlab("Year")+ theme_classic()+ theme(#legend.position=c(0.97,1.15),legend.justification=c(1,1), legend.position="none", legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_colour_manual(name = "Region",values = cols)+ scale_x_continuous(breaks=seq(1980,2015,5)) psme.rp ## SD of RP pipo.sd<-ggplot(nd.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_smooth(aes(fill=region),se=F)+ scale_color_manual(values=cols)+ scale_fill_manual(name = "Region",values = cols)+ xlab("Year")+ ylab("SD of recruitment probability")+ theme(legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) pipo.sd psme.sd<-ggplot(nd.ps.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_smooth(aes(fill=region),se=F)+ scale_color_manual(values=cols)+ scale_fill_manual(name = "Region",values = cols)+ xlab("Year")+ ylab("SD of recruitment probability")+ theme(legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) psme.sd plot_grid(NULL,pipo.sd,pipo.rp,NULL,psme.sd,psme.rp,rel_widths=c(0.08,0.6,1),labels=c("A","","","B","","")) ################################## ## 50% threshold models ## pipo datpi$q50<-ifelse(datpi$region=="CO",0.0080,ifelse(datpi$region=="NR",0.0071,ifelse(datpi$region=="CA",0.0053,0.0025))) datpi$PA50<-ifelse(datpi$density>=datpi$q50,1,0) names(datpi) vars<-c(7,19,28,33,35) # simplified final model ## binomial m.pi<-gbm.step(data=datpi, gbm.x = vars, gbm.y = 44, #response is "PA50" family = "bernoulli", tree.complexity = 3, learning.rate = 0.002, bag.fraction = 0.75, fold.vector=datpi[,39],n.folds=length(unique(datpi[,39]))) #fol vector is "cv.num" gbm.plot(m.pi) summary(m.pi) ## plot partial dependency plots plotdata.tsf <- plot.gbm(m.pi, i.var = c(1), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") tsf<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.43, label="21%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.1,0.465))+ geom_rug(data=datpi,aes(time_since_fire,PA50),sides="b") plotdata.dss <- plot.gbm(m.pi, i.var = c(2), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dss<-ggplot(plotdata.dss,aes(dist_PIPO_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(88,0.43, label="24%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.1,0.465))+ geom_rug(data=datpi,aes(dist_PIPO_seed_source,PA50),sides="b") plotdata.vwc <- plot.gbm(m.pi, i.var = c(4), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vwc<-ggplot(plotdata.vwc,aes(vwc_mean_driestM,y))+geom_line()+ ylab("Probability")+ xlab("Soil moist. driest month (vwc)")+ geom_text(aes(0.085,0.43, label="9%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.1,0.465))+ geom_rug(data=datpi,aes(vwc_mean_driestM,PA),sides="b") plotdata.vpd<- plot.gbm(m.pi, i.var = c(5), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") vpd<-ggplot(plotdata.vpd,aes(vpd_mean678,y))+geom_line()+ ylab("Probability")+ xlab("Mean summer VPD (z)")+ geom_text(aes(1.8,0.43, label="17%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.1,0.465))+ geom_rug(data=datpi,aes(vpd_mean678,PA50),sides="b") vpd plotdata.dnbr<- plot.gbm(m.pi, i.var = c(3), n.trees = m.pi$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(800,0.43, label="29%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.1,0.465))+ geom_rug(data=datpi,aes(dnbr,PA50),sides="b") pipo.pd<-plot_grid(tsf,dss,dnbr,vwc,vpd,nrow=2,vjust = 1.2) pipo.pd ## psme datps$q50<-ifelse(datps$region=="CA",0.0042,ifelse(datps$region=="NR",0.0129,0.0024)) datps$PA50<-ifelse(datps$density>=datps$q50,1,0) names(datps) vars<-c(7,20,28,32,34) # simplified model final ## binomial m.ps<-gbm.step(data=datps, gbm.x = vars, gbm.y = 44, #response is "PA50" family = "bernoulli", tree.complexity = 3, learning.rate = 0.001, bag.fraction = 0.75, fold.vector=datps[,39],n.folds=length(unique(datps[,39]))) #fold vector is "cv.num" gbm.plot(m.ps, write.title = FALSE) summary(m.ps) plotdata.tsf <- plot.gbm(m.ps, i.var = c(1), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tsf.ps<-ggplot(plotdata.tsf,aes(time_since_fire,y))+geom_line()+ ylab("Probability")+ xlab("Time since fire (yr)")+ geom_text(aes(20,0.35, label="33%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.06,0.42))+ geom_rug(data=datps,aes(time_since_fire,PA50),sides="b") plotdata.dss <- plot.gbm(m.ps, i.var = c(2), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dss.ps<-ggplot(plotdata.dss,aes(dist_PSME_seed_source,y))+geom_line()+ ylab("Probability")+ xlab("Dist. seed source (m)")+ geom_text(aes(85,0.35, label="19%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9), plot.margin=margin(t = 5.5, r = 8, b = 5.5, l = 5.5, unit = "pt"))+ scale_y_continuous(limits=c(0.06,0.42))+ geom_rug(data=datps,aes(dist_PSME_seed_source,PA50),sides="b") plotdata.vwc <- plot.gbm(m.ps, i.var = c(5), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") vwc.ps<-ggplot(plotdata.vwc,aes(vwc_mean345,y))+geom_line()+ ylab("Probability")+ xlab("Spring soil moist. (vwc)")+ geom_text(aes(0.33,0.35, label="21%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.06,0.42))+ geom_rug(data=datps,aes(vwc_mean345,PA50),sides="b") plotdata.tskin<- plot.gbm(m.ps, i.var = c(4), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") tskin.ps<-ggplot(plotdata.tskin,aes(tskin_max,y))+geom_line()+ ylab("Probability")+ xlab(add_degree("Max. skin temp. "))+ geom_text(aes(58,0.35, label="18%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.06,0.42))+ geom_rug(data=datps,aes(tskin_max,PA50),sides="b") tskin.ps plotdata.dnbr<- plot.gbm(m.ps, i.var = c(3), n.trees = m.ps$gbm.call$best.trees, return.grid = TRUE,type="response") dnbr.ps<-ggplot(plotdata.dnbr,aes(dnbr,y))+geom_line()+ ylab("Probability")+ xlab("Fire severity (dNBR)")+ geom_text(aes(800,0.35, label="9%"),size=3.4,check_overlap = TRUE)+ theme_classic()+ theme(axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_y_continuous(limits=c(0.06,0.42))+ geom_rug(data=datps,aes(dnbr,PA50),sides="b") psme.pd<-plot_grid(tsf.ps,vwc.ps,dss.ps,dnbr.ps,tskin.ps,nrow=2,vjust = 1.2) psme.pd ## plot just partial dependency both species ### plot_grid(pipo.pd,psme.pd,nrow=2,labels="AUTO",vjust=1.1) ########################################################### ## predict for 1 yr postfire and constant dss, severity ####################### pred.brt<-function(sites,dat,nd.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PIPO_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean_driestM<-rep(dat.test$vwc_mean_driestM) nd$vpd_mean678<-rep(dat.test$vpd_mean678) nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.pi, nd, n.trees=m.pi$gbm.call$best.trees, type="response") nd$site<-i nd.sum<-rbind(nd.sum,nd) } assign("nd.fin.cur", nd.sum, envir = .GlobalEnv) } nd.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_PIPO_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.pi$site) pred.brt(sites,dat.c.pi,nd.sum) ### summarize by region ### nd.reg.cur<- nd.fin.cur %>% group_by(region,year) %>% summarize(pr=mean(pr)) nd.reg.cur<-data.frame(nd.reg.cur) head(nd.reg.cur) ## get sd of annual recruitment probability by region nd.reg.cur<-nd.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.reg.cur<-data.frame(nd.reg.cur) #psme pred.brt<-function(sites,dat,nd.ps.sum){ for(i in sites){ dat.test<-dat[dat$site==i,] nd<-expand.grid(time_since_fire=1,dist_PSME_seed_source=50, dnbr=412, year=dat.test$year) nd$vwc_mean345<-rep(dat.test$vwc_mean345) nd$tskin_max<-rep(dat.test$tskin_max) nd$region<-rep(dat.test$region) nd$pr<-predict.gbm(m.ps, nd, n.trees=m.ps$gbm.call$best.trees, type="response") nd$Site<-i nd.ps.sum<-rbind(nd.ps.sum,nd) } assign("nd.ps.fin.cur", nd.ps.sum, envir = .GlobalEnv) } nd.ps.sum<-setNames(data.frame(matrix(ncol = 9, nrow = 0)), c("time_since_fire", "dist_psme_seed_source", "dnbr", "year","vwc_mean_driestM","vpd_mean678","region","pr","site")) sites<-levels(dat.c.ps$site) pred.brt(sites,dat.c.ps,nd.ps.sum) ### summarize by region ### nd.ps.reg.cur<- ddply(nd.ps.fin.cur,.(region,year),summarize,pr=mean(pr)) ## get sd of annual recruitment probability by region nd.ps.reg.cur<-nd.ps.reg.cur %>% group_by(region) %>% mutate(sd.pr=rollapply(pr, width=10, sd, align='center', partial=F, fill=NA)) nd.ps.reg.cur<-data.frame(nd.ps.reg.cur) ### RP plots cols <- c("Both"="#0570b0","PIPO"="#74a9cf","PSME"="#bdc9e1","CA"="#CC79A7","CO"="#000000","NR"="#009E73","SW"="#E69F00")#new ones pipo.rp<-ggplot(nd.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability")+xlab("Year")+ theme_classic()+ scale_colour_manual(values = cols)+ theme(legend.position=c(0.97,1.04),legend.justification=c(1,1), legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5)) pipo.rp psme.rp<-ggplot(nd.ps.reg.cur,aes(year,pr,color=region))+geom_line()+ ylab("Recruitment probability")+xlab("Year")+ theme_classic()+ theme(#legend.position=c(0.97,1.15),legend.justification=c(1,1), legend.position="none", legend.key.height = unit(0.4,"cm"), legend.title=element_blank(), axis.text=element_text(color=1), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_colour_manual(name = "Region",values = cols)+ scale_x_continuous(breaks=seq(1980,2015,5)) psme.rp ## SD of RP pipo.sd<-ggplot(nd.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_smooth(aes(fill=region),se=F)+ scale_color_manual(values=cols)+ scale_fill_manual(name = "Region",values = cols)+ xlab("Year")+ ylab("SD of recruitment probability")+ theme(legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) pipo.sd psme.sd<-ggplot(nd.ps.reg.cur,aes(year,sd.pr,color=region))+ theme_classic()+ geom_smooth(aes(fill=region),se=F)+ scale_color_manual(values=cols)+ scale_fill_manual(name = "Region",values = cols)+ xlab("Year")+ ylab("SD of recruitment probability")+ theme(legend.position="none", legend.title=element_blank(), axis.text.x = element_text(size=8,color=1), axis.text.y = element_text(size=8,color=1), axis.title = element_text(size=9))+ scale_x_continuous(breaks=seq(1980,2015,5),labels=c("1980","","1990","","2000","","2010","")) psme.sd plot_grid(NULL,pipo.sd,pipo.rp,NULL,psme.sd,psme.rp,rel_widths=c(0.08,0.6,1),labels=c("A","","","B","",""))