Supplement to Stamation et al. (2020) - Endang Species Res 41:373-383 - https://doi.org/10.3354/esr01031


Supplement 3

## R script for theoretical population model for SRW ##

# Constants to be used in the model

Surv2r=0.988; seSurv2r=0.001; # Survival in females 2+ years (Brando, 2013)

Surv1c=0.87; seSurv1c=0.17; # Survival first year (Carroll, 2016)

Surv2c=0.95; seSurv2c=0.05; # Survival 2-8 years (Carroll, 2016)

Mature=9; # Age at which females reach maturity (Carroll, 2016)

# Function to simulate life history for whales (births, sex and deaths).

WhaleBirths <- function(Born, Mature=FALSE, End=2018){

  if(Mature){

    Whale <- data.frame(Year=1980:End, Born=Born, Sex=1, Alive=1, Mature=NA, Calf=NA)

    Whale$Mature <- ifelse(Whale$Year >= Born+9, 1, 0)

  } else {

    Whale <- data.frame(Year=1980:End, Born=Born, Sex=NA, Alive=NA, Mature=NA, Calf=NA)

    Whale$Sex <- rbinom(1,1,0.5)

  }

  # print(c(Born, Whale$Sex[1]))

  if(Whale$Sex[1]==1){

    if(!Mature){

      Whale$Alive[which(Whale$Year==Whale$Born)] <- 1

      Whale$Mature[which(Whale$Year==Whale$Born)] <- 0

      if(which(Whale$Year==Whale$Born) < nrow(Whale)){

        Whale$Alive[which(Whale$Year==Whale$Born)+1] <- rbinom(1,1,Surv1c)

        if(Whale$Alive[which(Whale$Year==Whale$Born)+1]==1){

          Whale$Mature[which(Whale$Year==Whale$Born)+1] <- 0

        }

      }

      if(which(Whale$Year==Whale$Born)+1 < nrow(Whale)){

        for(i in (which(Whale$Year==Whale$Born)+2):min(which(Whale$Year==Whale$Born)+7, nrow(Whale))){

          Whale$Alive[i] <- rbinom(1,1,Surv2c)*Whale$Alive[i-1]

          if(Whale$Alive[i]==1){

            Whale$Mature[i] <- 0

          }

        }

      }

    }

    if(Whale$Born[1]+8 <= max(Whale$Year)){

      j <- which(Whale$Year==Whale$Born+8)

      while(j <= nrow(Whale)){

        Whale$Alive[j] <- rbinom(1,1,Surv2r)*Whale$Alive[j-1]

        if(Whale$Alive[j]==1){

          Whale$Mature[j] <- 1

        }

        j <- j + 1

      }

    }

    if(max(Whale$Mature, na.rm=TRUE)==1){

    Whale$Calf[which(Whale$Year==min(Whale$Year[which(Whale$Mature==1)])):nrow(Whale)] <- 0

      Whale$Calf[which(Whale$Year==min(Whale$Year[which(Whale$Mature==1)]))] <- 1

      Birth1 <- CalfIntData$Time %*% rmultinom(1, size=1, prob=CalfIntData$Prob)

      while(which.max(Whale$Calf*Whale$Year)+Birth1 <= nrow(Whale)){

        Whale$Calf[which.max(Whale$Calf*Whale$Year)+Birth1] <- 1

        Birth1 <- CalfIntData$Time %*% rmultinom(1, size=1, prob=CalfIntData$Prob)

      }

      Whale$Calf[Whale$Alive==0] <- NA

    }

  } else{

    ny <- which(Whale$Year >= Whale$Born)

    surv <- c(1, Surv1c, rep(Surv2c, 7), rep(Surv2r, nrow(Whale)-9))[1:length(ny)]

    live <- rbinom(length(ny), 1, surv)

    kk <- 1

    for(k in surv){

      Whale$Alive[k] <- prod(live[1:kk])

      }

  }

  return(Whale)

}

# Function to collate the data from the life history simulator for individual whales

WhaleTotal <- function(Years, End=2018){

  Births <- data.frame(Year=Years, Sex=1, Calves=NA, Death=NA)

  i <- 1

  while(nrow(Births) >= i){

    # print(i)

    data <- WhaleBirths(Births$Year[i], ifelse(length(Years)>=i, TRUE, FALSE), End)

    Births$Sex[i] <- data$Sex[1]

    Births$Calves[i] <- ifelse(Births$Sex[i]==1, sum(data$Calf, na.rm=TRUE), 0)

    Births$Death[i] <- max(data$Year*data$Alive, na.rm=TRUE)

    if(Births$Calves[i] > 0){

      Births <- rbind(Births, data.frame(Year=data$Year[which(data$Calf==1)], Sex=NA, Calves=NA, Death=NA))

    }

    i <- i + 1

  }

  return(Births)

}

Between <- function(x, a, b){

  sum(a <= x & b >= x)

}

# Function to collate individuals life histories into the group as whole over time.

WhaleSim <- function(data, Initial=1973:1978, End=2018, nsim=1000){

  Yearly <- array(NA, c(length(1980:End), 5, nrow(data)),

                  dimnames=list(NULL, c('Year', 'Mean', 'Median', 'LB','UB'), NULL))

  Yearly[,1,] <- 1980:End

  for(i in 1:nrow(data)){

    nCalf <- integer(length = nsim)

    nPop <- matrix(NA, nrow=nsim, ncol=nrow(Yearly))

    for(j in 1:nsim){

      Years <- sample(Initial, size=data$Pop0[i], replace=TRUE)

      Try <- WhaleTotal(Years, End=End)

      nCalf[j] <- sum(Try$Calves)

      for(jj in 1:nrow(Yearly)){

        nPop[j,jj] <- sum(Try$Year <= Yearly[jj,1,i] & Try$Death >= Yearly[jj,1,i])

      }

    }

    data$Mean[i] <- mean(nCalf)

    data[i,3:5] <- quantile(nCalf, probs=c(0.5, 0.025, 0.975))

    Yearly[,2,i] <- t(colMeans(nPop))

    Yearly[,3:5,i] <- t(apply(nPop, 2, MedCI))

  }

  return(list(data, Yearly))

}