Annotation of ext/extinction_risk_1.R, Revision 1.9
1.7 hako 1: # extinction_risk_1.R
1.9 ! hako 2: # $Id: extinction_risk_1.R,v 1.8 2015/08/02 12:41:02 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.8 hako 38: # ext1(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE,
39: # formatted = TRUE)
1.1 hako 40: #
41: # Arguments:
42: # dat data.frame of 2 variables: Time and Population size
43: # t a time period of interest
44: # ne a lower (extinction) threshold of population size
1.8 hako 45: # alpha 1 - confidence level
46: # verbose If set to FALSE, give the ML estimate of the probability
47: # of extinction within a specific time period t. If set to
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,
53: # the ML estimate of the probability of extinction within
54: # a specific time period t (P),
55: # a lower alpha % confidence limit of parameter *
56: # (lower.CL.*), and
57: # an upper alpha % confidence limit of parameter *
58: # (upper.CL.*).
59: # formatted If set to TRUE, give the result by formatted output.
60: # If set to FALSE, give a list of estimates.
1.1 hako 61: #
62: # References:
1.8 hako 63: # R. Lande and S. H. Orzack. Extinction dynamics of age-structured
64: # populations in a fluctuating environment. Proceedings of the
65: # National Academy of Sciences, 85(19):7418-7421, 1988.
66: # B. Dennis, P. L. Munholland, and J. M. Scott. Estimation of
67: # growth and extinction parameters for endangered species.
68: # Ecological Monographs, 61:115-143, 1991.
69: # H. Hakoyama (in preparation).
1.1 hako 70: #
71:
1.4 hako 72: ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE, formatted = TRUE) {
1.1 hako 73: yr <- ts(dat[, 1], start = c(dat[, 1][1])) # Year
74: ps <- ts(dat[, 2], start = c(dat[, 1][1])) # Population size
75: complete <- complete.cases(yr, ps)
76: yr <- yr[complete]
77: ps <- ps[complete]
78: tau <- diff(yr) # time intervals, \tau_i = t_i - t_{i-1}
1.8 hako 79: delta.log.n <- diff(log(ps))
1.1 hako 80: q <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\}
81: tq <- sum(tau)
1.8 hako 82: mu <- sum(delta.log.n) / tq # ML estimate of growth rate
83: s <- (1 / q) * sum((delta.log.n - mu * tau)^2 / tau) # ML estimate of variance
1.1 hako 84: us <- q * s / (q - 1) # an unbiased estimate of variance
85: xd <- log(ps[length(ps)] / ne)
1.8 hako 86:
87: lower.CL.mu <- mu - qt(1 - alpha / 2, q - 1) * sqrt(us / tq)
88: upper.CL.mu <- mu + qt(1 - alpha / 2, q - 1) * sqrt(us / tq)
89: lower.CL.s <- q * s / qchisq(1 - alpha / 2, q - 1)
90: upper.CL.s <- q * s / qchisq(alpha / 2, q - 1)
91:
92: w <- function(mu, xd, s, t) (mu * t + xd) / sqrt(s * t)
93: z <- function(mu, xd, s, t) (- mu * t + xd) / sqrt(s * t)
94: pr <- function(w, z) {
95: if(z < 35) {
96: pnorm(-w) + exp((z^2 - w^2) / 2) * pnorm(-z)
97: } else {
1.9 ! hako 98: pnorm(-w) + exp(- w^2 / 2) * (sqrt(2) / (2 * sqrt(pi))) * (1 / z - 1 / z^3 + 3 / z^5 - 15 / z^7 + 15 * 7 / z^9 - 15 * 7 * 9 / z^11 + 15 * 7 * 9 * 11 / z^13)
1.8 hako 99: }
1.1 hako 100: }
1.8 hako 101:
102: ww <- w(mu, xd, s, t)
103: zz <- z(mu, xd, s, t)
104: pp <- pr(ww, zz)
105:
106: c.limit <- function(tq, q, t, est, a, width = 10) {
107: t.obs <- sqrt((q - 1) * tq / (q * t)) * est
108: df <- q - 1
109: const <- sqrt(tq / t)
110: f <- function(x) {
111: pt(t.obs, df, const * x) - a
112: }
113: d.est <- est / width + 1
114: uniroot(f, c(- d.est + est, d.est + est), extendInt = "yes")$root
1.1 hako 115: }
1.8 hako 116:
117: confidence.interval <- function(mu, xd, s, t, tq, q, alpha) {
118: den1 <- sqrt(s * t)
119: w.est <- (mu * t + xd) / den1
120: z.est <- (- mu * t + xd) / den1
121: lower.CL.w <- c.limit(tq, q, t, w.est, 1 - alpha / 2)
122: upper.CL.w <- c.limit(tq, q, t, w.est, alpha / 2)
123: lower.CL.z <- c.limit(tq, q, t, z.est, 1 - alpha / 2)
124: upper.CL.z <- c.limit(tq, q, t, z.est, alpha / 2)
125: lower.CL.P <- pr(upper.CL.w, lower.CL.z)
126: upper.CL.P <- pr(lower.CL.w, upper.CL.z)
127: c(lower.CL.P, upper.CL.P)
1.1 hako 128: }
1.8 hako 129:
130: CL.P <- confidence.interval(mu, xd, s, t, tq, q, alpha)
131: lower.CL.P <- CL.P[[1]]
132: upper.CL.P <- CL.P[[2]]
133:
1.1 hako 134: if (verbose == TRUE) {
135: results <- list(ne = ne, t = t, verbose = verbose,
1.8 hako 136: n = q + 1,
1.6 hako 137: xd = xd,
1.1 hako 138: Growth.rate = mu,
1.8 hako 139: lower.CL.mu = lower.CL.mu,
140: upper.CL.mu = upper.CL.mu,
1.1 hako 141: Variance = s,
1.8 hako 142: lower.CL.s = lower.CL.s,
143: upper.CL.s = upper.CL.s,
144: # Unbiased.variance = us,
145: Extinction.probability = pp,
146: lower.CL.P = lower.CL.P,
147: upper.CL.P = upper.CL.P)
1.4 hako 148: if (formatted == TRUE) {
149: class(results) <- "ext1"
150: }
1.1 hako 151: return(results)
152: } else {
153: results <- list(ne = ne, t = t, verbose = verbose,
1.8 hako 154: Extinction.probability = pp)
1.4 hako 155: if (formatted == TRUE) {
156: class(results) <- "ext1"
157: }
1.1 hako 158: return(results)
159: }
160: }
161:
162: print.ext1 <- function(obj, digits = 5) {
163: if (obj$verbose == TRUE) {
164: output.est <- data.frame(
165: c(formatC(obj$Growth.rate, digits = digits),
166: formatC(obj$Variance, digits = digits),
1.8 hako 167: formatC(obj$xd, digits = digits),
168: formatC(obj$n, digits = digits),
1.1 hako 169: formatC(obj$Extinction.probability, digits = digits)),
1.8 hako 170: c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",
171: formatC(obj$upper.CL.mu, digits = digits),")", sep = ""),
172: paste("(",formatC(obj$lower.CL.s, digits = digits),", ",
173: formatC(obj$upper.CL.s, digits = digits),")", sep = ""),
174: "-",
175: "-",
176: paste("(",formatC(obj$lower.CL.P, digits = digits),", ",
177: formatC(obj$upper.CL.P, digits = digits),")", sep = "")))
1.1 hako 178: dimnames(output.est) <- list(
179: c("Growth rate:",
180: "Variance:",
1.8 hako 181: "Current population size, xd = log(N(q) / ne):",
182: "Sample size, n = q + 1:",
1.1 hako 183: paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
184: c("Estimate","95% confidence interval"))
185: print(output.est)
186: } else {
187: output.est <- data.frame(
188: c(formatC(obj$Extinction.probability, digits = digits)))
189: dimnames(output.est) <- list(
190: c(paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
191: c("Estimate"))
192: print(output.est)
193: }
194: }
195: #
196: # Examples
197: # Yellowstone grizzly bears (from Dennis et al., 1991)
1.4 hako 198: # 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),
1.1 hako 199: # 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))
200: # The probability of extinction (of decline to population size 1) within 100 years
1.8 hako 201: # ext1(dat, t = 100)
1.1 hako 202: # The probability of decline to 10 individuals within 100 years
1.8 hako 203: # ext1(dat, t = 100, ne = 10, verbose = 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