# Proportion of people alive in 1945 that are still alive

In demography, we like to use life tables to estimate the probability that someone born in 1945 (say) is still alive nowadays.  But another interesting quantity might be the probability that someone alive in 1945 is still alive nowadays.

The main difference is that we do not know when that person, alive in 1945, was born. Someone who was old in 1945 is very unlikely still alive in 2017. To compute those probabilities, we can use datasets from http://www.mortality.org/hmd/. More precisely, we need both death and birth data. I assume that datasets (text files) were downloaded (it is necessary to register – for free – to get the data).

D=read.table("FRDeaths_1x1.txt",skip=1,header=TRUE) B=read.table("FRBirths.txt",skip=1,header=TRUE)

In the death dataset, there is a “110+” for people older than 110 years. For convenience, let us cap our observations at 110 years old,

D$Age=as.numeric(as.character(D$Age)) D$Age[is.na(D$Age)]=110

Consider now a first function that will return, for people born in 1930 (say) two informations

• the number of people (here, let us consider women only) born in 1930 (from the birth database)
• the number of death of people of age 0 in 1930, people of age 1 in 1931, people of age 2 in 1932, etc…

The code is simple

+   date2=as.Date(paste(data$date_death_y,"-",data$date_death_m,"-",data$date_death_d,sep=""),"%Y-%m-%d") + idx=which(!(is.na(date1)|is.na(date2))) + date1=date1[idx] + date2=date2[idx] + itg=try(age<-age_years(date1,date2),silent=TRUE) + if(inherits(itg, "try-error")) age=trunc((date2-date1)/365.25) + w=weekdays(date2) + T=table(age,w) + Tab=matrix(0,106,7) + for(i in 1:nrow(T)) if(as.numeric(rownames(T)[i])<106) Tab[as.numeric(rownames(T)[i]),]=T[i,] + return(Tab) + } + D <- lapply( seq_len(nrow( go_through)),count_birthday) + T=D[[1]] + for(s in 2:length(D)) T=T+D[[s]] + return(T) + } If we run that function on the three files > D1=TABLE_AGE_DAY("ssdm1") |========================================| 100% > D2=TABLE_AGE_DAY("ssdm2") |========================================| 100% > D3=TABLE_AGE_DAY("ssdm3") |========================================| 100% we can visualize not percentages, as on the figure above, but counts > D=D1+D2+D3 > colnames(D)= c("Sun","Thu","Mon","Tue","Wed","Sat","Fri") > D=D1[, c("Sun","Mon","Tue","Wed","Thu","Fri","Sat")] and we have here (I remove the Saturday to get a better output) > D[,1:6] Sun Mon Tue Wed Thu Fri [1,] 2843 2888 2943 3020 2979 3038 [2,] 2007 1866 1918 1974 1990 2137 [3,] 1613 1507 1532 1530 1515 1613 [4,] 1322 1256 1263 1259 1207 1330 [5,] 1155 1061 1092 1128 1112 1171 [6,] 1067 985 950 1082 1009 1055 [7,] 1129 901 915 954 941 1044 [8,] 1026 927 944 935 911 1005 [9,] 1029 1012 871 908 939 998 [10,] 1093 1011 974 958 928 1018 [11,] 1106 1031 1019 1036 1087 1122 [12,] 1289 1219 1176 1215 1141 1292 [13,] 1618 1455 1487 1484 1466 1633 [14,] 2121 2000 1900 1941 1845 2138 [15,] 2949 2647 2519 2499 2524 2748 [16,] 4488 3885 3798 3828 3747 4267 [17,] 5709 4612 4520 4422 4443 5005 [18,] 7280 5618 5400 5271 5344 5986 [19,] 8086 6172 5833 5820 6004 6628 [20,] 8389 6507 6166 6055 6430 6955 [21,] 8794 7038 6794 6628 6841 7572 [22,] 8578 6528 6512 6472 6757 7342 [23,] 8345 6750 6483 6469 6714 7338 [24,] 8361 6859 6589 6623 6854 7369 [25,] 8398 6974 6892 6766 6964 7613 [26,] 8432 7210 7012 7175 7343 7801 [27,] 8757 7641 7526 7352 7674 7950 [28,] 9190 8041 7843 7851 7940 8268 [29,] 9495 8409 8555 8400 8469 8934 [30,] 9876 9041 9015 9166 9106 9641 [31,] 10567 9952 9506 9634 9770 10212 [32,] 11417 10428 10402 10275 10455 11169 [33,] 11992 11306 11124 11095 11243 11749 [34,] 12665 12327 11760 12025 12137 12443 [35,] 13629 13135 13179 13037 12968 13724 [36,] 14560 14009 13927 13822 14105 14436 [37,] 15660 14990 15013 15009 15101 15700 [38,] 16749 16504 16148 16091 15912 16863 [39,] 17815 17760 17519 17144 17553 17943 [40,] 19366 19057 18918 18517 18760 19604 [41,] 20770 20458 20154 20339 20349 21238 [42,] 21962 22194 22020 21499 21690 22347 [43,] 23803 23922 23701 23681 23437 24227 [44,] 25685 26133 25559 25209 25287 26115 [45,] 27506 28110 27363 27042 27272 28228 [46,] 29366 29744 29555 29245 29678 30444 [47,] 31444 32193 31817 31504 31753 32302 [48,] 33452 34719 33529 33954 33441 34618 [49,] 36186 37150 36005 36064 36226 37138 [50,] 38401 39244 38813 38465 38506 39884 [51,] 40331 41830 41168 41110 40937 42014 [52,] 43181 44351 43975 43949 43579 44734 [53,] 45307 47134 46522 46149 46089 47286 [54,] 47996 49441 49139 48678 48629 49903 [55,] 50635 52424 51757 51433 51477 52550 [56,] 53509 55337 54556 54482 54406 55906 [57,] 55703 58482 58016 57400 57097 58758 [58,] 59016 61453 60652 61024 60557 62473 [59,] 62475 65651 64169 63824 63829 65592 [60,] 66621 69185 68885 68217 68752 69963 [61,] 69759 73144 72421 71784 71745 73414 [62,] 80346 84253 83044 83177 82416 83833 [63,] 86851 90059 89002 88985 89245 90334 [64,] 91839 95465 94602 93985 94154 96195 [65,] 98461 102846 101348 101328 101306 103170 [66,] 104569 108722 107768 107711 107729 109350 [67,] 111230 115477 114418 114743 113935 116356 [68,] 116999 122053 120727 120342 119782 122926 [69,] 123695 128339 127184 126822 126639 129037 [70,] 129956 136123 134555 135120 133842 137390 [71,] 137984 142964 141316 142855 141419 143620 [72,] 145132 150708 148407 149345 149448 151910 [73,] 152877 157993 155861 156349 155924 158725 [74,] 159109 164652 162722 163499 163157 165744 [75,] 165848 172121 170730 170482 170585 173431 [76,] 172457 179036 177185 177328 177392 180215 [77,] 179936 185015 183223 183932 183237 186663 [78,] 185900 191053 189986 189730 189639 193038 [79,] 191498 196694 194246 194810 195246 197812 [80,] 195505 201289 199684 199561 198968 203226 [81,] 199031 204927 202204 202622 202951 205792 [82,] 201589 207928 204929 204001 204396 208224 [83,] 201665 206743 205194 204676 205256 207980 [84,] 200965 205653 203422 202393 203422 206012 [85,] 197445 202692 199498 199730 200075 201728 [86,] 192324 195961 193589 194754 193800 196102 [87,] 183732 188063 185153 186104 186021 188176 [88,] 174258 177474 175822 176078 176761 177449 [89,] 163180 166706 162810 164367 164281 166436 [90,] 149169 151738 150148 150212 150535 152435 [91,] 134218 136866 134959 134922 135027 136381 [92,] 118936 121106 119591 119509 119793 120998 [93,] 102734 104955 102944 102865 103345 104776 [94,] 87418 88885 88023 86963 87546 87872 [95,] 72023 72698 72151 71579 71530 72287 [96,] 56985 58238 57478 57319 57163 57615 [97,] 44447 45058 44607 44469 43888 44868 [98,] 33457 34132 33022 33409 33454 33642 [99,] 24070 24317 24305 24089 24020 24383 [100,] 17165 17295 16755 17115 16957 17207 [101,] 11799 12125 11709 11816 11824 11719 [102,] 7714 7741 7959 7691 7648 7633 [103,] 5024 5012 4822 4792 4882 4916 [104,] 2987 3101 2978 3049 3093 2906 [105,] 1781 1894 1811 1756 1734 1834 So clearly, for young people, the number of deaths is rather small… And to visualize it, as above, we can use > P=D/apply(D,1,sum)*100 > range(P) [1] 12.34857 17.59386 > dP=trunc((P-min(P))/(max(P)+.01-min(P))*11) > library(RColorBrewer) > CLR=rev(brewer.pal(name="RdYlBu", 11)) > plot(0:1,0:1,ylim=c(55,110),xlim=c(-1,7)) > for(i in 1:106){ + for(j in 1:7){ + rect(j-1,108-i,j,107-i,col=CLR[dP[i,j]]) + }} > text(rep(-.5,106),107.5-1:106,0:105,cex=.4) As above, we observe a strong difference among weekdays for the date of death for young people (below 30) which disappear after (even if there is still a sunday effect) # Démographie québécoise Cet automne, j’étais allé voir Mummy, de Xavier Dolan. Et j’avoue avoir adoré le film ! J’ai tout aimé ! La musique, le cadrage, le rythme, les acteurs (même si j’avais du mal à croire Anne Dorval qui restera toujours pour moi la maman des Parent, surtout que j’ai passé mon temps à Montréal à croiser Daniel Brière – l’acteur, pas le joueur des Avalanches – qui habitait à côté de la maison, et qui faisait souvent le marché en même temps que moi… pour moi les deux étaient attachés à tout jamais…. et la mère de la série et celle du film de Xavier Dolan sont assez différentes). Le seul reproche que j’aurais pu faire, ce sont les sous-titres (qui nous étaient imposés, en France). Au Québec, on avait pris l’habitude de vivre sans sous-titres, et on a appris une langue en la parlant, au quotidien. Par exemple, j’ai appris ce que voulait dire câlisser sans avoir à ouvrir un dictionnaire, ou sans demander de traduction. De même que j’ai découvert qu’il existait des variantes, comme décâlisser (et j’ai même fini par comprendre le sens). Mais je n’avais jamais eu besoin de visualiser la traduction du mot. Voir des traductions de mots que j’avais fini par découvrir et comprendre m’a déstabilisé. Je reverrais d’autant mieux le film en DVD, sans les sous-titres qui m’étaient imposés au cinéma. C’est amusant car pendant les vacances, j’en ai profité pour finir Magasin Général de Régis Loisel et Jean-Louis Tripp. Et j’ai eu plaisir à retrouver des expressions québécoises tout au long de la lecture. Même si l’histoire se déroule à Notre Dame du Lac, situé sur le bas Saint Laurent (même si le nom ne semble plus exister depuis quelques années), on est loin des dialogues que l’on peut entendre en région, et on a une version agréable de ce qu’on pourrait entendre à Montréal… En fait, comme cela est indiqué (trop) discrètement, c’est Jimmy Beaulieu (dont j’avais déjà souligné le travail admirable dans un précédant billet) qui a fait la “traduction” des dialogues, pour avoir (comme le disent les premières pages du livres) des dialogues en québécois “qui soient compréhensibles des deux côtés de l’Atlantique”. Ce qui montre bien qu’un français de France peut comprendre le québécois (moyennant un peu de bonne volonté). Je crois que j’aurais bien aimé avoir des sous-titres pour Mummy dans la même langue que celle parlée dans Magasin Général. En lisant en particulier les deux derniers tomes de Magasin Général (oui, j’avais un peu de retard), et les histoires de natalité, je me suis souvenu du travail que l’on avait fait il y a un peu plus d’un an avec Julie, qui était venu faire un stage à l’UQAM, et avec qui on avait découvert la démographie du Lac Saint Jean (certes, par rapport à Magasin Général, on est de l’autre côté du Saint Laurent). les slides (de la soutenance de stage) sont en ligne, # Reinterpreting Lee-Carter Mortality Model Last week, while I was giving my crash course on R for insurance, we’ve been discussing possible extensions of Lee & Carter (1992) model. If we look at the seminal paper, the model is defined as follows # Mais que s’est-il passé pendant la Première Guerre Mondiale? La réponse courte est que des gens sont morts. Beaucoup. Cela étant dit, on ne dit pas grand chose. On peut comparer les pyramides des âges pour mieux comprendre ce qui a pu se passer. Juste avant la guerre (en 1913), la pyramide des âges ressemblait à ça (en utilisant les données de mortality.org) > EXPO <- read.table( + "http://freakonometrics.free.fr/Exposures-France.txt", header=TRUE,skip=2) > EM=EXPO$Male
> EF=EXPO$Female > Y= EXPO$Year
> A= EXPO$Age > I=which(A=="110+") > base=data.frame(Female=EF,Male=EM,Y=Y,Ages=A) > base=base[-I,] > France1913=base[base$Y==1913,]
> France1919=base[base$Y==1919,] > France1913$Ages=as.numeric(
+ as.character(France1913$Ages)) > France1919$Ages=as.numeric(
+ as.character(France1919$Ages)) > France1913=France1913[,c("Male","Female", + "Ages")] > library(pyramid) > plot(c(0,100), c(0,100), type="n", + frame=FALSE, axes=FALSE, xlab="", ylab="", + main="Pyramide des Ages, France, 1913") > pyramidf(France1913, frame=c(10, 75, 0, 90), + Clab="", Lcol="skyblue", Rcol="pink", + Cstep=10, Laxis=0:4*60000, AxisFM="d") En revanche, juste après la guerre (en 1919), la pyramide des âges des âges ressemblait à celle là Continue reading Mais que s’est-il passé pendant la Première Guerre Mondiale? # Men set to live as long as women by 2030? A few months ago, in Men set to live as long as women, figures show, it was mentioned that (in the U.K.) the gap between male and female life expectancy is closing and men could catch up by 2030, according to an adviser for the Office for National Statistics. (the slides are available online http://cass.city.ac.uk/…). # Your Life in Weeks This week, I discovered a picture on http://waitbutwhy.com/, which represent a (so-called) typical human life, in weeks, I found that interesting. But the first problem is that I don’t understand the limit, below: 90 years, that’s not the average life length. That’s not what you should expect to live when you get born. The second problem is that it cannot be as static as it might seem, when you look at the picture. I mean, life expectancy at age 0 is not the same as life expectancy at age 30, or 50. So I did try to make an animated graph, using prospective life tables. Here a code to generate life tables, at different period, for a French population (I distinguish, here male and female) library(demography) france.LC1 <- lca(fr.mort,adjust="e0",series="female",years=c(1900,2100)) france.fcast <- forecast(france.LC1,h=100) L2 <- lifetable(france.fcast) ex2=L2$ex
L1=lifetable(fr.mort,series="female")
ex1=L1$ex exF=cbind(ex1,ex2) france.LC1 <- lca(fr.mort,adjust="e0",series="male",years=c(1900,2100)) france.fcast <- forecast(france.LC1,h=100) L2 <- lifetable(france.fcast) ex2=L2$ex
L1=lifetable(fr.mort,series="male")
ex1=L1$ex exM=cbind(ex1,ex2) Y=colnames(exF) Based on those lifetables, we can extract remaining life expectancy, at various ages (say, for instance 50, 51, 52, etc), for someone born on some given year (say 1950). Based on those expected remaining lifetimes, we can plot picture=function(yearborn=1950,age=50){ k=which(Y==yearborn) M=diag(exM[,k+0:100]) F=diag(exF[,k+0:100]) par(mfrow=c(1,2)) va=0:(52*100-1) plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","light blue","white")[1+ (va>=age*52)*1+(va>(age+M[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week", ylab="Age",main=paste("Man, born on ",yearborn, ", age ",age,sep="")) axis(1) axis(2) plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","pink","white")[1+ (va>=age*52)*1+(va>(age+F[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week", ylab="Age",main=paste("Woman, born on ",yearborn, ", age ",age,sep="")) axis(1) axis(2)} For instance, if we want the graph above, for someone age 30, born in 1980, we use picture(1980,30) Now, if we run a code to get an animated gif, we can get, for someone born in 1950, and for someone born in 2000 Now, if I could get historical datasets, with the average time spent in schools, ages of retirement, etc, I guess I could add it on the graph. But that’s another story… # Smoothing mortality rates This morning, I was working with Julie, a student of mine, coming from Rennes, on mortality tables. Actually, we work on genealogical datasets from a small region in Québec, and we can observe a lot of volatiliy. If I borrow one of her graph, we get something like Since we have some missing data, we wanted to use some Generalized Nonlinear Models. So let us see how to get a smooth estimator of the mortality surface. We will write some code that we can use on our data later on (the dataset we have has been obtained after signing a lot of official documents, and I guess I cannot upload it here, even partially). DEATH <- read.table( "http://freakonometrics.free.fr/Deces-France.txt", header=TRUE) EXPO <- read.table( "http://freakonometrics.free.fr/Exposures-France.txt", header=TRUE,skip=2) library(gnm) D=DEATH$Male
E=EXPO$Male A=as.numeric(as.character(DEATH$Age))
Y=DEATH$Year I=(A<100) base=data.frame(D=D,E=E,Y=Y,A=A) subbase=base[I,] subbase=subbase[!is.na(subbase$A),]

The first idea can be to use a Poisson model, where the mortality rate is a smooth function of the age and the year, something like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s(x,t)}])$that can be estimated using

library(mgcv)
regbsp=gam(D~s(A,Y,bs="cr")+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regbsp,newdata=data.frame(A=a,Y=y,E=1))
vX=trunc(seq(0,99,length=41))
vY=trunc(seq(1900,2005,length=41))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is here

It is also possible to extract the average value of the years, which is the interpretation of the $a_x$ coefficient in the Lee-Carter model,

predAx=function(a) mean(predict(regbsp,newdata=data.frame(A=a,
Y=seq(min(subbase$Y),max(subbase$Y)),E=1)))
plot(seq(0,99),Vectorize(predAx)(seq(0,99)),col="red",lwd=3,type="l")

We have the following smoothed mortality rate

Recall that the Lee-Carter model is

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}a_x+b_x\cdot k_t}])$

where parameter estimates can be obtained using

regnp=gnm(D~factor(A)+Mult(factor(A),factor(Y))+offset(log(E)),
data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regnp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The (crude) mortality surface is

with the following $a_x$ coefficients.

plot(seq(1,99),coefficients(regnp)[2:100],col="red",lwd=3,type="l")

Here we have a lot of coefficients, and unfortunately, on a smaller dataset, we have much more variability. Can we smooth our Lee-Carter model ? To get something which looks like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s_a(x)+s_b(x)\cdot s_k(t)}])$

Actually, we can, and the code is rather simple

library(splines)
knotsA=c(20,40,60,80)
knotsY=c(1920,1945,1980,2000)
regsp=gnm(D~bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3),
bs(subbase$Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))+
offset(log(E)),data=subbase, family=quasipoisson)
BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) BpY=bs(seq(min(subbase$Y),max(subbase$Y)),knots=knotsY,Boundary.knots= range(subbase$Y),degre=3)
predmodel=function(a,y)
predict(regsp,newdata=data.frame(A=a,Y=y,E=1)) v
Z=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is now

and again, it is possible to extract the average mortality rate, as a function of the age, over the years,

BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) Ax=BpA%*%coefficients(regsp)[2:8] plot(seq(0,99),Ax,col="red",lwd=3,type="l") We can then play with the smoothing parameters of the spline functions, and see the impact on the mortality surface knotsA=seq(5,95,by=5) knotsY=seq(1910,2000,by=10) regsp=gnm(D~bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3), bs(Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))
+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regsp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

We now have to use those functions our our small data sample ! That should be fun….

# Job for life ? Bishop of Rome ?

The job of Bishop of Rome – i.e. the Pope – is considered to be a life-long commitment. I mean, it usually was. There have been 266 popes since 32 A.D. (according to http://oce.catholic.com/…): almost all popes have served until their death. But that does not mean that they were in the job for long… One can easily extract the data from the website,

> L2=scan("http://oce.catholic.com/index.php?title=List_of_Popes",what="character")
> index=which(L2=="</td><td>Reigned")
> X=L2[index+1]
> Y=strsplit(X,split="-")

But one should work a little bit because sometimes, there are inconsistencies, e.g. 911-913 and then 913-14, so we need some more lines. Further, we can extract from this file the years popes started to reign, the year it ended, and the length, using those functions

> diffyears=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=1}
+ if(length(x)==2){s=diff(as.numeric(x))}
+ return(s)}
> whichyearsbeg=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=as.numeric(x)}
+ if(length(x)==2){s=as.numeric(x)[1]}
+ return(s)}
> whichyearsend=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=as.numeric(x)}
+ if(length(x)==2){s=as.numeric(x)[2]}
+ return(s)}

On our file, we have

> Years=unlist(lapply(Y,whichyearsbeg))
> YearsB=c(Years[1:91],752,Years[92:length(Years)])
> YearsB[187]=1276
> Years=unlist(lapply(Y,whichyearsend))
> YearsE=c(Years[1:91],752,Years[92:length(Years)])
> YearsE[187]=1276
> YearsE[266]=2013
> YearsE[122]=914
> W=unlist(lapply(Y,diffyears))
> W=c(W[1:91],1,W[92:length(W)])
> W[W==-899]=1
> which(is.na(W))
[1] 187 266
> W[187]=1
> W[266]=2013-2005

If we plot it, we have the following graph,

> plot(YearsB,W,type="h")

and if we look at the average length, we have the following graph,

> n=200
> YEARS = seq(0,2000,length=n)
> Z=rep(NA,n)
> for(i in 2:(n-1)){
+ index=which((YearsB>YEARS[i]-50)&(YearsE<YEARS[i]+50))
+ Z[i] = mean(W[index])}
> plot(YEARS,Z,type="l",ylim=c(0,30))
> n=50
> YEARS = seq(0,2000,length=n)
> Z=rep(NA,n)
> for(i in 2:(n-1)){
+ index=which((YearsB>YEARS[i]-50)&(YearsE<YEARS[i]+50))
+ Z[i] = mean(W[index])}
> lines(YEARS,Z,type="l",col="grey")

which does not reflect mortality improvements that have been observed over two millenniums. It might related to the fact that the average age at time of election has  increased over time (for instance, Benedict XVI was elected at 78 – one of the oldest to be elected). Actually, serving a bit more than 7 years is almost the median,

> mean(W>=7.5)
[1] 0.424812

(42% of the Popes did stay at least 7 years in charge) or we can look at the histogram,

> hist(W,breaks=0:35)

Unfortunately, I could not find more detailed database (including the years of birth for instance) to start a life-table of Popes.

# Les députés sont-ils à l’image de la population

Beaucoup de choses ont été écrites sur le fait que les députés ne sont pas vraiment le reflet de la population, que ce soit en terme de profession, de sexe, d’origine, d’age, etc. La liste pourrait être longue. Il y a plusieurs mois, j’avais commencé à regarder le profil des députés, par age. En effet, le site http://assemblee-nationale.fr/ permet d’accéder à des données sur tous les députés, depuis la Révolution. Y compris leur date de naissance. En croisant ces données avec des données de population, par exemple via http://www.mortality.org/, on peut comparer la répartition des ages des députés, avec la répartition des ages de la population.

Pour les amateurs, le code pour récupérer les données (ou au moins les dates de naissance des députés) ressemble à

N=2002
URL=paste("http://www.assemblee-nationale.fr/
regle_nom=est&Nom=&departement=&
choixdate=intervalle&D%C3%A9butMin=01%2F01%2F",
N,"&FinMin=31%2F12%2F",N,"&Dateau=&legislature=",
s,"&choixordre=chrono&Rechercher=
Lancer+la+recherche",sep="")
HTML=scan(URL,what="character")

k=which(HTML=="class=\"titre\">Né")
vHTML=HTML[k:length(HTML)]
vk=which(substr(vHTML,1,7)==">&nbsp;")
liste=vHTML[vk]
naissance=liste[seq(1,length(liste),by=2)]
NAISSANCE=as.Date(substr(naissance,8,17),
"%d/%m/%Y")

Maintenant, pour être tout à fait honnête, je ne suis pas certain de ce qui est vraiment renvoyé, et j’ai des doutes que cela correspondent réellement à la requête faite. En effet, même si je demande à avoir la liste des députés après l’élection, j’ai trop de monde… mais peut-être est-ce du aux décès éventuels, et il est possible que l’ensemble des députés qui ont siégé pendant la mandature apparaissent dans le résultat de la requête.

Sur la figure suivante on voit, sur plusieurs élections depuis plus de 100 ans, comment les deux distributions se déforment, avec en rouge la distribution de l’age des députés, et en bleu, la distribution de la population française, dans son ensemble (population de plus de 18 ans)

Si on veut tout suivre sur un graphique, au lieu de se regarder une animation, on peut représenter les différents quantiles (10%, 25%, 75% et 90%, retenus sur la population de plus de 18 ans, et l’age médian, au centre), avec la population française l’année de l’élection,

et l’ensemble des élus au parlement,

Si on veut faciliter la comparaison, on peut se contenter de visualiser l’évolution des ages moyens,

ou encore, du ratio (en % de différence) entre l’age moyen des députés, et celui de l’ensemble de la population.

Sur ce graphique, on voit que depuis 30 ans, l’age moyen des députés croit plus vite que celui de la population: la population français vieilli, mais moins que ses députés… La gérontocratie perdure donc en France. En espérant que cela ne débouche pas sur le clash générationnel que l’on semble observer ces temps-ci au Québec…

# Date of death, birthday and Elvis Presley

10 days ago, a study published on http://www.annalsofepidemiology.org/ mentioned that “Death has a preference for birthdays” (as claimed in the title). The conclusion of the paper is that, in general, birthdays do not evoke a postponement mechanism but appear to end up in a lethal way more frequently than expected (“anniversary reaction”). Well, this is not new, and several previous articles have mentioned that point, e.g. Angermeyer et al. (1987).

I found the idea interesting since in demography, there is a large literature trying to extrapolate death rates from discrete to continuous time. Extrapolation are usually extremely smooth. But none of them integrate that aspect of mortality precisely on the birthday. The problem is that it is rather difficult to say something since datasets with individual observations are rare, online.

But yesterday, @coulmont sent me a tweet mentioning a website. I do not know if this is legal (even if some explanations are given), but I will mention courtesy of http://ssdmf.info/. It is a so-called Social Security Death Master File, containing individual informations about deaths in the US, as well as geographic information (as described on http://www.ssa.gov/), for people having a social security number.

With R, it is possible to work on those files (even they are huge, with tens of millions observations). For instance, we can check who is inside.

> elvis=scan("ssdm2",skip=22371720,n=1,what="character",sep=",")
> elvis
[1] " 409522002PRESLEY         ELVIS     0800197701081935  "

If you believe that Elvis is dead, you might agree that this database can be accurate (or at least, not too bad). And further, we can see here how to read the result: Elvis was born on January 8, 1935 (8 last digits), and died on August 16, 1977 (8 digits before). Obviously here, there are some problems with the dataset (we do not have the day of the death of Elvis). So here, we remove all the observations that do not give us proper dates. Then, the idea is to assume that the person died in 2000 (or any year since the point is to focus on days and months). Then, we count the number of days between the day of death and the birthday in 2001 (that would have been after) and the one in 2000 (that was either before or after the death), so that we can derive the number of days after the birthday,

dates=substr(base,66,81)
death=as.Date(substr(dates,1,8),"%m%d%Y")
birth=as.Date(substr(dates,9,16),"%m%d%Y")
indice=is.na(death)|is.na(birth)
mean(indice)
mdeath=substr(dates,1,2)
ddeath=substr(dates,3,4)
mbirth=substr(dates,9,10)
dbirth=substr(dates,11,12)
indice=which(ddeath!="00")
birth1=as.Date(paste(mbirth[indice],
dbirth[indice],"2000",sep=""),"%m%d%Y")
birth2=as.Date(paste(mbirth[indice],
dbirth[indice],"2001",sep=""),"%m%d%Y")
death=as.Date(paste(mdeath[indice],ddeath[indice],
"2000",sep=""),"%m%d%Y")
k=length(indice)
diffday=cbind((as.numeric(death-birth1))[1:k],
(as.numeric(death-birth2))[1:k])
DIFF=apply(diffday,1,function(x) {min(x[x>=0])})

What we have here is the number of days following the previous birthday. If we look at the distribution of that number of days, we obtain

counts=table(DIFF)
plot(as.numeric(names(counts)),
as.numeric(counts))
counts["0"]/(mean(counts[100:200]))
> counts["0"]/(mean(counts[100:200]))
0
1.121261

Thus, the death excess on the day of birth was around 12%, which is rather close to the one obtained from the Swiss mortality statistics 1969–2008 (in Ajdacic-Gross et al. (2012)). Note that here, we just play with a small subset of the entire dataset,

That database is probably extremely interesting, except that it suffers a huge selection bias, since only dead people are in that database. So it might be useless if we wish to study life expectancy of people named Bill versus people named Georges (that was something I wanted to investigate initially). But we’ll see what else we can do with it (since Ewen have been able to write some code to go through that huge dataset).

# La pyramide des âges, au Canada

L’autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la “pyramide” des âges (évoquée par exemple surhttp://www.lapresse.ca/). Bon, j’ai des données assez proches, tirées dehttp://www.mortality.org/).

> pop=read.table(
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,] On peut faire assez facilement des pyramides des ages, par année, avec le code suivant, > library(plotrix) > seuils=seq(0,110,by=10) > pop$tranche=cut(pop$Age,seuils, right = FALSE) > an=1955 > base=pop[pop$Year==an,]
> women=aggregate(base$Female, + list(base$tranche),sum)[,2]
> men=aggregate(base$Female, + list(base$tranche),sum)[,2]
> nom=as.character(unique(pop$tranche)) > pyramid.plot(men/sum(men)*100, + women/sum(women)*100,labels=nom,gap=2, + lxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red"), + rxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red")) Au lieu de faire des barres horizontales, on peut cumuler, par tranche d’age, et regarder l’évolution dans le temps des proportions “moins de 20 ans” (bleu), “entre 20 et 60 ans” (mauve) et “plus de 60 ans” (rouge). > seuils=c(0,20,60,110) > pop$tranche2=cut(pop$Age,seuils, right = FALSE) > YEAR=1921:2008 > totaux=matrix(NA,length(YEAR),3) > for(i in 1:length(YEAR)){ + base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total, + list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu’après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l’isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d’interprétation). A suivre…

# Longevity and mortality dynamics with R

Following the previous post on life contingencies and actuarial models in life insurance, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The second part of the talk (on Actuarial models with R) will be dedicated to longevity and mortality. A complete set of slides can be downloaded from the blog, but again, only some part will be presented.

As mentioned earlier, the codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides above can be downloaded from here, and datasets are the following,

> DEATH <- read.table(
+ "http://freakonometrics.free.fr/Deces-France.txt",
+ "http://freakonometrics.free.fr/Exposures-France.txt",
+ header=TRUE,skip=2)

For additional resources, I will use Rob Hyndman‘s package on demography, Heather Turner and David Firth’s package on generalized nonlinear models (e.g. the slides of the short course Heather gave in Rennes at the UseR! conference in 2009), as well as functions developed by JPMorgan’s LifeMetrics (functions are  fully documented in the LifeMetrics Technical Document). All those functions can be obtained using

> library(demography)
> library(gnm)
> source("http://freakonometrics.free.fr/fitModels.R")

# Births and week-ends, in France

This week, I have seen on the internet (sorry, I cannot find proper references) the graph produced here on the right: which birthday is most likely ? The fact that I have no further information is important, since I do not know in which country such a graph was obtained. At least, I know it should not be France…

In France, I have already mentioned that there is a strong week-end effect: nowadays, there is 25% less deliveries during week-ends than during the week. Calot (1981) observed already that there were less deliveries on Sundays. This has been confirmed more recently, e.g. in http://www.lepoint.fr/ or http://www.prepabl.fr/, with a significant difference between week days, and week-ends. Here  is the number of birth per day, over 40 years, with in blue the average trend during the week, and in red, during week-ends,

naissance=read.table(
"http://freakonometrics.free.fr/naissanceFR2.txt")
attach(naissance)
date=as.Date(date)
plot(date, nbre,cex=.5)
t2=as.POSIXlt(date)
jour=t2$wday X=naissance$date
Y=naissance$nbre J=jour df=data.frame(X,Y,J) library(splines) regs=lm(Y~bs(X,df=20),data=df[jour%in%c(0,6),]) Yp=predict(regs,newdata=df) lines(X,Yp,col="red",lwd=3) regs=lm(Y~bs(X,df=20),data=df[jour%in%1:5,]) Yp=predict(regs,newdata=df) lines(X,Yp,col="blue",lwd=3) If we look at the evolution of the ratio week-ends over weeks days, we have the following graph t2=as.POSIXlt(date) jour=t2$wday
jour=jour[1:(1982*7)]
nbre2=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbre2[1:5+7*(i-1)]=nbre[1:5+7*(i-1)]*taux
nbre2[6:7+7*(i-1)]=nbre[6:7+7*(i-1)]
nbre2[1:7+7*(i-1)]=
mean(nbre[1:7+7*(i-1)])/mean(nbre2[1:7+7*(i-1)])*
nbre2[1:7+7*(i-1)]
}
nbretaux=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbretaux[1:7+7*(i-1)]=taux
}
plot(date[1:length(nbre2)],nbretaux)
X= date[1:length(nbre2)]
Y=nbretaux
library(splines)
reg=lm(Y~bs(X,df=20))
Yp=predict(reg)
lines(X,Yp,col="red",lwd=3)

In the beginning of the 70’s, during week-ends, there were 5% less deliveries, but 25% less around 2000. It is then possible to produce the same kind of graphs as the one above, per year of birth. And here, we clearly observe the importance of the week end effect (maybe also because of color choice)

naissance=read.csv(
"http://freakonometrics.free.fr/naissanceFR.csv",
sep=";")
M=as.matrix(naissance[,3:ncol(naissance)])
BIRTH=as.vector(t(M))
YEAR=rep(1968:2005,each=12*31)
MONTH=rep(rep(1:12,each=31),38)
DAY=rep(1:31,12*38)
X=NA
for(y in 1968:2005){
sbase=base[YEAR==y,]
X=c(X,sbase$BIRTH/sum(sbase$BIRTH,
na.rm=TRUE))
}
base=data.frame(YEAR,MONTH,DAY,
BIRTH,BIRTHDAYPROB=X[-1])

m1=min(base$BIRTHDAYPROB,na.rm=TRUE) m2=max(base$BIRTHDAYPROB,na.rm=TRUE)
y=1980
colr=rev(heat.colors(100))
sbase=base[YEAR==y,]
plot(0:1,0:1,col="white",xlim=c(-1,12),
ylim=c(-31,1),axes=FALSE,xlab=
paste("Naissance en",y,sep=" "),ylab="")
for(x in 1:nrow(sbase)){
a=sbase$MONTH[x];b=sbase$DAY[x]
polygon(c(a-.9,a-.9,a-.1,a-.1),-c(b-.9,b-.1,
b-.1,b-.9),col=colr[(sbase\$BIRTHDAYPROB[x]-m1)/
(m2-m1)*100],border=NA)
}
text((1:12)-.5,.5,c("J","F","M","A","M","J","J",
"A","S","O","N","D"),cex=.7)
text(-.5,-(1:31)+.5,1:31,cex=.7)