ext/extinction_risk_1.R - diff

Return to extinction_risk_1.R CVS log

Up to [local] / ext

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

version 1.22, 2015/08/23 11:43:24

version 1.23, 2017/12/20 11:37:39

Line 2

# $Id$

# Author: Hiroshi Hakoyama <hako@affrc.go.jp>

# Redistribution and use in source and binary forms, with or without

# modification, are permitted provided that the following conditions

Line 32

# of population size based on the Wiener-drift process model.

# An estimator for confidence interval of extinction risk

# developed by Hakoyama (w-z method, in preparation) is implemented,

# which is better than the estimator from the ordinary Delta method.

# which is better than the estimator by the ordinary Delta method.

# Usage:

# 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)

# better.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

Line 50

# the ML estimate of variance (s),

# Current population size, xd = log(N(q) / ne),

# Sample size, n = q + 1,

# AIC for the distribution of X = log(N),

# AIC for the distribution of N,

# the ML estimate of the probability of extinction within

# a specific time period t (P),

# a lower alpha % confidence limit of parameter *

# a lower (1 - alpha) % confidence limit of parameter *

# (lower.CL.*), and

# an upper alpha % confidence limit of parameter *

# an upper (1 - alpha) % confidence limit of parameter *

# (upper.CL.*).

# exact.CL If set to FALSE, give a CI by the w-z method.

# better.CL If set to FALSE, give a CI by the w-z method.

# If set to TRUE, give an asymptotically exact CI.

# If set to TRUE, give a better CI.

# n.sim The number of iteration for the asymptotically exact CI.

# n.sim The number of iteration for the better CI method.

# 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.

# qq.plot If set to TRUE, give a QQ plot for log(n_i/n_{i-1})/sqrt(tau_i) ~ N(0, 1).

# formatted If set to TRUE, give the result by formatted output.

# If set to FALSE, give a list of estimates.

# References:

Line 71

# B. Dennis, P. L. Munholland, and J. M. Scott. Estimation of

# growth and extinction parameters for endangered species.

# Ecological Monographs, 61:115-143, 1991.

# H. Hakoyama (in preparation).

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) {

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 <- ts(dat[, 1], start = c(dat[, 1][1])) # Year

yr <- dat[, 1] # Year

ps <- ts(dat[, 2], start = c(dat[, 1][1])) # Population size

ps <- dat[, 2] # Population size

complete <- complete.cases(yr, ps)

yr <- yr[complete]

ps <- ps[complete]

tau <- diff(yr) # time intervals, \tau_i = t_i - t_{i-1}

ti <- yr - yr[1] # elapsed time, ti = \{t_0 = 0, t_1, \dots, t_q\}

delta.log.n <- diff(log(ps))

tau <- diff(ti) # time intervals, \tau_i = t_i - t_{i-1}

qq <- length(yr) - 1 # yr = \{t_0, t_1, \dots, t_q\}

delta.log.ps <- diff(log(ps))

tq <- sum(tau)

qq <- length(yr) - 1

mu <- sum(delta.log.n) / tq # MLE of growth rate

tq <- ti[length(ti)]

s <- sum((delta.log.n - mu * tau)^2 / tau) / qq # MLE of variance

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

us <- qq * s / (qq - 1) # an unbiased estimate of variance

xd <- log(ps[length(ps)] / ne)

l <- - log(2 * pi * s) * qq / 2 - sum((delta.log.n - mu)^2) / (2 * s)

lnl <- - sum(log(ps[-1] * sqrt(2 * tau * pi))) - (qq / 2) * log(s) - sum((1 / tau) * (delta.log.ps - mu * tau)^2) / (2 * s)

AIC.X <- - 2 * l + 4 # AIC for the distribution of X = log(N)

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.n)

qqnorm(delta.log.ps / tau)

qqline(delta.log.n)

qqline(delta.log.ps / tau)

}

lower.CL.mu <- mu - qt(1 - alpha / 2, qq - 1) * sqrt(us / tq)

Line 105 ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05,

Line 107 ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05,

zz <- Z(mu, xd, s, tt)

pp <- Pr(ww, zz)

if (exact.CL == TRUE) {

if (better.CL == TRUE) {

CL.P <- AsymptoticallyExactConfidenceInterval(mu, xd, s, tt, tq, qq, alpha, n.sim, alpha.sim)

CL.P <- BetterConfidenceInterval(mu, xd, s, tt, tq, qq, alpha, n.sim, alpha.sim)

} else {

CL.P <- ConfidenceInterval(mu, xd, s, tt, tq, qq, alpha)

}

Line 116 ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05,

Line 118 ext1 <- function(dat, tt = 100, ne = 1, alpha = 0.05,

if (verbose == TRUE) {

results <- list(ne = ne, tt = tt, verbose = verbose,

AIC.X = AIC.X,

alpha = alpha,

AIC = AIC,

n = qq + 1,

xd = xd,

Growth.rate = mu,

Line 159 Pr <- function(w, z) {

Line 162 Pr <- function(w, z) {

}

#-------------------------------------------------------------------

ConfidenceInterval <- function(mu, xd, s, tt, tq, qq, alpha) {

den1 <- sqrt(s * tt)

den.1 <- sqrt(s * tt)

w.est <- (mu * tt + xd) / den1

w.est <- (mu * tt + xd) / den.1

z.est <- (- mu * tt + xd) / den1

z.est <- (- mu * tt + xd) / den.1

df <- qq - 1

const.1 <- sqrt((qq - 1) * tq / (qq * tt))

const.2 <- sqrt(tq / tt)

FindCL <- function(tq, qq, tt, est, a, width = 10) {

t.obs <- const.1 * est

f <- function(x) pt(t.obs, df, const.2 * x) - a

d.est <- est / width + 1

uniroot(f, c(- d.est + est, d.est + est), extendInt = "yes")$root

Line 203 RejectionRate <- function(mu.obs, xd, s.obs, tt, tq, q

Line 206 RejectionRate <- function(mu.obs, xd, s.obs, tt, tq, q

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, verbose = FALSE) {

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))

complete <- complete.cases(CI.dist[1, ])

CL <- CI.dist[1, ][complete]

n.estimables <- length(CL)

n.rejects <- sum(P.obs < CL)

if (n.estimables > 0) {

binom <- binom.test(n.rejects, n.estimables, p = gam / 2)

} else {

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) {

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))

complete <- complete.cases(CI.dist[2, ])

CL <- CI.dist[2, ][complete]

n.estimables <- length(CL)

n.rejects <- sum(P.obs > CL)

if (n.estimables > 0) {

binom <- binom.test(n.rejects, n.estimables, p = gam / 2)

} else {

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) {

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(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)

gg <- res[[1]]

if (res[[2]] > alpha.sim) {

pp <- res[[2]]

cat("No need of using the better CI method.", "\n")

if (verbose == TRUE) cat(alpha, gg, "\n")

return(ci.function(mu.obs, xd, s.obs, tt, tq, qq, alpha))

ai <- alpha^2 / gg

}

p.value <- 0

gamma.lower <- alpha

i <- 0

while (pp < alpha.sim) {

while (i < min.iteration || p.value < alpha.sim) {

i <- i + 1

res <- RejectionRate(mu.obs, xd, s.obs, tt, tq, qq, ai, alpha, n.sim, ci.function)

res <- RejectionRate.lower(mu.obs, xd, s.obs, tt, tq, qq, gamma.lower, alpha, n.sim, ci.function)

gg <- res[[1]]

ff <- res[[1]]

pp <- res[[2]]

p.value <- res[[2]]

if (verbose == TRUE) cat(ai, gg, "\n")

if (verbose == TRUE) cat("gamma =", gamma.lower, "f_hat(gamma) =", ff, "p-value =", p.value, "\n")

ai <- ai * alpha / gg

gamma.lower <- gamma.lower - (ff - alpha / 2) / i

}

if (verbose == TRUE) {

c(res[[3]], res[[4]], ai, i)

p.value <- 0

} else {

gamma.upper <- alpha

c(res[[3]], res[[4]])

i <- 0

while (i < min.iteration || p.value < 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

}

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]])

}

#-------------------------------------------------------------------

print.ext1 <- function(obj, digits = 5) {

Line 236 print.ext1 <- function(obj, digits = 5) {

Line 284 print.ext1 <- function(obj, digits = 5) {

formatC(obj$Variance, digits = digits),

formatC(obj$xd, digits = digits),

formatC(obj$n, digits = digits),

formatC(obj$AIC.X, digits = digits),

formatC(obj$AIC, 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 = ""),

Line 252 print.ext1 <- function(obj, digits = 5) {

Line 300 print.ext1 <- function(obj, digits = 5) {

"Variance:",

"Current population size, x_d = log(N(q) / n_e):",

"Sample size, n = q + 1:",

"AIC for the distribution of X = log(N):",

"AIC for the distribution of N:",

paste("Probability of decline to", obj$ne, "within", obj$tt, "years:")),

c("Estimate","95% confidence interval"))

c("Estimate", paste(as.character(100 * (1 - obj$alpha)), "% CI")))

print(output.est)

} else {

output.est <- data.frame(

Line 276 print.ext1 <- function(obj, digits = 5) {

Line 324 print.ext1 <- function(obj, digits = 5) {

# 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)

# ext1(dat, t = 100, ne = 1, verbose = TRUE, better.CL = TRUE)