Données Personnelles et Assurance

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)

Continue reading Données Personnelles et Assurance

Actuarial Pricing Game

Within a few days, we will launch the fourth actuarial pricing game.

  1. People have to register for that new game. People have to fill a form online before February 28th to get a training dataset. As in our previous versions, players can be individuals or in a team (identified as one player).
  2. When registration closes, players will receive a training dataset (training.csv). The training dataset will be based on two years of data on household policies. There will be information about the insurance policy from underwriters, as well as, claims data. They will also receive a pricing dataset (pricing.csv) with only underwriters information and no claims. The goal is to propose a premium for all records in the pricing dataset. Note that those registered will also receive additional information and template files after registration closes.
  3. Before  April 9th, 2018, players will have to provide two files: a pricing function (either a formula.py Python function, or a formula.R R function) and a dataset with prices prices.csv). The pricing function should be an understandable explicit function which, once applied to the training dataset, will yield premium prices for each record. There will be limits on the complexity of this function. To ensure this, players will receive instructions on what can be used to construct formulas either in Python or R.
  4. Before May 14th, 2018, players are asked again to submit two files: a dataset with prices (second_prices.csv) and a pricing model (either a model.py Python file, or a model.R R file). For this round there are no restrictions on what can be used in the model, as long as, by using the model file submitted the prices dataset (second_prices.csv) can be reproduced. Players will again receive further instruction on this after registration.

Unicode, UTF-8… la base de l’encodage

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 !

De la communication scientifique…

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.

Seven Questions Raised by the Bootstrap

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

  • why do maximum likelihood estimated distributions tend to be short-tailed ?
  • why does the delta method tend to underestimate standard errors ?
  • why are cross-validation estimators so variable ?
  • what is a correct confidence interval ?
  • what is a good nonparametric pivotal quantity ?
  • what are computationally efficients ways to bootstrap ?

Je me permets d’en rajouter une septième

  • why generate samples of size n, why not try smaller sizes ?

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.

 

Course on “Big Data for Economics”

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

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