Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Margin of error, and comparing proportions in the same sample

Irecently tried to answer a simple question, asked by @adelaigue. Actually, I thought that the answer would be obvious… but it is a little bit more compexe than what I thought. In a recent survey about elections in Brazil, it was mentionned in a French newspapper that “Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7%” (i.e. proportions obtained from the survey). It is also mentioned that “la marge d’erreur du sondage est de 2,2% ” i.e. the margin of error is 2.2%, which means (for the journalist) that there is a “grande probabilité que les 2 candidats soient à égalité” (there is a “large probability” to have equal proportions).
Usually, in sampling theory, we look at the margin of error of a single proportion. The idea is that the variance of https://latex.codecogs.com/gif.latex?%20\widehat{p}, obtained from a sample of size https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png is

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

thus, the standard error is

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

The standard 95% confidence interval, derived from a Gaussian approximation of the binomial distribution is

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

The largest value is obtained when p is 1/2, and then we have a worst case confidence interval (an upper bound) which is

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

So with a margin of error https://perso.univ-rennes1.fr/arthur.charpentier/latex/m205.png means that https://perso.univ-rennes1.fr/arthur.charpentier/latex/m206.png. Hence, with a 5% margin of error, it means that n=400. While 2.2% means that n=2000:
> 1/.022^2
[1] 2066.116
Classically, we compare proportions between two samples: surveys at two different dates, surveys in different regions, surveys paid by two different newpapers, etc. But here, we wish to compare proportions within the same sample. This has been consider in an “old” paper published in 1993 in the American Statistician,

It contains nice figures to illustrate the difference between the standard approach,

and the one we would like to study here.

This point is mentioned in the book by Kish, survey sampling (thanks Benoit for the reference),


Let https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png denote empirical frequencies we have obtained from the sample, based on https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png observations. Then since

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

and

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

we have

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

Thus, a natural margin of error on the difference between the two proportion is here

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

which is here 4 points
> n=2000
> p1=46.8/100
> p2=42.7/100
> 1.96*sqrt((p1+p2)-(p1-p2)^2)/sqrt(n)
[1] 0.04142327
Which is exactly the difference we have here ! Hence, the probability of reaching such a value is quite small (2%)
> s=sqrt(p1*(1-p1)/n+p2*(1-p2)/n+2*p1*p2/n)
> (p1-p2)/s
[1] 1.939972
> 1-pnorm(p1-p2,mean=0,sd=sqrt((p1+p2)-(p1-p2)^2)/sqrt(n))
[1] 0.02619152

Actually, we can compare the three margin of errors we have so far,

  • the upper bound
https://perso.univ-rennes1.fr/arthur.charpentier/latex/m208.png
  • the “average one”
https://perso.univ-rennes1.fr/arthur.charpentier/latex/m209.png

where

https://perso.univ-rennes1.fr/arthur.charpentier/latex/m212.png
  • the more accurate one we just obtained,
https://perso.univ-rennes1.fr/arthur.charpentier/latex/m213.png

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/m214.png.
> p=seq(0,.5,by=.01)
> ic1=rep(1.96/sqrt(4*n),length(p))
> ic2=1.96*sqrt(p*(1-p))/sqrt(n)
> delta=.01
> ic31=1.96*sqrt(2*p-delta^2)/sqrt(n)
> delta=.2
> ic32=1.96*sqrt(2*p-delta^2)/sqrt(n)
> plot(p,ic32,type=”l”,col=”blue”)
> lines(p,ic31,col=”red”)
> lines(p,ic2)
> lines(p,ic1,lty=2)
So on the graph below, the dotted line is the standard upper bound, the plain line in black being a more accurate one when the probability is https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png (the x-axis). The red line is the true margin of error with a large difference between candidates (20 points) and the blue line with a small difference (1 point).


Remark: an alternative is to consider a chi-square test, comparering two multinomial distributions, with probabilities https://perso.univ-rennes1.fr/arthur.charpentier/latex/m215.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/m216.png where https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png is the average proportion, i.e. 44.75%. Then

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

i.e.  https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png=3.71
> p=(p1+p2)/2
> (x2=n*((p1-p)^2/p+(p2-p)^2/p))
[1] 3.756425
> 1-pchisq(x2,df=1)
[1] 0.05260495
Under the null hypothesis, https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png should have a chi-square distribution, with one degree of freedom (since the average is fixed here). Here the probability to reach that level is around 5% (which can be compared with the 2% we add before).

So finally, I would think that here, stating that there is a “large probability” is not correct…

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à….

Copules et risques corrélés

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.copula-density-proj_m.jpg

J’avais promis dans un commentaire que je mettrais bientôt en ligne un survey sur les copules…. après plusieurs mois de retard, le document est en ligne [pdf], et il est sobrement intitulé copules et risques multiples. Toutes les remarques et critiques sont les bienvenues ! Il s’agit d’un chapitre pour un livre dont je mettrais la référence en ligne ultérieurement. A priori, ce document devrait servir de base pour le cours qui sera donné dans un mois au CIRM (mentionné ici, dans le cadre des Journées d’Études Statistiques).

Quelques heures d’effervescence sur Twitter

Pour ceux qui dormaient entre 11 heures et 15 heures (à Montréal, soit 17 heures – 21 heures en France), hier, Bernard Clavel est mort. Auparavant, j’allais sur le site de Libé, ou du Monde, pour me tenir informé de ce qui se passe en dehors de mon bureau (ou de ma cuisine),

Maintenant je suis sur Twitter… et autant l’avouer, j’ai essayé des outils comme TweetDeck mais au bout de 2 heures, force est de constater que j’ai arrêté l’application car j’en avais marre d’être harcelé par ces fenêtres…
Sur Twitter, tout le monde a quelque chose à dire sur tout, et ne peut s’empêcher de vouloir le dire à tout le monde…. Par exemple “Bernard Clavel” ne faisait pas franchement parti des mots tapé régulièrement sur Twitter…. Mais à l’annonce de sa mort, tout le monde veut diffuser l’information à son réseau…. Ça allait tellement vite que j’étais obligé de faire des backups des pages internets tous les quarts d’heure (sinon l’historique disparaît, ce qui limite l’intérêt des outils statistiques). Pour être précis, ma recherchait était

http://twitter.com/#search?q=bernard%20clavel

et pendant près de 4 heures, ça allait du tweet original, au retweet, au message recopié sans citation, etc…. J’ai récupéré ainsi plus de 400 tweets, mais sans l’horaire…. Finalement, j’ai découvert

http://search.tweetreports.com/search_results.php?search=bernard+clavel&sear=Y

Pendant l’heure qui suit les premiers tweets, c’est l’effervescence, puis ça retombe… avec un pic à 20 heures, visiblement certains l’apprenant aux informations du soir…

Manifestement la loi gamma (la date initiale étant le premier tweet annonçant le décès) n’est pas adaptée…

Moralité, pour poursuivre le billet de la semaine dernière (ici), on peut tenter un mélange, en l’occurrence un mélange gamma/exponentielle…. la première loi gamma permet d’intégrer le délai de transmission de l’information au réseau, alors que pour le second, il semble que ce soit une annonce au 20 heures qui ait déclenché la soudaine envie de tweeter..

Allez, le jour où je maîtrise mieux les outils linguistiques, je regarde en plus quand les tweets anglophone et hispaniques sont arrivés…..

Journées d’Etudes Statistiques, 2010

Cette année, les Journées d’Études Statistiques seront organisées du lundi 15 novembre au vendredi 19 (ici pour les inscriptions) au CIRM, à Marseille, et portera cette année sur les approches statistiques du risque.

Les journées sont organisées par Jean Jacques Droesbeke, Gilbert Saporta et Christine Thomas-Agnan. Et parmi les intervenants, on retrouvera Patrice Bertail, Anne Laure Fougères, Stéphane Loisel, Jean Marc Tallon, ou encore Eric Parent. Les journées de 2008 portaient sur les variables latentes et les mélanges (ici) et celles de 2006 sur les méthodes nonparamétriques en régression. Les contributions sont d’ailleurs régulièrement publiées. A suivre donc…

Too large datasets for regression ? What about subsampling….

recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results…

  • Time to run a regression, as a function of the number of observations

Here, I generate a dataset as follows

http://freakonometrics.hypotheses.org/files/2016/11/largesamp01.png

and we fit

http://freakonometrics.hypotheses.org/files/2016/11/largesamp02.png

where http://freakonometrics.hypotheses.org/files/2016/11/largesamp03.png is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply

> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
+1*(X2=="A")-2*(X2=="B")-5*(X2=="C")
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)

We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png) of the dataset.

> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25

Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,

> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86

Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it will be more than two timess longer,

Convexity can simply be interpreted as “too large datasets take time, and too small too…”. Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets….

  • Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?

Here, we have datasets with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines. The questions is how long will it take if we subdived into http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png subsamples (of equal size), and run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regressions ?

> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))

[…]

  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31

On the graph below, we have the time (y-axis, here on a log scale) it took to run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, as function of http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png (x-axis), including the time it took to run the regression on a dataset of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png which is the concentration of dots on the left (i.e. http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=1), both on the 6 regressors – in black – and with a strepwise procedure – in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, when http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png is not too large, i.e. less than 10 or 15, it is not longer than the regression on http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines… But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?

  • What about the quality of the output ?

Here, we consider only one dataset, with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=100,000 lines (just to make it run a bit faster). And http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=20 subsets. Recall that the generated dataset is from

http://freakonometrics.hypotheses.org/files/2016/11/largesamp01.png

and we fit

http://freakonometrics.hypotheses.org/files/2016/11/largesamp02.png

Here, we plot here http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png and a confidence interval, defined as

http://freakonometrics.hypotheses.org/files/2016/11/largesamp08.png

The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png segments on the right are the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators (each on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png).

We can see that we have much more volatility on those http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators, but the average (horizontal doted lines) are not so bad… The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the http://freakonometrics.hypotheses.org/files/2016/11/largesamp09.png the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we’ve seen earlier) than running only one on the whole dataset…. but it provides better estimates. So the tradeoff is not that simple, and maybe running several regressions on huge datasets can be a proper alternative.

Detecting distributions with infinite mean

In a post I published a few month ago (in French, here, based on some old paper,there), I mentioned a statistical procedure to test if the underlying distribution https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO02.pngof an i.i.d. sample https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO01.png had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO03.png=200 observations generated from a Pareto distribution

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

(here https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, as a start)

> a=2
> X=runif(200)^(-1/a)

Here, we will use the package developped by Mathieu Ribatet,

> library(RFA)
  • Using Generalized Pareto Distribution (and LR test)

A first idea is to fit a GPD distribution on observations that exceed a threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO11.png>1.
Since we would like to test https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png (against the assumption that the expected value is finite, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO07.png), it is natural to consider the likelihood ratio

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

Under the null hypothesis, the distribution of https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO09.png should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let  us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill’s estimator) to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO10.png, and thus, to fit a GPD on the https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.pnglargest values. And then to plot all that on a graph (like the Hill plot)

> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }

Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model (https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png the likelihood ratio statistics has a chi-square distribution).
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png largest values (and a confidence interval for the estimator),

If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance), we have

i.e. for those two models, we clearly reject the assumption of infinite mean (even if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png gets too close from 1, we should consider thresholds large enough). On the other hand, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),

  • Using Hill’s estimator

An alternative could be to use Hill’s estimator (with Alexander McNeil’s package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill’s estimator,

> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)

Again, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we clearly rejct the assumption of infinite mean,

and similarly, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance)

Here the test is more robust than the one based on the GPD. And if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean), again we accept https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png,

Note that if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=0.7, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,

(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

Tiens voilà, la pluie. Ah! quel sale temps. Où est-il l’été ? l’été où est-il ?

Dans Les cafards, de Jo Nesbø, Harry Hole (qui vient d’Oslo pour ceux qui n’ont pas dévoré la série) s’entend dire à un moment qu’à Bangkok, les discussions quotidiennes portent rarement sur “la pluie et le beau temps” pour la bonne raison que le climat est assez prévisible à Bangkok (par contre on parle quotidiennement du trafic routier). J’avais évoqué ici l’idée reçue que nous avions, sur le fait que le temps à Rennes n’est pas si changeant que ça. Et je me suis rendu compte depuis que nous sommes à Montréal que le temps pouvait vraiment changer d’un jour sur l’autre (et encore, paraît-il, je n’ai rien vu).
Alors tout d’abord, pour ceux qui en douteraient, la température à Rennes, et à Montréal, ça n’est pas vraiment la même chose,

pour Montréal, alors que pour Rennes on obtient (je précise mais je pense qu’on aurait deviné sans)

Autrement dit, on a une dispersion presque deux fois plus grande à Montréal. Si regarde également les variations (min/max) dans la journée, on obtient, pour Montréal les courbes suivantes,

avec enrouge, la valeur moyenne (sur 15 ans) du maximum observé dans la journée, et enbleu, la valeur moyenne du minimum. A Rennes, on retrouve le fait que la dispersion est relativement faible,

En fait, les choses sont encore pire quand on regarde les quantiles, i.e. le pire minimum observé dans 10% des cas (les plus froids) et le pire maximum observé dans 10% des cas (les plus chauds);

à  Montréal, alors qu’à Rennes, on est beaucoup plus resserrés,

Mais ces pires de cas ne sont pas obtenus dans la même journée, on ne peut pas en dire grand chose. Si on veut aller un peu plus loin, on peut regarder la dispersion dans la journée (i.e. différence entre le minimum et le maximum). On note qu’il est assez stable à Montréal (toujours d’environ 10 degrés)

alors qu’à Rennes, la différence est de 5 degrés l’hiver, contre 10 degré l’été,

Mais surtout, si on regarde la variation d’un jour sur l’autre, danns 90% des cas, on varie de +/- 8 degrés, sur la température moyenne journalière,

avec pas mal de cas un peu extrêmes, avec des variation de 10 ou 15 degré d’un jour sur l’autre sur la température moyenne, alors qu’à Rennes, on a plutôt +/- 5 degrés (avec très peu d’évènements extrêmes, car j’ai tronqué les ordonnées afin de mieux voir au centre)

Bref, je ne suis pas prêt de m’arrêter de parler longtemps de la pluie et du beau temps !

EM and mixture estimation

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

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

In order to describe how EM algorithm works, assume first that both https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.pngare perfectly known, and the mixture parameter is the only one we care about.

  • The simple model, with only one parameter that is unknown

Here, the likelihood is

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

so that we write the log likelihood as

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

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (that cannot be observed), taking value when https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM07.png is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and 0 if it is drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png. More generally (especially in the case we want to extend our model to 3, 4, … mixtures), https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png.
With that notation, the likelihood becomes

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

and the log likelihood

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

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png is known, denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png. Then I can predict the value of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png) for all observations,

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

So I can inject those values into my log likelihood, i.e. in

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

having maximum (no need to run numerical tools here)

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

that will be denoted https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM17.png. And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngis the best predictor of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM19.png given my observations (as well as my belief in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png). Then comes a maximization (M) step, where using https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png, I can estimate probability https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png.

  • A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM30.png. Recall that we had https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png.
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png, instead of being in the segment https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM31.png was in https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM32.png, then we could have considered mean and standard deviations of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=0, and similarly on the subset of observations such that https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=1.
But we can’t. So what can be done is to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png as the weight we should give to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and similarly, 1-https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.pngwould be weights given to observation i when estimating parameters of https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png.
So we set, as before

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

and then

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

and for the variance, well, it is a weighted mean again,

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

and this is it.

  • Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png and then https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }
Then I get https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png, or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png. So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?… next time, maybe….

Optimization and mixture estimation

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised…

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-01.png proportional to

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

The true model is https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-03.png, and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+  sd2 <- param[5]
+  -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 – p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+  }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+   fn = logvrainegX)
Actually, since I might have identification problems, I take either https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-06.png, depending whether https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-07.png or https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-08.png is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: theblack vertical segments. Obviously, when the sample is not enough heterogeneous (i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-09.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical  starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again… so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad…. but the probability to be far away from the tru value is not small at all… except when the difference between the two means exceeds 3…
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad… and the two algorithm perform similarly (maybe the first one being a bit closer to the trueparameter).

La tarification avec SAS

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble  que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

DATA contrats;
SET lib.contrats;
lnexpo = log(expo);
run;

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1_Param
           Type3=Genmod1_Var
           Modelfit=Genmod1_InfoModele; 
MODEL nbsin = ageconducteur /
                  dist = poisson   
                  link = log   
                  offset = lnexpo 
                  type3;
RUN; QUIT;

La sortie SAS a alors l’allure suivante

                                  The GENMOD Procedure
                    Critère pour évaluer la qualité de l'ajustement
              Critère                   DF          Valeur       Valeur/DF
              Deviance                63E3      26872.5334          0.4237
              Scaled Deviance         63E3      26872.5334          0.4237
              Pearson Chi-Square      63E3      73275.5362          1.1553
              Scaled Pearson X2       63E3      73275.5362          1.1553
              Log Likelihood                   -18474.2667

       Algorithm converged.
                       Analyse des résultats estimés de paramètres

                              Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

                         Statistiques LR pour Analyse de Type 3
                      Source           DF      Khi 2    Pr > Khi 2
                      ageconducteur     1     148.72        <.0001

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

PROC GAM DATA = base;
MODEL nbsin = spline(ageconducteur) / dist = Poisson;
OUTPUT OUT=gam PREDICTED; 
RUN; QUIT;
PROC SORT DATA = gam NODUPKEY; BY age_cond; RUN; QUIT;

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;
RUN;

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
PROC MEANS DATA = base NOPRINT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo
y_uni_nbsin;
RUN; QUIT;

ensuite, on fait un modèle GLM, et  un modèle GAM, et on récupère les sorties

PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
      gam;
BY age_cond;
RUN;

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1bis_Param
           Type3=Genmod1bis_Var
           Modelfit=Genmod1bis_InfoModele; 
MODEL nbsin = ageconducteur /
                 dist = poisson  
                 link = log      
                 offset = lnexpo 
                 type3           
                 scale = deviance;
RUN; QUIT;

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

                      Analyse des résultats estimés de paramètres

                                   Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

%MACRO CALCUL_AIC_BIC(infomodel=, param=);
    DATA _null_;
    SET &infomodel.;
    IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
    IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
    RUN;
    DATA _null_;
    SET &param.  end=fin;
    RETAIN nb_df 0;
    nb_df = nb_df + df;
    IF fin THEN CALL SYMPUT("k", nb_df);
    RUN;
    DATA genmod_aic_bic;
    SET &param.;
    FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
    Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
    IF Parameter = "Intercept";
    KEEP Loglike K N AIC_CALC BIC_CALC;
    RUN;
    DATA genmod_aic_bic;
    SET genmod_aic_bic;
    Loglike = &loglike.;
    K = &k.;
    N = %eval(&n_etoile. + &k.);
    AIC_CALC = 2 * Loglike + 2 * K;
    BIC_CALC = 2 * Loglike + K * log(N);
    RUN;
    PROC PRINT DATA = genmod_aic_bic;
    RUN; QUIT;
%MEND CALCUL_AIC_BIC;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);

Qui peut m’aider à comprendre les sorties de SAS ?

Je m’étais promis que j’évoquerais une bizarrerie rencontrée avec SAS lors d’une formation…. Écrire ce billet permettra à ceux qui auraient des éléments d’explication de poster un commentaire.
Pour cela, comparons une régression logistique faite avec deux outils différents, sous SAS,

  • avec la procédure logistique

Le code pour faire une régression logistique ressemble à ça

PROC LOGISTIC DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / selection=stepwise;
RUN; QUIT;

ce qui donne la sortie suivante (je passe l’introduction pour insister sur les coefficients)

                                 The LOGISTIC Procedure

                   Analyse des estimations de la vraisemblance maximum
                                                     Erreur         Khi 2
  Paramètre                    DF    Estimation         std       de Wald    Pr > Khi 2

  Intercept                     1        1.7833      0.0676      696.9022        <.0001
  sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
  age_soc      1_AGESOC_-60     1        0.4578      0.0667       47.1020        <.0001
  fract_paiemt Annuel           1        0.6021      0.0997       36.4862        <.0001
  fract_paiemt Mensuel          1       -0.5410      0.0842       41.2342        <.0001
  • avec la procédure genmod (car la régression logistique est un glm)

On peut faire exactement la même chose (théoriquement) en ajustement un modèle GLM,

PROC GENMOD DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / dist = binomial;
RUN;

et la sortie ressemble à ça

                                  The GENMOD Procedure
                      Analyse des résultats estimés de paramètres

                                                  Erreur      Wald 95Limites
Paramètre                     DF   Estimation   standard      de confiance %       Khi 2
Intercept                      1       1.5073     0.1501     1.2131     1.8014    100.85
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000       .
age_soc        1_AGESOC_-60    1       0.9156     0.1334     0.6542     1.1771     47.10
age_soc        Z_AGESOC_+60    0       0.0000     0.0000     0.0000     0.0000       .
fract_paiemt   Annuel          1       0.6634     0.1770     0.3165     1.0104     14.05
fract_paiemt   Mensuel         1      -0.4798     0.1510    -0.7759    -0.1838     10.09
fract_paiemt   Semestriel      0       0.0000     0.0000     0.0000     0.0000       .
Scale                          0       1.0000     0.0000     1.0000     1.0000
  • comparaison des deux sorties

Si on regarde l’impact du sexe par exemple, dans la première sortie on peut lire

sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
alors que dans la seconde sortie, on a
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000

On dira ce qu’on veut, mais moi je trouve cette différence troublante…. Dans la seconde sortie, le coefficient vaut le double de l’autre….
Alors SAS semble s’y retrouver car si on lui demande d’afficher le score prédit pour un individu au hasard (le premier de la base par exemple), les prédictions sont très proches,

                           fract_                                 proba1_       proba1_
  Obs  sexe_soc   age_soc  paiemt      SPOCAM  proba1_logit        
1      Homme          71  Annuel         0      0.10242637    0.10241302

Si quelqu’un sait interpréter ce qui est fait avec cette procédure logistique (car R donne la même chose que la sortie GLM), je suis preneur…..


“Le prix de l’or bat de nouveaux records”… et alors ?

Histoire de changer un peu, un rapide billet d’économie… bien que l’économie monétaire ne soit pas mon domaine de prédilection, loin de là. Mes compétences en la matière (comme sur beaucoup de notions en économie) se limitent à la lecture (mais aussi une intense relecture) d’Obélix et Compagnie

je dois avouer que je suis surpris par les commentaires de certains journalistes économiques, qui nous expliquent que “l’or a atteint des niveaux historiques” (ici, même si ça date un peu) ou que “le prix de l’or bat de nouveaux records” (ce matin cette fois, ). On peut aussi lire ici ou . On peut s’en convaincre en jetant un œil aux graphiques ici, , ou encore. Sans vouloir être méchant, si je transpose ca a une discussion que je pourrais avoir avec mes enfants, ça donnerait: “On a mis 150 euros sur un livret A le jour de votre naissance… c’est incroyable, mais tous les 15 jours, on bat des records: vous êtes de plus en plus riche !”

En fait, effectivement, le prix de l’or atteint des records… mais qu’est-ce que je peux faire de ce “prix historique” quand on sait qu’à coté, tout augmente (comme dirait ma grand mère)… On pourrait normaliser par rapport a autre chose, par exemple un indice boursier (car on parle d’investissement), ou comme le font plutôt les économistes, avec quelque chose en rapport avec la masse monétaire. Non ?
Faisons les calculs: les séries du prix de l’or se trouvent ici par exemple (avec des séries journalières pour les années les plus récentes), et pour les données permettant d’actualiser, on pourra aller ici. Pour faire simple, je me suis contenté des données mensuelles en ligne ici, et pour la masse monétaire, M3 serait pas mal (je renvoie ici pour les détails), mais la Réserve Fédérale américaine ne le publie plus depuis 2006. J’aurais bien aimé aller ici, mais c’est payant. Finalement je suis allé , en extrapolant M3 sur les années récentes (de 2006 à 2010) en utilisant le taux de croissance observé sur M2. Les données ne sont pas corrigées des variations saisonnières.

>  OR = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/prix-or-mensuel.csv",
+ sep=";",header=TRUE)
>  tOR=as.vector(t(OR[,2:13]))
>  VOR=tOR[1:511]

car j’avais mis des 0 pour les mois non observés de 2010 (on les supprime tout simplement). Sinon on retrouve la courbe qui indique qu’effectivement, on a atteint des records…

>    TOR=ts(VOR, frequency = 12, start = c(1968, 1))
>    plot(TOR,col="red")

So what ? Regardons un peu notre masse monétaire qui permettra de ramener le prix de l’or dans le passé au prix réel de l’argent,

>    M = read.table("https://perso.univ-rennes1.fr/
arthur.charpentier/M1M2M3.txt")
>    M2=as.numeric(as.character(rev(M$V4)))
>    M3=as.numeric(as.character(rev(M$V5)))
>    M2b=M2[560:length(M2)]/M2[559]
>    M3[560:length(M3)]=M3[559]*M2b
>    TM3=ts(M3,frequency=12,start=c(1959,1))
>    plot(TM3,col="blue")

Bon, j’ai fait une approximation rapide, mais c’est vraiment histoire de voir ce qui se passe si on normalise le prix de l’or. Ceux qui veulent aller plus loin peuvent améliorer le code. Et en corrigeant pour travailler à quantité d’argent équivalent, on obtient des choses un peu différentes (j’aime les euphémismes)

>    TM3b=ts(M3[(12*9+1):length(M3)],
+    frequency = 12, start = c(1968, 1))
>    TORb=TOR/TM3b
>    plot(TORb,col="purple")
>    abline(h=TORb[length(TORb)],lty=2)

On a certes une hausse relativement importante depuis un peu moins de 10 ans du prix de l’or, mais le niveau n’est en rien comparable avec la crise qui s’est observé en 1980 (évoqué ici, , ou encore). Fort heureusement, on ne nous tient pas le même discours sur le PIB, nous expliquant qu’il ne cesse de battre, mois après mois, de nouveaux records ! On notera que ce graphique se rapproche d’autres, en ligne, ici ou
Pour ceux qui ne sont pas convaincu que l’or puisse être relié à une grandeur monétaire, on peut voir l’or comme une alternative à un placement en bourse, et comparer à un indice boursier (par exemple le SP500 pour rester aux États Unis),

>    SP=read.table("http://ichart.finance.yahoo.com/
table.csv?s=%5EGSPC&a=00&b=1&c=1968&d=08
&e=18&f=2010&g=m&ignore=.csv",
+       sep=",",header=TRUE)
>    tSP=rev(SP$Close)
>    TSP=ts(tSP, frequency = 12, start = c(1968, 1))
>    plot(TSP,col="red")

>    TORc=TOR/TSP
>    plot(TORc,col="purple")
>    abline(h=TORc[length(TORc)],lty=2)

Bref, quelle que soit la quantité retenue pour normaliser, on se rend compte que les niveaux atteints ces derniers mois sont très loin de ceux observés en 1980, voire dans les années 80. Pas de quoi en faire la une de la presse….

Marches aléatoires en dimension 1, 2 ou 3, ou histoire éthylique

Revenons à nos histoires d’alcoolique un instant (les précédentes sont en ligne iciet ). L’avant dernière fois (ici), nous avions parlé de la marche aléatoire, et du temps mis avant de trouver la sortie de la rue de la soif.

Et j’avais évoqué la probabilité la probabilité qu’il ne récupère pas les clés (i.e. qu’il ne repasse pas par l’origine) s’il peut sortir un jour de la rue. En probabilité, on appelle cela la probabilité de repasser par l’origine et cela vérifie une propriété étrange dès lors que l’on change de dimension… On va suivre notre alcoolique et noter https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy01.png sa position à la date https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy02.png, puis poser

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

autrement dit la probabilité de se retrouver à sa position initiale à la date n (et de retrouver ainsi ses clés), le https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png pouvait être en n’importe quelle dimension, on posera

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

correspondant à la probabilité de premier retour à la position initiale, et enfin

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

qui est la probabilité de revenir un jour au point de départ, c’est à dire https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png… C’est cette dernière grandeur qui nous intéresse….

  • De la probabilité de récupérer ses clés dans la rue

La rue est ici infinie (contrairement à l’autre fois, où on pouvait en sortir, soit vers la place de Lices, soit vers la place Sainte Anne). Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy08.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. S’il va à droite ou à gauche avec probabilités https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl99.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl98.png respectivement, notons que https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy11.png, car il faut un nombre de pas pairs pour revenir à l’origine, et

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

Il est possible de montrer que

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

et en particulier

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

On peut relier les https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy15.png et les https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy16.png en notant que

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

Or https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy16.png constitue une distribution de probabilité, dont on peut définir la fonction génératrice

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

On notera que que l’on peut écrire

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

soit, en utilisant l’expression donnée auparavant

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

Soit finalement

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

ce qui se traduit par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy21.png si la série https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy22.png converge. Aussi,

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

donc si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy26.png, notre alcoolique repassera presque sûrement par le point d’origine, alors que si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl97.png, il ne repasse pas par le bar où il a oublié ses clés avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy27.png|.

  • De la probabilité de récupérer ses clés sur une place

La place est ici infinie de part et d’autre. Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy09.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. Il a alors 1 chance sur 4 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite. On peut noter que là encore
https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy11.png, car il faut un nombre de pas pairs pour revenir à l’origine, et

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

Or on peut réécrire cette somme en notant que

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

et que

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

aussi,

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

que l’on peut approcher (via la formule de Stirling) sous le forme

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

La série https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy22.png est alors divergente, et donc https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethyl31.png. Aussi, même sur une place, notre alcoolique reviendra presque sûrement à sa place initiale (et ce une infinité de fois, car une fois qu’il y est revenu et qu’il en repart, il repassera presque sûrement à cette position, etc).

  • De la probabilité de récupérer ses clés dans l’espace

Bon, cette fois notre alcoolique passe en dimension 3…. j’admets, il faut de l’imagination. On s’intéresse à la probabilité qu’une marche aléatoire sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy10.png issue de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png repasse un jour par https://perso.univ-rennes1.fr/arthur.charpentier/latex/ethy04.png. Il a alors 1 chance sur 6 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite, d’aller en haut ou d’aller en bas.

Il y a désormais 6 voisins, et en raisonnant comme avant (le maillage est un peu plus complexe), on obtient que

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

Je laisse les amateurs de combinatoires faire des calculs ou des approximations, mais on peut montrer que

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

La série \sum u_n est alors convergente, et

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

qui vaut un peu plus de 1/3. Autrement dit, on a seulement 1 chance sur 3 de revenir à l’origine et de retrouver ses clés…
Bref, en dimension 1 ou 2, on repasse presque sûrement une infinité de fois à l’origine, mais en dimension 3, on a 2 chances sur 3 de ne jamais y repasser…. J’ai toujours trouvé ça étonnant….

Ramadan, entre le sexe et la faim…

Dans mon précédant billet (ici) je volais un dessin de Martin Vidberg qui parlait de google et de statistiques… Google trend est un outil magique pour noter des tendances et faire des statistiques simples et marrantes. Par exemple, puisque le ramadan touche à sa fin, on peut regarder ce qui motivent les algériens pendant le ramadan. Manifestement, le sexe n’est pas une priorité, que ce soit en anglais,

ou en français,

De manière générale, le mot clé le plus tapé sous google a tendance à se faire tout petit pendant la période du ramadan, alias porno,

En fait, on cherche un peu moins les femmes sur le net,

Par contre, visiblement, on a faim, et on veut bien manger: cuisine explose pendant le ramadan,

On peut aussi regarder les autres pays à forte population musulmane, et observer des tendances similaires, comme le Maroc,

les Émirats Arabes Unis,

ou encore le Pakistan,

(dans ce dernier cas, on observe un fort retour du mot clé sexune fois le ramadan passé… une espèce de manque virtuel). Bref, avec Google, on peut vraiment observer des choses intéressantes… non ?

"sendo l'intento mio scrivere cosa utile a chi la intende…"