Life table calculation.

LifeTable

##e.g. data from kmi in LOS_toy_example.R
##non-infected only
#data <- data[-173,]

## (death) count at each event time
nDx <- table(data$exit)

## event times
times <- as.numeric(names(table(data$exit)))

## population just before each event time  
nKx <- max(cumsum(nDx))-c(0, cumsum(nDx)[-length(nDx)])

## hazard of death at each event time
## e.g. age-specific death rate
##"bathtub-shaped" hazard cf. reliability theory
nMx <- nDx/nKx

## residual lifetime
plot(lt$times, lt$ex, type="b")
plot(lt$times, lt$times+lt$ex, type="b")
abline(a=0,b=1)

## exponential hazard (Gompertz) for >30 days
#plot(times[times>30], log(nMx[times>30]), type="b")
#plot(times[times<=30], log(nMx[times<=30]), type="b")
#plot(times, nMx, type="b")
#fit <- lm(log(nMx[times>30])~times[times>30])
#fit <- lm(log(nMx[times<=30])~times[times<=30])
#abline(a=fit$coefficients[1], b=fit$coefficients[2])
a <- fit$coefficients[1]
b <- fit$coefficients[2]


et <- (log(1+(b/nMx))-(0.5*((1+(nMx/b))^(-2))))/b
#et_hat <- (log(b/nMx)-0.5)/b
lines(times, et, type="b", col="red")


lifeTable <- function(times, nMx){
  ## simple lifetable using Keyfitz and Flieger separation factors and 
  ## exponential tail of death distribution (to close out life table)
  ## http://www.lho.org.uk/Download/Public/7656/1/tech_supp_3.pdf
  ## call: lt <- lifeTable(times, nMx)
  
  #b0 <- 0.07
  #b1 <- 1.7

  nmax <- length(times)
  n <- c(diff(times), 999)          # width of the intervals
  
  ## Average proportion of the year lived by those who die
  nax <- n/2                # midway point. default to .5 of interval
  #nax[1] <- b0 + b1 *nMx[1]    # from Keyfitz & Flieger(1968)
  #nax[2] <- 1.5
  nax[nmax] <- 1/nMx[nmax] 	  # e_x at open age interval
  
  ## probability of dying
  nqx <- (n*nMx) / (1 + (n-nax)*nMx)
  nqx<-ifelse(nqx > 1, 1, nqx) # necessary for high nMx
  nqx[nmax] <- 1.0
  
  ## probability of surviving
  npx <- 1-nqx
  
  lx <- c(1,cumprod(npx))   # (cumulative) survivorship lx
  lx <- lx[1:length(nMx)]
  
  ## number of death during interval
  ndx <- lx * nqx 
  
  ## Number of person years lived through the interval
  nLx <- n*lx - nax*ndx       # equivalent to n*l(x+n) + (n-nax)*ndx
  nLx[nmax] <- lx[nmax]*nax[nmax]
  
  ## Total number of person years lived after the interval
  Tx <- rev(cumsum(rev(nLx)))
  
  ## Expectation of life
  ex <- ifelse(lx[1:nmax] > 0, Tx/lx[1:nmax], NA)
  
  
  lt <- data.frame(times=times, nax=nax, nmx=nMx, nqx=nqx, lx=lx, ndx=ndx, nLx=nLx, Tx=Tx, ex=ex)
  return(lt)
}
## END FUNCTION ##


##e.g. data from kmi in LOS_toy_example.R
##non-infected only
#data <- data[-173,]

## (death) count at each event time
nDx <- table(data$exit)

## event times
times <- as.numeric(names(table(data$exit)))

## population just before each event time  
nKx <- max(cumsum(nDx))-c(0, cumsum(nDx)[-length(nDx)])

## hazard of death at each event time
## e.g. age-specific death rate
##"bathtub-shaped" hazard cf. reliability theory
nMx <- nDx/nKx

## residual lifetime
plot(lt$times, lt$ex, type="b")
plot(lt$times, lt$times+lt$ex, type="b")
abline(a=0,b=1)

## exponential hazard (Gompertz) for >30 days
#plot(times[times>30], log(nMx[times>30]), type="b")
#plot(times[times<=30], log(nMx[times<=30]), type="b")
#plot(times, nMx, type="b")
#fit <- lm(log(nMx[times>30])~times[times>30])
#fit <- lm(log(nMx[times<=30])~times[times<=30])
#abline(a=fit$coefficients[1], b=fit$coefficients[2])
a <- fit$coefficients[1]
b <- fit$coefficients[2]


et <- (log(1+(b/nMx))-(0.5*((1+(nMx/b))^(-2))))/b
#et_hat <- (log(b/nMx)-0.5)/b
lines(times, et, type="b", col="red")


lifeTable <- function(times, nMx){
  ## simple lifetable using Keyfitz and Flieger separation factors and 
  ## exponential tail of death distribution (to close out life table)
  ## http://www.lho.org.uk/Download/Public/7656/1/tech_supp_3.pdf
  ## call: lt <- lifeTable(times, nMx)
  
  #b0 <- 0.07
  #b1 <- 1.7

  nmax <- length(times)
  n <- c(diff(times), 999)          # width of the intervals
  
  ## Average proportion of the year lived by those who die
  nax <- n/2                # midway point. default to .5 of interval
  #nax[1] <- b0 + b1 *nMx[1]    # from Keyfitz & Flieger(1968)
  #nax[2] <- 1.5
  nax[nmax] <- 1/nMx[nmax] 	  # e_x at open age interval
  
  ## probability of dying
  nqx <- (n*nMx) / (1 + (n-nax)*nMx)
  nqx<-ifelse(nqx > 1, 1, nqx) # necessary for high nMx
  nqx[nmax] <- 1.0
  
  ## probability of surviving
  npx <- 1-nqx
  
  lx <- c(1,cumprod(npx))   # (cumulative) survivorship lx
  lx <- lx[1:length(nMx)]
  
  ## number of death during interval
  ndx <- lx * nqx 
  
  ## Number of person years lived through the interval
  nLx <- n*lx - nax*ndx       # equivalent to n*l(x+n) + (n-nax)*ndx
  nLx[nmax] <- lx[nmax]*nax[nmax]
  
  ## Total number of person years lived after the interval
  Tx <- rev(cumsum(rev(nLx)))
  
  ## Expectation of life
  ex <- ifelse(lx[1:nmax] > 0, Tx/lx[1:nmax], NA)
  
  
  lt <- data.frame(times=times, nax=nax, nmx=nMx, nqx=nqx, lx=lx, ndx=ndx, nLx=nLx, Tx=Tx, ex=ex)
  return(lt)
}
## END FUNCTION ##