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 !