Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@cjabradshaw
cjabradshaw committed Oct 6, 2021
1 parent 60b324f commit cb03a14
Showing 1 changed file with 319 additions and 0 deletions.
319 changes: 319 additions & 0 deletions WombatPopModel.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,319 @@
##########################################################################################################################################
## common wombat (Vombatus ursinus) demographic model
##
## Corey Bradshaw
## [email protected]
## Flinders University, September 2021
##########################################################################################################################################

## functions
# beta distribution shape parameter estimator function
estBetaParams <- function(mu, var) {
alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
beta <- alpha * (1 / mu - 1)
return(params = list(alpha = alpha, beta = beta))
}

AICc <- function(...) {
models <- list(...)
num.mod <- length(models)
AICcs <- numeric(num.mod)
ns <- numeric(num.mod)
ks <- numeric(num.mod)
AICc.vec <- rep(0,num.mod)
for (i in 1:num.mod) {
if (length(models[[i]]$df.residual) == 0) n <- models[[i]]$dims$N else n <- length(models[[i]]$residuals)
if (length(models[[i]]$df.residual) == 0) k <- sum(models[[i]]$dims$ncol) else k <- (length(models[[i]]$coeff))+1
AICcs[i] <- (-2*logLik(models[[i]])) + ((2*k*n)/(n-k-1))
ns[i] <- n
ks[i] <- k
AICc.vec[i] <- AICcs[i]
}
return(AICc.vec)
}

delta.AIC <- function(x) x - min(x) ## where x is a vector of AIC
weight.AIC <- function(x) (exp(-0.5*x))/sum(exp(-0.5*x)) ## Where x is a vector of dAIC
ch.dev <- function(x) ((( as.numeric(x$null.deviance) - as.numeric(x$deviance) )/ as.numeric(x$null.deviance))*100) ## % change in deviance, where x is glm object

linreg.ER <- function(x,y) { # where x and y are vectors of the same length; calls AICc, delta.AIC, weight.AIC functions
fit.full <- lm(y ~ x); fit.null <- lm(y ~ 1)
AIC.vec <- c(AICc(fit.full),AICc(fit.null))
dAIC.vec <- delta.AIC(AIC.vec); wAIC.vec <- weight.AIC(dAIC.vec)
ER <- wAIC.vec[1]/wAIC.vec[2]
r.sq.adj <- as.numeric(summary(fit.full)[9])
return(c(ER,r.sq.adj))
}

## source
source("matrixOperators.r")


##############################
## VOMBATUS (ursinus) (VU)
## sources: Roger et al. 2011 Popul Ecol 53:215-227

# mass
VU.mass <- 25 # Vombatus ursinus (Saran et al. 2011 Pacific Conservation Biology 17:310-319)

## predicted rm (from Henneman 1983 Oecologia 56:104-108)
## log10rm = 0.6914 - 0.2622*log10m (mass in g)
VU.rm.pred <- 10^(0.6914 - (0.2622*log10(VU.mass*1000)))
VU.lm.pred <- exp(VU.rm.pred)

## theoretical population density for mammalian herbivores based on body size (Damuth 1981; Freeland 1990)
## log10D = 4.196 − 0.74*(log10m)
VU.D.pred <- (10^(4.196 - (0.74*log10(VU.mass*1000))))/2 # divided by 2 for females only
VU.D.pred # animals/km2

## max age
## non-volant birds & mammals (Healy K et al. 2014 PRSB)
## log10ls = 0.89 + 0.13log10m (mass in grams; ls = years)
VU.age.max <- round(10^(0.89 + (0.13*log10(VU.mass*1000))), 0)
VU.age.max <- 26 # reset based on McIlroy 2008

## age vector
VU.age.vec <- 0:VU.age.max

## fertility
## total fecundity from Allainé et al. 1987 (Oecologia)
## lnF = 2.719 - 0.211lnM (all mammals)
VU.F.pred <- exp(2.719 - (0.211*log(VU.mass*1000)))/2 # divided by 2 for females
VU.F.pred <- 1/2/2 # inter-birth interval of 2 years / 2 for daughters only (McIlroy 1995)

## age at primiparity
## lnalpha = 0.214 + 0.263*lnM (https://dx.doi.org/10.1093%2Fgerona%2F62.2.149)
VU.alpha <- ceiling(exp(-1.34 + (0.214*log(VU.mass*1000))))
VU.alpha <- 2 # reset according to Roger et al. (2011)

## define m function with age
## pouch young per year
VU.pypy <- 0.5

## proportion of females breeding each year
VU.pfbr <- 0.84

## proportion young that are female
VU.sr <- 0.5

## m vector
VU.m.vec <- c(0, 0.3*VU.F.pred, 0.9*VU.F.pred, rep(VU.F.pred, 24))
VU.m.sd.vec <- 0.05*VU.m.vec
plot(VU.age.vec, VU.m.vec, type="b", pch=19, xlab="age (yrs)", ylab="m")

# fit sigmoidal function
# logistic power function y = a / (1+(x/b)^c)
VU.m.dat <- data.frame(VU.age.vec, VU.m.vec)
param.init <- c(0.2, 2, -3)
VU.fit.logp <- nls(VU.m.vec ~ a / (1+(VU.age.vec/b)^c),
data = VU.m.dat,
algorithm = "port",
start = c(a = param.init[1], b = param.init[2], c = param.init[3]),
trace = TRUE,
nls.control(maxiter = 1000, tol = 1e-05, minFactor = 1/1024))
VU.fit.logp.summ <- summary(VU.fit.logp)
plot(VU.age.vec, VU.m.vec, type="b", pch=19, xlab="age (yrs)", ylab="m")
VU.age.vec.cont <- seq(0,max(VU.age.vec),1)
VU.pred.p.m <- coef(VU.fit.logp)[1] / (1+(VU.age.vec.cont/coef(VU.fit.logp)[2])^coef(VU.fit.logp)[3])
VU.pred.p.mm <- ifelse(VU.pred.p.m > 1, 1, VU.pred.p.m)
VU.pred.p.mm[2] <- 0
lines(VU.age.vec.cont, VU.pred.p.mm,lty=2,lwd=3,col="red")

## survival
## mean adult survival (McCarthy et al. 2008 Am Nat)
## ln{-ln[s(t)]} = ln(a) + bln(M) + ln (t)
ln.a.s <- -0.5; b.s <- -0.25
VU.s.tran <- ln.a.s + b.s*log(VU.mass*1000) + log(1)
VU.s.ad.yr <- exp(-exp(VU.s.tran))

# Siler hazard h(x) (Gurven et al. 2007)
a1 <- 1 - (1.05*VU.s.ad.yr) # initial infant mortality rate (also known as αt)
b1 <- 3.1 # rate of mortality decline (also known as bt)
a2 <- 1 - VU.s.ad.yr # age-independent mortality (exogenous mortality due to environment); also known as ct
a3 <- 0.1e-04 # initial adult mortality rate (also known as βt)
b3 <- 0.02 # rate of mortality increase
longev <- VU.age.max
x <- seq(0,longev,1) # age vector
h.x <- a1 * exp(-b1*x) + a2 + a3 * exp(b3 * x) # Siler's hazard model
plot(x,h.x,pch=19,type="l")
plot(x,log(h.x),pch=19,type="l")
l.x <- exp((-a1/b1) * (1 - exp(-b1*x))) * exp(-a2 * x) * exp(a3/b3 * (1 - exp(b3 * x))) # Siler's survival (proportion surviving) model
init.pop <- 10000
lx <- round(init.pop*l.x,0)
len.lx <- length(lx)
dx <- lx[1:(len.lx-1)]-lx[2:len.lx]
qx <- dx/lx[1:(length(lx)-1)]
VU.Sx <- c(0.99*VU.s.ad.yr, 1 - qx)
plot(x, VU.Sx, pch=19, type="l", xlab="age (years)", ylab="Sx")
VU.s.sd.vec <- 0.05*VU.Sx

## create matrix
VU.popmat <- matrix(data = 0, nrow=VU.age.max+1, ncol=VU.age.max+1)
diag(VU.popmat[2:(VU.age.max+1),]) <- VU.Sx[-(VU.age.max+1)]
VU.popmat[VU.age.max+1,VU.age.max+1] <- 0
VU.popmat[1,] <- VU.pred.p.mm
colnames(VU.popmat) <- c(0:VU.age.max)
rownames(VU.popmat) <- c(0:VU.age.max)
VU.popmat.orig <- VU.popmat ## save original matrix

## matrix properties
max.lambda(VU.popmat.orig) ## 1-yr lambda
VU.lm.pred
max.r(VU.popmat.orig) # rate of population change, 1-yr
VU.ssd <- stable.stage.dist(VU.popmat.orig) ## stable stage distribution
plot(VU.age.vec, VU.ssd, type="l", pch=19, xlab="age (yrs)", ylab="ssd")
R.val(VU.popmat.orig, VU.age.max) # reproductive value
VU.gen.l <- G.val(VU.popmat.orig, VU.age.max) # mean generation length

## initial population vector
area <- 500*500 # km × km
VU.pop.found <- round(area*VU.D.pred, 0) # founding population size (estimated density * 100 × 100 km region [10,000 km2])
VU.init.vec <- VU.ssd * VU.pop.found

#################
## project
## set time limit for projection in 1-yr increments
yr.st <- 1
#************************
yr.end <- round(40*VU.gen.l, 0) # set projection end date
#************************
t <- (yr.end - yr.st)

VU.tot.F <- sum(VU.popmat.orig[1,])
VU.popmat <- VU.popmat.orig
yr.vec <- seq(yr.st,yr.end)

## set population storage matrices
VU.n.mat <- matrix(0, nrow=VU.age.max+1,ncol=(t+1))
VU.n.mat[,1] <- VU.init.vec

## set up projection loop
for (i in 1:t) {
VU.n.mat[,i+1] <- VU.popmat %*% VU.n.mat[,i]
}

VU.n.pred <- colSums(VU.n.mat)
yrs <- seq(yr.st, yr.end, 1)
plot(yrs, log10(VU.n.pred),type="l",lty=2,pch=19,xlab="year",ylab="log10 N")

# compensatory density feedback
VU.K.max <- 1*VU.pop.found
VU.K.vec <- c(1, VU.K.max/2, 0.75*VU.K.max, VU.K.max)
VU.red.vec <- c(1,0.985,0.951,0.873)
plot(VU.K.vec, VU.red.vec,pch=19,type="b")
VU.Kred.dat <- data.frame(VU.K.vec, VU.red.vec)

# logistic power function a/(1+(x/b)^c)
VU.param.init <- c(1, 2*VU.K.max, 2)
VU.fit.lp <- nls(VU.red.vec ~ a/(1+(VU.K.vec/b)^c),
data = VU.Kred.dat,
algorithm = "port",
start = c(a = VU.param.init[1], b = VU.param.init[2], c = VU.param.init[3]),
trace = TRUE,
nls.control(maxiter = 1000, tol = 1e-05, minFactor = 1/1024))
VU.fit.lp.summ <- summary(VU.fit.lp)
plot(VU.K.vec, VU.red.vec, pch=19,xlab="N",ylab="reduction factor")
VU.K.vec.cont <- seq(1,2*VU.pop.found,1)
VU.pred.lp.fx <- coef(VU.fit.lp)[1]/(1+(VU.K.vec.cont/coef(VU.fit.lp)[2])^coef(VU.fit.lp)[3])
lines(VU.K.vec.cont, VU.pred.lp.fx, lty=3,lwd=3,col="red")

VU.a.lp <- coef(VU.fit.lp)[1]
VU.b.lp <- coef(VU.fit.lp)[2]
VU.c.lp <- coef(VU.fit.lp)[3]

## compensatory density-feedback deterministic model
## set population storage matrices
VU.n.mat <- matrix(0, nrow=VU.age.max+1, ncol=(t+1))
VU.n.mat[,1] <- VU.init.vec
VU.popmat <- VU.popmat.orig

## set up projection loop
for (i in 1:t) {
VU.totN.i <- sum(VU.n.mat[,i])
VU.pred.red <- as.numeric(VU.a.lp/(1+(VU.totN.i/VU.b.lp)^VU.c.lp))
diag(VU.popmat[2:(VU.age.max+1),]) <- (VU.Sx[-(VU.age.max+1)])*VU.pred.red
VU.popmat[VU.age.max+1,VU.age.max+1] <- 0
VU.popmat[1,] <- VU.pred.p.mm
VU.n.mat[,i+1] <- VU.popmat %*% VU.n.mat[,i]
}

VU.n.pred <- colSums(VU.n.mat)
plot(yrs, VU.n.pred, type="l",lty=2,pch=19,xlab="year",ylab="N")
abline(h=VU.pop.found, lty=2, col="red", lwd=2)

## stochatic projection with density feedback
## set storage matrices & vectors
iter <- 100
itdiv <- iter/10

VU.n.sums.mat <- matrix(data=NA, nrow=iter, ncol=(t+1))
VU.s.arr <- VU.m.arr <- array(data=NA, dim=c(t+1, VU.age.max+1, iter))

for (e in 1:iter) {
VU.popmat <- VU.popmat.orig

VU.n.mat <- matrix(0, nrow=VU.age.max+1,ncol=(t+1))
VU.n.mat[,1] <- VU.init.vec

for (i in 1:t) {
# stochastic survival values
VU.s.alpha <- estBetaParams(VU.Sx, VU.s.sd.vec^2)$alpha
VU.s.beta <- estBetaParams(VU.Sx, VU.s.sd.vec^2)$beta
VU.s.stoch <- rbeta(length(VU.s.alpha), VU.s.alpha, VU.s.beta)

# stochastic fertilty sampler (gaussian)
VU.fert.stch <- rnorm(length(VU.popmat[,1]), VU.pred.p.mm, VU.m.sd.vec)
VU.m.arr[i,,e] <- ifelse(VU.fert.stch < 0, 0, VU.fert.stch)

VU.totN.i <- sum(VU.n.mat[,i], na.rm=T)
VU.pred.red <- VU.a.lp/(1+(VU.totN.i/VU.b.lp)^VU.c.lp)

diag(VU.popmat[2:(VU.age.max+1),]) <- (VU.s.stoch[-(VU.age.max+1)])*VU.pred.red
VU.popmat[VU.age.max+1,VU.age.max+1] <- 0
VU.popmat[1,] <- VU.m.arr[i,,e]
VU.n.mat[,i+1] <- VU.popmat %*% VU.n.mat[,i]

VU.s.arr[i,,e] <- VU.s.stoch * VU.pred.red

} # end i loop

VU.n.sums.mat[e,] <- ((as.vector(colSums(VU.n.mat))/VU.pop.found))

if (e %% itdiv==0) print(e)

} # end e loop

VU.n.md <- apply(VU.n.sums.mat, MARGIN=2, median, na.rm=T) # mean over all iterations
VU.n.up <- apply(VU.n.sums.mat, MARGIN=2, quantile, probs=0.975, na.rm=T) # upper over all iterations
VU.n.lo <- apply(VU.n.sums.mat, MARGIN=2, quantile, probs=0.025, na.rm=T) # lower over all iterations

par(mfrow=c(1,3))
plot(yrs,VU.n.md,type="l", main = "", xlab="year", ylab="pN1", lwd=2, ylim=c(0.95*min(VU.n.lo),1.05*max(VU.n.up)))
lines(yrs,VU.n.lo,lty=2,col="red",lwd=1.5)
lines(yrs,VU.n.up,lty=2,col="red",lwd=1.5)

VU.s.add <- VU.m.add <- rep(0, VU.age.max+1)
for (m in 1:iter) {
VU.s.add <- rbind(VU.s.add, VU.s.arr[ceiling(VU.gen.l):(t+1),,m])
VU.m.add <- rbind(VU.m.add, VU.m.arr[ceiling(VU.gen.l):(t+1),,m])
}
VU.s.add <- VU.s.add[-1,]
VU.m.add <- VU.m.add[-1,]

VU.s.md <- apply(VU.s.add, MARGIN=2, median, na.rm=T) # mean s over all iterations
VU.s.up <- apply(VU.s.add, MARGIN=2, quantile, probs=0.975, na.rm=T) # upper over all iterations
VU.s.lo <- apply(VU.s.add, MARGIN=2, quantile, probs=0.025, na.rm=T) # lower over all iterations

plot(VU.age.vec,VU.s.md,type="l", main = "", xlab="age", ylab="s", lwd=2, ylim=c(0.95*min(VU.s.lo),1.05*max(VU.s.up)))
lines(VU.age.vec,VU.s.lo,lty=2,col="red",lwd=1.5)
lines(VU.age.vec,VU.s.up,lty=2,col="red",lwd=1.5)

VU.m.md <- apply(VU.m.add, MARGIN=2, median, na.rm=T) # mean s over all iterations
VU.m.up <- apply(VU.m.add, MARGIN=2, quantile, probs=0.975, na.rm=T) # upper over all iterations
VU.m.lo <- apply(VU.m.add, MARGIN=2, quantile, probs=0.025, na.rm=T) # lower over all iterations

plot(VU.age.vec,VU.m.md,type="l", main = "", xlab="age", ylab="m", lwd=2, ylim=c(0.95*min(VU.m.lo),1.05*max(VU.m.up)))
lines(VU.age.vec,VU.m.lo,lty=2,col="red",lwd=1.5)
lines(VU.age.vec,VU.m.up,lty=2,col="red",lwd=1.5)
par(mfrow=c(1,1))

0 comments on commit cb03a14

Please sign in to comment.