Tag Archives: Fisher

La guerre des étoiles : distinguer le signal du bruit

La grande difficulté dans la modélisation et la construction de modèles prédictifs est de réussir à distinguer le signal et le bruit (pour reprendre le titre du classique de Nate Silver). La réponse statistique est la notion de significativité, et la recherche des ‘étoiles’ dans les sorties de régression. Avec l’explosion du nombre de données, il est devenu crucial de faire cette distinction, de savoir quelles sont les interactions qui sont significatives.

Approches historiques de cette notion de « significativité »

Le débat sur la significativité est ancien, même si sa formulation s’est faite historiquement dans des termes assez vagues. Par exemple, dès 1710, le médecin et mathématicien John Arbuthnot s’était interrogé sur le ratio du nombre de naissance de garçons et de filles, se demandant si la différence était « statistiquement significative ». Pour être plus précis, en utilisant des statistiques sur près de 90 ans, il avait noté “There seems no more probable Cause to be assigned in Physics for this Equality of the Births, than that in our ’first Parents Seed there were at first formed an equal Number of both Sexes.John Arbuthnot pose le premier la question en terme probabilistes. Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un « test de significativité », (là encore avec notre terminologie actuelle). Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après-midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste, en se demandant s’il est « extrêmement probable » qu’il y ait une différence entre les deux mesures. Il avait alors introduit, le premier, un test de comparaison entre des valeurs moyennes : ayant noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important, il en conclue que les séries sont significativement différentes. En 1885, Francis Edgeworth avait repris ces idées, en comparant la taille des criminels, et la taille des gens ordinaires. Mais il faudra attendre les travaux de William Gosset (plus connu sous le nom Student), de Karl Pearson et surtout de Ronald Fisher pour avoir une définition plus rigoureuse de la significativité.

La notion de significativité est cruciale dans la construction de modèles prédictifs. En assurance automobile, l’âge du conducteur est une variable « significative » quand il s’agit d’expliquer la fréquence de sinistres. Formellement, cela signifie que l’âge et la fréquence de sinistre sont corrélées, et que cette corrélation, notée R, est « significativement non-nulle ». En 1921, Ronald Fisher, en proposant la construction d’un intervalle de confiance, propose ainsi un test de significativité. Comme il l’écrit, « from these values, we obtain the difference 0.0471 ± [0.0142] which might well be regarded as significant ». Sur la Figure 1, la figure de droite correspond au cas où la valeur est significativement non nulle (car l’intervalle de confiance à 95% ne contient pas 0) alors que la figure de gauche correspond au cas où la valeur est non-significative (car l’intervalle de confiance à 95% contient 0). Ronald Fisher passera alors plusieurs années à formaliser et expliquer cette idée de significativité statistique.

Figure 1 : distribution théorique de l’estimateur de la corrélation, R, et intervalle de confiance.

Tests, prise de décision et erreurs

Une des contributions majeures de la statistique des années 1920 a été de formaliser la prise de décision. La Figure 2 montre ainsi le mécanisme binaire de la prise de décision, et les deux types d’erreur : rejeter à tort une hypothèse, ou accepter à tort une hypothèse.

Figure 2 : schéma de la prise de décision.

Déclarer qu’une variable, ou une différence, est significative peut conduire à deux types d’erreur : la déclarer comme non-significative, alors qu’elle l’était (prendre un signal pour du bruit) ; et la déclarer comme significative, alors qu’elle ne l’était (prendre du bruit pour un signal). Quand on construit un test, on va essayer de contrôler la probabilité de commettre de telles erreurs.

Mais de faibles taux d’erreur ne veulent pas forcément dire qu’un test est efficace, et le principal danger des tests médicaux est lié à une mauvaise interprétation de ces taux d’erreur. Supposons qu’une maladie touche une personne sur mille, dans une population (selon le test, i.e. 1 personne sur 1000 est ‘positive’). Supposons aussi qu’il soit relative « efficace » au sens où 90% des cas ‘positifs’ sont effectivement atteints, et le test est négatif dans 99% des cas quand est n’est pas touché par la maladie. Ces chiffres sont bien au-dessus de la plupart des principaux tests couramment utilisés. Si le test est effectué sur 10000 personnes, le test sera positif sur 10 personnes (en moyenne), parmi lesquelles 9 sont effectivement malade, mais 1 est saine. A côté, le test sera négatif pour 9990 personnes. Et parmi ces 9990 négatifs, 9890 sont effectivement saines, et une centaine est pourtant malade. Donc au final, sur les 109 personnes malades, 9 ont été détectées, mais 100 sont passés inaperçues, soit un peu plus de 90% des malades ! Ce test supposé efficace ne l’est peut-être pas tant que ça.

Figure 3 : calculs des taux d’erreurs

Le concept de « p-value » et le mythe des 5%

Le concept de p-value est lié justement à l’erreur de rejeter, à tort une hypothèse. Si Ronald Fisher en parle abondement, il faut toutefois attribuer la paternité du concept à son collègue, Karl Pearson. Ce dernier, pour définir la notion de « significativité », utilise la formulation suivante : «  P=.1227 or the odds are now only 8 to 1 against a system of deviations as improbable as or more improbable than this one ». Pour reprendre l’exemple de la corrélation, Pearson nous dit qu’il y a une chance sur 8 pour obtenir une valeur aussi improbable. Comme l’illustre la Figure 4, la p-value est la probabilité d’avoir une statistique aussi grande ou aussi petite que celle obtenue sur l’échantillon, si effectivement la corrélation était nulle. Ou encore p=P{|R| > r | H0],

Figure 4 : distribution théorique de l’estimateur de la corrélation sous l’hypothèse où cette dernière serait nulle et probabilité que |R| dépasse la valeur empirique, observée.

Sur la figure de gauche, la p-value correspondant à l’aire rouge est de l’ordre de 10% alors que sur la figure de droite l’aire est de 0.1%. En fait, si la p-value est inférieur à 5%, alors 0 n’est pas dans l’intervalle de confiance à 95% de l’estimateur de corrélation. Donner la p-value est alors suffisant pour juger de la significativité d’une statistique.

C’est dans le chapitre 4 de son livre que Ronald Fisher pose les bases de la pratique (toujours en vigueur aujourd’hui) des tests statistiques : « if p is between 10% and 90% there is certainly no reason to suspect the hypothesis tested. If it is below 2% it is strongly indicated that the hypothesis fails to account for the whole of facts ». Autrement dit, avec une p-value supérieure à 10%, on peut accepter notre hypothèse (souvent que la variable n’est pas significative) et avec une p-value inférieure à 2%, on va la rejeter (la variable sera alors significative, si on fait un test de nullité d’une corrélation). Mais entre les deux ? Ronald Fisher clôt sans vraiment s’en rendre compte le débat, en affirmant « we whall not often be astray if we draw a conventional line at 5% ». Le choix – aujourd’hui dogmatique – d’un seuil à 5% repose sur l’idée de rejeter à tort une hypothèse avec une chance sur 20, ou bien correspond au fait de s’éloigner de 2 écart-types de la moyenne d’une loi normale centrée-réduite (ce qui arrive avec une chance sur 22).

Cette règle des 5% – une chance sur 20 – est encore en vigueur aujourd’hui. C’est elle que l’on utilise dans tous les modèles économétrique, en regardant les « étoiles » associées à chaque des variables explicatives (3 étoiles si la p-value est inférieure à 0.1%, 1 étoile si elle est inférieure à 5%, et rien au-delà de 10%). Si la -value est une mesure continue de la distance entre l’hypothèse que l’on cherche à tester et les données, ces étoiles ont instauré des seuils, malheureusement supposé infranchissables. Avec cette méthode, comme le notait Guelman (2015), «it seems impressive to see multiple independent findings that are statistically significant but with enough effort it is possible to find statistical significance anywhere. »

Pour illustrer la difficulté de la prise de décision, Power (2014) prend un exemple tiré de la survie à des séances de chimiothérapie, pour guérir un cancer.

Figure 5: chimiothérapie, données Williams et al. (1987).

Dans un “research hospital” le taux de survie est de 86%, alors qu’il est de 81% dans un “non-research hospital”. Un test du chi-deux, mesurant l’indépendance entre le choix de l’hôpital et la survie, donne une statistique de test de 0.5635, correspondant à une p-value de 45.29%. Autrement dit, on devrait être indifférent à l’endroit où on va se faire traiter, puisque l’hôpital est non significativement corrélé à la survie. Mais lorsque la vie est en jeu, est-on prêt à dire qu’une différence de 5 points (86% vs. 81%) n’est pas “statistiquement significative” ? Même avec une p-value proche de 50% !

Le problème des tests multiples

Si cette méthode pour juger de la significativité d’une variable dans un modèle prédictif a eu beaucoup de succès, il convient d’insister sur ses limites quand on se retrouver face à des données massives. Si on cherche à expliquer une variable (comme la sinistralité) par 100 variables, possiblement indépendantes de notre variable d’intérêt, l’analyse précédant devrait nous pousser à retenir 5 variables comme « significativement » corrélées avec la sinistralité.

Figure 6 : Exemple de séries significativement (très) corrélées, données Vigen (2015)

Le problème des tests multiples est d’autant plus important en imagerie médicale. Sur la Figure 7, 100 échantillons Gaussiens sont simulés, de moyenne nulle, et on teste la nullité des moyennes (seules les p-values sont présentées). Imaginons des tests sur des pixels d’une image d’IRM par exemple. Sur une image obtenue sur un individu sain, si des tests sont effectués, pixel par pixel, 10% des points présenteront une p-value inférieure à 10%, 5% une p-value inférieure à 5%, etc.

Figure 7 : 100 tests de nullité de moyenne, sur 100 échantillons simulés (de moyenne nulle).

En l’occurrence, la vision Bayésienne de la prise de décision (décrite dans Greenland & Poole (2013)) peut s’avérer plus juste, en plus d’offrir une interprétation plus claire.

Les alternatives pour juger de la significativité

Les p-values sont un outil important, central, dans la construction de modèles prédictifs. Et la pratique des tests est souvent un exercice intéressant : si la p-value est petite, on se sent conforté dans l’idée que la variable est significative. Mais si elle est grande (disons excède 10%), on ose parfois se dire qu’on n’a pas eu de chance… Dans certains cas, on peut même s’offrir la chance de faire le test sur un autre échantillon et – avec un peu de chance – la p-value sera plus faible.

Comme le note Briggs (2013), les p-value “encourage magical thinking […] they focus attention on the unobservable“. Même avec un regard critique, l’utilisation des p-value est dangereuse. Et l’utiliser de manière automatique, dans un algorithme d’apprentissage, l’est encore plus. Sans bon sens, on verra des variables extrinsèques exotiques utilisées dans les modèles prédictifs (pour reprendre la terminologie de Cass & Shell (1983). A quand un assureur qui utiliserait la pointure de pieds, les résultats au brevet des collèges ou la couleur de la boite à lettres dans son tarif d’assurance auto ?

Références

Arbuthnot, J. (1710). An argument for Divine Providence, taken from the constant regularity observed in the births of both sexes. Royal Society’s Philosophical Transactions.

Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle.

Briggs, W. (2013) Everything Wrong With P-Values Under One Roof. http://wmbriggs.com/…

Cass, D. & Shell, K. (1983). Do Sunspots Matter?. Journal of Political Economy 91 (21): 193–228.

Fisher, R. (1925). Statistical Methods for Research Workers. Oliver & Boyd.

Greenland S. & Poole C. (2013). Living with P values: resurrecting a Bayesian perspective on frequentist statistics. Epidemiology, 24: 62-8.

Guelman, A. (2015). P-values and statistical practice. andrewgelman.com/…

Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Australian & New Zealand Journal of Statistics, 28.

Powers, P.R. (2014). Acts of God and Man Ruminations on Risk and Insurance. ColumbiaUniversity Press.

Silver, N. (2012) The Signal and the Noise: Why So Many Predictions Fail — but Some Don’t, Penguin Press

Vigen, T (2015). Spurious Correlations http://tylervigen.com/…

Calcul(s) d’information de Fisher

La semaine passée, on avait fait quelques calculs pour obtenir l’information de Fisher pour des lois classiques. Je voulais juste remettre au propre les calculs pour les lois à plusieurs paramètres. Pour la loi Gamma,

la log-vraisemblance s’écrit

de telle sorte que

Ici, même pas besoin de prendre une espérance car la Hessienne est constante

Continue reading Calcul(s) d’information de Fisher

Box plot, Fisher’s style

In a recent issue of Significance, I discovered an interesting – and amuzing – figure, about some box & beard plot, in Dr Fisher’s casebook: Beard the statistician in his den.

In French, the box plot (introduced by John Tukey, not George Box, as discussed in a previous post) is popular under the name boîte à moustaches (box with a mustache, for a simple translation).

> set.seed(2)
> x=rnorm(500)
> boxplot(x,horizontal=TRUE,axes=FALSE)
> axis(1)

I was wondering if it was possible to reproduce that Fisher’s style of box plot, with a beard (Ronald Fisher was famous for his beard). Technically, it is not complicated,

> boxplot(x,horizontal=TRUE,xlim=c(-1,1.3),axes=FALSE)
> axis(1)
> Q=quantile(x,c(.25,.75))
> y=cut(x[(x>=Q[1])&(x<=Q[2])],seq(Q[1],Q[2],length=11))
> tb=table(y)
> u=seq(Q[1],Q[2],length=11)
> umid=(u[1:10]+u[2:11])/2
> for(i in 1:10) segments(umid[i],1-.2,umid[i],1-.2-tb[i]/20,lwd=3)

Yes, I kept the mustaches here, but ploting a simple box, one could easily draw a goatee, or a chin strip

> rect(Q[1],1-.2,Q[2],1+.2)
> segments(median(x),1-.2,median(x),1+.2,lwd=2)

The only problem, is that between the first and the third quartiles, the distribution is much more ‘uniform‘ than it might look on the figure in Significance. Unless, perhaps, if add more bars on the left of the first quartile, and on the right of the third one,

> du=diff(umid)[1]
> y=cut(x,seq(Q[1]-du*8,Q[2]+du*8,length=11+16))
> tb=table(y)
> u=seq(Q[1]-du*8,Q[2]+du*8,length=11+16)
> umid=(u[1:26]+u[2:27])/2
> for(i in 1:8) segments(umid[i],1,umid[i],1-.2-tb[i]/20,lwd=3)
> for(i in 19:26) segments(umid[i],1,umid[i],1-.2-tb[i]/20,lwd=3)

Of course, we can try to add a smile on that face, below the mustache,

> vu=seq(-2.5,2.5,by=.02)
> vv=dnorm(vu)
> lines(vu,1-vv*4,col="red",lty=2)
> d=density(x,bw=.1)
> lines(d$x,1-d$y*4,col="red")

So far, I am not convinced. I mean, it should be possible to add something on the box plot,  to add more information. But I still wonder what… maybe using colors?

De la significativité (statistique), suite

Suite à mon rapide billet sur l’article donnant les conclusions d’une étude des effets de la e-cigarette (et la comparaison avec le patch pour arrêter de fumer), j’avais entendu cette phrase “ça veut dire quoi ‘statistiquement non significatif’ ?“. Comme c’est une très bonne question, qui plus est relativement complexe, j’ai pris un peu de temps pour y répondre… Il ne s’agit pas de complexité mathématique, en s’entend car tous les cours de statistique abordent cette notion, mais plus d’une notion complexe à vulgariser. On va revenir un peu sur l’histoire des tests de significativité et sur les travaux de Ronald Fisher (pour comprendre l’origine des 5% utilisé comme seuil lors de l’utilisation de la -value), mais il est possible d’aller voir directement les exemples.

  • Un peu d’histoire

Si on veut commencer par la préhistoire, les travaux les plus anciens qui posent la question de la significativité sont ceux de John Arbuthnot, en 1710. En prenant 82 années de données de natalité, à Londres, il avait observé que chacune des années, il y avait eu plus de naissances de garçons que de filles (évoqué dans Hacking (1975)).

Il s’était demandé, avec la terminologie d’il y a plus de 300 ans, si cette différence pouvait être jugée comme ‘statistiquement significative. Je pense que l’on peut utiliser le terme ‘statistiquement‘ car John Arbuthnot pose la question en termes probabilistes, très clairement,

Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un ‘test de significativité‘, là encore avec notre terminologie actuelle. Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste (que l’on verra formalisé bien plus tard), en se demandant s’il est ‘extrêmement probable‘ qu’il y ait une différence entre les deux mesures. Il avait alors fait un test de comparaison des valeurs moyennes, et noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important. A la fin du XIXème siècle, Francis Edgeworth reprendra cette idée (sur la comparaison de moyennes) lorsqu’il écrira en 1885,

où on retrouve les premières idées qui serviront de base aux travaux ultérieurs. Car pour bien comprendre l’origine de la notion de significativité en statistique, il faut relire les travaux de Ronald Fisher des années 20, mais aussi de Student (William Gosset) et Karl Pearson. Ronald Fisher s’interroge,  dans son papier de 1921, sur la valeur d’un coefficient de corrélation (c’est dans ce papier qu’il introduit formelle sa transformation avec une tangente hyperbolique pour faire un test de nullité d’un coefficient de corrélation), où il note

Aussi, dans le papier, une valeur (ici une corrélation) peut être jugée comme ‘significative‘ (c’est à dire ‘significativement non nulle‘). Encore une fois, l’idée n’est pas nouvelle, Edgeworth aurait dit ‘non accidentelle‘ mais Fisher pose ici la terminologie, en parlant de  ‘significativité‘. L’année suivante, en 1922, il parle de l’idée de tester la significativité, dans le cadre d’un modele de régression (on est toujours dans le contexte de la corrélation)  ou en 1924, on retrouve l’idée de tester la significativité, Mais c’est surtout dans son livre Statistical Methods for Research Workers que l’on va trouver une formalisation claire (et très illustrée) de cette idée, et surtout la connecter à la notion de -value (on ne peut pas parler de significatité sans parler de -value). Pour l’anecdote, dans l’ouvrage de Ronald Fisher, on y parle de -value (oui, avec une majuscule). Cela dit, la paternité de la -value est à attribuer davantage à Pearson, par exemple dans son papier datant de 1900, ou encore, sur un autre exemple avec en plus l’interprétation usuelle en terme de significativité (sans toutefois utiliser ce terme), puisque Pearson dit qu’on a 1 chance sur 8 d’avoir une valeur aussi improbable. Plus loin, il dit aussi qu’arriver avec 1 chance sur 70 peut être qualifier de très improbable (un long débat suivra sur le seuil à partir duquel, effectivement, un événement peut être qualifié d’improbable) Mais revenons au livre de Ronald Fisher, et plus spécifiquement au chapitre 4, qui propose précisément d’étudier la ‘significativité de la différence entre deux moyennes

On y retrouve le test de Pearson basé sur la statistique du chi-deux, comme test d’ajustement (goodness of fit) et d’indépendance. et la comparaison de la -value avec (ce qui va correspondre aux mythiques) 5%. Il parle plus loin d’une ‘ligne conventionnelle‘ que l’on tracera à 5%, Il faut d’ailleurs un lien entre une -value à 5% (correspondant à un événement qui survient avec 1 chance sur 20) et le fait de s’éloigner de 2 écart-types de la moyenne pour une loi normale (qui survient avec 1 chances sur 22) Tout au long du livre, on retrouve évoqué le chiffre mythique 5% qui restera à la postérité. Le plus intéressant pour apprendre comme juger la significativité d’une valeur, on peut regarder les exemples, comme l’exemple 11 l’exemple 12, l’exemple 27, ou encore l’exemple 32 (pour en prendre juste quelques uns) (on retrouve l’utilisation de 2 écart-types pour une loi normale). La discussion amorcée dans ce livre sera poursuivie dans son papier de 1926, qui commence avec Le papier est passionnant, car il pose les bases de ce que font tous les statisticiens depuis 80 ans (il n’est pas surprenant de voir cet article reproduit dans le second volume des Breakthroughs in Statistics), où le seuil de 5% est encore utilisé. Il ne cesse tout au long de l’article de justifier cette règle de ‘un sur vingt‘, tout en autorisant d’autres valeurs, comme 2%, ou 10% (Fisher semble moins orthodoxe que certains utilisateurs de la notion de significativité) Il faudra ensuite attendre les travaux de Jerzy Neyman et Egon Pearson. qui vont axiomatiser et formaliser les tests statistiques, mais ça dépasse (de loin) l’objectif de mon billet d’aujourd’hui. Donc on va repousser à plus tard la présentation de leurs apports. Pour illustrer tout ce qu’on vient de voir, considérons deux exemples simples.

  • Petit exemple, tester l’indépendance

Comme Ronald Fisher, commençons par ce qui me semble le plus simple car on n’a pas besoin de parler modèle probabiliste ici. On se demande si deux variables sont indépendantes, ou pas. Deux cas peuvent se présenter:

  • on a des variables continues, et on veut voir si elles sont non-corrélées (oui, pour commencer on va simplifier, et dire que si elles sont non-corrélées… ça va se traduire par une notion que l’on peut assimiler – rapidement – à une notion d’indépendance. Promis, d’ici 6 mois il y aura des dizaines de billets sur le sujet sur le blog puisque je redonne mon cours sur la dépendance cet hiver). Je laisse les plus motivés jouer avec les données qui traînent sur le blog pour voir si la taille et le poids sont corrélées, ou utiliser un exemple que j’avais l’habitude de présenter dans mes cours de modèles de régression, où on se demande si – chez les femmes – le tour de poitrine et le tour de taille sont corrélés.
  • on a des variables catégorielles, par exemple la couleur des yeux, et la couleur des cheveux. Et on se demande si les deux sont couleurs sont indépendantes, ou pas.

Regardons ce dernier exemple, avec un petit jeu de 25 observations.

> X=c(rep("Brun",8+5),rep("Blond",9+3))
> Y=c(rep("Noisette",8),rep("Bleu",5),rep("Noisette",3),rep("Bleu",9))

On peut résumer l’information dans un tableau de contingence,

> T=as.matrix(table(X,Y))
> T
       Y
X       Bleu Noisette
  Blond    9        3
  Brun     5        8

On va essayer de comparer ce tableau avec ce qu’on devait avoir si les deux variables étaient indépendantes (car c’est ce que l’on souhaite tester : on veut savoir si le tableau que l’on a est significativement différent de ce qu’on aurait avec des variables indépendantes). Il faut ici ressortir la définition de l’indépendance entre deux variables aléatoires, pour tout  et . Par exemple   peut être avoir des cheveux blonds, et  avoir des yeux bleus. Pour rappel, on avait en tout 25 personnes ayant ces deux caractéristiques en même temps, donc empiriquement, on avait observé une probabilité de 9/25 soit 36% (c’est le terme de gauche). Pour le terme de droite, on notera que , c’est la probabilité d’avoir les cheveux blonds. Empiriquement, on en avait 9+3 sur 25, soit 48%. Et  c’est la probabilité d’avoir les yeux bleus. Or on avait 9+5 personnes avec ces caractéristiques, soit 56%. Si les deux variables étaient indépendantes, avec la formule précédente, la probabilité d’avoir à la fois les cheveux blonds et les yeux bleus serait 48% x 56%, soit 26.88%. A comparer aux 36% observés. En fait, au lieu de comparer des pourcentages, on peut aussi comparer des nombres: comme il y a 25 personnes en tout, 26.88% signifie que dans notre groupe, 6.72 personnes devraient avoir les cheveux blonds et les yeux bleus. On peut s’amuser à faire ces calculs pour tous les   et .

> Tind=chisq.test(T)$expected
> Tind
       Y
X       Bleu Noisette
  Blond 6.72     5.28
  Brun  7.28     5.72

Visuellement, on peut mettre ça en forme,

(sur le graphique de droite, les traits noirs correspondent aux nombres sous hypothèse d’indépendance, et en couleur, les valeurs observées).  Bon, maintenant, il faut réfléchir un peu… On se demande si nos deux tableaux sont proches ou pas. Ou, pour utiliser notre mot de la journée, significativement proches. Alors là, ça va devenir technique le temps d’un rapide paragraphe. Dans chaque case du tableau, on compte des gens. Sur nos  personnes, on peut s’attendre à en observer un nombre aléatoire, noté ,  qui suit une loi binomiale (c’est la définition de la loi binomiale). On notera que quand  est grand, la loi normale est alors proche de la loi Gaussienne. Dans chaque case, on a des lois normales (non indépendantes, car les nombres sont contraints, par exemple au total, il faut  personnes), et en prenant les carrés des lois normales (centrées et réduites), on va obtenir une loi du chi-deux (je renvoie vers un précédant billet pour les aspects techniques). Prendre le carré est intéressant, et naturel, car ça fait penser à une distance Euclidienne (usuelle). On va alors poser

où  est le nombre de personnes observées, de caractéristiques (jointes)   et . et  est ce qu’on aurait si les variables étaient indépendantes (calculé auparavant), soit

avec les notations naturelles. Ici,

> sum((T-Tind)^2/Tind)

[1] 3.380994

Bon… on en fait quoi maintenant de ce 3.38 ? C’est là qu’on utilise le petit résultat évoqué auparavant, qui nous dit que, si   est grand, et que les variables sont effectivement indépendantes, alors  va suivre une loi du chi-deux (en l’occurrence à 1 degré de liberté, pour tenir compte des diverses contraintes). On va alors se demander s’il est possible, ou vraisemblable, d’avoir 3.38 avec une loi du chi-deux.

> pchisq(3.380994,df=1)
[1] 0.9340477
> 1-pchisq(3.380994,df=1)
[1] 0.06595227

Aussi, une loi du chi-deux à un degré de liberté a 93.4% chances d’être inférieure à la valeur obtenue, ou encore

 avec 

C’est cette grandeur que l’on appelle -value et que l’on va essayer de comparer au seuil désormais mythique de 5%,

  • si cette probabilité est inférieure à 5%, on va rejeter l’hypothèse d’indépendance
  • si cette probabilité est supérieure à 5%, on va retenir l’hypothèse d’indépendance (pour faire simple)
> chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 3.381, df = 1, p-value = 0.06595

Les plus attentifs devraient me faire remarquer que c’est vaseux mon histoire : tout ce que je raconte est valide si  est grand, et dire que 25 c’est grand… c’est limite ! Tout d’abord, je dirais que c’est ce qui est utilisé dans un paquet d’études en médecine (peut-être un peu plus, mais guère plus). Ensuite je dirais que, malheureusement, la plupart des résultats que l’on voit dans les cours de statistiques sont de ce genre, à savoir des résultats asymptotiques, valides seulement si   est grand (c’est pareil en économie, ou des tests asymptotiques sont parfois évoqués avec 25 années d’observations). Mais heureusement, en s’accrochant un peu, on peut s’affranchir de cette hypothèse. A condition de faire un peu de bootstrap. En fait, l’idée est très très simple. On peut utiliser des générateurs de nombres aléatoires pour générer des tableaux de contingences pour lesquels les variables sont indépendantes. Car générer des variables aléatoires indépendantes, c’est ce qui est le plus simple. En fait, ici, on va se contenter de mélanger les variables,

> s=1
> set.seed(s)
> Xs=sample(X)
> Ys=sample(Y)

En faisant une permutation, je m’assure d’avoir autant de personnes dont la couleur des cheveux est   que dans la population initiale (je ne fais que permuter) et de personnes dont la couleur des yeux est . Par contre, comme je mélange indépendamment, je suis certain que les variables sont indépendantes,

> Ts=as.matrix(table(Xs,Ys))
> Ts
       Ys
Xs      Bleu Noisette
  Blond    7        5
  Brun     7        6

Générer un tableau de contingence, c’est juste pour illustrer, car l’idée est de générer mille (ou plus encore) tableaux de contingence obtenus en simulant des variables indépendantes, et de voir ce que cette statistique du chi-deux donne

> Q=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ Xs=sample(X)
+ Ys=sample(Y)
+ Ts=as.matrix(table(Xs,Ys))
+ Q[s]=sum((Ts-Tind)^2/Tind)}

Si on regarde les valeurs obtenues

>  table(Q) 
Q
0.0509906759906761  0.337162837162837   1.06560106560107   1.92411754911755 
               298                262                200                133 
  3.38099400599401   4.81185481185481    6.9971694971695   9.00037462537462 
                60                 33                  9                  3 
  11.9141275391275    14.489676989677 
                 1                  1

autrement dit, dans 29.8% des scénarios, on avait une distance (un statistique) valant 0.05, dans 26.2% des scénarios, on avait 0.33, et dans 20% des scénarios, on avait 1.06. Etc. Si on cumule on obtient

>  cumsum(table(Q))/10
0.0509906759906761  0.337162837162837   1.06560106560107   1.92411754911755 
              29.8               56.0               76.0               89.3 
  3.38099400599401   4.81185481185481    6.9971694971695   9.00037462537462 
              95.3               98.6               99.5               99.8 
  11.9141275391275    14.489676989677 
              99.9              100.0

Pour rappel, on avait eu

>  sum((T-Tind)^2/Tind) 
[1] 3.380994

qui peut effectivement être obtenu avec des variables indépendantes. Mais seulement dans 6% des scénarios. En fait, dans 95.3% des scénarios, on a eu une valeur inférieure ou égale à celle obtenue sur notre échantillon. Et dans 4.7% des scénarios, on a strictement dépassé cette valeur. On pourrait parler de  -value obtenue par bootstrap, ou rééchantillonage. On voit que l’on flirte ici avec le seuil des 5%… Visusellement, on a

> plot(c(0,as.numeric(names(table(Q)))),c(0,
+ cumsum(table(Q))/10),type="s",xlab="",ylab="")
> abline(v=sum((T-Tind)^2/Tind),lty=2)
> lines(seq(0,15,by=.01),pchisq(seq(0,15,by=.01),1)*100,col="red")

Ici, on est bien embêté pour conclure quoi que ce soit… Ce qui est amusant, c’est qu’avec un échantillon plus grand (trois fois plus grand par exemple), et les mêmes pourcentages, il n’y a plus aucune ambiguïté,

> T=T*3
> T
       Y
X       Bleu Noisette
  Blond   27        9
  Brun    15       24
>  Tind=chisq.test(T)$expected
> Tind Y X Bleu Noisette Blond 20.16 15.84 Brun 21.84 17.16

que l’on peut encore visualiser ci-dessous,

Ici, on obtient une -value de l’ordre de 0.1%,

>  chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 10.143, df = 1, p-value = 0.001449

On peut aussi faire du bootstrap si on trouve que 75 observations, c’est trop faible, mais la conclusion est la même,

  • Petit exemple, sur la comparaison de moyenne de deux groupes

Cette fois, on a deux populations, et on se demande si les deux groupes sont ‘comparables‘. On peut penser à comparer le nombre de personnes guéries au sein de deux populations, une qui a eu un placebo, et l’autre un médicament. Cela dit, cet exemple simple est en fait compliqué car ‘guérie’ n’est peut-être pas une notion claire. Prenons plutôt des quantités que l’on peut mesurer. Disons la taille des gens, et demandons nous si dans un groupe, les gens sont significativement plus grands que dans l’autre. Je vais tricher un peu… je vais prendre deux populations que je sais être différentes, a priori. Un groupes d’hommes, et un groupes de femmes. On va noter  et  ces deux échantillons, car il va falloir formaliser un minimum.

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)] 
> Davis=Davis[order(Davis$height),]
> attach(Davis)
> set.seed(1)
> X=sample(height[sex=="F"],size=6)
> Y=sample(height[sex=="M"],size=5)

On a ici 5 personnes dans le premier groupe et 6 personnes dans le second. Et on se demande si les tailles sont significativement différentes entre les deux groupes,

> (mx=mean(X))
[1] 165.5
> (sx=sd(X))
[1] 4.679744
> (my=mean(Y))
[1] 178
> (sy=sd(Y))
[1] 7.968689

On veut savoir si l’écart entre 165,5 cm et 178 cm est significativement différent de zéro (on pourrait aussi se demander s’il est significativement positif, mais ça serait une autre histoire). On va faire une hypothèse ici : on va supposer que nos échantillons sont obtenus comme des tirages de variables aléatoires, indépendantes,  et . Le modèle qu’on va supposer est que la taille des gens suit une loi normale, conditionnellement à leur sexe,

Si on veut visualiser nos données, on peu utiliser

> u=seq(150,200,by=.25)
> v1=dnorm(u,mx,sx)
> v2=dnorm(u,my,sy)
> plot(u,v1,col="red",type="l",ylim=c(0,.09),axes=FALSE)
> axis(1)
> lines(u,v2,col="blue")
> points(X,rep(.075,length(X)),col="red")
> points(Y,rep(.075,length(Y)),col="blue")
> abline(v=mx,lty=2,col="red")
> abline(v=my,lty=2,col="blue")

 

Le plus simple aurait été de supposer que les deux lois ont la même variance, mais quand on regarde les densités sur le graphique ci-dessus, il faut rester crédible. Bref, on va utiliser le test de Welch. L’idée est d’utiliser un petit exercice de probabilité. Soient  et  des variables aléatoires, Gaussiennes (comme supposé juste avant de faire un dessin) et indépendantes. Posons alors

et

(avec les quantités similaires pour le second échantillon) alors si , et si on note

la variable aléatoire

va suivre une loi de Student (connue depuis les travaux de William Gosset), avec

degrés de liberté. On notera qu’il ne s’agit pas ici d’un résultat asymptotique. Par contre, il repose intégralement sur une hypothèse – forte – de normalité des tailles des individus. L’idée est alors la même qu’auparavant : on va calculer cette statistique sur notre échantillon,

> nx=length(X)
> ny=length(Y)
> s=sqrt(sx^2/nx+sy^2/ny)
> T=(mx-my)/s
> T
[1] -3.091371

et se demander s’il est vraisemblable d’avoir une telle valeur avec une loi de Student. Le nombre de degrés de liberté est ici

> d=(sx^2/nx+sy^2/ny)^2/((sx^2/nx)^2/(nx-1)+(sy^2/ny)^2/(ny-1))
> d
[1] 6.218682

On pourrait alors calculer

  où 

mais ça serait un test pour savoir si la différence est significativement positive, ce qui n’est pas ce qu’on avait dit que l’on cherchait : on veut ici significativement non-nulle.

  où 

Si on veut visualiser un peu tout ça, on peut utiliser

> u=seq(-4,4,by=.01)
> v=dt(u,df=d)
> plot(u,v,type="l")
> abline(v=T,col="red")
> u=seq(-4,T,length=100)
> polygon(c(u,rev(u)),c(dt(u,df=d),rep(0,100)),col="red",border=NA)

Pour calculer ces probabilités, on utilise

> pt(T,d)+(1-pt(-T,d))
[1] 0.02038577

et comme la loi de Student est symétrique par rapport à l’origine,

> pt(T,d)*2           
[1] 0.02038577

Aussi, la -value est de l’ordre de 2%, ce qui est suffisamment petit pour dire que, oui, nos moyennes sont différentes. C’est d’ailleurs ce que renvoie la commande

> t.test(X,Y)

	Welch Two Sample t-test

data:  X and Y
t = -3.0914, df = 6.219, p-value = 0.02039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.31039  -2.68961
sample estimates:
mean of x mean of y 
    165.5     178.0

On avait vu tout à l’heure que l’idée de la valeur de 5% était liée au scénario qui se répète avec 1 chance sur 20 : si on pouvait générer d’autres échantillons, indépendants, alors 1 fois sur 20 on devrait atteindre une statistique de test de l’ordre de celle observée. Pour générer des échantillons vraiment indépendants, on avait utilisé du bootstrap, en rééchantillonnant dans notre population de test. Ici, on a bien plus que nos 11 observations, on en a en réalité 200. On peut se demande ce qui se serait passé si on avait constitué d’autres échantillons,

> Diff=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ X=sample(height[sex=="F"],size=6)
+ Y=sample(height[sex=="M"],size=5) 
+ Diff[s]=(mean(X)-mean(Y))/(sqrt(sd(X)^2/nx+sd(Y)^2/ny))}
> hist(Diff,col="light blue",border="white")

Sur aucun échantillon, la moyenne du premier groupe n’a dépassé la moyenne du second. Alors qu’on aurait pu penser en avoir une vingtaine (2% des mille scénarios générés). Ceci vient peut-être du fait que l’hypothèse de normalité pour les tailles peut être trop restrictive…

  • Pour aller plus loin

Le billet aura été un peu long, mais on aura eu le temps de voir en détails d’où vient cette notion de significativité statistique, et ses liens avec la -value. On aura même eu le temps de voir sur deux exemples l’importance des hypothèses sous-jacentes, en particulier si le test est asymptotique, ou s’il repose sur des hypothèses fortes de lois. Mais ce n’était qu’une courte mise en bouche. Pour aller plus loin, je suggère la lecture de Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle, M.V. (2008) Sound and fury: McCloskey and significance testing in economics, Liberman, M. (2013). “Significance”, in 1885 and today et le passionnant Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Et sur les -value, Goodman, S. (2002A Dirty Dozen:Twelve P-Value Misconceptions. Je devrais aussi citer Ziliak, S. & McCloskey, D. (2008). The cult of statistical significance, mais je ne l’ai pas encore lu. Par contre, Xi’an en parlait sur son blog. Et si j’étais courageux, je tenterais une lecture bayésienne de ce billet…. promis, dès que je trouve un peu de temps !

Logistic regression and categorical covariates

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
     Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 ** 
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543    
X3D          0.08004    0.34916   0.229 0.818672    
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality . Which means that for someone with characteristics , we predict the following probability

where  denotes the cumulative distribution function of the logistic distribution

For someone with characteristics , we predict the following probability

For someone with characteristics , we predict the following probability

(etc.) Here, if we accept  (against ), it means that modality  cannot be considerd as different from .

A natural idea can be to change the reference modality, and to look at the -values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k])
+ reg = glm(Y~X1+X2+X3,family=binomial,data=db)
+ M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4]
+ } 
> M
             A            B            C            D            E
A           NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05
B 4.771853e-04           NA 4.841204e-04 0.0009474491 4.743636e-01
C 9.225428e-01 0.0004841204           NA 0.7506242347 9.194193e-05
D 8.186724e-01 0.0009474491 7.506242e-01           NA 1.730589e-04
E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589           NA

and if we simply want to know if the -value exceeds – or not – 5%, we get the following,

> M.TF = M>.05
> M.TF
      A     B     C     D     E
A    NA FALSE  TRUE  TRUE FALSE
B FALSE    NA FALSE FALSE  TRUE
C  TRUE FALSE    NA  TRUE FALSE
D  TRUE FALSE  TRUE    NA FALSE
E FALSE  TRUE FALSE FALSE    NA

The first column is obtained when  is the reference, and then, we see which parameter should be considered as null. The interpretation is the following:

  •  and  are not different from 
  •  is not different from 
  •  and  are not different from 
  •  and  are not different from 
  •  is not different from 

Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand)

> library(car)
> db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 ** 
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *  
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories,  and , let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.39439    1.02791  -4.275 1.91e-05 ***
X1           0.51378    0.11138   4.613 3.97e-06 ***
X2           0.20807    0.07234   2.876  0.00402 ** 
X3bisBE      1.94905    0.36852   5.289 1.23e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 398.31  on 996  degrees of freedom
AIC: 406.31

Number of Fisher Scoring iterations: 7

Here, all the categories are significant. So we do have a proper model.

Une application du test de Fisher, partie 2

Poursuivons un peu la discussion d’hier. En fait, si on poursuit, on arrive à regrouper des régions ensemble…. mais sans réelle cohérence car des régions proches au niveau électorale peuvent être éloignées géographiquement.

Si on veut que les régions regroupées soient proches, il faut une base avec des informations spatiales, comme savoir si deux départements se touchent, ou pas. Pour cela, on utilise la base suivante, qui contient

  • le numéro du département i
  • le numéro du département j
  • une variable indicatrice indiquant si les département se touche
  • la distance (en km) entre les centroïdes des départements
> france=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/departements-france.csv",
+ header=TRUE,sep=";")

On va nettoyer un peu, en enlevant les départements d’outre-mer,

> france0=france; i200=c(201,202)
> france0=france0[-c(which(france0$depi%in%i200),
+ which(france0$depj%in%i200)),]
> france0$ij=100*france0$depi+france0$depj

La base ressemble à ça,

> head(france0)
depi depj Cij Dij   ij
1   77   75   0  51 7775
2   78   75   0  39 7875
3   91   75   0  40 9175
4   92   75   1  10 9275
5   93   75   1  12 9375
6   94   75   1  13 9475

Plusieurs stratégies sont possible: regarde des départements voisins et voir si on peut les regrouper. Ou regarder des départements proches électoralement, et voir s’il sont proches géographiquement. Par exemple, si on ordonne les coefficients obtenus lors de la régression,

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)
> COEF=coefficients(summary(reg))
> dbcoef=data.frame(COEF[,1:3],Dpt=rownames(COEF))
> rownames(dbcoef)=rank(COEF[,1])
> COEFtrie=dbcoef[as.character(1:94),]
> COEFtrie[1:6,]
Estimate Std..Error   t.value    Dpt
1 35.69903  0.3246129 109.97414 DptD06
2 36.46840  0.3244806 112.39010 DptD67
3 36.67594  0.3896297  94.13025 DptD68
4 37.39283  0.3202112 116.77552 DptD83
5 39.92023  0.3921762 101.79156 DptD74
6 42.61811  0.6182921  68.92877 DptD10

On peut se demander si les quatre premiers départements pourraient être regroupés ensemble. Une des modalités sert de modalité de référence, et on regarde les trois autres par rapport à cette dernière

> sbelection$Dpt=relevel(sbelection$Dpt,"D67")
> reg0=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> COEF0=coefficients(summary(reg0))
> COEF0[c(1,7,67,82),]
Estimate Std. Error     t value   Pr(>|t|)
(Intercept) 36.4684029  0.3244806 112.3900971 0.00000000
DptD06      -0.7693736  0.4589784  -1.6762740 0.09369320
DptD68       0.2075400  0.5070493   0.4093093 0.68231515
DptD83       0.9244243  0.4558759   2.0277980 0.04258819

On fait ensuite une analyse de la variance,

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68","D83")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F   Pr(>F)
1  36114 2281072960
2  36111 2280191044  3    881917 4.6556 0.002954 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On rejette ici l’hypothèse que les quatre départements se comportent de la même manière. Mais si on tente de regrouper les trois premières, cette fois, on accepte l’hypothèse

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280476723
2  36111 2280191044  2    285679 2.2621 0.1041

Sauf que sur les trois départements, seulement deux se touchent,

> france0[france0$ij%in%c(6768,667,668),]
depi depj Cij Dij   ij
3510    6   67   0 533  667
3552   67   68   1  95 6768
3605    6   68   0 439  668

On peut essayer de parcourir les sorties de la régression, et regarder si on a des voisins. Pour chaque département, on regarde si le voisin est un voisin géographique,

> departement=as.numeric(substr(COEFtrie$Dpt,5,6))
> deptij=departement[1:(nrow(COEFtrie)-1)]*100+
> departement[(2):(nrow(COEFtrie))]
> voisins=france0[france0$ij%in%deptij,]
> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)
> COEFtrie$voisin=NA
> v=voisinage[as.character(deptij)]
> COEFtrie$voisin[1:(nrow(COEFtrie)-1)]=
> v[is.na(v)==FALSE]
> COEFtrie[1:10,]
Estimate Std..Error   t.value    Dpt voisin
1  35.69903  0.3246129 109.97414 DptD06      0
2  36.46840  0.3244806 112.39010 DptD67      1
3  36.67594  0.3896297  94.13025 DptD68      0
4  37.39283  0.3202112 116.77552 DptD83      0
5  39.92023  0.3921762 101.79156 DptD74      0
6  42.61811  0.6182921  68.92877 DptD10      0
7  42.77423  0.4397101  97.27826 DptD01      0
8  43.53395  0.4448212  97.86843 DptD84      0
9  44.38759  0.3897553 113.88579 DptD85      0
10 44.65361  0.4577710  97.54575 DptD51      1

On note qu’il y a plusieurs voisins. On peut ensuite faire un test de Fisher. S’il est positif, on fusionne les classes. Pour la régression, mais aussi dans la base des distances et des voisinages,

>  if(COEFtrie$voisin[i]==1){
+  sbelection$Dpt[sbelection$Dpt==substr(
+ COEFtrie$Dpt[i+vsn],4,6)]=substr(COEFtrie$Dpt[i],4,6)
+  france0$depi[france0$depi==as.numeric(
+  substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))
+  france0$depi[france0$depj==as.numeric(
+ substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))

Dans la base des départements, il y aura plusieurs paires. Si A et B ont fusionné, il suffit que C touche A ou B pour que C soit voisin des régions fusionnées. C’est pour cela que l’on utilise

> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)

pour savoir si les départements sont voisins. On peut alors répéter plusieurs fois l’opération.

Par exemple, on note que 28, 45 et 89 peuvent fusionner,

> sbelection$Dpt=relevel(sbelection$Dpt,"D45")
+ reg10=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ COEF10=coefficients(summary(reg10))
+ COEF10[c(1,28,88),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 45.9303788  0.4170114 110.1417762 0.0000000
DptD28       0.5643212  0.6594028   0.8558065 0.3921105
DptD89       0.9423212  0.7016026   1.3430982 0.1792486
+ sbelection$DptGroupe=sbelection$Dpt
+ sbelection$DptGroupe[sbelection$Dpt%in%
+ c("D28","D89")]="D45"
+ model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
+ model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280314007
2  36111 2280191044  2    122963 0.9737 0.3777

Géographiquement, ce sont effectivement des voisins

Mais peu de départements voisins peuvent effectivement être fusionnés.

Une application du test de Fisher, partie 1

Un petit (?) billet pour revenir sur un point évoqué mardi dernier en cours, sur l’utilisation du test de Fisher (et des techniques d’analyse de variance). Je vais juste évoquer rapidement un vieux billet sur l’histoire du 5% utilisé lorsque l’on calcule la -value (ce qui complètera la réponse faite en classe mardi lorsque la question a été posée). Mardi dernier, j’avais évoqué rapidement le test d’une contrainte linéaire sur les paramètres. On se demande ici si les deux modèles suivants sont vraiment différents, entre

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

et sa version contrainte

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

On cherche alors à tester l’hypothèse

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

contre l’hypothèse alternative

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

Pour ça, on peut utiliser ce qu’a proposé Ronald Fisher, à savoir d’utiliser la somme des carrés des résidus (qui correspondent – à un facteur multiplicatif près – à l’estimateur de la variance du bruit), et il a posé

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

Je parlais tout à l’heure de choix de modèle… on notera que l’on en est pas très loin car on peut réécrire cette expression en faisant apparaître le  des régressions, contraintes et noncontraintes.

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

Fisher a montré que sous l’hypothèse H0,

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

i.e. la loi de la statistique de Fisher est précisément une loi de Fisher, Pour aller plus loin sur l’analyse de la variance (ANOVA, analysis of variance), je renvois àFaraway (2005).

J’avais dit oralement qu’une application intéressant du test de nullité de plusieurs coefficient dans une régression était son utilisation pour des variables catégorielles, et pour faire du regroupement de modalité.

Considérons un cas d’école. On dispose des résultats des votes dans plusieurs villes, en France, lors des dernières élections présidentielles (ces bases sont en ligne sur http://www.data.gouv.fr/), opposant François Hollande (avec un H) et Nicolas Sarkozy (avec un S),

> election=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/
205cf26f15a43974fdb7d3b25397613d.csv",
+ sep=";",header=TRUE,dec=",")
> head(election)
Code.du.dpt Inscrits Votants VoixH X.VoixExpH VoixS X.VoixExpS
1           1      592     506   195      41.05   280      58.95
2           1      215     175    82      52.23    75      47.77
3           1     8208    6447  3009      50.41  2960      49.59
4           1     1152     949   340      38.16   551      61.84
5           1      105      83    35      44.30    44      55.70
6           1     1702    1461   630      46.46   726      53.54

On va nettoyer un peu la base

> sbelection=election[election$Code.du.département
+ %in%as.character(1:95),]

car il y a des départements d’outre mers et la Corse, qui a une codification particulière. Bref, on retient les 95 départements métropolitain, moins le 20ème qui est la Corse. On va commencer par créer une variable de département qui soit un facteur qualitatif. Pour l’instant, le département est un numéro, entre 1 et 95 mais ce n’est aucunement un nombre !

> sbelection$Dpt=(paste("D",
> sbelection$Code.du.département,sep=""))
> I=as.numeric(as.character(sbelection$Code.du.département))<10
> sbelection$Dpt[I]=(paste("D0",
+ sbelection$Code.du.département,sep=""))
> sbelection$Dpt=factor(sbelection$Dpt)

On peut alors faire une régression, où on peut essayer de modéliser le taux de vote pour François Hollande (par exemple) en fonction de variables explicatives, et en particulier, on veut voir s’il n’y a pas un effet spatial. Bref, on va régresser sur le département. Et d’ailleurs que sur le département, inutile de mettre d’autre variable pour illustrer.

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)

Par contre, je fais une régression pondérée, car les villes (et villages) ne font pas du tout la même taille. C’est exactement la même idée que de faire une moyenne pondérée pour avoir le taux moyen dans une région (on rediscutera des régressions pondérées en fin de cours).

> COEF=coefficients(summary(reg))

La sortie ressemble à

> COEF
Estimate Std. Error   t value Pr(>|t|)
DptD01 42.77423  0.4397101  97.27826        0
DptD02 52.39423  0.4572016 114.59765        0
DptD03 56.89434  0.5468760 104.03518        0
DptD04 51.06544  0.7807748  65.40355        0
DptD05 50.70525  0.8403979  60.33481        0
DptD06 35.69903  0.3246129 109.97414        0
DptD07 53.44185  0.5564489  96.04089        0
DptD08 51.86319  0.6360901  81.53434        0
DptD09 64.67228  0.8013389  80.70528        0
DptD10 42.61811  0.6182921  68.92877        0

[…]

DptD89 46.87270  0.5642231  83.07476        0
DptD90 50.47561  0.9063735  55.68964        0
DptD91 53.41748  0.3172310 168.38670        0
DptD92 49.54305  0.2850597 173.79888        0
DptD93 65.29062  0.3338954 195.54217        0
DptD94 56.48855  0.3196643 176.71211        0
DptD95 53.89382  0.3344373 161.14778        0

On peut d’ailleurs faire une prédiction, par département (en prenant comme prédiction celle obtenue pour la première ville du département)

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)

On peut ensuite essayer de faire un petite dessin pour visualiser cette prédiction,

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)
> nomdpt=read.table(
+ "http://freakonometrics.blog.free.fr/
public/data/departements.txt",sep=";")
> nomdpt$nom=tolower(nomdpt$V2)
> dpt=as.character(c(paste("0",
+ 1:9,sep=""),c(10:19,21:95)))
> noms=as.character(nomdpt[nomdpt$V3%in%dpt,2])
> couleurliste=rev(rainbow(n=100))
> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur=couleurliste[trunc((PredDept-35)*3)]

C’est un peu long (et compliqué) car il faut identifier le bon département sur la carte R (qui est codée en lettres). D’ailleurs la fonction suivante ne marche pas

> match=match.map(francemap,dpt)
> francemap$names[which(is.na(match.map(
+ francemap,dpt[35]))==FALSE)]
[1] "Indre-et-Loire" "Indre"

car en l’occurrence, seule l’Indre aurait du être obtenue. On va faire la correspondance à la main,

> match=rep(NA,length(francemap$names))
> for(i in 1:length(francemap$names)){
+ if(sum(dpt==francemap$names[i])>0)
+ {match[i]=which(dpt==francemap$names[i])}
+ }

Une fois les départements identifiés, on colorie,

> color=couleur[match]
> map(database="france", fill=TRUE, col=color)

C’est joli, mais comme toujours c’est un peu compliqué de faire des cartes…

A titre d’information, les couleurs correspondent aux valeurs suivantes (avec en abscisse, le taux de vote obtenu par François Hollande)

Bref, on a beaucoup de modalité car il y a 95 département. Mais on devrait pouvoir regrouper. Si on ne sait pas trop par où commencer, prenons un département au hasard, disons le 35 (i.e. Rennes),

> k=34

On peut ensuite regarder les estimations des coefficients, ou mieux, les intervalles de confiances des différents estimateurs, et voir ceux qui chevauche celui de Rennes (on va prendre des intervalles de confiance à 70% histoire d’avoir des chances que les deux coefficients soient vraiment proches)

> alpha=.7
> CI.COEF=confint(reg,level=alpha)

On peut visualiser ces intervalles, avec une bande horizontale pour Rennes,

> plot(c(1:19,21:95),COEF[,1],ylim=c(35,
+ 60),pch=3,cex=.5)
> segments(c(1:19,21:95),CI.COEF[,1],
+ c(1:19,21:95),CI.COEF[,2])
> lk=which((CI.COEF[,1]<CI.COEF[k,2])&
+ (CI.COEF[,2]>CI.COEF[k,1]))
> segments(c(1:19,21:95)[lk],CI.COEF[lk,1],
+ c(1:19,21:95)[lk],CI.COEF[lk,2],col="red",lwd=2)

et en rouge, les intervalles qui la chevauchent

Ça y est, on a une liste potentielle de six départements à regrouper… On va prendre comme modalité de référence le département qui nous intéresse, et voir si les coefficients des autres départements sont (ou pas) significatifs.

> sbelection$Dpt=relevel(sbelection$Dpt,"D35")
> reg35=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)
> COEF35=coefficients(summary(reg35))
> COEF35[c(1,lk[lk<k]+1,lk[lk>k]),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 55.66226653  0.3275407 169.9399801 0.0000000
DptD11       0.59043789  0.6272106   0.9413710 0.3465210
DptD32       0.96845567  0.7855154   1.2328920 0.2176241
DptD36      -0.03007518  0.7414655  -0.0405618 0.9676455
DptD62       0.56238714  0.4247897   1.3239189 0.1855384
DptD75      -0.06226653  0.4098013  -0.1519432 0.8792326
DptD81      -0.11970183  0.6067040  -0.1972986 0.8435950

On a confirmation que ces six départements sont non significativement (individuellement) différents de celui de référence (et donc on accepte l’hypothèse de nullité du coefficient). Mais peut-on pour autant les regrouper? C’est là que l’on peut faire un test de Fisher. Pour cela on fait deux régressions. La première est celle avec un regroupement de départements en une modalité unique (les sept sont désormais ensemble)

> sbelection$DptGroupe=sbelection$Dpt
> sbelection$DptGroupe[sbelection$Dpt%
+ in%substr(rownames(COEF[lk,]),4,6)]=paste("D",k+1,sep="")
> model0=lm(X.VoixExpH~DptGroupe,
+ weights=Votants,data=sbelection)

et la seconde est celle que l’ont vient de faire (qui sera notre hypothèse alternative, i.e. il ne faut pas regrouper les départements ensemble car ils sont différents – ou au moins l’un est différent)

> model1=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)

On peut alors faire le test de Fisher (en continuant à corriger de l’effet taille de la population, et en faisant des sommes pondérées)

> SCR0=sum(sbelection$Votants*residuals(model0)^2)
> SCR1=sum(sbelection$Votants*residuals(model1)^2)
> dl0=model0$df.residual
> dl1=model1$df.residual
> (F=(SCR0-SCR1)/(dl0-dl1)/SCR1*dl1)
[1] 0.917397

On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value

> 1-pf(F,dl0-dl1,dl1)
[1] 0.4809411

On a fait (sans le savoir) une analyse de la variance, qui peut se faire plus simplement (?) en utilisant la fonction suivante

> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36117 2280538612
2  36111 2280191044  6    347568 0.9174 0.4809

Et on accepte ici l’hypothèse que ces sept départements peuvent être regroupé en un seul.

Enfin, si statistiquement, il y a du sens à regrouper les modalités, on peut se demander si ce regroupement de régions géographique à du sens,

> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur="red"
> match=match.map(francemap,dpt)
> color=couleur[match]
> map(dataLa première  estbase="france", fill=TRUE, col=color)

Les régions que l’on souhaite associer à Rennes sont fort éloignées, et les regrouper n’a peut-être pas de sens… (à suivre donc).

Fisher-Tippett theorem with an historical perspective

A couple of weeks ago, Rafael asked me if I had something on the history of extreme value theory. Since I will get back to fundamental results about extremes in my course, I promised I will write down a short post on all that.
To start from the beginning, in 1928, Ronald Fisher and Leonard Tippett formulated the three types of limiting distributions for the maximum term of a random sample (Fisher & Tippett (1928)). The problem was to characterize function http://freakonometrics.blog.free.fr/public/perso5/ext-01.gif such that

http://freakonometrics.blog.free.fr/public/perso5/ext-2.gif

where http://freakonometrics.blog.free.fr/public/perso5/ext-3.gif where http://freakonometrics.blog.free.fr/public/perso5/ext-4.gif‘s are i.i.d. with cumulative distribution function http://freakonometrics.blog.free.fr/public/perso5/ext-5.gif. They had supporting arguments, but no (rigorous) proof. Nevertheless, the obtained that the only possible types for G were

http://freakonometrics.blog.free.fr/public/perso5/ext-6.gif

i.e. Fréchet type (Pareto-type tails), or

http://freakonometrics.blog.free.fr/public/perso5/ext-7.gif

i.e. Weibull type (bounded distribution type), or

http://freakonometrics.blog.free.fr/public/perso5/ext-8.gif

i.e. Gumbel type (exponential-type tails). Emil Gumbel has been intensively using the so-called Gumbel distribution on river flows, since (as he explained in 1958), “it seems that the rivers know the theory. It only remains to convince the engineers of the validity of this analysis“.
Independently of that work (published in 1928), Maurice Fréchet considered in 1927 (in Sur la loi de probabilité de l’écart maximum) possible limits of

http://freakonometrics.blog.free.fr/public/perso5/ext-9.gif

and obtained only http://freakonometrics.blog.free.fr/public/perso5/ext-10.gif as possible limit. Richard von Mises gave in 1936 sufficient, but not necessary conditions for their (max) domain of attraction, i.e. characterization of function http://freakonometrics.blog.free.fr/public/perso5/ext-11.gif such that the maxima converges to some specific function http://freakonometrics.blog.free.fr/public/perso5/ext-01.gif (von Mises (1936)). E.g. he noticed that a sufficient condition on http://freakonometrics.blog.free.fr/public/perso5/ext-11.gifto be in the (max) domain of attraction of the Gumbel distribution is that

http://freakonometrics.blog.free.fr/public/perso5/ext-13.gif

Then in 1943, Boris Gnedenko gave a complete characterization of those three types, with a complete characterization for two of them (heavy tails, i.e. Fréchet type and bounded support, i.e. Weibull) but his necessary and sufficient condition was based on a function that was not explicitly defined (see Gnedenko (1943)). Laurens de Haan in the 70’s derived checkable condition for Gumbel’s type.
Boris Gnedenko proved (in Section 4 of his paper) that F is the (max) domain of attraction of http://freakonometrics.blog.free.fr/public/perso5/ext-10.gif if and only if http://freakonometrics.blog.free.fr/public/perso5/ext-16.gif is regularly varying at infinity, with index http://freakonometrics.blog.free.fr/public/perso5/ext-17.gif (even if the term “regular variation” was not mentioned in the paper). Similar results were derived to characterize functions in the (max) domain of attraction of Weibull. For the (max) domain of attraction of http://freakonometrics.blog.free.fr/public/perso5/ext-18.gif, Boris Gnedenko obtained that a necessary and sufficient condition was that there exists a function http://freakonometrics.blog.free.fr/public/perso5/ext-19.gif such http://freakonometrics.blog.free.fr/public/perso5/ext-19.gif goes to 0 at infinity and

http://freakonometrics.blog.free.fr/public/perso5/ext-20.gif

Several papers have discussed what function http://freakonometrics.blog.free.fr/public/perso5/ext-19.gif could be e.g. David Mejzler in 1949 (in Russian, but see also his 1965 paper), and Laurens de Hann in 1970 and 1971 (following the dramatic flood in the Netherlands in 1953, researchers in the Netherlands have focuses on dikes, and extreme value applications).

Mejzler’s idea was to work on quantiles, and not on the cumulative distribution function. I.e. define

http://freakonometrics.blog.free.fr/public/perso5/ext-21.gif

Then a necessary and sufficient condition for F to be in the (max) domain of attraction of http://freakonometrics.blog.free.fr/public/perso5/ext-18.gif is that

http://freakonometrics.blog.free.fr/public/perso5/ext-23.gif

Laurens de Haan proved in 1971 that function http://freakonometrics.blog.free.fr/public/perso5/ext-19.gif can be – in general – given by

http://freakonometrics.blog.free.fr/public/perso5/ext-25.gif

And in 1976, Laurens de Haan obtained a three-type convergence working on quantile function http://freakonometrics.blog.free.fr/public/perso5/ext-26.gif (with a much shorter proof).
There have been many many papers extending Fisher-Tippett’s theorem, e.g. on non-independent sequences, like exchangeable ones (in a paper by Simeon Berman in 1962, or on stationary Gaussian sequences in 1964).

Fisher-Tippett theorem and limiting distribution for the maximum

Tomorrow, we will discuss Fisher-Tippett theorem. The idea is that there are only three possible limiting distributions for normalized versions of the maxima of i.i.d. samples http://freakonometrics.blog.free.fr/public/perso5/max-00.gif. For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. http://freakonometrics.blog.free.fr/public/perso5/max-09.gif on the unit interval. Let http://freakonometrics.blog.free.fr/public/perso5/max-10.gif and http://freakonometrics.blog.free.fr/public/perso5/max-11.gif. Then, for all http://freakonometrics.blog.free.fr/public/perso5/max-12.gif and http://freakonometrics.blog.free.fr/public/perso5/max-13.gif,

http://freakonometrics.blog.free.fr/public/perso5/max-14.gif

i.e. the limiting distribution of the maximum is Weibull’s.

set.seed(1)
s=1000000
n=100
M=matrix(runif(s),n,s/n)
V=apply(M,2,max)
bn=1
an=1/n
U=(V-bn)/an
hist(U,probability=TRUE,,col="light green",
xlim=c(-7,1),main="",breaks=seq(-20,10,by=.25))
u=seq(-10,0,by=.1)
v=exp(u)
lines(u,v,lwd=3,col="red")

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function http://freakonometrics.blog.free.fr/public/perso5/max-05.gif. Let http://freakonometrics.blog.free.fr/public/perso5/max-06.gif and http://freakonometrics.blog.free.fr/public/perso5/max-07.gif, then

http://freakonometrics.blog.free.fr/public/perso5/max-08.gif

which means that the limiting distribution is Fréchet’s.

set.seed(1)
s=1000000
n=100
M=matrix((runif(s))^(-1/2),n,s/n)
V=apply(M,2,max)
bn=0
an=n^(1/2)
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(0,7),main="",breaks=seq(0,max(U)+1,by=.25))
u=seq(0,10,by=.1)
v=dfrechet(u,shape=2)
lines(u,v,lwd=3,col="red")

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function http://freakonometrics.blog.free.fr/public/perso5/max-01.gif. Let http://freakonometrics.blog.free.fr/public/perso5/max-02.gif and http://freakonometrics.blog.free.fr/public/perso5/max-03.gif, then

http://freakonometrics.blog.free.fr/public/perso5/max-04.gif

i.e. the limiting distribution for the maximum is Gumbel’s distribution.

library(evd)
set.seed(1)
s=1000000
n=100
M=matrix(rexp(s,1),n,s/n)
V=apply(M,2,max)
(bn=qexp(1-1/n))
log(n)
an=1
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Consider now a Gaussian http://freakonometrics.blog.free.fr/public/perso5/max-17.gif sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

http://freakonometrics.blog.free.fr/public/perso5/max-15.gif

as http://freakonometrics.blog.free.fr/public/perso5/max-16.gif. Let http://freakonometrics.blog.free.fr/public/perso5/max-18.gif and http://freakonometrics.blog.free.fr/public/perso5/max-19.gif. Then we can get

http://freakonometrics.blog.free.fr/public/perso5/max-20.gif

as http://freakonometrics.blog.free.fr/public/perso5/max-21.gif. I.e. the limiting distribution of the maximum of a Gaussian sample is Gumbel’s. But what we do not see here is that for a Gaussian sample, the convergence is extremely slow, i.e., with 100 observations, we are still far away from Gumbel distribution,

and it is only slightly better with 1,000 observations,

set.seed(1)
s=10000000
n=1000
M=matrix(rnorm(s,0,1),n,s/n)
V=apply(M,2,max)
(bn=qnorm(1-1/n,0,1))
an=1/bn
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Even worst, consider lognormal observations. In that case, recall that if we consider (increasing) transformation of variates, we are in the same domain of attraction. Hence, since http://freakonometrics.blog.free.fr/public/perso5/max-22.gif, if

http://freakonometrics.blog.free.fr/public/perso5/max-23.gif

then

http://freakonometrics.blog.free.fr/public/perso5/max-24.gif

i.e. using Taylor’s approximation on the right term,

http://freakonometrics.blog.free.fr/public/perso5/max-25.gif

This gives us normalizing coefficients we should use here.

set.seed(1)
s=10000000
n=1000
M=matrix(rlnorm(s,0,1),n,s/n)
V=apply(M,2,max)
bn=exp(qnorm(1-1/n,0,1))
an=exp(qnorm(1-1/n,0,1))/(qnorm(1-1/n,0,1))
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,40,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Confidence interval for predictions with GLMs

Consider a (simple) Poisson regression https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss01.gif. Given a sample https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss02.gif where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss03.gif, the goal is to derive a 95% confidence interval for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss04.gif given https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss05.gif, where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss04.gif is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below

> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1$fit,pch=4,lwd=3)

i.e.

Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif denote the maximum likelihood estimator of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss07.gif. Then
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss40.gif
where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss101.gif is Fisher information of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif (from standard maximum likelihood theory). Recall that
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss13.gif
where computation of those values is based on the following calculations
http://freakonometrics.blog.fre<br /><br /> e.fr/public/latex/poiss21.gif
In the case of the log-Poisson regression
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss36.gif
Let us get back to our initial problem.

  • confidence interval for the linear combination

A first idea to get a confidence interval for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss49.gif is to get a confidence interval for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss100.gif (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss40.gif

thus, an approximation for the variance matrix of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif will be based on https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss45.gif, obtained by plugging estimators of the parameters.
Then, since https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss47.gif has a normal distribution, centered on https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss49.gif, with variance https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss102.gif where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/Poiss110.gif is the variance of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/poiss06.gif. All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response"))
> i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05

Based on those values, it is easy to derive the standard deviation for the linear combination,

> x=30
> P2=predict(r,type="link",se.fit=TRUE,
+ newdata=data.frame(speed=x))
> P2
$fit
1
5.046034

$se.fit
[1] 0.05747075

$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

  • delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> P1
$fit
1
155.4048

$se.fit
1
8.931232

$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101
  • bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,

Tests statistiques et p value

Mercredi, nous avons abordé le concept de http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value en cours (avant de parler tout à l’heure de tests asymptotiques). Intuitivement, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité que la statistique de test prenne une valeur au moins aussi extrême que celle qui a été observée, si l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif est vraie.

  • Fisher vs. Neyman (et Pearson)

Dans l’approche de Fisher, on suppose que http://freakonometrics.blog.free.fr/public/perso2/testp-01.gif et on cherche à tester http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif. On dispose pour cela d’un échantillon http://freakonometrics.blog.free.fr/public/perso2/testp-02.gif.
On choisit une statistique http://freakonometrics.blog.free.fr/public/perso2/testp-03.gif telle que plus http://freakonometrics.blog.free.fr/public/perso2/testic-5.gif est grande, plus on s’éloigne de http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif (on calcule une distance à l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif, quelque chose comme ça). On calcule alors la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme

http://freakonometrics.blog.free.fr/public/perso2/testp-07.gif

et on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif si http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif est trop petite. La notion de “plus petite” est liée à la comparaison avec un seuil prédéfini, souvent 5% (cf plus loin pour une discussion).
C’est un peu heuristique, mais c’est comme cela que ça fonctionne.
Dans l’approche de Neyman-Pearson, on cherche à testerhttp://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contre une hypothèse alternative, du genre http://freakonometrics.blog.free.fr/public/perso2/testp-11.gif. Si http://freakonometrics.blog.free.fr/public/perso2/testp-12.gif, on va construire une région de rejet de la forme

http://freakonometrics.blog.free.fr/public/perso2/testp-13.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est choisi en fonction d’une probabilité d’erreur de première espèce souhaitée a priori, i.e.

http://freakonometrics.blog.free.fr/public/perso2/testp-15.gif

Par exemple, considérons un échantillon suivant une loi http://freakonometrics.blog.free.fr/public/perso2/testp-16.gif. On dispose de http://freakonometrics.blog.free.fr/public/perso2/testp-17.gif observations pour tester http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif. L’échantillon est le suivant

> set.seed(4)
> (X=rnorm(10))
 [1]  0.2167549 -0.5424926  0.8911446  0.5959806  1.6356180  0.6892754
 [7] -1.2812466 -0.2131445  1.8965399  1.7768632

La statistique proposée par Fisher (je ne discuterais pas le choix optimal de la statistique – si quelqu’un a des références, les commentaires sont ouvertes) est la moyenne normalisée,

http://freakonometrics.blog.free.fr/public/perso2/testp-21.gif

On sait que sous H_0, cette statistique suit une loi normale centrée réduite.  La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value se calcule alors simplement

> (T=sqrt(10)*mean(X))
[1] 1.791523
 
> 2*(1-pnorm(T,mean=0,sd=0)))
[1] 0.07320942

Graphiquement, cela se représente de la manière suivante.

Cette probabilité étant faible, on aurait tendance vouloir rejeter http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.
Dans l’optique de Neyman Pearson, si on teste http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif contre http://freakonometrics.blog.free.fr/public/perso2/testp-22.gif, la région critique s’obtient à l’aide du rapport de vraisemblance

http://freakonometrics.blog.free.fr/public/perso2/testp-25.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso2/testp-26.gif

ce qui se traduit par

http://freakonometrics.blog.free.fr/public/perso2/testp-27.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est tel que

http://freakonometrics.blog.free.fr/public/perso2/testp-29.gif

Or sous H_0, la somme suit une loi normale centrée de variance http://freakonometrics.blog.free.fr/public/perso2/testp-30.gif. Ici, http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est donné par

> qnorm(1-alpha,mean=0,sd=sqrt(10))
[1] 5.201484

Or pour rappel, la somme vaut ici

> sum(X)
[1] 5.665293

On est dans la région de rejet: on rejette alors http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.

  • p-value pour les tests usuels de R

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est intéressante en pratique car elle permet à l’utilisateur de ne pas rentrer trop dans les détails de la construction du test, et de la loi de la statistique considérée.
Par exemple (comme celui abordé ici), pour un test d’égalité de proportion, sous R, la commande est

> prop.test(nx,n,.5)
 
1-sample proportions test with continuity correction
 
data:  nx out of n, null probability 0.5
X-squared = 0.05, df = 1, p-value = 0.823
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.3204804 0.7617145
sample estimates:
p
0.55

Et de manière générale, tous les tests sous R renvoient une http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value.

  • p ne signifie pas puissance

(mais probabilité). Dans un test simple,http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contrehttp://freakonometrics.blog.free.fr/public/perso2/testh13.gif, la puissance est la probabilité de rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif sous l’hypothèse où http://freakonometrics.blog.free.fr/public/perso2/H1.gif est vraie (rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif avec raison si on veut). C’est donc une fonction de  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.
Pour un test multiple, avec une hypothèse alternative de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-36.gif, on définit alors une fonction puissance pour toute valeur  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.  On a alors des graphiques de fonction puissance, comme ceux évoqués ici. Bref, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est une valeur, alors que la puissance sera une fonction, définie sur http://freakonometrics.blog.free.fr/public/perso2/testp-40.giflorsque H_1 est une hypothèse de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-41.gif.

  • L’interprétation de p

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité d’obtenir (au moins) la statistique obtenue pour l’échantillon si http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie. Classiquement, on rejette http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif si la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est inférieure au seuil http://freakonometrics.blog.free.fr/public/perso2/testh04.gif .
Certains interprètent la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme la probabilité que http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, http://freakonometrics.blog.free.fr/public/perso2/testp-55.gif. Ce qui n’a pas de sens si l’on s’en tient au formalisme énoncé en cours: http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, ou pas.
Ce n’est pas non plus la probabilité de rejeter, à tort, http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif , http://freakonometrics.blog.free.fr/public/perso2/testp-57.gif, correspondant à la puissance du test (cf la discussion dans le paragraphe précédant).

  • Pourquoi 5% ?

En pratique, on prend toujours comme valeur critique 5%, avec la règle décision simple suivante:

  • http://freakonometrics.blog.free.fr/public/perso2/testp-65.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-61.gif on accepte http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • http://freakonometrics.blog.free.fr/public/perso2/testp-66.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-60.gif on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-33.gif

C’est Fisher qui semble être à l’origine de ce 5%, en 1925 “The value for which P=0.05, or 1 in 20, is 1.96 or nearly 2; it is convenient to take this point as a limit in judging whether a deviation ought to be considered significant or not. Deviations exceeding twice the standard deviation are thus formally regarded as significant. Using this criterion we should be led to follow up a false indication only once in 22 trials, even if the statistics were the only guide available.” On retrouve ce “un sur vingt” un an plus tard, “Either there is something in the treatment, or a coincidence has occurred such as does not occur more than once in twenty trials”.
Voilà ce que l’on retient souvent de la lecture des conseils de Fisher.
Maintenant, il lui arrive aussi, à l’occasion, de penser que 4% ce n’est pas significatif, e.g. “P is between .02 and .05. The result must be judged significant, though barely so; in view of the data we cannot ignore the possibility that on this field, and in conjunction with the other manures used, nitrate of soda has conserved the fertility better than sulphate of ammonia; the data do not, however, demonstrate this point beyond the possibility of doubt.” Ou parfois pense que 8% est significatif “P=.089. Thus a larger value of 2 would be obtained by chance only 8.9 times in a hundred, from a series of values in random order. There is thus some reason to suspect that the distribution of rainfall in successive years is not wholly fortuitous, but that some slowly changing cause is liable to affect in the same direction the rainfall of a number of consecutive years.”

Statistiques, de la théorie à la pratique (partie 2)

Il y a une quinzaine de jours (ici) on avait vu comment construire des estimateurs pour l’indice de queue de la loi de Pareto. On avait 4 valeurs (différentes). Maintenant, on peut se demander lequel devrait être le “meilleur“.

  • propriétés de l’estimateur du maximum de vraisemblance

Maintenant que l’on a un peu avancé dans le cours, on peut continuer l’étude de la distribution des couts de sinistres (étude commencée ici).

http://freakonometrics.blog.free.fr/public/perso2/pareto-100.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/pareto-101.gif

de telle sorte que

http://freakonometrics.blog.free.fr/public/perso2/pareto-102.gif

et

http://freakonometrics.blog.free.fr/public/perso2/pareto-204.gif

Aussi l’information de Fisher est http://freakonometrics.blog.free.fr/public/perso2/pareto-203.gif.
La variance aymptotique de http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gif est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-206.gif
  • propriétés de l’estimateur de la méthode des moments

Nous avions noté que

http://freakonometrics.blog.free.fr/public/perso2/pareto-207.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-208.gif, alors http://freakonometrics.blog.free.fr/public/perso2/pareto-209.gif. Aussi, comme le théorème central limite limite

http://freakonometrics.blog.free.fr/public/perso2/pareto-210.gif

http://freakonometrics.blog.free.fr/public/perso2/pareto-211.gif

la delta method permet d’écrire, en posant http://freakonometrics.blog.free.fr/public/perso2/pareto-311.gif,

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso2/correc-latex-pareto2.gif

La variance asymptotique est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-310.gif
  • propriétés de l’estimateur de la méthode de la médiane

Pour le dernier estimateur, le théorème de Glivenko Cantelli sur les fonctions quantiles permet d’écrire

http://freakonometrics.blog.free.fr/public/perso2/pareto-309.gif

i.e. pour la médiane,

http://freakonometrics.blog.free.fr/public/perso2/pareto-308.gif

où la variance asymptotique (pour la médiane) peut se simplifier sous la forme

http://freakonometrics.blog.free.fr/public/perso2/pareto-307.gif

Là encore, en utilisant la delta method (je peux renvoyer ici pour une application numérique de cette méthode), avec

http://freakonometrics.blog.free.fr/public/perso2/pareto-306.gif

i.e. http://freakonometrics.blog.free.fr/public/perso2/pareto-305.gif, on en déduit que

http://freakonometrics.blog.free.fr/public/perso2/pareto-304.gif

La variance asymptotique de notre estimateur est alors (en enlevant le facteur http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gifau dénominateur – car il est commun à tous),

http://freakonometrics.blog.free.fr/public/perso2/pareto-303.gif

Pour trois de nos estimateurs (je n’ai pas fait le quatrième, ça serait hors programme), on a pu obtenir facilement des variances asymptotiques. Graphiquement, on a les évolutions suivantes pour les variances asymptotiques de nos trois estimateurs, en fonction de la valeur de http://freakonometrics.free.fr/blog/latex/information09.gif, avec le maximum de vraisemblance (en rouge), celui de la méthode des moments (en bleu) et celui basé sur la médiane (en mauve)

Manifestement le premier estimateur (le maximum de vraisemblance) a une variance asymptotique plus faible que les deux autres (mais on l’avait vu en cours, c’est ce que garantit l’inégalité de Cramér-Rao). Autrement dit, l’estimateur du maximum de vraisemblance est probablement le plus fiable….
En utilisant ces résultats mentionnés, on peut aussi construire des intervalles de confiance asymptotiques pour la valeur du paramètre http://freakonometrics.free.fr/blog/latex/information09.gif, à l’aide de nos trois estimateurs, i.e. pour un intervalle de confiance à 90%,

http://freakonometrics.blog.free.fr/public/perso2/pareto-302.gif

où les estimations des variances asymptotiques sont obtenues en substituant http://freakonometrics.blog.free.fr/public/perso2/pareto-205.gifà  http://freakonometrics.free.fr/blog/latex/information09.gif dans les formules au dessus.
Si on considère que l’estimateur du maximum de vraisemblance est le meilleur que l’on puisse utiliser, on obtient comme intervalle de confiance pour http://freakonometrics.free.fr/blog/latex/information09.gif,

2.845093-1.64/sqrt(100)*2.845093
[1] 2.378498
 
2.845093+1.64/sqrt(100)*2.845093
[1] 3.311688

L’étape d’après sera d’introduire des tests. On pourra ainsi se demander si l’on peut accepter l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/pareto-301.gif. Mais c’est une autre histoire… à suivre donc….

Information de Fisher

Nous avons défini dans le cours l’information de Fisher par

http://freakonometrics.free.fr/blog/latex/information01.gif

où http://freakonometrics.free.fr/blog/latex/information02.gif. On verra la semaine prochaine que cette information peut s’écrire (moyennant quelques conditions de régularité)

http://freakonometrics.free.fr/blog/latex/information03.gif

Dans un cadre multivarié, on a que http://freakonometrics.free.fr/blog/latex/information04.gif avec

http://freakonometrics.free.fr/blog/latex/information05.gif

ou encore

http://freakonometrics.free.fr/blog/latex/information06.gif

(avec une écriture analogue au cas univarié en utilisant la matrice Hessienne).
Étant donné un échantillon, l’information de Fisher empirique est la pente de la log-vraisemblance en https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/information09.gif, i.e. dans le cas d’échantillons de 20 tirages de pile/face (avec une pièce non truquée)

http://freakonometrics.blog.free.fr/public/perso2/info-fisher.gif
X=sample(c("PILE","FACE"),size=20,prob=c(.5,.5),replace=TRUE)
X=(X=="PILE")*1
vecP=seq(0,1,by=.01)
logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
L=lapply(vecP,logL)
plot(vecP,L,type="l",col="blue",lwd=2,xlab="",ylab="")
abline(v=.5,col="red",lty=2)

En répétant plusieurs milliers de fois, on peut étudier la distribution de cette pente

En note que le score est centré, et ici la variance est

> var(SCORE)
[1] 80.00168
> 20/(.5*(1-.5))
[1] 80
20/(.5*(1-.5))

(ce qui correspond effectivement à l’information de Fisher). Mais cette informationest un concept plus général qu’une variance…
Bon, maintenant essayons de creuser un peu cette histoire d’information (c’est Ronald Fisher qui a – le premier – utilisé ce nom en 1925, bien avant les travaux de Shannon et Wiener sur des mesures logarithme d’information dans un contexte de théorie de la communication).

On cherche a quantifier le mesure d’information d’un échantillon http://freakonometrics.free.fr/blog/latex/information07.gif correspondant à une collection de variables aléatoires http://freakonometrics.free.fr/blog/latex/information08.gif, i.i.d. de loi http://freakonometrics.free.fr/blog/latex/information10.gif, et de densité/loi de probabilité http://freakonometrics.free.fr/blog/latex/information11.gif.

Une hypothèse est qu’avoir http://freakonometrics.free.fr/blog/latex/information12.gif observations apporte deux fois plus d’information qu’en avoir seulement http://freakonometrics.free.fr/blog/latex/information13.gif. On a alors quelque chose qui pourrait s’écrire

http://freakonometrics.free.fr/blog/latex/information14.gif

Aussi on suppose que pas d’observation signifie pas d’information… J’ai l’impression que les bayésiens ne seraient pas d’accord, mais c’est une autre histoire…

On suppose aussi que l’information dépend de la loi  http://freakonometrics.free.fr/blog/latex/information11.gif, et pas de l’échantillon http://freakonometrics.free.fr/blog/latex/information07.gif.

On peut aussi noter (c’est une idée peut être plus simple à comprendre dans le contexte de la communication) que la compression ne peut que détériorer l’information. Autrement dit, si on utiliser http://freakonometrics.blog.free.fr/public/perso2/infofi-01.gif et non pas http://freakonometrics.blog.free.fr/public/perso2/infori-02.gif, l’information est moins bonne, i.e.

http://freakonometrics.blog.free.fr/public/perso2/infofi-03.gif

Et le seul cas pour lequel on n’a pas de perte d’information est le cas d’une statistique exhaustive.

On notera aussi qu’une information doit être positive (au pire on n’a pas d’information).

Mais pour relier cette information a une propriété sur la variance de notre estimateur, il faut creuser un peu. Considérons un estimateur sans biais http://freakonometrics.free.fr/blog/latex/information15.gif de http://freakonometrics.free.fr/blog/latex/information09.gif, i.e.

http://freakonometrics.free.fr/blog/latex/infromation16.gif

ou encore (en supposant que l’échantillon soit réduit à une unique observation)

http://freakonometrics.free.fr/blog/latex/information17.gif

soit, en notant que

http://freakonometrics.free.fr/blog/latex/information18.gif

(car on a une densité)

http://freakonometrics.free.fr/blog/latex/information19.gif

En différenciant par rapport à http://freakonometrics.free.fr/blog/latex/information09.gif, on peut en déduire que

http://freakonometrics.free.fr/blog/latex/information20.gif

soit (en admettant que l’on puisse intervertir intégrale et dérivée)

http://freakonometrics.free.fr/blog/latex/information21.gif

Or on peut réécrire

http://freakonometrics.free.fr/blog/latex/information22.gif

ce qui donne

http://freakonometrics.free.fr/blog/latex/information23.gif

ou encore

http://freakonometrics.free.fr/blog/latex/information24.gif

Et là, on invoque l’inégalité de Cauchy-Schwarz afin de garantir que

http://freakonometrics.free.fr/blog/latex/information25.gif

On reconnait alors deux espérances, i.e.

http://freakonometrics.free.fr/blog/latex/information26.gif

(pour le terme de gauche, on a utilisé le fait qu’on avait un estimateur sans biais). Bref, il existe plein d’application de cette information, cette dernière inégalité étant probablement la plus importante.
Pour aller un peu plus loin, je peux renvoyer vers des notes de cours de Tran Viet Chi, en ligne ici. Et pour ceux qui veulent aller vraiment plus loin, John Pierce Robert Ash et surtout Solomon Kullback ont écrit des livres entiers sur statistiques et mesure d’information….

Khi de la poule ou de l’oeuf…

Samedi soir, Djalil et Florent se demandaient ce que Karl Pearson avait inventé en premier, la loi du khi (ou du chi), ou la loi du khi-deux (pensant que la loi du chi-deux était relativement naturelle, et que la loi du khi découlait de la loi du khi-deux). La loi du khi admet pour densité

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

et est aussi parfois définie à partir de la fonction Gamma incomplète, i.e. sa fonction de répartition est

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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/gamma-reg-1.png

avec

https://perso.univ-rennes1.fr/arthur.charpentier/latex/gamma-reg-2.png
et classiquement,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/gamma-reg-3.png

Pour la loi du khi-deux, sa densité est

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

dont la fonction de répartition est

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

Bref, on peut définir ces deux lois complètement indépendamment de leurs applications: si on considère des variables indépendantes https://perso.univ-rennes1.fr/arthur.charpentier/latex/Xi-normaux.png, alors

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

suit une loi du khi-deux, alors que

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

suit une loi du khi. Bref, intuitivement, on définit d’abord la loi du chi-deux avant d’introduire la loi du khi, si l’on part de l’interprétation. Et si l’on reprend le texte fondamental de Fisher, ici, datant de 1922, on ne retrouve que le khi-deux

Mais dans “The lady drinking tea“, David Salsburg suggère que Karl Pearson avait introduit le khi avant le khi-deux… Si l’on y regarde de plus près, Karl Pearson a proposé le “test du chi-deux” en 1900 (ici), avec l’idée que l’exposant qui apparait dans la loi normale pouvait être noté, comme on le retrouve ci-dessous

On voit aussi interprété dans ce texte le khi, et pas seulement le khi-carré.

Effectivement, si l’on regarde ici par exemple, la loi du chi-deux est reliée à l’étude du volume d’un ellipsoïde. Mais ce n’est pas la première apparition de cette variation dite du khi-deux. Si on y regarde de plus près, par exemple ici, ou , on notera que Karl Pearson avait défini, dès 1896 aussi bien la loi du khi que la loi du khi-deux, dans le cadre de ce qui s’appelle les distribution de Pearson de type III,

Autrement dit, Karl Pearson a commencé par définir les deux lois simultanément, via les fonction gamma incomplète, et autres équations d’analyse fonctionnelle, et c’est seulement plus tard que l’utilisation de la loi du khi-deux l’a emporté, de part son interprétation relativement élégante (somme de carrés de lois normales centrées réduites) et le test d’indépendance (ou d’adéquation) du khi-deux.
Mais en 1931, Karl Pearson renvoie à un papier plus ancien où l’on trouve pour la première fois la loi du khi-deux, suggérant d’ailleurs d’appeler cette loi la loi de Helmert, en mémoire à Friedrich Robert Helmert, qui avait également beaucoup travailler sur les moindres carrés (dont il est naturel d’arriver à des sommes de carrés de lois normales, ici), qui aurait trouvé cette loi en 1875. Mais on peut également remonter aux travaux d’Ernst Karl Abbe en 1863 qui avait également mentionné la loi.

Mais c’est de Karl Pearson que l’on a hérité la notation et le nom du khi-deux.

Some historical remarks on extreme values

I will start here a short post on extreme values, with some historical perspective. In a recent paper (in French), I mentioned the use of the Pareto distribution as a standard model for extremes, but if reinsurers have been using the Pareto distribution for a long time (see here e.g.), the oldest mathematical models when dealing with extreme value should be related to work on maximum values in finite samples.

  • The work of Ronald Fisher and Leonard Tippett

Leonard Henry Tippett, a former student of Karl Pearson published in Biometrika a note on extremes, in 1925. The goal was “the determination of the distribution of the range and the extremes for a large number of samples“. In 1925, everyone was looking for the Gaussian distribution everywhere, and Leonard Tippett observed that the distribution of the largest value did not have a Gaussian distribution.
A few years after, a joint work with Ronald Fisher was presented to the Cambridge Philosophical Society. The starting point was the idea of “stability” (even if the term did not appear explicitely in their work): the limiting distribution the maximum should be of the “same type” as the underlying distribution. Thus, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-01.png stands for the cumulative distribution function, it should satisfy functional equation

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

Solutions of that functional equation will give all possible limiting distributions. Thus, Fisher and Tippett obtained three possible limits,

  • solutions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-03.png, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-04.png
  • solutions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-06.png with https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-07.png (i.e. finite lower bound for the support), i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-08.png
  • solutions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-10.png if https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-11.png (i.e. finite upper bound for the support), i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-09.png

Based on those possible limiting distributions, Fisher and Tippett wanted to derive what has been called later on the “domain of attraction” of those distributions.

  • The work of Maurice Fréchet, at the same time

In 1926, Maurice Fréchet wrote a paper on “la loi de probabilité de l’écart maximum“. That paper, as well as the one by Fisher and Tippett (wrote at the same time), investigated asymptotic limits. Both obtained functional equations, but only Maurice Fréchet understood the importance of the stability concept, pointed out by Paul Levy in the context of sums. Thus, Maurice Fréchet introduced the concept of what is called now “max-stability“. But Fréchet solve only functional equation https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-05.png. The point is that Fréchet studied absolute values of errors, i.e. strictly positive random variables. Thus, Maurice Fréchet considered distribution

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

wherehttps://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-92.png is an arbitrary positive constant. The “2” comes from the fact that Fréchet considered errors with respect to the median. But he did not introduced that new distribution function, he also proved that the distribution appears as a limit when the underlying distribution of the https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-13.png‘s has an algebraic behavior at infinity, i.e. equivalent to https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-90.png, for some https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-91.png. I.e. he proved that Pareto-type tailed distibutions where in the domain of attraction of the Fréchet distribution.

  •  Later on, the work of Emil Gumbel

In 1932, Emil Gumbel gave a talk in France on the “âge limite“. But as he wrote it “on peut donc supposer que la distribution de l’âge limite – c’est à dire la probabilité que la probabilité de cet âge ait une valeur donnée – soit Gaussienne“. But a few years after, he read about Fisher’s work, and observed also that “la distribution d’une valeur extrêmes peut être représentée pour un nombre suffisant d’observations par la formule doublement exponentielle, pourvu que la distribution initiale se comporte asymptotiquement comme une exponentielle. La formule devient rigoureuse si la distribution initiale est exponentielle“, as he wrote in 1935. Thus, as Fréchet proved that Pareto type distribution were in the max-domain of attraction of Fréchet’s distribution, Gumbel obtained that exponential type distributions were in the max-domain of attraction of Gumbel’s distribution. He also introduced the term “distribution de type exponentiel
For Emil Gumbel, it was natural to study the logarithmic derivative of the distribution, since it is the mortality rate in demography (area that Emil Gumbel studied previously). As he mentioned “d’un point de vue théorique, il est intéressant de noter que M. Fréchet a construit une distribution initiale d”une variable aléatoire pour laquelle la valeur absolue de la dérivée logarithmique diminue sans limite“. But since it was not a valuable property for practical applications, he decided that “nous nous bornerons au traitement des données de type exponentiel“. Emil Gumbel always tried to relate his work on extremes and what he did on demograpy.
For instance in 1937, he wrote a paper on “les centennaires” that can also be related to the work of Bortkiewicz on rare events. He also applied his work on radioactivity, and hydrology.
In the 30’s, hydrographs as Hazen or Graszberger introduced the concept of “yearly maximum” of
a river level. They actually proposed to look for actuarial models to study decennial or centennial floods.  But they only used the lognormal distribution to model yearly maxima. In 1936, French hydrologist Aimé Coutagne met Emil Gumbel (who was teaching at the ISFA, in Lyon). At that time, Emil Gumbel was looking for possible applications (outside demography) for his doubly exponential distribution. As as pointed out by Aimé, “sa formule devait être applicable au cas des crues; c’est à dire des plus grands débits, problème analogue à celui des plus grands âges“. Not only Gumbel’s distribution gave better empirical results, but also it came with a theoritical justification.

  • Gumbel’s distribution properties

Consider the Gumbel distribution, with location and scale parameters alpha and beta respectively, i.e.

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

Note that the associated quantile function is

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

with mean

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

and variance

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-44.png
  • The work of Waloddi Weibull

Waloddi Weibull, a Swedish physict proposed a distribution in 1939, to represent the distribution of breaking strength of materials. He used it in the 50’s in reliability concept. Actually, Weibull appeared late in the story of extremes, since Fréchet, Fisher and Tippett mentioned it already in the mid-20’s.

  • From the central limit theorem (on the average) to Fisher-Tippett theorem (on the maxima)

In order to visualize those two theorem, consider the following animation, where samples of 20 exponential variables are generated. From those 20 values, we plot the maximum in blue, and the average in red, on top. Just below, be rescale those points by considering https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-16.png, and below again, https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-15.png}. When then look at the position of https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-14.png and the one of the mean of https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-15.png. We then build an histogram to visualize the distribution of the rescaled maximum (in blue) and the rescale average (in red).

For those who might be busy, after 1000 generations of samples, we obtain the following histograms (below), including the Gaussian distribution below (i.e. the average of exponential variables looks Gaussian, even with only 20 observations, actually the Gaussian distribution is only asymptotic, i.e. we should consider samples of size 2000), and the maximum over 20 observations of exponential variables (on top) looks like a Gumbel distribution (actually, here it is the exact distribution, and it is the asymptotic distribution for exponential type variables).

  • The GEV distribution

The unified expression of those three distributions is call the GEV distribution. The generalized extreme value distribution has cumulative distribution function

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

for https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-21.png, where https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-22.png is the location parameter, https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-23.png the scale parameter and https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-24.png the shape parameter. Note that the expected value is
https://perso.univ-rennes1.fr/arthur.charpentier/latex/ext-30.png