[BACK]Return to extinction_risk_1.R CVS log [TXT][DIR] Up to [local] / ext

Diff for /ext/extinction_risk_1.R between version 1.23 and 1.24

version 1.23, 2017/12/20 11:37:39 version 1.24, 2017/12/23 00:18:12
Line 76 
Line 76 
   
 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) {  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    yr <- dat[, 1] # Year
   ps <- dat[, 2] # Population size    n  <- dat[, 2] # Population size
   complete <- complete.cases(yr, ps)    complete <- complete.cases(yr, n)
   yr <- yr[complete]    yr <- yr[complete]
   ps <- ps[complete]    n  <-  n[complete]
   ti <- yr - yr[1] # elapsed time, ti = \{t_0 = 0, t_1, \dots, t_q\}    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}    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    qq <- length(yr) - 1
   tq <- ti[length(ti)]    tq <- ti[length(ti)]
   mu <- log(ps[length(ps)] / ps[1]) / tq # MLE of growth rate    nq <- n[length(n)]
   s <- sum((delta.log.ps - mu * tau)^2 / tau) / qq # MLE of variance    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    us <- qq * s / (qq - 1) # an unbiased estimate of variance
   xd <- log(ps[length(ps)] / ne)    xd <- log(nq / 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)    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)    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)  # \log(n_i/n_{i-1})/\sqrt{\tau_i} \sim N(0, 1)
   if (qq.plot == TRUE) {    if (qq.plot == TRUE) {
     qqnorm(delta.log.ps / tau)      qqnorm(delta.log.n / tau)
     qqline(delta.log.ps / tau)      qqline(delta.log.n / tau)
   }    }
   
   lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq)    lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq)
Line 120  ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05, 
Line 121  ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05, 
     results <- list(ne = ne, tt = tt, verbose = verbose,      results <- list(ne = ne, tt = tt, verbose = verbose,
       alpha = alpha,        alpha = alpha,
       AIC = AIC,        AIC = AIC,
       n = qq + 1,        sample.size = qq + 1,
         nq = nq,
       xd = xd,        xd = xd,
       Growth.rate = mu,        Growth.rate = mu,
       lower.CL.mu = lower.CL.mu,        lower.CL.mu = lower.CL.mu,
Line 189  randomCI <- function(mu.obs, xd, s.obs, tt, tq, qq, al
Line 191  randomCI <- function(mu.obs, xd, s.obs, tt, tq, qq, al
   ci.function(mu, xd, s, tt, tq, qq, alpha)    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) {  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))    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))    CI.dist <- replicate(n.sim, randomCI(mu.obs, xd, s.obs, tt, tq, qq, alpha, ci.function))
Line 219  RejectionRate.lower <- function(mu.obs, xd, s.obs, tt,
Line 204  RejectionRate.lower <- function(mu.obs, xd, s.obs, tt,
   binom <- list(estimate = NaN, p.value = NaN)    binom <- list(estimate = NaN, p.value = NaN)
   }    }
   c(binom$estimate[[1]], binom$p.value)    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) {  RejectionRate.upper <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, gam, n.sim, ci.function = ConfidenceInterval) {
Line 235  RejectionRate.upper <- function(mu.obs, xd, s.obs, tt,
Line 219  RejectionRate.upper <- function(mu.obs, xd, s.obs, tt,
   binom <- list(estimate = NaN, p.value = NaN)    binom <- list(estimate = NaN, p.value = NaN)
   }    }
   c(binom$estimate[[1]], binom$p.value)    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) {  BetterConfidenceInterval <- function(mu.obs, xd, s.obs, tt, tq, qq, alpha, n.sim, alpha.sim, ci.function = ConfidenceInterval, verbose = FALSE, min.iteration = 5) {
Line 244  BetterConfidenceInterval <- function(mu.obs, xd, s.obs
Line 227  BetterConfidenceInterval <- function(mu.obs, xd, s.obs
     return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha))      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)    rr.lower <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function)
   if (res[[2]] > alpha.sim) {    rr.upper <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, alpha, alpha, n.sim, ci.function)
         cat("No need of using the better CI method.", "\n")    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))          return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha))
   }    }
   
   p.value <- 0    p.value.lower <- 0
   gamma.lower <- alpha    gamma.lower <- alpha
   i <- 0    i <- 0
    while (i < min.iteration || p.value < alpha.sim) {     while (i < min.iteration || p.value.lower < alpha.sim) {
         i <- i + 1          i <- i + 1
     res <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower, alpha, n.sim, ci.function)      rr.lower <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower, alpha, n.sim, ci.function)
     ff <- res[[1]]      f.lower <- rr.lower[[1]]
     p.value <- res[[2]]      p.value.lower <- rr.lower[[2]]
     if (verbose == TRUE) cat("gamma =", gamma.lower, "f_hat(gamma) =", ff, "p-value =", p.value, "\n")      if (verbose == TRUE) cat("gamma =", gamma.lower, "f_hat(gamma) =", f.lower, "p-value =", p.value.lower, "\n")
     gamma.lower <- gamma.lower - (ff - alpha / 2) / i      gamma.lower <- gamma.lower - (f.lower - alpha / 2) / i
   }    }
   
   p.value <- 0    p.value.upper <- 0
   gamma.upper <- alpha    gamma.upper <- alpha
   i <- 0    i <- 0
    while (i < min.iteration || p.value < alpha.sim) {     while (i < min.iteration || p.value.upper < alpha.sim) {
         i <- i + 1          i <- i + 1
     res <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, gamma.upper, alpha, n.sim, ci.function)      rr.upper <- RejectionRate.upper(mu.obs, xd, s.obs, tt, tq, qq, gamma.upper, alpha, n.sim, ci.function)
     ff <- res[[1]]      f.upper <- rr.upper[[1]]
     p.value <- res[[2]]      p.value.upper <- rr.upper[[2]]
     if (verbose == TRUE) cat("gamma =", gamma.upper, "f_hat(gamma) =", ff, "p-value =", p.value, "\n")      if (verbose == TRUE) cat("gamma =", gamma.upper, "f_hat(gamma) =", f.upper, "p-value =", p.value.upper, "\n")
     gamma.upper <- gamma.upper - (ff - alpha / 2) / i      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]])  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]])
Line 282  print.ext1 <- function(obj, digits = 5) {
Line 266  print.ext1 <- function(obj, digits = 5) {
     output.est <- data.frame(      output.est <- data.frame(
     c(formatC(obj$Growth.rate, digits = digits),      c(formatC(obj$Growth.rate, digits = digits),
       formatC(obj$Variance, digits = digits),        formatC(obj$Variance, digits = digits),
         formatC(obj$nq, digits = digits),
       formatC(obj$xd, 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$AIC, digits = digits),
       formatC(obj$Extinction.probability, digits = digits)),        formatC(obj$Extinction.probability, digits = digits)),
     c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",      c(paste("(",formatC(obj$lower.CL.mu, digits = digits),", ",
Line 293  print.ext1 <- function(obj, digits = 5) {
Line 278  print.ext1 <- function(obj, digits = 5) {
       "-",        "-",
       "-",        "-",
       "-",        "-",
         "-",
       paste("(",formatC(obj$lower.CL.P, digits = digits),", ",        paste("(",formatC(obj$lower.CL.P, digits = digits),", ",
       formatC(obj$upper.CL.P, digits = digits),")", sep = "")))        formatC(obj$upper.CL.P, digits = digits),")", sep = "")))
     dimnames(output.est) <- list(      dimnames(output.est) <- list(
       c("Growth rate:",        c("Growth rate:",
         "Variance:",          "Variance:",
         "Current population size, x_d = log(N(q) / n_e):",          "Current population size, nq:",
         "Sample size, n = q + 1:",          "xd = ln(nq / ne):",
           "Sample size, q + 1:",
         "AIC for the distribution of N:",          "AIC for the distribution of N:",
         paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")),          paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")),
       c("Estimate", paste(as.character(100 * (1 - obj$alpha)), "% CI")))        c("Estimate", paste(as.character(100 * (1 - obj$alpha)), "% CI")))
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.24