R/lifetable.R

In robjhyndman/demography: Forecasting Mortality, Fertility, Migration and Population Data

# Lifetable functions
# Produce lifetable from mortality rates

#' Construct lifetables from mortality rates
#'
#' Computes period and cohort lifetables from mortality rates for multiple years.
#'
#' For period lifetables, all years and all ages specified are included in the tables. For cohort lifetables,
#' if \code{ages} takes a scalar value, then the cohorts are taken to be of that age in each year contained in \code{years}.
#' But if \code{ages} is a vector of values, then the cohorts are taken to be of those ages in the first year contained in \code{years}.
#'
#' For example, if \code{ages=0} then lifetables of the birth cohorts for all years in \code{years} are computed. On the other hand,
#' if \code{ages=0:100} and \code{years=1950:2010}, then lifetables of each age cohort in 1950 are computed.
#'
#' In all cases, \eqn{q_x = m_x/(1+[(1-a_x)m_x])}{qx = mx/(1 + ((1-ax) * mx))} as per Chiang (1984).
#'
#' Warning: the code has only been tested for data based on single-year age groups.
#'
#' @param data Demogdata object such as obtained from \code{\link{read.demogdata}},
#' \code{\link{forecast.fdm}} or \code{\link{forecast.lca}}.
#' @param series Name of series to use.  Default is the first series in \code{data[["rate"]]}.
#' @param years Vector indicating which years to include in the tables.
#' @param ages Vector indicating which ages to include in table.
#' @param max.age Age for last row. Ages beyond this are combined.
#' @param type Type of lifetable: \code{period} or \code{cohort}.
#'
#' @return Object of class \dQuote{lifetable} containing the following components:
#' \item{label}{Name of region from which data are taken.}
#' \item{series}{Name of series}
#' \item{age}{Ages for lifetable}
#' \item{year}{Period years or cohort years}
#' \item{mx}{Death rate at age x.}
#' \item{qx}{The probability that an individual of exact age x will die before exact age x+1.}
#' \item{lx}{Number of survivors to exact age x.  The radix is 1.}
#' \item{dx}{The number of deaths between exact ages x and x+1.}
#' \item{Lx}{Number of years lived between exact age x and exact age x+1.}
#' \item{Tx}{Number of years lived after exact age x.}
#' \item{ex}{Remaining life expectancy at exact age x.}
#' Note that the lifetables themselves are not returned, only their components. However, there is a print method that constructs (and returns)
#' the lifetables from the above components.
#'
#' @seealso \code{\link{life.expectancy}}
#' @author Heather Booth, Leonie Tickle, Rob J Hyndman, John Maindonald and Timothy Miller
#' @references
#' Chiang CL. (1984) \emph{The life table and its applications}. Robert E Krieger Publishing Company: Malabar.
#'
#' Keyfitz, N, and Caswell, H. (2005) \emph{Applied mathematical demography}, Springer-Verlag: New York.
#'
#' Preston, S.H., Heuveline, P., and Guillot, M. (2001) \emph{Demography: measuring and modeling population processes}. Blackwell
#'
#' @examples
#' france.lt <- lifetable(fr.mort)
#' plot(france.lt)
#' lt1990 <- print(lifetable(fr.mort,year=1990))
#'
#' france.LC <- lca(fr.mort)
#' france.fcast <- forecast(france.LC)
#' france.lt.f <- lifetable(france.fcast)
#' plot(france.lt.f)
#'
#' # Birth cohort lifetables, 1900-1910
#' france.clt <- lifetable(fr.mort,type="cohort",age=0, years=1900:1910)
#'
#' # Partial cohort lifetables for 1950
#' lifetable(fr.mort, years=1950)
#' @keywords models
#' @export
lifetable <- function(data, series=names(data$rate)[1], years=data$year, ages=data$age,
  max.age=min(100,max(data$age)), type=c("period","cohort"))
{
  if(!is.element("demogdata",class(data)))
    stop("data must be a demogdata object")
  if(data$type != "mortality")
    stop("data must be a mortality object")
  type <- match.arg(type)
  if(!is.el(series,names(data$rate)))
    stop(paste("Series",series,"not found"))
  if(is.na(sum(match(years,data$year))))
    stop("Years not present in data")
  sex <- series
  if(sex!="female" & sex!="male" & sex!="total")
  {
    if(is.element("model",names(data)))
      sex <- names(data$model)[4]
  }
  na <- length(data$age)
  if(na > 4)
    agegroup <- data$age[na-1]-data$age[na-2]
  else
    stop("Insufficient age groups")

  if(type=="period")
  {
    max.age <- min(max.age,max(ages))
    data <- extract.ages(data,ages,combine.upper=FALSE)
    if(max.age < max(ages))
      data <- extract.ages(data,min(ages):max.age,combine.upper=TRUE)
    data <- extract.years(data,years=years)
    mx <- get.series(data$rate,series)
    n <- length(years)
    p <- nrow(mx)
    rownames(mx) <- ages <- data$age
    colnames(mx) <- years
    qx <- lx <- dx <- Lx <- Tx <- ex <- rx <- mx*NA
    rownames(rx) <- ages - 1
    rownames(rx)[ages==0] <- "B"
    for(i in 1:n)
    {
      ltable <- lt(mx[,i],min(ages),agegroup,sex)
      nn <- length(ltable$qx)
      qx[1:nn,i] <- ltable$qx
      lx[1:nn,i] <- ltable$lx
      dx[1:nn,i] <- ltable$dx
      Lx[1:nn,i] <- ltable$Lx
      Tx[1:nn,i] <- ltable$Tx
      ex[1:nn,i] <- ltable$ex
      rx[1:nn,i] <- ltable$rx
    }
  }
  else if(type=="cohort" & length(ages)>1) # multiple ages, single year.
  {
    data <- extract.ages(data,min(ages):max.age,combine.upper=TRUE)
    data <- extract.years(data,years=seq(min(years),max(data$year),by=1))
    n <- length(data$year)
    p <- length(data$age)
    cmx <- matrix(NA,p,p)
    rownames(cmx) <- data$age
    colnames(cmx) <- paste(min(years)," age ",data$age,sep="")
    qx <- dx <- Tx <- lx <- Lx <- ex <- cmx
    cohort <- match(ages,data$age)
    cohort <- cohort[!is.na(cohort)]
    if(length(cohort)==0)
      stop("No data available")
    if(min(data$age[cohort]+n-1) < max.age)
      warning("Insufficient data for other lifetables. Try reducing max.age")
    for(coh in cohort)
    {
      if(data$age[coh]+n-1 > max.age)
      {
        subdata <- extract.ages(data,data$age[coh]:max.age,combine.upper=TRUE)
        mx <- get.series(subdata$rate,series)
        p <- nrow(mx)
        for (j in 1:p)
          cmx[coh+j-1,coh] <- mx[j,j]
        ltable <- lt(cmx[coh+(1:p)-1,coh], data$age[coh], agegroup, sex=sex)
        p <- length(ltable$lx)
        lx[coh+(1:p)-1,coh] <- ltable$lx
        Lx[coh+(1:p)-1,coh] <- ltable$Lx
        ex[coh+(1:p)-1,coh] <- ltable$ex
        qx[coh+(1:p)-1,coh] <- ltable$qx
        dx[coh+(1:p)-1,coh] <- ltable$dx
        Tx[coh+(1:p)-1,coh] <- ltable$Tx
      }
    }
    mx <- cmx
    # Retain columns in required cohort
    mx <- mx[,cohort,drop=FALSE]
    lx <- lx[,cohort,drop=FALSE]
    Lx <- Lx[,cohort,drop=FALSE]
    ex <- ex[,cohort,drop=FALSE]
    qx <- qx[,cohort,drop=FALSE]
    dx <- dx[,cohort,drop=FALSE]
    Tx <- Tx[,cohort,drop=FALSE]
    rx <- NULL
  }
  else #single age, multiple years.
  {
    data <- extract.years(data,years=seq(min(years),max(data$year),by=1))
    data <- extract.ages(data,ages:max.age,combine.upper=TRUE)
    n <- length(data$year)
    p <- length(data$age)
    ny <- length(years)
    cmx <- matrix(NA,p,ny)
    rownames(cmx) <- data$age
    colnames(cmx) <- paste(years," age ",ages,sep="")
    qx <- dx <- Tx <- lx <- Lx <- ex <- cmx
    minage <- max.age
    for(i in 1:ny)
    {
      subdata <- extract.years(data,years=seq(years[i],max(data$year),by=1))
      upperage <- min(ages+length(subdata$year)-1)
      minage <- min(minage,upperage)
      if(upperage >= max.age)
      {
        mx <- get.series(subdata$rate,series)
        p <- nrow(mx)
        for (j in 1:p)
          cmx[j,i] <- mx[j,j]
        ltable <- lt(cmx[,i],ages, agegroup, sex=sex)
        p <- length(ltable$lx)
        lx[1:p,i] <- ltable$lx
        Lx[1:p,i] <- ltable$Lx
        ex[1:p,i] <- ltable$ex
        qx[1:p,i] <- ltable$qx
        dx[1:p,i] <- ltable$dx
        Tx[1:p,i] <- ltable$Tx
      }
    }
    mx <- cmx
    rx <- NULL
  }

  return(structure(list(age=ages,year=years, mx=mx,qx=qx,lx=lx,dx=dx,Lx=Lx,Tx=Tx,ex=ex,rx=rx,
        series=series, type=type, label=data$label),class="lifetable"))
}



lt <- function (mx, startage = 0, agegroup = 5, sex)
{
  # Omit missing ages
  if (is.na(mx[1]))
    mx[1] <- 0
  firstmiss <- (1:length(mx))[is.na(mx)][1]
  if (!is.na(firstmiss))
    mx <- mx[1:(firstmiss - 1)]
  nn <- length(mx)
  if (nn < 1)
    stop("Not enough data to proceed")

  # Compute width of each age group
  if (agegroup == 1)
    nx <- c(rep(1, nn - 1), Inf)
  else if (agegroup == 5) # First age group 0, then 1-4, then 5-year groups.
    nx <- c(1, 4, rep(5, nn - 3), Inf)
  else stop("agegroup must be either 1 or 5")

  if (agegroup == 5 & startage > 0 & startage < 5)
    stop("0 < startage < 5 not supported for 5-year age groups")

  if (startage == 0) # for single year data and the first age (0) in 5-year data
  {
    if (sex == "female")
    {
      if (mx[1] < 0.107)
        a0 <- 0.053 + 2.8 * mx[1]
      else a0 <- 0.35
    }
    else if (sex == "male")
    {
      if (mx[1] < 0.107)
        a0 <- 0.045 + 2.684 * mx[1]
      else a0 <- 0.33
    }
    else # if(sex == "total")
    {
      if (mx[1] < 0.107)
        a0 <- 0.049 + 2.742 * mx[1]
      else a0 <- 0.34
    }
  }
  else if (startage > 0)
    a0 <- 0.5
  else
    stop("startage must be non-negative")
  if (agegroup == 1)
  {
    if (nn > 1)
      ax <- c(a0, rep(0.5, nn - 2), Inf)
    else
      ax <- Inf
  }
  else if (agegroup == 5 & startage == 0)
  {
    if (sex == "female")
    {
      if (mx[1] < 0.107)
        a1 <- 1.522 - 1.518 * mx[1]
      else
        a1 <- 1.361
    }
    else if (sex == "male")
    {
      if (mx[1] < 0.107)
        a1 <- 1.651 - 2.816 * mx[1]
      else
        a1 <- 1.352
    }
    else # sex == "total"
    {
      if (mx[1] < 0.107)
        a1 <- 1.5865 - 2.167 * mx[1]
      else a1 <- 1.3565
    }
    ax <- c(a0, a1, rep(2.6, nn - 3), Inf)
    ### ax=2.5 known to be too low esp at low levels of mortality
  }
  else # agegroup==5 and startage > 0
  {
    ax <- c(rep(2.6, nn - 1), Inf)
    nx <- c(rep(5, nn))
  }
  qx <- nx * mx/(1 + (nx - ax) * mx)
  # age <- startage + cumsum(nx) - 1
  # if (max(age) >= 75) {
  #    idx <- (age >= 75)
  #   ax[idx] <- (1/mx + nx - nx/(1 - exp(-nx * mx)))[idx]
  #  qx[idx] <- 1 - exp(-nx * mx)[idx]
  #    }
  #qx[qx > 1] <- 1  ################  NOT NEEDED IN THEORY

  #plot(qx)  #### TO CHECK RESULT RE QX>1

  qx[nn] <- 1
  if (nn > 1)
  {
    lx <- c(1, cumprod(1 - qx[1:(nn - 1)]))
    dx <- -diff(c(lx, 0))
  }
  else
    lx <- dx <- 1
  Lx <- nx * lx - dx * (nx - ax)
  Lx[nn] <- lx[nn]/mx[nn]
  Tx <- rev(cumsum(rev(Lx)))
  ex <- Tx/lx
  if (nn > 2)
    rx <- c(Lx[1]/lx[1], Lx[2:(nn - 1)]/Lx[1:(nn - 2)], Tx[nn]/Tx[nn-1])
  else if (nn == 2)
    rx <- c(Lx[1]/lx[1], Tx[nn]/Tx[nn - 1])
  else
    rx <- c(Lx[1]/lx[1])
  if (agegroup == 5)
    rx <- c(0, (Lx[1] + Lx[2])/5 * lx[1], Lx[3]/(Lx[1]+Lx[2]),
            Lx[4:(nn - 1)]/Lx[3:(nn - 2)], Tx[nn]/Tx[nn-1])
  result <- data.frame(ax = ax, mx = mx, qx = qx, lx = lx,
      dx = dx, Lx = Lx, Tx = Tx, ex = ex, rx = rx, nx = nx)
  return(result)
}


#' Plot life expectancy from lifetable
#'
#' plots life expectancy for each age and each year as functional time series.
#'
#' @param x Output from \code{\link{lifetable}}.
#' @param years Years to plot. Default: all available years.
#' @param main Main title.
#' @param xlab Label for x-axis.
#' @param ylab Label for y-axis.
#' @param ... Additional arguments passed to \code{\link[rainbow]{plot.fds}}.
#'
#' @seealso \code{\link{life.expectancy}}, \code{\link{lifetable}}.
#' @author Rob J Hyndman
#' @examples
#' france.lt <- lifetable(fr.mort)
#' plot(france.lt)
#'
#' france.LC <- lca(fr.mort)
#' france.fcast <- forecast(france.LC)
#' france.lt.f <- lifetable(france.fcast)
#' plot(france.lt.f,years=2010)
#' @export
#' @keywords models
plot.lifetable <- function(x,years=x$year,main,xlab="Age",ylab="Expected number of years left",...)
{
  if(x$type != "period")
    stop("Currently only period lifetables can be plotted.")
  # Extract years
  idx <- match(years,x$year)
  idx <- idx[!is.na(idx)]
  idx <- idx[idx <= ncol(x$ex)]
  if(length(idx)==0)
    stop("Year not available")
  years <- x$year[idx]
  ny <- length(years)

  if(missing(main))
  {
    main <- paste("Life expectancy:",x$label,x$series)
    if(ny>1)
      main <- paste(main," (",min(years),"-",max(years),")",sep="")
    else
      main <- paste(main," (",years,")",sep="")
  }

  plot(fts(x$age,x$ex[,idx],start=years[1],frequency=1),main=main,ylab=ylab,xlab=xlab,...)
}

#' @rdname plot.lifetable
#' @export
lines.lifetable <- function(x,years=x$year,...)
{
  if(x$type != "period")
    stop("Currently only period lifetables can be plotted.")
  # Extract years
  idx <- match(years,x$year)
  idx <- idx[!is.na(idx)]
  idx <- idx[idx <= ncol(x$ex)]
  if(length(idx)==0)
    stop("Year not available")
  years <- x$year[idx]
  ny <- length(years)

  lines(fts(x$age,x$ex[,idx],start=x$year[1],frequency=1),...)
}

#' @export
as.data.frame.lifetable <- function(x, ...) {
  ny <- ncol(x$ex)
  if(x$type=="period") {
    year <- x$year
  } else {
    year <- as.numeric(colnames(x$ex))
  }
  outlist <- vector(length=ny,mode="list")
  for(i in seq(ny)) {
    idx2 <- !is.na(x$mx[,i])
    if(sum(idx2)>0) {
      outlist[[i]] <- data.frame(year[i], as.numeric(rownames(x$ex)), x$mx[,i],x$qx[,i],x$lx[,i],x$dx[,i],x$Lx[,i],x$Tx[,i],x$ex[,i])[idx2,]
      colnames(outlist[[i]]) <- c("year","x","mx","qx","lx","dx","Lx","Tx","ex")
    }
  }
  out <- do.call("rbind", outlist)
  rownames(out) <- NULL
  return(out)
}

#' @export
print.lifetable <- function(x,digits=4,...)
{
  out <- as.data.frame(x)
  if(x$type=="period") {
    cat("Period ")
  } else {
    cat("Cohort ")
  }
  cat(paste("lifetable for",x$label,":",x$series,"\n\n"))
  print(round(out, digits = digits))
  invisible(out)
}

# Compute expected age from single year mortality rates
get.e0 <- function(x,agegroup,sex,startage=0)
{
  lt(x, startage, agegroup, sex)$ex[1]
}

# Compute expected ages for multiple years
#' Estimate life expectancy from mortality rates
#'
#' All three functions estimate life expectancy from \code{lifetable}.
#' The function \code{flife.expectancy} is primarily designed for forecast life expectancies and will optionally
#' produce prediction intervals. Where appropriate, it will package the results as a forecast object
#' which makes it much easier to product nice plots of forecast life expectancies.
#' The \code{e0} function is a shorthand wrapper for \code{flife.expectancy} with \code{age=0}.
#'
#' @return Time series of life expectancies (one per year), or a forecast object of life expectancies (one per year).
#'
#' @param data Demogdata object of type \dQuote{mortality} such as obtained from \code{\link{read.demogdata}},
#' or an object of class \code{fmforecast} such as the output from  \code{\link{forecast.fdm}} or \code{\link{forecast.lca}},
#' or an object of class \code{fmforecast2} such as the output from \code{\link{forecast.fdmpr}}.
#' @param series Name of mortality series to use. Default is the first demogdata series in data.
#' @param years Vector indicating which years to use.
#' @param type Either \code{period} or \code{cohort}.
#' @param age Age at which life expectancy is to be calculated.
#' @param max.age Maximum age for life table calculation.
#' @param PI If TRUE, produce a prediction interval.
#' @param nsim Number of simulations to use when computing a prediction interval.
#' @param ... Other arguments passed to \code{simulate} when producing prediction intervals.
#'
#' @seealso \code{\link{lifetable}}
#' @author Rob J Hyndman
#' @examples
#' plot(life.expectancy(fr.mort),ylab="Life expectancy")
#'
#' france.LC <- lca(fr.mort,adjust="e0",years=1950:1997)
#' france.fcast <- forecast(france.LC,jumpchoice="actual")
#' france.e0.f <- life.expectancy(france.fcast)
#'
#' france.fdm <- fdm(extract.years(fr.mort,years=1950:2006))
#' france.fcast <- forecast(france.fdm)
#' \dontrun{
#'   e0.fcast <- e0(france.fcast,PI=TRUE,nsim=200)
#'   plot(e0.fcast)}
#'
#' life.expectancy(fr.mort,type='cohort',age=50)
#'
#' @keywords models
#' @export
life.expectancy <- function(data,series=names(data$rate)[1],years=data$year,
    type=c("period","cohort"), age=min(data$age), max.age=min(100,max(data$age)))
{
  type <- match.arg(type)
  if(!is.el(series,names(data$rate)))
    stop(paste("Series",series,"not found"))
  if(is.null(max.age))
    max.age <- min(100,max(data$age))
  if(age > max.age | age > max(data$age))
    stop("age is greater than maximum age")
  else if(age < min(data$age))
    stop("age is less than minimum age")
  if(type=="period")
    data.lt <- lifetable(data,series,years,type=type,max.age=max.age)$ex
  else
   data.lt <- lifetable(data,series,years,type=type,ages=age,max.age=max.age)$ex
  idx <- match(age,rownames(data.lt))
  #if(sum(is.na(data.lt[idx,]))>0
  # max(data.lt[idx,]) > 1e9)
  #    warning("Some missing or infinite values in the life table calculation.\n  These can probably be avoided by setting max.age to a lower value.")

  return(ts(data.lt[idx,],start=years[1],frequency=1))
}

#' @rdname life.expectancy
#' @export
flife.expectancy <- function(data, series=NULL, years=data$year,
    type=c("period","cohort"), age, max.age=NULL, PI=FALSE, nsim=500, ...)
{
  type <- match.arg(type)
  if(is.element("fmforecast",class(data)))
  {
    if(data$type != "mortality")
       stop("data not a mortality object")
    hdata <- list(year=data$model$year, age=data$model$age,
         type=data$type, label=data$model$label, lambda=data$lambda)
    hdata$rate <- list(data$model[[4]])
    if(min(hdata$rate[[1]],na.rm=TRUE) < 0 | !is.null(data$model$ratio)) # Transformed
      hdata$rate <- list(InvBoxCox(hdata$rate[[1]],data$lambda))
    if(type=="cohort")
    {
      hdata$year <- c(hdata$year,data$year)
      hdata$rate <- list(cbind(hdata$rate[[1]],data$rate[[1]]))
    }
    names(hdata$rate) <- names(data$model)[4]
    if(!is.null(data$model$pop))
    {
      hdata$pop = list(data$model$pop)
      names(hdata$pop) <- names(hdata$rate)
      if(type=="cohort") # Add bogus population for future years
      {
        n <- ncol(hdata$pop[[1]])
        h <- length(hdata$year)-n
        hdata$pop[[1]] <- cbind(hdata$pop[[1]],matrix(rep(hdata$pop[[1]][,n],h),nrow=nrow(hdata$pop[[1]]),ncol=h))
      }
    }
    class(hdata) <- "demogdata"
    # Fix missing values. Why are they there?
    hdata$rate[[1]][is.na(hdata$rate[[1]])] <- 1-1e-5
    if(is.null(max.age))
      max.age <- max(data$age)
    if(missing(age))
      age <- min(hdata$age)

    x <- stats::window(life.expectancy(hdata,type=type,age=age,max.age=max.age),end=max(data$model$year))
    xf <- life.expectancy(data,years=years,type=type,age=age,max.age=max.age)
    if(type=="cohort")
    {
      xf <- ts(c(stats::window(x,start=max(data$model$year)-max.age+age+1, extend=TRUE),xf),end=max(time(xf)))
      if(sum(!is.na(xf)) > 0)
        xf <- stats::na.omit(xf)
      else
        xf <- stop("Not enough data to continue")
      if(min(time(x)) > max(data$model$year)-max.age+age)
        x <- NULL#ts(NA,end=min(time(xf))-1)
      else
        x <- stats::window(x,end=max(data$model$year)-max.age+age)
    }

    out <- structure(list(x=x,mean=xf,method="FDM model"),class="forecast")
    if(is.element("lca",class(data$model)))
      out$method = "LC model"
    else if(!is.null(data$product))
      out$method = "Coherent FDM model"
    if(PI) # Compute prediction intervals
    {
      e0calc <- (!is.element("product",names(data$rate)) & !is.element("ratio",names(data$rate)) & is.null(data$model$ratio))
      if(is.null(data$product) & is.null(data$var) & is.null(data$kt.f))
        warning("Incomplete information. Possibly this is from a coherent\n  model and you need to pass the entire object.")
      else
      {
        sim <- simulate(data,nsim,adjust.model.var=FALSE,...)
        if(type=="cohort") # Add actual rates for first few years
        {
          usex <- length(x)*any(!is.na(x))
          ny <- length(data$model$year) - usex
          sim2 <- array(NA,c(dim(sim)[1],dim(sim)[2]+ny,dim(sim)[3]))
          sim2[,(ny+1):dim(sim2)[2],] <- sim
          hrates <- hdata$rate[[1]][,usex + (1:ny)]
          sim2[,1:ny,] <- array(rep(hrates,dim(sim)[2]),c(dim(sim)[1],ny,dim(sim)[3]))
          sim <- sim2
          rm(sim2)
        }
        if(e0calc)
        {
          e0sim <- matrix(NA,dim(sim)[2],dim(sim)[3])
          simdata <- data
          if(type=="cohort")
            simdata$year <- min(time(out$mean))-1 + 1:dim(sim)[2]
          for(i in 1:dim(sim)[3])
          {
            simdata$rate <- list(as.matrix(sim[,,i]))
            names(simdata$rate) <- names(data$rate)[1]
            e0sim[,i] <- life.expectancy(simdata,type=type,age=age,max.age=max.age)
          }
          e0sim <- e0sim[1:length(xf), , drop = FALSE]
          if(is.element("lca",class(data$model)))
            out$level <- data$kt.f$level
          else
            out$level <- data$coeff[[1]]$level
          out$lower <- ts(apply(e0sim,1,quantile,prob=0.5 - out$level/200,na.rm=TRUE))
          out$upper <- ts(apply(e0sim,1,quantile,prob=0.5 + out$level/200,na.rm=TRUE))
          stats::tsp(out$lower) <- stats::tsp(out$upper) <- stats::tsp(out$mean)
        }
        out$sim <- sim
      }
    }
    return(out)
  }
  else if(is.element("fmforecast2",class(data)))
  {
    if(data[[1]]$type != "mortality")
      stop("data not a mortality object")
    if(is.null(series))
      series <- names(data)[1]
    if(is.element("product",names(data))) # Assume coherent model
    {
      if(missing(age))
        age <- min(data[["product"]]$age)
      if(is.null(max.age))
        max.age <- max(data[["product"]]$age)
      if(series=="total")
      {
        # Compute forecast total mortality rates
        # Assume sex ratios of populations stay the same as last k years
        ny <- length(data[["product"]]$model$year)
        h <- length(data[["product"]]$year)
        k <- 5
        sex.ratio <- (data$female$model$pop / data$male$model$pop)[,ny-(1:k)+1]
        sex.ratio <- apply(sex.ratio, 1, median)
        frate <- data$female$rate$female
        mrate <- data$male$rate$male
        totalrate <- frate/(1+1/sex.ratio) + mrate/(1+sex.ratio)
        # Construct demogdata object with estimated total rates
        # Use last year of population as approximation
        total <- demogdata(totalrate,
                   pop=matrix(data$female$model$pop[,ny] +
                        data$male$model$pop[,ny],
                              nrow=max.age+1,ncol=h),
                   ages=data[["product"]]$age,
                   years=data[["product"]]$year,
                   type="mortality",
                   label=data$product$model$label,
                   name="total",
                   lambda=0)
        # Compute total life expectancy
        out <- flife.expectancy(total, PI=FALSE, age=age, max.age=max.age, type=type)
        if(PI)
        {
          # Create forecast object
          if(is.element("ts",class(out)))
          {
            out <- structure(list(mean=out, method="Coherent FDM model"), class='forecast')
            frate <- exp(data$female$model$female)
            mrate <- exp(data$male$model$male)
            htotalrate <- frate/(1+1/sex.ratio) + mrate/(1+sex.ratio)
            historictotal <- demogdata(htotalrate,
                   pop=data$female$model$pop + data$male$model$pop,
                   ages=data$female$model$age,
                   years=data$female$model$year,
                   type="mortality",
                   label=data$product$model$label,
                   name="total",
                   lambda=0)
            out$x <- life.expectancy(historictotal)
          }
          # Generate simulated products and ratios
          prodsim <- flife.expectancy(data$product,nsim=nsim,PI=TRUE,age=age,max.age=max.age,type=type)
          data$ratio[["female"]]$model$ratio <- TRUE # To avoid PI calculation on flife.expectancy
          ratiosim <- flife.expectancy(data$ratio[["female"]],nsim=nsim,PI=TRUE,age=age,max.age=max.age,type=type)
          # Simulated future rates
          fsim <- prodsim$sim * ratiosim$sim
          msim <- prodsim$sim / ratiosim$sim
          sim <- fsim/(1+1/sex.ratio) + msim/(1+sex.ratio)
          dimsim <- dim(sim)
          simdata <- total
          nnn <- dim(total$rate[[1]])
          useyears <- (1:nnn[2]) + (dimsim[2]-nnn[2])
          e0sim <- matrix(NA,nnn[2],dimsim[3])
          for(i in 1:dimsim[3])
          {
            simdata$rate <- list(as.matrix(sim[,useyears,i]))
            names(simdata$rate) <- names(total$rate)[1]
            fl <- flife.expectancy(simdata,type=type,age=age,max.age=max.age)
            if(inherits(fl, "ts"))
              e0sim[,i] <- fl
            else if(inherits(fl, "forecast"))
              e0sim[,i] <- fl$mean
            else
              stop("No idea what's going on here.")
          }
          out$level <- data$product$coeff[[1]]$level
          out$lower <- ts(apply(e0sim,1,quantile,prob=0.5 - out$level/200))
          out$upper <- ts(apply(e0sim,1,quantile,prob=0.5 + out$level/200))
          stats::tsp(out$lower) <- stats::tsp(out$upper) <- stats::tsp(out$mean)
        }
      }
      else
      {
        if(missing(age))
          age <- min(data[[series]]$age)
        if(is.null(max.age))
          max.age <- max(data[[series]]$age)
        out <- flife.expectancy(data[[series]],PI=FALSE,age=age,max.age=max.age,type=type)
        out$method <- "Coherent FDM model"
        if(PI)
        {
          # Generate simulated products and ratios
          prodsim <- flife.expectancy(data$product,nsim=nsim,PI=TRUE,age=age,max.age=max.age,type=type)
          data$ratio[[series]]$model$ratio <- TRUE # To avoid PI calculation on flife.expectancy
          ratiosim <- flife.expectancy(data$ratio[[series]],nsim=nsim,PI=TRUE,age=age,max.age=max.age,type=type)
          # Simulated future rates
          sim <- prodsim$sim * ratiosim$sim
          dimsim <- dim(sim)
          simdata <- data[[series]]
          nnn <- dim(data[[series]]$rate[[1]])
          useyears <- (1:nnn[2]) + (dimsim[2]-nnn[2])
          e0sim <- matrix(NA,nnn[2],dimsim[3])
          for(i in 1:dim(sim)[3])
          {
            simdata$rate <- list(as.matrix(sim[,useyears,i]))
            names(simdata$rate) <- names(data[[series]]$rate)[1]
            fl <- flife.expectancy(simdata,type=type,age=age,max.age=max.age)
            if(inherits(fl, "ts"))
              e0sim[,i] <- fl
            else if(inherits(fl, "forecast"))
              e0sim[,i] <- fl$mean
            else
              stop("No idea what's going on here.")
          }
          #browser()
          out$level <- data$product$coeff[[1]]$level
          out$lower <- ts(apply(e0sim,1,quantile,prob=0.5 - out$level/200))
          out$upper <- ts(apply(e0sim,1,quantile,prob=0.5 + out$level/200))
          stats::tsp(out$lower) <- stats::tsp(out$upper) <- stats::tsp(out$mean)
        }
      }
    }
    else
      out <- flife.expectancy(data[[series]],PI=PI,nsim=nsim,max.age=max.age,type=type,age=age)
    return(out)
  }
  else
  {
    if(!is.element("demogdata",class(data)))
      stop("data must be a demogdata object")
    if(data$type != "mortality")
      stop("data must be a mortality object")
    if(is.null(series))
      series <- names(data$rate)[1]
    if(missing(age))
      age <- min(data$age)
    return(life.expectancy(data,series=series,years=years,type=type,age=age,max.age=max.age))
  }
}

#' @rdname life.expectancy
#' @export
e0 <- function(data, series=NULL, years=data$year,
    type=c("period","cohort"), max.age=NULL,PI=FALSE, nsim=500, ...)
{
  flife.expectancy(data, series=series, years=years,age=0,
    type=type, max.age=max.age,PI=PI,nsim=nsim,...)
}