The Encyclopedia of Educational Research, Measurement, and Evaluation (edited by Bruce B. Frey) is out…
Bruce kindly asked me to write an entry on the central limit theorem in the encyclopedia…
The Encyclopedia of Educational Research, Measurement, and Evaluation (edited by Bruce B. Frey) is out…
Bruce kindly asked me to write an entry on the central limit theorem in the encyclopedia…
Depuis plusieurs semaines (pour ne pas dire plusieurs mois), je discute avec Delphine Cocteau-Senn et Rodolphe Bigot (dialogue de juristes avec un actuaire) sur la protection des données personnelles en assurance. Tout avait commencé il y 15 mois, avec une table ronde organisée à Amiens, avec un resto la veille, les discussions pendant la table ronde, puis des dizaines et des dizaines d’échanges de messages, de coup de fils, de restos…. Rodolphe Bigot m’avait sollicité pour une autre table ronde sur le droit de l’assurance à Caen, et de mon côté, j’avais demandé à Delphine Cocteau-Senn d’intervenir dans la formation Actuariat – Data Science de l’Institut des Actuaires, depuis l’an dernier. Bref, on a eu de très nombreux échanges, et quand on voit les prises de positions ici et là (pour ne pas parler de brassage d’air), on s’était dit que nos discussions pourraient éclairer les débats qui agitent actuaires et juristes.
C’était d’autant plus intéressant, je pense, qu’on a vraiment appris, mutuellement, de ces échanges. Echanges entre un actuaire et des juristes, tous les trois universitaires. Depuis plusieurs semaines, on tente de mettre par écrit ces échanges, et le résultat est aujourd’hui en ligne. La forme retenue est une présentation juridique de plusieurs points (en lien avec des notions évoquées dans le Réglement Général sur la Protection des Données – n° 2016/679), quelques questions actuarielles, et une ébauche de réponse. On a longuement hésité sur la forme, mais on a gardé la forme de la discussion, probablement plus agréable à lire qu’un texte écrit à trois. Je mets ici l’introduction (le document complet est en ligne)
Vendredi, je serais à la Défense toute la journée pour un séminaire “Analytics and Actuarial Pricing” chez BNP Paribas Cardif.
Slides are available online.
Within a few days, we will launch the fourth actuarial pricing game.
Un petit billet rapide, sans aucune prétention. Juste pour revenir sur des trucs bizarres sur l’encodage… J’avais commencé à m’y intéresser il y a 3 ans, alors que je discutais avec des utilisateurs de R en Inde et au Japon. Mais revenons à la base… Si on regarde la page “unicode” de Wikipedia (ou plus précisément, si on regarde le code source, on peut avoir des informations (dans les métadonnées) : en l’occurrence, on voit que la page est encodée en UTF-8
C’est cohérent car de plus en plus de documents en ligne semblent encodé dans ce format. Si on revient à la source, l’encodage c’est lié au passage par les bytes 0/1. Plein de sites expliquent (en gros) la philosophie. Quand j’étais petit, on nous expliquait les caractères ASCII, qui utilisent 7 bits ce qui permet de disposer 128 caractères uniquement, numérotés de 0 à 127. Genre “A” était le 65ème caractère, et “a” le 97ème. Il y avait des symboles exotiques (comme “@” qui est le 64ème) ou la ponctuation, mais pas d’accents, comme “à”. Si on veut les lettres avec des accents, il faut plus que 128 caractères. Est-alors arrivée ISO/CEI 8859, qui proposait d’encoder les caractères sur 8 bits (et pas 7). On parle aussi de ISO 8859-1 plus populaire sous le nom “latin-1” très populaires dans les pays latins, avec le “à“, et puis on a continué d’enrichir pour arriver au ISO 8859-15 (latin-9) avec dedans le symbole “€”. Pour les langues asiatiques, il a fallu enrichir encore. Bref, c’est un jeu sans fin. Dans les années 90, est arrivée la norme ISO 10646 correspondant au “jeu universel de caractères” ou UCS pour Universal Character Set (en anglais), qui donnera Unicode. Je crois que les puristes font la différence entre les jeux de caractères (l’ensemble de caractères auxquels on attribue à chacun un point de code unique, la liste des caractères ASCII par exemple) et l’encodage (qui est la façon de représenter ce point de code en mémoire). Mais on va faire simple…
Là où ça se complique c’est que les OS utilisaient (par défaut) des encodages différents. Sous Windows, UTF-16 et Unicode sont majoritairement utilisé. Sous Mac OS,il y avait un truc appelé MacRoman mais il semble que UTF-8 l’emporte aujourd’hui. Enfin, sous linux, historiquement c’était du latin-1 par défaut, mais progressivement l’UTF-8 semble prendre le dessus. Bref, tous ceux qui naviguent entre les OS le savent : les accents c’est l’enfer. Et le principal soucis c’est qu’il est impossible, quand on récupère un fichier d’avoir son encodage a priori. Regardons le fichier suivant BaseANSI.txt
On peut tenter de regarder s’il y a des sortes de méta-données dans le document…
mais ce n’est pas le cas. Le document est un document brut, avec juste les données, et aucune métadonnées. Rien sur l’encodage (en particulier, car c’est ce qui nous intéresse aujourd’hui). Pour l’anecdote, si au lieu d’ouvrir le fichier avec mon bloc-note sous Windows, on tente sur un mac, on a
Et quand je passe par le terminal DOS sous Windows, c’est la même chose (ce qui est un peu rassurant)
Regardons un autre fichier BaseUNICODE.txt. Si vous voyez une différence avec l’autre, dites le moi…
Pourtant il y en a une : la preuve, dans un terminal, on peut lire les accents du fichier,
Un dernier ? Regardons BaseUTF8.txt,Là encore, on a l’impression d’avoir le même fichier. Et pourtant non.
Commençons par récupérer les trois fichiers (via R)
download.file(url="http://freakonometrics.free.fr/BaseUTF8.txt", destfile="BaseUTF8.txt") download.file(url="http://freakonometrics.free.fr/BaseANSI.txt", destfile="BaseANSI.txt") download.file(url="http://freakonometrics.free.fr/BaseUNICODE.txt", destfile="BaseUNICODE.txt") |
Techniquement, ces documents sont différents, même si dans le bloc note, on ne voyait rien : il n’ont pas la même taille !
file.size("BaseUNICODE.txt") [1] 322 file.size("BaseANSI.txt") [1] 163 file.size("BaseUTF8.txt") [1] 171 |
On peut essayer de les lire, pour voir…
read.table("BaseANSI.txt",header=TRUE,sep=";") Département No Habitants 1 Ain 01 631877 2 Ardèche 07 324209 3 Bouches-du-Rhône 13 2016622 4 Corse-du-Sud 2A 152730 5 Côte-dOr;21;533147\nCôtes-dArmor 22 598357 |
On a un petit soucis avec l’apostrophe, mais rien de grave
read.table("BaseANSI.txt",header=TRUE,sep=";",quote="\"") Département No Habitants 1 Ain 01 631877 2 Ardèche 07 324209 3 Bouches-du-Rhône 13 2016622 4 Corse-du-Sud 2A 152730 5 Côte-d'Or 21 533147 6 Côtes-d'Armor 22 598357 |
R arrive parfaitement a lire (sur mon ordinateur du moins) cette première base de données. Tentons une seconde
read.table("BaseUTF8.txt",header=TRUE,sep=";",quote="\"") ï..DÃ.partement No Habitants 1 Ain 01 631877 2 Ardèche 07 324209 3 Bouches-du-Rhône 13 2016622 4 Corse-du-Sud 2A 152730 5 Côte-d'Or 21 533147 6 Côtes-d'Armor 22 598357 |
On a manifestement un soucis, un soucis d’encodage. Le soucis c’est que, d’ordinaire, on ne sait pas comment le fichier a été encodé. Et comme on l’a vu, il n’y a pas de métadonnées qui nous donne l’encodage. Comment faire ? Heureusement, il existe une fonction sous R qui devine l’encodage. Littéralement.
library(readr) guess_encoding("BaseUTF8.txt", n_max = 1000) # A tibble: 3 x 2 encoding confidence 1 UTF-8 1.00 2 ISO-8859-1 0.620 3 ISO-8859-2 0.430 |
Autrement dit, le best-guess est ici un encodage UTF-8. Essayons pour voir,
read.table("BaseUTF8.txt",header=TRUE,sep=";",encoding="UTF-8",quote="\"") X.U.FEFF.Département No Habitants 1 Ain 01 631877 2 Ardèche 07 324209 3 Bouches-du-Rhône 13 2016622 4 Corse-du-Sud 2A 152730 5 Côte-d'Or 21 533147 6 Côtes-d'Armor 22 598357 |
Ça marche ! ou disons que ça marche presque. J’ai un soucis dans mes libellés de variables. Mais disons que c’est anecdotique ici.
Tentons le dernier fichier
read.table("BaseUNICODE.txt",header=TRUE,sep=";",quote="\"") ÿþD 1 NA 2 NA 3 NA 4 NA 5 NA 6 NA 7 NA 8 NA |
Comme ça ne marche pas, essayons d’avoir un best-guess sur l’encodage
guess_encoding("BaseUNICODE.txt", n_max = 1000) # A tibble: 3 x 2 encoding confidence 1 UTF-16LE 1.00 2 ISO-8859-1 0.530 3 ISO-8859-2 0.270 |
On en tient un !
read.table("BaseUNICODE.txt",header=TRUE,encoding="UTF-16LE") ÿþD 1 \n 2 B 3 \n 4 C 5 \n 6 C 7 \n 8 C |
Damned, ça ne marche pas… Heureusement, Ewen est venu à ma rescousse. La solution est relativement étrange
read.table("BaseUNICODE.txt",header=TRUE,fileEncoding="UTF-16LE", sep=";", quote="") Département No Habitants 1 Ain 01 631877 2 Ardèche 07 324209 3 Bouches-du-Rhône 13 2016622 4 Corse-du-Sud 2A 152730 5 Côte-d'Or 21 533147 6 Côtes-d'Armor 22 598357 |
Oui, on n’utilise pas l’encodage via “encoding” mais “fileEncoding” ! Subtil non… Cela dit, si on voulait être cohérent on devrait utiliser la fonction du package pour lire la base, non ? Car dans le package “readr” il y a une fonction “read_table”
read_table("BaseUNICODE.txt",locale=locale(encoding = "UTF-16LE")) Error in guess_types_(datasource, tokenizer, locale, n = guess_max) : Incomplete multibyte sequence |
Sauf que quand on la lance, on a un message d’erreur… Bref, l’encodage c’est compliqué. C’est ce qui est expliqué dans le forum de discussion du package cela dit. Je ne suis pas sûr que mon billet serve à quoi que ce soit, mais je voulais garder une trace de tout ça !
Vendredi dernier, j’étais convié à une matinée de discussion autour de la communication scientifique, par des universitaires, à destination des doctorants de la région. Muriel Lefebvre a fait une présentation captivante sur la communication scientifique, et comme promis, je vais mettre par écrit quelques notes prises (je pourrais préciser que j’ai rajouté des réflexions personnelles, mélées à mes notes, car je n’ai pas réussi à déméler). L’idée était de questionner la communication scientifique, et ses pratiques diverses. Le premier point était de noter la variété des formes, la variété des objectifs (écrire permet de réfléchir, de communiquer, de valider, etc), et la variété des publics. En effet, on peut écrire pour soi, pour les pairs (présents, mais aussi passés, en revisitant des études), et pour d’autres sphères. Pour illustrer ce dernier point, Muriel Lefebvre a proposé de revenir sur la rose des vents proposée par Callon, Courtial & Penan (1993), que je reprends sommairement
(bon, le géomètre qui est en moi bout, car une rose des vents est une étoile à 32 directions, pour les 32 directions du vent, car 32=25 or ici on n’a que 5 branches…. mais peu importe). La “communication scientifique” est un objet difficile à identifier parce qu’elle est confrontée à ces 5 axes, qui sont de nature très différente.
En haut à gauche, on retrouve les deux branches classiques du monde académique : la communication liée à l’enseignement, et celle liée à la recherche académique. Dans ces deux sphères, l’enseignant-chercheur est légitime parce que c’est dans son nom : il est enseignant et chercheur. Le terme “média” est à entendre au sens large, il s’agit de toucher un grand public. La “sphère économique” est liée à d’autres termes importants aujourd’hui, comme l’innovation, ou la valorisation. Pour ce dernier, si on pense naturellement aux applications industrielles, on peut toutefois rajouter des activités commerciales au sens large, comme cela peut se voir en droit. Enfin, les “institutions” pourrait être un terme fourre-tour, dès lors que l’enseignant-chercheur est appelé à partager son expertise (par exemple dans des discussions en lien avec la bioéthique).
Il y a aussi une confusion sur les termes, où les termes de publication, de médiation, de vulgarisation, sont utilisés dans tous les sens. Historiquement, la communication scientifique passe par l’écriture, qui légitime le travail. La publication permet de diffuser des travaux, tout en servant d’instrument de validation, d’assurer la paternité d’une découverte, de valoriser un travail, et elle sert aujourd’hui de base pour évaluer la recherche. Cette publication scientifique passe par un processus éditorial et une évaluation par les pairs. La chaîne de l’écriture scientifique est un processus non linéaire, largement décrit par Hans Reichenbach et Karl Popper, mais qui historiquement, restait dans la sphère académique. La vulgarisation (du latin vulgaris) vise à transmettre un contenu (scientifique) d’une manière aisément compréhensible par un large public. Son apogée a été au XIXème siècle, et a été basé sur une relation fondamentalement hiérarchique, du savant vers l’ignorant, du sacré vers le profane. La transmission est alors unidirectionnelle. Ce terme a aujourd’hui une connotation négative, on lui préferera celui d’instruction, ou d’éducation populaire. Cette vulgarisation n’est alors pas une fin en soi, mais une méthode. Puis dans les années 60/70, on voit arriver l’idée de médiation scientifique. Il s’agit alors de mettre du lien, de dialoguer, dans une époque qui remet en question l’autorité de la science, qui n’associe plus science et progrès aussi aisément. Ce dialogue qui doit être rétabli entre science et société se retrouve plus généralement entre la culture et la société. C’est l’époque des maisons de la culture, nées de l’idée d’ouvrir et de rentre accessible l’art (et la culture) à tout le monde. Par la suite, dans les années 80, la notion de culture va s’élargir, en reconnaissant des formes plurielles et diverses. Et la science va bénéficier de ce mouvement. C’est l’idée de “(re)mettre la science en culture“. La culture scientifique est alors une sous partie de la culture commune.
Ces changements vont se traduire par des nouvelles formes de publicisation de la science, pour reprendre la terminologie de Jürgen Habermas. Il s’agit pour le scientifique de rendre public ses travaux, d’être visible. La publicisation se fait par les revues savantes traditionnelles, mais aussi les revues électroniques, les archives ouvertes (hal, arxiv, etc), les réseaux sociaux académiques (research gate, academia, etc), qui restent tous dans la sphère académique traditionnelle. Mais il y a aujourd’hui également de la publicisation vers un publique non-académique, qui peut être atteint par les carnets et blogs de recherche, ou par des chaînes youtube. Depuis 1982, et la redéfinition du rôle de l’enseignant-chercheur, l’importance de la diffusion n’a cessé d’augmenter.
Les expériences de vulgarisation/médiation ne cessent de se multiplier, allant de la fête de la science, à la nuit européenne des chercheurs ou MT180. L’objectif n’est plus très clair, même si les institutions académiques semblent valoriser ces activités, allant d’un souhait de partager des savoirs à une forme d’injonction de visibilité. Pour conclure, Muriel Lefebvre est longuement revenue sur cette double injonction, de qualité et de visibilité. Les fontières semblent se brouiller, Il y a une demande d’humanisation (voire une personnalisation) de la recherche.
On a ensuite eu une table ronde, et ma prise de notes en a pâti (on s’est un peu écouté parler). Il y a eu de nombres interrogations sur la légitimité pour faire de la médiation, et sur l’ouverture des revues savantes. Mais une intervention a aussi souligné l’émergence d’un nouvel acteur depuis quelques années dans la médiation : les associations. Le point de départ semble avoir été les associations en lien avec les personnes malades, atteintes du HIV, dans les années 80, qui ont voulu comprendre. Ces associations peuvent jouer un rôle important dans la médiation. L’histoire amusante est un article de Muriel Lefebvre, intitulé The Circulation of Scientific Articles in the Sphere of Web-Based Media: Citation Practices, Communities of Interests and Local Ties. Dans cet article, elle s’interroge sur la vie d’un article paru en 2012, dans PLoS, intitulé “Freshwater Killer Whales”: Beaching Behavior of an Alien Fish to Hunt Land Birds, et plus particulièrement qui mentionne cet article. Pour information, cet article étudie un gros poisson, la silure (Silurus glanis)
et ce poisson semble avoir tendance à manger des oiseaux, comme ces pigeons (Columbia livia). L’article illustre avec quelques photos (les plus courageux pour aussi trouver en ligne des vidéos).
La question n’était pas de savoir qui dans les revues scientifiques citent cet article (comme peut le faire la bibliométrie traditionnelle) mais plutôt de comprendre au sens large qui a pu mentionner cette étude (publiée en accès ouvert). En fait, cet article a été énormément mentionné sur des sites internet d’associations de pêcheurs (comme seriouslyfish), qui peuvent relayer une information publiée par des chercheurs, pour la diffuser auprès d’un public beaucoup plus large. En fait, le processus de la production scientifique est complexe, car ces associations peuvent aussi faire remonter un problème auprès de scientifiques, qui prendront alors le temps de creuser, et de comprendre. Il s’agit alors véritablement d’échanges. Plus généralement aujourd’hui, on retrouvera beaucoup d’association de malades, qui tentent de comprendre la recherche, mais qui tentent aussi de sensibiliser (et donc de financer). Il n’y a plus de hiérarchie, entre sachants et profanes, mais au contraire une vraie médiation jouée par le chercheur.
Pour le cours de magistère, suite à la première partie du cours, présentant les méthodes de simulation, faite par Magali Fromont, nous allons travailler sur le boostrap, en partant de l’article de Bradly Efron, Six Questions Raised by the Bootstrap. L’article date un peu (il a été publié en 1990) mais il pose six questions qui permettront de revenir sur le fonctionnement du bootstrap (et de la statistique mathématique).
Je me permets d’en rajouter une septième
Ce dernier point pose la question du subsampling. Je renvoie aux notes de cours de Charles Geyer qui présente ce concept. Je recommande à tous la lecture de What Teachers Should Know about the Bootstrap. Enfin, pour un ouvrage de référence plus moderne que l’article de Bradley Efron, je ne peux que recommander les chapitres 10 et 11 de Computer Age Statistical Inference.
In July, I will give a lecture at the 2018 edition of the Summer School at the UB School of Economics, in Barcelona. It will be a four day crash course. Registration should be opened.
Lecture 1: Introduction : Why Big Data brings New Questions
Lecture 2: Simulation Based Techniques & Bootstrap
Lecture 3: Loss Functions : from OLS to Quantile Regression
Lecture 4: Nonlinearities and Discontinuities
Lecture 5: Cross-Validation and Out-of-Sample diagnosis
Lecture 6: Variable and model selection
Lecture 7: New Tools for Classification Problems
Lecture 8: New Tools for Time Series & Forecasting
Ce vendredi, toute la matinée, je participe à une table ronde sur le thème communication scientifique : nouveaux visages, nouvelles frontières, sur le campus de Beaulieu, à l’université de Rennes. J’essayerais de prendre des notes, et de faire un billet, en fin de semaine.
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)) |