#### Appendix S2: R code for constructing and analyzing matrices with nonbreeders #### ############################################# ### ### ### Model of nonbreeder dynamics ### ### ### ############################################# # Parameters -------------------------------------------------------------- # census = "pre"/"post" - type of census to be used in building and analyzing matrix # am = youngest possible age at maturity # no.pre = number of prebreeder classes # no.b = number of breeder classes # no.n = number of nonbreeder classes. Usually no.n=no.b ### Survival probabilities: # o.surv = first-year survival if census="post." If pre-breeding census is used first-year survival is included in fec. # i.surv = survival probability of immatures (vector, or single number if all are equal) # p.surv = survival probability of prebreeders (vector, or single number if all are equal) # b.surv = survival probability of breeders (vector, or single number if all are equal) # n.surv = survival probability of nonbreeders (vector, or single number if all are equal) ## Note that if survival is added as a single number the terminal age class is automatically pooled. ### Survival subsets: # i.surv.sub = subsetting of immature survival (e.g. c(1,2,2) means that immature survival has one value for ages 1, one for ages 2-3, and one for ages 4-5) # p.surv.sub = subsetting of prebreeder survival # b.surv.sub = subsetting of breeder survival # n.surv.sub = subsetting of nonbreeder survival # Note that subsets only need to be specified if parameter input is a vector (not if it is a single value) ### Probabilities of becoming a breeder: # pb.am = probability of individuals entering breeding population at age am # pb.pre = probability of prebreeders entering breeding population at next time step (vector, or single number if no.pre=1 or all are equal) # pb.b = probability of breeders being breeders again next time step (vector, or single number if no.b=1 or all are equal) # pb.n = probability of nonbreeders entering breeding population at next time step (vector, or single number if no.n=1 or all are equal) # old = the probability of an individual entering the "old" class at or after age no.pre. For prebreeder probability of breeding will then be (1-old)*pb.pre, and prob. of staying in the prebreeder class is (1-old)(1-pb.pre). For breeders prob of breeding is (1-old)*pb.b # old = 0 means that it goes back to the "normal" model (without senescence) ### Breeding probability subsets: # pb.pre.sub = subsetting of prebreeding transition probabilities (see above) # pb.b.sub = subsetting of breeding transition probabilities # pb.n.sub = subsetting of nonbreeding transition probabilities # Note that in the "old" (senescence) system only the "old" class is programmed as nonbreeders. Earlier ones are prebreeders. # Also note, in this code lambda is calculated based on the immature, young nonbreeder, and breeder population, but not the senescent nonbreeder class (but the whole matrix is included in model runs). # This is not strictly necessary, as results will be the same if the whole matrix is used. # However, including the nonbreeders gives a reducible matrix which causes dependence on initial values for some measures (like absolute population size). # We have therefore set up the code in this way to allow for easy extensions. # pre.pool = TRUE or FALSE. If FALSE, individiuals from last PB class that do not breed enter NB class (instead of staying in PB) # fec = vector of offspring production (including survival to age 1 if prebreeding census) or single number if no.b=1 or all the same # fec.var = variance in offspring production (vector or single number corresponding to fec) # fec.sub = subsetting of fecundity if fec is a vector # LESLIE - TRUE/FALSE, whether or not a Leslie matrix should be constructed and analyzed. Can only be true if age (or if no.n=no.b=1). # classif - "time"/"age", whether breeding probability depends on time spent in class, or age # time.ages - number of adult age classes in Leslie matrix for "recent breeding history" model. If set to zero the number of adult age classes in the Leslie matrix will be set equal to the sum of nonbreeder and breeder classes. # ignoN.model - TRUE/FALSE, testing effects of not including nonbreeders in a model in which only breeders are seen, but individuals are identified # Note that if ignoN.model=TRUE max.age will automatically be set equal to b.last (so all individuals die afther that age). # sens - TRUE/FALSE, Should sensitivity analysis be run? # freq.dep = "none"/"surv"/"fec", Is there frequency dependence? If "surv" or "fec" nonbreeders (and prebreeders if present) have a negative effect on breeder survival or fecundity, respectively. # Realized bsurv = bsurv/(1+NB/B) or realized fec = fec/(1+NB/B) # If freq.dep is not none, use b.surv.base or fec.base, instead of b.surv and fec. # startvec - initial values to use if there is frequency dependence # times - number of time steps to let projection run before testing for equilibrium Nonbreeders1 <- function(census=c("pre", "post"), am, no.pre, no.b, no.n, o.surv=NA, i.surv, i.surv.sub=NULL, p.surv, p.surv.sub=NULL, b.surv=NA, b.surv.base=NA, b.surv.sub=NULL, n.surv, n.surv.sub=NULL, fec=NA, fec.base=NA, fec.sub=NULL, fec.var, pb.am, pb.pre, pb.pre.sub=NULL, pb.b, pb.b.sub=NULL, pb.n, pb.n.sub=NULL, old=0, pre.pool=TRUE, Leslie=FALSE, classif=c("time", "age"), time.ages=0, ignoN.model=FALSE, sens=FALSE, freq.dep=c("none", "surv", "fec"), startvec, times){ if (missing(no.n)) no.n <- no.b if (no.pre==1 & match.arg(classif)=="age") stop("classif must be time when no.pre=1") if (no.n==1 & no.b==1 & match.arg(classif)=="age") stop("classif must be time when no.b=1 and no.n=1") if (match.arg(classif)=="age" & no.pre>no.b) stop("no.pre cannot be greater than no.b when age") if (match.arg(classif)=="age" &! no.n==no.b &! old>0) stop("no.n must be equal to no.b when age (unless old)") if (match.arg(classif)=="age" & pre.pool==TRUE & old==0 &! no.pre==no.b) stop("Cannot pool prebreeders when model is age-classified, unless no.pre=no.b, or old") if (ignoN.model==TRUE & match.arg(census)=="pre") stop("ignoN.model assumes post-breeding census") if (old>0 & Leslie==TRUE) stop("Cannot use Leslie when old") if (old>0 & no.pre==0) stop("Must have prebreeders when old") if (length(i.surv)>1 & is.null(i.surv.sub)) stop("Missing i.surv.sub") if (length(p.surv)>1 & is.null(p.surv.sub)) stop("Missing p.surv.sub") if (length(b.surv)>1 & is.null(b.surv.sub)) stop("Missing b.surv.sub") if (length(n.surv)>1 & is.null(n.surv.sub)) stop("Missing n.surv.sub") if (length(pb.pre)>1 & is.null(pb.pre.sub)) stop("Missing pb.pre.sub") if (length(pb.n)>1 & is.null(pb.n.sub)) stop("Missing pb.n.sub") if (length(pb.b)>1 & is.null(pb.b.sub)) stop("Missing pb.b.sub") if (length(fec)>1 & is.null(fec.sub)) stop("Missing fec.sub") if (no.b>0 & match.arg(freq.dep)=="none" & (is.null(b.surv) | is.null(pb.b) | is.na(b.surv[1]) | is.na(pb.b[1]))) stop("Missing breeder parameter") if (no.n>0 & match.arg(freq.dep)=="none" & (is.null(n.surv) | is.null(pb.n) | is.na(n.surv[1]) | is.na(pb.n[1]))) stop("Missing nonbreeder parameter") if (no.pre>0 & match.arg(freq.dep)=="none" & (is.null(p.surv) | is.null(pb.pre) | is.na(p.surv[1]) | is.na(pb.pre[1]))) stop("Missing prebreeder parameter") if (match.arg(freq.dep)=="surv" & is.na(b.surv.base[1])) stop("Missing b.surv.base") if (match.arg(freq.dep)=="fec" & is.na(fec.base)) stop("Missing fec.base") if (match.arg(freq.dep)!="none" & ignoN.model==TRUE) stop("freq.dep must be none when ignoN.model is true") if (match.arg(freq.dep)!="none" & sens==TRUE) stop("freq.dep must be none when sens is true") if (match.arg(freq.dep)=="fec" & is.na(b.surv.base[1])==FALSE) stop("b.surv.base should not be used when freq.dep=fec") if (match.arg(freq.dep)=="surv" & is.na(fec.base)==FALSE) stop("fec.base should not be used when freq.dep=surv") if (match.arg(census)=="post" & is.na(o.surv)) stop("Need o.surv if post-breeding census") if ((is.null(i.surv.sub)==FALSE & sum(i.surv.sub)!=(am-1)) | (is.null(p.surv.sub)==FALSE & sum(p.surv.sub)!=no.pre) | (is.null(b.surv.sub)==FALSE & sum(b.surv.sub)!=no.b) | (is.null(n.surv.sub)==FALSE & sum(n.surv.sub)!=no.n)) stop("Number of parameters in subclass doesn't match number of classes") if ((is.null(pb.pre.sub)==FALSE & sum(pb.pre.sub)!=no.pre) | (is.null(pb.b.sub)==FALSE & sum(pb.b.sub)!=no.b) | (is.null(pb.n.sub)==FALSE & sum(pb.n.sub)!=no.n)) stop("Number of parameters in subclass doesn't match number of classes") require(popbio) if (match.arg(freq.dep)=="fec") {fec <- fec.base} if (match.arg(freq.dep)=="surv") {b.surv <- b.surv.base} # Subsets if (is.null(i.surv.sub)) {i.surv.sub <- am-1} if (match.arg(census)=="post") {if (am==1) {i.surv.sub <- 1} else {i.surv.sub <- c(1, i.surv.sub)}} i.s.start <- c(1,cumsum(i.surv.sub[-length(i.surv.sub)])+1) i.s.end <- cumsum(i.surv.sub) if (is.null(p.surv.sub)) p.surv.sub <- no.pre p.s.start <- c(1, cumsum(p.surv.sub[-length(p.surv.sub)])+1) p.s.end <- cumsum(p.surv.sub) if (is.null(b.surv.sub)) b.surv.sub <- no.b b.s.start <- c(1, cumsum(b.surv.sub[-length(b.surv.sub)])+1) b.s.end <- cumsum(b.surv.sub) if (is.null(n.surv.sub)) n.surv.sub <- no.n n.s.start <- c(1, cumsum(n.surv.sub[-length(n.surv.sub)])+1) n.s.end <- cumsum(n.surv.sub) if (is.null(pb.pre.sub)) pb.pre.sub <- no.pre pb.pre.start <- c(1, cumsum(pb.pre.sub[-length(pb.pre.sub)])+1) pb.pre.end <- cumsum(pb.pre.sub) if (is.null(pb.b.sub)) pb.b.sub <- no.b pb.b.start <- c(1, cumsum(pb.b.sub[-length(pb.b.sub)])+1) pb.b.end <- cumsum(pb.b.sub) if (is.null(pb.n.sub)) pb.n.sub <- no.n pb.n.start <- c(1, cumsum(pb.n.sub[-length(pb.n.sub)])+1) pb.n.end <- cumsum(pb.n.sub) if (is.null(fec.sub)) fec.sub <- no.b fec.start <- c(1, cumsum(fec.sub[-length(fec.sub)])+1) fec.end <- cumsum(fec.sub) # Parameter vectors if (match.arg(census)=="pre") i.surv.v <- rep(i.surv, i.surv.sub) if (match.arg(census)=="post") {if (am>1) {i.times <- c(i.surv.sub[1]-1, i.surv.sub[-1]) i.surv.v <- c(o.surv, rep(i.surv, i.times[i.times>0]))} else{i.surv.v <- o.surv}} p.surv.v <- rep(p.surv, p.surv.sub) if (ignoN.model==FALSE){ b.surv.v <- rep(b.surv, b.surv.sub) n.surv.v <- rep(n.surv, n.surv.sub)} else { b.surv.v <- c(rep(b.surv, c(b.surv.sub[-length(b.surv.sub)], b.surv.sub[length(b.surv.sub)]-1)), 0) n.surv.v <- c(rep(n.surv, c(n.surv.sub[-length(n.surv.sub)], n.surv.sub[length(n.surv.sub)]-1)), 0) b.surv.sub <- c(b.surv.sub[-length(b.surv.sub)], b.surv.sub[length(b.surv.sub)]-1, 1) n.surv.sub <- c(n.surv.sub[-length(n.surv.sub)], n.surv.sub[length(n.surv.sub)]-1, 1) b.s.start <- c(1, cumsum(b.surv.sub[-length(b.surv.sub)])+1) b.s.end <- cumsum(b.surv.sub) n.s.start <- c(1, cumsum(n.surv.sub[-length(n.surv.sub)])+1) n.s.end <- cumsum(n.surv.sub) if (b.surv[length(b.surv)]!=0) b.surv <- c(b.surv, 0) if (n.surv[length(n.surv)]!=0) n.surv <- c(n.surv, 0) max.age <- no.b+am # This is max.age counted in number of classes, not in years (i.e. 0-year-olds are counted as 1) } fec.v <- rep(fec, fec.sub) fec.var.v <- rep(fec.var, fec.sub) pb.pre.v <- rep(pb.pre, pb.pre.sub) pb.b.v <- rep(pb.b, pb.b.sub) pb.n.v <- rep(pb.n, pb.n.sub) ### PRE-BREEDING CENSUS ------------------------------------------------------------------------------------------------------------------------------------- if (match.arg(census)=="pre"){ # Define matrix positions for easy programming i.last <- am-1 p1 <- am p.last <- am+no.pre-1 b1 <- am+no.pre b.last <- am+no.pre+no.b-1 nb1 <- am+no.pre+no.b nb.last <- am+no.pre+no.b+no.n-1 # List of parameters (subsets) for survival and transitions par.sub.names <- c( sapply(1:length(i.surv.sub), function(x) paste("i.surv",x, sep="")), sapply(1:length(p.surv.sub), function(x) paste("p.surv",x, sep="")), sapply(1:length(b.surv.sub), function(x) paste("b.surv",x, sep="")), sapply(1:length(n.surv.sub), function(x) paste("n.surv",x, sep="")), "pb.am", sapply(1:length(pb.pre.sub), function(x) paste("pb.pre",x, sep="")), sapply(1:length(pb.b.sub), function(x) paste("pb.b",x, sep="")), sapply(1:length(pb.n.sub), function(x) paste("pb.n",x, sep="")), "old") if (am==1) par.sub.names <- par.sub.names[-which(par.sub.names=="pb.am")] # Parameter values, separated by survival and fecundity if (am>1) {par.vals <- c(i.surv, p.surv, b.surv, n.surv, pb.am, pb.pre, pb.b, pb.n, old) } else {par.vals <- c(i.surv, p.surv, b.surv, n.surv, pb.pre, pb.b, pb.n, old)} if (am>1) {par.vals.f <- fec} else {par.vals.f <- c(fec, pb.am)} # Create matrix while keeping track of all parameters ----------------------------------------------------------- # Set up lists for storing indices Ind.list.surv <- lapply(1:length(par.sub.names), function(x) matrix(NA, nrow=nb.last*nb.last, ncol=3, dimnames=list(NULL,c("x", "y", "Inv")))) names(Ind.list.surv) <- par.sub.names Ind.list.fec <- lapply(1:(length(fec.sub)+1), function(x) matrix(NA, nrow=nb.last*nb.last, ncol=3, dimnames=list(NULL,c("x", "y", "Inv")))) names(Ind.list.fec) <- c(sapply(1:length(fec.sub), function(x) paste("fec", x, sep="")), "pb.am") if (am>1) Ind.list.fec <- Ind.list.fec[-which(names(Ind.list.fec)=="pb.am")] # Set up temporary matrix for holding indices temp.indmat <- matrix(NA, nrow=nb.last, ncol=3) # Survival of immatures if (am>2) {temp.indmat[1:(i.last-1),] <- cbind(2:i.last, 1:(i.last-1), 0) for (i in 1:length(i.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", i, sep=""))]][1:i.surv.sub[i],] <- temp.indmat[i.s.start[i]:i.s.end[i],]} } # Transitions from immature to nonbreeder if no.pre=0 if (no.pre==0 & am>1) { Ind.list.surv$pb.am[1,] <- c(nb1, i.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(nb1, i.last, 0)} # Survival and probability of prebreeders staying in PB (or immatures entering). Also PB entering NB if pre.pool=FALSE if (no.pre >0) { if (old==0){ if (pre.pool==FALSE) { if (match.arg(classif)=="time") { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb1, p.last, 1) } else { if (no.n > no.pre) { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1+no.pre, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb1+no.pre, p.last, 1) } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb.last, p.last, 1) }} } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(p.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(p.last, p.last, 1) }} else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1, p.last, 0) Ind.list.surv$old[1,] <- c(nb1, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][2,] <- c(p.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(p.last, p.last, 1) Ind.list.surv$old[2,] <- c(p.last, p.last, 1) } if (am>1) { Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(p1, i.last, 0) Ind.list.surv$pb.am[2,] <- c(p1, i.last, 1)} if (no.pre>1) { temp.indmat[,] <- NA temp.indmat[1:(no.pre-1),] <- cbind((p1+1):(p.last), p1:(p.last-1), 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } temp.indmat[1:(no.pre-1),3] <- 1 for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }} # If age matters (both B and NB age-classified): if (no.b>1 & match.arg(classif)=="age"){ # Transitions from PB (or juv) to B if (am==1 & no.pre>0){ if (no.pre==no.b) { temp.indmat[,] <- NA temp.indmat[1:(no.pre),] <- cbind(c((b1+1):(b.last), b.last), p1:(p.last), 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } Ind.list.surv$old[min(which(is.na(Ind.list.surv$old))),] <- c(b.last, p.last, 1) } else{ temp.indmat[,] <- NA temp.indmat[1:(no.pre),] <- cbind((b1+1):(am+no.pre^2), p1:(p.last), 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }} else{ if (no.pre==no.b) { temp.indmat[,] <- NA temp.indmat[1:no.b,] <- cbind(c((b1+1):b.last, b.last), p1:p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) if (am>1) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am))),] <- c(b1, i.last, 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } if (old>0) {Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(b.last, p.last, 1) } } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) if (am>1) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- c(b1, i.last, 0) if (no.pre>0){ temp.indmat[,] <- NA temp.indmat[1:no.pre,] <- cbind(c((b1+1):(am+no.pre*2)), p1:p.last, 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }}} # Transitions from NB to B if (old==0) { temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((b1+1):b.last, b.last), nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][1:n.surv.sub[i],] <- temp.indmat[n.s.start[i]:n.s.end[i],] } for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][1:pb.n.sub[i],] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } } # Transitions from B to NB if (old==0) { temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((nb1+1):nb.last, nb.last), b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][1:b.surv.sub[i],] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][1:pb.b.sub[i],] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } else{ Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][,1]))),] <- c(nb1, b.last, 0) Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(nb1, b.last, 0) # Transitions from B to PB in old system temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((p1+1):p.last, p.last), b1:b.last, 0) for (i in 1:(length(b.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:(length(pb.b.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(p.last, b.last, 1) } } # If time matters for B or NB (or they are both 1) if (match.arg(classif)=="time"){ # Transitions from PB (or immature) to B if (am==1 & no.pre!=0){ temp.indmat[,] <- NA temp.indmat[1:no.pre,] <- cbind(b1, p1:p.last, 0) for (i in 1:(length(p.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:(length(pb.pre.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } } if (am>1) { Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- c(b1, i.last, 0) if (no.pre>0){ temp.indmat[,] <- NA temp.indmat[1:no.pre,] <- cbind(b1, p1:p.last, 0) for (i in 1:(length(p.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:(length(pb.pre.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }} # Transitions from NB to B temp.indmat[,] <- NA temp.indmat[1:(no.n),] <- cbind(b1, nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1])))+n.surv.sub[i]-1),] <- temp.indmat[n.s.start[i]:n.s.end[i],] } for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1])))+pb.n.sub[i]-1),] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } # Transitions from B to NB temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(nb1, b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } # Survival and probability of breeders staying in B if (no.b>1){ temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((b1+1):b.last, b.last), b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } else {Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][,1]))),] <- c(b.last, b.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", length(pb.b.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", length(pb.b.sub), sep=""))]][,1]))),] <- c(b.last, b.last, 0) } if (old>0) {Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(b.last, b.last, 1) } # Survival and probability of nonbreeders staying in NB if (no.n>1){ temp.indmat[,] <- NA temp.indmat[1:(no.n),] <- cbind(c((nb1+1):nb.last, nb.last), nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1])))+n.surv.sub[i]-1),] <- temp.indmat[n.s.start[i]:n.s.end[i],] } temp.indmat[1:no.n,3] <- 1 temp.indmat[no.n,] <- NA for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1])))+pb.n.sub[i]-1),] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } if (old==0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))),] <- c(nb.last, nb.last, 1) } else{ Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", length(n.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", length(n.surv.sub), sep=""))]][,1]))),] <- c(nb.last, nb.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][,1]))),] <- c(nb.last, nb.last, 1) } # Fecundity if (am>1) { temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(1, b1:b.last, 0) for (i in 1:length(fec.sub)){Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][1:fec.sub[i],] <- temp.indmat[fec.start[i]:fec.end[i],]} } else{ temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(b1, b1:b.last, 0) for (i in 1:length(fec.sub)){Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][1:fec.sub[i],] <- temp.indmat[fec.start[i]:fec.end[i],]} Ind.list.fec$pb.am[min(which(is.na(Ind.list.fec$pb.am[,1]))):(min(which(is.na(Ind.list.fec$pb.am[,1])))+no.b-1),] <- temp.indmat[1:no.b,] if (no.pre==0){ temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(nb1, b1:b.last, 0) for (i in 1:length(fec.sub)){Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][min(which(is.na(Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[fec.start[i]:fec.end[i],]} temp.indmat[1:no.b,3] <- 1 Ind.list.fec$pb.am[min(which(is.na(Ind.list.fec$pb.am[,1]))):(min(which(is.na(Ind.list.fec$pb.am[,1])))+no.b-1),] <- temp.indmat[1:no.b,] } else{ temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(p1, b1:b.last, 0) for (i in 1:length(fec.sub)){Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][min(which(is.na(Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.fec[[which(names(Ind.list.fec)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[fec.start[i]:fec.end[i],]} temp.indmat[1:no.b,3] <- 1 Ind.list.fec$pb.am[min(which(is.na(Ind.list.fec$pb.am[,1]))):(min(which(is.na(Ind.list.fec$pb.am[,1])))+no.b-1),] <- temp.indmat[1:no.b,] }} # Create separate matrices for each of the parameters par.mats.s <- lapply(1:length(par.sub.names), function(x) matrix(NA, nrow=nb.last, ncol=nb.last)) names(par.mats.s) <- paste(par.sub.names, ".mat", sep="") for (i in 1:length(par.sub.names)) { par.mats.s[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,2])] <- par.vals[i] par.mats.s[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,2])] <- (1-par.vals[i])} par.mats.f <- lapply(1:(length(fec.sub)+1), function(x) matrix(NA, nrow=nb.last, ncol=nb.last)) names(par.mats.f) <- paste(names(Ind.list.fec), ".mat", sep="") for (i in 1:(length(names(Ind.list.fec)))) { par.mats.f[[i]][cbind(Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==0,1], Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==0,2])] <- par.vals.f[i] par.mats.f[[i]][cbind(Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==1,1], Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==1,2])] <- 1-par.vals.f[i]} # Replace NA values with 1, or 0 if they are all NA (this makes it possible to multiply them together properly) unl.Ind.s <- lapply(1:length(par.mats.s), function(x) unlist(par.mats.s[x])) surv.NAs <- lapply(1:(nb.last*nb.last), function(x) sapply(1:length(par.mats.s), function(y) is.na(unl.Ind.s[[y]][x]))) surv.all.NAs <- lapply(1:length(surv.NAs), function(x) all(surv.NAs[[x]]==TRUE)) for (i in 1:length(par.sub.names)) {par.mats.s[[i]][which(surv.all.NAs==TRUE)] <- 0 par.mats.s[[i]][is.na(par.mats.s[[i]])] <- 1 } unl.Ind.f <- lapply(1:length(par.mats.f), function(x) unlist(par.mats.f[x])) fec.NAs <- lapply(1:(nb.last*nb.last), function(x) sapply(1:length(par.mats.f), function(y) is.na(unl.Ind.f[[y]][x]))) fec.all.NAs <- lapply(1:length(fec.NAs), function(x) all(fec.NAs[[x]]==TRUE)) for (i in 1:(length(fec.sub)+1)) {par.mats.f[[i]][which(fec.all.NAs==TRUE)] <- 0 par.mats.f[[i]][is.na(par.mats.f[[i]])] <- 1 } # Multiply together all matrices to create a single survival (and transition) matrix and a fecundity matrix for (i in 2:length(par.mats.s)) {if (i==2) {NMs <- par.mats.s[[1]]} NMs <- NMs*par.mats.s[[i]]} for (i in 2:length(par.mats.f)) {if (i==2) {NMf <- par.mats.f[[1]]} NMf <- NMf*par.mats.f[[i]]} NMsurv <- NMs NMfec <- NMf # Add together survival and fecundity matrix into total population projection matrix NM <- NMsurv+NMfec # Row names Classes <- which(c(am-1, no.pre, no.b, no.n)>0) name.let <- list("i%d", "p%d", "b%d", "nb%d") name.num <- list(1:(am-1), 1:no.pre, 1:no.b, 1:no.n) Nclasses <- c(unlist(mapply(sprintf, name.let[Classes], name.num[Classes]))) rownames(NM) <- Nclasses # Frequency dependent case ------------------------------------------------------------------------------------------------ if (match.arg(freq.dep)!="none"){ # Create matrix to record population projection popsize <- matrix(nrow=nb.last, ncol=times+1) # Insert initial values popsize[,1] <- startvec # Define function for updating frequency dependent matrix fdMat <- function(tt, popsize=popsize){ if (old>0 | no.pre>0) {noNB <- sum(popsize[c(p1:p.last, nb1:nb.last),tt-1])} else { noNB <- sum(popsize[c(nb1:nb.last),tt-1]) } noB <- sum(popsize[c(b1:b.last), tt-1]) if (freq.dep=="fec"){ # Calculate breeder survival adjusted by nonbreeder density fec.fd <- fec/(1+noNB/noB) # Insert fec.fd into matrix for (i in 1:(length(fec.sub))) {par.mats.f[[i]][cbind(Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==0,1], Ind.list.fec[[i]][is.na(Ind.list.fec[[i]][,1])==FALSE & Ind.list.fec[[i]][,3]==0,2])] <- fec.fd[i]} for (i in 2:length(par.mats.f)) {if (i==2) {NMf <- par.mats.f[[1]]} NMf <- NMf*par.mats.f[[i]]} NMfd <- NMsurv+NMf} if (freq.dep=="surv"){ # Calculate breeder survival adjusted by nonbreeder density b.surv.fd <- b.surv/(1+noNB/noB) # Insert b.surv.fd into matrix for (i in grep("b.surv", par.sub.names)) {par.mats.s[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,2])] <- b.surv.fd[i-min(grep("b.surv", par.sub.names))+1] par.mats.s[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,2])] <- (1-b.surv.fd[i-min(grep("b.surv", par.sub.names))+1])} for (i in 2:length(par.mats.s)) {if (i==2) {NMs <- par.mats.s[[1]]} NMs <- NMs*par.mats.s[[i]]} NMfd <- NMs+NMfec} return(NMfd)} # Project population over "times" time steps for (tt in 2:(times+1)){ NMfd <- fdMat(tt, popsize) # Project population at next time step newpopsize <- NMfd %*% popsize[,tt-1] popsize[,tt] <- newpopsize } NMfd.eq <- fdMat(times+1, popsize) rownames(NMfd.eq) <- rownames(NM) if (old>0){populationsize <- popsize[1:b.last,] NMfd.equil <- NMfd.eq[-nrow(NMfd.eq), -nrow(NMfd.eq)]} else { populationsize <- popsize NMfd.equil <- NMfd.eq} phat <- populationsize[,times]/sum(populationsize[,times]) phat.test <- (NMfd.equil %*% phat) / norm(NMfd.equil %*% phat) test <- phat.test - phat if (any(abs(test)>0.000000001)) stop("not at equilibrium") # Find deterministic growth rate eig_info <- eigen(NMfd.equil) ev <- which.max(Re(eig_info$values)) eig_val <- Re(eig_info$values[ev]) Lambda.test <- (colSums(populationsize)[times+1]/colSums(populationsize)[times-9])^(1/10) if (abs(Lambda.test-eig_val)>0.000000001) stop("not at equilibrium") } else { # Find stable stage structure and reproductive values ------------------------------- # For "old" system use matrix without senescent class. if (old>0) {NM <- NM[-nrow(NM), -nrow(NM)] NMsurv <- NMsurv[-nb1, -nb1] NMfec <- NMfec[-nb1, -nb1]} # Define function Eig.calc <- function(m.matrix){ # Stable stage structure eig_info <- eigen(m.matrix) ev <- which.max(Re(eig_info$values)) eig_val <- Re(eig_info$values[ev]) u_vals <- Re((eig_info$vectors[,ev])/sum(eig_info$vectors[,ev])) # Reproductive values eig_info2 <- eigen(t(m.matrix)) ev2 <- which(round(eig_val, digits=5)==round(eig_info2$values, digits=5)) v_vals <- Re(eig_info2$vectors[,ev2])/sum(u_vals*Re(eig_info2$vectors[,ev2])) # Test eigenvalues uvector.test <- m.matrix%*%u_vals - eig_val*u_vals if (any(uvector.test>10^-4)) stop("problem with u_vals") vvector.test <- t(v_vals)%*%m.matrix - eig_val*t(v_vals) if (any(vvector.test>10^-4)) stop("problem with v_vals") # if (any(vvector.test>10^-4)) v_vals <- NA # With the old system there will very occasionally be a problem with the v-values (only shows up when running many iterations). This line lets it continue running anyway. Use with care. return(list(eig_val, u_vals, v_vals))} # Run function on NM matrix NM.info <- Eig.calc(NM) eig_val <- NM.info[[1]] u_vals <- NM.info[[2]] v_vals <- NM.info[[3]] if (old>0) Nclasses <- Nclasses[-length(Nclasses)] info.mat <- data.frame(Nclasses, u_vals, v_vals) # Caculate sensitivities ------------------------------------------------ if (sens==TRUE){ # Matrix element level sensitivities tot.NMsens <- sensitivity(NM, zero=TRUE) # Create matrix to store lower level sensitivities in sens.mat <- matrix(nrow=length(names(Ind.list.surv))+length(names(Ind.list.fec)), ncol=1, dimnames=list(c(par.sub.names, names(Ind.list.fec)))) for (i in 1:length(names(Ind.list.surv))) { if (old>0){Ind.list.surv[[i]] <- Ind.list.surv[[i]][-c(which(apply(Ind.list.surv[[i]], 1, function(x) any(x==nb.last))==TRUE), nb.last*nb.last),] } Ind.no <- length(which(is.na(Ind.list.surv[[i]][,1])==FALSE)) if (Ind.no>0) { temp.sens.ind <- rep(NA, Ind.no) for (j in 1:Ind.no) { if (Ind.list.surv[[i]][j, 3]==0){ temp.sens.ind[j] <- tot.NMsens[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]*NMsurv[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]/par.vals[i] } else { temp.sens.ind[j] <- -1*tot.NMsens[Ind.list.surv[[i ]][j, 1], Ind.list.surv[[i]][j, 2]]*NMsurv[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]/(1-par.vals[i]) } } sens.mat[i] <- sum(temp.sens.ind) }} for (i in 1:length(names(Ind.list.fec))) { Ind.no <- length(which(is.na(Ind.list.fec[[i]][,1])==FALSE)) if (Ind.no>0) { temp.sens.ind <- rep(NA, Ind.no) for (j in 1:Ind.no) { if (Ind.list.fec[[i]][j, 3]==0){ temp.sens.ind[j] <- tot.NMsens[Ind.list.fec[[i]][j, 1], Ind.list.fec[[i]][j, 2]]*NMfec[Ind.list.fec[[i]][j, 1], Ind.list.fec[[i]][j, 2]]/par.vals.f[i] } else { temp.sens.ind[j] <- -1*tot.NMsens[Ind.list.fec[[i]][j, 1], Ind.list.fec[[i]][j, 2]]*NMfec[Ind.list.fec[[i]][j, 1], Ind.list.fec[[i]][j, 2]]/(1-par.vals.f[i]) } } sens.mat[length(names(Ind.list.surv))+i] <- sum(temp.sens.ind) }} } else {tot.NMsens <- NA sens.mat <- NA} # Calculate demographic variance ------------------------------------------ # Remove values associated with old individuals if (old>0){u_vals_dv <- u_vals[-nb1] v_vals_dv <- v_vals[-nb1] last.cov <- b.last} else {u_vals_dv <- u_vals v_vals_dv <- v_vals last.cov <- nb.last} # Variance in survival term: surv.var <- apply(NMsurv, 2, function(x) sum(x*(1-x)* v_vals_dv^2)) surv.cov <- fec.var.term <- fec.cov.term <- vector("numeric", length=length(surv.var)) # Covariance in survival/transition term: first.cov <- max(am-1,1) # There are no covarances for immature survival terms for (i in first.cov:last.cov){if (length(which(NMsurv[,i]>0))>1){ rows <- combn(which(NMsurv[,i]>0),2) surv.cov[i] <- sum(apply(rows, 2, function(x) -prod(NMsurv[x,i])*prod(v_vals_dv[x])))} else surv.cov[i] <- 0} if (am>1){ # Variance in reproduction term if am>1: fec.var.term[b1:b.last] <- fec.var.v*v_vals_dv[1]^2 # Demographic variance if am>1: DemVar <- sum(u_vals_dv*(surv.var+2*surv.cov+fec.var.term)) } else{ # Variance in reproduction terms if am=1: if(no.b>1) fec.probs <- t(apply(NMfec[, b1:b.last], 1, function(x) x/fec.v)) else fec.probs <- NMfec[, b1]/fec.v # Check if probabilities sum to one if (no.b>1){ if (all(colSums(fec.probs)>0.9999999999 & colSums(fec.probs) < 1.0000000001)==FALSE) stop("Problem with colSums of fec.probs") fec.rows <- which(fec.probs[,1]>0) fec.var.term[b1:b.last] <- sapply(1:no.b, function(x) sum((rep(fec.v[x],2)*fec.probs[fec.rows, x]*(1-fec.probs[fec.rows,x]) + rep(fec.var.v[x],2)*fec.probs[fec.rows,x]^2)*v_vals_dv[fec.rows]^2)) # Covariance in reproduction terms if am=1 fec.cov.term[b1:b.last] <- sapply(1:ncol(fec.probs), function(x) prod(fec.probs[fec.rows,x])*(fec.var.v[x]-fec.v[x])*prod(v_vals_dv[fec.rows])) } else { if (sum(fec.probs)<0.99999 | sum(fec.probs)>1.00001) stop("Problem with sum of fec.probs") fec.rows <- which(fec.probs>0) fec.var.term[b1] <- sum((fec.v*fec.probs[fec.rows]*(1-fec.probs[fec.rows]) + fec.var.v*fec.probs[fec.rows]^2)*v_vals_dv[fec.rows]^2) fec.cov.term[b1] <- prod(fec.probs[fec.rows])*(fec.var.v-fec.v)*prod(v_vals_dv[fec.rows]) } DemVar <- sum(u_vals_dv*(surv.var+2*surv.cov+fec.var.term+2*fec.cov.term)) } # Construct and analyze Leslie matrix ------------------------------------------ if (Leslie==TRUE){ # Not intended for use with old # In "recent breeding history" system. First create expanded matrix with both age and breeding structure. # Note that some stage combinations are impossible even if they have demographic parameters in the expanded matrix (entering these stages is impossible, so those parameters will never be used). if (match.arg(classif)=="time" & no.pre==0 & (no.b>1 | no.n>1)) { agecl.matrix <- matrix(0, nrow=nrow(NM)+(no.n+no.b)*(no.n+no.b-1), ncol=nrow(NM)+(no.n+no.b)*(no.n+no.b-1)) rownames(agecl.matrix) <- c(1:(am-1), rep(am:(am-1+no.n+no.b), each=(no.n+no.b))) #Rownames indicate ages colnames(agecl.matrix) <- c(row.names(NM)[1:(am-1)], rep(row.names(NM)[am:dim(NM)[1]], times=(no.n+no.b))) #Colnames indicate breeding class # If time.ages=0 number of age classes is set equal to sum of nonbreeder and breeder classes max.age.class <- ifelse(time.ages==0, am+no.n+no.b-1, am+time.ages-1) # Fecundity agecl.matrix[1:(am-1),1:nrow(NM)] <- NM[1:(am-1),] age.blocks <- seq(nrow(NM)+1, nrow(agecl.matrix), no.n+no.b) for (i in age.blocks){ agecl.matrix[1:(am-1),i:(i+(no.n+no.b-1))] <- NM[1:(am-1), am:nrow(NM)]} # Survival agecl.matrix[am:nrow(NM), 1:(am-1)] <- NM[am:(nrow(NM)),1:(am-1)] for (i in age.blocks){ agecl.matrix[i:(i+(no.n+no.b-1)), (i-(no.n+no.b)):(i-1)] <- NM[am:nrow(NM), am:nrow(NM)] } agecl.matrix[age.blocks[length(age.blocks)]:(age.blocks[length(age.blocks)]+(no.n+no.b-1)),age.blocks[length(age.blocks)]:(age.blocks[length(age.blocks)]+(no.n+no.b-1))] <- NM[am:nrow(NM), am:nrow(NM)] # Analyze the expanded matrix agecl.info <- Eig.calc(agecl.matrix) agecl.u_vals <- agecl.info[[2]] # Construct Leslie matrix directly from expanded matrix (using stable age distribution) umage <- t(apply(agecl.matrix, 1, function(x) x*agecl.u_vals)) L.matrixB.t <- matrix(NA, ncol=max.age.class, nrow=max.age.class) for (i in 1:(max.age.class-1)){for (j in 1:(max.age.class-1)){ L.matrixB.t[i,j] <- sum(umage[which(rownames(agecl.matrix)==i), which(rownames(agecl.matrix)==j)]) }} for (i in 1:(max.age.class-1)){L.matrixB.t[i, max.age.class] <- sum(umage[which(rownames(agecl.matrix)==i), which(rownames(agecl.matrix)>=max.age.class)])} for (j in 1:(max.age.class-1)){L.matrixB.t[max.age.class, j] <- sum(umage[which(rownames(agecl.matrix)>=max.age.class), which(rownames(agecl.matrix)==j)])} L.matrixB.t[max.age.class, max.age.class] <- sum(umage[which(rownames(agecl.matrix)>=max.age.class), which(rownames(agecl.matrix)>=max.age.class)]) u.agecl.m <- rep(NA, max.age.class) for (i in 1:(max.age.class-1)){u.agecl.m[i] <- sum(agecl.u_vals[which(rownames(agecl.matrix)==i)])} u.agecl.m[max.age.class] <- sum(agecl.u_vals[which(rownames(agecl.matrix)>=max.age.class)]) L.matrixB <- t(apply(L.matrixB.t, 1, function(x) x/u.agecl.m)) # Calculate fecundity for the Leslie model L.fec <- c(rep(0, am-1), rep(NA, max.age.class-am)) if(max.age.class>am){ for (i in 1:(max.age.class-am)){L.fec[am-1+i] <- sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.b-1)]*fec.v)/sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.n+no.b-1)])}} else { L.fec[am] <- sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.b-1)]*fec.v)/sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.n+no.b-1)])} L.fec[max.age.class] <- sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.b-1))]*fec.v)/sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.n+no.b-1))]) L.fec.var <- c(rep(0, am-1), rep(NA, max.age.class-am)) if(max.age.class>am){ for (i in 1:(max.age.class-am)){L.fec.var[am-1+i] <- sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.b-1)]*fec.var.v+agecl.u_vals[(c(am, age.blocks)[i]+no.b):(c(am, age.blocks)[i]+no.n+no.b-1)]*fec.v^2)/sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.n+no.b-1)])}} else { L.fec.var[am] <- sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.b-1)]*fec.var.v+agecl.u_vals[(c(am, age.blocks)[1]+no.b):(c(am, age.blocks)[1]+no.n+no.b-1)]*fec.v^2)/sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[i]+no.n+no.b-1)]) } L.fec.var[max.age.class] <- sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.b-1))]*fec.var.v+agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)]+no.b, function(x) x:(x+no.n-1))]*fec.v^2)/sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.n+no.b-1))]) # Calculate survival for the Leslie model if (am>1) {j.term <- i.surv.v} else j.term <- NULL L.surv <- c(j.term, rep(NA, max.age.class-am)) if(max.age.class>am){for (i in 1:(max.age.class-am)){L.surv[am-1+i] <- sum((agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.b-1)])*b.surv.v+(agecl.u_vals[(c(am, age.blocks)[i]+no.b):(c(am, age.blocks)[i]+no.n+no.b-1)])*n.surv.v)/sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.n+no.b-1)])} } else { L.surv[am] <- sum((agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[i]+no.b-1)])*b.surv.v+(agecl.u_vals[(c(am, age.blocks)[1]+no.b):(c(am, age.blocks)[1]+no.n+no.b-1)])*n.surv.v)/sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.n+no.b-1)]) } L.surv[max.age.class] <- sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.b-1))]*b.surv.v+(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)]+no.b, function(x) x:(x+no.n-1))])*n.surv.v)/sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.n+no.b-1))]) # Construct Leslie matrix using calculated parameters L.matrix <- matrix(0, ncol=max.age.class, nrow=max.age.class) L.matrix[1,] <- L.fec if (am>1) {L.matrix[cbind(2:max.age.class, 1:(max.age.class-1))] <- L.surv[-max.age.class]} L.matrix[max.age.class, max.age.class] <- L.surv[max.age.class] # Covariance between fecundity and survival nsurvterm <- rep(NA, max.age.class-am+1) bsurvterm <- rep(NA, max.age.class-am+1) if(max.age.class>am){ for (i in 1:(max.age.class-am)) nsurvterm[i] <- sum((n.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[(c(am, age.blocks)[i]+no.b):(c(am, age.blocks)[i]+no.n+no.b-1)])/sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.n+no.b-1)])} else { nsurvterm[1] <- sum((n.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[(c(am, age.blocks)[1]+no.b):(c(am, age.blocks)[1]+no.n+no.b-1)])/sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.n+no.b-1)])} nsurvterm[max.age.class-am+1] <- sum((n.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)]+no.b, function(x) x: (x+no.n-1))])/sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.n+no.b-1))]) if(max.age.class>am){ for (i in 1:(max.age.class-am)) bsurvterm[i] <- sum((b.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.b-1)])/sum(agecl.u_vals[c(am, age.blocks)[i]:(c(am, age.blocks)[i]+no.n+no.b-1)])} else { bsurvterm[1] <- sum((b.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.b-1)])/sum(agecl.u_vals[c(am, age.blocks)[1]:(c(am, age.blocks)[1]+no.n+no.b-1)]) } bsurvterm[max.age.class-am+1] <- sum((b.surv.v-L.surv[am:nrow(L.matrix)])*agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x: (x+no.b-1))])/sum(agecl.u_vals[sapply(c(am, age.blocks)[(max.age.class-am+1):(length(age.blocks)+1)], function(x) x:(x+no.n+no.b-1))]) LScov <- c(rep(0, am-1), (-L.fec[am:nrow(L.matrix)])*nsurvterm + (sum(fec.v*u_vals[b1:b.last]/sum(u_vals[b1:b.last]))-L.fec[am:nrow(L.matrix)])*bsurvterm) } else { # Stable age distribution u.breed <- u_vals[b1:b.last] u.non <- u_vals[nb1:nb.last] if (no.pre>0) {u.pre <- c(u_vals[p1:p.last],rep(0, no.b-no.pre))} else u.pre <- 0 if (am>1) {u.juv <- u_vals[1:(am-1)]} else u.juv <- 0 u.age.ad <- u.breed+u.non+u.pre if (am>1) {u.age <- c(u.juv, u.age.ad)} else u.age <- u.age.ad # sorting classes by age if (am>1) {j.age.index <- lapply(1:i.last, function(z) z)} else {j.age.index <- NULL} if (no.pre>0) {p.age.index <- lapply(am:p.last, function(x) c(x, p.last+x-am+1, b.last+x-am+1))} else {p.age.index <- NULL} if (no.b > no.pre) {a.age.index <- lapply(1:(no.b-no.pre), function(y) c(p.last+no.pre+y, b.last+no.pre+y) )} else {a.age.index <- NULL} age.index <- c(j.age.index, p.age.index, a.age.index) # Construct Leslie matrix directly from NM matrix (using stable age distribution) uNM <- t(apply(NM, 1, function(x) x*u_vals)) L.matrixB.t <- matrix(NA, ncol=length(u.age), nrow=length(u.age)) for (i in 1:length(u.age)){ for (j in 1:length(u.age)){ L.matrixB.t[i,j] <- sum(uNM[age.index[[i]], age.index[[j]]]) } } L.matrixB <- t(apply(L.matrixB.t, 1, function(x) x/u.age)) # Calculate fecundity for the Leslie model L.fec <- c(rep(0, am-1), u.breed*fec.v/u.age.ad) L.fec.var <- c(rep(0, am-1), (u.breed*fec.var.v+(u.pre+u.non)*fec.v^2)/u.age.ad) # Calculate survival for the Leslie model if (no.pre>0) {p.term <- c(p.surv.v,rep(0, no.b-no.pre))*u.pre} else p.term <- 0 if (am>1) {j.term <- i.surv.v} else j.term <- NULL L.surv <- c(j.term, (b.surv.v*u.breed+n.surv.v*u.non+p.term)/u.age.ad) # Construct Leslie matrix using calculated parameters L.matrix <- matrix(0, ncol=am-1+no.b, nrow=am-1+no.b) L.matrix[1,] <- L.fec if (am>1 | no.b>1) {L.matrix[cbind(2:(am-1+no.b), 1:(am-2+no.b))] <- L.surv[-(am-1+no.b)]} if (am==1 & no.b==1) {L.matrix[1,1] <- L.matrix[1,1]+L.surv[1]} else { L.matrix[(am-1+no.b), (am-1+no.b)] <- L.surv[am-1+no.b]} # Covariance between fecundity and survival LScov <- c(rep(0, am-1), (-L.fec[am:nrow(L.matrix)])*(n.surv.v-L.surv[am:nrow(L.matrix)])*u.non/u.age.ad+(sum(fec.v*u_vals[b1:b.last]/sum(u_vals[b1:b.last]))-L.fec[am:nrow(L.matrix)])*(b.surv.v-L.surv[am:nrow(L.matrix)])*u.breed/u.age.ad) } # Check that the two Leslie matrices are the same if (any(abs(L.matrix-L.matrixB)>0.000001)) stop("Problem with L.matrix") # Growth rate, stable stage distribution and reproductive value L.info <- Eig.calc(L.matrix) L.eig_val <- L.info[[1]] L.u_vals <- L.info[[2]] L.v_vals <- L.info[[3]] # Demographic variance L.DemVar <- sum(L.u_vals*(L.surv*(1-L.surv)*c(L.v_vals[-1]^2,L.v_vals[dim(L.matrix)[1]]^2)+L.fec.var*L.v_vals[1]^2)) L.DemVar.cov <- sum(L.u_vals*(L.surv*(1-L.surv)*c(L.v_vals[-1]^2,L.v_vals[dim(L.matrix)[1]]^2)+L.fec.var*L.v_vals[1]^2+2*LScov*L.v_vals[1]*c(L.v_vals[-1],L.v_vals[dim(L.matrix)[1]]))) } else {L.DemVar <- NA L.DemVar.cov <- NA L.eig_val <- NA L.info <- NA L.u_vals <- NA L.v_vals <- NA} ignoN.eig_val <- NA ignoN.u_vals <- NA ignoN.v_vals <- NA} } ### POST-BREEDING CENSUS ------------------------------------------------------------------------------------------------------------------------------------- if (match.arg(census)=="post"){ # List of parameters (subsets) for survival and transitions par.sub.names <- c( sapply(1:length(i.surv.sub), function(x) paste("i.surv",x, sep="")), sapply(1:length(p.surv.sub), function(x) paste("p.surv",x, sep="")), sapply(1:length(b.surv.sub), function(x) paste("b.surv",x, sep="")), sapply(1:length(n.surv.sub), function(x) paste("n.surv",x, sep="")), "pb.am", sapply(1:length(pb.pre.sub), function(x) paste("pb.pre",x, sep="")), sapply(1:length(pb.b.sub), function(x) paste("pb.b",x, sep="")), sapply(1:length(pb.n.sub), function(x) paste("pb.n",x, sep="")), "old", sapply(1:length(fec.sub), function(x) paste("fec",x, sep=""))) fec.v <- c(fec.v, fec.v[length(fec.v)]) fec.var.v <- c(fec.var.v, fec.var.v[length(fec.var.v)]) # Parameter values if (am==1) {par.vals <- c(o.surv, p.surv, b.surv, n.surv, pb.am, pb.pre, pb.b, pb.n, old, fec)} else { par.vals <- c(o.surv, i.surv, p.surv, b.surv, n.surv, pb.am, pb.pre, pb.b, pb.n, old, fec)} # Redefine am, since first class now consists of 0-year-olds (post-breeding census) am <- am+1 # Define matrix positions for easy programming i.last <- am-1 p1 <- am p.last <- am+no.pre-1 b1 <- am+no.pre b.last <- am+no.pre+no.b-1 nb1 <- am+no.pre+no.b nb.last <- am+no.pre+no.b+no.n-1 # Create matrix while keeping track of all parameters ----------------------------------------------------------- # Set up lists for storing indices Ind.list.surv <- lapply(1:length(par.sub.names), function(x) matrix(NA, nrow=nb.last*nb.last, ncol=3, dimnames=list(NULL,c("x", "y", "Inv")))) names(Ind.list.surv) <- par.sub.names # Set up temporary matrix for holding indices temp.indmat <- matrix(NA, nrow=nb.last, ncol=3) # Survival of immatures if (am>2) {temp.indmat[1:(i.last-1),] <- cbind(2:i.last, 1:(i.last-1), 0) for (i in 1:(length(i.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", i, sep=""))]][1:i.surv.sub[i],] <- temp.indmat[i.s.start[i]:i.s.end[i],]} } # Transitions from immature to nonbreeder if no.pre=0 if (no.pre==0) { Ind.list.surv$pb.am[1,] <- c(nb1, i.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(nb1, i.last, 0)} # Survival and probability of prebreeders staying in PB (or immatures entering). Also PB entering NB if pre.pool=FALSE if (no.pre >0) { if (old==0){ if (pre.pool==FALSE) { if (match.arg(classif)=="time") { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb1, p.last, 1) } else { if (no.n > no.pre) { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1+no.pre, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb1+no.pre, p.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(p1, i.last, 0) Ind.list.surv[[which(names(Ind.list.surv)=="pb.am")]][2,] <- c(p1, i.last, 1) } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(nb.last, p.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(p1, i.last, 0) Ind.list.surv[[which(names(Ind.list.surv)=="pb.am")]][2,] <- c(p1, i.last, 1) }} } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(p.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(p.last, p.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(p1, i.last, 0) Ind.list.surv[[which(names(Ind.list.surv)=="pb.am")]][2,] <- c(p1, i.last, 1) }} else { Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][1,] <- c(nb1, p.last, 0) Ind.list.surv$old[1,] <- c(nb1, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", length(p.surv.sub), sep=""))]][2,] <- c(p.last, p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", length(pb.pre.sub), sep=""))]][1,] <- c(p.last, p.last, 1) Ind.list.surv$old[2,] <- c(p.last, p.last, 1) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(p1, i.last, 0) Ind.list.surv[[which(names(Ind.list.surv)=="pb.am")]][2,] <- c(p1, i.last, 1) } if (no.pre>1) { temp.indmat[,] <- NA temp.indmat[1:(no.pre-1),] <- cbind((p1+1):(p.last), p1:(p.last-1), 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } temp.indmat[1:(no.pre-1),3] <- 1 for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }} # If age matters (both B and NB age-classified): if (no.b>1 & match.arg(classif)=="age"){ # Transitions from PB (or juv) to B if (no.pre==no.b) { temp.indmat[,] <- NA temp.indmat[1:no.b,] <- cbind(c((b1+1):b.last, b.last), p1:p.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- c(b1, i.last, 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } if (old>0) {Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(b.last, p.last, 1) } } else { Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- c(b1, i.last, 0) if (no.pre>0){ temp.indmat[,] <- NA temp.indmat[1:no.pre,] <- cbind(c((b1+1):(am+no.pre*2)), p1:p.last, 0) for (i in 1:length(p.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:length(pb.pre.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } }} # Transitions from NB to B if (old==0) { temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((b1+1):b.last, b.last), nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][1:n.surv.sub[i],] <- temp.indmat[n.s.start[i]:n.s.end[i],] } for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][1:pb.n.sub[i],] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } } # Transitions from B to NB if (old==0) { temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((nb1+1):nb.last, nb.last), b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][1:b.surv.sub[i],] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][1:pb.b.sub[i],] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } else{ Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][,1]))),] <- c(nb1, b.last, 0) Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(nb1, b.last, 0) # Transitions from B to PB in old system temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((p1+1):p.last, p.last), b1:b.last, 0) for (i in 1:(length(b.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:(length(pb.b.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(p.last, b.last, 1) } } # If time matters for B or NB (or they are both 1) if (match.arg(classif)=="time"){ # Transitions from PB (or immature) to B Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(b1, i.last, 0) Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- c(b1, i.last, 0) if (no.pre>0){ temp.indmat[,] <- NA temp.indmat[1:no.pre,] <- cbind(b1, p1:p.last, 0) for (i in 1:(length(p.surv.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[p.s.start[i]:p.s.end[i],] } for (i in 1:(length(pb.pre.sub))){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[pb.pre.start[i]:pb.pre.end[i],] } } # Transitions from NB to B temp.indmat[,] <- NA temp.indmat[1:(no.n),] <- cbind(b1, nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1])))+n.surv.sub[i]-1),] <- temp.indmat[n.s.start[i]:n.s.end[i],] } for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1])))+pb.n.sub[i]-1),] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } # Transitions from B to NB temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(nb1, b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } temp.indmat[1:no.b,3] <- 1 for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } # Survival and probability of breeders staying in B if (no.b>1){ temp.indmat[,] <- NA temp.indmat[1:(no.b),] <- cbind(c((b1+1):b.last, b.last), b1:b.last, 0) for (i in 1:length(b.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- temp.indmat[b.s.start[i]:b.s.end[i],] } for (i in 1:length(pb.b.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[pb.b.start[i]:pb.b.end[i],] } } else {Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", length(b.surv.sub), sep=""))]][,1]))),] <- c(b.last, b.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", length(pb.b.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", length(pb.b.sub), sep=""))]][,1]))),] <- c(b.last, b.last, 0) } if (old>0) {Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(b.last, b.last, 1) } # Survival and probability of nonbreeders staying in NB if (no.n>1){ temp.indmat[,] <- NA temp.indmat[1:(no.n),] <- cbind(c((nb1+1):nb.last, nb.last), nb1:nb.last, 0) for (i in 1:length(n.surv.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1])))+n.surv.sub[i]-1),] <- temp.indmat[n.s.start[i]:n.s.end[i],] } temp.indmat[1:no.n,3] <- 1 temp.indmat[no.n,] <- NA for (i in 1:length(pb.n.sub)){Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1])))+pb.n.sub[i]-1),] <- temp.indmat[pb.n.start[i]:pb.n.end[i],] } if (old==0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))),] <- c(nb.last, nb.last, 1) } else{ Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", length(n.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", length(n.surv.sub), sep=""))]][,1]))),] <- c(nb.last, nb.last, 0) Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", length(pb.n.sub), sep=""))]][,1]))),] <- c(nb.last, nb.last, 1) } # Fecundity temp.indmat[,] <- NA if (old>0){temp.indmat[1:(b.last-am+2),] <- cbind(1, (am-1):b.last, 0)} else { temp.indmat[1:(nb.last-am+2),] <- cbind(1, (am-1):nb.last, 0)} Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("i.surv", length(i.surv.sub), sep=""))]][,1]))),] <- c(1, am-1, 0) for (i in 1:length(b.surv.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("b.surv", i, sep=""))]][,1])))+b.surv.sub[i]-1),] <- c(temp.indmat[(b.s.start[i]+1+no.pre):(b.s.end[i]+1+no.pre),])} for (i in 1:length(pb.b.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.b", i, sep=""))]][,1])))+pb.b.sub[i]-1),] <- temp.indmat[(pb.b.start[i]+1+no.pre):(pb.b.end[i]+1+no.pre),]} if (old==0) {for (i in 1:length(n.surv.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("n.surv", i, sep=""))]][,1])))+n.surv.sub[i]-1),] <- temp.indmat[(n.s.start[i]+1+no.pre+no.b):(n.s.end[i]+1+no.pre+no.b),]} for (i in 1:length(pb.n.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.n", i, sep=""))]][,1])))+pb.n.sub[i]-1),] <- temp.indmat[(pb.n.start[i]+1+no.pre+no.b):(pb.n.end[i]+1+no.pre+no.b),]} } Ind.list.surv$pb.am[min(which(is.na(Ind.list.surv$pb.am[,1]))),] <- temp.indmat[1,] if (no.pre>0){ for (i in 1:length(p.surv.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("p.surv", i, sep=""))]][,1])))+p.surv.sub[i]-1),] <- temp.indmat[(p.s.start[i]+1):(p.s.end[i]+1),]} for (i in 1:length(pb.pre.sub)) {Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1]))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("pb.pre", i, sep=""))]][,1])))+pb.pre.sub[i]-1),] <- temp.indmat[(pb.pre.start[i]+1):(pb.pre.end[i]+1),]} } if (match.arg(classif)=="time"){ Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+no.pre),] <- temp.indmat[1:(1+no.pre),] Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+fec.sub[1]-1),] <- temp.indmat[(2+no.pre):(1+no.pre+fec.sub[1]),] if (old==0) {Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+no.n-1),] <- temp.indmat[(2+no.pre+no.b):(1+no.pre+no.b+no.n),]} if (length(fec.sub)>1) {for (i in 2:length(fec.sub)) Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[(1+no.pre+fec.start[i]):(1+no.pre+fec.end[i]),]} } if (match.arg(classif)=="age"){ Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))),] <- c(1, am-1, 0) Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+fec.sub[1]-1),] <- temp.indmat[(2+no.pre):(1+no.pre+fec.sub[1]),] if (length(fec.sub)>1) {for (i in 2:length(fec.sub)) Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[(1+no.pre+fec.start[i]):(1+no.pre+fec.end[i]),]} if (old==0){ Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+fec.sub[1]-1),] <- temp.indmat[(2+no.pre+no.b):(1+no.pre+no.b+fec.sub[1]),] if (length(fec.sub)>1) {for (i in 2:length(fec.sub)) Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[(1+no.pre+no.b+fec.start[i]):(1+no.pre+no.b+fec.end[i]),]}} if (no.pre>0){ Ind.list.surv$fec1[(min(which(is.na(Ind.list.surv$fec1[,1])))):(min(which(is.na(Ind.list.surv$fec1[,1])))+fec.sub[1]-1),] <- temp.indmat[2:(1+fec.sub[1]),] if (length(fec.sub)>1) {for (i in 2:length(fec.sub)) Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))):(min(which(is.na(Ind.list.surv[[which(names(Ind.list.surv)==paste("fec", i, sep=""))]][,1])))+fec.sub[i]-1),] <- temp.indmat[(1+fec.start[i]):(1+fec.end[i]),]} } if (old>0) { Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(1, p.last, 1) Ind.list.surv$old[min(which(is.na(Ind.list.surv$old[,1]))),] <- c(1, b.last, 1) } } # Create separate matrices for each of the parameters par.mats <- lapply(1:length(par.sub.names), function(x) matrix(NA, nrow=nb.last, ncol=nb.last)) names(par.mats) <- paste(par.sub.names, ".mat", sep="") for (i in 1:length(par.sub.names)) { par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,2])] <- par.vals[i] par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,2])] <- (1-par.vals[i])} # Replace NA values with 1, or 0 if they are all NA (this makes it possible to multiply them together properly) unl.Ind.s <- lapply(1:length(par.mats), function(x) unlist(par.mats[x])) surv.NAs <- lapply(1:(nb.last*nb.last), function(x) sapply(1:length(par.mats), function(y) is.na(unl.Ind.s[[y]][x]))) surv.all.NAs <- lapply(1:length(surv.NAs), function(x) all(surv.NAs[[x]]==TRUE)) for (i in 1:length(par.sub.names)) {par.mats[[i]][which(surv.all.NAs==TRUE)] <- 0 par.mats[[i]][is.na(par.mats[[i]])] <- 1 } # Multiply together all matrices to create a single matrix for (i in 2:length(par.mats)) {if (i==2) {NMs <- par.mats[[1]]} NMs <- NMs*par.mats[[i]]} NM <- NMs # Row names Classes <- which(c(am-1, no.pre, no.b, no.n)>0) name.let <- list("i%d", "p%d", "b%d", "nb%d") name.num <- list(1:(am-1), 1:no.pre, 1:no.b, 1:no.n) Nclasses <- c(unlist(mapply(sprintf, name.let[Classes], name.num[Classes]))) rownames(NM) <- Nclasses # Frequency dependent case ------------------------------------------------------------------------------------------------ if (match.arg(freq.dep)!="none"){ # Create matrix to record population projection popsize <- matrix(nrow=nb.last, ncol=times+1) # Insert initial values popsize[,1] <- startvec # Define function for updating frequency dependent matrix fdMat <- function(tt, popsize=popsize){ if (old>0 | no.pre>0) {noNB <- sum(popsize[c(p1:p.last, nb1:nb.last),tt-1])} else { noNB <- sum(popsize[c(nb1:nb.last),tt-1]) } noB <- sum(popsize[c(b1:b.last), tt-1]) if (freq.dep=="fec"){ # Calculate breeder survival adjusted by nonbreeder density fec.fd <- fec/(1+noNB/noB) # Insert fec.fd into matrix for (i in grep("fec", par.sub.names)) {par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,2])] <- fec.fd[i-min(grep("fec", par.sub.names))+1] par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,2])] <- (1-fec.fd[i-min(grep("fec", par.sub.names))+1])} for (i in 2:length(par.mats)) {if (i==2) {NMs <- par.mats[[1]]} NMs <- NMs*par.mats[[i]]} NMfd <- NMs} if (freq.dep=="surv"){ # Calculate breeder survival adjusted by nonbreeder density b.surv.fd <- b.surv/(1+noNB/noB) # Insert b.surv.fd into matrix for (i in grep("b.surv", par.sub.names)) {par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==0,2])] <- b.surv.fd[i-min(grep("b.surv", par.sub.names))+1] par.mats[[i]][cbind(Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,1], Ind.list.surv[[i]][is.na(Ind.list.surv[[i]][,1])==FALSE & Ind.list.surv[[i]][,3]==1,2])] <- (1-b.surv.fd[i-min(grep("b.surv", par.sub.names))+1])} for (i in 2:length(par.mats)) {if (i==2) {NMs <- par.mats[[1]]} NMs <- NMs*par.mats[[i]]} NMfd <- NMs} return(NMfd)} # Project population over "times" time steps for (tt in 2:(times+1)){ NMfd <- fdMat(tt, popsize) # Project population at next time step newpopsize <- NMfd %*% popsize[,tt-1] popsize[,tt] <- newpopsize } NMfd.eq <- fdMat(times+1, popsize) rownames(NMfd.eq) <- rownames(NM) if (old>0){populationsize <- popsize[1:b.last,] NMfd.equil <- NMfd.eq[-nrow(NMfd.eq), -nrow(NMfd.eq)]} else { populationsize <- popsize NMfd.equil <- NMfd.eq} phat <- populationsize[,times]/sum(populationsize[,times]) phat.test <- (NMfd.equil %*% phat) / norm(NMfd.equil %*% phat) test <- phat.test - phat if (any(abs(test)>0.000000001)) stop("not at equilibrium") # Find deterministic growth rate eig_info <- eigen(NMfd.equil) ev <- which.max(Re(eig_info$values)) eig_val <- Re(eig_info$values[ev]) Lambda.test <- (colSums(populationsize)[times+1]/colSums(populationsize)[times-9])^(1/10) if (abs(Lambda.test-eig_val)>0.000000001) stop("not at equilibrium") } else { # Find stable stage structure and reproductive values ------------------------------- # For "old" system use matrix without senescent class if (old>0) {NM <- NM[-nrow(NM), -nrow(NM)]} # Define function Eig.calc <- function(m.matrix){ # Stable stage structure eig_info <- eigen(m.matrix) ev <- which.max(Re(eig_info$values)) eig_val <- Re(eig_info$values[ev]) u_vals <- Re((eig_info$vectors[,ev])/sum(eig_info$vectors[,ev])) # Reproductive values eig_info2 <- eigen(t(m.matrix)) ev2 <- which(round(eig_val, digits=5)==round(eig_info2$values, digits=5)) v_vals <- Re(eig_info2$vectors[,ev2])/sum(u_vals*Re(eig_info2$vectors[,ev2])) # Test eigenvalues uvector.test <- m.matrix%*%u_vals - eig_val*u_vals if (any(uvector.test>10^-4)) stop("problem with u_vals") vvector.test <- t(v_vals)%*%m.matrix - eig_val*t(v_vals) if (any(vvector.test>10^-4)) stop("problem with v_vals") # if (any(vvector.test>10^-4)) v_vals <- NA # With the old system there will very occasionally be a problem with the v-values (only shows up when running many iterations). This line lets it continue running anyway. Use with care. return(list(eig_val, u_vals, v_vals))} # Run function on NM matrix NM.info <- Eig.calc(NM) eig_val <- NM.info[[1]] u_vals <- NM.info[[2]] v_vals <- NM.info[[3]] if (old>0) Nclasses <- Nclasses[-length(Nclasses)] info.mat <- data.frame(Nclasses, u_vals, v_vals) # Caculate sensitivities ------------------------------------------------ if (sens==TRUE){ # Matrix element level sensitivities tot.NMsens <- sensitivity(NM, zero=TRUE) # Create matrix to store lower level sensitivities in sens.mat <- matrix(nrow=length(names(Ind.list.surv)), ncol=1, dimnames=list(par.sub.names)) for (i in 1:length(names(Ind.list.surv))) { Ind.no <- length(which(is.na(Ind.list.surv[[i]][,1])==FALSE)) if (Ind.no>0) { temp.sens.ind <- rep(NA, Ind.no) for (j in 1:Ind.no) { if (Ind.list.surv[[i]][j, 3]==0){ temp.sens.ind[j] <- tot.NMsens[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]*NMs[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]/par.vals[i] } else { temp.sens.ind[j] <- -1*tot.NMsens[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]*NMs[Ind.list.surv[[i]][j, 1], Ind.list.surv[[i]][j, 2]]/(1-par.vals[i]) } } sens.mat[i] <- sum(temp.sens.ind) }} } else {tot.NMsens <- NA sens.mat <- NA} } # Construct and analyze Leslie matrix ------------------------------------------ if (Leslie==TRUE){ # Not intended for use with old # In "recent breeding history" system. First create expanded matrix with both age and breeding structure. if (match.arg(classif)=="time" & no.pre==0 & (no.b>1 | no.n>1)) { agecl.matrix <- matrix(0, nrow=nrow(NM)+(no.n+no.b)*(no.n+no.b-1), ncol=nrow(NM)+(no.n+no.b)*(no.n+no.b-1)) rownames(agecl.matrix) <- c(1:(am-1), rep(am:(am-1+no.n+no.b), each=(no.n+no.b))) #Rownames indicate ages colnames(agecl.matrix) <- c(row.names(NM)[1:(am-1)], rep(row.names(NM)[am:dim(NM)[1]], times=(no.n+no.b))) #Colnames indicate breeding class max.age.class <- ifelse(time.ages==0, am+no.n+no.b-1, am+time.ages-1) # Fecundity agecl.matrix[1:(am-1),1:nrow(NM)] <- NM[1:(am-1),] age.blocks <- seq(nrow(NM)+1, nrow(agecl.matrix), no.n+no.b) for (i in age.blocks){ agecl.matrix[1:(am-1),i:(i+(no.n+no.b-1))] <- NM[1:(am-1), am:nrow(NM)]} # Survival agecl.matrix[am:nrow(NM), 1:(am-1)] <- NM[am:(nrow(NM)),1:(am-1)] for (i in age.blocks){ agecl.matrix[i:(i+(no.n+no.b-1)), (i-(no.n+no.b)):(i-1)] <- NM[am:nrow(NM), am:nrow(NM)] } agecl.matrix[age.blocks[length(age.blocks)]:(age.blocks[length(age.blocks)]+(no.n+no.b-1)),age.blocks[length(age.blocks)]:(age.blocks[length(age.blocks)]+(no.n+no.b-1))] <- NM[am:nrow(NM), am:nrow(NM)] # Analyze the expanded matrix agecl.info <- Eig.calc(agecl.matrix) agecl.u_vals <- agecl.info[[2]] # Construct Leslie matrix directly from expanded matrix (using stable age distribution) umage <- t(apply(agecl.matrix, 1, function(x) x*agecl.u_vals)) L.matrixB.t <- matrix(NA, ncol=max.age.class, nrow=max.age.class) for (i in 1:(max.age.class-1)){L.matrixB.t[i, max.age.class] <- sum(umage[which(rownames(agecl.matrix)==i), which(rownames(agecl.matrix)>=max.age.class)])} for (j in 1:(max.age.class-1)){L.matrixB.t[max.age.class, j] <- sum(umage[which(rownames(agecl.matrix)>=max.age.class), which(rownames(agecl.matrix)==j)])} L.matrixB.t[max.age.class, max.age.class] <- sum(umage[which(rownames(agecl.matrix)>=max.age.class), which(rownames(agecl.matrix)>=max.age.class)]) u.agecl.m <- rep(NA, max.age.class) for (i in 1:(max.age.class-1)){u.agecl.m[i] <- sum(agecl.u_vals[which(rownames(agecl.matrix)==i)])} u.agecl.m[max.age.class] <- sum(agecl.u_vals[which(rownames(agecl.matrix)>=max.age.class)]) L.matrixB <- t(apply(L.matrixB.t, 1, function(x) x/u.agecl.m)) L.matrix <- L.matrixB } else { # Stable age distribution u.breed <- u_vals[b1:b.last] u.non <- u_vals[nb1:nb.last] if (no.pre>0) {u.pre <- c(u_vals[p1:p.last],rep(0, no.b-no.pre))} else u.pre <- 0 u.juv <- u_vals[1:(am-1)] u.age.ad <- u.breed+u.non+u.pre u.age <- c(u.juv, u.age.ad) # sorting classes by age j.age.index <- lapply(1:i.last, function(z) z) if (no.pre>0) {p.age.index <- lapply(am:p.last, function(x) c(x, p.last+x-am+1, b.last+x-am+1))} else {p.age.index <- NULL} if (no.b > no.pre) {a.age.index <- lapply(1:(no.b-no.pre), function(y) c(p.last+no.pre+y, b.last+no.pre+y) )} else {a.age.index <- NULL} age.index <- c(j.age.index, p.age.index, a.age.index) # Construct Leslie matrix directly from NM matrix (using stable age distribution) uNM <- t(apply(NM, 1, function(x) x*u_vals)) L.matrix.t <- matrix(NA, ncol=length(u.age), nrow=length(u.age)) for (i in 1:length(u.age)){ for (j in 1:length(u.age)){ L.matrix.t[i,j] <- sum(uNM[age.index[[i]], age.index[[j]]]) } } L.matrixB <- t(apply(L.matrix.t, 1, function(x) x/u.age)) # Calculate fecundity for the Leslie model (including survival) ifelse(no.pre==0, L.fec <- c(rep(0, am-2), i.surv.v[max(!is.na(i.surv.v))]*pb.am*fec.v[1], (b.surv.v*pb.b.v*u.breed+n.surv.v*pb.n.v*u.non)*fec.v[-1]/u.age.ad), L.fec <- c(rep(0, am-2), i.surv.v[max(!is.na(i.surv.v))]*pb.am*fec.v[1], (b.surv.v*pb.b.v*u.breed+n.surv.v*pb.n.v*u.non+p.surv.v*pb.pre.v*u.pre)*fec.v[-1]/u.age.ad)) # Calculate survival for the Leslie model if (no.pre>0) {pre.term <- p.surv.v*u.pre} else pre.term <- 0 if (am>2) {L.surv <- c(i.surv.v[-length(i.surv.v)], i.surv.v[max(!is.na(i.surv.v))], (b.surv.v*u.breed+n.surv.v*u.non+pre.term)/u.age.ad)} else { L.surv <- c(i.surv.v[max(!is.na(i.surv.v))], (b.surv.v*u.breed+n.surv.v*u.non+pre.term)/u.age.ad) } # Construct Leslie matrix using calculated parameters L.matrix <- matrix(0, ncol=am-1+no.b, nrow=am-1+no.b) L.matrix[1,] <- L.fec L.matrix[cbind(2:(am-1+no.b), 1:(am-2+no.b))] <- L.surv[-(am-1+no.b)] L.matrix[(am-1+no.b), (am-1+no.b)] <- L.surv[am-1+no.b] # Check that the two Leslie matrices are the same if (any(abs(L.matrix-L.matrixB)>0.000001)) stop("Problem with L.matrix") } # Growth rate, stable stage distribution and reproductive value L.info <- Eig.calc(L.matrix) L.eig_val <- L.info[[1]] L.u_vals <- L.info[[2]] L.v_vals <- L.info[[3]] L.DemVar <- NA L.DemVar.cov <- NA } else {L.DemVar <- NA L.DemVar.cov <- NA L.eig_val <- NA L.info <- NA L.u_vals <- NA L.v_vals <- NA} # Model of breeding population ------------------------------------------ # All breeders are seen, but nonbreeders are not. All offspring are marked. if (ignoN.model==TRUE){ # 1 means breeder, 2 means nonbreeder # Define function for calculating probability of different paths through life cycle. path.probs <- function(cvec){ prob <- ifelse (cvec[1]==1, prod(i.surv.v)*pb.am, prod(i.surv.v)*(1-pb.am)) for (x in 1:(length(cvec)-1)) { iprobs <- matrix(NA, nrow=2, ncol=2) iprobs[1,1] <- b.surv.v[x]*pb.b.v[x] iprobs[1,2] <- b.surv.v[x]*(1-pb.b.v[x]) iprobs[2,1] <- n.surv.v[x]*pb.n.v[x] iprobs[2,2] <- n.surv.v[x]*(1-pb.n.v[x]) prob <- prob*iprobs[cvec[x], cvec[x+1]] } prob <- ifelse(cvec[length(cvec)]==1, prob*pb.b.v[length(cvec)]*b.surv.v[length(cvec)], prob*pb.n.v[length(cvec)]*n.surv.v[length(cvec)]) return(prob) } # All possible paths onetwo <- list(1:2) combs <- rep(onetwo, max.age-am+1) all.combs <- cbind(matrix(2, ncol=am-1), expand.grid(combs)) #Columns of twos make numbering of columns easier for variable am # Calculating survival estimate up to focal age surv.to.age <- vector("numeric", length=max.age-am+1) for (i in am:max.age) { if (i==am) {focal.probs <- prod(i.surv.v)*pb.am} else{ if (i==(am+1)) {focal.probs <- prod(i.surv.v)*pb.am*b.surv.v[1]*pb.b.v[1]+prod(i.surv.v)*(1-pb.am)*n.surv.v[1]*pb.n.v[1]} else { focal.combs <- expand.grid(rep(onetwo, i-am)) focal.probs <- apply(focal.combs, 1, function(x) path.probs(x))}} if (i