Demo codes to calculate mortality rates, then fit with flexsurv

I generated some normally distributed random lifespan, fit with flexsurv for Gomertz and Weibull. Interestingly, it fits Weibull better than Gompertz. 


source("lifespan.r")

lifespan = round(rnorm(1000)*10, 0) + 30

tb = calculate.s(lifespan)
head(tb)
tb$ds=NA; tb$dt=NA

tb$dt[1] = tb$s[1]
tb$ds[1] = 1 - tb$s[1]
tb$mortality.rate[1] = tb$ds[1] / tb$dt[1]

for( j in 2:length(tb[,1])) {
  tb$ds[j] =  tb$s[j-1] - tb$s[j] 
  tb$dt[j] = -tb$t[j-1] + tb$t[j]
  tb$mortality.rate[j] = tb$ds[j] / ( tb$s[j] * tb$dt[j])
}
plot( tb$s ~ tb$t)
plot( log10(tb$mortality.rate) ~ tb$t ) #linear for Gompertz
plot( log10(tb$mortality.rate) ~ log10(tb$t ) ) #linear for Weibull

tb2 = tb2[-length(tb2[1,]), ]

summary(lm(log10(tb2$mortality.rate) ~ tb2$t, na.rm=T ))

summary(lm(log10(tb$mortality.rate) ~ log10(tb$t), na.rm=T ))

log(m)~t is the linear Gompertz plot.  It can be seen that Gompertz is better suited for the early life stage.

log(m)~log(t) is the linear Weibull plot. Interesting, normally distributed lifespan is better fitted with Weibull model, especially for the late life stage. 

Not surprsingly, Weibull model give better likelihood fitting.