############################################################################################################# # # Disentangling the effects of habitat suitability, dispersal and fragmentation on river fish distribution # # Radinger J. and Wolter C. # # R-script for plots: i) boxplots of habitat suitability (HS), disperal probability (DI) and barrier effects (BE) # contrasting absence and presence sites. ii) Graphical display of interaction effects of habitat, dispersal and # barriers on the presences of fish species. iii) Development of standardized effect sizes of habitat, dispersal # and barriers over nine modelled yearly time steps # ############################################################################################################# library(lme4) library(arm) library(reshape2) library(Hmisc) # Resultsdir results_dir <- "path/to/habitat_dispersal_barriers_raw_data.csv" # load dataframe and previousely generated GLMM model results habitat_dispersal_barriers_raw_data <- read.csv(paste(results_dir,"habitat_dispersal_barriers_raw_data.csv",sep="")) glmm_over_years <- readRDS(paste(results_dir,"glmm_over_years.rds",sep="")) #Created with previous script ####################### Prepare dataframe ########################## habitat_dispersal_barriers_raw_data$SiteID <- as.factor(habitat_dispersal_barriers_raw_data$SiteID) habitat_dispersal_barriers_raw_data$variable <- gsub("point_fidimo_corrected","dispersal",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$variable <- gsub("point_barrier_effect_corrected","barriers",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$variable <- gsub("point_MaxEnt","habitat",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$year <- unlist(lapply(strsplit(as.character(habitat_dispersal_barriers_raw_data$variable), "_"),"[",2)) habitat_dispersal_barriers_raw_data$variable <- unlist(lapply(strsplit(as.character(habitat_dispersal_barriers_raw_data$variable), "_"),"[",1)) habitat_dispersal_barriers_raw_data$variable <- factor(habitat_dispersal_barriers_raw_data$variable, levels=c("habitat","dispersal","barriers")) df_habitat_dispersal_barriers <- dcast(habitat_dispersal_barriers_raw_data, species+SiteID+presabs~variable,mean) df_habitat_dispersal_barriers$habitat_centered <- scale(df_habitat_dispersal_barriers$habitat,scale=FALSE,center=median(df_habitat_dispersal_barriers$habitat)) df_habitat_dispersal_barriers$dispersal_centered <- scale(df_habitat_dispersal_barriers$dispersal,scale=FALSE,center=median(df_habitat_dispersal_barriers$dispersal)) df_habitat_dispersal_barriers$barriers_centered <- scale(df_habitat_dispersal_barriers$barriers,scale=FALSE,center=median(df_habitat_dispersal_barriers$barriers)) ################################################################################################ ################################### Boxplots for HS, DI, BE #################################### cairo_pdf(paste(results_dir,"02_presabs_boxplots.pdf",sep=""),width=10,height=4) par(fig=c(0,.33,0,1),xpd=TRUE,mar=c(3.8,3.5,1,0),mgp=c(1.8,0.3,0)) boxplot((df_habitat_dispersal_barriers$habitat)~df_habitat_dispersal_barriers$presabs, #log="y", names=c("Absence","Presence"), ylab="Habitat suitability (HS)", axes=FALSE) axis(1,tcl=0.3,at=c(1,2), tick=FALSE, labels=c("Absence","Presence")) axis(2,tcl=0.3) yrange <- par("usr")[3:4] ypos <- yrange[1]-diff(yrange)/8 segments(1,ypos,2,ypos) segments(1,ypos,1,ypos+diff(yrange)/50) segments(2,ypos,2,ypos+diff(yrange)/50) text(1.5,ypos+diff(yrange)/50,"<") text(1.5,ypos-diff(yrange)/50,"***") text(0,yrange[2],"(a)") par(fig=c(0.33,0.66,0,1),new=TRUE,xpd=TRUE) boxplot(log(df_habitat_dispersal_barriers$dispersal+1)~df_habitat_dispersal_barriers$presabs, #log="y", names=c("Absence","Presence"), ylab="Species dispersal (DI)", yaxt="n", axes=FALSE) axis(1,tcl=0.3,at=c(1,2), tick=FALSE, labels=c("Absence","Presence")) axis(2,at=log(c(1,2,5,10,20,50,100)),labels=c(1,2,5,10,20,50,100),tcl=0.3) yrange <- par("usr")[3:4] ypos <- yrange[1]-diff(yrange)/8 segments(1,ypos,2,ypos) segments(1,ypos,1,ypos+diff(yrange)/50) segments(2,ypos,2,ypos+diff(yrange)/50) text(1.5,ypos+diff(yrange)/50,"<") text(1.5,ypos-diff(yrange)/50,"***") text(0,yrange[2],"(b)") par(fig=c(0.66,0.99,0,1),new=TRUE,xpd=TRUE) boxplot((df_habitat_dispersal_barriers$barriers)~df_habitat_dispersal_barriers$presabs, #log="y", names=c("Absence","Presence"), ylab="Barrier effects (BE)", ann=FALSE, axes=FALSE) axis(1,tcl=0.3,at=c(1,2), tick=FALSE, labels=c("Absence","Presence")) axis(2,tcl=0.3,at=seq(-3.5,0,0.5),labels=FALSE) axis(2,tcl=0.3,at=seq(-3,0,1),tick=FALSE) yrange <- par("usr")[3:4] ypos <- yrange[1]-diff(yrange)/8 segments(1,ypos,2,ypos) segments(1,ypos,1,ypos+diff(yrange)/50) segments(2,ypos,2,ypos+diff(yrange)/50) text(1.5,ypos+diff(yrange)/50,"<") text(1.5,ypos-diff(yrange)/50,".") text(0,yrange[2],"(c)") dev.off() plot.new() ################################################################################################ ################################### Simulations for displaying interactions #################### # Model mod_m1 <- glmer(presabs~barriers_centered+dispersal_centered+habitat_centered+barriers_centered*habitat_centered+dispersal_centered*habitat_centered+dispersal_centered*barriers_centered+(barriers_centered+habitat_centered+dispersal_centered|species), family=binomial,data=df_habitat_dispersal_barriers) # Simulated model coefficients (Approach adapted from Burton et. al 2013, Early maternal experience shapes offspring performance in the wild. Ecology, 94, 618-626) n_sim <- 10000 sim_mod_m1 <- sim(mod_m1, n_sim)@fixef colnames(sim_mod_m1) <- names(coef(mod_m1)$species) ########### Interaction HS*BE ###################### length_HS_seq <- 100 HS_seq <- seq(from=min(df_habitat_dispersal_barriers$habitat), to=1, #to=max(df_habitat_dispersal_barriers$habitat), length=length_HS_seq) HS_matrix <- matrix(HS_seq, nr = n_sim, nc = length_HS_seq, byrow = TRUE) length_DI_seq <- 100 DI_seq <- seq(from=min(df_habitat_dispersal_barriers$dispersal), to=0.6, #to=quantile(df_habitat_dispersal_barriers$dispersal,probs=0.8), length=length_DI_seq) DI_matrix <- matrix(DI_seq, nr = n_sim, nc = length_DI_seq, byrow = TRUE) length_BE_seq <- 100 BE_seq <- seq(from=min(df_habitat_dispersal_barriers$barriers), to=quantile(df_habitat_dispersal_barriers$barriers,probs=1), length=length_BE_seq) BE_matrix <- matrix(BE_seq, nr = n_sim, nc = length_BE_seq, byrow = TRUE) #fixed variables BE_low <- quantile(df_habitat_dispersal_barriers$barriers,probs=c(0.99)) BE_high <- quantile(df_habitat_dispersal_barriers$barriers,probs=c(0.01)) HS_high <- 0.6 HS_low <- 0.2 DI_high <- quantile(df_habitat_dispersal_barriers$dispersal,probs=c(0.99)) DI_low <- quantile(df_habitat_dispersal_barriers$dispersal,probs=c(0.01)) ### Calculated model results based on simulated betas ### ########### Interaction HS*BE ###################### # for DI fixed to median, varying HS and 2-fixed barriers (min, max) presabs_modelled_BE_high_HS <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (BE_high -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (BE_high -median(df_habitat_dispersal_barriers$barriers)) * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal)) * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (BE_high-median(df_habitat_dispersal_barriers$barriers)) * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal))) presabs_modelled_BE_low_HS <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (BE_low -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (BE_low -median(df_habitat_dispersal_barriers$barriers)) * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal)) * (HS_matrix-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (BE_low-median(df_habitat_dispersal_barriers$barriers)) * (median(df_habitat_dispersal_barriers$dispersal) - median(df_habitat_dispersal_barriers$dispersal))) ########### Interaction BE*DI ###################### # for HS fixed to median, varying DI and 2-fixed barriers (min, max) presabs_modelled_BE_high_DI <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (BE_high -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (BE_high -median(df_habitat_dispersal_barriers$barriers)) * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (BE_high-median(df_habitat_dispersal_barriers$barriers)) * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal))) presabs_modelled_BE_low_DI <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (BE_low -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (BE_low -median(df_habitat_dispersal_barriers$barriers)) * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) * (median(df_habitat_dispersal_barriers$habitat)-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (BE_low-median(df_habitat_dispersal_barriers$barriers)) * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal))) ########### Interaction HS*DI ###################### # for BE fixed to median, varying HS and 2-fixed barriers (min, max) presabs_modelled_HS_high_DI <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (0 -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (HS_high-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (0 -median(df_habitat_dispersal_barriers$barriers)) * (HS_high-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) * (HS_high-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (0-median(df_habitat_dispersal_barriers$barriers)) * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal))) presabs_modelled_HS_low_DI <- plogis(sim_mod_m1[,"(Intercept)"] + sim_mod_m1[,"barriers_centered"] * (0 -median(df_habitat_dispersal_barriers$barriers)) + sim_mod_m1[,"dispersal_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) + sim_mod_m1[,"habitat_centered"] * (HS_low-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:habitat_centered"] * (0 -median(df_habitat_dispersal_barriers$barriers)) * (HS_low-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"dispersal_centered:habitat_centered"] * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal)) * (HS_low-median(df_habitat_dispersal_barriers$habitat)) + sim_mod_m1[,"barriers_centered:dispersal_centered"] * (0-median(df_habitat_dispersal_barriers$barriers)) * (DI_matrix - median(df_habitat_dispersal_barriers$dispersal))) ################## ## Plots of interactions ### cairo_pdf(paste(results_dir,"04_GLMM_interactions.pdf",sep=""),width=10,height=4) # create x and y axes par(fig=c(0,.33,0,1),xpd=TRUE,mar=c(3.8,3.5,1,0),mgp=c(1.8,0.3,0)) plot(0, 0, xlim = c(0,1), ylim = c(0, 1), type= "n", xlab= "HS", ylab= "Probability of presence", axes=FALSE) axis(1,tcl=0.3) axis(2,tcl=0.3) yrange <- par("usr")[3:4] xrange <- par("usr")[1:2] text(xrange[1]-diff(xrange)*0.185,1,"(a)") # Draw polygons of 95% CI polygon(c(HS_seq,rev(HS_seq)), c(apply(presabs_modelled_BE_high_HS, 2, quantile, 0.025), rev(apply(presabs_modelled_BE_high_HS, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) polygon(c(HS_seq,rev(HS_seq)), c(apply(presabs_modelled_BE_low_HS, 2, quantile, 0.025), rev(apply(presabs_modelled_BE_low_HS, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) # plot fitted values & 95% CI for barriers max (0 barrier effect) lines(HS_seq , apply(presabs_modelled_BE_high_HS, 2, mean), lty=1, lwd=2, col="black") lines(HS_seq , apply(presabs_modelled_BE_high_HS, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(HS_seq , apply(presabs_modelled_BE_high_HS, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # plot fitted values & 95% CI for barriers min (maximum barrier effect) lines(HS_seq , apply(presabs_modelled_BE_low_HS, 2, mean), lty=4, lwd=2, col="black") lines(HS_seq , apply(presabs_modelled_BE_low_HS, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(HS_seq , apply(presabs_modelled_BE_low_HS, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # legend legend(0.5*xrange[2],0.2,c("high BE","low BE"), lty=c(1,4), lwd=2, bty="n", seg.len=1.5) # create x and y axes par(fig=c(0.33,0.66,0,1),new=TRUE,xpd=TRUE) plot(0, 0, xlim = c(0,0.6), ylim = c(0, 1), type= "n", xlab= "DI", ylab= "Probability of presence", axes=FALSE) axis(1,tcl=0.3) axis(2,tcl=0.3) yrange <- par("usr")[3:4] xrange <- par("usr")[1:2] text(xrange[1]-diff(xrange)*0.185,1,"(b)") # Draw polygons of 95% CI polygon(c(DI_seq,rev(DI_seq)), c(apply(presabs_modelled_HS_high_DI, 2, quantile, 0.025), rev(apply(presabs_modelled_HS_high_DI, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) polygon(c(DI_seq,rev(DI_seq)), c(apply(presabs_modelled_HS_low_DI, 2, quantile, 0.025), rev(apply(presabs_modelled_HS_low_DI, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) # plot fitted values & 95% CI for high HS (e.g. 0.6) lines(DI_seq , apply(presabs_modelled_HS_high_DI, 2, mean), lty=1, lwd=2, col="black") lines(DI_seq , apply(presabs_modelled_HS_high_DI, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(DI_seq , apply(presabs_modelled_HS_high_DI, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # plot fitted values & 95% CI for low HS (eg. 0.1) lines(DI_seq , apply(presabs_modelled_HS_low_DI, 2, mean), lty=4, lwd=2, col="black") lines(DI_seq , apply(presabs_modelled_HS_low_DI, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(DI_seq , apply(presabs_modelled_HS_low_DI, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # legend legend(0.5*xrange[2],0.2,c("high HS","low HS"), lty=c(1,4), lwd=2, bty="n", seg.len=1.5) # create x and y axes par(fig=c(0.66,0.99,0,1),new=TRUE,xpd=TRUE) plot(0, 0, xlim = c(0,0.6), ylim = c(0, 1), type= "n", xlab= "DI", ylab= "Probability of presence", axes=FALSE) axis(1,tcl=0.3) axis(2,tcl=0.3) yrange <- par("usr")[3:4] xrange <- par("usr")[1:2] text(xrange[1]-diff(xrange)*0.185,1,"(c)") # Draw polygons of 95% CI polygon(c(DI_seq,rev(DI_seq)), c(apply(presabs_modelled_BE_high_DI, 2, quantile, 0.025), rev(apply(presabs_modelled_BE_high_DI, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) polygon(c(DI_seq,rev(DI_seq)), c(apply(presabs_modelled_BE_low_DI, 2, quantile, 0.025), rev(apply(presabs_modelled_BE_low_DI, 2, quantile, 0.975))), border=NA, col=adjustcolor("grey70",alpha=0.5)) # plot fitted values & 95% CI for high HS (e.g. 0.6) lines(DI_seq , apply(presabs_modelled_BE_high_DI, 2, mean), lty=1, lwd=2, col="black") lines(DI_seq , apply(presabs_modelled_BE_high_DI, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(DI_seq , apply(presabs_modelled_BE_high_DI, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # plot fitted values & 95% CI for low HS (eg. 0.1) lines(DI_seq , apply(presabs_modelled_BE_low_DI, 2, mean), lty=4, lwd=2, col="black") lines(DI_seq , apply(presabs_modelled_BE_low_DI, 2, quantile, 0.025), lty=1, lwd=0.25, col="gray28") lines(DI_seq , apply(presabs_modelled_BE_low_DI, 2, quantile, 0.975), lty=1, lwd=0.25, col="gray28") # legend legend(0.5*xrange[2],0.2,c("high BE","low BE"), lty=c(1,4), lwd=2, bty="n", seg.len=1.5) dev.off() plot.new() ################################################################################################ ######## Plot development of standardized parameter coefs (Betas) over years ################## plot_estimates_df <- do.call("rbind",glmm_over_years) plot_estimates_df$varname <- unlist(lapply(strsplit(rownames(plot_estimates_df), "\\."),"[",2)) row.names(plot_estimates_df) <- NULL plot_estimates_df <-droplevels(plot_estimates_df[!plot_estimates_df$varname %in% c("(Intercept)", "sd_(Intercept)|species"),]) plot_estimates_df$varname <- factor(plot_estimates_df$varname, levels=c("DI","HS","BE")) ### Start Plot cairo_pdf(paste(results_dir,"05_Standardized_estimates_over_years.pdf",sep=""),width=5,height=5) #### Start with DI par(new=F,mar=c(3.8,3.5,1,0),mgp=c(1.5,0.3,0)) plot(x=plot_estimates_df$Time[plot_estimates_df$varname=="DI"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="DI"], type="l", lty="dotted", col="grey40", axes=FALSE, ylim=c(0,15), xlab="", ylab="") points(x=plot_estimates_df$Time[plot_estimates_df$varname=="DI"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="DI"], pch=16, col="white", cex=3.5) par(new=T) errbar(x=plot_estimates_df$Time[plot_estimates_df$varname=="DI"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="DI"], yplus=plot_estimates_df$BootCI_upr[plot_estimates_df$varname=="DI"], yminus=plot_estimates_df$BootCI_lwr[plot_estimates_df$varname=="DI"], axes=FALSE, ylim=c(0,15), xlab="", ylab="", cex=0) points(x=plot_estimates_df$Time[plot_estimates_df$varname=="DI"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="DI"], pch=24, bg="white") ##### HS lines(x=plot_estimates_df$Time[plot_estimates_df$varname=="HS"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="HS"], lty="dotted", col="grey40") points(x=plot_estimates_df$Time[plot_estimates_df$varname=="HS"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="HS"], pch=16, col="white", cex=3.5) par(new=T) errbar(x=plot_estimates_df$Time[plot_estimates_df$varname=="HS"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="HS"], yplus=plot_estimates_df$BootCI_upr[plot_estimates_df$varname=="HS"], yminus=plot_estimates_df$BootCI_lwr[plot_estimates_df$varname=="HS"], axes=FALSE, ylim=c(0,15), xlab="", ylab="", cex=0) points(x=plot_estimates_df$Time[plot_estimates_df$varname=="HS"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="HS"], pch=21, bg="white") #### BE lines(x=plot_estimates_df$Time[plot_estimates_df$varname=="BE"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="BE"], lty="dotted", col="grey40") points(x=plot_estimates_df$Time[plot_estimates_df$varname=="BE"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="BE"], pch=16, col="white", cex=3.5) par(new=T) errbar(x=plot_estimates_df$Time[plot_estimates_df$varname=="BE"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="BE"], yplus=plot_estimates_df$BootCI_upr[plot_estimates_df$varname=="BE"], yminus=plot_estimates_df$BootCI_lwr[plot_estimates_df$varname=="BE"], axes=FALSE, ylim=c(0,15), xlab="", ylab="", cex=0) points(x=plot_estimates_df$Time[plot_estimates_df$varname=="BE"], y=plot_estimates_df$Estimates[plot_estimates_df$varname=="BE"], pch=22, bg="white") ### Axis axis(1,tcl=0.3,at=c(1:9), labels=c(1:9)) axis(2,tcl=0.3,at=c(0:14),labels=FALSE) axis(2,at=seq(0,14,2),labels=seq(0,14,2),tick=FALSE) title(xlab="Time (years)",ylab=expression(paste("Standardized Parameter Estimates ",beta[i])),outer=FALSE) ### Legend points(x=5.5,y=14,pch=24,bg="white") points(x=5.5,y=13,pch=21,bg="white") points(x=5.5,y=12,pch=22,bg="white") text(x=5.7,y=14,labels=expression(paste("Species dispersal ",beta[DI])),adj=0) text(x=5.7,y=13,labels=expression(paste("Habitat suitability ",beta[HS])),adj=0) text(x=5.7,y=12,labels=expression(paste("Barrier effects ",beta[BE])),adj=0) dev.off()