Annotation of ext/extinction_risk_1.R, Revision 1.22
1.7 hako 1: # extinction_risk_1.R
1.22 ! hako 2: # $Id: extinction_risk_1.R,v 1.21 2015/08/16 13:59:30 hako Exp $
1.1 hako 3: #
4: # Author: Hiroshi Hakoyama <hako@affrc.go.jp>
5: # Copyright (c) 2013-2015 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
6: #
7: # Redistribution and use in source and binary forms, with or without
8: # modification, are permitted provided that the following conditions
9: # are met:
10: # 1. Redistributions of source code must retain the above copyright
11: # notice, this list of conditions and the following disclaimer.
12: # 2. Redistributions in binary form must reproduce the above copyright
13: # notice, this list of conditions and the following disclaimer in the
14: # documentation and/or other materials provided with the distribution.
15: #
16: # THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY
17: # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
18: # THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
19: # PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
20: # AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21: # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22: # NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23: # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
24: # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
25: # STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26: # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
27: # ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28: #
29: # Description:
1.8 hako 30: # Estimates the demographic parameters and the probability of
31: # extinction within a specific time period, t, from a time series
32: # of population size based on the Wiener-drift process model.
33: # An estimator for confidence interval of extinction risk
34: # developed by Hakoyama (w-z method, in preparation) is implemented,
35: # which is better than the estimator from the ordinary Delta method.
1.1 hako 36: #
37: # Usage:
1.20 hako 38: # ext1(dat, tt = 100, ne = 1, alpha = 0.05, verbose = FALSE,
1.15 hako 39: # exact.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, qq.plot = FALSE, formatted = TRUE)
1.1 hako 40: #
41: # Arguments:
1.13 hako 42: # dat data.frame of 2 variables: Time and Population size
1.20 hako 43: # tt a time period of interest
44: # ne a lower (extinction) threshold of population size (ne >= 1)
1.8 hako 45: # alpha 1 - confidence level
1.13 hako 46: # verbose If set to FALSE, give the ML estimate of the probability
1.20 hako 47: # of extinction within a specific time period tt. If set to
1.8 hako 48: # TRUE, give a more verbose output:
49: # the ML estimate of the growth rate (mu),
50: # the ML estimate of variance (s),
51: # Current population size, xd = log(N(q) / ne),
52: # Sample size, n = q + 1,
1.19 hako 53: # AIC for the distribution of X = log(N),
1.8 hako 54: # the ML estimate of the probability of extinction within
55: # a specific time period t (P),
56: # a lower alpha % confidence limit of parameter *
57: # (lower.CL.*), and
58: # an upper alpha % confidence limit of parameter *
59: # (upper.CL.*).
1.10 hako 60: # exact.CL If set to FALSE, give a CI by the w-z method.
61: # If set to TRUE, give an asymptotically exact CI.
62: # n.sim The number of iteration for the asymptotically exact CI.
63: # alpha.sim Level of significance to decide the convergence of iteration.
1.13 hako 64: # qq.plot If set to TRUE, give a QQ plot for delta.log.n.
1.14 hako 65: # formatted If set to TRUE, give the result by formatted output.
1.8 hako 66: # If set to FALSE, give a list of estimates.
1.1 hako 67: # References:
1.8 hako 68: # R. Lande and S. H. Orzack. Extinction dynamics of age-structured
69: # populations in a fluctuating environment. Proceedings of the
70: # National Academy of Sciences, 85(19):7418-7421, 1988.
71: # B. Dennis, P. L. Munholland, and J. M. Scott. Estimation of
72: # growth and extinction parameters for endangered species.
73: # Ecological Monographs, 61:115-143, 1991.
74: # H. Hakoyama (in preparation).
1.1 hako 75: #
76:
1.20 hako 77: ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05, verbose = FALSE, exact.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, qq.plot = FALSE, formatted = TRUE) {
1.1 hako 78: yr <- ts(dat[, 1], start = c(dat[, 1][1])) # Year
79: ps <- ts(dat[, 2], start = c(dat[, 1][1])) # Population size
80: complete <- complete.cases(yr, ps)
81: yr <- yr[complete]
82: ps <- ps[complete]
83: tau <- diff(yr) # time intervals, \tau_i = t_i - t_{i-1}
1.8 hako 84: delta.log.n <- diff(log(ps))
1.20 hako 85: qq <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\}
1.1 hako 86: tq <- sum(tau)
1.20 hako 87: mu <- sum(delta.log.n) / tq # MLE of growth rate
88: s <- sum((delta.log.n - mu * tau)^2 / tau) / qq # MLE of variance
89: us <- qq * s / (qq - 1) # an unbiased estimate of variance
1.1 hako 90: xd <- log(ps[length(ps)] / ne)
1.20 hako 91: l <- - log(2 * pi * s) * qq / 2 - sum((delta.log.n - mu)^2) / (2 * s)
92: AIC.X <- - 2 * l + 4 # AIC for the distribution of X = log(N)
1.17 hako 93:
1.13 hako 94: if (qq.plot == TRUE) {
95: qqnorm(delta.log.n)
96: qqline(delta.log.n)
97: }
98:
1.20 hako 99: lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq)
100: upper.CL.mu <- mu + qt(1 - alpha / 2, qq - 1) * sqrt(us / tq)
101: lower.CL.s <- qq * s / qchisq(1 - alpha / 2, qq - 1)
102: upper.CL.s <- qq * s / qchisq(alpha / 2, qq - 1)
103:
104: ww <- W(mu, xd, s, tt)
105: zz <- Z(mu, xd, s, tt)
106: pp <- Pr(ww, zz)
1.8 hako 107:
1.20 hako 108: if (exact.CL == TRUE) {
109: CL.P <- AsymptoticallyExactConfidenceInterval(mu, xd, s, tt, tq, qq, alpha, n.sim, alpha.sim)
1.10 hako 110: } else {
1.20 hako 111: CL.P <- ConfidenceInterval(mu, xd, s, tt, tq, qq, alpha)
1.10 hako 112: }
113:
1.8 hako 114: lower.CL.P <- CL.P[[1]]
115: upper.CL.P <- CL.P[[2]]
116:
1.1 hako 117: if (verbose == TRUE) {
1.20 hako 118: results <- list(ne = ne, tt = tt, verbose = verbose,
119: AIC.X = AIC.X,
120: n = qq + 1,
1.6 hako 121: xd = xd,
1.1 hako 122: Growth.rate = mu,
1.8 hako 123: lower.CL.mu = lower.CL.mu,
124: upper.CL.mu = upper.CL.mu,
1.1 hako 125: Variance = s,
1.8 hako 126: lower.CL.s = lower.CL.s,
127: upper.CL.s = upper.CL.s,
128: # Unbiased.variance = us,
129: Extinction.probability = pp,
130: lower.CL.P = lower.CL.P,
131: upper.CL.P = upper.CL.P)
1.4 hako 132: if (formatted == TRUE) {
133: class(results) <- "ext1"
134: }
1.1 hako 135: return(results)
136: } else {
1.20 hako 137: results <- list(ne = ne, tt = tt, verbose = verbose,
1.8 hako 138: Extinction.probability = pp)
1.4 hako 139: if (formatted == TRUE) {
140: class(results) <- "ext1"
141: }
1.1 hako 142: return(results)
143: }
144: }
145:
1.20 hako 146: #-------------------------------------------------------------------
147: W <- function(mu, xd, s, tt) (mu * tt + xd) / sqrt(s * tt)
148: #-------------------------------------------------------------------
149: Z <- function(mu, xd, s, tt) (- mu * tt + xd) / sqrt(s * tt)
150: #-------------------------------------------------------------------
151: Pr <- function(w, z) {
152: if(z < 35) {
153: pnorm(-w) + exp((z^2 - w^2) / 2) * pnorm(-z)
154: } else {
155: pnorm(-w) + exp(- w^2 / 2) * (sqrt(2) / (2 * sqrt(pi))) *
156: (1 / z - 1 / z^3 + 3 / z^5 - 15 / z^7
157: + 105 / z^9 - 945 / z^11 + 10395 / z^13)
158: }
159: }
160: #-------------------------------------------------------------------
161: ConfidenceInterval <- function(mu, xd, s, tt, tq, qq, alpha) {
162: den1 <- sqrt(s * tt)
163: w.est <- (mu * tt + xd) / den1
164: z.est <- (- mu * tt + xd) / den1
1.22 ! hako 165: df <- qq - 1
! 166: const.1 <- sqrt((qq - 1) * tq / (qq * tt))
! 167: const.2 <- sqrt(tq / tt)
! 168: FindCL <- function(tq, qq, tt, est, a, width = 10) {
! 169: t.obs <- const.1 * est
! 170: f <- function(x) pt(t.obs, df, const.2 * x) - a
! 171: d.est <- est / width + 1
! 172: uniroot(f, c(- d.est + est, d.est + est), extendInt = "yes")$root
! 173: }
1.20 hako 174: lower.CL.w <- FindCL(tq, qq, tt, w.est, 1 - alpha / 2)
175: upper.CL.w <- FindCL(tq, qq, tt, w.est, alpha / 2)
176: lower.CL.z <- FindCL(tq, qq, tt, z.est, 1 - alpha / 2)
177: upper.CL.z <- FindCL(tq, qq, tt, z.est, alpha / 2)
178: lower.CL.P <- Pr(upper.CL.w, lower.CL.z)
179: upper.CL.P <- Pr(lower.CL.w, upper.CL.z)
180: c(lower.CL.P, upper.CL.P)
181: }
182: #-------------------------------------------------------------------
183: randomCI <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function = ConfidenceInterval) {
184: mu <- rnorm(1, mean = mu.obs, sd = sqrt(s.obs / qq))
185: s <- rchisq(1, df = qq - 1) * s.obs / qq
186: ci.function(mu, xd, s, tt, tq, qq, alpha)
187: }
188: #-------------------------------------------------------------------
189: RejectionRate <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, gam, n.sim, ci.function = ConfidenceInterval) {
190: CI.dist <- replicate(n.sim, randomCI(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function))
191: P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt))
192: ci <- ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha)
193: complete <- complete.cases(CI.dist[1, ], CI.dist[2, ])
194: lower <- CI.dist[1, ][complete]
195: upper <- CI.dist[2, ][complete]
196: n.estimables <- length(lower)
197: n.rejects <- (sum(P.obs < lower) + sum(P.obs > upper))
198: if (n.estimables > 0) {
199: binom <- binom.test(n.rejects, n.estimables, p = gam)
200: } else {
201: binom <- list(estimate = NaN, p.value = NaN)
202: }
203: c(binom$estimate[[1]], binom$p.value, ci)
204: }
205: #-------------------------------------------------------------------
1.22 ! hako 206: AsymptoticallyExactConfidenceInterval <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, n.sim, alpha.sim, ci.function = ConfidenceInterval, verbose = FALSE) {
1.20 hako 207: P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt))
208: if (P.obs == 0 || P.obs == 1) {
1.21 hako 209: return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha))
1.20 hako 210: }
211: res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function)
212: gg <- res[[1]]
213: pp <- res[[2]]
1.22 ! hako 214: if (verbose == TRUE) cat(alpha, gg, "\n")
! 215: ai <- alpha^2 / gg
! 216: i <- 0
1.20 hako 217: while (pp < alpha.sim) {
1.22 ! hako 218: i <- i + 1
1.20 hako 219: res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, ai, alpha, n.sim, ci.function)
220: gg <- res[[1]]
221: pp <- res[[2]]
1.22 ! hako 222: if (verbose == TRUE) cat(ai, gg, "\n")
! 223: ai <- ai * alpha / gg
! 224: }
! 225: if (verbose == TRUE) {
! 226: c(res[[3]], res[[4]], ai, i)
! 227: } else {
! 228: c(res[[3]], res[[4]])
1.20 hako 229: }
230: }
231: #-------------------------------------------------------------------
1.1 hako 232: print.ext1 <- function(obj, digits = 5) {
233: if (obj$verbose == TRUE) {
234: output.est <- data.frame(
235: c(formatC(obj$Growth.rate, digits = digits),
236: formatC(obj$Variance, digits = digits),
1.8 hako 237: formatC(obj$xd, digits = digits),
238: formatC(obj$n, digits = digits),
1.20 hako 239: formatC(obj$AIC.X, digits = digits),
1.1 hako 240: formatC(obj$Extinction.probability, digits = digits)),
1.8 hako 241: c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",
242: formatC(obj$upper.CL.mu, digits = digits),")", sep = ""),
243: paste("(",formatC(obj$lower.CL.s, digits = digits),", ",
244: formatC(obj$upper.CL.s, digits = digits),")", sep = ""),
245: "-",
246: "-",
1.16 hako 247: "-",
1.8 hako 248: paste("(",formatC(obj$lower.CL.P, digits = digits),", ",
249: formatC(obj$upper.CL.P, digits = digits),")", sep = "")))
1.1 hako 250: dimnames(output.est) <- list(
251: c("Growth rate:",
252: "Variance:",
1.20 hako 253: "Current population size, x_d = log(N(q) / n_e):",
1.8 hako 254: "Sample size, n = q + 1:",
1.19 hako 255: "AIC for the distribution of X = log(N):",
1.20 hako 256: paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")),
1.1 hako 257: c("Estimate","95% confidence interval"))
258: print(output.est)
259: } else {
260: output.est <- data.frame(
261: c(formatC(obj$Extinction.probability, digits = digits)))
262: dimnames(output.est) <- list(
1.20 hako 263: c(paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")),
1.1 hako 264: c("Estimate"))
1.13 hako 265: print(output.est)
1.1 hako 266: }
267: }
268: #
269: # Examples
270: # Yellowstone grizzly bears (from Dennis et al., 1991)
1.21 hako 271: dat <- data.frame(Year = c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987),
272: Population = c(44, 47, 46, 44, 46, 45, 46, 40, 39, 39, 42, 44, 41, 40, 33, 36, 34, 39, 35, 34, 38, 36, 37, 41, 39, 51, 47, 57, 47))
1.1 hako 273: # The probability of extinction (of decline to population size 1) within 100 years
1.8 hako 274: # ext1(dat, t = 100)
1.1 hako 275: # The probability of decline to 10 individuals within 100 years
1.8 hako 276: # ext1(dat, t = 100, ne = 10, verbose = TRUE)
1.13 hako 277: # with QQ-plot
278: # ext1(dat, t = 100, ne = 10, verbose = TRUE, qq.plot = T)
1.10 hako 279: # The probability of decline to 1 individuals within 100 years with an asymptotically exact confidence interval
280: # ext1(dat, t = 100, ne = 1, verbose = TRUE, exact.CL = TRUE)
![[BACK]](/icons/back.gif){kind=icon}
Return to extinction_risk_1.R CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / ext