############################################################# ## ## Hargreaves & Eckert ECOLOGY LETTERS ## ## Matrix models to estimate population viability, taking dormant seeds into account ## Results of this script are in Fig S5 ## code for Fig 5 at end of script ############################################################# ## MATRIX MODEL BACKGROUND ##### # 3 life history stages: s1 (1styear seeds), s2 (dormant seeds), P (plants) # s1 can make s2 or P # s2 can make s2 or P # P can make s1 tmatrix=matrix(c(1,2,3,4,5,6,7,8,9),ncol=3,nrow=3) #matrix: tmatrix #1,1 - p s1 -> s1 = 0 #2,1 - p s2 -> s1 = 0 #3,1 - p P -> s1 = mean LTRS #1,2 - p s1 -> s2 = dormancy fraction #2,2 - p s2 -> s2 = dormancy fraction #3,2 - p P -> s2 = 0 #1,3 - p s1 -> P = prop Em #2,3 - p s2 -> P = prop Em 2nd year #3,3 - p P -> P = 0 #basic model w example (ie made up) matrix state=c(0,0,100) #starting state = 0 seeds, 0 dormant seeds, 100 adults tmatrix=matrix(c(0,0,10,.5,0,0,0,.5,.8),ncol=3,nrow=3) #matrix: tmatrix for(i in 1:100){ #for each of the 100 plants mem_state=state state=state %*% tmatrix #treating all individuals/seeds as individuals (i.e., steadstate) - but could also base it on adult growth growth=sum(state)/sum(mem_state) print(growth) } ## DATA ################ NKp <- read.csv(file="~/Desktop/Hargreaves&EckertEcolLett data/HargreavesEcolLett data prop<FNK.csv", header=T); NKp$X <- NULL; summary(NKp) NKp$site <- factor(NKp$site, levels=c('L','M','HE','A1','A2')) NKp$Source <- factor(NKp$Source, levels=c('L','M','H','E','A','MT','F')) NKp$year <- as.factor(NKp$year) #make sure year is factor summary(NKp) #for matrix models need following: 1) prop 1st year seeds that emerge (s1s2), 2) seeds per emerged plant (Ps1), 3) prop seeds that stay dormant (use prop that emerge in 2nd year) #1) prop 1st year seeds that emerge (s1s2) - use pEm across all plots (ie include TR and OTC plots) NKp$pEm <- NKp$nEm / NKp$nseed #use nEm (not nEmzt): all high-elevation sites checked at same phenological stage every year NKEm <- aggregate(NKp[c('pEm')], by=NKp[c('site','Source')], mean, na.rm=T); head(NKEm) #2) seeds per emerged plant (Ps1) - distinguish between TR (control) and OTC (warmed) plots NKp$seedperEm <- NKp$seedperplot / NKp$nEm NKav <- aggregate(NKp[c('seedperEm')], by=NKp[c('site','treat','Source')], mean, na.rm=T); NKav #3) prop seeds that stay dormant (use prop that emerge in 2nd year) NKd <- read.csv(file="~/Desktop/Hargreaves&EckertEcolLett data/HargreavesEcolLett data dormancyNK.csv", header=T); NKd$X <- NULL #EXPLANATION OF VARIABLES (variables used from NKp are the same) summary(NKd) #transect=elevational transect (NK=Nakiska, HB=Hailstone Butte) #site=site seeds were transplanted into. L=low elevation, M=mid-elevation (the range centre), H=high-elevation (within 100 masl of absolute high-elevation range edge), HE=high-elevation range edge, A=above the range #yrplanted=year seeds were planted #plot=plot, marked with central wooden stake, with multiple subplots (1 subplot per source population) #Source=source population: low-elevation from that transet (L), mid-elevation from that transect (M), high-elevation (within 100 masl of absolute high-elevation range edge) from that transect (H), edge (absolute high-elevation range edge E), high-elevation (2000m) populations on neighbouring Moose Mountain (MT) and Fortress Mountain (F), and seeds produced by plants transplanted above the range the previous year (A) #year=year plants were active (ie year seeds germinated) #yrEm=year of emergence. Non-dormant seeds germinate in yrEm=1. Plants that emerge the second summer after they were planted yrEm=2 etc. #treat=dormant #nseed.orig=# of seeds planted in the subplot (usually 25) #nseed.rem=# of seeds remaining in the subplot after previous years of potential emergence. When yrEm=1, nseed.rem=nseed.orig #nEm=# of plants that emerged in the subplot that year (cannot be > nseed.rem) #nFlwr=# of plants in the subplot that flowered that year (cannot be > nEm) #nR=# of plants in the subplot that produced viable seed (cannot be > nFlwr) #avLn1new=average # of primary stem leaf nodes (not counting the cotelydon) plants in that subplot had the first time they were seen. Lowest possible value = 0 if all plants were first seen when they had only cotyledons #avJnew=average Julien day (Jan 1st = 1) plants were first seen #excl.str=stricter criteria for excluding subplots that might have been contaminated by seed from previous years. =excl if even 1 fruit opened in the subplot & was old enough to have dropped seed. NKd[is.na(NKd$nseed.rem),] NKd$pEm.fromrem <- NKd$nEm / NKd$nseed.rem #all from low & mid sites in yrs 3+ - result of missing a dormancy check in an early year summary(NKd) NKyEm2 <- NKd[NKd$yrEm==2 & (NKd$site=='HE' | NKd$site=='A1' | NKd$site=='A2') & NKd$Source!='A',]; NKyEm2 <- droplevels(NKyEm2); summary(NKyEm2) NKd <- aggregate(NKyEm2[c('pEm.fromrem')], by=NKyEm2[c('site','Source')], mean, na.rm=T); NKd colnames(NKd) <- c('site','Source','pEmyr2') #get overall dormancy fraction for highest three sites overalld <- aggregate(NKyEm2[c('pEm.fromrem')], by=NKyEm2[c('transect')], mean, na.rm=T); overalld #0.0724 # HB prop data ### HBp <- read.csv(file="~/Desktop/Hargreaves&EckertEcolLett data/HargreavesEcolLett data prop<FHB.csv", header=T); HBp$X <- NULL; summary(HBp) HBp$site <- factor(HBp$site, levels=c('L','M','H','HE','A')) HBp$Source <- factor(HBp$Source, levels=c('L','M','H','E','A','MT','F')) HBp$year <- as.factor(HBp$year) #make sure year is factor summary(HBp) #for matrix models need following: 1) prop 1st year seeds that emerge (s1s2), 2) seeds per emerged plant (Ps1), 3) prop seeds that stay dormant (use prop that emerge in 2nd year) #1) prop 1st year seeds that emerge (s1s2) - use pEm across all plots (ie include TR and OTC plots) HBp$pEm <- HBp$nEm / HBp$nseed #use nEm (not nEmzt): all high-elevation sites checked at same phenological stage every year HBEm <- aggregate(HBp[c('pEm')], by=HBp[c('site','Source')], mean, na.rm=T); head(HBEm) #2) seeds per emerged plant (Ps1) - distinguish between TR (control) and OTC (warmed) plots HBp$seedperEm <- HBp$seedperplot / HBp$nEm HBav <- aggregate(HBp[c('seedperEm')], by=HBp[c('site','treat','Source')], mean, na.rm=T); HBav #3) prop seeds that stay dormant (use prop that emerge in 2nd year) HBd <- read.csv(file="~/Desktop/Hargreaves&EckertEcolLett data/HargreavesEcolLett data dormancyHB.csv", header=T); HBd$X <- NULL summary(HBd) HBd[is.na(HBd$nseed.rem),] HBd$pEm.fromrem <- HBd$nEm / HBd$nseed.rem #all from mid sites in yrs 3 - result of missing a dormancy check in an early year summary(HBd) #get dormant fraction in second year for highest sites only HByEm2 <- HBd[HBd$yrEm==2 & (HBd$site=='H' | HBd$site=='HE' | HBd$site=='A') & HBd$Source!='A',]; head(HByEm2) HBd <- aggregate(HByEm2[c('pEm.fromrem')], by=HByEm2[c('site','Source')], mean, na.rm=T); HBd colnames(HBd) <- c('site','Source','pEmyr2') summary(HBd) #get overall dormancy fraction for top three sites overalld <- aggregate(HByEm2[c('pEm.fromrem')], by=HByEm2[c('transect')], mean, na.rm=T); overalld #0.032 ## NK matrix models ########## #create dummy data frame for values to go into df <- data.frame(transect='NK', site='test', Source='test', treat='test', s1s2=1, seedperEm_Ps1=1, pEm_s1P=1, lambda=1); df head(NKav) trn='NK' #PROCEDURE: execute a line of set up code (eg s='HE'; S='L'; t='TR'), then run model except the write.csv line, repeat for next line of set up code (eg s='A1'; S='L' ; t='TR') etc. Once have run model for all lines of set up code can write the csv) #low seeds s='HE'; S='L'; t='TR' s='A1'; S='L' ; t='TR' s='A2'; S='L' ; t='TR' #mid seeds s='HE'; S='M'; t='TR' #fluctuates btwn 2 values s='A1'; S='M'; t='TR' s='A2'; S='M'; t='TR' #high seeds ### s='HE'; S='H'; t='TR' #still fluctuating s='A1'; S='H'; t='TR' s='A2'; S='H'; t='TR' #edge seeds ### s='HE'; S='E'; t='TR' s='A2'; S='E'; t='TR' #>1! but only had E seeds a A2 in best year so take with grain of salt #MT seeds ### s='HE'; S='MT'; t='TR' s='A1'; S='MT'; t='TR' s='A2'; S='MT'; t='TR' #at NKA2 >1! #F seeds ### s='HE'; S='F'; t='TR' s='A1'; S='F'; t='TR' s='A2'; S='F'; t='TR' #OTC NK s='HE'; S='M'; t='OTC' s='HE'; S='H'; t='OTC' s='HE'; S='MT'; t='OTC' s='A2'; S='M'; t='OTC' s='A2'; S='H'; t='OTC' s='A2'; S='MT'; t='OTC' #matrix runs NK ### state=c(0,0,100) #starting state = 0 seeds, 0 dormant seeds, 100 adults #define matrix s1s1 = s2s1 = Ps2 = pP = 0 #impossible transitions Ps1 = NKav$seedperEm[NKav$site==s & NKav$Source==S & NKav$treat==t]; Ps1 #s1s2 = s2s2 = NKd$pEmyr2[NKd$site==s & NKd$Source==S]; s1s2 #option to use dormancy fraction for that source at that site s1s2 = s2s2 = .072 #using overall dormancy fraction from highest 3 sites (USED in paper) s1P = NKEm$pEm[NKEm$site==s & NKEm$Source==S]; s1P #use emergence fraction across treatments s2P = s1P tmatrix=matrix(c(s1s1,s2s1,Ps1,s1s2,s2s2,Ps2,s1P,s2P,pP),ncol=3,nrow=3); tmatrix #matrix: #run for(i in 1:500){ #repeat 100 times mem_state=state #record of starting state state=state %*% tmatrix #multiply matrix by vector state, returns new vector with new number of 1st year seeds, dormant seeds, plants #treating all individuals/seeds as individuals (i.e., steadstate) - but could also base it on adult growth growth=sum(state)/sum(mem_state) # sum(state) = pop size in time step 2 (adults + seeds), sum(mem_state) = pop size in time step 1 (adults + seeds). print(growth) } #make vector new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda=growth[1]); new #new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda= 0.2072922); new #if pop goes to zero before end of run #add vector to data frame df <- rbind(df, new) ; df #df$pEmfrom<-NULL df$pEmfrom <- 'TR&OTCplots at same site as P1s' #use if used overall odormancy fraction #df$dormancy <- 'by siteXsource'; #if used different dormancy fraction for each site x source combination #write.csv(df, "~/Desktop/???") ## HB matrix models ########## #HB create dummy data frame for values to go into df <- data.frame(transect='HB', site='test', Source='test', treat='test', s1s2=1, seedperEm_Ps1=1, pEm_s1P=1, lambda=1); df head(HBav) trn='HB' #low seeds s='H'; S='L' ; t='TR' s='HE'; S='L' ; t='TR' s='A'; S='L' ; t='TR' #mid seeds s='H'; S='M'; t='TR' s='HE'; S='M'; t='TR' s='A'; S='M'; t='TR' #high seeds ## s='H'; S='H'; t='TR' s='HE'; S='H'; t='TR' s='A'; S='H'; t='TR' #edge seeds s='HE'; S='E'; t='TR' #MT seeds s='H'; S='MT'; t='TR' s='HE'; S='MT'; t='TR' s='A'; S='MT'; t='TR' #F seeds s='H'; S='F'; t='TR' s='HE'; S='F'; t='TR' s='A'; S='F'; t='TR' #OTC HB s='H'; S='M'; t='OTC' #get oscillating pop size bc plants have very low chance of reproducing but then make tones of seed. but no oscillating if use local dormancy s='H'; S='H'; t='OTC' s='H'; S='MT'; t='OTC' s='A'; S='M'; t='OTC' s='A'; S='H'; t='OTC' s='A'; S='MT'; t='OTC' #matrix runs HB ### state=c(0,0,100) #starting state = 0 seeds, 0 dormant seeds, 100 adults #define matrix s1s1 = s2s1 = Ps2 = pP = 0 #impossible transitions Ps1 = HBav$seedperEm[HBav$site==s & HBav$Source==S & HBav$treat==t]; Ps1 s1s2 = s2s2 = 0.032 #use single dormancy fraction from top three sites only (new recalculated value) #s1s2 = s2s2 = HBd$pEmyr2[HBd$site==s & HBd$Source==S]; s1s2 #let d vary by sitex source s1P = HBEm$pEm[HBEm$site==s & HBEm$Source==S]; s1P #use emergence fraction across treatments s2P = s1P tmatrix=matrix(c(s1s1,s2s1,Ps1,s1s2,s2s2,Ps2,s1P,s2P,pP),ncol=3,nrow=3); tmatrix #matrix: #run for(i in 1:1000){ #repeat 100 times mem_state=state #record of starting state state=state %*% tmatrix #multiply matrix by vector state, returns new vector with new number of 1st year seeds, dormant seeds, plants #treating all individuals/seeds as individuals (i.e., steadstate) - but could also base it on adult growth growth=sum(state)/sum(mem_state) # sum(state) = pop size in time step 2 (adults + seeds), sum(mem_state) = pop size in time step 1 (adults + seeds). print(growth) } #make vector new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda=growth[1]); new new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda= 0.84255); new #if pop goes to zero before end of run #if dormancy = zero then lambda really does oscillate #add vector to data frame df <- rbind(df, new) ; df df$pEmfrom <- 'TR&OTCplots at same site as P1s' # df$dormancy <- 'by siteXsource' df$dormancy <- 'av across HEAandsource' #write.csv(df, "~/Desktop/???.csv") #ones to check - L source at HE with sitexsource dormancy (d=0). get stable oscillations. take mean (3.964665+ 3.968773)/2 ## redo OTC treatment with Emergence values from range edge - ie would warming the growing season make site sustainable ##### #create dummy data frame for values to go into df <- data.frame(transect='test', site='test', Source='test', treat='test', s1s2=1, seedperEm_Ps1=1, pEm_s1P=1, lambda=1, pEmfrom='test'); df #OTC HB A using Em from HB-H trn='HB'; s='A'; S='M'; t='OTC'; sEm='H' #doesnt quite stabilize s='A'; S='H'; t='OTC'; sEm='H' s='A'; S='MT'; t='OTC'; sEm='H' state=c(0,0,100) #starting state = 0 seeds, 0 dormant seeds, 100 adults #define matrix s1s1 = s2s1 = Ps2 = pP = 0 #impossible transitions Ps1 = HBav$seedperEm[HBav$site==s & HBav$Source==S & HBav$treat==t]; Ps1 s1s2 = s2s2 = 0.032 #dormancy fraction, unchanging s1P = HBEm$pEm[HBEm$site==sEm & HBEm$Source==S]; s1P #use emergence fraction from H site #s1P = HBEm$pEm[HBEm$site==s & HBEm$Source==S]; s1P #use emergence fraction from H site s2P = s1P tmatrix=matrix(c(s1s1,s2s1,Ps1,s1s2,s2s2,Ps2,s1P,s2P,pP),ncol=3,nrow=3); tmatrix #matrix: #run for(i in 1:1000){ #repeat 100 times mem_state=state #record of starting state state=state %*% tmatrix #multiply matrix by vector state, returns new vector with new number of 1st year seeds, dormant seeds, plants #treating all individuals/seeds as individuals (i.e., steadstate) - but could also base it on adult growth growth=sum(state)/sum(mem_state) # sum(state) = pop size in time step 2 (adults + seeds), sum(mem_state) = pop size in time step 1 (adults + seeds). print(growth) } new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda=growth[1], pEmfrom=sEm); new new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda=1.932889, pEmfrom=sEm); new #add vector to data frame df <- rbind(df, new) ; df #OTC NK A2 trn='NK'; s='A2'; S='M'; t='OTC'; sEm='HE' s='A2'; S='H'; t='OTC'; sEm='HE' s='A2'; S='MT'; t='OTC'; sEm='HE' state=c(0,0,100) #starting state = 0 seeds, 0 dormant seeds, 100 adults #define matrix s1s1 = s2s1 = Ps2 = pP = 0 #impossible transitions Ps1 = NKav$seedperEm[NKav$site==s & NKav$Source==S & NKav$treat==t]; Ps1 s1s2 = s2s2 = .072 #dormancy fraction, unchanging from top 3 sites s1P = NKEm$pEm[NKEm$site==sEm & NKEm$Source==S]; s1P #use emergence fraction from H site #s1P = HBEm$pEm[HBEm$site==s & HBEm$Source==S]; s1P #use emergence fraction from H site s2P = s1P tmatrix=matrix(c(s1s1,s2s1,Ps1,s1s2,s2s2,Ps2,s1P,s2P,pP),ncol=3,nrow=3); tmatrix #matrix: #run for(i in 1:500){ #repeat 100 times mem_state=state #record of starting state state=state %*% tmatrix #multiply matrix by vector state, returns new vector with new number of 1st year seeds, dormant seeds, plants #treating all individuals/seeds as individuals (i.e., steadstate) - but could also base it on adult growth growth=sum(state)/sum(mem_state) # sum(state) = pop size in time step 2 (adults + seeds), sum(mem_state) = pop size in time step 1 (adults + seeds). print(growth) } #make vector new <- data.frame(transect=trn, site=s, Source=S, treat=t, s1s2=s1s2, seedperEm_Ps1=Ps1, pEm_s1P=s1P, lambda=growth[1], pEmfrom=sEm); new #new <- data.frame(site=s, Source=S, treat=t, s1s2=s1s2, lambda= 0.1474157); new #if pop goes to zero before end of run do by hand #add vector to data frame df <- rbind(df, new) ; df #write.csv(df, "~/Desktop/???.csv") ## FIGURE S5 ##### nrr <- read.csv(file="~/Desktop/Hargreaves&EckertEcolLett data/HargreavesEcolLett lambda from matrix models.csv", header=T) nrr$site <- factor(nrr$site, levels=c('H','HE','A1','A','A2')) nrr$Source <- factor(nrr$Source, levels=c('L','M','H','E','MT','F')) summary(nrr) table(nrr$site, nrr$dormancyfrom) bwplot(lambda ~ Source|site, data=nrr, subset=transect=='NK') #make shorter code for dormancy nrr$dorm <- as.factor(ifelse(nrr$dormancy=='av across HEAandsource', 'av', 'sXS')) nrr <- nrr[with(nrr, order(transect, pEmfrom, dorm, Source, site)), ] summary(nrr) quartz(width=4.5,height=4.8, family="Helvetica") outp=1.4; inp=1.2; #cex of outer (black) point & inner coloured point cexax=.8; cexlb=1; cexst=.75; #font size of axis labels, transect label, stats #general par(mfrow=c(2,2), oma=c(2.5,2.5,1,0), bty='l', pch=20); marNK=c(1.5,.7,.5,0); marHB=c(1.5,0,.5,.7) ltyh=1 #line type for horizontal ref line at 1 (ie sustainability if all seeds germinated) #axes & legend xmin=.75; xmaxNK=3.6; xmaxHB=3.5; abxNK=1.3; abxHB=2.3; lymin=.05; lymax=4.3; yl=3 #log scale ab=1.2 #line width for vertical range limit abline leg=15 #ymax for legend in ltf panel - set to0 high if dont need #lines NK1=1; HB1=1; n=3; nHB=3 # sets lengths of lines. since analysis is only top three sites just have lines for those l =1.1 #linewidth for lines connecting points for each source, set to 0 if dont want lines #points colL='firebrick3'; colM='olivedrab3'; colH='deepskyblue2'; colE='cyan'; colMT=colH; colF=colH; colOTC='orange' #colour by elevation pchL=19; pchM=19; pchH=19; pchE=19; pchMT=17; pchF=15 #if want shape to vary by elevation j=.05; jL=-2; jM=-1; jH=0; jE=1; jMT=2; jF=3 #amount by which to jitter points - even spacing #NK Lifetime fitness - lambda from matrix model ### x=c(1,2,3); ylim=c(lymin,lymax); xlimNK=c(xmin, xmaxNK); par(mar=marNK) xL=x+jL*j; yL=nrr$lambda[nrr$Source=='L' & nrr$transect=='NK' & nrr$dorm=='av'] plot(xL, yL, pch=pchL, cex=outp, ylim=ylim, xlim=xlimNK, xaxt="n", ylab='', xlab="", las=1) abline(h=1, col='grey45', lwd=ab, lty=ltyh) # line showing lambda=1 = get in early so points are on top abline(v=abxNK, col='grey45', lwd=ab, lty = 2) # line showing RL xM=x+jM*j; yM=nrr$lambda[nrr$Source=='M' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='TR'] xH=x+jH*j; yH=nrr$lambda[nrr$Source=='H' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='TR'] xE=c(1,3)+jE*j; yE=nrr$lambda[nrr$Source=='E' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='TR'] xMT=x+jMT*j; yMT=nrr$lambda[nrr$Source=='MT' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='TR'] xF=x+jF*j; yF=nrr$lambda[nrr$Source=='F' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='TR'] lines(xL[NK1:n], yL[NK1:n], col=colL, lwd=l); lines(xM[NK1:n], yM[NK1:n], col=colM, lwd=l); lines(x=xH[NK1:n], yH[NK1:n], col=colH, lwd=l) #lines for L,M,H lines(xMT, yMT, col=colMT, lwd=l); lines(xF, yF, col=colF, lwd=l) #lines for extra sources points(xL, yL, col=colL, pch=pchL, cex=inp) points(x=xM, yM, pch=pchM, cex=outp); points(x=xM, yM, col=colM, pch=pchM, cex=inp) points(x=xH, y=yH, pch=pchH, cex=outp); points(x=xH, y=yH, col=colH, pch=pchH, cex=inp) points(xE, yE, pch=pchE, cex=outp); points(x=xE, y=yE, col=colE, pch=pchE, cex=inp) points(x=xMT, y=yMT, pch=pchMT, cex=outp); points(x=xMT, y=yMT, col=colMT, pch=pchMT, cex=inp) points(x=xF, y=yF, pch=pchF, cex=outp); points(x=xF, y=yF, col=colF, pch=pchF, cex=inp) axis(side=1, at=x, labels=NA) ls=.5; mtext(side=1, at=x, line=ls, text=c('edge','above1','above2'), cex=cexax) mtext('Lambda from matrix models', side=2, line=1.9, cex=cexlb, at=0) mtext('NK transect', side=3, line=0, cex=cexlb) lines(x=c(abxNK,abxNK), y=c(3.8,lymax), lwd=5, col='white') ylb=4; text(x=xmin, y=ylb, 'Natural', cex=cexlb+.1, pos=4) # Lambda HB par(mar=marHB); xlimHB=c(xmin,xmaxHB) xL=x+jL*j; yL=nrr$lambda[nrr$Source=='L' & nrr$transect=='HB' & nrr$dorm=='av'] plot(xL, yL, pch=pchL, cex=outp, ylim=ylim, xlim=xlimHB, xaxt="n", yaxt="n", ylab='', xlab="", las=1) abline(h=1, col='grey45', lwd=ab, lty=ltyh) # line showing lambda=1 = get in early so points are on top abline(v=abxHB, col='grey45', lwd=ab, lty = 2) # line showing RL xM=x+jM*j; yM=nrr$lambda[nrr$Source=='M' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='TR'] xH=x+jH*j; yH=nrr$lambda[nrr$Source=='H' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='TR'] xE=2+jE*j; yE=nrr$lambda[nrr$Source=='E' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='TR'] xMT=x+jMT*j; yMT=nrr$lambda[nrr$Source=='MT' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='TR'] xF=x+jF*j; yF=nrr$lambda[nrr$Source=='F' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='TR'] lines(xL[HB1:n], yL[HB1:n], col=colL, lwd=l); lines(xM[HB1:n], yM[HB1:n], col=colM, lwd=l); lines(x=xH[HB1:n], yH[HB1:n], col=colH, lwd=l) #lines for L,M,H lines(xMT, yMT, col=colMT, lwd=l); lines(xF, yF, col=colF, lwd=l) #lines for extra sources points(xL, yL, col=colL, pch=pchL, cex=inp) points(x=xM, yM, pch=pchM, cex=outp); points(x=xM, yM, col=colM, pch=pchM, cex=inp) points(x=xH, y=yH, pch=pchH, cex=outp); points(x=xH, y=yH, col=colH, pch=pchH, cex=inp) points(xE, yE, pch=pchE, cex=outp); points(x=xE, y=yE, col=colE, pch=pchE, cex=inp) points(x=xMT, y=yMT, pch=pchMT, cex=outp); points(x=xMT, y=yMT, col=colMT, pch=pchMT, cex=inp) points(x=xF, y=yF, pch=pchF, cex=outp); points(x=xF, y=yF, col=colF, pch=pchF, cex=inp) axis(side=1, at=x, labels=NA, cex.axis=cexax) axis(side=2, labels=NA,tck=-.02) mtext(side=1, at=x, line=ls, text=c('high','edge','above'), cex=cexax) mtext('HB transect', side=3, line=0, cex=cexlb) #legend - first switch back to normal scale xleg=c(rep(2.8,6)); yleg=c(seq(2.6, to=4.2, length.out=6)-.1); points(x=xleg, y=yleg, pch=c(pchL,pchM,pchH,pchE,pchMT,pchF), cex=c(rep(outp,6))) points(x=xleg, y=yleg, pch=c(pchL,pchM,pchH,pchE,pchMT,pchF), cex=c(rep(inp,5)), col=c(colL,colM,colH,colE,colMT,colF)) text(x=xleg+.04, y=yleg-.05, labels=c('Low','Mid', 'High', 'Edge','MT-High','FT-high'), pos=4, cex=cexst) #NK OTCs - lambda from matrix model ### x=c(1,2.7); x2=3.3; outpO=outp+.2; inpO=inp-.1; ylim=c(lymin,lymax); xlimNK=c(xmin, xmaxNK); par(mar=marNK); xM=x+jM*j; yM=nrr$lambda[nrr$Source=='M' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] plot(xM, yM, pch=pchM, cex=outpO, ylim=ylim, xlim=xlimNK, xaxt="n", ylab='', xlab="", las=1, col=colOTC) abline(h=1, col='grey45', lwd=ab, lty=ltyh) # line showing lambda=1 = get in early so points are on top abline(v=abxNK, col='grey45', lwd=ab, lty = 2) # line showing RL xMb=x2+jM*j; yMb=nrr$lambda[nrr$Source=='M' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xMb, yMb, pch=pchM, cex=outpO, col=colOTC) xH=x+jH*j; yH=nrr$lambda[nrr$Source=='H' & nrr$transect=='NK' & nrr$dorm=='av'& nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] points(xH, yH, pch=pchH, cex=outpO, col=colOTC) xHb=x2+jH*j; yHb=nrr$lambda[nrr$Source=='H' & nrr$transect=='NK' & nrr$dorm=='av'& nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xHb, yHb, pch=pchH, cex=outpO, col=colOTC) xMT=x+jMT*j; yMT=nrr$lambda[nrr$Source=='MT' & nrr$transect=='NK' & nrr$dorm=='av'& nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] points(xMT, yMT, pch=pchMT, cex=outpO, col=colOTC) xMTb=x2+jMT*j; yMTb=nrr$lambda[nrr$Source=='MT' & nrr$transect=='NK' & nrr$dorm=='av' & nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xMTb, yMTb, pch=pchMT, cex=outpO, col=colOTC) lines(c(xM,xMb), c(yM,yMb), col=colM, lwd=l); lines(c(xH,xHb), c(yH,yHb), col=colH, lwd=l); lines(c(xMT,xMTb), c(yMT,yMTb), col=colMT, lwd=l); points(x=xM, yM, col=colM, pch=pchM, cex=inpO); points(x=xMb, yMb, col=colM, pch=pchM, cex=inpO) points(x=xH, y=yH, col=colH, pch=pchH, cex=inpO); points(x=xHb, y=yHb, col=colH, pch=pchH, cex=inpO) points(x=xMT, y=yMT, col=colMT, pch=pchMT, cex=inpO); points(x=xMTb, y=yMTb, col=colMT, pch=pchMT, cex=inpO) axis(side=1, at=c(1,3), labels=NA, cex.axis=cexax) mtext(side=1, at=c(1,3), line=ls, text=c('edge','above2'), cex=cexax) pEad=.1; lpE=1.3; mtext(side=1, at=c(x-pEad,x2+pEad), line=lpE, text=c('','pEm from','pEm from'), cex=cexax-.2) mtext(side=1, at=c(x-pEad,x2+pEad), line=lpE+.5, text=c('','above2','edge'), cex=cexax-.2) lines(x=c(abxNK,abxNK), y=c(3.8,lymax), lwd=5, col='white') text(x=xmin, y=ylb, 'Warmed', cex=cexlb+.1, pos=4, col=colOTC) # Lambda HB par(mar=marHB); xlimHB=c(xmin,xmaxHB) xM=x+jM*j; yM=nrr$lambda[nrr$Source=='M' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] plot(xM, yM, pch=pchM, cex=outpO, ylim=ylim, xlim=xlimNK, xaxt="n",yaxt="n", ylab='', xlab="", las=1, col=colOTC) abline(h=1, col='grey45', lwd=ab, lty=ltyh) # line showing lambda=1 = get in early so points are on top abline(v=abxNK, col='grey45', lwd=ab, lty = 2) # line showing RL xMb=x2+jM*j; yMb=nrr$lambda[nrr$Source=='M' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xMb, yMb, pch=pchM, cex=outpO, col=colOTC) xH=x+jH*j; yH=nrr$lambda[nrr$Source=='H' & nrr$transect=='HB' & nrr$dorm=='av'& nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] points(xH, yH, pch=pchH, cex=outpO, col=colOTC) xHb=x2+jH*j; yHb=nrr$lambda[nrr$Source=='H' & nrr$transect=='HB' & nrr$dorm=='av'& nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xHb, yHb, pch=pchH, cex=outpO, col=colOTC) xMT=x+jMT*j; yMT=nrr$lambda[nrr$Source=='MT' & nrr$transect=='HB' & nrr$dorm=='av'& nrr$treat=='OTC' & nrr$pEmfrom=='TR&OTCplots at same site as P1s'] points(xMT, yMT, pch=pchMT, cex=outpO, col=colOTC) xMTb=x2+jMT*j; yMTb=nrr$lambda[nrr$Source=='MT' & nrr$transect=='HB' & nrr$dorm=='av' & nrr$pEmfrom!='TR&OTCplots at same site as P1s'] points(xMTb, yMTb, pch=pchMT, cex=outpO, col=colOTC) lines(c(xM,xMb), c(yM,yMb), col=colM, lwd=l); lines(c(xH,xHb), c(yH,yHb), col=colH, lwd=l); lines(c(xMT,xMTb), c(yMT,yMTb), col=colMT, lwd=l); points(x=xM, yM, col=colM, pch=pchM, cex=inpO); points(x=xMb, yMb, col=colM, pch=pchM, cex=inpO) points(x=xH, y=yH, col=colH, pch=pchH, cex=inpO); points(x=xHb, y=yHb, col=colH, pch=pchH, cex=inpO) points(x=xMT, y=yMT, col=colMT, pch=pchMT, cex=inpO); points(x=xMTb, y=yMTb, col=colMT, pch=pchMT, cex=inpO) axis(side=2, labels=NA,tck=-.02) axis(side=1, at=c(1,3), labels=NA, cex.axis=cexax) mtext(side=1, at=c(1,3), line=.5, text=c('high','above'), cex=cexax) mtext(side=1, at=c(x-pEad,x2+pEad), line=lpE, text=c('','pEm from','pEm from'), cex=cexax-.2) mtext(side=1, at=c(x-pEad,x2+pEad), line=lpE+.5, text=c('','above','edge'), cex=cexax-.2) mtext('Site', side=1, line=3, at=.6, cex=cexlb)