CSEGriskfigure <- function(data, te = 100, absolutethresh = FALSE, threshold = 0.1, datalogged = FALSE, silent = FALSE, return.model = FALSE, CI.method = "hessian", CI.sim = 1000) {
if (!(CI.method %in% c("hessian", "parametric", "innovations", "none"))) {
stop("Stopped in CSEGriskfigure(): Allowed values of 'CI.method' are \"none\", \"hessian\", \"parametric\", and \"innovations\".\n", call. = FALSE)
}
if (!silent) {
if (datalogged) {
cat("Analysis assumes that data and threshold are already logged")
} else {
cat("Data and threshold are assumed to be unlogged")
}
cat("\n")
if (absolutethresh) {
cat("Using an absolute threshold for the extinction threshold")
cat("\n")
cat("Threshold is ")
cat(threshold)
cat("\n")
if (datalogged) thresh <- exp(threshold)
}
else {
cat("Using an percentage threshold for the extinction threshold")
cat("\n")
cat("Threshold is ")
cat((1 - threshold) * 100)
cat(" percent decline")
cat("\n")
}
}
par(mfrow = c(3, 2))
# set up pi
pi <- get("pi", pos = "package:base")
# a=read.table(data, skip=dataskip)
a <- as.matrix(data)
if (dim(a)[2] != 2) stop("Stopped in CSEGriskfigure() because this function requires a data file with 2 columns; time in first column.", call. = FALSE)
if (datalogged) a[, 2] <- exp(a[, 2]) # plot unlogged
nyr <- length(a[, 1])
# estimate parameters
dat <- a[, 2]
dat <- log(dat) # because I exp-ed it above if logged
kem <- MARSS(dat, model = list(U = matrix("U")), silent = TRUE)
if (kem$par$Q != 0) kem <- MARSSparamCIs(kem)
# PANEL #1 plot the data
plot(a[, 1], a[, 2],
type = "p", bty = "L", xaxp = c(a[1, 1], a[nyr, 1], nyr - 1),
xlab = "", ylab = "Pop. Estimate", ylim = c(.9 * min(a[, 2], na.rm = TRUE), 1.1 * max(a[, 2], na.rm = TRUE)), xlim = c(a[1, 1] - 1, a[nyr, 1] + 1)
)
lines(a[, 1], exp(kem$states[1, ]), col = "red")
title(
paste("u est = ", format(kem$par$U, digits = 2), " (95% CIs ", format(kem$par.lowCI[["U"]], digits = 2), ", ", format(kem$par.upCI[["U"]], digits = 2), ")", "\n", "Q est = ", format(kem$par$Q, digits = 2)),
cex.main = .9
)
# PANEL #2 plot the cdf of extinction
tyrs <- 1:te
kal.cdf <- matrix(0, nrow = te)
N0 <- exp(kem$states[1, length(dat)])
if (!absolutethresh) thresh <- threshold * N0 else thresh <- threshold
xd <- log(N0) - log(thresh)
kal.u <- as.vector(kem$par$U)
kal.Q <- as.vector(kem$par$Q) # to create scalar from 1x1 matrix
if (kal.u <= 0) p.ever <- 1 else p.ever <- exp(-2 * kal.u * xd / kal.Q)
if (N0 <= thresh) {
kal.cdf <- rep(1, te)
} else {
if (is.finite(exp(2 * xd * abs(kal.u) / kal.Q))) {
sec.part <- exp(2 * xd * abs(kal.u) / kal.Q) * pnorm((-xd - abs(kal.u) * tyrs) / sqrt(kal.Q * tyrs))
} else {
sec.part <- 0
}
kal.cdf <- p.ever * (pnorm((-xd + abs(kal.u) * tyrs) / sqrt(kal.Q * tyrs)) + sec.part)
}
plot(tyrs, kal.cdf, xlab = "time steps into future", ylab = "probability to hit threshold", ylim = c(0, 1), bty = "l", type = "l")
title(paste("Prob. to hit ", format(thresh, digits = 1)))
if (kem$par$Q != 0) {
# Add CIs to cdf of quasi-extinction plot if requested
nsim <- CI.sim
plotCI <- TRUE
if (CI.method == "hessian") {
plotCI <- FALSE
kem <- MARSShessian(kem)
# The R and Q are chol transformed so = sqrt(R) and sqrt(Q) respectively
boot.params <- try(mvtnorm::rmvnorm(n = nsim, mean = kem$parMean, sigma = kem$parSigma, method = "chol"), silent = TRUE)
if (!inherits(boot.params, "try-error")) {
plotCI <- TRUE
boot.params <- boot.params[boot.params[, which(substr(colnames(boot.params), 1, 1) == "Q")] > 0, ]
boot.params[, which(substr(colnames(boot.params), 1, 1) == "Q")] <- boot.params[, which(substr(colnames(boot.params), 1, 1) == "Q")]^2
boot.params[, which(substr(colnames(boot.params), 1, 1) == "R")] <- boot.params[, which(substr(colnames(boot.params), 1, 1) == "R")]^2
}
}
if (CI.method == "parametric") {
boot.params <- t(MARSSboot(kem, nboot = 1000, output = "parameters", sim = "parametric")$boot.params)
}
if (CI.method == "innovations") {
boot.params <- t(MARSSboot(kem, nboot = 1000, output = "parameters", sim = "nonparametric")$boot.params)
}
if (CI.method != "none" & plotCI) {
nsim <- dim(boot.params)[1]
kal.cdf.cis <- matrix(0, nrow = nsim, ncol = te)
for (j in 1:nsim) {
u <- boot.params[j, which(substr(colnames(boot.params), 1, 1) == "U")]
Q <- boot.params[j, which(substr(colnames(boot.params), 1, 1) == "Q")]
if (u <= 0) {
p.ever <- 1
} else {
p.ever <- exp(-2 * u * xd / Q)
}
if (is.finite(exp(2 * xd * abs(u) / Q))) {
sec.part <- exp(2 * xd * abs(u) / Q) * pnorm((-xd - abs(u) * tyrs) / sqrt(Q * tyrs))
} else {
sec.part <- 0
}
kal.cdf.cis[j, ] <- p.ever * pnorm((-xd + abs(u) * tyrs) / sqrt(Q * tyrs)) + sec.part
}
CIs.ex <- apply(kal.cdf.cis, 2, quantile, c(0.0275, .125, .875, .9725), na.rm = TRUE)
lines(tyrs, CIs.ex[1, ], col = "red")
lines(tyrs, CIs.ex[4, ], col = "red")
lines(tyrs, CIs.ex[2, ], col = "green", lty = 2)
lines(tyrs, CIs.ex[3, ], col = "green", lty = 2)
legend("topleft", c("95% CI", "75% CI", "mean"), col = c("red", "green", "black"), lty = c(1, 2, 1), bty = "n")
}
# Panel 3 Now make plot of pdf of the time to extinction
kal.pdf <- matrix(nrow = te * 5)
tyrs <- 1:(te * 5)
if (N0 <= thresh) {
kal.pdf[1] <- 1
} else {
for (i in 1:(te * 5)) {
kal.pdf[i] <-
xd * (1 / sqrt(2 * pi * kal.Q * tyrs[i]^3)) *
exp((-(xd - abs(kal.u) * tyrs[i])^2) / (2 * kal.Q * tyrs[i]))
} # end i loop
}
i <- tyrs
aa <- max(kal.pdf)
i.left <- 1:which(kal.pdf == aa)
i.min <- max(which(kal.pdf[i.left] < aa * .001))
i.right <- which(kal.pdf == aa):max(tyrs)
if (min(kal.pdf[i.right]) < (aa * .001)) {
i.max <- which(kal.pdf == aa) - 1 + min(which(kal.pdf[i.right] < aa * .001))
} else {
i.max <- max(tyrs)
}
plot(tyrs[i.min:i.max], kal.pdf[i.min:i.max], xlab = "time steps into future", ylab = "probability to hit threshold", bty = "l", type = "l")
title(main = "PDF of time to threshold \n given it IS reached")
# Panel 4 Now make plot of probability of hitting different thresholds in te yrs
thshes <- seq(0.5 * thresh, 2 * thresh, (2 * thresh - 0.5 * thresh) / 100)
kal.thresh <- matrix(0, nrow = length(thshes))
for (i in 1:length(thshes)) {
xd <- log(N0) - log(thshes[i])
if (kal.u <= 0) {
p.ever <- 1
} else {
p.ever <- exp(-2 * kal.u * xd / kal.Q)
}
if (N0 <= thshes[i]) {
kal.thresh[i] <- 1
} else {
kal.thresh[i] <-
p.ever * pnorm((-xd + abs(kal.u) * te) / sqrt(kal.Q * te))
}
+exp(2 * xd * abs(kal.u) / kal.Q) * pnorm((-xd - abs(kal.u) * te) / sqrt(kal.Q * te))
} # end i loop
plot(thshes, kal.thresh, xlab = "Number of ind. at Ne", ylab = "probability to hit threshold", bty = "l", ylim = c(0, 1), type = "l")
abline(v = thresh, col = "red")
text(thresh * 1.1, max(kal.thresh) * .9, "90% threshold")
title(paste("Prob. of hitting threshold", " in ", te, " time steps"))
N <- matrix(0, ncol = te, nrow = 10)
for (i in 1:10) {
N[i, ] <- N0 * cumprod(exp(rnorm(te, kal.u, sqrt(kal.Q))))
}
graphics::matplot(t(N), type = "l", bty = "L", main = "Sample projections", xlab = "time steps into the future", ylab = "N")
# Finally make TMU plot
CSEGtmufigure(N = nyr, u = kal.u, s2p = kal.Q, make.legend = FALSE)
# KW
par(mfrow = c(1, 1))
if (return.model) {
return(kem)
}
} else {
cat("Extinction PDF plots not possible because Q estimate is 0.\n")
if (return.model) {
return(kem)
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.