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.
No comments:
Post a Comment