Mortality rate with different interval lengths for the different observations

[maurolepore]

Contributed by Gabriel Arellano (via @)

# Mortality rate calculation using different interval lengths 
# for the different observations. Based on Kubo et al. 2000.
mortality_rate_Kubo <- function(data)
{
  # Approach taken by Kubo et al., using different
  # intervals for each observation, and the same
  # mortality rate for all the individuals:
  # Journal of Tropical Ecology (2000) 16:753-756.
  # The mortality rate is the value that makes
  # Equation 3 of Kubo et al. 2000 equal to 0.
  Eq3 <- function(mortality.rate, data)
  {
    alive <- which(data$status == "A")
    dead <- which(data$status == "D")
    Ti <- data[alive, "t"]
    Tj <- data[dead, "t"]
    x = -sum(Ti) + sum((Tj * exp(-mortality.rate * Tj))/(1 - exp(-mortality.rate * Tj)))
    abs(x)
  }
  
  # It is a single-parameter optimization problem.
  # The upper limit of the instantaneous mortality rate is +Inf in theory,
  # but at low mortality rates it should be almost equal to the truly annual
  # mortality rate, which is constrained within [0, 1]. In our system the
  # mortality rates are typically low, so often an upper limit of 1 is reasonable.
  # In some cases (certain species, seedling data with very short intervals, ...)
  # this upper limit should be increased to allow mortality rates > 1.
  
  # Definition of the upper limit, as 10 times the "typical" mortality rate:
  N = sum(data$status %in% c("A", "D"))
  S = sum(data$status == "A")
  mean.t = mean(data$t[data$status %in% c("A", "D")])
  U = 10 * ((log(N) - log(S))/mean.t)
  
  # (Very loosely) constrained optimization:
  mortality.rate <- optimize(f = Eq3, lower = 0, upper = U, data = data)$minimum
  return(mortality.rate)
}


if(FALSE) # example -- if FALSE it will not run by doing source
{
  # Comparison between methods:
  N = 1000
  data <- data.frame(id = 1:N, 
                     status = sample(c("A", "D"), size = N, prob = c(0.6, 0.4), replace = TRUE), 
                     t = runif(N, min = 1, max = 8))
  
  # Typical calculation of the mortality rate,
  # assuming the average time for all the observations:
  N = sum(data$status %in% c("A", "D"))
  S = sum(data$status == "A")
  mean.t = mean(data$t[data$status %in% c("A", "D")])
  (log(N) - log(S))/mean.t
  
  # Alternative calculation
  # (should be very similar).
  mortality_rate_Kubo(data)
  
  # Using extremely high mortality rates (>1)
  data2 <- data
  data2$t <- data2$t / 100
  
  N = sum(data2$status %in% c("A", "D"))
  S = sum(data2$status == "A")
  mean.t = mean(data2$t[data2$status %in% c("A", "D")])
  (log(N) - log(S))/mean.t
  
  mortality_rate_Kubo(data2)
}