Poisson ou Binomiale ? Application sur un problème de démographie. Les slides sont en ligne (slides 21) et la vidéo aussi (slides 21)
Tag Archives: espérance de vie
Prédire l’espérance de vie (à la naissance)
Pour finir les exemples de modèles de prévision, on a parlé de calculs d’espérance de vie (à la naissance). Le point de départ des modèles démographiques, ce sont les tables de mortalité, dont je parlais dans le précédant billet. Mais cette vision est très biaisée, car elle suppose qu’il n’y aura pas d’amélioration des conditions de vie. Pour faire les choses proprement, on utilise des données plus complètes, avec des nombres de décès en fonction de l’âge x, mais aussi de la date t.
url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt" download.file(url,"FRD.txt") url="http://freakonometrics.free.fr/FranceExposures_1x1.txt" download.file(url,"FRE.txt") DECES=read.table("FRD.txt",skip = 3,header=TRUE) EXPOS=read.table("FRE.txt",skip = 3,header=TRUE) |
On notera D_{x,t} le nombre de décès observé, et E_{x,t} l’exposition. De telle sorte que la probabilité de décéder dans l’année, pour une personne d’âge x à la date t est q_{x,t}=D_{x,t}/E_{x,t}. Classiquement, ces données sont stockées dans une matrice pour la visualisation, et dans une base de données pour les régressions
annees = 1816:2015 ages = 0:110 MDF=matrix(DECES$Female,length(ages),length(annees)) MDH=matrix(DECES$Male,length(ages),length(annees)) MEF=matrix(EXPOS$Female,length(ages),length(annees)) MEH=matrix(EXPOS$Male,length(ages),length(annees)) MEF[MEF<=0]=1 MEH[MEH<=0]=1 QF=MDF/MEF QH=MDH/MEH QF[QF<=0]=NA QH[QH<=0]=NA |
(il faut bricoler un peu pour éviter les problèmes de valeurs nulles car (i) on fait des ratios (ii) on passera ensuite au logarithme). On peut visualiser cette surface, en fonction de x et t.
persp(log(QF)) |
ou (plus visuel)
library(rgl) persp3d(ages,annees,log(QH),col="light blue") |
Pour modèliser l’évolution de q_{x,t} on peut s’inspirer du modèle de Lee & Carter (1992) qui suppose que \log(q_{x,t})=A_x+B_x\cdot K_t. Les \mathbf{A}=(A_0,A_1,\cdots,A_{110}) c’est en quelques sortes l’évolution moyenne de \log(q_{x,t}). \mathbf{K}=(K_{1816},K_{1817},\cdots,K_{2015}) va permettre de comprendre les gains de conditions de vie, la baisse de la probabilité de décéder dans l’année. Ces gains ne sont pas uniformes, et on utilise \mathbf{B}=(B_0,B_1,\cdots,B_{110}) pour faire dépendre les gains de l’âge.
Pour estimer les paramères, \mathbf{A}, \mathbf{B} et \mathbf{K}, on peut être tentés d’utiliser un modèle binomial, au sens où D_{x,t}\sim \mathcal{B}(E_{x,t},q_{x,t}), qui est le modèle de base en assurance-vie. Autrement dit, ici D_{x,t}\sim \mathcal{B}(E_{x,t},\exp[A_x+B_x\cdot K_t]).
Une autre piste est d’utiliser la loi des petits nombres, qui dit que si la probabilité est faible (et c’est le cas quand on s’intéresse à la probabilité de mourir dans l’année), la loi binomiale peut être approchée par une loi de Poisson. Autrement dit, D_{x,t}\sim \mathcal{P}(E_{x,t}\cdot q_{x,t}), i.e. D_{x,t}\sim \mathcal{P}(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t]). On reconnait ici une régression de Poisson, avec comme variables explicatives l’âge x, l’année t, et l’exposition apparait comme variable offset. Le seul soucis est que ce n’est pas une régression linéaire. On a ici un modèle non linéaire, puisque \mathbb{E}[D_{x,t}]=(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t]). On a à la fois une somme et un produit. Bref, on va partir vers un modèle nonlinéaire généralisé.
base=data.frame(DH=DECES$Male,DF=DECES$Female, EH=EXPOS$Male,EF=EXPOS$Female, age=rep(ages,length(annees)), annee=rep(annees,each=length(ages))) idx=which((base$DH<=0)|(base$DF<=0)|(base$EH<=0)|(base$EF<=0)) base=base[-idx,] library(gnm) reg=gnm( DH ~ offset(log(EH)) + as.factor(age) + Mult(as.factor(age),as.factor(annee)), family = poisson(link="log"), data=base) |
On a nos estimateurs \hat{\mathbf{A}}, \hat{\mathbf{B}} et \hat{\mathbf{K}}.
Ax=reg$coefficients[2:111] Bx=reg$coefficients[112:222] Kt=reg$coefficients[223:length(reg$coefficients)] |
On va passer sous silence la non identifiabilité du modèle. On peut représenter nos trois jeux de coefficients. Tout d’abord les \hat{\mathbf{A}} qui nous donnent l’évolution moyenne,
plot(ages[-1],Ax) |
On peut aussi tracer \hat{\mathbf{K}} qui donne l’évolution temporelle.
plot(annees,-Kt) |
Encore une fois, le modèle n’est pas identifiable. Si (\hat{\mathbf{A}},\hat{\mathbf{B}},\hat{\mathbf{K}}) est solution, il en est de même pour (\hat{\mathbf{A}},-\hat{\mathbf{B}},-\hat{\mathbf{K}}). Bref, le sens d’évolution n’a pas de sens. On peut représenter -\hat{\mathbf{K}} qui présente l’avantage d’être décroissant (et de décrire l’amélioration des conditions de vie). Enfin, traçons aussi -\hat{\mathbf{B}}
plot(ages,-Bx) |
La difficulté est que pour prévoir l’espérance de vie, on a besoin de calculer q_{t,x} pour des valeurs importantes de t (non encore observée). Par exemple, on pourrait avoir envie de connaître q_{50,2020} (pour quelqu’un né en 1970). On a envie d’utiliser q_{50,2020}=\exp(\hat{A}_{50}+\hat{B}_{50}\hat{K}_{2020}). Le soucis c’est que \hat{K}_{2020} ne fait pas partie des grandeurs estimées \hat{\mathbf{K}}.
L’idée est de construire un prédicteur \tilde{K}_{2020} à partir de \{\hat{K}_{1816},\cdots,\hat{K}_{2015}\}. Classiquement (et c’est l’idée originale de Lee & Carter (1992)) on peut tenter un modèle ARIMA. Ici, la série est un peu trop régulière. On peut tenter soit un modèle exponentiel, soit un modèle linéaire (sur la série originale \hat{\mathbf{K}}, au delà de 1950)
idx=which(annees>=1950) regk=lm(log(Kt[idx])~annees[idx]) futur=2016:2125 prevKt=exp(regk$coefficients[1]+regk$coefficients[2]*futur) plot(annees,Kt,xlim=c(1815,2125),ylim=range(c(Kt,prevKt))) points(futur,prevKt,col="red") regk2=lm(Kt[idx]~annees[idx]) prevKt2=(regk2$coefficients[1]+regk2$coefficients[2]*futur) points(futur,prevKt2,col="blue") |
On peut alors construire la série des prédictions sur les données passées, \hat{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\hat{K}_{t}], et sur les données futures \tilde{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\tilde{K}_{t}].
Si on reste sur les données passées, c’est assez simple, avec ici les probabilités de décès pour 1880
BASE=data.frame( x=base$age, t=base$annee, indice_x=1+base$age, indice_t=base$annee-1815) BASE$Ax=Ax[BASE$indice_x] BASE$Bx=Bx[BASE$indice_x] BASE$Kt=Kt[BASE$indice_t] cste=reg$coefficients[1] BASE$Qpred=exp(cste+BASE$Ax+BASE$Bx*BASE$Kt) plot(BASE$x[BASE$t==1880],BASE$Qpred[BASE$t==1880], log="y") |
On peut faire la même chose pour le futur (ici pour l’année 2050) avec les deux modèles
BASE2=data.frame( x=rep(ages,length(futur)), t=rep(futur,each=length(ages)), indice_x=1+rep(ages,length(futur)), indice_t=rep(futur,each=length(ages))-2015) BASE2$Ax=Ax[BASE2$indice_x] BASE2$Bx=Bx[BASE2$indice_x] BASE2$Kt1=prevKt[BASE2$indice_t] BASE2$Kt2=prevKt2[BASE2$indice_t] cste=reg$coefficients[1] BASE2$Qpred1=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt1) BASE2$Qpred2=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt2) head(BASE2) plot(BASE2$x[BASE2$t==2050],BASE2$Qpred1[BASE2$t== 2050],log="y") |
pour la prévision exponentielle, et
pour la prévision linéaire On va ensuite bricoler pour créer le vecteur \mathbf{q} qui contiendra soit \hat{q}_{h,t+h}, soit \tilde{q}_{h,t+h}. En effet, pour faire un calcul d’espérance de vie, on ne regarde pas à un instant donnée, mais on suit une “cohorte” dans le temps. Par exemple, pour une personne née en 1968, on a la probabilité de décès (dans l’année) suivante
annee = 1968 vt=annee+ages vx=ages vq=rep(NA,length(vx)) sbase=data.frame(x=vx,t=vt) for(i in 1:nrow(sbase)){ if(sbase$t[i] <= 2015) {vq[i]=BASE[(BASE$x==sbase$x[i])&(BASE$t==sbase$t[i]),"Qpred"]} if(sbase$t[i] > 2015) {vq[i]=BASE2[(BASE2$x==sbase$x[i])&(BASE2$t==sbase$t[i]),"Qpred2"]} } plot(ages,vq,log="y") abline(v=2015-1968+.5,lty=2) |
avec à gauche les valeurs estimées par notre modèle (on aurait pu utiliser les valeurs brutes), et à droite, les valeurs prédites (puisqu’on parle de probabilités de décès en 2030, entre autres).
Pour calculer l’espérance de vie à la naissance, on utilise le code suivant
sum(cumprod(exp(-vq[1:110]))) [1] 77.62047 |
On peut alors faire une petite fonction, afin de visualiser l’évolution de cette espérance de vie
Ex=function(annee = 1968,past="obs",fut="2"){ vt=annee+ages vx=ages sbase=data.frame(x=vx,t=vt,vq=NA) for(i in 1:nrow(sbase)){ if(sbase$t[i] <= 2015) {sbase$vq[i]=BASE[(BASE$x==sbase$x[i])&(BASE$t==sbase$t[i]),paste0("Q",past)]} if(sbase$t[i] > 2015) {sbase$vq[i]=BASE2[(BASE2$x==sbase$x[i])&(BASE2$t==sbase$t[i]),paste0("Qpred",fut)]} } vP = cumprod(exp(-(sbase$vq[1:110]))) sum(vP)} |
Aussi
ANNEES=1930:2010 E2=Vectorize(Ex)(ANNEES) plot(ANNEES,E2) |
Si on regarde en variation, on voit qu’il y a (en gros) une variation de 0.25 d’année en année (d’où le fameux “on gagne un trimestre d’espérance de vie tous les ans”)
plot(ANNEES[-1],diff(E2),ylim=c(0,1)) |
Comme on plafonne ici les âges (à 110 ans) on comprend que ce gain ne peut que dimunuer, ce qu’on observe ici.
Par contre, si on prend une prévision qui conserve la variation exponentielle de K_t, on obtient
E1=Vectorize(function(a) Ex(a,"obs","1"))(ANNEES) plot(ANNEES,E1) |
ce qui n’est pas réaliste. Encore moins en considérant la variation de cette courbe (beaucoup trop optimiste)
plot(ANNEES[-1],diff(E1),ylim=c(0,1)) |
Les Espérances de Vie (à la naissance)
Ce soir, nous avions une discussion passionnante avec @FrancoisGourio @msolignac et @Anne__Lavigne au sujet d’un article paru dans Le Monde
Comme souvent, seul un point très précis est retenu d’une longue étude sur la démographie française (ici l’espérance de vie à la naissance) publié par l’INSEE. Au delà du débat sur les chiffres, j’ai pu comprendre ce que l’INSEE mettait derrière le terme “espérance de vie à la naissance” à savoir “la durée de vie moyenne d’une génération fictive qui connaîtrait tout au long de son existence les conditions de mortalité par âge de l’année considérée“. En fait, je n’avais pas compris jusqu’à ce soir ce que cachait cette idée de “génération fictive soumise aux conditions de mortalité de l’année“. En formalisant un peu, ça veut dire que les calculs sont fait de manière transversale, et pas du tout prospective. Car pour faire simple, des données de mortalité, ce sont des données sous la forme suivante,
A différentes dates, disons qu’on fait un recensement. On a alors le nombre de personnes en vie, et décédée à un âge donné. Par exemple, p_{x,t} désigne ici la probabilité qu’une personne d’âge x survive jusqu’à son prochain anniversaire, l’année t. Un peu de calculs permet de voire que
e= \sum_{x=0}^{120} \prod_{n=0}^x p_{n,t}peut être interprété comme une espérance de vie à la naissance. C’est la version discrète de la formule
\mathbb{E}[T]=\int_0^\infty \mathbb{P}(T>t)dt(quand on travaille sur une variable aléatoire positive) pour la somme, et le produit est juste lié au fait que
\mathbb{P}[T>t]=\mathbb{P}[T>t\vert T>t-1]\cdot \mathbb{P}[T>t-1\vert T>t-2] \cdots(avec des mots, la probabilité d’atteindre 60 ans, c’est la probabilité de survivre à son 1er anniversaire, puis à son second, …, puis à son 59ème). Bref, on a toutes ces probabilités de survie dans nos données. Des données, on peut en trouver par exemple sur mortality.org.
> library(demography) > france <- hmd.mx("FRACNP", username, password, "France")
La méthode décrite par l’INSEE ressemble à la suivante: 
On considère toutes les probabilités de survie brute, et on cumule,
> kPx_obs=1-france$rate$total[,"2013"] > sum(cumprod(kPx_obs)) [1] 81.24053
(je ne fais pas ici de distinction Homme/Femme, les plus curieux pourront creuser par eux-mêmes). On note que si en 2013, il y a eu une grosse grippe, qui a touché tout le monde, l’espérance de vie en sera forcement impactée ! On n’a pas ici de correction ici, pour lisser les décès tout au long de l’année, mais c’est possible de le faire.
Une autre stratégie peut être d’ajuster un modèle sur les données de mortalité, par exemple le modèle de Lee & Carter (1992), qui dit que la probabilité qu’une personne d’âge x décède avant son prochain anniversaire, l’année t, est de la forme
1-p_{x,t}=\exp[a_x+b_x\cdot \kappa_t]
On peut ajuster ce modèle et considérer alors
On peut utiliser une version transversale comme auparavant,
> france.LC <- lca(france,adjust="e0",years=1950:2013) > france.fcast <- forecast(france.LC,jumpchoice="actual",h=150) > Ax=france.LC$ax > Bx=france.LC$bx > Kt=c(france.LC$kt,france.fcast$kt$mean) > years=c(1950:2013,2014:2163)
L’espérance de vie à la naissance (telle qu’on l’avait définie auparavant) est alors
> kPx_model=1-exp(Ax+Bx*Kt[years==2013]) > sum(cumprod(kPx_model)) [1] 84.27459
J’avoue être ici surpris par l’écart avec le précédant. Mais pourquoi pas…
Enfin, autre méthode, basée sur un modèle prospectif. L’idée est de noter que le vieillissement d’un individu, dans une cohorte donnée, se fait en diagonale, et pas verticalement, dans le diagramme de Lexis
On utilise ici le fait que quelqu’un né en 2015 autant 40 ans en 2055, et qu’il faut prévoir la probabilité de décéder l’année de ses 40 ans dans 40 ans. En intégrant les gains sur l’espérance de vie, on obtient cette fois l’espérance de vie à la naissance
> kPx_pred=1-exp(Ax+Bx*Kt[years%in 13:(2012+length(Ax))]) > sum(cumprod(kPx_pred)) [1] 90.15392
pour quelqu’un né en 2013. Cette fois, on a un gain vraiment significatif, correspondant à des améliorations des conditions de vie avec le temps, et la baisse de la mortalité.
Je reste toutefois surpris que l’INSEE continue à utiliser une version transversale (et pas du tout prospective) dans ses calculs d’espérance de vie, même si, comme le notait @FrancoisGourio l’espérance de vie à la naissance est une statistique internationale publiée par tous les pays, et donc doit être calculée de façon simple (ou pas trop complexe). Et je renvoie vers un vieux billet pour ceux qui veulent aller plus sur ces calculs prospectifs, et le fameux “tous les ans, on gagne un trimestre d’espérance de vie“. Ou encore un article sur les tables de mortalité, paru dans Risques en 2007.