Annotation of ext/extinction_risk_1.R, Revision 1.4
1.4 ! hako 1: # extinction_risk_1.R, ver. 1.6 2015/6/24
1.1 hako 2: #
3: # Author: Hiroshi Hakoyama <hako@affrc.go.jp>
4: # Copyright (c) 2013-2015 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
5: #
6: # Redistribution and use in source and binary forms, with or without
7: # modification, are permitted provided that the following conditions
8: # are met:
9: # 1. Redistributions of source code must retain the above copyright
10: # notice, this list of conditions and the following disclaimer.
11: # 2. Redistributions in binary form must reproduce the above copyright
12: # notice, this list of conditions and the following disclaimer in the
13: # documentation and/or other materials provided with the distribution.
14: #
15: # THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY
16: # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
17: # THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
18: # PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
19: # AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20: # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
21: # NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22: # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
23: # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
24: # STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25: # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
26: # ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27: #
28: # Description:
29: # Estimates the demographic parameters and the probability of extinction within
30: # a specific time period, t, from a time series of population size based on the
31: # Wiener-drift process model (Dennis et al., 1991).
32: #
33: # Usage:
34: # ext1(dat, t = 100, ne = 1, verbose = FALSE)
35: #
36: # Arguments:
37: # dat data.frame of 2 variables: Time and Population size
38: # t a time period of interest
39: # ne a lower (extinction) threshold of population size
40: # verbose If set to FALSE, give the ML estimate of the probability of extinction
41: # within a specific time period t. If set to TRUE, give a more verbose output:
42: # the ML estimate of the growth rate (mu),
43: # the ML estimate of variance (s),
44: # the unbiased estimate of variance (us),
45: # the ML estimate of the probability attaining the threshold ne (psi),
46: # the ML estimate of the conditional mean time to reach the threshold, given
47: # the threshold is reached (theta),
48: # the ML estimate of the conditional probability of extinction within a specific
49: # time period t, given the threshold is reached (G),
50: # the ML estimate of the probability of extinction within a specific time period t (G*psi), and
51: # a lower 95 % confidence limit of parameter * (Low.CL.*), and
52: # a higher 95 % confidence limit of parameter * (High.CL.*)
1.4 ! hako 53: # formatted If set to TRUE, give the result by formatted output. If set to FALSE, give a list of estimates.
1.1 hako 54: #
55: # References:
56: # R. Lande and S. H. Orzack. Extinction dynamics of age-structured populations
57: # in a fluctuating environment. Proceedings of the National Academy of
58: # Sciences, 85(19):7418-7421, 1988.
59: # B. Dennis, P. L. Munholland, and J. M. Scott. Estimation of growth and
60: # extinction parameters for endangered species. Ecological Monographs,
61: # 61:115-143, 1991.
62: #
63:
1.4 ! hako 64: ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE, formatted = TRUE) {
1.1 hako 65: yr <- ts(dat[, 1], start = c(dat[, 1][1])) # Year
66: ps <- ts(dat[, 2], start = c(dat[, 1][1])) # Population size
67: complete <- complete.cases(yr, ps)
68: yr <- yr[complete]
69: ps <- ps[complete]
70: tau <- diff(yr) # time intervals, \tau_i = t_i - t_{i-1}
71: w <- diff(log(ps)) # W_i = \log(N(t_i)/N(t_{i-1})) = X(t_i) - X(t_{i-1})
72: q <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\}
73: tq <- sum(tau)
74: mu <- sum(w) / tq # ML estimate of growth rate
75: s <- (1 / q) * sum((w - mu * tau)^2 / tau) # ML estimate of variance
76: us <- q * s / (q - 1) # an unbiased estimate of variance
77: xd <- log(ps[length(ps)] / ne)
78: psi <- function(xd, mu, s) {
79: # The probability that the process will attain the threshold
80: ifelse (mu <= 0, 1, exp(- 2 * mu * xd / s))
81: }
82: var.for.psi <- function(xd, mu, s) {
83: ifelse (mu <= 0, 0, (4 * xd^2 / s) * (2 *(q -1) * mu^2 / (q^2 * s)))
84: }
85: Low.CL.psi <- function(xd, mu, s) {
86: ifelse (mu <= 0, 1, exp(- (2 * mu * xd / s) - qnorm(1 - alpha / 2) * sqrt(var.for.psi(xd, mu, s))))
87: }
88: High.CL.psi <- function(xd, mu, s) {
89: ifelse (mu <= 0, 1, exp(- (2 * mu * xd / s) + qnorm(1 - alpha / 2) * sqrt(var.for.psi(xd, mu, s))))
90: }
91: theta <- xd / abs(mu)
92: Low.CL.mu <- mu - qt(1 - alpha / 2, q - 1) * sqrt(us / tq)
93: High.CL.mu <- mu + qt(1 - alpha / 2, q - 1) * sqrt(us / tq)
94: Low.CL.s <- q * s / qchisq(1 - alpha / 2, q - 1)
95: High.CL.s <- q * s / qchisq(alpha / 2, q - 1)
96: var.mu <- s / tq
97: var.s <- 2 * s^2 * (q - 1) / q^2
98: erfc <- function (x) 2 * pnorm(-sqrt(2) * x)
99: dGpsidmu <- - (xd / s) *
100: exp(- (2 * mu * xd) / s) *
101: erfc((- mu * t + xd) / sqrt(2 * s * t))
102: dGpsids <- (exp(- (mu * t + xd)^2 / (2 * s * t)) * t * xd) /
103: (sqrt(2 * pi) * (s * t)^(3/2)) +
104: ((mu * xd) / s^2) * exp(- (2 * mu * xd) / s) *
105: erfc((- mu * t + xd) / sqrt(2 * s * t))
106: G <- function(t, xd, mu, s) {
107: # G = Pr[T <= t]: See Eq. (16) and Appendix of Dennis et al. (1991).
108: # Note that there is a typo in Eq. (16) of Dennis et al. (1991).
109: a <- xd / sqrt(s)
110: b <- abs(mu) / sqrt(s)
111: y <- (b * t - a) / sqrt(t)
112: z <- (b * t + a) / sqrt(t)
113: d0 <- 0.2316419
114: d1 <- 0.319381530
115: d2 <- - 0.356563782
116: d3 <- 1.781477937
117: d4 <- - 1.821255978
118: d5 <- 1.330274429
119: qz <- 1 / (1 + d0 * z)
120: g <- ifelse(z >= 4, pnorm(y) + dnorm(y) *
121: (1 -
122: (1 / z^2) +
123: (1 * 3 / z^4) -
124: (1 * 3 * 5 / z^6) +
125: (1 * 3 * 5 * 7 / z^8) -
126: (1 * 3 * 5 * 7 * 9 / z^10) +
127: (1 * 3 * 5 * 7 * 9 * 11 / z^12) -
128: (1 * 3 * 5 * 7 * 9 * 11 * 13 / z^14)
129: ) / z,
130: pnorm(y) +
131: dnorm(y) * (d1 * qz + d2 * qz^2 + d3 * qz^3 + d4 * qz^4 + d5 * qz^5))
132: gg <- pnorm((-xd + abs(mu) * t) / sqrt(s * t)) + exp(2 * xd * abs(mu) / s) *
133: pnorm((-xd - abs(mu) * t) / sqrt(s * t))
134: ifelse(is.nan(gg), g, gg)
135: }
136: H <- function(t, xd, mu, s) {
137: log(G(t, xd, mu, s) * psi(xd, mu, s) / (1 - G(t, xd, mu, s) * psi(xd, mu, s)))
138: }
139: Gpsi <- G(t, xd, mu, s) * psi(xd, mu, s)
140: var.H <- var.mu * (dGpsidmu / (Gpsi * (1 - Gpsi)))^2 +
141: var.s * (dGpsids / (Gpsi * (1 - Gpsi)))^2
142: hh <- H(t, xd, mu, s)
143: Low.CL.Gpsi <- 1 / (1 + exp(- hh + qnorm(1 - alpha / 2) * sqrt(var.H)))
144: High.CL.Gpsi <- 1 / (1 + exp(- hh - qnorm(1 - alpha / 2) * sqrt(var.H)))
145: if (verbose == TRUE) {
146: results <- list(ne = ne, t = t, verbose = verbose,
147: Growth.rate = mu,
148: Low.CL.mu = Low.CL.mu,
149: High.CL.mu = High.CL.mu,
150: Variance = s,
151: Low.CL.s = Low.CL.s,
152: High.CL.s = High.CL.s,
153: Low.CL.psi = Low.CL.psi(xd, mu, s),
154: High.CL.psi = High.CL.psi(xd, mu, s),
155: Unbiased.variance = us,
156: Probability.attaining.the.threshold = psi(xd, mu, s),
157: Conditional.mean.extinction.time = theta,
158: Conditional.extinction.probability = G(t, xd, mu, s),
159: Extinction.probability = G(t, xd, mu, s) * psi(xd, mu, s),
160: Low.CL.Gpsi = Low.CL.Gpsi,
161: High.CL.Gpsi = High.CL.Gpsi)
1.4 ! hako 162: if (formatted == TRUE) {
! 163: class(results) <- "ext1"
! 164: }
1.1 hako 165: return(results)
166: } else {
167: results <- list(ne = ne, t = t, verbose = verbose,
168: Extinction.probability = G(t, xd, mu, s) * psi(xd, mu, s))
1.4 ! hako 169: if (formatted == TRUE) {
! 170: class(results) <- "ext1"
! 171: }
1.1 hako 172: return(results)
173: }
174: }
175:
176: print.ext1 <- function(obj, digits = 5) {
177: # cat("\nExtinction risk\n")
178: if (obj$verbose == TRUE) {
179: output.est <- data.frame(
180: c(formatC(obj$Growth.rate, digits = digits),
181: formatC(obj$Variance, digits = digits),
182: formatC(obj$Unbiased.variance, digits = digits),
183: formatC(obj$Probability.attaining.the.threshold, digits = digits),
184: # formatC(obj$Conditional.mean.extinction.time, digits = digits),
185: # formatC(obj$Conditional.extinction.probability, digits = digits),
186: formatC(obj$Extinction.probability, digits = digits)),
187: c(paste("(",formatC(obj$Low.CL.mu, digits = digits),", ",
188: formatC(obj$High.CL.mu, digits = digits),")", sep = ""),
189: paste("(",formatC(obj$Low.CL.s, digits = digits),", ",
190: formatC(obj$High.CL.s, digits = digits),")", sep = ""),
191: paste("(",formatC(obj$Low.CL.s, digits = digits),", ",
192: formatC(obj$High.CL.s, digits = digits),")", sep = ""),
193: paste("(",formatC(obj$Low.CL.psi, digits = digits),", ",
194: formatC(obj$High.CL.psi, digits = digits),")", sep = ""),
195: # "-",
196: paste("(",formatC(obj$Low.CL.Gpsi, digits = digits),", ",
197: formatC(obj$High.CL.Gpsi, digits = digits),")", sep = "")))
198: dimnames(output.est) <- list(
199: c("Growth rate:",
200: "Variance:",
201: "Unbiased variance:",
202: "Probability attaining the threshold:",
203: # "Conditional mean extinction time:",
204: # "Conditional extinction probability:",
205: paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
206: c("Estimate","95% confidence interval"))
207: print(output.est)
208: } else {
209: output.est <- data.frame(
210: c(formatC(obj$Extinction.probability, digits = digits)))
211: dimnames(output.est) <- list(
212: c(paste("Probability of decline to", obj$ne, "within", obj$t, "years:")),
213: c("Estimate"))
214: print(output.est)
215: }
216: }
217: #
218: # Examples
219: # Yellowstone grizzly bears (from Dennis et al., 1991)
1.4 ! hako 220: # 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 221: # 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))
222: # The probability of extinction (of decline to population size 1) within 100 years
223: # ext1(dat, 100)
224: # The probability of decline to 10 individuals within 100 years
225: # ext1(dat, 100, 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