Category Archives: Statistics

Statistics

Une courte histoire des expériences randomisées

5 Comments

On a beaucoup parlé ces derniers mois des expériences randomisées. Cet automne pour commencer, lorsque le prix Nobel d’économie a été attribué à Esther Duflo, Michael Kremer et Abhijit Banerjee. Plus récemment, lors de la polémique, en France, où, dans une conférence donnée le 13 février 2020, intitulée contre la méthode, un médecin marseillais affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ».

Les expériences randomisées

Quand on se demande l’effet d’un traitement (par rapport à un autre), la bonne méthode consiste à se poser la question « que ce serait-il passé si la personne avait choisi l’autre option ? » – ce que John Stuart Mill parlait de la « méthode des différences ». On parle ici d’inférence causale car on ne peut pas savoir ce qu’aurait donné l’autre option, qui n’a pas été retenue. L’effet causal d’une option, ou d’un traitement, n’est jamais observé car il s’agit de la différence potentielle entre deux alternatives. La méthode la plus simple consiste à prendre deux patients proches (à défaut d’être identiques), et de voir comment les deux répondent aux deux options, la proximité se mesurant en fonction de covariables (comme l’âge, le genre, etc), et les patients n’ayant pas choisi les options proposées. Il y a trois points importants dans cette approche : (1) la notion de groupe de contrôle, autrement dit il faut qu’il y ait des personnes « non traitées » (2) la notion d’aveuglement, ou d’absence d’information, autrement dit, les patients ne doivent pas savoir quel traitement ils ont eu (3) la notion d’affectation aléatoire, autrement dit, le traitement doit être choisi de manière totalement indépendante des patients. Ce sont des points qui opposent fondamentalement les expériences randomisées aux données observationnelles, comme les oppose par exemple Paul Rosenbaum dans Observation and Experiment. Si la méthode d’expérience randomisée a été formalisée par Ronald Fisher en 1934, nous allons revenir ici sur deux moments importants dans l’histoire de cette méthodologie. La première est l’utilisation de la randomisation pour faire des tests autour de 1890, en particulier avec les travaux de Peirce puis Richet, sur le spiritualisme et la transmission de pensée (qui ont très probablement influencé Fisher). Et la seconde est l’expérience qui a été menée en 1954 aux États-Unis, pour tester l’efficacité d’un vaccin contre la polio (et qui montre au passage que faire des expériences randomisées dans un contexte de maladies infectieuses a tout son sens).

Continue reading Une courte histoire des expériences randomisées →

Statistics

Computing AIC on a Validation Sample

10 Comments

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 →

Statistics

On a toujours besoin de petits poids (chez soi)

Leave a comment

Après avoir mis en ligne le billet d’hier sur les manifestants, j’ai passé la journée à me demander “pourquoi ai-je mis des poids tenant compte d’un effet taille dans ma régression” ?. Sur le moment ça me semblait très naturel, et plus la journée avançait, plus je me mettais à douter…

Bon, pour rappel quand on fait une régression classique, sans poids, et sans constante car on veut un coefficient de proportionnalité, par “moindres carrés“, on résout

ce qui donne

La dernière écriture montre que l’on fait une moyenne pondérée des pentes pour chacune des observations… Graphiquement, ça donne

les traits pointillés permettent de montrer ou se trouve le point moyen, avec la moyenne sur toutes les villes pour les syndicats et pour la police, pour les deux dates. Les flèches indiquent les pentes moyennes, c’est à dire le nombre total de manifestants en France selon les syndicats divisé par le nombre total de manifestants en France selon la police. Bref, la droite de régression “classique” n’a pas pour pente le ratio moyen…

Une alternative usuelle est d’introduire des poids, et on cherche à résoudre

dont la solution a la forme

Hier j’avais proposé de prendre des poids correspondant à la taille de la ville (pour intégrer un effet médiatique – les médias ne parlant parfois que des nombres de manifestants dans les très grosses villes), mais on pourrait proposer, sans le même genre d’idée que de prendre (et qui est représenté ci-dessous)

Toutefois, l’idée la plus pertinente aurait dû être http://freakonometrics.blog.free.fr/public/perso/wls08.png ce qui donne comme estimateur par moindres carrés pondérés

http://freakonometrics.blog.free.fr/public/perso/wls09.png

c’est à dire le ratio moyen national… Bref, on est colinéaire avec nos vecteurs tracés depuis le début

Autrement dit, il aurait été plus malin ici d’attribuer des petits poids aux grosses villes pour avoir un coefficient interprétable….

Mais au débat du débat (probablement légitime) sur le choix de la régression, ce qui est troublant c’est que plus les poids donnés aux grosses villes sont importants, plus la pente est importante… autrement dit, on modèle quadratique pourrait être plus adapté, ou au moins quelque chose de non linéaire, et de strictement croissant. Sur le graphique suivant, on voit l’estimation par moindres carrées locaux, et par splines,

ce qui pourrait nous conforter dans l’idée que “les grandes villes ont des plus grands trottoirs“(ramené à la taille de la chaussée !), ce qui peut surprendre un peu.. Si on regarde l’estimation des pentes obtenues par ces deux méthodes, on a, pour les régressions locales

(avec en trait horizontal en pointillés les pentes des flèches, i.e. les taux observés nationalement), et par splines

On retrouve l’idée que dans les petites villes, le comptage est plus précis (des deux côtés) avec un facteur de 2 et 3 pour les deux manifestations respectivement. En mars dernier, on avait ensuite un facteur croissant passant de 2 à 4 en fonction de la taille de la manifestation. Et pour la dernière manifestation en revanche, c’est un peu n’importe quoi (si je peux me permettre).

Bref, pour revenir au point qui m’intéressait au départ, on note qu’il ne semble pas y avoir de méthode scientifique permettant de passer du nombre estimé par la police à celui proposé par les syndicats, car sinon il y aurait une relative stabilité entre les manifestations… stabilité que l’on n’a pas (quelle que soit le type de régression considéré…).

Statistics

Lillie versus Kolmogorov-Smirnov

Leave a comment

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 là), 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

Statistics

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

Leave a comment

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 (là 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…

Statistics

Mais qui connaît la loi logistique ?

Leave a comment

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 , là ou encore là).

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 là, 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…