## Population Projections for LUCO7 using vital rates from ## Monitoring data in Saddle, Forest, and Red sites ## load variables yrspan <- 50 nptotal <- 42 npsite <- 14 dimtotal <- 15 dimsite <- 5 sites <- 3 runs <- 1000 tmax <- 50 bootreps <- 500 ## load data files corrdata <- read.table("MScorrdata.csv", header=TRUE, sep=",") NPvitals <- read.table("MSvitals.csv", header=TRUE, sep=",") FPvitals <- read.table("FPvitals.csv", header=TRUE, sep=",") RPvitals <- read.table("RPvitals.csv", header=TRUE, sep=",") NPvrmeans <- NPvitals$mean NPvrvars <- NPvitals$var vrtypes <- NPvitals$type vrmins <- NPvitals$min vrmaxs <- NPvitals$max FPvrmeans <-FPvitals$mean FPvrvars <- FPvitals$var RPvrmeans <- RPvitals$mean RPvrvars <- RPvitals$var ## load popbio package library(popbio) ############################################################## ## define the transition matrices MSmxdef <- function(vrs) { matrix(c(vrs[2]*vrs[7], 0, 0, vrs[6]*vrs[2]*vrs[7], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[2]*(1-vrs[7]), 0, 0, vrs[6]*(1-vrs[7])*(vrs[2]^(9/12))*vrs[1], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[3]*(1-vrs[8]), vrs[4]*(1-vrs[9]), vrs[5]*(1-vrs[11]), vrs[4]*(1-vrs[13]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[3]*vrs[8], vrs[4]*(vrs[9]-vrs[10]), vrs[5]*(vrs[11]-vrs[12]), vrs[4]*(vrs[13]-vrs[14]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[4]*vrs[10], vrs[5]*vrs[12], vrs[4]*vrs[14], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[16]*vrs[21], 0, 0, vrs[20]*vrs[16]*vrs[21], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[16]*(1-vrs[21]), 0, 0, vrs[20]*(1-vrs[21])*(vrs[16]^(9/12))*vrs[15], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[17]*(1-vrs[22]), vrs[18]*(1-vrs[23]), vrs[19]*(1-vrs[25]), vrs[18]*(1-vrs[27]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[17]*vrs[22], vrs[18]*(vrs[23]-vrs[24]), vrs[19]*(vrs[25]-vrs[26]), vrs[18]*(vrs[27]-vrs[28]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[18]*vrs[24], vrs[19]*vrs[26], vrs[18]*vrs[28], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[30]*vrs[35], 0, 0, vrs[34]*vrs[30]*vrs[35], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[30]*(1-vrs[35]), 0, 0, vrs[34]*(1-vrs[35])*(vrs[30]^(9/12))*vrs[29], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[31]*(1-vrs[36]), vrs[32]*(1-vrs[37]), vrs[33]*(1-vrs[39]), vrs[32]*(1-vrs[41]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[31]*vrs[36], vrs[32]*(vrs[37]-vrs[38]), vrs[33]*(vrs[39]-vrs[40]), vrs[32]*(vrs[41]-vrs[42]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, vrs[32]*vrs[38], vrs[33]*vrs[40], vrs[32]*vrs[42]), nrow = 15, byrow = TRUE) } ## Forest site matrix definition Fmxdef <- function(vrs) { matrix(c(vrs[2]*vrs[7], 0, 0, vrs[6]*vrs[2]*vrs[7], 0, vrs[2]*(1-vrs[7]), 0, 0, vrs[6]*(1-vrs[7])*(vrs[2]^(9/12))*vrs[1], 0, 0, vrs[3]*(1-vrs[8]), vrs[4]*(1-vrs[9]), vrs[5]*(1-vrs[11]), vrs[4]*(1-vrs[13]), 0, vrs[3]*vrs[8], vrs[4]*(vrs[9]-vrs[10]), vrs[5]*(vrs[11]-vrs[12]), vrs[4]*(vrs[13]-vrs[14]), 0, 0, vrs[4]*vrs[10], vrs[5]*vrs[12], vrs[4]*vrs[14]), nrow = 5, byrow = TRUE) } ## Red site matrix definition Rmxdef <- function(vrs) { matrix(c(vrs[16]*vrs[21], 0, 0, vrs[20]*vrs[16]*vrs[21], 0, vrs[16]*(1-vrs[21]), 0, 0, vrs[20]*(1-vrs[21])*(vrs[16]^(9/12))*vrs[15], 0, 0, vrs[17]*(1-vrs[22]), vrs[18]*(1-vrs[23]), vrs[19]*(1-vrs[25]), vrs[18]*(1-vrs[27]), 0, vrs[17]*vrs[22], vrs[18]*(vrs[23]-vrs[24]), vrs[19]*(vrs[25]-vrs[26]), vrs[18]*(vrs[27]-vrs[28]), 0, 0, vrs[18]*vrs[24], vrs[19]*vrs[26], vrs[18]*vrs[28]), nrow = 5, byrow = TRUE) } ## Saddle site matrix definition Smxdef <- function(vrs) { matrix(c(vrs[30]*vrs[35], 0, 0, vrs[34]*vrs[30]*vrs[35], 0, vrs[30]*(1-vrs[35]), 0, 0, vrs[34]*(1-vrs[35])*(vrs[30]^(9/12))*vrs[29], 0, 0, vrs[31]*(1-vrs[36]), vrs[32]*(1-vrs[37]), vrs[33]*(1-vrs[39]), vrs[32]*(1-vrs[41]), 0, vrs[31]*vrs[36], vrs[32]*(vrs[37]-vrs[38]), vrs[33]*(vrs[39]-vrs[40]), vrs[32]*(vrs[41]-vrs[42]), 0, 0, vrs[32]*vrs[38], vrs[33]*vrs[40], vrs[32]*vrs[42]), nrow = 5, byrow = TRUE) } ############################################################ ## make the multi-site M12 matrix MScorrE <- cor(corrdata, use="pairwise.complete.obs", method="spearman") for(i in 1:84){ for(j in 1:84){ if(is.na(MScorrE[i,j])){MScorrE[i,j]=0} } } diag(MScorrE)=1 eigen(MScorrE)$values MScorrC <- matrix(MScorrE[1:42,1:42], nrow=42) MScorrB <- matrix(MScorrE[43:84,1:42], nrow=42) MScorrM <- matrix(0, nrow=yrspan*nptotal, ncol=yrspan*nptotal) for(i in 1:yrspan){ for(j in 1:yrspan){ if(i==j){mxpiece=MScorrC; expo=1} if(i>j){expo=i-j; mxpiece=(MScorrB)^expo} if(ij){expo=i-j; mxpiece=(FcorrB)^expo} if(ij){expo=i-j; mxpiece=(RcorrB)^expo} if(ij){expo=i-j; mxpiece=(ScorrB)^expo} if(imx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen001PrExt[b,t] <- Scen001PrExt[b,t]+1; break} } Scen001sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ########################################################## ##Scenario 002: Multi-site, No Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3,750,15,32,13,1,2735,117,24,72,6) weights <- c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0) nx <- 30 Scen002PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen002sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) MSm <- c(MSm[(nptotal+1):MSsz],rnorm(nptotal)) MSnewxy <- MSznew%*%MSm MSoldxy <- MSnewxy for(y in 1:nptotal){ index = round(100*pnorm(MSnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- MSmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:15){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1||g==6||g==11){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen002PrExt[b,t] <- Scen002PrExt[b,t]+1; break} } Scen002sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################### ##Scenario 007: Multi-site, Full Predation, Seed Bank Estimate, QET Adult n0 <- c(278,24,23,24,3,750,15,32,13,1,2735,117,24,72,6) weights <- c(0,0,1,1,0,0,0,1,1,0,0,0,1,1,0) nx <- 10 Scen007PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen007sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) MSm <- c(MSm[(nptotal+1):MSsz],rnorm(nptotal)) MSnewxy <- MSznew%*%MSm MSoldxy <- MSnewxy for(y in 1:nptotal){ index = round(100*pnorm(MSnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- MSmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:15){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1||g==6||g==11){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen007PrExt[b,t] <- Scen007PrExt[b,t]+1; break} } Scen007sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################### ##Scenario 008: Multi-site, Full Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3,750,15,32,13,1,2735,117,24,72,6) weights <- c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0) nx <- 30 Scen008PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen008sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) MSm <- c(MSm[(nptotal+1):MSsz],rnorm(nptotal)) MSnewxy <- MSznew%*%MSm MSoldxy <- MSnewxy for(y in 1:nptotal){ index = round(100*pnorm(MSnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- MSmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:15){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1||g==6||g==11){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen008PrExt[b,t] <- Scen008PrExt[b,t]+1; break} } Scen008sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################### ##Scenario 013: Multi-site, Real Predation, Seed Bank Estimate, QET Adult n0 <- c(278,24,23,24,3,750,15,32,13,1,2735,117,24,72,6) weights <- c(0,0,1,1,0,0,0,1,1,0,0,0,1,1,0) nx <- 10 Scen013PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen013sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) MSm <- c(MSm[(nptotal+1):MSsz],rnorm(nptotal)) MSnewxy <- MSznew%*%MSm MSoldxy <- MSnewxy for(y in 1:nptotal){ index = round(100*pnorm(MSnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- MSmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:15){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1||g==6||g==11){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen013PrExt[b,t] <- Scen013PrExt[b,t]+1; break} } Scen013sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################### ##Scenario 014: Multi-site, Real Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3,750,15,32,13,1,2735,117,24,72,6) weights <- c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0) nx <- 30 Scen014PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen014sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) MSm <- c(MSm[(nptotal+1):MSsz],rnorm(nptotal)) MSnewxy <- MSznew%*%MSm MSoldxy <- MSnewxy for(y in 1:nptotal){ index = round(100*pnorm(MSnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- MSmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:15){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1||g==6||g==11){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2||g==7||g==12){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3||g==8||g==13){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4||g==9||g==14){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen014PrExt[b,t] <- Scen014PrExt[b,t]+1; break} } Scen014sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 019: Forest, No Predation, Seed Bank Estimate, QET Adult n0 <- c(278,24,23,24,3) weights <- c(0,0,1,1,0) nx <- 3 Scen019PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen019sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen019PrExt[b,t] <- Scen019PrExt[b,t]+1; break} } Scen019sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 020: Forest, No Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3) weights <- c(1,0,0,0,0) nx <- 10 Scen020PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen020sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen020PrExt[b,t] <- Scen020PrExt[b,t]+1; break} } Scen020sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 025: Forest, Full Predation, Seed Bank Estimate, QET Adult n0 <- c(278,24,23,24,3) weights <- c(0,0,1,1,0) nx <- 3 Scen025PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen025sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen025PrExt[b,t] <- Scen025PrExt[b,t]+1; break} } Scen025sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 026: Forest, Full Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3) weights <- c(1,0,0,0,0) nx <- 10 Scen026PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen026sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen026PrExt[b,t] <- Scen026PrExt[b,t]+1; break} } Scen026sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 031: Forest, Real Predation, Seed Bank Estimate, QET Adult n0 <- c(278,24,23,24,3) weights <- c(0,0,1,1,0) nx <- 3 Scen031PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen031sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen031PrExt[b,t] <- Scen031PrExt[b,t]+1; break} } Scen031sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 032: Forest, Real Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(278,24,23,24,3) weights <- c(1,0,0,0,0) nx <- 10 Scen032PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen032sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=14) Fm <- c(Fm[(npsite+1):Fsz],rnorm(npsite)) Fnewxy <- Fznew%*%Fm Foldxy <- Fnewxy for(y in 1:14){ index = round(100*pnorm(Fnewxy[y])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- Fmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen032PrExt[b,t] <- Scen032PrExt[b,t]+1; break} } Scen032sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 037: Red, No Predation, Seed Bank Estimate, QET Adult n0 <- c(750,15,32,13,1) weights <- c(0,0,1,1,0) nx <- 3 Scen037PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen037sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Rm <- c(Rm[(npsite+1):Rsz],rnorm(npsite)) Rnewxy <- Rznew%*%Rm Roldxy <- Rnewxy for(y in 15:28){ index = round(100*pnorm(Rnewxy[y-14])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Rmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen037PrExt[b,t] <- Scen037PrExt[b,t]+1; break} } Scen037sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 038: Red, No Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(750,15,32,13,1) weights <- c(1,0,0,0,0) nx <- 10 Scen038PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen038sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Rm <- c(Rm[(npsite+1):Rsz],rnorm(npsite)) Rnewxy <- Rznew%*%Rm Roldxy <- Rnewxy for(y in 15:28){ index = round(100*pnorm(Rnewxy[y-14])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Rmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen038PrExt[b,t] <- Scen038PrExt[b,t]+1; break} } Scen038sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 043: Red, Full Predation, Seed Bank Estimate, QET Adult n0 <- c(750,15,32,13,1) weights <- c(0,0,1,1,0) nx <- 3 Scen043PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen043sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Rm <- c(Rm[(npsite+1):Rsz],rnorm(npsite)) Rnewxy <- Rznew%*%Rm Roldxy <- Rnewxy for(y in 15:28){ index = round(100*pnorm(Rnewxy[y-14])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Rmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen043PrExt[b,t] <- Scen043PrExt[b,t]+1; break} } Scen043sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 044: Red, Full Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(750,15,32,13,1) weights <- c(1,0,0,0,0) nx <- 10 Scen044PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen044sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Rm <- c(Rm[(npsite+1):Rsz],rnorm(npsite)) Rnewxy <- Rznew%*%Rm Roldxy <- Rnewxy for(y in 15:28){ index = round(100*pnorm(Rnewxy[y-14])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Rmxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen044PrExt[b,t] <- Scen044PrExt[b,t]+1; break} } Scen044sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 055: Saddle, No Predation, Seed Bank Estimate, QET Adult n0 <- c(2735,117,24,72,6) weights <- c(0,0,1,1,0) nx <- 3 Scen055PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen055sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen055PrExt[b,t] <- Scen055PrExt[b,t]+1; break} } Scen055sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 056: Saddle, No Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(2735,117,24,72,6) weights <- c(1,0,0,0,0) nx <- 10 Scen056PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen056sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = NPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen056PrExt[b,t] <- Scen056PrExt[b,t]+1; break} } Scen056sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 061: Saddle, Full Predation, Seed Bank Estimate, QET Adult n0 <- c(2735,117,24,72,6) weights <- c(0,0,1,1,0) nx <- 3 Scen061PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen061sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen061PrExt[b,t] <- Scen061PrExt[b,t]+1; break} } Scen061sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 062: Saddle, Full Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(2735,117,24,72,6) weights <- c(1,0,0,0,0) nx <- 10 Scen062PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen062sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = FPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen062PrExt[b,t] <- Scen062PrExt[b,t]+1; break} } Scen062sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 067: Saddle, Real Predation, Seed Bank Estimate, QET Adult n0 <- c(2735,117,24,72,6) weights <- c(0,0,1,1,0) nx <- 3 Scen067PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen067sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen067PrExt[b,t] <- Scen067PrExt[b,t]+1; break} } Scen067sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ################################################################# ##Scenario 068: Saddle, Real Predation, Seed Bank Estimate, QET Seed Bank n0 <- c(2735,117,24,72,6) weights <- c(1,0,0,0,0) nx <- 10 Scen068PrExt <- matrix(0, nrow=bootreps, ncol=tmax) Scen068sLam <- matrix(0, nrow=bootreps, ncol=runs) for(b in 1:bootreps){ for(x in 1:runs){ nt <- n0 for(t in 1:tmax){ vrs <- rep(0, each=42) Sm <- c(Sm[(npsite+1):Ssz],rnorm(npsite)) Snewxy <- Sznew%*%Sm Soldxy <- Snewxy for(y in 29:42){ index = round(100*pnorm(Snewxy[y-28])) if(index==0){index=1} if(index==100){index=99} vrs[y] = RPbetas[index,y] } mx <- Smxdef(vrs) mx[mx<0] <-0 nt2 <- mx%*%nt for(g in 1:5){ if(nt[g]<=50){q<-runif(round(nt[g])) if(g==1){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g]))} else{if(g==2){ nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q<=mx[g+1,g]); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>mx[g+1,g]&q<=(mx[g+1,g]+mx[g+2,g]))} else{if(g==3){ nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q<=mx[g,g]); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>mx[g,g]&q<=(mx[g,g]+mx[g+1,g])); nt2[g+2]<-nt2[g+2]-(nt[g]*mx[g+2,g])+sum(q>(mx[g,g]+mx[g+1,g])&q<=(mx[g,g]+mx[g+1,g]+mx[g+2,g]))} else{if(g==4){ nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q<=mx[g-1,g]); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>mx[g-1,g]&q<=(mx[g-1,g]+mx[g,g])); nt2[g+1]<-nt2[g+1]-(nt[g]*mx[g+1,g])+sum(q>(mx[g-1,g]+mx[g,g])&q<=(mx[g-1,g]+mx[g,g]+mx[g+1,g]))} else{nt2[g-2]<-nt2[g-2]-(nt[g]*mx[g-2,g])+sum(q<=mx[g-2,g]); nt2[g-1]<-nt2[g-1]-(nt[g]*mx[g-1,g])+sum(q>mx[g-2,g]&q<=(mx[g-2,g]+mx[g-1,g])); nt2[g]<-nt2[g]-(nt[g]*mx[g,g])+sum(q>(mx[g-2,g]+mx[g-1,g])&q<=(mx[g-2,g]+mx[g-1,g]+mx[g,g]))} } } } } } nt <- nt2 if(sum(weights*nt)<=nx){ Scen068PrExt[b,t] <- Scen068PrExt[b,t]+1; break} } Scen068sLam[b,x] <- (sum(nt)/sum(n0))^(1/tmax) } } ############################################################ (CDF001 <- (cumsum(Scen001PrExt))/1000) (CDF002 <- (cumsum(Scen002PrExt))/1000) (CDF003 <- (cumsum(Scen003PrExt))/1000) (CDF004 <- (cumsum(Scen004PrExt))/1000) (CDF005 <- (cumsum(Scen005PrExt))/1000) (CDF006 <- (cumsum(Scen006PrExt))/1000) (CDF007 <- (cumsum(Scen007PrExt))/1000) (CDF008 <- (cumsum(Scen008PrExt))/1000) (CDF009 <- (cumsum(Scen009PrExt))/1000) (CDF010 <- (cumsum(Scen010PrExt))/1000) (CDF011 <- (cumsum(Scen011PrExt))/1000) (CDF012 <- (cumsum(Scen012PrExt))/1000) (CDF013 <- (cumsum(Scen013PrExt))/1000) (CDF014 <- (cumsum(Scen014PrExt))/1000) (CDF015 <- (cumsum(Scen015PrExt))/1000) (CDF016 <- (cumsum(Scen016PrExt))/1000) (CDF017 <- (cumsum(Scen017PrExt))/1000) (CDF018 <- (cumsum(Scen018PrExt))/1000) (CDF019 <- (cumsum(Scen019PrExt))/1000) (CDF020 <- (cumsum(Scen020PrExt))/1000) (CDF021 <- (cumsum(Scen021PrExt))/1000) (CDF022 <- (cumsum(Scen022PrExt))/1000) (CDF023 <- (cumsum(Scen023PrExt))/1000) (CDF024 <- (cumsum(Scen024PrExt))/1000) (CDF025 <- (cumsum(Scen025PrExt))/1000) (CDF026 <- (cumsum(Scen026PrExt))/1000) (CDF027 <- (cumsum(Scen027PrExt))/1000) (CDF028 <- (cumsum(Scen028PrExt))/1000) (CDF029 <- (cumsum(Scen029PrExt))/1000) (CDF030 <- (cumsum(Scen030PrExt))/1000) (CDF031 <- (cumsum(Scen031PrExt))/1000) (CDF032 <- (cumsum(Scen032PrExt))/1000) (CDF033 <- (cumsum(Scen033PrExt))/1000) (CDF034 <- (cumsum(Scen034PrExt))/1000) (CDF035 <- (cumsum(Scen035PrExt))/1000) (CDF036 <- (cumsum(Scen036PrExt))/1000) (CDF037 <- (cumsum(Scen037PrExt))/1000) (CDF038 <- (cumsum(Scen038PrExt))/1000) (CDF039 <- (cumsum(Scen039PrExt))/1000) (CDF040 <- (cumsum(Scen040PrExt))/1000) (CDF041 <- (cumsum(Scen041PrExt))/1000) (CDF042 <- (cumsum(Scen042PrExt))/1000) (CDF043 <- (cumsum(Scen043PrExt))/1000) (CDF044 <- (cumsum(Scen044PrExt))/1000) (CDF045 <- (cumsum(Scen045PrExt))/1000) (CDF046 <- (cumsum(Scen046PrExt))/1000) (CDF047 <- (cumsum(Scen047PrExt))/1000) (CDF048 <- (cumsum(Scen048PrExt))/1000) (CDF055 <- (cumsum(Scen055PrExt))/1000) (CDF056 <- (cumsum(Scen056PrExt))/1000) (CDF057 <- (cumsum(Scen057PrExt))/1000) (CDF058 <- (cumsum(Scen058PrExt))/1000) (CDF059 <- (cumsum(Scen059PrExt))/1000) (CDF060 <- (cumsum(Scen060PrExt))/1000) (CDF061 <- (cumsum(Scen061PrExt))/1000) (CDF062 <- (cumsum(Scen062PrExt))/1000) (CDF063 <- (cumsum(Scen063PrExt))/1000) (CDF064 <- (cumsum(Scen064PrExt))/1000) (CDF065 <- (cumsum(Scen065PrExt))/1000) (CDF066 <- (cumsum(Scen066PrExt))/1000) (CDF067 <- (cumsum(Scen067PrExt))/1000) (CDF068 <- (cumsum(Scen068PrExt))/1000) (CDF069 <- (cumsum(Scen069PrExt))/1000) (CDF070 <- (cumsum(Scen070PrExt))/1000) (CDF071 <- (cumsum(Scen071PrExt))/1000) (CDF072 <- (cumsum(Scen072PrExt))/1000) CDFs <- matrix(c(CDF001, CDF002, CDF003, CDF004, CDF005, CDF006, CDF007, CDF008, CDF009, CDF010, CDF011, CDF012, CDF013, CDF014, CDF015, CDF016, CDF017, CDF018, CDF019, CDF020, CDF021, CDF022, CDF023, CDF024, CDF025, CDF026, CDF027, CDF028, CDF029, CDF030, CDF031, CDF032, CDF033, CDF034, CDF035, CDF036, CDF037, CDF038, CDF039, CDF040, CDF041, CDF042, CDF043, CDF044, CDF045, CDF046, CDF047, CDF048, CDF055, CDF056, CDF057, CDF058, CDF059, CDF060, CDF061, CDF062, CDF063, CDF064, CDF065, CDF066, CDF067, CDF068, CDF069, CDF070, CDF071, CDF072), nrow=50) write.csv(CDFs, "CDFs.csv", row.names=FALSE, col.names=FALSE) ######################################################################## Lam001 <- c(mean(Scen001sLam), mean(Scen001sLam)+((sd(Scen001sLam)/sqrt(1000))*1.96), mean(Scen001sLam)-((sd(Scen001sLam)/sqrt(1000))*1.96)) Lam002 <- c(mean(Scen002sLam), mean(Scen002sLam)+((sd(Scen002sLam)/sqrt(1000))*1.96), mean(Scen002sLam)-((sd(Scen002sLam)/sqrt(1000))*1.96)) Lam003 <- c(mean(Scen003sLam), mean(Scen003sLam)+((sd(Scen003sLam)/sqrt(1000))*1.96), mean(Scen003sLam)-((sd(Scen003sLam)/sqrt(1000))*1.96)) Lam004 <- c(mean(Scen004sLam), mean(Scen004sLam)+((sd(Scen004sLam)/sqrt(1000))*1.96), mean(Scen004sLam)-((sd(Scen004sLam)/sqrt(1000))*1.96)) Lam005 <- c(mean(Scen005sLam), mean(Scen005sLam)+((sd(Scen005sLam)/sqrt(1000))*1.96), mean(Scen005sLam)-((sd(Scen005sLam)/sqrt(1000))*1.96)) Lam006 <- c(mean(Scen006sLam), mean(Scen006sLam)+((sd(Scen006sLam)/sqrt(1000))*1.96), mean(Scen006sLam)-((sd(Scen006sLam)/sqrt(1000))*1.96)) Lam007 <- c(mean(Scen007sLam), mean(Scen007sLam)+((sd(Scen007sLam)/sqrt(1000))*1.96), mean(Scen007sLam)-((sd(Scen007sLam)/sqrt(1000))*1.96)) Lam008 <- c(mean(Scen008sLam), mean(Scen008sLam)+((sd(Scen008sLam)/sqrt(1000))*1.96), mean(Scen008sLam)-((sd(Scen008sLam)/sqrt(1000))*1.96)) Lam009 <- c(mean(Scen009sLam), mean(Scen009sLam)+((sd(Scen009sLam)/sqrt(1000))*1.96), mean(Scen009sLam)-((sd(Scen009sLam)/sqrt(1000))*1.96)) Lam010 <- c(mean(Scen010sLam), mean(Scen010sLam)+((sd(Scen010sLam)/sqrt(1000))*1.96), mean(Scen010sLam)-((sd(Scen010sLam)/sqrt(1000))*1.96)) Lam011 <- c(mean(Scen011sLam), mean(Scen011sLam)+((sd(Scen011sLam)/sqrt(1000))*1.96), mean(Scen011sLam)-((sd(Scen011sLam)/sqrt(1000))*1.96)) Lam012 <- c(mean(Scen012sLam), mean(Scen012sLam)+((sd(Scen012sLam)/sqrt(1000))*1.96), mean(Scen012sLam)-((sd(Scen012sLam)/sqrt(1000))*1.96)) Lam013 <- c(mean(Scen013sLam), mean(Scen013sLam)+((sd(Scen013sLam)/sqrt(1000))*1.96), mean(Scen013sLam)-((sd(Scen013sLam)/sqrt(1000))*1.96)) Lam014 <- c(mean(Scen014sLam), mean(Scen014sLam)+((sd(Scen014sLam)/sqrt(1000))*1.96), mean(Scen014sLam)-((sd(Scen014sLam)/sqrt(1000))*1.96)) Lam015 <- c(mean(Scen015sLam), mean(Scen015sLam)+((sd(Scen015sLam)/sqrt(1000))*1.96), mean(Scen015sLam)-((sd(Scen015sLam)/sqrt(1000))*1.96)) Lam016 <- c(mean(Scen016sLam), mean(Scen016sLam)+((sd(Scen016sLam)/sqrt(1000))*1.96), mean(Scen016sLam)-((sd(Scen016sLam)/sqrt(1000))*1.96)) Lam017 <- c(mean(Scen017sLam), mean(Scen017sLam)+((sd(Scen017sLam)/sqrt(1000))*1.96), mean(Scen017sLam)-((sd(Scen017sLam)/sqrt(1000))*1.96)) Lam018 <- c(mean(Scen018sLam), mean(Scen018sLam)+((sd(Scen018sLam)/sqrt(1000))*1.96), mean(Scen018sLam)-((sd(Scen018sLam)/sqrt(1000))*1.96)) Lam019 <- c(mean(Scen019sLam), mean(Scen019sLam)+((sd(Scen019sLam)/sqrt(1000))*1.96), mean(Scen019sLam)-((sd(Scen019sLam)/sqrt(1000))*1.96)) Lam020 <- c(mean(Scen020sLam), mean(Scen020sLam)+((sd(Scen020sLam)/sqrt(1000))*1.96), mean(Scen020sLam)-((sd(Scen020sLam)/sqrt(1000))*1.96)) Lam021 <- c(mean(Scen021sLam), mean(Scen021sLam)+((sd(Scen021sLam)/sqrt(1000))*1.96), mean(Scen021sLam)-((sd(Scen021sLam)/sqrt(1000))*1.96)) Lam022 <- c(mean(Scen022sLam), mean(Scen022sLam)+((sd(Scen022sLam)/sqrt(1000))*1.96), mean(Scen022sLam)-((sd(Scen022sLam)/sqrt(1000))*1.96)) Lam023 <- c(mean(Scen023sLam), mean(Scen023sLam)+((sd(Scen023sLam)/sqrt(1000))*1.96), mean(Scen023sLam)-((sd(Scen023sLam)/sqrt(1000))*1.96)) Lam024 <- c(mean(Scen024sLam), mean(Scen024sLam)+((sd(Scen024sLam)/sqrt(1000))*1.96), mean(Scen024sLam)-((sd(Scen024sLam)/sqrt(1000))*1.96)) Lam025 <- c(mean(Scen025sLam), mean(Scen025sLam)+((sd(Scen025sLam)/sqrt(1000))*1.96), mean(Scen025sLam)-((sd(Scen025sLam)/sqrt(1000))*1.96)) Lam026 <- c(mean(Scen026sLam), mean(Scen026sLam)+((sd(Scen026sLam)/sqrt(1000))*1.96), mean(Scen026sLam)-((sd(Scen026sLam)/sqrt(1000))*1.96)) Lam027 <- c(mean(Scen027sLam), mean(Scen027sLam)+((sd(Scen027sLam)/sqrt(1000))*1.96), mean(Scen027sLam)-((sd(Scen027sLam)/sqrt(1000))*1.96)) Lam028 <- c(mean(Scen028sLam), mean(Scen028sLam)+((sd(Scen028sLam)/sqrt(1000))*1.96), mean(Scen028sLam)-((sd(Scen028sLam)/sqrt(1000))*1.96)) Lam029 <- c(mean(Scen029sLam), mean(Scen029sLam)+((sd(Scen029sLam)/sqrt(1000))*1.96), mean(Scen029sLam)-((sd(Scen029sLam)/sqrt(1000))*1.96)) Lam030 <- c(mean(Scen030sLam), mean(Scen030sLam)+((sd(Scen030sLam)/sqrt(1000))*1.96), mean(Scen030sLam)-((sd(Scen030sLam)/sqrt(1000))*1.96)) Lam031 <- c(mean(Scen031sLam), mean(Scen031sLam)+((sd(Scen031sLam)/sqrt(1000))*1.96), mean(Scen031sLam)-((sd(Scen031sLam)/sqrt(1000))*1.96)) Lam032 <- c(mean(Scen032sLam), mean(Scen032sLam)+((sd(Scen032sLam)/sqrt(1000))*1.96), mean(Scen032sLam)-((sd(Scen032sLam)/sqrt(1000))*1.96)) Lam033 <- c(mean(Scen033sLam), mean(Scen033sLam)+((sd(Scen033sLam)/sqrt(1000))*1.96), mean(Scen033sLam)-((sd(Scen033sLam)/sqrt(1000))*1.96)) Lam034 <- c(mean(Scen034sLam), mean(Scen034sLam)+((sd(Scen034sLam)/sqrt(1000))*1.96), mean(Scen034sLam)-((sd(Scen034sLam)/sqrt(1000))*1.96)) Lam035 <- c(mean(Scen035sLam), mean(Scen035sLam)+((sd(Scen035sLam)/sqrt(1000))*1.96), mean(Scen035sLam)-((sd(Scen035sLam)/sqrt(1000))*1.96)) Lam036 <- c(mean(Scen036sLam), mean(Scen036sLam)+((sd(Scen036sLam)/sqrt(1000))*1.96), mean(Scen036sLam)-((sd(Scen036sLam)/sqrt(1000))*1.96)) Lam037 <- c(mean(Scen037sLam), mean(Scen037sLam)+((sd(Scen037sLam)/sqrt(1000))*1.96), mean(Scen037sLam)-((sd(Scen037sLam)/sqrt(1000))*1.96)) Lam038 <- c(mean(Scen038sLam), mean(Scen038sLam)+((sd(Scen038sLam)/sqrt(1000))*1.96), mean(Scen038sLam)-((sd(Scen038sLam)/sqrt(1000))*1.96)) Lam039 <- c(mean(Scen039sLam), mean(Scen039sLam)+((sd(Scen039sLam)/sqrt(1000))*1.96), mean(Scen039sLam)-((sd(Scen039sLam)/sqrt(1000))*1.96)) Lam040 <- c(mean(Scen040sLam), mean(Scen040sLam)+((sd(Scen040sLam)/sqrt(1000))*1.96), mean(Scen040sLam)-((sd(Scen040sLam)/sqrt(1000))*1.96)) Lam041 <- c(mean(Scen041sLam), mean(Scen041sLam)+((sd(Scen041sLam)/sqrt(1000))*1.96), mean(Scen041sLam)-((sd(Scen041sLam)/sqrt(1000))*1.96)) Lam042 <- c(mean(Scen042sLam), mean(Scen042sLam)+((sd(Scen042sLam)/sqrt(1000))*1.96), mean(Scen042sLam)-((sd(Scen042sLam)/sqrt(1000))*1.96)) Lam043 <- c(mean(Scen043sLam), mean(Scen043sLam)+((sd(Scen043sLam)/sqrt(1000))*1.96), mean(Scen043sLam)-((sd(Scen043sLam)/sqrt(1000))*1.96)) Lam044 <- c(mean(Scen044sLam), mean(Scen044sLam)+((sd(Scen044sLam)/sqrt(1000))*1.96), mean(Scen044sLam)-((sd(Scen044sLam)/sqrt(1000))*1.96)) Lam045 <- c(mean(Scen045sLam), mean(Scen045sLam)+((sd(Scen045sLam)/sqrt(1000))*1.96), mean(Scen045sLam)-((sd(Scen045sLam)/sqrt(1000))*1.96)) Lam046 <- c(mean(Scen046sLam, na.rm=TRUE), mean(Scen046sLam, na.rm=TRUE)+((sd(Scen046sLam, na.rm=TRUE)/sqrt(1000))*1.96), mean(Scen046sLam, na.rm=TRUE)-((sd(Scen046sLam, na.rm=TRUE)/sqrt(1000))*1.96)) Lam047 <- c(mean(Scen047sLam), mean(Scen047sLam)+((sd(Scen047sLam)/sqrt(1000))*1.96), mean(Scen047sLam)-((sd(Scen047sLam)/sqrt(1000))*1.96)) Lam048 <- c(mean(Scen048sLam), mean(Scen048sLam)+((sd(Scen048sLam)/sqrt(1000))*1.96), mean(Scen048sLam)-((sd(Scen048sLam)/sqrt(1000))*1.96)) Lam055 <- c(mean(Scen055sLam), mean(Scen055sLam)+((sd(Scen055sLam)/sqrt(1000))*1.96), mean(Scen055sLam)-((sd(Scen055sLam)/sqrt(1000))*1.96)) Lam056 <- c(mean(Scen056sLam), mean(Scen056sLam)+((sd(Scen056sLam)/sqrt(1000))*1.96), mean(Scen056sLam)-((sd(Scen056sLam)/sqrt(1000))*1.96)) Lam057 <- c(mean(Scen057sLam), mean(Scen057sLam)+((sd(Scen057sLam)/sqrt(1000))*1.96), mean(Scen057sLam)-((sd(Scen057sLam)/sqrt(1000))*1.96)) Lam058 <- c(mean(Scen058sLam), mean(Scen058sLam)+((sd(Scen058sLam)/sqrt(1000))*1.96), mean(Scen058sLam)-((sd(Scen058sLam)/sqrt(1000))*1.96)) Lam059 <- c(mean(Scen059sLam), mean(Scen059sLam)+((sd(Scen059sLam)/sqrt(1000))*1.96), mean(Scen059sLam)-((sd(Scen059sLam)/sqrt(1000))*1.96)) Lam060 <- c(mean(Scen060sLam), mean(Scen060sLam)+((sd(Scen060sLam)/sqrt(1000))*1.96), mean(Scen060sLam)-((sd(Scen060sLam)/sqrt(1000))*1.96)) Lam061 <- c(mean(Scen061sLam), mean(Scen061sLam)+((sd(Scen061sLam)/sqrt(1000))*1.96), mean(Scen061sLam)-((sd(Scen061sLam)/sqrt(1000))*1.96)) Lam062 <- c(mean(Scen062sLam), mean(Scen062sLam)+((sd(Scen062sLam)/sqrt(1000))*1.96), mean(Scen062sLam)-((sd(Scen062sLam)/sqrt(1000))*1.96)) Lam063 <- c(mean(Scen063sLam), mean(Scen063sLam)+((sd(Scen063sLam)/sqrt(1000))*1.96), mean(Scen063sLam)-((sd(Scen063sLam)/sqrt(1000))*1.96)) Lam064 <- c(mean(Scen064sLam), mean(Scen064sLam)+((sd(Scen064sLam)/sqrt(1000))*1.96), mean(Scen064sLam)-((sd(Scen064sLam)/sqrt(1000))*1.96)) Lam065 <- c(mean(Scen065sLam), mean(Scen065sLam)+((sd(Scen065sLam)/sqrt(1000))*1.96), mean(Scen065sLam)-((sd(Scen065sLam)/sqrt(1000))*1.96)) Lam066 <- c(mean(Scen066sLam, na.rm=TRUE), mean(Scen066sLam, na.rm=TRUE)+((sd(Scen066sLam, na.rm=TRUE)/sqrt(1000))*1.96), mean(Scen066sLam, na.rm=TRUE)-((sd(Scen066sLam, na.rm=TRUE)/sqrt(1000))*1.96)) Lam067 <- c(mean(Scen067sLam), mean(Scen067sLam)+((sd(Scen067sLam)/sqrt(1000))*1.96), mean(Scen067sLam)-((sd(Scen067sLam)/sqrt(1000))*1.96)) Lam068 <- c(mean(Scen068sLam), mean(Scen068sLam)+((sd(Scen068sLam)/sqrt(1000))*1.96), mean(Scen068sLam)-((sd(Scen068sLam)/sqrt(1000))*1.96)) Lam069 <- c(mean(Scen069sLam), mean(Scen069sLam)+((sd(Scen069sLam)/sqrt(1000))*1.96), mean(Scen069sLam)-((sd(Scen069sLam)/sqrt(1000))*1.96)) Lam070 <- c(mean(Scen070sLam), mean(Scen070sLam)+((sd(Scen070sLam)/sqrt(1000))*1.96), mean(Scen070sLam)-((sd(Scen070sLam)/sqrt(1000))*1.96)) Lam071 <- c(mean(Scen071sLam), mean(Scen071sLam)+((sd(Scen071sLam)/sqrt(1000))*1.96), mean(Scen071sLam)-((sd(Scen071sLam)/sqrt(1000))*1.96)) Lam072 <- c(mean(Scen072sLam), mean(Scen072sLam)+((sd(Scen072sLam)/sqrt(1000))*1.96), mean(Scen072sLam)-((sd(Scen072sLam)/sqrt(1000))*1.96)) Lams <- matrix(c(Lam001, Lam002, Lam003, Lam004, Lam005, Lam006, Lam007, Lam008, Lam009, Lam010, Lam011, Lam012, Lam013, Lam014, Lam015, Lam016, Lam017, Lam018, Lam019, Lam020, Lam021, Lam022, Lam023, Lam024, Lam025, Lam026, Lam027, Lam028, Lam029, Lam030, Lam031, Lam032, Lam033, Lam034, Lam035, Lam036, Lam037, Lam038, Lam039, Lam040, Lam041, Lam042, Lam043, Lam044, Lam045, Lam046, Lam047, Lam048, Lam055, Lam056, Lam057, Lam058, Lam059, Lam060, Lam061, Lam062, Lam063, Lam064, Lam065, Lam066, Lam067, Lam068, Lam069, Lam070, Lam071, Lam072), byrow=TRUE, nrow=66) write.csv(Lams, "Lams.csv", row.names=FALSE, col.names=FALSE)