Category Archives: Statistics

Modéliser et prévoir, entre ‘avoir faux’ et ‘avoir pas de chance’

Depuis plusieurs heures, il y a une phrase qui revient sans cesse, et qui m’agace au plus haut point, sur “le modèle est faux…“. C’est une phrase que j’ai entendu pas plus tard qu’au séminaire ce midi, lorsqu’un collègue a dit que très clairement les modèles étaient faux puisqu’ils n’avaient pas prédit le vainqueur des élections américaines. Par exemple, le jour des élections, Huffington Post annonçait que Donald Trump avait 1,6% de chances d’être élu…

et pourtant, il a gagné. C’est donc pourtant que le modèle était “faux”, c’est évident, non ?

Pourtant, tout ceux qui ont fait un minimum de statistique dans leurs vies le savent… les évènements rares, ça peut se produire,

> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "trump"

Ce type de raisonnement me dérange car il dénote une mauvaise compréhense de ce qu’est une probabilité. Autre exemple: supposons que j’annonce que ma belle mère a une chance sur un million de mourrir lorsqu’elle prendra la route pendant les prochaines vacances, et qu’elle meurt effectivement cet hiver d’un accident de voiture. Comme pour les élections américaines, on sera surpris, sous le choc. Mais est-ce que ça veut dire que mon modèle était faux ? Non…Si j’avais annoncé qu’elle avait 85% de ne pas passer les vacances en prenant la voiture, j’aurais eu un modèle des plus douteux… Le problème ici c’est qu’on avance des probabilités sur des évènements qu’on ne pourra observer qu’une fois. On ne peut pas le tester suffisement pour savoir s’il est juste, en moyenne.

Cela dit, c’est un problème assez profond, car il est très dur de savoir ce qui relève effectivement de la malchance… C’est en tous cas comme ça que je comprends le strip ci-dessous

Et je trouve dommage que sous prétexte qu’un évènement rare se réalise, on blâme la personne en charge de la prévision… Il y a manifestement une incompréhension profonde de ce que peut être le hasard. C’est triste…

Ce que la courbe ROC (et l’AUC) ne raconte pas

En préparant une intervention pour mardi prochain, j’épluchais les résultats renvoyés pour un exercice, et j’ai eu un résultat assez étrange avec un modèle de classification. J’avais donné la même base cet automne à l’ensae, et j’avais donc près d’une trentaine d’autres modèles, pour comparer (disons plutôt que sur la même base de test, j’ai une trentaine de prévisions). Les observations noires sont celles obtenues cet automne (le trait correspond aux meilleurs AUC sur la base de test), et les observations rouges sont celles que j’ai obtenu pour l’intervention de mardi (là encore, le trait vertical correspond aux meilleurs modèles, au sens du AUC), sur une observation de la base de test,

Ce sont les probabilités prédites (mais j’ai enlevé l’échelle).

Pour presque toutes mes observations, les poids rouges sont bien au dessus des autres… Mais cela n’enlève en rien le fait que l’AUC obtenu (pour les deux modèles rouges) est très bon. C’est effectivement un résultat important (et connu) : le critère AUC (mais plus généralement la courbe ROC en entier) n’indique en rien si la valeur prédite est bonne, ou pas. Il nous dit juste si l’ordre obtenu est correct. Si les valeurs les plus importantes sont effectivement les valeurs pour lesquelles on a un 1, l’AUC sera très bon.

C’est ce qu’on peut observer sur le petit exemple ci-dessous. Considérons un modèle logistique simulé assez simple,

> n=1e3
> set.seed(1)
> x1=rnorm(n)
> x2=runif(n)
> u=-3+x2+x1
> p=exp(u)/(1+exp(u))
> y=rbinom(n,prob=p,size=1)
> library(ROCR)
> df=data.frame(y,x1,x2)
> mean(df$y)
[1] 0.116
> reg=glm(y~.,data=df,family=binomial)
> p=predict(reg,type="response")
> mean(p)
[1] 0.116
> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")

L’AUC est ici

> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

et la courbe ROC est la suivante,

> plot(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="blue")

Supposons que l’on change la constante du modèle logistique,

> reg$coefficients[1]=0

Dans ce cas, notre prévision est assez mauvaise, car la probabilité moyenne prédite est ici

> u=reg$coefficients[1]+reg$coefficients[2]*
+ df$x1+reg$coefficients[3]*df$x2
> p=exp(u)/(1+exp(u))
> mean(p)
[1] 0.6060676

(on est loin des 11,6% de 1 dans la base). Pourtant le AUC est bon

> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")
> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

(c’est en fait la même valeur qu’auparavant, ce qui a du sens puisque la courbe ROC est identique)

> lines(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="red")

Autrement dit, ces outils, classiquement utilisés pour juger la qualité d’un classifieur, ne permettent en aucun cas de dire que la probabilité prédite a du sens. Ces critères permettent juste de dire qu’on identifie assez bien les personnes qui ont le plus de chance d’avoir la réponse 1. Ce qui n’est pas si mal… mais c’est un autre problème que celui d’avoir une probabilité qui soit pertinente.

Estimer une Proportion sur une Question Sensible, par Enquête

Comme le rapportait Slate hier (dans un article qui, malheureusement ne mérite pas d’être lu), un joli sujet a été proposé au baccalauréat au lycée français de Pondichéry, en mathématiques.

Ce sujet est d’autant plus intéressant qu’il revient sur une méthode assez classique pour questionner les gens sur des questions sensibles (ici sur le téléchargement pirate, mais on peut poser la même question sur la fraude, ou des pratiques sexuelles).

Comme le notait , cette technique avait été introduite initialement par Stanley Warner, voilà un demi-siècle,

Cette technique est particulièrement intéressante, tant d’un point de vu intellectuel que pratique. Plusieurs sur Twitter me faisaient remarquer que, contrairement à ce qui était raconté dans Slate, cette technique est effectivement utilisée.

Quand j’enseignais cette technique, je prenais la version simple, avec une pièce de monnaie. Supposons qu’une proportion  de personnes fraudent le fisc. Et qu’une proportion  aime Prince, va à la piscine une fois par mois, a un diplôme de l’enseignement supérieur, etc. Peut importe, la seconde proportion doit pouvoir être estimée sans problème avec une autre enquête (voire par le recensement). Bref, on a un estimateur de . Ici, les personnes interrogées tirent une pièce et répondent soit à la question “avez-vous déjà fraudé au Fisc?” s’ils tombent sur face, et à la question “aimez-vous Prince?” s’ils tombent sur pile. On pourrait dire qu’ils tirent une pièce dans leur tête, mais il faut s’assurer que la probabilité de tomber sur face est effectivement 1/2.

La magie est de noter que la proportion de personnes qui fraudent le fisc peut s’écrire

On peut avoir, par une autre enquête un estimateur du second terme, . Et en regardant des personnes qui ont répondu “oui” à une question inconnue, on a un estimateur pour

Et par différence, on en déduit un estimateur de . Les maths sont dans l’article de Stanley Warner, mais il faut reconnaître que la méthode est vraiment jolie…. Dans le cas du sujet du bac, la mise en oeuvre est un peu différente, mais l’idée est un peu la même….

(si certains ont des exemples pratiques d’application, les commentaires sont ouverts, ou sinon par mail pour joindre des documents en pièce jointe)

Radial Graphs for Time Series

On How to: Weather Radials, there was a nice visualisation of temperatures. Since I am too old fashioned for ggplot2, I wanted to reproduce a similar graph with the old plot style.

Assume that daily temperature is in a vector X (e.g. temperature in Montréal, QC, in 2009). To get a radial plot, use

> n=length(X)
> theta=seq(0,1-1/n,length=n)*2*pi
> r=30+X
> plot(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",xlab="",ylab="",axes=FALSE)
> for(t in 1:n){
+   if(X[t]>0) CL=rgb(0,0,1,.4)
+   if(X[t]<0) CL=rgb(1,0,0,.4)
+   if(X[t]==0) CL="white"
+   segments((30+X[t])*cos(pi/2-theta[t]),(30+X[t])*sin(pi/2-theta[t]),30*cos(pi/2-theta[t]),30*sin(pi/2-theta[t]),col=CL)
+ }
> for(r in 10*seq(0,6)) lines(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",col="light blue")

Bonferroni, les anniversaires et les terroristes (ou pas)

Le paradoxe des anniversaires est un exercice de probabilité que l’on découvre assez vite. Mais assez souvent, on l’oublie, alors qu’il est en réalité assez profond. Le problème est le suivant: quel est le nombre de personnes qu’il faut réunir pour avoir (au moins) une chance sur 2 pour que deux personnes aient leur anniversaire le même jour. Si on calcule la probabilité que deux personnes dans un groupe de n fêtent leur anniversaire le même jour, on peut calculer la probabilité que personne ne soit né le même jour, parmi n personnes. Un calcul rapide montre qu’avec 23 personnes, la probabilité qu’au moins deux personnes fête son anniversaire le même jour est 50.73%. Autrement dit, il suffit d’une groupe de 23 personnes (c’était la taille d’une classe d’élèves, du moins à mon époque…).

L’intuition, pour comprendre pourquoi nous avons l’impression que ce nombre est petit est que, naturellement, si on a un événement qui survient avec probabilité p, et qu’on répète l’expérience n fois, la probabilité que rien ne se produise est (1-p)n. Aussi, en première approximation (et j’insiste sur ce point parce que c’est avec ce genre d’argument fallacieux qu’on conclue que la probabilité d’avoir une catastrophe nucléaire dépasse 100%, conduisant à du “certitude statistique“), on en déduit que la probabilité qu’un événement survienne est de l’ordre de np. Sauf qu’il, on compte des paires. Et avec n personnes, on constitue environ n2/2 paires. Et donc la probabilité devient proche de n2p/2… ce qui est bien plus grand que np. C’est grossièrement ce qui explique que notre intuition se sent trompée : avec seulement 23 personnes, on peut constituer plus de 250 paires. La probabilité que toutes les paires correspondent à des dates d’anniversaire différente n’est pas si petit qu’il y parait.

Jeff Ullman parle ainsi de principe de Bonferroni,

“if you look in more places for interesting patterns than your amount of data will support, you are bound to find crap”

Et cette intuition faussée quand on travaille sur des paires se retrouve dans les techniques utilisées par les services de renseigngement, tel que le montre Anand Rajaraman dans Mining Massive DatasetsSupposons qu’il existe un milliard d’individus qui pourraient potentiellement être des terroristes. Que les gens, quels qu’il soient vont à l’hôtel 1 jour sur 100. Un hôtel donné peut accueillir 100 personnes. On suppose qu’il existe 100,000 personnes, de sorte qu’en moyenne ils soient plein tous les jours. Et supposons qu’on regarde sur un historique de 3 ans, i.e. 1000 jours pour avoir des calculs ronds. Supposons que pour traquer les terroristes, on cherche les gens qui ont séjourné, ensemble, dans un hôtel, le même jour, au moins deux fois… On imagine que c’est effectivement rare, et donc suspicieux…

La probabilité que deux personnes données dorment à l’hôtel le même jour est faible, 1/100×1/100, soit 0,0001. Et comme il y a 100,000 hôtel, la probabilité que les deux personnes descendent, en plus, dans le même hôtel, il faut diviser 0,0001 par 100,000. Ce qui donne 10-9 . Mais ça, c’est la probabilité pour qu’ils se croisent une fois ! Si on regarde la probabilité qu’ils se croisent deux fois, en trois ans, on arrive à une probabilité de 10-18 Pour ceux qui ont du mal, ça doit faire environ

0,00000000000000001 %

Et encore. Ici il peut en plus se croiser la seconde fois dans un autre hôtel. La probabilité qu’il se croisent dans le même hôtel que la première fois devient encore plus petite !

Supposons maintenant qu’un organisme de surveillance avisé décider de lister les paires de personnes qui, sur des paires de jours, se sont retrouvé dans le même hôtel. Ça doit être suffisamment rare pour être suspicieux. Rappelons que quand on considère n individus, le nombre de paires possibles est de l’ordre de n2/2. Autrement dit, avec un million de personnes, le nombre de paires est de l’ordre de 5×1017. Et pareil, le nombre de paires de jours est de l’ordre de 5×105. Aussi, le nombre d’événements suspicieux que l’on observe, sur ces trois années est de l’ordre de

(5×1017 ) x (5×10) x10-18 = 250000

Autrement dit, 1/4 de millions de paires de personnes sont des terroristes potentiels ! Soit 1/2 millions d’individus.

Maintenant, supposons qu’il y ait 10 terroristes. S’ils sont assez stupides pour se retrouver dans les mêmes hôtels, les mêmes jours, ça veut dire que sur les 500000 personnes qui vont être surveillées, 10 sont effectivement des terroristes, et on a 499990 faux positifs… Enfin, “faux positifs”… ils ont forcément quelque chose d’autre à se reprocher, non ?

Dynamique de la Pyramide des Ages

Très joli billet sur blog.revolutionanalytics.com avec un code de @kyle_e_walker permettant, très simplement (moyennant une inscription pour avoir une clé permettant d’utiliser l’API du census) de construire une pyramide des âges dynamiques.

> devtools::install_github('walkerke/idbr')
> library(idbr)
> library(ggplot2)
> library(animation)
> library(dplyr)
> library(ggthemes)
> idb_api_key("mykey1239F2f324zf9GGZgege32R2ii4")

On importe alors les données pour les hommes et les femmes,

> male <- idb1('FR', 2010:2050, sex = 'male') %>%
+   mutate(POP = POP * -1,
+   SEX = 'Male')
> female <- idb1('FR', 2010:2050, sex = 'female') %>% mutate(SEX = 'Female')

et on stocke le tout

> france <- rbind(male, female) %>%
+   mutate(abs_pop = abs(POP))

Ensuite, on crée l’animation,

> saveGIF({
+   
+   for (i in 2010:2050) {
+     
+     title <- as.character(i)
+     
+     year_data <- filter(france, time == i)
+     
+     g1 <- ggplot(year_data, aes(x = AGE, y = POP, fill = SEX, width = 1)) +
+       coord_fixed() + 
+       coord_flip() +
+       annotate('text', x = 98, y = -800000, 
+       label = 'Data: US Census Bureau IDB; idbr R package', size = 3) + 
+       geom_bar(data = subset(year_data, SEX == "Female"), stat = "identity") +
+       geom_bar(data = subset(year_data, SEX == "Male"), stat = "identity") +
+       scale_y_continuous(breaks = seq(-1000000, 1000000, 500000),
+       labels = paste0(as.character(c(seq(1, 0, -0.5), c(0.5, 1))), "m"), 
+       limits = c(min(france$POP), max(france$POP))) +
+       theme_economist(base_size = 14) + 
+       scale_fill_manual(values = c('#ff9896', '#d62728')) + 
+       ggtitle(paste0('Population structure of France, ', title)) + 
+       ylab('Population') + 
+       xlab('Age') + 
+       theme(legend.position = "bottom", legend.title = element_blank()) + 
+       guides(fill = guide_legend(reverse = TRUE))
+     print(g1) 
+   }
+ }, movie.name = 'france_pyramid.gif', interval = 0.1, ani.width = 700, ani.height = 600)

Et le résultat est vraiment joli, non ?

Mortality by Weekday and Age

A few days ago, I did mention on Twitter a nice graph, with

My colleague Jean-Philippe was extremely sceptical, so I tried to reproduce that graph. The good thing is that we have the Social Security Death Master File, for data in the US. To be more specific, I have three big files on my hard drive, and in order to reproduce that graph, we’ll load the data by chunks. But before, because we have the day of birth, and the day of death, I need a function to compute the age. So here it is

> age_years <- function(earlier, later)
+ {
+   lt <- data.frame(earlier, later)
+   age <- as.numeric(format(lt[,2],format="%Y")) - as.numeric(format(lt[,1],format="%Y"))
+   dayOnLaterYear <- ifelse(format(lt[,1],format="%m-%d")!="02-29",
+                            as.Date(paste(format(lt[,2],format="%Y"),"-",format(lt[,1],format="%m-%d"),sep="")),
+                            ifelse(as.numeric(format(later,format="%Y")) %% 400 == 0 | as.numeric(format(later,format="%Y")) %% 100 != 0 & as.numeric(format(later,format="%Y")) %% 4 == 0,
+                                   as.Date(paste(format(lt[,2],format="%Y"),"-",format(lt[,1],format="%m-%d"),sep="")),
+                                   as.Date(paste(format(lt[,2],format="%Y"),"-","02-28",sep=""))))
+   age[which(dayOnLaterYear > lt$later)] <- age[which(dayOnLaterYear > lt$later)] - 1
+   age
+ }

from github.com/nzcoops. Now, it is possible to create a similar table, based on that huge file (we have almost 50 million observations)

> cols <- c(1,9,20,4,15,15,1,2,2,4,2,2,4,2,5,5,7)
> noms_col <- c ("code","ssn","last_name","name_suffix","first_name","middle_name","VorPCode","date_death_m","date_death_d","date_death_y","date_birth_m","date_birth_d","date_birth_y","state","zip_resid","zip_payment","blanks")
> library(LaF)

> TABLE_AGE_DAY=function(temp = "ssdm3"){
+ ssn <- laf_open_fwf( temp,column_widths = cols,column_types=rep("character",length(cols) ),column_names = noms_col,trim = TRUE)
+ object.size(ssn)
+ go_through <- seq(1,nrow(ssn),by = 1e05 )
+ if(go_through[ length(go_through)] != nrow( ssn)) go_through <- c(go_through,nrow( ssn))
+ go_through <- cbind(go_through[-length(go_through)],c(go_through[-c(1,length(go_through)) ]-1,go_through [ length(go_through)]))
+ go_through
+ 
+ pb <- txtProgressBar(min = 0, max = nrow( go_through), style = 3)
+ count_birthday <- function(s){
+   #print(s)
+   setTxtProgressBar(pb, s)
+   data <- ssn[ go_through[s,1]:go_through[s,2],c("date_death_y","date_death_m","date_death_d",
+                                                  "date_birth_y","date_birth_m","date_birth_d")]
+   date1=as.Date(paste(data$date_birth_y,"-",data$date_birth_m,"-",data$date_birth_d,sep=""),"%Y-%m-%d")
+   date2=as.Date(paste(data$date_death_y,"-",data$date_death_m,"-",data$date_death_d,sep=""),"%Y-%m-%d")
+   idx=which(!(is.na(date1)|is.na(date2)))
+   date1=date1[idx]
+   date2=date2[idx]
+   itg=try(age<-age_years(date1,date2),silent=TRUE)
+   if(inherits(itg, "try-error")) age=trunc((date2-date1)/365.25)
+   w=weekdays(date2)
+   T=table(age,w)
+   Tab=matrix(0,106,7)
+   for(i in 1:nrow(T)) if(as.numeric(rownames(T)[i])<106) Tab[as.numeric(rownames(T)[i]),]=T[i,]
+   return(Tab)
+ }
+ D <- lapply( seq_len(nrow( go_through)),count_birthday) 
+ T=D[[1]]
+ for(s in 2:length(D)) T=T+D[[s]]
+ return(T)
+ }

If we run that function on the three files

> D1=TABLE_AGE_DAY("ssdm1")
|========================================| 100%
> D2=TABLE_AGE_DAY("ssdm2")
|========================================| 100%
> D3=TABLE_AGE_DAY("ssdm3")
|========================================| 100%

we can visualize not percentages, as on the figure above, but counts

> D=D1+D2+D3
> colnames(D)=
c("Sun","Thu","Mon","Tue","Wed","Sat","Fri")
> D=D1[,
c("Sun","Mon","Tue","Wed","Thu","Fri","Sat")]

and we have here (I remove the Saturday to get a better output)

> D[,1:6]
          Sun    Mon    Tue    Wed    Thu    Fri
  [1,]   2843   2888   2943   3020   2979   3038
  [2,]   2007   1866   1918   1974   1990   2137
  [3,]   1613   1507   1532   1530   1515   1613
  [4,]   1322   1256   1263   1259   1207   1330
  [5,]   1155   1061   1092   1128   1112   1171
  [6,]   1067    985    950   1082   1009   1055
  [7,]   1129    901    915    954    941   1044
  [8,]   1026    927    944    935    911   1005
  [9,]   1029   1012    871    908    939    998
 [10,]   1093   1011    974    958    928   1018
 [11,]   1106   1031   1019   1036   1087   1122
 [12,]   1289   1219   1176   1215   1141   1292
 [13,]   1618   1455   1487   1484   1466   1633
 [14,]   2121   2000   1900   1941   1845   2138
 [15,]   2949   2647   2519   2499   2524   2748
 [16,]   4488   3885   3798   3828   3747   4267
 [17,]   5709   4612   4520   4422   4443   5005
 [18,]   7280   5618   5400   5271   5344   5986
 [19,]   8086   6172   5833   5820   6004   6628
 [20,]   8389   6507   6166   6055   6430   6955
 [21,]   8794   7038   6794   6628   6841   7572
 [22,]   8578   6528   6512   6472   6757   7342
 [23,]   8345   6750   6483   6469   6714   7338
 [24,]   8361   6859   6589   6623   6854   7369
 [25,]   8398   6974   6892   6766   6964   7613
 [26,]   8432   7210   7012   7175   7343   7801
 [27,]   8757   7641   7526   7352   7674   7950
 [28,]   9190   8041   7843   7851   7940   8268
 [29,]   9495   8409   8555   8400   8469   8934
 [30,]   9876   9041   9015   9166   9106   9641
 [31,]  10567   9952   9506   9634   9770  10212
 [32,]  11417  10428  10402  10275  10455  11169
 [33,]  11992  11306  11124  11095  11243  11749
 [34,]  12665  12327  11760  12025  12137  12443
 [35,]  13629  13135  13179  13037  12968  13724
 [36,]  14560  14009  13927  13822  14105  14436
 [37,]  15660  14990  15013  15009  15101  15700
 [38,]  16749  16504  16148  16091  15912  16863
 [39,]  17815  17760  17519  17144  17553  17943
 [40,]  19366  19057  18918  18517  18760  19604
 [41,]  20770  20458  20154  20339  20349  21238
 [42,]  21962  22194  22020  21499  21690  22347
 [43,]  23803  23922  23701  23681  23437  24227
 [44,]  25685  26133  25559  25209  25287  26115
 [45,]  27506  28110  27363  27042  27272  28228
 [46,]  29366  29744  29555  29245  29678  30444
 [47,]  31444  32193  31817  31504  31753  32302
 [48,]  33452  34719  33529  33954  33441  34618
 [49,]  36186  37150  36005  36064  36226  37138
 [50,]  38401  39244  38813  38465  38506  39884
 [51,]  40331  41830  41168  41110  40937  42014
 [52,]  43181  44351  43975  43949  43579  44734
 [53,]  45307  47134  46522  46149  46089  47286
 [54,]  47996  49441  49139  48678  48629  49903
 [55,]  50635  52424  51757  51433  51477  52550
 [56,]  53509  55337  54556  54482  54406  55906
 [57,]  55703  58482  58016  57400  57097  58758
 [58,]  59016  61453  60652  61024  60557  62473
 [59,]  62475  65651  64169  63824  63829  65592
 [60,]  66621  69185  68885  68217  68752  69963
 [61,]  69759  73144  72421  71784  71745  73414
 [62,]  80346  84253  83044  83177  82416  83833
 [63,]  86851  90059  89002  88985  89245  90334
 [64,]  91839  95465  94602  93985  94154  96195
 [65,]  98461 102846 101348 101328 101306 103170
 [66,] 104569 108722 107768 107711 107729 109350
 [67,] 111230 115477 114418 114743 113935 116356
 [68,] 116999 122053 120727 120342 119782 122926
 [69,] 123695 128339 127184 126822 126639 129037
 [70,] 129956 136123 134555 135120 133842 137390
 [71,] 137984 142964 141316 142855 141419 143620
 [72,] 145132 150708 148407 149345 149448 151910
 [73,] 152877 157993 155861 156349 155924 158725
 [74,] 159109 164652 162722 163499 163157 165744
 [75,] 165848 172121 170730 170482 170585 173431
 [76,] 172457 179036 177185 177328 177392 180215
 [77,] 179936 185015 183223 183932 183237 186663
 [78,] 185900 191053 189986 189730 189639 193038
 [79,] 191498 196694 194246 194810 195246 197812
 [80,] 195505 201289 199684 199561 198968 203226
 [81,] 199031 204927 202204 202622 202951 205792
 [82,] 201589 207928 204929 204001 204396 208224
 [83,] 201665 206743 205194 204676 205256 207980
 [84,] 200965 205653 203422 202393 203422 206012
 [85,] 197445 202692 199498 199730 200075 201728
 [86,] 192324 195961 193589 194754 193800 196102
 [87,] 183732 188063 185153 186104 186021 188176
 [88,] 174258 177474 175822 176078 176761 177449
 [89,] 163180 166706 162810 164367 164281 166436
 [90,] 149169 151738 150148 150212 150535 152435
 [91,] 134218 136866 134959 134922 135027 136381
 [92,] 118936 121106 119591 119509 119793 120998
 [93,] 102734 104955 102944 102865 103345 104776
 [94,]  87418  88885  88023  86963  87546  87872
 [95,]  72023  72698  72151  71579  71530  72287
 [96,]  56985  58238  57478  57319  57163  57615
 [97,]  44447  45058  44607  44469  43888  44868
 [98,]  33457  34132  33022  33409  33454  33642
 [99,]  24070  24317  24305  24089  24020  24383
[100,]  17165  17295  16755  17115  16957  17207
[101,]  11799  12125  11709  11816  11824  11719
[102,]   7714   7741   7959   7691   7648   7633
[103,]   5024   5012   4822   4792   4882   4916
[104,]   2987   3101   2978   3049   3093   2906
[105,]   1781   1894   1811   1756   1734   1834

So clearly, for young people, the number of deaths is rather small…

And to visualize it, as above, we can use

> P=D/apply(D,1,sum)*100
> range(P)
[1] 12.34857 17.59386
> dP=trunc((P-min(P))/(max(P)+.01-min(P))*11)
> library(RColorBrewer)
> CLR=rev(brewer.pal(name="RdYlBu", 11))

> plot(0:1,0:1,ylim=c(55,110),xlim=c(-1,7))
> for(i in 1:106){
+   for(j in 1:7){
+  rect(j-1,108-i,j,107-i,col=CLR[dP[i,j]])
+   }}
> text(rep(-.5,106),107.5-1:106,0:105,cex=.4)

As above, we observe a strong difference among weekdays for the date of death for young people (below 30) which disappear after (even if there is still a sunday effect)

La guerre des étoiles : distinguer le signal du bruit

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

Approches historiques de cette notion de « significativité »

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

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

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

Tests, prise de décision et erreurs

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

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

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

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

Figure 3 : calculs des taux d’erreurs

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

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

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

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

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

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

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

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

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

Le problème des tests multiples

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

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

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

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

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

Les alternatives pour juger de la significativité

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

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

Références

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

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

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

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

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

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

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

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

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

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

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

Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that  has a similar distribution… so there is no reason to test for a Gamma model for  before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if  has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution  is Gaussian but also the non-conditional distribution of .

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of  is also Gaussian

> library(nortest)
> ad.test(Y)

	Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

	Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here  is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value}
 
> plot(2:100,Vectorize(pv)(2:100),type="l")
> abline(h=.05,col="red")

So here, it could be possible to have also a Gaussian distribution, for . As least to accept that assumption, statistically.

In the context of a Poisson regression, it is well know that it’s not possible to have at the same time  that is Poisson distributed (that’s a Poisson regression) and also  that is Poisson distributed. That simply comes from the fact that

while

and because of the conditional Poisson distribution, then

Thus,

So  cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for .

More generally, it is very difficult to have a distribution family for   that is also the distribution of the non-conditional variable . In the context of a finite mixture ( takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference).

Hence, as a keep saying, before running a regression model on with some given family, it is never a good idea to check if the non-conditional distribution  has the same distribution. Because there is no reason, usually, to remain in the same family.

Confidence Regions for Parameters in the Simplex

Consider here the case where, in some parametric inference problem, parameter  is a point in the Simplex,

For instance, consider some regression, on compositional data,

> library(compositions)
>  data(DiagnosticProb)
>  Y=DiagnosticProb[,"type"]-1
>  X=DiagnosticProb[,c("A","B","C")]
>  model = glm(Y~ilr(X),family=binomial)
>  b = ilrInv(coef(model)[-1],orig=X)
>  as.numeric(b)
[1] 0.3447106 0.2374977 0.4177917

We can visualize that estimator on the simplex, using

>  tripoint=function(s){
+    p=s/sum(s)
+    abc2xy(matrix(p,1,3))
+  }

>  lab=LETTERS[1:3]
>  xl=c(-.1,1.25)
>  yl=c(-.1,1.15)
>  library(trifield)
>  A=abc2xy(matrix(c(1,0,0),1,3)) 
>  B=abc2xy(matrix(c(0,1,0),1,3))
>  C=abc2xy(matrix(c(0,0,1),1,3)) 
>  plot(0:1,0:1,col="white",
+  xlim=xl,ylim=yl,xlab="",ylab="",axes=FALSE)
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

If we want to compute a ‘confidence region’, we can either use Bayesian models (with a Dirichlet distribution as prior distribution), or use bootstrap. We will use here the second idea

>  MB=matrix(NA,1e4,2)
>  for(sim in 1:1e4){
+    idx=sample(1:nrow(DiagnosticProb),
+    size=nrow(DiagnosticProb),replace=TRUE)
+  Y=DiagnosticProb[idx,"type"]-1
+  X=DiagnosticProb[idx,c("A","B","C")]
+  model = glm(Y~ilr(X),family=binomial)
+  MB[sim,]=tripoint(as.numeric(
+    ilrInv(coef(model)[-1],orig=X)))}

To get some ‘confidence region’, we can then use the bagplot, to get either a region where 50% of the boostraped estimators are, or 95%,

>  library(aplpack)
> P1=bagplot(MB[,1],MB[,2], factor =1.96, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 
> P2=bagplot(MB[,1],MB[,2], factor =0.67, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 

Then we can easily plot those two regions,

>  plot(0:1,0:1,col="white")
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  polygon(P1$hull.loop,col="yellow",border=NA)
>  polygon(P2$hull.loop,col="orange",border=NA)
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

 

Collectionner les Figurines Lego avec ses Ami(e)s

Dans mon précédant billet, je regardais combien de figurines Lego je devais acheter pour avoir toute la série, avec ma fille. Pour la série de 16, il fallait, en moyenne, acheter 54 sachets, à 3 euros pièces. Dans les commentaires, @Guillaume et @Gaelle me faisaient noter que les échanges, ça pouvait aider.

Bon, tout d’abord j’aime pas qu’on donne des conseils à mes filles pour avoir raison contre moi ! Parce que oui, c’est sûr qu’en s’y mettant à plusieurs, le coût individuel sera plus faible… Ensuite, parce que les maths me semblent plus compliqués…. Cela dit, en simulant des parties, ça tombe assez vite…

Pour commencer, je fais une petite fonction (histoire d’alléger mes codes) pour incrémenter mes comptages de figurines, à chaque fois que j’en achète une

> f=function(mat,i){
+   mat[i]=mat[i]+1
+   return(mat)
+ }

Ensuite, on peut avoir une petite fonction, correspondant à une simulation, qui me renvoie le nombre de figurine de chaque type que j’ai du acheter pour avoir les 16 une fois. Ou plus. En fait, je veux avoir en tout k fois toutes les figurines (ce qui permet de partager parmi k ami(e)s)

> tirage=function(k=1){
+ test=FALSE
+ cpte=0
+ M=rep(0,16)
+ while(!test){
+   no=sample(1:16,size=1)
+   M=f(M,no)
+   cpte=cpte+1
+   if(min(M)==k) test=TRUE
+ }
+ return(list(compte=M,
+             nbre=cpte))
+ }

Maintenant on peut faire des simulations. Par exemple, pour retrouver le fait qu’il faut, en moyenne, acheter 54 sachet pour que ma fille (seule) ait toutes les figurines

> N=unlist(lapply(1:10000,function(x) tirage(1)$nbre))
> mean(N)
[1] 54.1854

Si maintenant on le fait avec mes deux filles,

> N=unlist(lapply(1:10000,function(x) tirage(2)$nbre))
> mean(N)
[1] 82.7961

ce qui est (beaucoup) moins que le double. Si au lieux de faire ça entre mes deux filles, mais que ma fille collectionne les figurines avec une amie à elle (et que les papas partagent les frais), on a en moyenne 41 sachets à acheter, chacun,

> mean(N)/2
[1] 41.39805

et pour trois, on passe à 36 sachets

> N=unlist(lapply(1:10000,function(x) tirage(3)$nbre))
> mean(N)/3
[1] 36.04363

etc

> V=unlist(lapply(1:10,function(k) mean(unlist(lapply(1:10000,function(x) tirage(k)$nbre)))))
> plot(1:10,V/1:10,type="b",pch=19)

On peut aller plus loin, en regardant les quantiles. Avec 4 ami(e)s, si on achète tous 45 sachets, on a 90% de chances que tout le monde ait les 16 figurines

D’ailleurs, si on regarde l’évolution du quantile (à 90%) sur la moyenne, on voit qu’il baisse vite avec la taille…. Avec 2 ami(e)s, il faut 35% de sachets en plus pour s’assurer d’avoir toutes les figurines dans 90% de cas (par rapport au nombre moyen) alors qu’il en faut juste 20% de plus avec 8 ami(e)s.

Maintenant, si quelqu’un a le courage (et le temps) de faire les maths, je suis preneur !

Collectionner les Figurines Lego

 

Tous ceux qui ont des enfants – ou qui jouent encore au Lego – le savent : régulièrement, Lego sort des séries de figurines, vendus en sachets individuels hermétiques, un peu comme les anciennes cartes à jouer. On est à la saison 15, et il y a 16 figurines par séries,

D’ordinaire, les enfants réclament un paquet, quand on va dans un magasin de jeux (genre pour acheter un cadeau pour un anniversaire) mais là, ma fille les a réclamé tous.

Bon, après le moment où je ris dans ma barbe, j’ai quand même fini par me demander combien ça pourrait coûter si on achetait des paquets, de manière à récupérer les 16 figurines.

La bonne stratégie pour faire le calcul me semble d’ordonner les personnage par ordre d’apparition,  et de modéliser le nombre de paquets achetés entre deux premières apparitions, comme indiqué sur le dessin ci-dessous

Si  désigne le nombre total de paquets à acheter, on note  le nombre de paquets à acheter, entre l’apparition de la ()ème nouvelle figurine, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première figurine), on a

On notera que, si  désigne la loi géométrique,
et à partir de là, les calculs sont simples, puisque

soit (par linéarité de l’espérance)

i.e.

ou encore

Le nombre moyen de paquets à acheter est alors
> sum(16/(1:16)) 
[1] 54.09166

Pour ceux qui ont oublié leur calculatrice, on peut se souvenir que

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)
> log(16)*16 
[1] 44.36142

qui diffère du résultat numérique car il manquait la constante d’Euler

https://latex.codecogs.com/gif.latex?gamma%20=%20lim_{n%20<br /> ightarrow%20infty%20}%20left(%201+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20-%20ln(n)%20<br /> ight)

i.e.

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)+gamma

Numériquement, on obtient, en prenant pour https://latex.codecogs.com/gif.latex?gamma 0.57721 (cf Google)

> (log(16)+ 0.57721)*16 
[1] 53.59678

Mais je m’écarte du sujet… à 3€, il faut débourser, en moyenne, 162€ pour avoir les 16 personnages, soit plus de 10€ (en moyenne) par personnage !

> sum(16/(1:16))*3 
[1] 162.275

On est ici au même prix que l’hotel Lego friends (oui, c’est une demande récurrente)

Par simulation, on peut confirmer cet ordre de grandeur (avec ici un million de simulation de la variable),

> rpaquet=function(i) sum(1+rgeom(16,prob=(1:16)/16))
> S=unlist(lapply(1:1e6,rpaquet))
> mean(3*S)
[1] 162.1709

mais surtout, on peut avoir des quantiles, par exemple

> median(S)*3
[1] 150

après avoir déboursé 150 euros, il y a 50% de chances que l’on ait déjà toutes les figurines.

> quantile(S,.9)*3
90% 
237

 

Et il y a 1 chance sur 10 pour qu’on dépasse les 237€. Mine de rien, pour 16 petites figurines, on dépasse là le prix du Falcon Millenium,

Je pense que les enfants ne sont pas prêt que je leur rachète des petites figurines !

Regression with Splines: Should we care about Non-Significant Components?

Following the course of this morning, I got a very interesting question from a student of mine. The question was about having non-significant components in a splineregression.  Should we consider a model with a small number of knots and all components significant, or one with a (much) larger number of knots, and a lot of knots non-significant?

My initial intuition was to prefer the second alternative, like in autoregressive models in R. When we fit an AR(6) model, it’s not really a big deal if most coefficients are not significant (but the last one). It’s won’t affect much the forecast. So here, it might be the same. With a larger number of knots, we should be able to capture small bumps that we’ll never capture with a smaller number.

Here is what a have with a small number of knots, and cubic splines

and with a larger number of knots

In order to understand what’s going on, consider a simple model, with the two splines above, in red

> set.seed(1)
> library(splines)
> x=seq(0,1,by=.01)
> v=bs(x,10)
> x2=v[,2]
> x10=v[,10]
> set.seed(1)
> y=1+3*x2+5*x10+rnorm(length(x))/4
> y_test=1+3*x2+5*x10+rnorm(length(x))/4

Note that here I have generated two sets of data, one to train a model, and one to test it.  Here, the data looks like that

> plot(x,y)

It is based on two splines,

> lines(df$x,1+3*x2+5*x10)

If we use a spline model with 10 degrees of freedom, we get

> df=data.frame(x,y)
> reg=lm(y~bs(x,10),data=df)
> summary(reg)
 
Coefficients:
            Estimate Std. Er t value Pr(>|t|)    
(Intercept)  0.91671 0.17068   5.371 6.08e-07 ***
bs(x, 10)1   0.20485 0.32696   0.627    0.533    
bs(x, 10)2   3.15593 0.22534  14.005  < 2e-16 ***
bs(x, 10)3   0.04847 0.25075   0.193    0.847    
bs(x, 10)4   0.09373 0.21597   0.434    0.665    
bs(x, 10)5   0.11624 0.22939   0.507    0.614    
bs(x, 10)6   0.24829 0.22293   1.114    0.268    
bs(x, 10)7  -0.06825 0.23498  -0.290    0.772    
bs(x, 10)8   0.19633 0.26241   0.748    0.456    
bs(x, 10)9   0.27557 0.26976   1.022    0.310    
bs(x, 10)10  4.78134 0.24116  19.826  < 2e-16 ***

which makes sense, from what we have generated. Indeed, most of the components are not significant, but the second and the tenth. We can actually test that all those components are null (at the same time)

> A=matrix(0,8,11)
> colnames(A)=names(coefficients(reg))
> A[1,2]=A[2,4]=A[3,5]=A[4,6]=A[5,7]=
+ A[6,8]=A[7,9]=A[8,10]=1
> b=rep(0,8)
> linearHypothesis(reg, A,b)
Linear hypothesis test
 
Hypothesis:
bs(x, 10)1 = 0
bs(x, 10)3 = 0
bs(x, 10)4 = 0
bs(x, 10)5 = 0
bs(x, 10)6 = 0
bs(x, 10)7 = 0
bs(x, 10)8 = 0
bs(x, 10)9 = 0
 
Model 1: restricted model
Model 2: y ~ bs(x, 10)
 
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     98 4.8766                           
2     90 4.6196  8   0.25701 0.6259  0.754

and yes, those coefficients are not significant.

> yp10=predict(reg)
> lines(df$x,yp10,col="red")

Continue reading Regression with Splines: Should we care about Non-Significant Components?

How Could Classification Trees Be So Fast on Categorical Variables?

I think that over the past months, I have been saying non-correct things about classification with categorical covariates. Because I never took time to look at it carefuly. Consider some simulated dataset, with a logistic regression,

> n=1e3
> set.seed(1)
> X1=runif(n)
> q=quantile(X1,(0:26)/26)
> q[1]=0
> X2=cut(X1,q,labels=LETTERS[1:26])
> p=exp(-.1+qnorm(2*(abs(.5-X1))))/(1+exp(-.1+qnorm(2*(abs(.5-X1)))))
> Y=rbinom(n,size=1,p)
> df=data.frame(X1=X1,X2=X2,p=p,Y=Y)

Here, we use some continuous covariate, except that is considered as not-observed. Instead, we have a categorical covariate with 26 categories. The (theoretical) relationship between the covariate and the probability is given below,

> vx1=seq(0,1,by=.001)
> vp=exp(-.1+qnorm(2*(abs(.5-vx1))))/(1+exp(-.1+qnorm(2*(abs(.5-vx1)))))
> plot(vx1,vp,type="l")

and the empirical probability, for each modality is

If we run a classification tree, we get

> library(rpart)
> tree=rpart(Y~X2,data=df)
> library(rpart.plot)
> prp(tree, type=2, extra=1)

To be more specific, the output is here

> tree
1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.3 0.302
    4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *
    5) X2=F,G,H,I 153  37.22876 0.4183007       *
  3) X2=A,B,C,D,E,S,T,U,V,W,X,Y,Z 501 109.61 0.67
    6) X2=B,C,D,E,S,T,U,V,W,X 385  90.38 0.623  *
    7) X2=A,Y,Z 116  14.50862 0.8534483         *

 

Note that it takes less than a second to get that output. So clearly, we did not look for all combinations between modalities. For the first node, there are like  possible groups, i.e.

> 67108864

It is big… not huge, but too big to try all combinations, since that’s only the first node, and we have to do it again on the two leaves, etc. Antoine (aka @ly_antoine) told me – while we were having a coffee after lunch today – the trick to get a fast algorithm, on categories. And as usual, the idea is very clever…

First, we need a function to compute Gini index

> gini=function(y,classe){
+    T=table(y,classe)
+    nx=apply(T,2,sum)
+    n=sum(T)
+    pxy=T/matrix(rep(nx,each=2),nrow=2)
+    omega=matrix(rep(nx,each=2),nrow=2)/n
+    g=-sum(omega*pxy*(1-pxy))
+    return(g)}

For the first node, the idea is very simple:

  • Compute empirical averages 
> cond_prob=aggregate(df$Y,by=list(df$X2),mean)
  • Then sort those values, ,
  • Based on that ordering, consider 
> Group_Letters=cond_prob[order(cond_prob$x),2]

  • Then consider (only)  possible partitions,

against 

> v_gini=rep(NA,26)
> for(v in 1:26){
+   CLASSE=df$X2 %in% Group_Letters[1:v]
+   v_gini[v]=gini(y=df$Y,classe=CLASSE)
+ }

If we plot them, we get

> plot(1:26,v_gini,type="b)

As for continuous variables, we seek for the maximum value, and then, we have our two groups,

> sort(Group_Letters[1:which.max(v_gini)])
 [1] F G H I J K L M N O P Q R

That’s exactly what we got with the tree function in R,

1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30

Now, consider the leaf on the left (for instance)

> sub_df=df[df$X2 %in% sort(Group_Letters[1:which.max(v_gini)]),]

Then use the same algorithm as before: sort the conditional means,

> cond_prob=aggregate(sub_df$Y,by=
+ list(sub_df$X2),mean)
> s_Group_Letters=cond_prob[order(cond_prob$x),2]

Then compute Gini indices based on groups obtained from that ordering,

> v_gini=rep(NA,length(sub_Group_Letters))
> for(v in 1:length(sub_Group_Letters)){
+   CLASSE=sub_df$X2 %in% s_Group_Letters[1:v]
+   v_gini[v]=gini(y=sub_df$Y,classe=CLASSE)
+ }

If we plot it, we get our two groups,

> plot(1:length(s_Group_Letters),v_gini,type="b")

And the first group is here

> sort(sub_Group_Letters[1:which.max(v_gini)])
[1] J K L M N O P Q R

Again, that’s exactly what we got with the R function

1) root 1000 249.90000 0.4900000  
  2) X2=F,G,H,I,J,K,L,M,N,O,P,Q,R 499 105.30 0.30
    4) X2=J,K,L,M,N,O,P,Q,R 346  65.12 0.25144  *

Clever, isn’t?