Return to extinction_risk_2.R CVS log

Up to [local] / ext

add Id

# extinction_risk_2.R, ver. 2.5 2014/2/12
# $Id: extinction_risk_2.R,v 1.2 2015/05/25 12:51:56 hako Exp $
#
# Author: Hiroshi Hakoyama <hako@affrc.go.jp>
# Copyright (c) 2013-2014 Hiroshi Hakoyama <hako@affrc.go.jp>, All rights reserved.
# 
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY 
# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 
# AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
# NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 
# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# Description:
# Estimates the demographic parameters and the probability of extinction within 
# a specific time period, t, from a time series of population size based on a 
# density-dependent population model with environmental and demographic
# stochasticity (Hakoyama and Iwasa, 2000).
#
# Usage:
# ext2(dat, gt = 3, t = 100, INDEX = FALSE, PS = 10^6, verbose = FALSE, accuracy = 0.25)
#
# Arguments:
# dat		data.frame of 2 variables: Time and Population size
# gt		generation time
# t			a time period of interest
# 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 per generation (r),
# the ML estimate of carrying capacity (K),
# the ML estimate of variance per generation (s),
# the ML estimate of the mean extinction time (met, year)
# the ML estimate of the probability of extinction per year (lambda)
# the ML estimate of the probability of extinction within a specific time period t (ext, year)
# INDEX		If data are an index of the population size, set to TRUE. If data are
# population size, set to FALSE. If INDEX is set to TRUE, give only r and ss, 
# and the probability of extinction time within a specific time period t based on 
# Te(r, PS, ss).
# PS		Population size from other information. Set PS for the data of 
# an index of the population size.
# accuracy	accuracy for the integrate function
#
# References:
#  H. Hakoyama and Y. Iwasa. Extinction risk of a density-dependent population
# estimated from a time series of population size. Journal of Theoretical Biology,
# 204:337-359, 2000.
#
ext2 <- function(dat,
                 gt = 3, t = 100, INDEX = FALSE, 
                 PS = 10^6, verbose = FALSE, accuracy = 0.25) {
  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)
  yr <- yr[complete]
  ps <- ps[complete]
  tau <- diff(yr) / gt # time intervals, \tau_i = (t_i - t_{i-1})/gt
  K <- mean(ps) # Carrying capacity
  y <- ps - K
  n <- length(ps)
  if (sum(y[1:length(y) - 1] * y[2:length(y)]) < 0) {
    print("The ML estimates do not exist, because the autocorrelation with delay tau is negative (see Hakoyama & Iwasa, 2000, p-357).")
    return()
  }
  if (all(tau == tau[1])) {
    f <- function(b) {
      (n - 1) / n * y[1]^2 * b * (1 - b^2) - b / n * sum((y[2:length(y)] - b *
      y[1:length(y) - 1])^2) + (1 - b^2) * sum(y[1:length(y) - 1] *
      (y[2:length(y)] - b * y[1:length(y) - 1]))
    }
    beta <- uniroot(f, c(0, 1), tol = 10^-9)$root
    alpha <- (y[1]^2 + sum((y[2:length(y)] -
    beta * y[1:length(y) - 1])^2) / (1 - beta^2)) / n
    r <- - log(beta) / tau[1] # growth rate per generation
  } else {
  	f <- function(b) {
  	  sum(- tau * b^(2 * tau - 1) / (1 - b^(2 * tau))) / n *
  	  (y[1]^2 + sum((y[2:length(y)] - b^tau * y[1:length(y) - 1])^2 / 
  	  (1 - b^(2 * tau)))) +
  	  sum(b^(tau - 1) * tau * (y[2:length(y)] - b^tau * y[1:length(y) - 1]) * 
  	  (b^tau * y[2:length(y)] - y[1:length(y) - 1]) / (1 - b^(2 * tau))^2)
    }
    beta <- uniroot(f, c(10^-300, 1 - 10^-9), tol = 10^-50)$root
    alpha <- (y[1]^2 + sum((y[2:length(y)] -
    beta^(tau) * y[1:length(y) - 1])^2 / (1 - beta^(2 * tau)))) / n 
    r <- - log(beta) # growth rate per generation
  }

  if (INDEX == FALSE) {
    s <- (2 * alpha * r - K) / K^2 # environmental variance
    if (s < 0) {
      print("Negative variance was estimated!")
      return()
    }
    met <- Te(r, K, s, a = accuracy) * gt # Mean extinction time (year)
    lambda <- 1 / met # Extinction probability per year
    ext <- ifelse(1 - exp(-t * lambda) > 0, 1 - exp(-t * lambda), t * lambda) # Extinction probability within t years based on exponential distribution or approximation for small lambda
  } else {
    ss <- (2 * alpha * r) / K^2 # s for INDEX data with large population size
    met2 <- Te(r, PS, ss, a = accuracy) * gt
    lambda2 <- 1 / met2
    ext2 <- ifelse(1 - exp(- t * lambda2) > 0, 1 - exp(-t * lambda2), t * lambda2)
  } # for INDEX data
  
  if (INDEX == FALSE) {
    if (verbose == TRUE) {
      results <- list(t = t, verbose = verbose, INDEX = INDEX, 
                      Growth.rate.per.generation = r,
                      Carrying.capacity = K,
                      Variance.per.generation = s, 
                      Mean.extinction.time = met, 
                      Extinction.probability.per.year = lambda,
                      Extinction.probability.within.t.years = ext)
      class(results) <- "ext2"
      return(results)
    } else {
      results <- list(t = t, verbose = verbose, INDEX = INDEX,
                    Extinction.probability.within.t.years = ext)
      class(results) <- "ext2"
      return(results)    
    }
  } else {
    if (verbose == TRUE) {
      results <- list(t = t, verbose = verbose, INDEX = INDEX, 
                      Growth.rate.per.generation = r, 
                      Variance.per.generation = ss,
                      Extinction.probability.within.t.years = ext2)
      class(results) <- "ext2"
      return(results)     
    } else {
      results <- list(t = t, verbose = verbose, INDEX = INDEX, 
                      Extinction.probability.within.t.years = ext2)
      class(results) <- "ext2"
      return(results)
    }
  }
}

print.ext2 <- function(obj, digits = 5) {
#  cat("\nExtinction risk\n")
  if ((obj$INDEX == FALSE) && (obj$verbose == TRUE)) {
    output.est <- data.frame(c(
      formatC(obj$Growth.rate.per.generation, digits = digits), 
      formatC(obj$Carrying.capacity, digits = digits), 
      formatC(obj$Variance.per.generation, digits = digits),
      formatC(obj$Mean.extinction.time, digits = digits),
      formatC(obj$Extinction.probability.per.year, digits = digits),
      formatC(obj$Extinction.probability.within.t.years, digits = digits)))
    dimnames(output.est) <- list(
      c("Growth rate per generation:",
        "Carrying capacity:",
        "Variance per generation:",
        "Mean extinction time (year):",
        "Extinction probability per year:",
        paste("Extinction probability within", obj$t, "years:")),
      c("Estimate"))
    print(output.est)
  }

  if ((obj$INDEX == TRUE) && (obj$verbose == TRUE)) {
    output.est <- data.frame(c(
      formatC(obj$Growth.rate.per.generation, digits = digits), 
      formatC(obj$Variance.per.generation, digits = digits),
      formatC(obj$Extinction.probability.within.t.years, digits = digits)))
    dimnames(output.est) <- list(
      c("Growth rate per generation:",
        "Variance per generation:",
        paste("Extinction probability within", obj$t, "years:")),
      c("Estimate"))
    print(output.est)
  }

  if (obj$verbose == FALSE) {
    output.est <- data.frame(
      c(formatC(obj$Extinction.probability, digits = digits)))
    dimnames(output.est) <- list(
      c(paste("Extinction probability within", obj$t, "years:")),
      c("Estimate"))
  	print(output.est)
  }
}

Te <- function(r, K, s2, x0 = K, a = 0.25, limit = 10) {
  # Mean extinction time, Te (generation time)
  R <- 2 * r / (s2 * K)
  DD <- 1 / s2
  f1 <- function(u, v) {
  	exp(- R * x0 * (1 / v - 0.5 * (tanh(0.5 * pi * sinh(u)) + 1)) + (R * (K + DD) + 1) * log(((x0 / v + DD) / (0.5 * x0 * (tanh(0.5 * pi * sinh(u)) + 1) + DD)))) * (0.5 * (pi * cosh(u)) / (1 + cosh(pi * sinh(u)))) / (1 + v * DD / x0)
  }
  T1 <- function() {
  	2 * DD * integrate(function(x) {
      sapply(x, function(x) {
        integrate(function(y) f1(x, y), 0, 1,
        stop.on.error = T, rel.tol = .Machine$double.eps^a)$value
      })
    }, -limit, limit, stop.on.error = T, rel.tol = .Machine$double.eps^a)$value
  }
  f2 <- function(u, v) {
  	exp(- R * (x0 * v - 0.5 * x0 * (tanh(0.5 * pi * sinh(u)) + 1)) + (R * (K + DD) + 1) * log(((x0 * v + DD) / (0.5 * x0 * (tanh(0.5 * pi * sinh(u)) + 1) + DD)))) * (0.5 * (pi * cosh(u)) / (1 + cosh(pi * sinh(u)))) / (v^2 + v * DD / x0)
  }
  T2 <- function() {
  	2 * DD * integrate(function(x) {
      sapply(x, function(x) {
        integrate(function(y) f2(x, y), 0.5 * (tanh(0.5 * pi * sinh(x)) + 1), 1,
        stop.on.error = T, rel.tol = .Machine$double.eps^a)$value
      })
    }, -limit, limit, stop.on.error = T, rel.tol = .Machine$double.eps^a)$value
  }
T1() + T2()
}

#
# 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))
# the probability of extinction (to decline to population size 0) within 100 years
# ext2(dat, gt = 10, t = 100)
# ext2(dat, gt = 10, t = 100, verbose = TRUE)
# CPUE data of the Japanese crucian carp in Lake Biwa (from Hakoyama and Iwasa, 2000)
# dat2 <- data.frame(Year = c(1955, 1956, 1957, 1958, 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, 1988, 1989), 
# CPUE = c(0.230711329, 0.287504747, 0.321044547, 0.248123271, 0.266040689, 0.276898734, 0.360330579, 0.365279529, 0.405311276, 0.415012942, 0.466610313, 0.373273942, 0.349548646, 0.245173745, 0.31368529, 0.320981211, 0.236671001, 0.263181412, 0.263037511, 0.346241458, 0.290079925, 0.250327654, 0.284950658, 0.253397633, 0.303960837, 0.35359857, 0.399908592, 0.320795504, 0.237847222, 0.260603205, 0.291603821, 0.301130524, 0.272430669, 0.221655329, 0.186635945))
# ext2(dat2, gt = 3, t = 100, INDEX = TRUE, PS = 10^6)