r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

生命表から平均寿命を計算する

1年生存率と5年生存率 - 驚異のアニヲタ社会復帰への道

に触発されて、5年生存率ではなく平均寿命やMSTを計算したくなった。

前回の記事を使って、簡単にコピペでデータを読み込めるようにする。

dat <- data.frame(numeric( 21 ))
dat $ age = c( 0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100 )
dat $ death = c( 4102,655,590,1802,3370,4170,5952,7469,10238,15754,28964,49579,62258,80829,120825,159362,174185,165385,127573,50503,9578 )
dat $ death_per_100k = c( 72.3,10.4,9.6,28,44.5,50.7,59.6,81.3,128.5,201.6,316.5,474.9,720,1045.2,1729.3,2953,4895.9,8626.9,14694.7,22969.8,35655.2 )
dat <- dat [,-1]

age.dr <- rep(dat$death_per_100k,each = 5)

death <- surv <- numeric(110)
surv[1] <- 1
for (i in 2:106){
      surv[i] <- surv[i-1]*(1 - age.dr[i-1]/100000)
}
surv[107:110] <- (3:0) * surv[106]/4 #110歳までに全員死亡することにする

plot(surv, type = "o", pch = 16) 

 

for(i in 1:109){
death[i] <- surv[i] - surv[i+1]
}
sum(death * (0:109))  #mean survival time
## [1] 81.78
library(survival)
## Loading required package: splines
s1 <- Surv(time = 0:109, event = rep(1, 110)) #これは、0.5~109.5の方がいいかも。
sfit <- survfit(s1~1, weights = death)
plot(sfit)

 

print(sfit, rmean = 'common') #rmean = AUC = mean survial time
## Call: survfit(formula = s1 ~ 1, weights = death)
## 
##    records      n.max    n.start     events     *rmean *se(rmean) 
##      110.0        1.0        1.0        1.0       81.8       14.5 
##     median    0.95LCL    0.95UCL 
##       85.0       62.0         NA 
##     * restricted mean with upper limit =  109

今回の場合平均寿命は81.8才、中央値は85才であった。