Tag Archives: test

Tests statistiques et intervalles de confiance

Comme je l’avais énoncé rapidement à la fin du chapitre sur les intervalles de confiance, il existe une forme de dualité entre tests et intervalles de confiance. Par exemple, si on considère un test sur la moyenne, de la forme http://freakonometrics.blog.free.fr/public/perso2/testic-01.gifcontre http://freakonometrics.blog.free.fr/public/perso2/testic-02.gif dans un modèle Gaussien, http://freakonometrics.blog.free.fr/public/perso2/test-ic03.gif, la région d’acceptation (de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif) sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/testic-05.gif
c’est à dire si

http://freakonometrics.blog.free.fr/public/perso2/testic-06.gif
Aussi, pour qu’une valeur http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif soit acceptée avec un niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif (erreur de première espèce associée à un test bilatéral), une condition nécessaire et suffisante est que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif appartienne à l’intervalle de confiance (symétrique) centré sur la moyenne empirique.
Il y a ainsi dualité entre deux concepts

  • prendre une valeur acceptée pour un test bilatéral de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif
  • prendre une valeur appartenant à l’intervalle de confiance symétrique de niveau 1-http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif

Considérons par exemple le test asymptotique du rapport de vraisemblance (évoqué ici), dont la région d’acceptation sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/tesp100.gif
Étant donné un échantillon http://freakonometrics.blog.free.fr/public/perso2/testh03.gif, l’intervalle de confiance pour  est

http://freakonometrics.blog.free.fr/public/perso2/tesp101.gif
Numériquement, reprenons l’exemple du jeu de pile ou face abordé ici.

> n=20; alpha=.05
> X=sample(0:1,size=n,replace=TRUE)
> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=0.5,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
>
> TLR=function(p0){2*(logL(pml)-logL(p0))<qchisq(1-alpha,df=1)}
> (IC=sapply(p,TLR))
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[73]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[97] FALSE FALSE FALSE FALSE FALSE
> c(p[min(which(IC==TRUE))],p[max(which(IC==TRUE))])
[1] 0.34 0.75


Graphiquement, on visualise l’intervalle de confiance (à 95%) entre les deux traits verticaux.
De manière duale, à partir d’un intervalle de confiance http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif pour , on peut construire un test. Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif, on accepte H_0,  sinon on rejette cette hypothèse. Par construction,

http://freakonometrics.blog.free.fr/public/perso2/testp-102.gif
Par exemple l’intervalle de confiance usuel – à 95% – pour une proportion est de la forme suivante

http://freakonometrics.blog.free.fr/public/perso2/testp107.gif
Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à cet intervalle, on accepte http://freakonometrics.blog.free.fr/public/perso2/testp109.gif (contre l’hypothèse bilatérale http://freakonometrics.blog.free.fr/public/perso2/testp101.gif).

> a=mean(X)-1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> b=mean(X)+1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> c(a,b)
[1] 0.3319638 0.7680362


La région entre les deux traits verticaux est la région d’acceptation.

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.”

Wald, score et rapport de vraisemblance

Vendredi, nous devrions aborder un peu la problématique des tests asymptotiques. Avant, rappelons que l’estimateur du maximum de vraisemblance de  vérifie la propriété asymptotique suivante.

Or on a une propriété intéressante sur les convergences de vecteurs Gaussiens (pour rester assez général). Soit  une suite de vecteurs aléatoires telle que  où . Alors

où  (pour un paramètre univarié, ça sera alors ). Mais ce résultat n’est intéressant que si  est connue. En fait, si on a la suite des variances vérifie , alors

Traduit sur notre estimateur du maximum de vraisemblance, cela signifie que

ou encore

Trois tests peut alors être mis en œuvre (et tant qu’à faire, on peut essayer de l’illustrer sur un exemple: encore et toujours du pile/face).

> set.seed(1)
> n=20
> X=sample(0:1,size=n,replace=TRUE)
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
> LL=sapply(p,logL)
> plot(p,LL,type="l",col="red",lwd=2)
> p0=.5
> points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
> abline(v=p0,lty=2)

On a alors 20 tirages de pile ou face, on a obtenu 11 piles, et on se demande si la pièce est “juste“.Mais avant toute chose, commençons par calculer l’estimateur du maximum de vraisemblance, ainsi que le score et l’information de Fisher. Dans ce modèle de variables de Bernoulli, on connaît des formes explicites de ces quantités. Mais ici, on va utiliser la méthode la plus générale qui soit, i.e. on va maximiser la log-vraisemblance, puis on va calculer le score et l’information de Fisher à la valeur de  que l’on souhaite tester.

> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=p0,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")
> Df2 = D(Df, "p")
> p=p0
> score=eval(Df)
> (IF=-eval(Df2))
[1] 80
> 1/(p0*(1-p0)/n)
[1] 80

On note que l’information de Fisher calculée par double différenciation de la log-vraisemblance est identique à la formule analytique. Pour rappels, on a

Trois tests (équivalent asymptotiquement comme on le verra en cours) sont alors possibles.

Tout d’abord le test de Wald propose de travailler sur la différence entre l’estimateur du maximum de vraisemblance, et la valeur que l’on cherche à tester (comme le montre le dessin ci-dessous). Cette différence doit être “petite” si  est la vraie valeur.On peut alors utiliser la statistique suivante

Asymptotiquement, en utilisant le résultat évoqué au début, cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> alpha=0.05
> (T=(pml-p0)^2*IF)
[1] 0.1999970
> T<qchisq(1-alpha,df=1)
[1] TRUE

Ici, on accepte l’hypothèse que la pièce n’est pas pipée.

Ensuite le test du rapport de vraisemblance propose de travailler sur les valeurs de la log-vraisemblance. Là encore, cette différence doit être “petite” si  est la vraie valeur. Graphiquement, on mesure une distance normalisée

La statistique est celle d’un test bilatéral,

On pose alors

(le 2 sera détaillé un peu en cours, c’est mis là de manière à avoir trois tests équivalents). Asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=2*(logL(pml)-logL(p0)))
[1] 0.2003347
> T<qchisq(1-alpha,df=1)
[1] TRUE

Là encore, on accepte l’hypothèse que la pièce n’est pas pipée.

Enfin, le test du score propose de travailler sur la pente de la log-vraisemblance en  . Cette pente doit être “petite” si  est la vraie valeur.

La pente correspond au score,

La statistique de test est alors ici

Et là encore, asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=slope^2/IF)
[1] 0.2
> T<qchisq(1-alpha,df=1)
[1] TRUE

Et le test accepte ici encore l’hypothèse que la pièce n’est pas pipée.

et c’est parti pour les élections…

Ce soir, @imparibus reprenait  une information du Point (ici). “Un sondage de l’institut Sensus diffusé jeudi crédite Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7% […]  la marge d’erreur du sondage est de 2,2%, a expliqué un reponsable de l’enquête [… ce qui donne une ] grande probabilité que les 2 candidats soient à égalité”. @imparibus affirmait que c’était faux,@adelaigue me demandant mon avis…
Bon, je suis flâté qu’Alexandre m’appelle à l’aide, mais je suis plutôt nul en sondage… Bref, pour faire une réponse intelligente et détaillée, il me faudra un peu plus de place que ce que ne me propose Twitter…
Soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial01.png la probabilité de voter pour R, et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial02.png la probabilité de voter pour S. On note https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial03.png la probabilité de voter pour un autre candidat, https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial04.png.
Si https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png sont les fréquences empiriques obtenues à partir de https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png personnes sondées, rappelons que

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

et

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

Aussi,

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

Si on veut tester https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial12.png (hypothèse https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png) car c’est finalement de ça qu’on parle ici, l’idée la plus naturelle serait d’introduire la statistique de test d’égalité des proportions,

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

qui suit – à la louche – une loi normale sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png.
C’est bien joli tout ça, mais je ne connais pas https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png ! En fait, on va essayer de s’aider de la petite phrase “la marge d’erreur du sondage est de 2,2%”
Alors si j’en crois le site d’ipsos (ici),  cette marge d’erreur correspond à 2000 personnes interrogées…
Or ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png =46% et et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png = 42.7%, avec

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

soit,

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

Si on compare à une loi normale, on a 2.6% de chances d’atteindre une telle valeur…
On peut aussi faire un test du chi-deux. En cas d’égalité des proportions, on devrait avoir https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png=44.75%
et si on utilise un test du chi-deux pour mesurer la distance entre les deux lois multinomiale (sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png et sous H1), on a comme distance du chi-deux

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

soit ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png=3.71 avec 2000 personnes interrogées. Sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png devrait suivre une loi du chi-deux à 1 degré de liberté, ce qui donne une p-value de 5.26%.
On peut aussi tenter un test de rapport de vraisemblance. Sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png, on a

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

alors que sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial23.png

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

(en utilisant la forme des lois multinomiales). Le rapport des vraisemblances s’écrit alors

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

soit ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial27.png= 10.436. La statistique de test est icihttps://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial27.png qui doit suivre, sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png une loi du chi-deux à un degré de liberté. La p-value associée est alors de l’ordre de 0,1%….
Bon, il semble qu’on puisse raffiner un peu les tests, en reprenant les idées de Nass (1959) et de Hoel (1938) – mais qui sont intéressant si on a peu de données (ce qui n’est pas vraiment le cas ici).
Bref, j’aurais tendance à être plutôt d’accord avec @imparibus. Maintenant, je suis nul en sondage, donc si quelqu’un a d’autres idées, je suis preneur ! Surtout que j’ai peut-être fait quelques erreurs de calculs ici ou là….

Sexe et naissances

Depuis quelques jours, nous faisons la présentation de la petite dernière à la famille et aux amis, et nous étions surpris du nombres de familles n’ayant que des filles (nous ne nous étions pas fait la réflexion jusqu’alors). La première piste était de remettre en cause l’indépendance des sexes entre enfants d’une même famille. La second théorie était de dire que les familles avec déjà deux petits gars devaient capituler rapidement, contrairement aux familles qui avaient des filles, qui pouvaient continuer… En gros, “les filles c’est plus simple à élever que les p’tits mecs”. Forcément, j’ai voulu mener l’enquête….

A partir du fichier enfant de l’enquête Education et Famille, EPCV d’octobre 2003, on peut récupérer quelques statistiques. Dans cette base, on dispose de 10027 enfants, issus de 4081 familles. Je ne reviens pas sur le sens de la famille disons que je vais travailler par identifiant du logement, i.e. les enfants habitant sous un même toit. Les personnes qui veulent en savoir plus iront sur le site de l’insee (ici).

  • De l’indépendance des sexes

Dans cette base, j’ai 50,52% de garçons et 49,48% de filles. 19% des enfants semblent être des enfants uniques, et dans 51,41% des cas il s’agit d’un garçon. Pour faire une étude sur la distribution des sexes, on note que les familles avec au moins deux enfants semblent représentatives de l’ensemble. Considérons alors les familles ayant eu au moins deux enfants, soit 3299 familles. La répartition se fait de la manière suivante, entre le sexe du premier enfant, et le sexe du second,

sexe Garçon (2) Fille (2)
Garçon (1) 857 801
Fille (1) 813 828

Si on fait un test d’indépendance du chi-deux, on observe qu’il y a indépendance entre le sexe des deux enfants,

>  chisq.test(tableau)
 Pearson's Chi-squared test with Yates' continuity correction
X-squared = 1.4346, df = 1, p-value = 0.231

Bref, on admettra qu’il y a indépendance entre les sexes au sein d’une famille….

  • Des familles dont les deux premiers enfants sont du même sexe

On notera (sur le tableau ci-dessus) qu’il y a plus de familles à deux garçons que de familles à deux filles. Parmi les 857 familles qui ont commencé par 2 garçons, 379 familles ont eu un troisième enfant, soit 44,22%, et 48,81% ont eu encore un garçon. Parmi les 828 familles qui ont commencé par 2 filles, 388 familles ont eu un troisième enfant, soit 46,86%, et d’ailleurs 50,77% ont encore eu une fille (cette statistique confirmant notre hypothèse d’indépendance).
Cette différence (entre les proportions de familles qui ont tenté un troisième enfant après les deux premiers du même sexe) ne semble pas négligeable, mais si on fait un test de Student, la différence n’est toutefois pas significative,

        Welch Two Sample t-test
t = -1.0858, df = 1680.417, p-value = 0.2777
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval: -0.07397060  0.02125328 
sample estimates: mean of x mean of y0.4422404 0.4685990

Si on compare aux familles qui ont eu un garçon et une fille, on notera que 42,68% de ces familles tentent le troisième enfant (contre 45.51% pour les familles qui ont eu deux premiers enfants du même sexe). Cet fois, la différence semble presque significative,

        Welch Two Sample t-test 
t = -1.637, df = 3292.66, p-value = 0.1017 
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval: -0.062203214  0.005596887 
sample estimates: mean of x mean of y 0.4268897 0.4551929

Bref, l’hypothèse comme quoi les familles avec uniquement des garçons sont moins nombreuses que les familles avec uniquement des filles ne semble pas tenir la route statistiquement….en revanche, il semble que les familles veulent à la fois avoir un garçon et une fille. Bon, c’est pas tout, mais j’ai une couche à changer… la prochaine fois je parlerais de sexe et de vaisselle…

PS: ce souhait d’aborder des sujets sur le sexe n’ont rien à voir avec le fait que, depuis que j’ai été référencé sur sexactu.com, la fréquentation sur mon blog a explosé… enfin, ça n’a duré qu’une journée, tout le monde s’est rendu compte de la supercherie…

Peut on faire l’économie du formalisme quand on parle d’extrêmes ?

Tous les blogs économiques saluent la parution en poche du joli petit livre de Daniel Zajdenweber, Economie des Extrêmes. En particulier, beaucoup de monde salue ce livre qui explique simplement des choses complexes…. Par exemple Alexandre dès 2001 “Passé le premier chapitre, un peu ardu, et qui nécessite du lecteur des connaissances de base en statistique et probabilités (notion de lois de probabilité, d’espérance, de variance…) qui décrit en termes littéraires les caractéristiques de ces lois, l’auteur applique ces résultats à un grand nombre de phénomènes concrets, et en tire les conséquences“. Mais peut-on parler d’économie des extrêmes sans être technique ?

Histoire que mon message ne soit pas déformé, je trouve passionnant ce petit livre introductif à la problématique des risques extrêmes (qui est un de mes dadas depuis quelques années) mais j’espère qu’il servira d’encouragement à une lecture d’ouvrages plus détaillés sur le sujet. Car la vulgarisation a des limites que l’on atteint vite quand on parle de sujets aussi complexes.

L’exemple que j’ai le plus étudié est celui des sinistres de perte d’exploitation (longuement évoqué par Daniel Zajdenweber dans son livre). Il y a quelques années j’avais utilisé cette partie du livre comme base pour faire un sujet d’examen pour le cours de “réassurance et grands risques” que je donnais alors à l’ENSAE1. Et malheureusement, mes compétences littéraires sont très limitées, donc je vais faire des maths. Dans le livre, la figure suivante est présentée,

qui correspond effectivement à la fonction tracée dès 1925 par Karl Gustav Hagstroem (j’avais souligné (ici) ses travaux précurseurs où l’intérêt de la loi de Pareto pour modéliser les très grands riques apparaissait pour la première fois). C’est en effet assez naturel: si on a une loi de Pareto, i.e.

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

alors on pourrait écrire, en passant au logarithme

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

Si on représente la version empirique, c’est à dire le nuage de points

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

alors pour une loi de Pareto, les points devraient être alignés suivant une droite, et la pente doit correspondre au paramètre de la fonction puissance. C’est visiblement l’idée exploitée ici.

Autrement dit, les pointillés ne sont un intervalle de confiance, mais juste un outils graphique pour se demander si la pente vaut 1, ou pas. Daniel Zajdenweber affirme que la pente doit ici être -1.

Le fait que la valeur soit unitaire ou pas a en effet un impact très important en terme d’assurabilité du risque de perte d’exploitation. Rappelons que pour une variable positive (et c’est le cas ici). Et si on a une telle loi de Pareto (de puissance unitaire), alors la prime pure d’un traité de réassurance, couvrant entre m et M s’écrit

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

soit

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

ce qui correspond aux calculs de Daniel Zajdenweber… Mais encore une fois “l’absence d’espérance mathématique de la distribution des sinistres” est une conclusion très forte sur laquelle on peut essayer de revenir.

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

aussi ici

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

autrement dit, l’espérance est finie si la pente est strictement plus grande (en valeur absolue) que 1. Si la pente est inférieure (ou égale) à 1, le risque n’est pas assurable ! Ce qui est une conclusion très très forte pour les assureurs.
J’ai donc demandé à la FFSA la base de données utilisée ici, et pour éviter des problèmes d’inflation des coûts de sinistres entre 1992 et 2000. Si je prends tous les sinsitres, on obtient l’ajustement de Pareto suivant

soit une pente (en valeur absolue de 1.47). Mais encore une fois, l’ajustement de Pareto se fait sur les grands sinistres. Hill a proposé un estimateur très populaire pour estimer ce coefficient, où on ne prend en compte que les k observations les plus grandes, et on regarde l’estimation de la pente du graphique de Pareto pour ces quelques valeurs. On représente alors l’estimation en fonction du nombre de grands sinistres, ou du seuil définissant les graphs sinistres. Numériquement, en posant

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

on peut écrire comme estimateur de la pente

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

soit, en simplifiant le numérateur,

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

tel que l’a construit Hill en 1975. Graphiquement, on a ici

Bref, la question est ce savoir si on atteint la valeur 1 pour les grands sinistres. Graphiquement, on a malgré tout envie de rejeter cette hypothèse.
Une solution peut être de faire un test statistique, basé sur de ratio de vraisemblance, comme le suggèrent Reiss & Thomas (2001) or Coles (2001). En fait, on peut même utiliser d’autres estimateurs que celui de Hill, comme celui obtenu en faisant un ajustement de loi GPD (Pareto généralisée
) sur la loi des Excès, ou une loi GEV sur des maximas par blocs (Generalized Extreme Value). On introduit alors la statistique de test suivante

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

et on regarde les p-value (ainsi que la correction de Bartlett à droite),

On peut aussi, plus simplement, estimer plusieurs coefficients de pentes pour des seuils différents, et regarder la borne supérieure de l’intervalle de confiance,

Bref, même si avec un des ajustements de loi GPD on hésite à retenir une pente unitaire, la plupart des tests rejettent cette hypothèse, et donc le risque de perte d’exploitation semble être assurable, d’espérance mathématique finie. Bref, les dessins c’est très bien pour faire passer une idée, mais ne retenir que ça pour en tirer des conclusions aussi fortes me laisser sceptique….

Choix multiples et théorie des tests

En poursuivant la première question du QCM de “probabilité & statistiques” (en ligne ici pour l’énoncé, et pour une ébauche de correction), on a l’occasion de voir une application intéressante (ou tout du moins selon mon point de vue) de la théorie des tests.

Je souhaite tester une hypothèse

H0:  la personne a répondu au hasard aux questions

(je spécifierais plus tard l’hypothèse alternative, de toutes les façons, on n’en a pas besoin pour construire la statistique de test.)
En répondant au hasard (i.e. sous l’hypothèse H0), la personne avait 1 chance sur 3 de gagner 1 point, et 2 chances sur 3 de perdre 0,5 point1. Bref, à une transformation affine prêt, pour chaque question, on a une loi de Bernoulli. Plus précisément le nombre de points gagnés (ou perdus) à la question i est

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-1.png

Aussi, la loi du nombre total de points est

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-2.png

or si l’on suppose que l’on réponde au hasard à toutes les questions, on peut suppose les variables indépendantes, et donc

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-3.png

et sans être trop prétentieux, on connaît de cette loi (en particulier sa fonction de répartition et sa fonction quantile, au moins numériquement en tous les cas). Bref, si F est cette fonction de répartition,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-6.png

ensuite on sait faire: c’est juste trouver la loi d’une transformation croissante de cette variable (linéaire en fait). Aussi

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-7.png

c’est à dire

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-9.png

Ben ça peut se représenter facilement avec R
> X=seq(-10,20,by=.1)
> G = function(x){pbinom((10+x)*2/3,size=20,prob=1/3)}
> Y=G(X)
> plot(X,Y,type=”l”)

Notons que je peux aussi avoir la fonction quantile de cette loi. On cherche q tel que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-10.png

C’est à dire tel que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-11.png

i.e.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-12.png

En  R, on obtient la fonction suivante
> G.inv = function(u){1.5*qbinom(u,size=20,prob=1/3)-10}
Voilà, j’ai donc la loi de la note (qui me servira de statistique de test) sous l’hypothèse H0. Bon, ça suffit presque à construire complètement le test. Il suffit de spécifier l’hypothèse alternative afin de construire les régions critiques. A priori je vois deux formes pour l’hypothèse alternative,

H1:  la personne n’a pas répondu au hasard aux questions

ce qui correspondra à un test bilatéral, incluant ceux qui ont des bonnes notes (et qui sont plus fort que le hasard), et ceux qui ont des très mauvaises notes (et qui sont donc moins fort que le hasard en quelque sorte). L’autre hypothèse alternative, qui a plus de sens à mes yeux serait

H1:  la personne est plus forte que le hasard

ce qui sous-entend qu’elle a un peu travaillé (ou qu’elle a raisonné un minimum).
Dans le premier cas, on cherche les deux seuils de la région de rejet
> (note.inf = G.inv(0.05))
-5.5
> (note.sup=G.inv(.975))
6.5 

ce qui donnera un intervalle d’acceptation à 95%. Dans le second cas, on prendra une région de rejet qui commencera à la note minimal (-10), et on montera jusqu’au quantile à 95% de G, soit
> (note.sup=G.inv(.95))

Moralité, il semblerait que pour 20% des élèves, on accepte l’hypothèse comme quoi ils ont répondu au hasard, et on atteint 30% si on prend l’autre hypothèse alternative. En fait, étrangement, si on regarde les réponses à la première question sur le calcul de la moyenne, presque qu’un tiers des élèves a répondu que s’ils cochaient au hasard, ils auraient 6,6667/20. Or ce tiers a précisément obtenu 6,83/20. Donc il semble qu’un tiers des personnes aient effectivement répondu hasard !
Si j’autorise la personne à choisir (au hasard) l’option de ne pas répondre, je me retrouve avec une loi dite multinomiale. On prendra par la suite la proba de ne pas répondre à la question comme paramètre, afin de voir comme le test évolue avec ce paramètre.
Pour décrire la loi mutinomiale

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-13.png

Bref, là aussi on peut obtenir la fonction de répartition de cette loi, ainsi que sa densité… au moins numériquement, car la somme de lois multinomiales indé
pendantes, ça peut vite devenir un peu compliqué.. Bon, je vais faire du numérique par la suite…
On va commencer par supposer que dans 10% des cas, on ne sait pas répondre, ou plutôt on ne répond rien… Dans ce cas, on note que cela ne change pas trop (au moins sur la région de rejet du test).

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

Si par contre on suppose que dans 1/4 des cas, on ne répond pas, la région rétrécit,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

et encore davantage si on suppose qu’on ne répond pas dans 50% des cas…

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

1 je ne prends pas en compte la possibilité qu’une personne puisse ne pas répondre…. on verra ça plus tard.