=================================================================== RCS file: /cvs/ext/extinction_risk_1.R,v retrieving revision 1.19 retrieving revision 1.20 diff -u -p -r1.19 -r1.20 --- ext/extinction_risk_1.R 2015/08/14 11:39:52 1.19 +++ ext/extinction_risk_1.R 2015/08/15 06:50:20 1.20 @@ -1,5 +1,5 @@ # extinction_risk_1.R -# $Id: extinction_risk_1.R,v 1.19 2015/08/14 11:39:52 hako Exp $ +# $Id: extinction_risk_1.R,v 1.20 2015/08/15 06:50:20 hako Exp $ # # Author: Hiroshi Hakoyama # Copyright (c) 2013-2015 Hiroshi Hakoyama , All rights reserved. @@ -35,16 +35,16 @@ # which is better than the estimator from the ordinary Delta method. # # Usage: -# ext1(dat, t = 100, ne = 1, alpha = 0.05, verbose = FALSE, +# ext1(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) # # 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 +# tt a time period of interest +# ne a lower (extinction) threshold of population size (ne >= 1) # alpha 1 - confidence level # verbose If set to FALSE, give the ML estimate of the probability -# of extinction within a specific time period t. If set to +# of extinction within a specific time period tt. If set to # TRUE, give a more verbose output: # the ML estimate of the growth rate (mu), # the ML estimate of variance (s), @@ -74,7 +74,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, qq.plot = FALSE, formatted = TRUE) { +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) { 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) @@ -82,118 +82,42 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v ps <- ps[complete] tau <- diff(yr) # time intervals, \tau_i = t_i - t_{i-1} delta.log.n <- diff(log(ps)) - q <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\} + qq <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\} tq <- sum(tau) - mu <- sum(delta.log.n) / tq # ML estimate of growth rate - s <- (1 / q) * sum((delta.log.n - mu * tau)^2 / tau) # ML estimate of variance - us <- q * s / (q - 1) # an unbiased estimate of variance + mu <- sum(delta.log.n) / tq # MLE of growth rate + s <- sum((delta.log.n - mu * tau)^2 / tau) / qq # MLE of variance + us <- qq * s / (qq - 1) # an unbiased estimate of variance xd <- log(ps[length(ps)] / ne) - l <- - (q / 2) * log(2 * pi * s) - sum((delta.log.n - mu)^2) / (2 * s) # ML for the distribution of x = log.n (not for n). - AIC <- - 2 * l + 4 # AIC for x + l <- - log(2 * pi * s) * qq / 2 - sum((delta.log.n - mu)^2) / (2 * s) + AIC.X <- - 2 * l + 4 # AIC for the distribution of X = log(N) 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) - upper.CL.s <- q * s / qchisq(alpha / 2, q - 1) + lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq) + upper.CL.mu <- mu + qt(1 - alpha / 2, qq - 1) * sqrt(us / tq) + lower.CL.s <- qq * s / qchisq(1 - alpha / 2, qq - 1) + upper.CL.s <- qq * s / qchisq(alpha / 2, qq - 1) - w <- function(mu, xd, s, t) (mu * t + xd) / sqrt(s * t) - z <- function(mu, xd, s, t) (- mu * t + xd) / sqrt(s * t) - pr <- function(w, z) { - 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 - + 105 / z^9 - 945 / z^11 + 10395 / z^13) - } - } + ww <- W(mu, xd, s, tt) + zz <- Z(mu, xd, s, tt) + pp <- Pr(ww, zz) - ww <- w(mu, xd, s, t) - zz <- z(mu, xd, s, t) - pp <- pr(ww, zz) - - c.limit <- function(tq, q, t, est, a, width = 10) { - t.obs <- sqrt((q - 1) * tq / (q * t)) * est - df <- q - 1 - const <- sqrt(tq / t) - f <- function(x) { - pt(t.obs, df, const * x) - a - } - d.est <- est / width + 1 - uniroot(f, c(- d.est + est, d.est + est), extendInt = "yes")$root - } - - confidence.interval <- function(mu, xd, s, t, tq, q, alpha) { - den1 <- sqrt(s * t) - w.est <- (mu * t + xd) / den1 - z.est <- (- mu * t + xd) / den1 - lower.CL.w <- c.limit(tq, q, t, w.est, 1 - alpha / 2) - upper.CL.w <- c.limit(tq, q, t, w.est, alpha / 2) - lower.CL.z <- c.limit(tq, q, t, z.est, 1 - alpha / 2) - upper.CL.z <- c.limit(tq, q, t, z.est, alpha / 2) - lower.CL.P <- pr(upper.CL.w, lower.CL.z) - upper.CL.P <- pr(lower.CL.w, upper.CL.z) - c(lower.CL.P, upper.CL.P) - } - - CI.rand <- function(mu.obs, xd, s.obs, t, tq, q, alpha) { - mu <- rnorm(1, mean = mu.obs, sd = sqrt(s.obs / q)) - s <- rchisq(1, df = q - 1) * s.obs / q - confidence.interval(mu, xd, s, t, tq, q, alpha) - } - - CI.sim <- function(n.sim, gam, alpha, mu.obs, xd, s.obs, t, tq, q) { - CI.dist <- replicate(n.sim, CI.rand(mu.obs, xd, s.obs, t, tq, q, alpha)) - P.obs <- pr(w(mu.obs, xd, s.obs, t), z(mu.obs, xd, s.obs, t)) - ci <- confidence.interval(mu.obs, xd, s.obs, t, tq, q, alpha) - complete <- complete.cases(CI.dist[1, ], CI.dist[2, ]) - lower <- CI.dist[1, ][complete] - upper <- CI.dist[2, ][complete] - n.estimables <- length(lower) - n.rejects <- (sum(P.obs < lower) + sum(P.obs > upper)) - if (n.estimables > 0) { - binom <- binom.test(n.rejects, n.estimables, p = gam) + if (exact.CL == TRUE) { + CL.P <- AsymptoticallyExactConfidenceInterval(mu, xd, s, tt, tq, qq, alpha, n.sim, alpha.sim) } else { - binom <- list(estimate = NaN, p.value = NaN) - } - c(binom$estimate[[1]], binom$p.value, ci) + CL.P <- ConfidenceInterval(mu, xd, s, tt, tq, qq, alpha) } - asymptotically.exact.confidence.interval <- function(n.sim, alpha.sim, mu.obs, xd, s.obs, t, tq, q, alpha) { - P.obs <- pr(w(mu.obs, xd, s.obs, t), z(mu.obs, xd, s.obs, t)) - if (P.obs == 0 || P.obs == 1) { - return(confidence.interval(mu.obs, xd, s.obs, t, tq, q, alpha)) - } - res <- CI.sim(n.sim, alpha, alpha, mu.obs, xd, s.obs, t, tq, q) - gg <- res[[1]] - pp <- res[[2]] - a <- alpha - while (pp < alpha.sim) { - res <- CI.sim(n.sim, alpha, a, mu.obs, xd, s.obs, t, tq, q) - gg <- res[[1]] - pp <- res[[2]] - a <- a * alpha / gg - } - c(res[[3]], res[[4]]) - } - - 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) - } lower.CL.P <- CL.P[[1]] upper.CL.P <- CL.P[[2]] if (verbose == TRUE) { - results <- list(ne = ne, t = t, verbose = verbose, - AIC = AIC, - n = q + 1, + results <- list(ne = ne, tt = tt, verbose = verbose, + AIC.X = AIC.X, + n = qq + 1, xd = xd, Growth.rate = mu, lower.CL.mu = lower.CL.mu, @@ -210,7 +134,7 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v } return(results) } else { - results <- list(ne = ne, t = t, verbose = verbose, + results <- list(ne = ne, tt = tt, verbose = verbose, Extinction.probability = pp) if (formatted == TRUE) { class(results) <- "ext1" @@ -219,6 +143,90 @@ ext1 <- function(dat, t = 100, ne = 1, alpha = 0.05, v } } +#------------------------------------------------------------------- +W <- function(mu, xd, s, tt) (mu * tt + xd) / sqrt(s * tt) +#------------------------------------------------------------------- +Z <- function(mu, xd, s, tt) (- mu * tt + xd) / sqrt(s * tt) +#------------------------------------------------------------------- +Pr <- function(w, z) { + 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 + + 105 / z^9 - 945 / z^11 + 10395 / z^13) + } +} +#------------------------------------------------------------------- +FindCL <- function(tq, qq, tt, est, a, width = 10) { + t.obs <- sqrt((qq - 1) * tq / (qq * tt)) * est + df <- qq - 1 + const <- sqrt(tq / tt) + f <- function(x) pt(t.obs, df, const * x) - a + d.est <- est / width + 1 + uniroot(f, c(- d.est + est, d.est + est), extendInt = "yes")$root +} +#------------------------------------------------------------------- +ConfidenceInterval <- function(mu, xd, s, tt, tq, qq, alpha) { + den1 <- sqrt(s * tt) + w.est <- (mu * tt + xd) / den1 + z.est <- (- mu * tt + xd) / den1 + lower.CL.w <- FindCL(tq, qq, tt, w.est, 1 - alpha / 2) + upper.CL.w <- FindCL(tq, qq, tt, w.est, alpha / 2) + lower.CL.z <- FindCL(tq, qq, tt, z.est, 1 - alpha / 2) + upper.CL.z <- FindCL(tq, qq, tt, z.est, alpha / 2) + lower.CL.P <- Pr(upper.CL.w, lower.CL.z) + upper.CL.P <- Pr(lower.CL.w, upper.CL.z) + c(lower.CL.P, upper.CL.P) +} +#------------------------------------------------------------------- +randomCI <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function = ConfidenceInterval) { + mu <- rnorm(1, mean = mu.obs, sd = sqrt(s.obs / qq)) + s <- rchisq(1, df = qq - 1) * s.obs / qq + ci.function(mu, xd, s, tt, tq, qq, alpha) +} +#------------------------------------------------------------------- +RejectionRate <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, gam, n.sim, ci.function = ConfidenceInterval) { + CI.dist <- replicate(n.sim, randomCI(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function)) + P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt)) + ci <- ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha) + complete <- complete.cases(CI.dist[1, ], CI.dist[2, ]) + lower <- CI.dist[1, ][complete] + upper <- CI.dist[2, ][complete] + n.estimables <- length(lower) + n.rejects <- (sum(P.obs < lower) + sum(P.obs > upper)) + if (n.estimables > 0) { + binom <- binom.test(n.rejects, n.estimables, p = gam) + } else { + binom <- list(estimate = NaN, p.value = NaN) + } + c(binom$estimate[[1]], binom$p.value, ci) +} +#------------------------------------------------------------------- +AsymptoticallyExactConfidenceInterval <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, n.sim, alpha.sim, ci.function = ConfidenceInterval) { + P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt)) + if (P.obs == 0 || P.obs == 1) { + return(c(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha), NaN, NaN)) + } + res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function) + if (is.nan(res[[1]])) return(c(NaN, NaN, alpha, 0)) + gg <- res[[1]] + pp <- res[[2]] + ai <- alpha + i <- 0 + while (pp < alpha.sim) { + i <- i + 1 + res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, ai, alpha, n.sim, ci.function) + if (is.nan(res[[1]])) return(c(NaN, NaN, ai, i)) + gg <- res[[1]] + pp <- res[[2]] + ai <- ai * alpha / gg + if (ai < 10^(-10) || ai > 1) return(c(NaN, NaN, ai, i)) + if (i > 20) return(c(NaN, NaN, ai, i)) + } + c(res[[3]], res[[4]], ai, i) +} +#------------------------------------------------------------------- print.ext1 <- function(obj, digits = 5) { if (obj$verbose == TRUE) { output.est <- data.frame( @@ -226,7 +234,7 @@ print.ext1 <- function(obj, digits = 5) { formatC(obj$Variance, digits = digits), formatC(obj$xd, digits = digits), formatC(obj$n, digits = digits), - formatC(obj$AIC, digits = digits), + formatC(obj$AIC.X, digits = digits), formatC(obj$Extinction.probability, digits = digits)), c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ", formatC(obj$upper.CL.mu, digits = digits),")", sep = ""), @@ -240,17 +248,17 @@ print.ext1 <- function(obj, digits = 5) { dimnames(output.est) <- list( c("Growth rate:", "Variance:", - "Current population size, xd = log(N(q) / ne):", + "Current population size, x_d = log(N(q) / n_e):", "Sample size, n = q + 1:", "AIC for the distribution of X = log(N):", - paste("Probability of decline to", obj$ne, "within", obj$t, "years:")), + paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")), c("Estimate","95% confidence interval")) print(output.est) } else { output.est <- data.frame( c(formatC(obj$Extinction.probability, digits = digits))) dimnames(output.est) <- list( - c(paste("Probability of decline to", obj$ne, "within", obj$t, "years:")), + c(paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")), c("Estimate")) print(output.est) }