R implementation of lo et al (2016) supplementary material

Lo_etal_2016.R

# First draft of this script by Markus Goehler, light edits by Tim Riffe
#################################### Preparation ###########################################

# Data from Lo et al supplementary material
# https://www.demographic-research.org/volumes/vol35/15/default.htm

mytable<- structure(
              list(x = structure(
                c(1L, 2L, 11L, 3L, 4L, 5L, 6L, 7L, 
                  8L, 9L, 10L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L), 
                 .Label = c("<1", "1-4", "10-14", "15-19", "20-24", "25-29", "30-34", "35-39", 
                         "40-44", "45-49", "5-9", "50-54", "55-59", "60-64", "65-69", 
                         "70-74", "75-79", "80-84", "85-89", "90+"), class = "factor"), 
               Dx = c(23L, 13L, 3L, 3L, 19L, 18L, 23L, 14L, 22L, 61L, 116L, 
                      180L, 219L, 268L, 326L, 503L, 735L, 868L, 849L, 770L), 
               Px = c(6318L,  24605L, 32625L, 40545L, 45060L, 41280L, 40100L, 35245L, 36455L, 
                      54745L, 64465L, 57980L, 49630L, 41975L, 31400L, 28380L, 24300L, 
                      16880L, 8460L, 4325L), 
               ax = c(0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 
                      0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 
                      0.5, 0.5), 
               nx = c(1L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
                      5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L)), class = "data.frame", row.names = c(NA, -20L))        

# get vectors for below
Dx       <-mytable$Dx
Px       <-mytable$Px
ax       <-mytable$ax
nx       <-mytable$nx
age      <-mytable$x
n        <-length(Px)
radix    <-1e5

################################### Chiang Life Expectancy #################################


############################ single elements 1 by 1 ########################################

mx.f<-function(Px,Dx){ 
  mx<-Dx/Px
  return(mx)
}

qx.f<-function(nx,mx,ax){
  qx<-(nx*mx)/(1+((1-ax)*nx*mx))
  qx[length(qx)]<-1
  return(qx)
}

px.f<-function(qx){
  px<-1-qx
  return(px)
}

lx.f<-function(px,radix){
  lx <- c(1,cumprod(px)) * radix
  lx <-lx[1:length(lx)-1]
  return(lx)
}

dx.f<-function(lx,qx) {
  dx<-lx*qx
  return(dx)
}
  

Lx.f<-function(lx,n,nx,ax,mx){
  Lx <- (lx[-n] + lx[-1])*(nx[-n]/2)
  Lx[1]<-(ax[1]*lx[1]+(1-ax[1])*lx[2])
  Lx[n] <- lx[n]/mx[n]
  return(Lx)
}

Tx.f<-function(Lx){
  Tx<-rev(cumsum(rev(Lx)))
  return(Tx)
}

ex.f<-function(Tx,lx){
  ex<-Tx/lx
  return(ex)
}

################################## all elements combined ###################################
lifetable<-function(n,Px,Dx,ax,nx,age,radix) {
  mx <-mx.f(Px,Dx) #death rates
  qx <-qx.f(nx,mx,ax) #probability of dying
  qx[length(qx)]<-1
  px <-px.f(qx) #probability of surviving
  lx <-lx.f(px,radix)
  dx <-dx.f(lx,qx)
  Lx <-Lx.f(lx,n,nx,ax,mx)
  Tx <-Tx.f(Lx)
  ex <-ex.f(Tx,lx)
  result<-data.frame(Age=age,
                     Px=Px,
                     Dx=Dx,
                     ax=ax,
                     nx=nx,
                     mx=mx,
                     qx=qx,
                     px=px,
                     lx=lx,
                     dx=dx,
                     Lx=Lx,
                     Tx=Tx,
                     ex=ex)
  return(result)
}

probe<-lifetable(n,Px,Dx,ax,nx,age,radix)



################################### Standard Chiang Variance #################################
qx<-probe$qx
lx<-probe$lx
ex<-probe$ex

############################ single elements 1 by 1 ########################################

var_qx.f<-function(qx,Dx){
  var_qx          <- qx^2*(1-qx)/Dx
  return(var_qx)
}

v1.f<-function(n,nx,ax,lx,ex,var_qx) {
  v1<-lx^2*((1-ax)*nx+c(ex[-1],0))^2*var_qx
  return(v1)
}



v2.f<-function(v1){
  v2             <- rev(cumsum(rev(v1)))
  return(v2)
}

var_chiang.f<-function(v2,lx){
  var_chiang     <- v2/lx^2
  return(var_chiang)
}

se_chiang.f<-function(var_chiang){
  se_chiang      <- sqrt(var_chiang)
  return(se_chiang)
}
################################## all elements combined ###################################
scv_getxxx<- function(n,qx,Dx,lx,nx,ax,ex,age){
  var_qx          <- var_qx.f(qx,Dx)
  v1              <- v1.f(n,nx,ax,lx,ex,var_qx)
  v2              <- v2.f(v1)
  var_chiang      <- var_chiang.f(v2,lx)
  se_chiang       <- se_chiang.f(var_chiang)
  result          <- data.frame(age=age,
                               Var_qx=var_qx,
                               v1=v1,
                               v2=v2,
                               var_chiang=var_chiang,
                               se_chiang=se_chiang)
  return(result) 
}
probe2<-scv_getxxx(n,qx,Dx,lx,nx,ax,ex,age)



################################## Adjusted Chiang Variance #################################
px<-probe$px
mx<-probe$mx
var_chiang<-probe2$var_chiang
se_chiang<-probe2$se_chiang

############################ single elements 1 by 1 ########################################

paw.f<-function(px){
  paw               <- rev(cumprod(rev(px[1:length(px)-1])))
  paw[length(paw)+1]<- 1
  return(paw)
}

chk_paw.f<-function(lx){
  chk_paw           <- lx[length(lx)]/lx 
  return(chk_paw)
}

term_adj.f<-function(chk_paw,mx,Px){
  term_adj          <- chk_paw^2/(mx[length(mx)]^3*Px[length(Px)])
  return(term_adj)
}

var_adj.f<-function(var_chiang,term_adj){
  var_adj           <- var_chiang+term_adj
  return(var_adj)
}

se_adj.f<-function(var_adj){
  se_adj            <-sqrt(var_adj)
  return(se_adj)
}

################################## all elements combined ###################################
acv_getxxx<- function (px,lx,mx,Px,var_chiang,se_chiang,age) {
  paw               <- paw.f(px)
  chk_paw           <- chk_paw.f(lx)
  term_adj          <- term_adj.f(chk_paw,mx,Px)
  var_adj           <- var_adj.f(var_chiang,term_adj)
  se_adj            <- se_adj.f(var_adj)
  result            <- data.frame(age=age,
                                  paw=paw,
                                  chk_paw=chk_paw,
                                  term_adj=term_adj,
                                  var_adj=var_adj,
                                  se_adj=se_adj)
 return(result)               
}

probe3<-acv_getxxx(px,lx,mx,Px,var_chiang,se_chiang,age)



################################## Adjusted Chiang with Population Error #################################
prop_error=0.05
critic=2
chk_paw<-probe3$chk_paw
var_adj<-probe3$var_adj
se_adj<-probe3$se_adj

############################ single elements 1 by 1 ########################################
term_pop.f<-function(Px,n,chk_paw,Dx,prop_error,critic){
  term_pop          <- (Px[n]*chk_paw/Dx[n])^2*(prop_error/critic)^2
  return(term_pop)
}

var_adj_pop.f<-function(var_adj,term_pop){
  var_adj_pop       <- var_adj + term_pop
  return(var_adj_pop)
}

se_adj_pop.f<-function(var_adj_pop){
  se_adj_pop        <- sqrt(var_adj_pop)
  return(se_adj_pop)
}

################################## all elements combined ###################################
acvpe_getxxx<-function(Px,chk_paw,Dx,prop_err,critic,n,var_adj,se_adj,age) {
  term_pop          <- term_pop.f(Px,n,chk_paw,Dx,prop_error,critic)
  var_adj_pop       <- var_adj_pop.f(var_adj,term_pop)
  se_adj_pop        <- se_adj_pop.f(var_adj_pop)
  result            <- data.frame (Age=age,
                                   term_pop=term_pop,
                                   var_adj_pop=var_adj_pop,
                                   se_adj_pop=se_adj_pop)
  return(result) 
}

probe4<-acvpe_getxxx(Px,chk_paw,Dx,prop_err,critic,n,var_adj,se_adj,age)



################## Standard Chiang with Negative Binomial Overdispersion #################################

############################ single elements 1 by 1 ########################################


Px_minus.f<-function(Px,nx,ax,Dx){
  Px_minus          <- Px/nx+(1-ax)*Dx
  return(Px_minus)
}

Px_plus.f<-function(Px_minus,Dx){
  Px_plus           <- Px_minus - Dx
  return(Px_plus)
}

VarNBQx.f<-function(Dx,Px_plus,Px_minus,n){
  VarNBQx           <- Dx*(1/(Px_plus*Px_minus))
  VarNBQx[n]        <- 0
  return(VarNBQx)
}  

v1.nb.f<-function(lx,ax,nx,ex,VarNBQx,n){
  v1<-lx^2*((1-ax)*nx+c(ex[-1],0))^2*VarNBQx
  return(v1)
}

v2.nb.f<-function(v1){
  v2                <- rev(cumsum(rev(v1)))
  return(v2)
}

var_NB.f<-function(v2,lx){
  var_NB            <- v2/lx^2 
  return(var_NB)
}

se_NB.f<-function(var_NB){
  se_NB           <- sqrt(var_NB)
  return(se_NB)
}

################################## all elements combined ###################################
scnbo_getxxx<- function(age,Px,Dx,ax,nx,lx,ex,n) {
  Px_minus          <- Px_minus.f(Px,nx,ax,Dx)
  Px_plus           <- Px_plus.f(Px_minus,Dx)
  VarNBQx           <- VarNBQx.f(Dx,Px_plus,Px_minus,n)
  v1                <- v1.nb.f(lx,ax,nx,ex,VarNBQx,n)
  v2                <- v2.nb.f(v1)
  var_NB            <- var_NB.f(v2,lx)
  se_NB             <- se_NB.f(var_NB)
  result            <- data.frame(age=age,
                                  Px_minus=Px_minus,
                                  Px_plus=Px_plus,
                                  VarNBQx=VarNBQx,
                                  v1=v1,
                                  v2=v2,
                                  var_NB=var_NB,
                                  se_NB=se_NB)
  return(result)
}

probe5<-scnbo_getxxx(age,Px,Dx,ax,nx,lx,ex,n)





################## Adjusted Chiang with NB Overdispersion #################################
Px_minus<-probe5$Px_minus
Px_plus<-probe5$Px_plus
var_NB<-probe5$var_NB

############################ single elements 1 by 1 ########################################

term_adj_nb.f<-function(mx,Px,n,chk_paw,Px_minus,Px_plus){
  term_adj        <- chk_paw^2/(mx[n]^3*Px[n])*(Px_minus[n]/Px_plus[n])
  return(term_adj)
}

var_adj_nb.f<-function(var_NB,term_adj){
  var_adj_nb      <- var_NB + term_adj
  return(var_adj_nb)
}

se_adj_nb.f<-function(var_adj_nb){
  se_adj_nb     <- sqrt(var_adj_nb)
  return(se_adj_nb)
}

################################## all elements combined ###################################
acnbo_getxxx<-function(age,mx,Px,n,chk_paw,Px_minus,Px_plus,var_NB){
  term_adj        <- term_adj_nb.f(mx,Px,n,chk_paw,Px_minus,Px_plus)
  var_adj_nb      <- var_adj_nb.f(var_NB,term_adj)
  se_adj_nb       <- se_adj_nb.f(var_adj_nb)
  result          <- data.frame(age=age,
                                term_adj=term_adj,
                                var_adj_nb=var_adj_nb,
                                se_adj_nb=se_adj_nb) 
  return(result) 
}

probe6<-acnbo_getxxx(age,mx,Px,n,chk_paw,Px_minus,Px_plus,var_NB)