Tag Archives: démographie

Data XXL et Démographie Participative

Petite virée parisienne pour trois jours. Je serais à Science Po mercredi pour travailler un peu, et à l’INED jeudi et vendredi pour la conférence Sciences XXL – Ce que l’abondance et la diversité des données numériques font aux sciences sociales. Jeudi, je serais discutant pour la session “outils”, et vendredi après-midi, nous allons clôture la conférence avec un rapide retour sur un travail joint avec Olivier et Ewen, sur ‘Les données collaboratives généalogiques au service de la démographie historique’. Le travail est très vaste, et nous présenterons quelques résultats de réconciliation entre les tables de mortalité historique de l’INED, et celle que nous avons reconstituées sur des données démographiques…. Par exemple avec les naissances de 1850

et celles de 1899

La suite…. vendredi (ou plus tard, on n’a pas fini de jouer avec ces données).

 

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(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> 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"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

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…

Combien de temps profite-t-on de ses grands parents ?

Ce Hier matin, je suis tombé un peu par hasard sur deux graphiques de l’INSEE (en France) avec l’age moyen des mères à l’accouchement, en fonction de rang de naissance de l’enfant, avec tout d’abord 1905-1965,

puis 1960-2000

Ces graphiques sont passionnants en soi – comme en ont témoigné pas mal de followers sur Twitter – mais ils m’ont fait m’interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n’a – malheureusement – pas accès au données complètes sur le site, mais on peut trouver d’autres donnéesintéressantes (en l’occurrence l’age moyen à la naissance).

> agenaissance=read.table("http://freakonometrics.blog.free.fr/
public/data/agenaissance.csv",header=TRUE,sep=",")
> agenaissance$Age=as.character(agenaissance$AGE)
> agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+
+ as.numeric(substr(agenaissance$Age,4,4))/10
> plot(agenaissance$ANNEE+.5,agenaissance$AGE,
+ type="l",lwd=2,col="blue")

Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus,

On peut alors aller en cran plus loin, en se demandant non pas quel était l’âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)

> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)-
+ agenaissance$AGE
> w=(trunc(agenaissance$NAIS.MERE-.5))
> rownames(agenaissance)=agenaissance$ANNEE
> a1=agenaissance[as.character(w),]$NAIS.MERE
> a2=agenaissance[as.character(w+1),]$NAIS.MERE
> p=agenaissance$NAIS.MERE-(w+.5)
> agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2
> agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5-
+ agenaissance$NAIS.GRD.MERE
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE
2000  2000 30.3 30,3     1970.2       1942.87        57.63
2001  2001 30.4 30,4     1971.1       1943.80        57.70
2002  2002 30.4 30,4     1972.1       1944.92        57.58
2003  2003 30.5 30,5     1973.0       1945.95        57.55
2004  2004 30.5 30,5     1974.0       1947.05        57.45
2005  2005 30.6 30,6     1974.9       1948.04        57.46
> plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE,
+ type="l",lwd=2,col="red")

Là encore, on peut visualiser l’âge de la grand-mère maternelle à la naissance

A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d’espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.

> Deces <- 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)
> Deces$Age <- as.numeric(as.character(Deces$Age))
> Deces$Age[is.na(Deces$Age)] <- 110
> Expo$Age <- as.numeric(as.character(Expo$Age))
> Expo$Age[is.na(Expo$Age)] <- 110
>  library(forecast)
>  library(demography)
>  YEAR <- unique(Deces$Year);nC=length(YEAR)
>  AGE  <- unique(Deces$Age);nL=length(AGE)
>  MUF  <- matrix(Deces$Female/Expo$Female,nL,nC)
>  POPF <- matrix(Expo$Female,nL,nC)
>  BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE,
+ years=YEAR, type="mortality",
+  label="France", name="Femmes", lambda=1)
> LCF <- lca(BASEF)
> LCFf<-forecast(LCF,h=100)
> A <- LCF$ax
> B <- LCF$bx
> K1 <- LCF$kt
> K2 <- K1[length(K1)]+LCFf$kt.f$mean
> K <- c(K1,K2)
> MU <- matrix(NA,length(A),length(K))
> for(i in 1:length(A)){
+ for(j in 1:length(K)){
+ MU[i,j] <- exp(A[i]+B[i]*K[j]) }}
> esp.vie = function(xentier,T){
+ s <- seq(0,99-xentier-1)
+ MUd <- MU[xentier+1+s,T+s-1898]
+ Pxt <- cumprod(exp(-diag(MUd)))
+ ext <- sum(Pxt)
+ return(ext) }
> EVIE = function(x,T){
+ x1 <- trunc(x)
+ x2 <- x1+1
+ return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) }
> agenaissance$EV=NA
> for(i in 1:100){
+ t <- 2006-i
+ agenaissance$EV[agenaissance$ANNEE==t]=
+ EVIE(x=agenaissance$age.GRD.MERE[
+ agenaissance$ANNEE==t],t) }
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE       EV
2000 30.3 30,3     1970.2       1942.87        57.63 29.13876
2001 30.4 30,4     1971.1       1943.80        57.70 29.17047
2002 30.4 30,4     1972.1       1944.92        57.58 29.39027
2003 30.5 30,5     1973.0       1945.95        57.55 29.52041
2004 30.5 30,5     1974.0       1947.05        57.45 29.72511
2005 30.6 30,6     1974.9       1948.04        57.46 29.80398

Autrement dit, sur la dernière ligne, l’espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d’environ 29.80 ans. On peut alors visualiser non seulement l’âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle,

> plot(agenaissance$ANNEE+.5,agenaissance$EV,
+ type="l",lwd=2,col="purple")

On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s’est stabilisé à une trentaine d’années. On peut aussi continuer, et remonter d’un cran (en refaisant tourner le code avec quelques modifications): on a alors l’âge (moyen) de son arrière grand mère à la naissance,

et la durée de vie (résiduelle) des arrière grand mères

On manque ici de données, mais il semble que l’on profite – en moyenne – environ 5 ans de son arrière grand mère. Maintenant on peut aussi s’interroger sur les limites de cette étude rapide. En particulier, de même qu’il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d’enfants et de petits-enfants a un impact sur la durée de vie résiduelle d’une personne (ou si on peut supposer l’indépendance comme on l’a fait ici).

Modélisation des tables de mortalité

En cours d’actuariat, nous commencerons à parler des tables de mortalité prospectives.

Côté théorie, je peux renvoyer vers quelques slides de Pierre Devolder, ici. Je renvoie ici pour un court papier de vulgarisation sur le sujet. Sur l’ajustement de lois de mortalité prospectives, je renvoie au très beau mémoire en ligne ici (ou  pour des choses très proches en anglais). Pour l’approche statique, je renvoie également ici pour une méthode simple d’ajustement de loi de Makeham. Pour les aspects dynamique, il y a des documents techniques ici ou .

Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et  (pour les décès). Il faut retravailler un peu les données, à cause d’une valeur “110+”, mais rien de plus simple,

> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Exposures_1x1.txt",header=TRUE)
> DEATH$Age=as.numeric(as.character(DEATH$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> DEATH$Age[is.na(DEATH$Age)]=110
> EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))
Message d'avis :
NAs introduits lors de la conversion automatique 
> EXPOSURE$Age[is.na(EXPOSURE$Age)]=110

Coté code, je renvoie à un précédant billet sur le sujet, ici, et  pour le package (non officiel) de démographie développé par Rob Hyndman.

Enfin, à la fin du cours, on parlera un peu de modèles de durée. J’essayerai de taper quelques choses dans les jours à venir, mais côté TD, on utilisera ce que fait Frédéric Planchet (qui connaît fort bien le sujet), et qui se trouve en ligne ici. Les données sont en ligne ici, avec un peu de code R pour retraiter les donnéesici et , et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier ().

Les codes R utilisés en TD sont les suivants (avec en plus les graphiques obtenus),

> D=DEATH[DEATH$Year==an,]
> E=EXPOSURE[EXPOSURE$Year==an,]
> MU = D[,3:5]/E[,3:5]
> plot(0:110,log(MU[,1]),type="l",col="red")
> lines(0:110,log(MU[,2]),col="blue")

Avec ci-dessus les log des taux de décès instantannés pour l’année 1986.

> PH=matrix(NA,111,111) # ligne X, colonne H
> PF=PH
> for(x in 0:110){
+ PH[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,2]))
+ PF[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,1]))
+ }
> x=0
> plot(1:111,PH[x+1,],ylim=c(0,1),type="l",col="blue")
> lines(1:111,PF[x+1,],col="red")

pour la fonction de survie à la naissance (ci-dessus) ou à 50 ans (ci-dessous),

> somme=function(x){sum(x,na.rm=TRUE)}
> EH=EF=rep(NA,111)
> EH=apply(PH,1,somme)
> EF=apply(PF,1,somme)
> plot(0:110,EH,type="l",col="blue")
> lines(0:110,EF,col="red")

Ce qui correspond aux espérances de vie résiduelles.

> for(k in 1:3){
+ Z=MU[,k]
+ Z[is.nan(Z)==TRUE]=NA
+ MU[,k]=Z
+ }
> AGE=unique(DEATH$Age)
> ANNEE=unique(DEATH$Year)
> MUF=matrix(MU[,1],length(AGE),length(ANNEE))
> MUH=matrix(MU[,2],length(AGE),length(ANNEE))
> persp(AGE[1:100],ANNEE,log(MUH[1:100,]),
+ theta=-30,col="light green",shade=TRUE)
> persp(AGE[1:100],ANNEE,log(MUH[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUH))

On a retourné le graph en mettant les années à l’envers, mais on peut aussi regarder la surface pour les femmes, ainsi que l’allure des courbes de niveau

> persp(AGE[1:100],ANNEE,log(MUF[1:100,107:1]),
+ theta=-30,col="light green",shade=TRUE)
> image(AGE,ANNEE,log(MUF))

Mais nous n’avons ici qu’une lecture statique des tables, et nous ne suivons pas vraiment un individu: en même temps que la personne vieilli, le temps passe…. Considérons une personne née en 1900, et regardons (ex post) l’évolution de sa probabilité de décès (en regardant non pas la probabilité de décéder à 18 d’une personne de 18 ans en 1900, mais en tenant compte du fait que la personne a fêté ses 18 ans en 1918).

> MUFP=MUHP=rep(NA,111)    # MU PROSPECTIF
> for(i in 1:111){
+     MUHP[i]=MUH[i,colonne+i]
+     MUFP[i]=MUF[i,colonne+i]
+ }
> plot(0:110,log(MUHP),col="purple",type="l")
> D=DEATH[DEATH$Year==annee.naissance,]
> E=EXPOSURE[EXPOSURE$Year==annee.naissance,]
> MU = D[,3:5]/E[,3:5]
> lines(0:110,log(MU[,2]),col="blue")

On peut aussi l’ajustement par un modèle à la Lee-Carter (via la fonction fournie par LifeMetrics)

> source("https://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
> DEATH=DEATH[DEATH$Age<90,] # virer age>90
> EXPOSURE=EXPOSURE[EXPOSURE$Age<90,] # virer age>90
> XV=unique(DEATH$Age)
> YV=unique(DEATH$Year)
> ETF=t(matrix(EXPOSURE[,3],length(XV),length(YV)))
> DTF=t(matrix(DEATH[,3],length(XV),length(YV)))
> ETH=t(matrix(EXPOSURE[,4],length(XV),length(YV)))
> DTH=t(matrix(DEATH[,4],length(XV),length(YV)))
> WA=matrix(1,length(YV),length(XV))
> LCF = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)
-1166100 -> -263940.2  ->-204175.2 
-188375.9 -> -185822.4  ->-184403.5 
-183223.3 -> -183069.8  ->-182350.3 
-181701.1 -> -181600.6  ->-181145.5 
-180729.6 -> -180671.5  ->-180371.8 
-180101.9 -> -180068.6  ->-179871.2 
> LCH = fit701(xv=XV,yv=YV,etx=ETH,dtx=DTH,wa=WA)
-2564403 -> -1028776  ->-953629.9 
-881334.2 -> -868279.3  ->-834413.6 
-804946.8 -> -802817.2  ->-786343 
-772775 -> -772007.5  ->-764501.7 
-758653.8 -> -758399  ->-755212 
-752931 -> -752851  ->-751620.2 
-750811.1 -> -750786.7  ->-750351.4 
> plot(LCF$x,LCF$beta1,type="l",col="red")
> lines(LCH$x,LCH$beta1,col="blue")

> plot(LCF$x,LCF$beta2,type="l",col="red")
> lines(LCH$x,LCH$beta2,col="blue")

> plot(LCF$y,LCF$kappa2,type="l",col="red")
> lines(LCH$y,LCH$kappa2,col="blue")

> MUF      = DTF/ETF
> MUH      = DTH/ETH
> MUFMODEL = MUF; MUHMODEL = MUH
> for(i in 1:nrow(MUF)){
+ for(j in 1:ncol(MUF)){
+ MUFMODEL[i,j] = exp(LCF$beta1[j]
+ +LCF$beta2[j]*LCF$kappa2[i])
+ MUHMODEL[i,j] = exp(LCH$beta1[j]
+ +LCH$beta2[j]*LCH$kappa2[i])
+ }}

soit, pour avoir un dessin,

> persp(LCH$y,LCH$x,log(MUH),
+ theta=45,col="light green",shade=TRUE)
> persp(LCH$y,LCH$x,log(MUHMODEL),
+ theta=45,col="light blue",shade=TRUE)

On peut ainsi comparer les taux bruts (en vert) avec les taux estimés (en bleu)

mais aussi regarder  se trouvent les erreurs le plus grandes,