Category Archives: Statistics

Computing AIC on a Validation Sample

This afternoon, we’ve seen in the training on data science that it was possible to use AIC criteria for model selection.

> library(splines)
> AIC(glm(dist ~ speed, data=train_cars, 
  family=poisson(link="log")))
[1] 438.6314
> AIC(glm(dist ~ speed, data=train_cars, 
  family=poisson(link="identity")))
[1] 436.3997
> AIC(glm(dist ~ bs(speed), data=train_cars, 
  family=poisson(link="log")))
[1] 425.6434
> AIC(glm(dist ~ bs(speed), data=train_cars, 
  family=poisson(link="identity")))
[1] 428.7195

And I’ve been asked why we don’t use a training sample to fit a model, and then use a validation sample to compare predictive properties of those models, penalizing by the complexity of the model.    But it turns out that it is difficult to compute the AIC of those models on a different dataset. I mean, it is possible to write down the likelihood (since we have a Poisson model) but I want a code that could work for any model, any distribution….

Hopefully, Heather suggested a very clever idea, using her package

And actually, it works well.

Continue reading Computing AIC on a Validation Sample

Lillie versus Kolmogorov-Smirnov

Toujours pour répondre à une question posée par mail (“ça quoi ça sert le test lillie sous R ? “), un court billet. Avant de parler de Lillie1, il faut revenir sur Kolmogorov-Smirnov. L’idée est de tester

https://perso.univ-rennes1.fr/arthur.charpentier/tex/KS-01b.png

contre

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-02.png

Je prends une forme paramétrique pour la fonction que l’on cherche à tester car c’est ce qui est souvent fait en pratique. Et surtout, ça permettra de mieux comprendre (je pense) l’idée du test lilllie. L’idée du test repose sur le théorème de Glivenko-Cantelli. En fait, Andreï Kolmogorov et Vladimir Smirnov ont montré que l’on pouvait encadrer la vitesse de convergence,

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-03.png

ce qui se montre en notant que https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-04.png converge en loi (si l’hypothèse H0 était la bonne) vers un pont brownien changé de temps par la fonction quantile https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-05.png.
La loi limite est appelée loi de Kolmogorv-Smirnov, et ne dépend pas de la loi sous-jacente. Par contre elle n’est qu’asymptotique. Pour récupérer la loi à distance finie, regardons les simulations suivantes. On notera en particulier que la loi de la différence ne dépend pas de la vraie loi des observations (même si c’est cette loi que l’on teste, comme toujours lorsque l’on cherche la loi sous l’hypothèse H0)
> n=20
> ns=100000
> SU=SN=SE=rep(NA,ns)
> for(i in 1:ns){
+ X=runif(n)
+ SU[i]=max(c(abs(rank(X)/n-X),abs((rank(X)-1)/n-X)))
+ X=rnorm(n,0,1)
+ SN[i]=max(c(abs(rank(X)/n-pnorm(X,0,1)),abs((rank(X)-1)/n-pnorm(X,0,1))))
+ X=rexp(n,1)
+ SE[i]=max(c(abs(rank(X)/n-pexp(X,1)),abs((rank(X)-1)/n-pexp(X,1))))
+ }
> plot(density(SN),col=”red”,lwd=2)
> lines(density(SU),col=”blue”,lwd=2)
> lines(density(SE),col=”green”,lwd=2)
avec la densité de la loi du maximum entre la fonction de répartition théorique et la fonction de répartition empirique enrouge pour l’échantillon Gaussien, enbleu pour l’échantillon uniforme et envert pour l’échantillon exponentiel. Peu importe la loi sous-jacente, la statistique de test est la même.

On peut en particulier récupérer les quantiles de cette lois pour n=20 observations (la loi du test dépend simplement du nombre d’observations).

> quantile(SE,.95)
95%
0.2937023
> quantile(SE,.90)
90%
0.264661
On peut d’ailleurs confirmer ces valeurs en regardant dans les tables (ici ou ), par exemple 0.294 à 95%,

Mais ce test ne sert que pour juger de l’adéquation à une loi précise (que nous avions noté https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-encore.png ici). En pratique, on ne connaît pas cette loi, mais on peut l’estimer, en particulier pour les lois paramétriques. On peut alors chercher à tester autre chose, comme

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-06.png

contre

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-07.png

qui s’apparente plus à un test d’ajustement à une famille de lois. Ce test étant différent du précédant, il n’y a pas de raison d’avoir la même loi pour la statistique de test (qui néanmoins, à  peut rester celle que nous avions considérée, à savoir la norme infinie entre les fonctions de répartition). Graphiquement par exemple, on cherchait aurparavant à tester l’ajustement de la loirouge, et maintenant on va tester l’ajustement de la loibleue,

qui va s’adapter à l’échantillon, contrairement au test initial,

Le code ici, dans le cas Gaussien par exemple (la loi sous-jacente n’est plus neutre dans ce cas) donne
> n=20
> ns=100000
> SN=SNML=rep(NA,ns)
> for(i in 1:ns){
+ X=rnorm(n,0,1)
+ SN[i]  =max(c(abs(rank(X)/n-pnorm(X,0,1)),
+               abs((rank(X)-1)/n-pnorm(X,0,1))))
+ SNML[i]=max(c(abs(rank(X)/n-pnorm(X,mean(X),sd(X))),
+               abs((rank(X)-1)/n-pnorm(X,mean(X),sd(X)))))
+ }
> plot(density(SN),col=”blue”,lty=2)
> lines(density(SNML),col=”red”,lwd=2)
Les densités sont sensiblement différentes,

et les quantiles aussi,
> quantile(SN,.95)
&nbsp
;     95%
0.2943206
> quantile(SNML,.95)
95%
0.1923295
Sous R, deux fonctions distinctes sont à utiliser. Par exemple l’ajustement à une loi https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-08.png donne
>  X=rnorm(20)
> mean(X)
[1] -0.553641
> sd(X)
[1] 0.8293167
>  ks.test(X,”pnorm”, mean=0,sd=1)

One-sample Kolmogorov-Smirnov test

data:  X
D = 0.3096, p-value = 0.03362
alternative hypothesis: two-sided
(on rejette ici l’ajustement à une loi normale centrée réduite) et l’ajustement à une loi https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-09.png donne
> lillie.test(X)

Lilliefors (Kolmogorov-Smirnov) normality test

data:  X
D = 0.1519, p-value = 0.2629
avec une p-value de l’ordre de 26%. Notons que cette p-value est interprétable, contrairement au test suivant, qui n’est pas valide,
>  ks.test(X,”pnorm”, mean=mean(X),sd=sd(X))

One-sample Kolmogorov-Smirnov test

data:  X
D = 0.1519, p-value = 0.69
alternative hypothesis: two-sided
Encore une fois, la statistique de test est la même, sauf que la loi de la statistique de test sous H0 n’est pas juste dans le dernier cas. En particulier, n sera beaucoup plus exigent dans le test car on a choisi la meilleure loi normale possible (à l’aide d’un critère de maximisation de vraisemblance). Sinon je renvoie ici pour un document sur le même sujet.
1 Lillie pour Hubert Lilliefors qui, le premier, a proposé ce test

… et d’où sort le modèle probit ?

Histoire de poursuivre le billet sur la régression logistique (ou logit, ici), je continue avec la régression probit. Chester Bliss (papier de 1934 en ligne ici) travaillait, en tant que biologiste, sur la conception de pesticides. Au lieu de travailler (comme habituellement) sur le pourcentage de nuisibles détruits, mais sur une grandeur qu’il a appelé “probability unit“, soit “prob…it“, qu’il condensa en probit. De manière arbitraire, un pourcentage de tués valant 0,0001 correspondait à une probit nulle, et un pourcentage de 0,9999 à une probit de 10.
Il avait d’ailleurs proposé une table de correspondance, permettant de passer du pourcentage au probit. Il proposait alors de relier le probit au logarithme de la dose de pesticides produits, ce qui devait aboutir – dû moins l’espérait-il – à une ligne droite. Il appela ce modèle un modèle probit.
David Finney ( pour un article rétrospectif) a poursuivi cette idée, toujours dans le domaine toxicologique. Dans son livre de 1952, on retrouve une table de convertion associée à cette grandeur introduite par Bliss

Mais pour aller plus loin, il faut remonter à un papier publié par Finney et Stevens en 1948, où ils eurent la bonne idée de renormer les grandeurs de Bliss en translatant tout de 5, autrement dit un pourcentage quasi-nul (0,0001) correspondait à une probit de -5, et un pourcentage de 0,9999 à une probit de 5. Aussi la probit de 0 correspondait au cas où le pourcentage valait 1/2. Ils notent ainsi que la probit Y associée à une proportion P est la grandeur qui vérifie

Voilà pour la petite histoire…

Mais qui connaît la loi logistique ?

Certains élèves ont pris peur en voyant apparaitre la loi logisitque dans le cours d’économétrie des variables qualitatives, car autant j’avais passé pas mal de temps lors des remises à niveau sur la loi normale, la loi de Student, la loi de Poisson…. jamais je n’avais présenté la loi logistiquue.. Mais qui connaît la loi logisitique ?

  • La loi logistique (en tant que loi)

En maths, on a introduit très tôt les fonctions de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zevzev.png

où K et r sont des réels positifs et a un réel quelconque, qui sont ce que l’on appelle les solutions (en temps continu) du modèle de Verhulst.

C’est Pierre François Verhulst qui a donné le nom de “courbe logistique“, en expliquant, en 1845  « nous donnerons le terme de logistique à cette courbe ». Comme toujours, il faut revenir aux racines étymologique, en notant que  « logistique » a même racine que logarithme, et que logistikos signifie “calcul” en grec (on peut retrouver plus d’information ici , ou encore ).

Ces fonctions ont été introduites par Pierre François Verhulst, qui était un élève d’Adolphe Quetelet, et qui dévait étudier un modèle d’évolution de population qui ne soit pas exponentielle. Cette courbe a été utilisée dans l’étude des populations à partir des années 20 par Raymond Pearl et Lowell Jacob Reed, mais qui ne créditent Verhulst de la paternité de la découverte qu’en 1922. En revanche, on ne retrouve le terme de “logistique” qu’en 1924 dans une correspondance entre George Yule et Reed. On notera que c’est Joseph Berkson qui défendra l’idée d’ajuster certaines courbes par une fonction logistique (modèle logit) plutôt que par la fonction de répartition de la loi de Gauss (modèle probit, introduit par Bliss en 1934, mais c’est une autre histoire que je raconterais sûrement un jour….).

En dynamique des populations, le modèle de Verhulst est un modèle de croissance proposé vers 1840. Pour reprendre l’explication de wikipedia sur le sujet (inspiré de plusieurs ouvrages sur le sujet dont de Handbook on logisitic distribution) Verhulst a proposé ce modèle en réponse au modèle de Malthus qui proposait un taux d’accroissement constant sans frein conduisant à une croissance exponentielle de la population. Le modèle de Verhulst imagine que le taux de natalité et le taux de mortalité sont des fonctions affines respectivement décroissante et croissante de la taille de la population. Autrement dit, plus la taille de la population augmente, plus son taux de natalité diminue et son taux de mortalité augmente.

On appelle y la taille de la population, m(y) le taux de mortalité et n(y) le taux de natalité. La taille de la population suit l’équation différentielle

https://perso.univ-rennes1.fr/arthur.charpentier/latex/585468azdf.png

Si m et n sont des fonction affines respectivement croissante et décroissante (ce qui une hypothèse un peu forte, mais qui a longtemps été utilisée comme on l’avait vu en cours d’actuariat en M1) alors n – m est une fonction affine décroissante. L’équation peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/325365ze4v.png

avec a et b deux réels positifs, si la croissance est positive quand y est proche de 0. En posant K=a/b, l’équation devient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/36854z6e5v46z.png

avec a > 0 et K > 0 . Aussi, histoire de formaliser, on recherche des fonctions strictement positives définies sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/r+.png et vérifiant le système https://perso.univ-rennes1.fr/arthur.charpentier/latex/5za365c4a.png et

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/54za6cf8.png

ce qui conduit à la solution dite logistique
https://perso.univ-rennes1.fr/arthur.charpentier/latex/az4c65.png
  • La cote au lieu de la probabilité

Bref, voilà pour la loi…. mais comment voit-on apparaître cette loi en économétrie ? Pour rappel, on cherche à modéliser la probabilité de survenant d’un évènement. Mais il existe d’autre outils que la probabilité en tant que telle. En particulier, la cote d’un évènement – ou odds in favor of an event” – est la quantité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/pipipiphn.png

où p est la probabilité de réalisation de cet évènement.
Cette grandeur est bien connue par les turfistes…. Bon, dans ce cadre (les paris), cette probabilité n’a rien d’une probabilité historique, mais il s’agit plutôt d’une probabilité risque neutre comme on dit en finance: il s’agit de la probabilité extraite des prix. Considérons un exemple simple (où on ne mise que sur le vainqueur)

  • Le parieur A mise 10 € sur le cheval numéro 1.
  • Le parieur B mise 5 € sur le cheval numéro 2.
  • Le parieur C mise 20 € sur le cheval numéro 3.
  • Le parieur D mise 15 € sur le cheval numéro 4.

On suppose que les sommes misées reflètent les croyances des joueurs (qui ont tous la même aversion pour le risques). Autrement dit, les joueurs pensent que le cheval 1 a deux fois plus de chances de gagner que le cheval 2. Bref, La cote d’un cheval est obtenue en divisant la somme des enjeux joués sur les autres chevaux par la somme des enjeux joués sur le cheval choisi (ou l’inverse, ça dépend de

la manière dont on souhaite le présenter), soit ici

  • Pour le cheval 1, une cote de [5+20+15]/10=4 (on dira 4 contre 1)
  • Pour
    le cheval 2, une cote de [10+20+15]/5=9 (on dira 9 contre 1)
  • Pour le cheval 3, une cote de [10+5+15]/20=3/2 (on dira 1.5 contre 1)
  • Pour le cheval 4, une cote de [10+5+20]/15=7/3

Ca existe également en football (comme en témoigne le ticket ci-contre où figurent les cotes). Bref, cette cote est une manière très naturelle de parler de probabilité. L’énorme avantage est qu’une probabilité est définie sur [0,1] alors qu’une cote peut prendre n’importe quelle valeur positive.

  • La régression logistique

On souhaite donc modéliser la cote à l’aide de variables explicatives. Mais comme le score peut être négative, il est assez “naturel” de prendre le logarithme (comme pour la régression Poissonnienne par exemple). Aussi, on peut comprendre l’écriture de la forme suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ojker%c3%b4bvjk%c3%b4.png

ou encore

https://perso.univ-rennes1.fr/arthur.charpentier/latex/hopvezuiogvp.png

Or la fonction qui apparaît à gauche est la fonction quantile de la loi logistique. Quod erat demonstrandi…

  • Un cas particulier de GLM

La transformation logistique apparaît très naturelle quand on utilise les GLM (modèles linéaires généralisés). En effet, il faut que la loi s’écrive sous la forme exponentielle
(comme j’en avais parlé ici ou , par exemple). Pour reprendre les notations du cours d’actuariat 2, il faut écrire la loi binomiale

https://perso.univ-rennes1.fr/arthur.charpentier/latex/khpok%c3%b9ikh%c3%b9.png

sous la forme^

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ezv65465.png

soit ici (on n’a pas de paramètre de surdispersion)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/6ze54c65.png

En notant que la loi binomiale se réécrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/azcac6ea46.png

ou encore

https://perso.univ-rennes1.fr/arthur.charpentier/latex/6c48ze684fg.png

On en déduit que le paramètre naturel est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zegzgh8796h8.png

et que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/az6f54ea68g54.png

(moyennant un peu de réécriture). Bref, on note une relation naturelle (on parlera de fonction de lien canonique) de la forme entre le paramètre usuel de la loi binomiale et le paramètre dit naturel de la loi exponentielle,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zegzgh8796h8.png

ou encore

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zevzg7h9.png

La fonction de répartition logistique apparaît là encore presque par enchantemant…. enfin, presque…