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))