Tag Archives: modèle

Justice “actuarielle”, algorithmes… et données

Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré surveillé et punir de Nicolas Bourgoin, qui contient un chapitre sur le sujet.

Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlait de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le coeur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie…

Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice.

Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable Y. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable X_3 et une autre X_4, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse \tilde{X}_2, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons X_1 le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par X_1. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)] Autrement dit, ma probabilité est une fonction croissante (linéaire) de X_3, et si je suis riche (ce qu’on note X_2=R), j’ai moins de chance d’être condamné, avec un facteur k\in[0,1], toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes.

Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe A ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.

 n=1000000
 r=.95
 k=1/3
 S=matrix(c(1,r,r,1),2,2)
 library(mnormt)
 set.seed(1)
 vectX=rmnorm(n,varcov=S)
 X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A"))
 X1=relevel(X1,"A")
 X2=cut(vectX[,2],breaks=c(-10,1,10),c("P","R"))
 X3=runif(n)
 X4=runif(n)
 Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P")))
 B0=data.frame(X1,X2,X3,X4,Y1)

On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement)

B1=data.frame(X1=X1[Y1==1],X2=X2[Y1==1],X3=X3[Y1==1],X4=X4[Y1==1])

Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées

 mean(B0$X1=="A")
[1] 0.158734
 mean(B1$X1=="A")
[1] 0.05016072
 t.test(B0$X1=="A",B1$X1=="A")
data:  B0$X1 == "A" and B1$X1 == "A"
t = 219.5, df = 1283500, p-value < 2.2e-16
95 percent confidence interval:
 0.1076038 0.1095428
sample estimates:
 mean of x  mean of y 
0.15873400 0.05016072

Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de X_3, qui était de 50% dans la population et de 58% dans la population condamnée

 t.test(B0$X3,B1$X3)
data:  B0$X3 and B1$X3
t = -163.46, df = 844940, p-value < 2.2e-16
95 percent confidence interval:
 -0.08455430 -0.08255058
sample estimates:
mean of x mean of y 
0.5000094 0.5835618

alors que pour X_4, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….

 t.test(B0$X4,B1$X4)
data:  B0$X4 and B1$X4
t = -0.81161, df = 48645, p-value = 0.417
95 percent confidence interval:
 -0.003869879  0.001603452
sample estimates:
mean of x mean of y 
0.5001326 0.5012659

Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de X_3, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.

 set.seed(2)
 n1=nrow(B1)
 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))
 B2=cbind(B1,Y2)
 reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

J’ai construit ici mon modèle permettant de calculer une probabilité de récidive. Le modèle est construit sur le niveau de diplôme X_1 (qui n’influence en rien, encore une fois, la probabilité de commettre un crime), la variable X_3 (qui effectivement joue un rôle) et la variable de bruit X_4.

 set.seed(2)
 n1=nrow(B1)
 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))
 B2=cbind(B1,Y2)
 reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

Je fais ici un score logistique car c’est le plus logique.

 summary(reg)
Coefficients:
            Estimate Std. Error z value  Pr(|z|)    
Inter.) -4.04753    0.11592 -34.916    2e-16 ***
X1C      0.61311    0.10914   5.618 1.93e-08 ***
X1B      0.68807    0.09977   6.896 5.34e-12 ***
X3       1.00610    0.08052  12.495    2e-16 ***
X4      -0.13562    0.07286  -1.861   0.0627 .

Premier élément à observer, la variable X_4 est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve X_1, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux B ou C) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné.

Plus intéressant, si je regarde la valeur du coefficient, en fonction de k, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante.

Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base.

On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens

qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative.

Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc.

Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances….

Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score Y d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches !

Faire parler les chiffres… n’importe comment

Cette fin de semaine, Martin Grandjean a mis en ligne un billet intéressant sur son blog, sur l’utilisation des statistiques (dans un but de propagande). L’exercice n’est pas nouveau, mais Martin soulève des questions, malheureusement importantes et complexes. Dans un paragraphe, intitulé “faire parler les chiffres… n’importe comment” (que j’ai repris comme titre, j’avoue avoir hésité avec “with great power comes great responsibility“), on retrouve l’analyse (rapide) d’un graphique, présenté ci-dessous. Le graphique parle de “statistiques” liées à des problèmes démographiques et d’immigration, sauf que comme le note Martin, il y a un “petit problème” : “ces statistiques n’en sont pas vraiment. Ou plutôt, il s’agit d’une extrapolation libre et “linéaire”“.

Je reviendrais sur ces points dans deux minutes, avant je voulais revenir sur une petite phrase que l’on retrouve un peu plus loin: “comment peut-on se permettre de prédire une tendance sur cinquante ans à partir d’une tendance observée sur dix ans ?“. Cette remarque est importante… car beaucoup de monde se pose cette question, et les statisticiens se doivent d’apporter une réponse. Ce problème m’avait été posé (plus ou moins) par une journaliste il y a quelques années maintenant, sur la solvabilité des compagnies d’assurance. Avec Solvabilité II, en Europe, les compagnies doivent calculer des VaR à 99.5%, c’est à dire liés à des événements ayant une période de retour de 200 ans. Avec au mieux 25 ans de données ! Dans ce cas, la nuance vient du fait que l’on cherche à estimer un quantile associé à une probabilité faible avec peu de données (et que l’on assimile rareté et période temporelle). C’est le même problème qui se pose en hydrologie, quand on veut construire une digue suffisamment haute pour que seul un événement millénaire provoque une inondation, mais que l’on doit prévoir avec 50 ans de données. Cela dit, quand un assureur fait de la mortalité prospective, il a souvent 20 ou 30 ans de données, et doit estimer la probabilité qu’un assuré de 25 ans soit encore en vie dans 50 voire 75 ans. C’est le genre de message que j’essayais de faire passer il y a quelques années dans une table ronde sur la pérennité des régimes de retraites, en France. Bref, cet exercice est un vrai exercice. Et il convient de bien comprendre le modèle pour comprendre ce qu’on fait ! Dire que c’est stupide ne pourra pas suffire.

Car quand on a peu de données (j’ai envie de dire que c’est le cas ici), c’est le modèle qui va imposer les conclusions: tous les modèles semblent valides quand on a peu de données, et les conclusions dépendront plus du modèle retenu que des données. C’est moins le cas quand on a des gros volumes de données (sans parler de big-data, je garde ça pour un autre billet).

Revenons sur le problème en question, car il est intéressant. On nous parle ici d'”extrapolation linéaire” dans la publicité. “Extrapolation” signifie que l’on va construire un modèle, et qu’on va l’utiliser pour faire une prévision. Et il y a “linéaire“. A priori, un modèle linéaire, ce n’est pas trop ambigu… sauf qu’on peut faire un modèle sur le niveau, ou sur le taux de croissance. Si on fait un modèle très simple – constant – sur le taux de croissance, on aura alors une croissance exponentielle, c’est aussi simple que ça. Donc le modèle peut être linéaire, ou exponentiel (on regardera un peu les deux).

Autre point. On a ici trois populations,

  • les suisses, 
  • les étrangers, 
  • la population totale, 

Le graphique est un peu trompeur… Il donne l’illusion que l’on extrapole la première variable  et la troisième, . Le danger est qu’il convient, dans ce cas, de contraindre un minimum les modèles, car . On pourrait très bien imaginer qu’en extrapolant, on ait . Donc il faudra faire un peu attention… Le plus simple sera probablement de construire des modèles sur  et  en s’assurant de la positivité de ces deux variables.

Dernier point. A-t-on le droit de faire des modèles indépendants. Peut-on supposer que les deux séries  et  évoluent indépendamment l’une de l’autre ? A priori non, et dans ce cas, c’est un peu plus complexe, car il faut un modèle bivarié. Bref, je trouve que ce petit exemple donne matière à réflexion sur ce qu’est la modélisation statistique. Et sur l’incroyable pouvoir qu’on les statisticiens et les économètres (car ceux sont eux qui construisent le modèle… “with great power comes great responsibility” comme dirait l’autre).

Car il faut que les choses soient claires. Il est impossible d’être neutre ! Ça n’a pas de sens… On peut construire plusieurs modèles, qui sembleront tous valides, et ensuite utiliser un critère de choix de modèles (les plus usuels étant peut-être le critère AIC d’Akaike, ou le BIC de Schwarz). Mais là aussi, je n’ai aucun moyen pour dire : ce modèle est le meilleur, et il est neutre. Non, il sera le meilleur pour un critère de choix, que je choisis.

Histoire d’être plus clair, essayons de faire des régressions. Le modèle le plus simple serait un modèle linéaire,

avec des résidus – a priori – non corrélés, . C’est que donnerait deux régressions linéaires sous Excel, pour ceux qui ne sont pas à l’aise avec la modélisation. Classiquement, ce modèle s’écrit

avec  (conditionnellement à notre variable explicative). Lançons nous… Bon, quand j’ai cherché les données, je n’ai pas trouvé 3 observations, mais une bonne trentaine pour la population résidente, et les étrangers. On va donc utiliser toutes les données pour faire notre projection (et éventuellement voir ce que donnerait une régression avec seulement 3 observations)

> D=read.table("http://freakonometrics.free.fr/suisse.csv",sep=",",
+ header=TRUE)
> tail(D)
      X      N1      N3      N2
27 2006 1554527 7508739 5954212
28 2007 1602093 7593494 5991401
29 2008 1669715 7701856 6032141
30 2009 1714004 7785806 6071802
31 2010 1766277 7870134 6103857
32 2011 1815994 7954662 6138668

Comme on a un modèle linéaire, on peut faire indifféremment l’estimation du modèle sur deux des trois séries. Par simplicité, modélisons des deux séries qui nous intéressent. L’estimation donne ici

> reg2=lm((N2/100000)~X,data=D)
> reg3=lm((N3/100000)~X,data=D)
> pred2=predict(reg2,newdata=data.frame(X=1980:2060))/100000
> pred3=predict(reg3,newdata=data.frame(X=1980:2060))/100000
> lines(1980:2060,pred2,col=COL[1],lwd=2)
> lines(1980:2060,pred3,col=COL[6],lwd=2)

Maintenant, il faut s’interroger sur la pertinence de notre modèle. Pour ma part, j’ai rarement vu des modèles linéaires pour modéliser des des évolutions de populations. Mais sur du court terme (une centaine d’années), pour quoi pas. Même avec un modèle de type Maltusien, c’est possible. Pour ma part, quand on travaille sur des données de comptage, j’aime beaucoup la régression de Poisson,

A la différence de notre premier modèle, on aura ici de l’hétéroscédasticité : plus la population grande, plus la variance (et donc le terme d’erreur) sera grande. C’est plus réaliste que le premier modèle. Rien de plus facile que d’estimer ce modèle

> reg2=glm((N2/100000)~X,data=D,family=poisson(link="identity"))
> reg1=glm((N1/100000)~X,data=D,family=poisson(link="identity"))

On peut en profiter d’ailleurs pour regarder un peu les intervalles de confiance (sur notre modèle) et relativiser un peu les conclusions

> pred2p=predict(reg2,newdata=data.frame(X=1980:2060),type="response",se.fit=TRUE)
> pred1p=predict(reg1,newdata=data.frame(X=1980:2060),type="response",se.fit=TRUE)

La taille de l’intervalle de confiance nous suggère ici d’être prudent avant d’affirmer quoi que ce soit… et encore, on a de la chance, on a ici 35 ans de données et pas 3 !

> I=which(D$X%in%c(1990,2000,2010))
> reg2=glm((N2/100000)~X,data=D[I,],family=poisson(link="identity"))
> reg1=glm((N1/100000)~X,data=D[I,],family=poisson(link="identity"))

Bon, classiquement, quand on fait une régression de Poisson, on utilise un lien exponentiel,

On va alors s’assurer que les deux comptages soient toujours positifs ! Ce qui est bien compte tenu de notre problème. Mais en contrepartie, on a la contrainte d’avoir une croissance (ou une décroissance) exponentielle de nos populations. On n’a pas trop le choix (car quand on y pense, avec un modèle linéaire, il y a forcément une date, avant laquelle, ou après laquelle, la population sera négative… ce qui est gênant)

> reg2=glm((N2/100000)~X,data=D,family=poisson)
> reg1=glm((N1/100000)~X,data=D,family=poisson)
> pred2p=predict(reg2,newdata=data.frame(X=1980:2060),type="response")
> pred1p=predict(reg1,newdata=data.frame(X=1980:2060),type="response")

(la courbe en pointillé en arrière plan est le modèle linéaire).

On continue ? On peut en fait avoir un modèle bivarié, pour dire que nos deux courbes évoluent ensemble. Il existe par exemple le modèle de Poisson à choc commun, de densité jointe

avec . Sous R, le package qui permettait de faire cette régression n’est plus disponible, mais peu importe, on peut le recoder, ça prendra trois minutes,

> f=function(Z,L){
+ si=function(i) choose(Z[1],i)*choose(Z[2],i)*gamma(i+1)*
+ (L[3]/(L[1]*L[2]))^i
+ s=Vectorize(si)(0:min(Z))
+ p=exp(-sum(L))*L[1]^Z[1]/gamma(Z[1]+1)*L[2]^Z[2]/gamma(Z[2]+1)*sum(s)
+ return(p)
+ }
> minuslogL=function(B){
+ h=function(x) exp(x)
+ logL=function(i) log(f(Y[i,],
+ c(h(B[1]+B[2]*X[i]),h(B[3]+B[4]*X[i]),h(B[5]+B[6]*X[i]))))
+ return(-sum(Vectorize(logL)(1:n)))
+ }
> optim(c(lm(log(Y[,1])~X)$coefficients,
+ lm(log(Y[,2])~X)$coefficients,0,0),minuslogL)->maxL
> Bstar=maxL$par
> Bstar
  (Intercept)             X   (Intercept)             X               
-4.343823e+01  2.303903e-02 -3.430377e+00  3.746384e-03 

 1.506016e-02 -9.743153e-04 
> predbiv=function(x,B=Bstar){
+ h=function(x) exp(x)
+ return(c(h(B[1]+B[2]*x)+h(B[5]+B[6]*x),
+          h(B[3]+B[4]*x)+h(B[5]+B[6]*x)))

Si on construit notre prévision, on est ici proche de notre précédant modèle,

On pourrait aussi utiliser un modèle de série temporelle: dire que l’on a une tendance linéaire, ou exponentielle, et dire que le bruit peut être modélisé par un modèle autorégressif (bivarié ou pas).

avec un modèle autorégressive pour le bruit,

ou de manière beaucoup plus générale,


i.e. un modèle autorégressif vectoriel,

On peut estimer ce modèle assez facilement,

> reg2=glm(N2~X,data=D,family=poisson)
> reg1=glm(N1~X,data=D,family=poisson)
> Z=cbind(residuals(reg2),residuals(reg1))
> library(vars)
> regvar=VAR(Z,p=1)

Si on regarde le modèle, on notera que la matrice d’autocorrélation est bien pleine, avec de la causalité dans tous les sens (ce qui devrait augmenter encore la taille de nos intervalles de confiance)

> summary(regvar)

VAR Estimation Results:

Estimation results for equation y1: 
=================================== 
y1 = y1.l1 + y2.l1 + const 

      Estimate Std. Error t value Pr(>|t|)    
y1.l1  1.13635    0.05191  21.891   <2e-16 ***
y2.l1 -0.11922    0.03411  -3.495   0.0016 ** 
const  0.73037    0.52316   1.396   0.1737    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimation results for equation y2: 
=================================== 
y2 = y1.l1 + y2.l1 + const 

      Estimate Std. Error t value Pr(>|t|)    
y1.l1  0.44698    0.12384   3.609  0.00118 ** 
y2.l1  0.83294    0.08139  10.234 5.76e-11 ***
const  0.87575    1.24811   0.702  0.48868    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Je continue ? Allez, un petit dernier exemple pour la route ! Regardons le petit modèle suivant,

Avec ce modèle, le ratio entre les deux populations sera toujours le même. Ce qui fera taire tous les débats… et pour l’estimer, il suffit de rajouter la contrainte dans les codes précédents (avec l’estimation par maximum de vraisemblance)

Maintenant, comme on avait des séries temporelles, je n’ai pas présenté de modèles nonparamétriques, car c’est délicat de faire de la prévision avec ces modèles. Par contre, si on avait des données individuelles, on aura pu s’en donner à cœur joie, à faire du lissage dans tous les sens ! Ce qui nous aurait donné encore plus de latitude dans notre modèle. Je passe aussi sous silence ici que mes modèles sont probablement trop simples. Je veux dire par là – je connais mal les procédures pour avoir la nationalité suisse – que quand ses enfants, et ses petits-enfants, sont nés en Suisse, j’ai du mal à me dire que la personne est un “étranger“. Au bout de quelques années, on peut imaginer un modèle dynamique où les étrangers passent de “étranger” à “Suisse“. Mais j’ai peut-être un biais canadien dans ma perception de la nationalité.

Ce que je voulais dire, c’est que faire des prévisions à partir d’un (relativement) petit jeu de données n’est jamais simple, et que le choix du modèle va complètement déterminer la visualisation que l’on veut avoir pour nos données, et du type de prévision que l’on cherche. Et ne me demandez pas d’être neutre, c’est impossible… Une fois qu’on a compris ce problème, on comprend que certains chercheurs avancent que plus de la moitié des études publiées (souvent sur des petits échantillons) sont fausses. Je ne dirais pas pour ma part qu’elles sont “fausses“. Juste qu’il est facile, quand on a peu de données, de faire passer le message que l’on veut ! Ici, c’est facile à démonter (Martin l’a fait en quelques lignes, et il m’en aura fallu un peu plus car j’ai besoin d’un peu plus de temps – et de formalisme – pour raconter mes histoires), et il s’agissait uniquement d’un dessin de propagande. Il faut garder en mémoire que beaucoup de décisions de santé publique sont basés sur ce genre d’exercice. Avec en plus des données rarement publiques, et entachées de beaucoup de bruit… On pensera aux études sur les OGM, sur les antennes relais, sur le tabagisme (évoqué dans un précédant billet), etc. Et pour démonter ce genre d’étude, cela prend beaucoup plus de temps !

Ordres d’un processus ARMA

Dans la méthodologie de Box & Jenkins, une étape qui arrive très rapidement est le choix des ordres du processus http://freakonometrics.blog.free.fr/public/perso6/ARMA__latex.gif, une fois que l’on validé l’hypothèse de stationnarité de la série, comme on l’a vu en cours la semaine passée. Considérons la série du trafic autoroutier,

> autoroute=read.table(
+"http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)

Un outils pratique de sélections des ordres http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif dans un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMAlatex.gif est la fonction d’autocorrélation étendue. La définition est donnée dans les notes de cours (Def. 223) à partir des statistiques proposées par Tsay & Tiao (1984),

>  EACF=eacf(A7)
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o x x x x x o o x  x  x  o
1 x x o o x x x o o o x  x  x  x
2 o x o o o o x o o o o  x  x  o
3 o x o x o o o o o o o  x  x  x
4 x x x x o o o o o o o  x  o  o
5 x x o o o x o o o o o  x  o  x
6 x x o o o x o o o o o  x  o  o
7 o x o o o x o o o o o  x  o  o
>  EACF
$eacf
           [,1]       [,2]        [,3]        [,4]
[1,]  0.6476234  0.2124105 -0.02413173 -0.24234535
[2,]  0.4889076  0.2797152 -0.02494135 -0.06094037
[3,]  0.1006541 -0.2285771  0.03514148  0.08000588
[4,]  0.1390240 -0.2788742  0.04386746  0.28307260
[5,] -0.5091680  0.3144899 -0.34572269  0.31450865
[6,] -0.4224571 -0.4877505  0.16054232 -0.09130728
[7,] -0.4731353 -0.4324857 -0.04847184 -0.10500350
[8,] -0.2129591 -0.4072901  0.09487899 -0.06493243
             [,5]        [,6]        [,7]        [,8]
[1,] -0.514330187 -0.61634046 -0.52314403 -0.28008661
[2,] -0.254912957 -0.28966664 -0.33963243 -0.21863077
[3,] -0.156624357 -0.01199786 -0.25116738  0.13079231
[4,] -0.183283544  0.03651508 -0.08711829  0.11626377
[5,] -0.190885091 -0.09786463 -0.09182557  0.08818875
[6,] -0.072904071  0.29271777 -0.09334712  0.01972648
[7,] -0.009873289  0.36909726  0.01698660 -0.03317456
[8,]  0.020485930  0.38342158  0.16981715 -0.02592442
             [,9]       [,10]       [,11]     [,12]
[1,] -0.082607058  0.15178887  0.56583179 0.8368975
[2,] -0.085026400  0.02731460  0.26158357 0.7844748
[3,]  0.001866994 -0.11658312  0.03038621 0.7026361
[4,]  0.025183793 -0.21608692  0.05660781 0.6674301
[5,] -0.006831894 -0.02514440 -0.07390257 0.4762563
[6,]  0.010058718 -0.03888613 -0.04382043 0.5338091
[7,]  0.032124905 -0.07022090 -0.04427400 0.4674165
[8,]  0.024189179 -0.20818201  0.01459933 0.4830369
          [,13]       [,14]
[1,]  0.5637439  0.17862571
[2,]  0.4530716  0.24413569
[3,]  0.2534178 -0.20160890
[4,]  0.2409861 -0.29462510
[5,] -0.1517324 -0.14763294
[6,] -0.1701182 -0.28771495
[7,] -0.1651515  0.05466457
[8,] -0.1403198 -0.04095030

$ar.max
[1] 8

$ma.ma
[1] 14

Les ronds dans la matrice désignent des valeurs non-significatives. Par défaut, le nombre de retards, pris en compte pour la composante autorégressive est faible, mais on peut l’augmenter.

> EACF=eacf(A7,13,13)
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12 13
0  x o o x x x x x o o x  x  x  o
1  x x o o x x x o o o x  x  x  x
2  o x o o o o x o o o o  x  x  o
3  o x o x o o o o o o o  x  x  x
4  x x x x o o o o o o o  x  o  o
5  x x o o o x o o o o o  x  o  x
6  x x o o o x o o o o o  x  o  o
7  o x o o o x o o o o o  x  o  o
8  o x o o o x o o o o o  x  o  o
9  x x o o x o x o o o o  x  o  o
10 x x o o x o x o o o o  x  o  o
11 x x o o x x x o o o x  x  o  o
12 x x o o o o o o o o o  o  o  o
13 x x o o o o o o o o o  o  o  o
> EACF$eacf
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,]  0.64  0.21 -0.02 -0.24 -0.51 -0.61 -0.52 -0.28 -0.08  0.15
[2,]  0.48  0.27 -0.02 -0.06 -0.25 -0.28 -0.33 -0.21 -0.08  0.02
[3,]  0.10 -0.22  0.03  0.08 -0.15 -0.01 -0.25  0.13  0.00 -0.11
[4,]  0.13 -0.27  0.04  0.28 -0.18  0.03 -0.08  0.11  0.02 -0.21
[5,] -0.50  0.31 -0.34  0.31 -0.19 -0.09 -0.09  0.08  0.00 -0.02
[6,] -0.42 -0.48  0.16 -0.09 -0.07  0.29 -0.09  0.01  0.01 -0.03
[7,] -0.47 -0.43 -0.04 -0.10  0.00  0.36  0.01 -0.03  0.03 -0.07
[8,] -0.21 -0.40  0.09 -0.06  0.02  0.38  0.16 -0.02  0.02 -0.20
[9,] -0.14 -0.50 -0.10 -0.06  0.11  0.38 -0.01 -0.02 -0.08 -0.20
[10,] -0.29  0.48  0.00  0.18  0.27 -0.12  0.41  0.00  0.16  0.15
[11,]  0.24  0.48 -0.04  0.02  0.46  0.10  0.37 -0.10  0.05  0.16
[12,] -0.59  0.49 -0.18  0.05  0.31 -0.32  0.34 -0.16  0.07  0.16
[13,] -0.31  0.31 -0.12  0.16 -0.11  0.04 -0.04  0.12 -0.09  0.00
[14,]  0.47  0.26  0.11  0.13 -0.01 -0.01  0.00  0.07 -0.03  0.02

On peut aussi visualiser graphiquement les différentes valeurs, avec en ordonnée les ordres autoréregressifs (http://freakonometrics.blog.free.fr/public/perso6/platex.gif) et en abscisse les ordres moyenne mobile (http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif).

> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> ceacf=matrix(as.numeric(cut(EACF$eacf,
+1,nrow(EACF$eacf),
+ ncol(EACF$eacf))
> for(i in 1:ncol(EACF$eacf)){
+ for(j in 1:nrow(EACF$eacf)){
+ polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
+ col=CL[ceacf[j,i]])
+ }}


Un bruit blanc est en bas à gauche (http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif nuls). Dans cette méthode, dite méthode des coins, on cherche un coin telle que dans le quadrant supérieur, les valeurs soient non-significatives (ternes sur le dessin ci-dessus). La figure ci-desssous correspond à l’analyse d’un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMA126latex.gif

Les fortes valeurs positives sont en bleu foncé, les fortes valeurs négatives sont en rouge foncé. On peut aussi regarder la fonction suivante, qui utilise une identification de modèle par MINIC (Minimum Information Criterion) à l’aide du critère de Shcwarz (BIC, ou SBC Schwarz’s Bayesian Criterion)

> armaselect(A7,nbmod=5)
p q      sbc
[1,] 12 1 1441.798
[2,] 13 0 1442.628
[3,] 12 0 1443.188
[4,] 12 2 1443.362
[5,] 14 0 1445.069

Enfin, une dernière fonction possible est évoquée dans la section 6.5. du livre deCryer & Chan (2008),

> ARMA.SELECTION=
+ armasubsets(A7,nar=14,nma=14,ar.method='ols')
> plot(ARMA.SELECTION)

basé sur le critère de Schwarz.
Je vais me répéter, mais ces méthodes ne sont que des outils, histoire d’avoir des pistes si on ne sait trop dans quelle direction partir. Et compte tenu de la saisonnalité de la série, je serais pour ma part parti sur un modèle avec http://freakonometrics.blog.free.fr/public/perso6/platex.gif=12, histoire de voir si on ne pourrait pas avoir un modèle simple, et facilement interprétable en plus…

  1. -3):3)/3,labels=1:6 []