男女別の生命表から平均寿命を計算する
厚生労働省の所にデータが有ったので、やり直す。
##男女別生命表
#データフレーム作成 厚生労働省H22国勢調査による完全生命表
dat <- data.frame(numeric( 226 ))
dat $ x = c( 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110 )
dat $ sexM1 = c( 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 )
dat $ lx = c( 100000,99790,99757,99734,99718,99707,99698,99690,99682,99676,99669,99663,99657,99650,99643,99633,99621,99608,99592,99573,99552,99528,99503,99477,99451,99425,99400,99373,99345,99314,99281,99245,99207,99168,99126,99081,99034,98983,98927,98866,98801,98731,98656,98575,98485,98387,98281,98168,98046,97911,97762,97599,97424,97237,97039,96826,96597,96351,96089,95808,95508,95184,94831,94451,94041,93604,93138,92637,92097,91513,90879,90183,89411,88552,87592,86520,85326,83996,82515,80861,79014,76960,74685,72167,69376,66284,62867,59134,55091,50771,46228,41531,36759,31993,27313,22818,18627,14840,11516,8686,6354,4496,3069,2016,1269,764,438,238,122,59,27,11,4,2,1,100000,99754,99716,99690,99672,99659,99647,99637,99628,99619,99612,99603,99594,99583,99570,99555,99536,99512,99482,99445,99401,99350,99293,99232,99169,99105,99041,98977,98911,98845,98777,98709,98638,98565,98489,98409,98325,98236,98139,98034,97918,97793,97656,97508,97346,97170,96978,96768,96538,96284,96006,95702,95370,95006,94608,94172,93694,93171,92601,91981,91308,90568,89764,88902,87980,86994,85938,84804,83588,82290,80904,79413,77807,76074,74199,72156,69928,67496,64851,61985,58902,55622,52169,48550,44767,40849,36852,32862,28943,25141,21495,18047,14876,12021,9506,7343,5529,4051,2880,1982,1317,842,517,303,169,90,45,21,9,4,1 )
dat $ ndx = c( 210,33,23,15,11,9,8,7,7,6,6,6,7,8,10,12,14,16,18,21,24,25,26,26,26,26,27,28,31,33,36,38,40,42,44,48,51,55,61,65,70,75,81,90,98,106,113,122,135,149,163,175,186,199,213,229,246,263,281,300,325,352,380,409,437,466,501,541,584,634,697,772,859,959,1072,1195,1330,1481,1655,1847,2054,2275,2517,2791,3093,3417,3733,4043,4321,4543,4697,4772,4766,4680,4495,4191,3787,3324,2830,2333,1858,1427,1054,746,505,326,200,116,63,32,15,7,3,1,0,246,37,26,18,13,11,10,9,8,8,8,10,11,13,15,19,24,30,37,44,51,57,61,63,64,64,64,66,67,68,68,70,73,76,80,84,89,96,106,115,126,137,149,162,176,192,210,231,254,278,304,332,363,398,436,478,523,570,620,673,739,804,862,922,986,1056,1134,1216,1298,1386,1490,1607,1732,1876,2043,2228,2432,2645,2866,3083,3279,3453,3619,3783,3918,3997,3990,3919,3802,3646,3448,3171,2855,2515,2163,1813,1479,1171,898,665,475,325,214,134,80,45,24,12,5,2,1 )
dat $ npx = c( 0.9979,0.99967,0.99977,0.99985,0.99989,0.99991,0.99992,0.99992,0.99993,0.99994,0.99994,0.99994,0.99993,0.99992,0.9999,0.99988,0.99986,0.99984,0.99981,0.99979,0.99976,0.99975,0.99974,0.99974,0.99974,0.99974,0.99973,0.99972,0.99969,0.99966,0.99964,0.99962,0.9996,0.99958,0.99955,0.99952,0.99948,0.99944,0.99939,0.99934,0.99929,0.99924,0.99918,0.99909,0.999,0.99892,0.99885,0.99876,0.99862,0.99847,0.99833,0.99821,0.99809,0.99796,0.99781,0.99764,0.99746,0.99727,0.99708,0.99687,0.9966,0.9963,0.99599,0.99566,0.99535,0.99502,0.99463,0.99416,0.99366,0.99308,0.99233,0.99144,0.99039,0.98917,0.98776,0.98619,0.98442,0.98237,0.97994,0.97716,0.974,0.97044,0.96629,0.96133,0.95542,0.94845,0.94063,0.93163,0.92157,0.91051,0.8984,0.8851,0.87036,0.85372,0.83542,0.81633,0.7967,0.77602,0.75427,0.73146,0.70758,0.68266,0.65672,0.62979,0.60194,0.57322,0.54373,0.51357,0.48285,0.45173,0.42035,0.38888,0.35753,0.32648,0.29596,0.99754,0.99963,0.99974,0.99982,0.99987,0.99989,0.9999,0.99991,0.99992,0.99992,0.99992,0.9999,0.99989,0.99987,0.99985,0.99981,0.99976,0.9997,0.99962,0.99955,0.99949,0.99943,0.99939,0.99936,0.99936,0.99936,0.99935,0.99934,0.99933,0.99932,0.99931,0.99929,0.99926,0.99923,0.99919,0.99915,0.9991,0.99902,0.99892,0.99882,0.99872,0.9986,0.99848,0.99834,0.99819,0.99802,0.99784,0.99762,0.99737,0.99711,0.99683,0.99653,0.99619,0.99581,0.99539,0.99493,0.99442,0.99388,0.99331,0.99268,0.9919,0.99112,0.99039,0.98963,0.98879,0.98786,0.98681,0.98566,0.98447,0.98315,0.98158,0.97977,0.97773,0.97534,0.97247,0.96913,0.96522,0.96081,0.9558,0.95026,0.94432,0.93792,0.93063,0.92207,0.91248,0.90215,0.89173,0.88074,0.86865,0.85497,0.83959,0.82431,0.80805,0.79078,0.77245,0.75305,0.73256,0.71095,0.68823,0.6644,0.63949,0.61351,0.58652,0.55858,0.52977,0.5002,0.46998,0.43925,0.40818,0.37696,0.34578 )
dat $ nqx = c( 0.0021,0.00033,0.00023,0.00015,0.00011,9e-05,8e-05,8e-05,7e-05,6e-05,6e-05,6e-05,7e-05,8e-05,1e-04,0.00012,0.00014,0.00016,0.00019,0.00021,0.00024,0.00025,0.00026,0.00026,0.00026,0.00026,0.00027,0.00028,0.00031,0.00034,0.00036,0.00038,4e-04,0.00042,0.00045,0.00048,0.00052,0.00056,0.00061,0.00066,0.00071,0.00076,0.00082,0.00091,0.001,0.00108,0.00115,0.00124,0.00138,0.00153,0.00167,0.00179,0.00191,0.00204,0.00219,0.00236,0.00254,0.00273,0.00292,0.00313,0.0034,0.0037,0.00401,0.00434,0.00465,0.00498,0.00537,0.00584,0.00634,0.00692,0.00767,0.00856,0.00961,0.01083,0.01224,0.01381,0.01558,0.01763,0.02006,0.02284,0.026,0.02956,0.03371,0.03867,0.04458,0.05155,0.05937,0.06837,0.07843,0.08949,0.1016,0.1149,0.12964,0.14628,0.16458,0.18367,0.2033,0.22398,0.24573,0.26854,0.29242,0.31734,0.34328,0.37021,0.39806,0.42678,0.45627,0.48643,0.51715,0.54827,0.57965,0.61112,0.64247,0.67352,0.70404,0.00246,0.00037,0.00026,0.00018,0.00013,0.00011,1e-04,9e-05,8e-05,8e-05,8e-05,1e-04,0.00011,0.00013,0.00015,0.00019,0.00024,3e-04,0.00038,0.00045,0.00051,0.00057,0.00061,0.00064,0.00064,0.00064,0.00065,0.00066,0.00067,0.00068,0.00069,0.00071,0.00074,0.00077,0.00081,0.00085,9e-04,0.00098,0.00108,0.00118,0.00128,0.0014,0.00152,0.00166,0.00181,0.00198,0.00216,0.00238,0.00263,0.00289,0.00317,0.00347,0.00381,0.00419,0.00461,0.00507,0.00558,0.00612,0.00669,0.00732,0.0081,0.00888,0.00961,0.01037,0.01121,0.01214,0.01319,0.01434,0.01553,0.01685,0.01842,0.02023,0.02227,0.02466,0.02753,0.03087,0.03478,0.03919,0.0442,0.04974,0.05568,0.06208,0.06937,0.07793,0.08752,0.09785,0.10827,0.11926,0.13135,0.14503,0.16041,0.17569,0.19195,0.20922,0.22755,0.24695,0.26744,0.28905,0.31177,0.3356,0.36051,0.38649,0.41348,0.44142,0.47023,0.4998,0.53002,0.56075,0.59182,0.62304,0.65422 )
dat $ mux = c( 0.06939,0.00045,0.00023,0.00019,0.00013,1e-04,9e-05,8e-05,7e-05,6e-05,6e-05,6e-05,6e-05,7e-05,9e-05,0.00011,0.00013,0.00015,0.00017,2e-04,0.00023,0.00025,0.00026,0.00026,0.00026,0.00026,0.00026,0.00028,3e-04,0.00032,0.00035,0.00037,0.00039,0.00041,0.00044,0.00046,5e-04,0.00054,0.00059,0.00064,0.00068,0.00073,0.00079,0.00087,0.00095,0.00104,0.00111,0.00119,0.00131,0.00145,0.0016,0.00173,0.00185,0.00198,0.00212,0.00228,0.00246,0.00264,0.00283,0.00302,0.00326,0.00355,0.00386,0.00418,0.0045,0.00482,0.00518,0.00561,0.0061,0.00663,0.0073,0.00812,0.0091,0.01025,0.01157,0.01308,0.01476,0.01669,0.01896,0.02162,0.02465,0.02809,0.03203,0.03671,0.04234,0.04909,0.05688,0.0658,0.07604,0.0875,0.10021,0.11432,0.13009,0.14808,0.16863,0.19138,0.21479,0.24009,0.26743,0.29696,0.32888,0.36337,0.40063,0.44089,0.4844,0.53141,0.58221,0.6371,0.69641,0.7605,0.82975,0.90457,0.98543,1.07279,1.16719,0.09375,0.00057,0.00026,0.00022,0.00015,0.00012,0.00011,1e-04,9e-05,8e-05,8e-05,9e-05,1e-04,0.00012,0.00014,0.00017,0.00021,0.00027,0.00034,0.00041,0.00048,0.00055,6e-04,0.00063,0.00064,0.00064,0.00065,0.00066,0.00067,0.00068,0.00069,7e-04,0.00073,0.00076,0.00079,0.00083,0.00088,0.00094,0.00103,0.00113,0.00123,0.00134,0.00146,0.00159,0.00173,0.00189,0.00207,0.00227,0.0025,0.00276,0.00303,0.00332,0.00364,0.004,0.0044,0.00485,0.00534,0.00586,0.00642,0.00701,0.00773,0.00853,0.00929,0.01003,0.01083,0.01172,0.01273,0.01385,0.01503,0.01629,0.01775,0.01947,0.02143,0.02367,0.02636,0.02955,0.03328,0.03759,0.04249,0.04802,0.05407,0.06056,0.0678,0.07629,0.08618,0.09717,0.10871,0.12062,0.13362,0.14839,0.16615,0.18378,0.2029,0.22364,0.24614,0.27056,0.29704,0.32578,0.35695,0.39077,0.42746,0.46726,0.51044,0.55729,0.60811,0.66325,0.72307,0.78796,0.85837,0.93474,1.01761 )
dat $ ex = c( 86.3,85.48,84.51,83.53,82.54,81.55,80.56,79.57,78.57,77.58,76.58,75.59,74.59,73.6,72.6,71.61,70.62,69.63,68.64,67.65,66.67,65.68,64.7,63.71,62.73,61.75,60.76,59.78,58.8,57.81,56.83,55.85,54.87,53.9,52.92,51.94,50.97,49.99,49.02,48.05,47.08,46.11,45.15,44.19,43.23,42.27,41.31,40.36,39.41,38.46,37.52,36.58,35.65,34.72,33.79,32.86,31.94,31.02,30.1,29.19,28.28,27.37,26.47,25.58,24.68,23.8,22.91,22.03,21.16,20.29,19.43,18.58,17.73,16.9,16.08,15.27,14.48,13.7,12.94,12.19,11.46,10.76,10.07,9.4,8.76,8.15,7.56,7.01,6.48,5.99,5.53,5.1,4.7,4.32,3.98,3.66,3.38,3.11,2.87,2.65,2.44,2.25,2.08,1.92,1.77,1.64,1.51,1.4,1.29,1.19,1.1,1.02,0.94,0.88,0.81,79.55,78.75,77.78,76.8,75.81,74.82,73.83,72.84,71.84,70.85,69.85,68.86,67.87,66.87,65.88,64.89,63.9,62.92,61.94,60.96,59.99,59.02,58.05,57.09,56.12,55.16,54.19,53.23,52.26,51.3,50.33,49.37,48.4,47.44,46.48,45.51,44.55,43.59,42.63,41.68,40.73,39.78,38.83,37.89,36.95,36.02,35.09,34.17,33.25,32.33,31.42,30.52,29.63,28.74,27.86,26.98,26.12,25.26,24.42,23.58,22.75,21.93,21.12,20.32,19.53,18.74,17.97,17.2,16.44,15.7,14.96,14.23,13.51,12.81,12.12,11.45,10.79,10.16,9.56,8.98,8.42,7.89,7.38,6.89,6.43,6,5.59,5.21,4.85,4.51,4.19,3.89,3.62,3.36,3.12,2.9,2.69,2.49,2.31,2.14,1.98,1.84,1.7,1.58,1.46,1.35,1.25,1.16,1.07,0.99,0.92 )
dat $ nLx = c( 99837,99773,99745,99726,99713,99703,99694,99686,99679,99673,99666,99660,99654,99647,99638,99627,99615,99600,99583,99563,99540,99516,99490,99464,99438,99413,99386,99359,99329,99297,99263,99226,99188,99147,99104,99058,99008,98955,98897,98834,98766,98694,98616,98531,98437,98335,98225,98108,97980,97838,97681,97512,97331,97139,96934,96713,96476,96221,95950,95660,95348,95010,94644,94249,93825,93374,92891,92371,91809,91201,90537,89804,88989,88081,87066,85934,84673,83269,81703,79954,78005,75842,73447,70796,67856,64602,61027,57138,52952,48515,43889,39148,34372,29642,25046,20693,16697,13138,10060,7479,5387,3749,2514,1620,999,588,329,175,87,41,18,7,3,1,0,99808,99733,99704,99681,99665,99653,99642,99632,99623,99615,99608,99599,99588,99577,99563,99546,99525,99498,99464,99423,99376,99322,99262,99200,99137,99073,99009,98944,98878,98811,98743,98674,98602,98527,98449,98367,98281,98188,98087,97977,97857,97725,97583,97428,97259,97075,96875,96655,96413,96147,95856,95538,95191,94810,94393,93936,93436,92890,92295,91649,90944,90172,89338,88446,87492,86472,85378,84203,82946,81605,80168,78620,76952,75149,73192,71058,68730,66192,63436,60461,57278,53910,50374,46672,42817,38854,34853,30894,27031,23303,19752,16436,13421,10734,8395,6407,4763,3441,2410,1632,1065,669,402,231,126,65,32,14,6,2,1 )
dat $ Tx = c( 8630132,8530295,8430522,8330776,8231051,8131338,8031636,7931942,7832256,7732577,7632904,7533238,7433578,7333924,7234277,7134639,7035012,6935397,6835798,6736215,6636652,6537112,6437596,6338106,6238642,6139204,6039791,5940405,5841046,5741717,5642419,5543156,5443930,5344743,5245596,5146492,5047435,4948426,4849471,4750574,4651740,4552974,4454279,4355663,4257132,4158695,4060360,3962135,3864027,3766047,3668210,3570528,3473016,3375685,3278546,3181612,3084899,2988423,2892202,2796252,2700592,2605244,2510234,2415590,2321342,2227516,2134143,2041252,1948881,1857072,1765872,1675335,1585531,1496542,1408461,1321395,1235461,1150788,1067519,985816,905862,827857,752015,678568,607772,539916,475314,414287,357149,304197,255682,211793,172646,138273,108631,83585,62892,46195,33057,22998,15519,10132,6383,3869,2250,1250,662,332,158,70,29,11,4,1,0,7955005,7855198,7755464,7655761,7556080,7456415,7356762,7257120,7157488,7057865,6958249,6858642,6759043,6659454,6559878,6460315,6360769,6261244,6161746,6062282,5962859,5863483,5764162,5664899,5565699,5466562,5367489,5268480,5169536,5070658,4971847,4873104,4774431,4675829,4577302,4478853,4380486,4282205,4184017,4085929,3987952,3890096,3792370,3694788,3597360,3500100,3403025,3306151,3209496,3113083,3016936,2921080,2825542,2730351,2635541,2541148,2447212,2353776,2260886,2168591,2076941,1985997,1895826,1806488,1718042,1630549,1544077,1458699,1374496,1291551,1209946,1129778,1051158,974206,899057,825865,754807,686077,619885,556449,495988,438711,384801,334427,287756,244938,206085,171231,140337,113307,90003,70252,53816,40395,29661,21266,14859,10096,6655,4245,2613,1548,879,478,247,121,56,24,10,4,1 )
dat <- dat [,-1]
plot(NULL, type = "n", xlim = c(0, 115), ylim = c(0,1), labels = NULL, ylab = "survival rate", xlab = "age")
lines(dat$x[dat$sexM1 == 1], dat$lx[dat$sexM1 == 1]/100000, type = "l", col = "blue") #年齢別生存確率
lines(dat$x[dat$sexM1 == 0], dat$lx[dat$sexM1 == 0]/100000, type = "l", col = "red")
#年の真ん中で死亡したことにするため、年齢に0.5を足すといいようだ。
sum( (dat$x[dat$sexM1 == 1] + 0.5) * dat$ndx[dat$sexM1 == 1] / 100000) #mean survival time M
sum( (dat$x[dat$sexM1 == 0] + 0.5) * dat$ndx[dat$sexM1 == 0] / 100000) #mean survival time F
library(survival)
s1 <- Surv(time = dat$x + 0.5, event = rep(1, length(dat$x)))
sfit <- survfit(s1 ~ sexM1, weights = ndx/10000, data = dat)
plot(sfit)
print(digits = 5, sfit, rmean = 'common') #rmean = AUC = mean survial time
median(rep(dat$x[dat$sexM1 == 1]+0.5, dat$ndx[dat$sexM1 ==1]))
median(rep(dat$x[dat$sexM1 == 0]+0.5, dat$ndx[dat$sexM1 ==0]))
ssp <- smooth.spline( dat$lx[dat$sexM1 == 1]/100000, dat$x[dat$sexM1 == 1])
predict(ssp, x=0.5)
ssp <- smooth.spline( dat$lx[dat$sexM1 == 0]/100000, dat$x[dat$sexM1 == 0])
predict(ssp, x=0.5)
> sum( (dat$x[dat$sexM1 == 1] + 0.5) * dat$ndx[dat$sexM1 == 1] / 100000) #mean survival time M
[1] 79.54972 男性の平均寿命はH22の平均寿命と一致
> sum( (dat$x[dat$sexM1 == 0] + 0.5) * dat$ndx[dat$sexM1 == 0] / 100000) #mean survival time F
[1] 86.30499 女性の平均寿命はH22の平均寿命と一致
> print(digits = 5, sfit, rmean = 'common') #rmean = AUC = mean survial time
Call: survfit(formula = s1 ~ sexM1, data = dat, weights = ndx/10000)
records n.max n.start events *rmean *se(rmean) median 0.95LCL 0.95UCL
sexM1=0 115 10.0004 10.0004 10.0004 86.302 4.1651 89.5 82.5 NA
sexM1=1 111 9.9997 9.9997 9.9997 79.552 4.5299 82.5 74.5 NA
* restricted mean with upper limit = 112.5
rmean=男性女性の平均寿命⇨H22の平均寿命と一致
median＝中央値はH22の報告と不一致
> median(rep(dat$x[dat$sexM1 == 1]+0.5, dat$ndx[dat$sexM1 ==1]))
[1] 82.5
> median(rep(dat$x[dat$sexM1 == 0]+0.5, dat$ndx[dat$sexM1 ==0]))
[1] 89.5
こちらでもmedian＝中央値はH22の報告と不一致
quantile関数ですべての種類をやっても不一致
予測線を書いてmedianを出す
> ssp <- smooth.spline( dat$lx[dat$sexM1 == 1]/100000, dat$x[dat$sexM1 == 1])
> predict(ssp, x=0.5)
$x
[1] 0.5
$y
[1] 82.60474
>
> ssp <- smooth.spline( dat$lx[dat$sexM1 == 0]/100000, dat$x[dat$sexM1 == 0])
> predict(ssp, x=0.5)
$x
[1] 0.5
$y
[1] 89.17287
こちらの方法で中央値もH22の報告と一致した
