=================================================================== RCS file: /cvs/ext/extinction_risk_1.R,v retrieving revision 1.12 retrieving revision 1.13 diff -u -p -r1.12 -r1.13 --- ext/extinction_risk_1.R 2015/08/06 13:09:31 1.12 +++ ext/extinction_risk_1.R 2015/08/08 03:38:28 1.13 @@ -1,5 +1,5 @@ # extinction_risk_1.R -# $Id: extinction_risk_1.R,v 1.12 2015/08/06 13:09:31 hako Exp $ +# $Id: extinction_risk_1.R,v 1.13 2015/08/08 03:38:28 hako Exp $ # # Author: Hiroshi Hakoyama # Copyright (c) 2013-2015 Hiroshi Hakoyama , All rights reserved. @@ -39,11 +39,11 @@ # exact.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, formatted = TRUE) # # Arguments: -# dat data.frame of 2 variables: Time and Population size -# t a time period of interest -# ne a lower (extinction) threshold of population size +# dat data.frame of 2 variables: Time and Population size +# t a time period of interest +# ne a lower (extinction) threshold of population size # alpha 1 - confidence level -# verbose If set to FALSE, give the ML estimate of the probability +# verbose If set to FALSE, give the ML estimate of the probability # of extinction within a specific time period t. If set to # TRUE, give a more verbose output: # the ML estimate of the growth rate (mu), @@ -60,6 +60,7 @@ # If set to TRUE, give an asymptotically exact CI. # n.sim The number of iteration for the asymptotically exact CI. # alpha.sim Level of significance to decide the convergence of iteration. +# qq.plot If set to TRUE, give a QQ plot for delta.log.n. # formatted If set to TRUE, give the result by formatted output. # If set to FALSE, give a list of estimates. # References: @@ -72,7 +73,7 @@ # H. Hakoyama (in preparation). # -ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE, exact.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, formatted = TRUE) { +ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE, exact.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, qq.plot = FALSE, formatted = TRUE) { yr <- ts(dat[, 1], start = c(dat[, 1][1])) # Year ps <- ts(dat[, 2], start = c(dat[, 1][1])) # Population size complete <- complete.cases(yr, ps) @@ -87,9 +88,14 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v us <- q * s / (q - 1) # an unbiased estimate of variance xd <- log(ps[length(ps)] / ne) - lower.CL.mu <- mu - qt(1 - alpha / 2, q - 1) * sqrt(us / tq) + if (qq.plot == TRUE) { + qqnorm(delta.log.n) + qqline(delta.log.n) + } + + lower.CL.mu <- mu - qt(1 - alpha / 2, q - 1) * sqrt(us / tq) upper.CL.mu <- mu + qt(1 - alpha / 2, q - 1) * sqrt(us / tq) - lower.CL.s <- q * s / qchisq(1 - alpha / 2, q - 1) + lower.CL.s <- q * s / qchisq(1 - alpha / 2, q - 1) upper.CL.s <- q * s / qchisq(alpha / 2, q - 1) w <- function(mu, xd, s, t) (mu * t + xd) / sqrt(s * t) @@ -98,7 +104,9 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v if(z < 35) { pnorm(-w) + exp((z^2 - w^2) / 2) * pnorm(-z) } else { - 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) + pnorm(-w) + exp(- w^2 / 2) * (sqrt(2) / (2 * sqrt(pi))) * + (1 / z - 1 / z^3 + 3 / z^5 - 15 / z^7 + + 105 / z^9 - 945 / z^11 + 10395 / z^13) } } @@ -174,10 +182,10 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v } if (exact.CL == TRUE) { - CL.P <- asymptotically.exact.confidence.interval(n.sim, alpha.sim, mu, xd, s, t, tq, q, alpha) - } else { - CL.P <- confidence.interval(mu, xd, s, t, tq, q, alpha) - } + CL.P <- asymptotically.exact.confidence.interval(n.sim, alpha.sim, mu, xd, s, t, tq, q, alpha) + } else { + CL.P <- confidence.interval(mu, xd, s, t, tq, q, alpha) + } lower.CL.P <- CL.P[[1]] upper.CL.P <- CL.P[[2]] @@ -239,17 +247,19 @@ print.ext1 <- function(obj, digits = 5) { dimnames(output.est) <- list( c(paste("Probability of decline to", obj$ne, "within", obj$t, "years:")), c("Estimate")) - print(output.est) + print(output.est) } } # # Examples # Yellowstone grizzly bears (from Dennis et al., 1991) - # 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), -# 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)) + # 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), + # 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)) # The probability of extinction (of decline to population size 1) within 100 years # ext1(dat, t = 100) # The probability of decline to 10 individuals within 100 years # ext1(dat, t = 100, ne = 10, verbose = TRUE) +# with QQ-plot +# ext1(dat, t = 100, ne = 10, verbose = TRUE, qq.plot = T) # The probability of decline to 1 individuals within 100 years with an asymptotically exact confidence interval # ext1(dat, t = 100, ne = 1, verbose = TRUE, exact.CL = TRUE) \ No newline at end of file