=================================================================== RCS file: /cvs/ext/extinction_risk_1.R,v retrieving revision 1.23 retrieving revision 1.24 diff -u -p -r1.23 -r1.24 --- ext/extinction_risk_1.R 2017/12/20 11:37:39 1.23 +++ ext/extinction_risk_1.R 2017/12/23 00:18:12 1.24 @@ -1,5 +1,5 @@ # extinction_risk_1.R -# $Id: extinction_risk_1.R,v 1.23 2017/12/20 11:37:39 hako Exp $ +# $Id: extinction_risk_1.R,v 1.24 2017/12/23 00:18:12 hako Exp $ # # Author: Hiroshi Hakoyama # Copyright (c) 2013-2017 Hiroshi Hakoyama , All rights reserved. @@ -76,26 +76,27 @@ ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05, verbose = FALSE, better.CL = FALSE, n.sim = 10000, alpha.sim = 0.05, qq.plot = FALSE, formatted = TRUE) { yr <- dat[, 1] # Year - ps <- dat[, 2] # Population size - complete <- complete.cases(yr, ps) + n <- dat[, 2] # Population size + complete <- complete.cases(yr, n) yr <- yr[complete] - ps <- ps[complete] + n <- n[complete] ti <- yr - yr[1] # elapsed time, ti = \{t_0 = 0, t_1, \dots, t_q\} tau <- diff(ti) # time intervals, \tau_i = t_i - t_{i-1} - delta.log.ps <- diff(log(ps)) + delta.log.n <- diff(log(n)) qq <- length(yr) - 1 tq <- ti[length(ti)] - mu <- log(ps[length(ps)] / ps[1]) / tq # MLE of growth rate - s <- sum((delta.log.ps - mu * tau)^2 / tau) / qq # MLE of variance + nq <- n[length(n)] + mu <- log(nq / n[1]) / 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) - lnl <- - sum(log(ps[-1] * sqrt(2 * tau * pi))) - (qq / 2) * log(s) - sum((1 / tau) * (delta.log.ps - mu * tau)^2) / (2 * s) + xd <- log(nq / ne) + lnl <- - sum(log(n[-1] * sqrt(2 * tau * pi))) - (qq / 2) * log(s) - sum((1 / tau) * (delta.log.n - mu * tau)^2) / (2 * s) AIC <- - 2 * lnl + 2 * 2 # AIC for the distribution of N (population size) # \log(n_i/n_{i-1})/\sqrt{\tau_i} \sim N(0, 1) if (qq.plot == TRUE) { - qqnorm(delta.log.ps / tau) - qqline(delta.log.ps / tau) + qqnorm(delta.log.n / tau) + qqline(delta.log.n / tau) } lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq) @@ -120,7 +121,8 @@ ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05, results <- list(ne = ne, tt = tt, verbose = verbose, alpha = alpha, AIC = AIC, - n = qq + 1, + sample.size = qq + 1, + nq = nq, xd = xd, Growth.rate = mu, lower.CL.mu = lower.CL.mu, @@ -189,23 +191,6 @@ randomCI <- function(mu.obs, xd, s.obs, tt, tq, qq, al 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) -} -#------------------------------------------------------------------- RejectionRate.lower <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, gam, n.sim, ci.function = ConfidenceInterval) { P.obs <- Pr(W(mu.obs, xd, s.obs, tt), Z(mu.obs, xd, s.obs, tt)) CI.dist <- replicate(n.sim, randomCI(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function)) @@ -219,7 +204,6 @@ RejectionRate.lower <- function(mu.obs, xd, s.obs, tt, binom <- list(estimate = NaN, p.value = NaN) } c(binom$estimate[[1]], binom$p.value) - # return {rejection.rate.est, p-value} } #------------------------------------------------------------------- RejectionRate.upper <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, gam, n.sim, ci.function = ConfidenceInterval) { @@ -235,7 +219,6 @@ RejectionRate.upper <- function(mu.obs, xd, s.obs, tt, binom <- list(estimate = NaN, p.value = NaN) } c(binom$estimate[[1]], binom$p.value) - # return {rejection.rate.est, p-value} } #------------------------------------------------------------------- BetterConfidenceInterval <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, n.sim, alpha.sim, ci.function = ConfidenceInterval, verbose = FALSE, min.iteration = 5) { @@ -244,34 +227,35 @@ BetterConfidenceInterval <- function(mu.obs, xd, s.obs return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha)) } - res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function) - if (res[[2]] > alpha.sim) { - cat("No need of using the better CI method.", "\n") + rr.lower <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function) + rr.upper <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function) + if (rr.lower[[2]] > alpha.sim && rr.upper[[2]] > alpha.sim) { +# cat("No need of using the better CI method.", "\n") return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha)) } - - p.value <- 0 + + p.value.lower <- 0 gamma.lower <- alpha i <- 0 - while (i < min.iteration || p.value < alpha.sim) { + while (i < min.iteration || p.value.lower < alpha.sim) { i <- i + 1 - res <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower, alpha, n.sim, ci.function) - ff <- res[[1]] - p.value <- res[[2]] - if (verbose == TRUE) cat("gamma =", gamma.lower, "f_hat(gamma) =", ff, "p-value =", p.value, "\n") - gamma.lower <- gamma.lower - (ff - alpha / 2) / i + rr.lower <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower, alpha, n.sim, ci.function) + f.lower <- rr.lower[[1]] + p.value.lower <- rr.lower[[2]] + if (verbose == TRUE) cat("gamma =", gamma.lower, "f_hat(gamma) =", f.lower, "p-value =", p.value.lower, "\n") + gamma.lower <- gamma.lower - (f.lower - alpha / 2) / i } - p.value <- 0 + p.value.upper <- 0 gamma.upper <- alpha i <- 0 - while (i < min.iteration || p.value < alpha.sim) { + while (i < min.iteration || p.value.upper < alpha.sim) { i <- i + 1 - res <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, gamma.upper, alpha, n.sim, ci.function) - ff <- res[[1]] - p.value <- res[[2]] - if (verbose == TRUE) cat("gamma =", gamma.upper, "f_hat(gamma) =", ff, "p-value =", p.value, "\n") - gamma.upper <- gamma.upper - (ff - alpha / 2) / i + rr.upper <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, gamma.upper, alpha, n.sim, ci.function) + f.upper <- rr.upper[[1]] + p.value.upper <- rr.upper[[2]] + if (verbose == TRUE) cat("gamma =", gamma.upper, "f_hat(gamma) =", f.upper, "p-value =", p.value.upper, "\n") + gamma.upper <- gamma.upper - (f.upper - alpha / 2) / i } c(ci.function(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower)[[1]], ci.function(mu.obs, xd, s.obs, tt, tq, qq, gamma.upper)[[2]]) @@ -282,8 +266,9 @@ print.ext1 <- function(obj, digits = 5) { output.est <- data.frame( c(formatC(obj$Growth.rate, digits = digits), formatC(obj$Variance, digits = digits), + formatC(obj$nq, digits = digits), formatC(obj$xd, digits = digits), - formatC(obj$n, digits = digits), + formatC(obj$sample.size, digits = digits), formatC(obj$AIC, digits = digits), formatC(obj$Extinction.probability, digits = digits)), c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ", @@ -293,13 +278,15 @@ print.ext1 <- function(obj, digits = 5) { "-", "-", "-", + "-", paste("(",formatC(obj$lower.CL.P, digits = digits),", ", formatC(obj$upper.CL.P, digits = digits),")", sep = ""))) dimnames(output.est) <- list( c("Growth rate:", "Variance:", - "Current population size, x_d = log(N(q) / n_e):", - "Sample size, n = q + 1:", + "Current population size, nq:", + "xd = ln(nq / ne):", + "Sample size, q + 1:", "AIC for the distribution of N:", paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")), c("Estimate", paste(as.character(100 * (1 - obj$alpha)), "% CI")))