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 ?

Evolution des Taux et Valeurs de Rentes

Dans mon billet, publié hier soir, sur Taux d’intérêt négatifs et explosion des valeurs de rentes, je montrais que les calculs de valeurs actuelles probables de rentes avec des taux de 5% (couramment utilisés il y a encore quelques années) ou -2% (les taux aujourd’hui sont faibles, voire négatifs) peut avoir un impact colossal. Mais qu’en est-il ‘pour de vrai’ ? Que se passe-t-il si j’actualise avec des taux ‘réellement’ observés, et pas des taux fixés arbitrairement.

Sur datamarket.com, on peut ainsi récupérer le taux à un an

Pourquoi un an? Je ne sais pas… Il me fallait en choisir un. C’est la plus courte maturité que j’ai pu trouver…. C’est un choix largement discutable… Il a un impact sur les valeurs numériques, probablement. Mais pas sur la tendance….

On commence par récupérer les données

> B=read.table("euro-yield-curves-daily-data.csv",sep=";",nrows=2925,header=TRUE)
> Y=as.numeric(as.character(B[,2]))
> Y[is.na(Y)]=as.numeric(paste("-",substr(as.character(B[is.na(Y),2]),2,nchar(as.character(B[is.na(Y),2]))),sep=""))
> D=as.Date(as.character(B[,1]),"%Y-%m-%d")
> YR=as.numeric(substr(as.character(B[,1]),1,4))
> plot(D,Y,type="l")
> abline(h=0,col="red")

Si on regarde uniquement la dernière année, depuis le 1er janvier 2015, les taux ont finalement très peu varié

en passant de 0.1% à -0.25% (soit une variation de l’ordre de 0.35 points, en 14 mois). Quel est l’impact sur les valeurs de rente ?

Reprenons le code d’hier (un peu modifié pour aller plus vite)

> file =paste("http://freakonometrics.free.fr/",
+ "HOM","-table-SPLx.csv",sep="")
> BH  = read.table(file,header=TRUE,sep=",")
> file =paste("http://freakonometrics.free.fr/",
+ "FEM","-table-SPLx.csv",sep="")
> BF  = read.table(file,header=TRUE,sep=",")
> PRIX=function(annee=2011,age,sexe="HOM",
+ taux=0.04, duree,C=1000){
+ if(sexe=="HOM") B=BH
+ if(sexe=="FEM") B=BF
+   an    = annee-age; if(an>2005){an=2005}
+   nom   = paste("X",an,sep="")
+   L     = B[,nom]
+   Q     = L[(age+1):length(L)]/L[(age+1)]
+   actualisation = (1+taux)^(0:min(duree,120-age))
+   prixsup = sum(Q[2:(min(duree,120-age)+1)]/
+   actualisation[2:(min(duree,120-age)+1)] )
+   prixinf = sum(Q[1:(min(duree,120-age))]/
+   actualisation[1:(min(duree,120-age))] )
+   return(C*c(prixsup,prixinf))}

On peut alors regarder l’évolution des rentes, avec les taux (et l’année, on va tenir compte aussi des gains d’espérance de vie)

Regardons pour un homme de 15 ans (l’idée était de regarder les rentes versées en cas d’accident corporel)

> F15=function(i) PRIX(annee=YEAR[i],age=15,sex="HOM",duree=150,taux=Y[i]/100)[1]
> rente15H=Vectorize(F15)(seq(1,length(Y),by=1))

ou pour une femme de 15 ans

> F15=function(i) PRIX(annee=YEAR[i],age=15,sex="FEM",duree=150,taux=Y[i]/100)[1]
> rente15F=Vectorize(F15)(seq(1,length(Y),by=1))

On peut aussi regarder un homme de 25 ans

> F25=function(i) PRIX(annee=YEAR[i],age=25,sex="HOM",duree=150,taux=Y[i]/100)[1]
> rente25H=Vectorize(F25)(seq(1,length(Y),by=1))

et une femme de 25 ans également

> F25=function(i) PRIX(annee=YEAR[i],age=25,sex="FEM",duree=150,taux=Y[i]/100)[1]
> rente25F=Vectorize(F25)(seq(1,length(Y),by=1))

Pour visualiser l’évolution des rentes, plaçons nous en base 100 en janvier 2008 (correspondant à l’époque où je rédigeais les annexes techniques)

> b15h=rente15H/rente15H[D==as.Date("2008-01-02","%Y-%m-%d")]
> b15f=rente15F/rente15F[D==as.Date("2008-01-02","%Y-%m-%d")]
> plot(D,b15f*100,col="red",type="l")
> lines(D,b15h*100,col="blue")

Aussi, une rente qui valait 100 en 2008 vaut aujourd’hui 350. Avec des taux passant de 3.5% à 0%. Et une croissance presque linéaire depuis 4 ans. Si on se limite aux 14 derniers mois (baisse de 0.35 points des taux), la valeur de la rente augmente, elle, de 10%

Si on regarde pour des jeunes de 25 ans (et plus 15 ans), les résultats sont relativement comparabless

> b25h=rente25H/rente25H[D==as.Date("2008-01-02","%Y-%m-%d")]
> b25f=rente25F/rente25F[D==as.Date("2008-01-02","%Y-%m-%d")]
> plot(D,b25f*100,col="red",type="l")
> lines(D,b25h*100,col="blue")

Une rente qui valait 100 en 2008 va aujourd’hui un peu moins de 350. Autrement dit, même en actualisant avec des taux de marché, on voit que la baisse des taux va avoir un impact très important sur les rentes. Même si la fréquence d’accidents graves baisse (ou disons reste stable), le coût des accidents corporels devrait continuer à augmenter dans les mois à venir… juste à cause des faibles taux d’intérêt (et de la valeur des rentes qui va exploser).

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)

Taux d’intérêt négatifs et explosion des valeurs de rentes

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrage évaluation du préjudice corporel, de Max le Roy. Comme l’indiquait le descriptif du livre (à l’époque) j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs,

“La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

En 2011, pour une mise à jour du livre, j’avais refait les calculs avec des tables mises à jours. Oui, à l’époque je trouvais que la priorité était d’arrêter d’utiliser les tables TD et TV qui dataient de 1988. Mais on était resté avec des taux à 5%. J’avais d’ailleurs été surpris (je peux le dire maintenant, il y a prescription) lorsque j’avais tenté de discuter le choix du taux d’intérêt, et qu’il m’avait été répondu “ah, ça change quelque chose ?“.

Oui, les taux d’actualisation impactent les valeurs des rentes. Et je n’avais pas idée à quel point. Le code utilisé pour construire les valeurs de rentes est assez simple.

> PRIX=function(annee=2011,age,sex="HOM",
+ taux=0.04, duree,C=1000){
+  file  = paste("http://freakonometrics.free.fr/",
+  sex,"-table-SPLx.csv",sep="")
+  B  = read.table(file,header=TRUE,sep=",")
+  an    = annee-age; if(an>2005){an=2005}
+  nom   = paste("X",an,sep="")
+  L     = B[,nom]
+  Q     = L[(age+1):length(L)]/L[(age+1)]
+  actualisation = (1+taux)^(0:min(duree,120-age))
+  prixsup = sum(Q[2:(min(duree,120-age)+1)]/
+  actualisation[2:(min(duree,120-age)+1)] )
+   prixinf = sum(Q[1:(min(duree,120-age))]/
+  actualisation[1:(min(duree,120-age))] )
+  return(C*c(prixsup,prixinf))}

Je renvois ici deux “prix” de rente, suivant que la rente est versée en fin ou en début d’année. Par exemple, avec un taux d’actualisation de 5% (comme dans le livre), pour un homme de 50 ans, si on considère une rente versée jusqu’à sa mort (j’ai mis ici une durée maximal de 150 années), la valeur actualisée espérée d’une rente (en tenant compte des probabilités de décès) de 1000 euros par an (encore une fois, tant que la personne est en vie) est de l’ordre de 17,000 euros.

> PRIX(age=50,sex="HOM",duree=150,taux=.05)
[1] 16435.31 17435.30

Avec un taux d’actualisation de 0%, on atteinte plus du double.

> PRIX(age=50,sex="HOM",duree=150,taux=0)
[1] 38822.48 39822.46

Et c’est encore pire pour les jeunes. Si on calcule pour plusieurs âges des rentes versées jusqu’au décès, avec différents niveaux de taux (5%, 2%, 0% et -2%, oui, je prends un taux négatif)

> vage=seq(15,65,by=5)
> vp5=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.05)[1])(vage)
> vp2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.02)[1])(vage)
> vp0=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0)[1])(vage)
> vpn2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=-0.02)[1])(vage)
> plot(vage,vpn2,type="b",col="red")
> lines(vage,vp0,type="b")
> lines(vage,vp2,type="b",col="blue")
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"))

on observe une explosion des rentes pour les jeunes (ce qui n’est pas non plus surprenant)

Pour un jeune de 20 ans, la valeur actualisée probable d’une rente est presque 10 fois plus élevée, si on passe d’un taux de 5% à -2%

> PRIX(age=15,sex="HOM",duree=150,taux=.05)
[1] 19427.94 20427.93
> PRIX(age=15,sex="HOM",duree=150,taux=-.02)
[1] 200778.3 201770.3

Ce qu’on peut visualiser sur le graphique ci-dessous (avec, par âge, 100 qui correspond à la valeur obtenue avec un taux de 5%)

> plot(vage,vpn2/vp5*100,type="b",col="red")
> lines(vage,vp0/vp5*100,type="b")
> lines(vage,vp2/vp5*100,type="b",col="blue")
> abline(h=100,lty=2)
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"),pch=1)

Passer de 5% à 2% donnait des rentes plus important, en passant de 100 à 200. En passant à 0% on passe de 100 à 400 pour les jeunes. Mais en passant à des taux négatifs, on passe de 100 à 1000 !

Dans un scénarios où les taux restent encore bas longtemps, on peut imaginer que les indemnités versées aux victimes de préjudices corporels vont exploser dans les mois à venir…

Spatial and Temporal Viz of Gas Price, in France

A great think in France, is that we can play with a great database with gas price, in all gas stations, almost eveyday. The file is rather big, so let’s make sure we have enough memory to run our codes,

> rm(list=ls())

To extract the data, first, we should extract the xml file, and then convert it in a more common R object (say a list)

> year=2014
> loc=paste("http://donnees.roulez-eco.fr/opendata/annee/",year,sep="")
> download.file(loc,destfile="oil.zip")

Content type 'application/zip' length 15248088 bytes (14.5 MB)

> unzip("oil.zip", exdir="./")
> fichier=paste("PrixCarburants_annuel_",year,
".xml",sep="")
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

We have a large dataset, with prices, for various types of gaz, for almost any gas station in France, almost every day, in 2014. It is a 1.4Gb list, with 11,064 elements (each of them being a gas station)

> length(l)
[1] 11064

There are two ways to look at the data. A first idea is to consider a gas station, and to extract the time series.

> time_series=function(no,type_gas="Gazole"){
+   prix=list()
+   date=list()
+   nom=list()
+   j=0
+   for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+         j=j+1
+  date[[j]]=as.character(l[[no]][[i]]["maj"])
+  prix[[j]]=as.character(l[[no]][[i]]["valeur"])
+  nom[[j]]=as.character(l[[no]][[i]]["nom"])
+       }}
+   }
+   id=which(unlist(nom)==type_gas)
+   n=length(id)
+   jour=function(j) as.Date(substr(date[[id[j]]],1,10),"%Y-%m-%d")
+   jour_heure=function(j) as.POSIXct(substr(date[[id[j]]],1,19), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   ext_y=function(j) substr(date[[id[j]]],1,4)
+   ext_m=function(j) substr(date[[id[j]]],6,7)
+   ext_d=function(j) substr(date[[id[j]]],9,10)
+   ext_h=function(j) substr(date[[id[j]]],12,13)
+   ext_mn=function(j) substr(date[[id[j]]],15,16)
+   prix_essence=function(i) as.numeric(prix[[id[i]]])/1000
+   base1=data.frame(indice=no,
+            id=l[[no]]$.attrs["id"],
+            adresse=l[[no]]$adresse,
+            ville=l[[no]]$ville,
+  lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+  lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+       cp=l[[no]]$.attrs["cp"],
+       saufjour=l[[no]]$ouverture["saufjour"], 
+       Y=unlist(lapply(1:n,ext_y)),
+       M=unlist(lapply(1:n,ext_m)),
+       D=unlist(lapply(1:n,ext_d)),
+       H=unlist(lapply(1:n,ext_h)),
+       MN=unlist(lapply(1:n,ext_mn)),
+    prix=unlist(lapply(1:n,prix_essence)))
+   
+   base1=base1[!is.na(base1$prix),]
+   
+   date_d=paste(year,"-01-01 12:00:00",sep="")
+   date_f=paste(year,"-12-31 12:00:00",sep="")
+   vecteur_date=seq(as.POSIXct(date_d, format =
+                 "%Y-%m-%d %H:%M:%S"),
+                    as.POSIXct(date_f, format = 
+                 "%Y-%m-%d %H:%M:%S"),by="days")
+   date=paste(base1$Y,"-",base1$M,"-",base1$D,
+   " ",base1$H,":",base1$MN,":00",sep="")
+   date_base=as.POSIXct(date, format = 
+                "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   idx=function(t) sum(vecteur_date[t]>=date_base)
+   vect_idx=Vectorize(idx)(1:length(vecteur_date))
+   P=c(NA,base1$prix)
+   base2=ts(P[1+vect_idx],
+         start=year,frequency=365)
+   list(base=base1,
+        ts=base2)
+ }

To get the time series, extrapolation is necessary, since we have here observation at irregular dates. Here, for instance, for the second gas station, we get

> plot(time_series(2)$ts,ylim=c(1,1.6),col="red")
> lines(time_series(2,"SP98")$ts,col="blue")

An alternative is to study gas price from a spatial perspective. Given a date, we want the price in all stations. As previously, we keep the last price observed, in each station,

> spatial=function(dt){
+   base=NULL
+   for(no in 1:length(l)){  
+     prix=list()
+     date=list()
+     j=0
+     for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+   j=j+1
+   date[[j]]=as.character(l[[no]][[i]]["maj"])
+       }}
+   }
+   n=j
+   D=as.Date(substr(unlist(date),1,10),"%Y-%m-%d")
+   k=which(D==D[which.max(D[D<=dt])])
+ if(length(k)>0){
+   B=Vectorize(function(i) l[[no]][[k[i]]])(1:length(k))
+ if("nom" %in%  rownames(B)){  
+   k=which(B["nom",]=="Gazole")
+   prix=as.numeric(B["valeur",k])/1000
+   if(length(prix)==0) prix=NA
+   base1=data.frame(indice=no,
+   lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+   lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+   gaz=prix)
+   base=rbind(base,base1)
+ }}}
+ return(base)}

For instance, for the 5th of May, 2014, we get the following dataset

> B=spatial(as.Date("2014-05-05"))

To visualize prices, consider only mainland France (excluding islands in the Pacific, or close to the Caribeans)

> idx=which((B$lon>(-10))&(B$lon<20)&
+ (B$lat>35)&(B$lat<55))
> B=B[idx,]
> Q=quantile(B$gaz,seq(0,1,by=.01),na.rm=TRUE)
> Q[1]=0
> x=as.numeric(cut(B$gaz,breaks=unique(Q)))
> CL=c(rgb(0,0,1,seq(1,0,by=-.025)),
+ rgb(1,0,0,seq(0,1,by=.025)))
> plot(B$lon,B$lat,pch=19,col=CL[x])

Red dots are the most expensive gas stations, that particular day.

If we add contours of the French regions, we get

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=CL[x])

 

We can also focus on some specific region, say the South of Brittany.

> library(OpenStreetMap)
> map <- openmap(c(lat= 48,   lon= -3),
+                c(lat= 47,   lon= -2))
> map <- openproj(map) 
> plot(map)
> points(B$lon,B$lat,pch=19,col=CL[x])

As we can see on that map, there are regions that are rather empty, where the closest gas station might be a bit far away. Actually, it is possible to add Voronoi sets on the map,

> dB=data.frame(lon=B$lon,lat=B$lat)
> idx=which(!duplicated(dB))
> dB=dB[idx,]

 

which could help to get the price of the closest gaz station.

> library(tripack)
> V <- voronoi.mosaic(dB$lon[id],dB$lat[id])
> plot(V,add=TRUE)

It is possible to plot each polygon with the color of the gaz station we add. Actually, it is a bit tricky, and I could not find a R function to to this. So I did it manually,

> plot(map)
> P <- voronoi.polygons(V)
> library(sp)
> point_in_i=function(i,point) point.in.polygon(point[1],point[2],P[[i]][,1],P[[i]][,2])
> which_point=function(i) which(Vectorize(function(j) point_in_i(i,c(dB$lon[id[j]],dB$lat[id[j]])))(1:length(id))>0)
> for(i in 1:length(P)) polygon(P[[i]],col=CL[x[id[which_point(i)]]],border=NA)

With this map, we can see that we have blue areas, i.e. all stations in a given area are cheap (because of competition), but in some places, a very expensive one is next to a very cheap one. I guess we should look closer at the dynamics… [to be continued….]

Le principe de mutualisation est remis en cause par les données

Interview sur le site de l’Institut Louis Bachelier.

En quoi la massification des données bouleverse les principes de l’assurance ?

La base de l’assurance est de regrouper des personnes afin de constituer une mutualité. Le prix des primes est calculé de sorte à couvrir le coût moyen des sinistres de cette mutualité.

L’arrivée de nouvelles données commence à remettre en cause ce principe de mutualisation, au profit d’assurances à la carte. La segmentation en fonction du profil de l’assuré est de plus en plus fine.

Différencier les produits et tarifs en fonction des caractéristiques du client, comme son âge ou son lieu d’habitation, n’est pas nouveau…

La segmentation a toujours existé. Un jeune conducteur, par exemple, paie davantage qu’un conducteur plus expérimenté. Mais elle peut s’effectuer aujourd’hui à un niveau beaucoup plus fin grâce aux données GPS ou celles des objets connectés notamment. La question est donc de savoir jusqu’où il est souhaitable d’aller ? Quel est le bon équilibre entre segmentation et mutualisation ?

Si on exclut l’idée de faire du bénéfice (ce qui n’est pas la vocation des mutuelles), l’assurance est « un «jeu à somme nulle ». Autrement dit, si les tarifs de certains clients baissent, d’autres augmenteront nécessairement, au risque d’exclure certaines personnes du marché, comme cela a pu être le cas pour les jeunes conducteurs en Irlande, avec des tarifs proposés très élevés.

Il s’agit donc d’une vraie question de société sur laquelle la recherche peut apporter des éléments de réponse.

Les assureurs ne seront-ils pas tentés de choisir uniquement les profils les moins risqués ?

Certains assureurs peuvent être tentés de chercher des niches, des segments peu risqués, et plus profitables. Mais la recherche de niche est un jeu dangereux. Il devient de plus en plus dur de distinguer le signal du bruit, et de s’assurer que la niche perdurera dans le temps. Il faut aussi tenir compte de biais de sélection dans la base.

Quels sont les autres risques d’une sur segmentation ?

En s’appuyant sur une plus large base d’informations, des compagnies vont mener des stratégies tarifaires assez agressives. Cela risque de générer une forte volatilité des primes d’assurance et d’entrainer d’importants mouvements chez les assurés.

Encourager la concurrence entre les compagnies d’assurance n’est-il pas souhaitable?

D’importants mouvements des assurés peuvent en tous cas soulever certaines problématiques. Aujourd’hui, les assureurs gardent souvent leurs clients plusieurs années et peuvent donc mutualiser les risques et les coûts associés dans le temps. Les jeunes conducteurs ont toujours plus d’accidents que la moyenne, mais avec les années d’expérience, leurs risques diminuent. Une compagnie peut donc accepter d’assurer un nouveau conducteur, à un tarif raisonnable, en misant sur le fait qu’elle le conservera comme client les années suivantes. Cette mutualisation dans le temps ne sera plus possible si les assurés deviennent plus volatiles.

Leurs comportement sont encore mal connus et compris. Au-delà de la tarification et de la création de produits, l’arrivée de nouvelles données peut aider à mieux les comprendre: pourquoi choisit-on un contrat plutôt qu’un autre ? Dans le cas des couvertures non-obligatoires, pourquoi décidons-nous de nous assurer ou pas ? Les données permettent de s’intéresser à l’assurance en tant qu’activité économique.

Les assureurs s’appuient sur des modèles mathématiques afin de prédire la sinistralité de leur portefeuille. L’ajout de nouvelles données peut-il permettre d’améliorer ces prédictions ?

Les algorithmes permettent déjà de prédire de façon très fine la probabilité qu’un sinistre survienne dans l’année. Les nouvelles données ne vont rien changer sur ce point tant que le cycle de production ne sera pas inversé. Par contre, elles seront certainement utiles sur le plan de la prévention. Grâce aux objets connectés, ou encore aux voitures intelligentes, il est désormais possible de détecter les comportements à risque. Des actions de prévention adéquates pourraient alors être mise en place.

Central Limit Theorem

The central limit theorem is a fundamental theorem of statistics. It prescribes that the sum of a sufficiently large number of independent and identically distributed random variables approximately follows a normal distribution.

History of the Central Limit Theorem

The term “central limit theorem” most likely traces back to Georg Pólya. As he recapitulated at the beginning of a paper published in 1920, it was “generally known that the appearance of the Gaussian probability density ” in a great many situations “can be explained by one and the same limit theorem” which plays “a central role in probability theory”. Laplace had discovered the essentials of this fundamental theorem in 1810, and with the designation “central limit theorem of probability theory” which was even emphasized in the paper’s title, Pólya gave it the name that has been in general use ever since.

In this paper of 1820, Laplace starts by proving the central limit theorem for some certain probability distributions. He then continues with arbitrary discrete and continuous distributions. But a more general (and rigorous) proof should be attributed to Siméon Denis Poisson. He also intuited that weaker version could easily be derived. As for Laplace, the main purpose of that Central Limit Theorem for Poisson was to be a tool in calculations, not so much to be a mathematical theorem in itself. Therefore, neither Laplace nor Poisson explicitly formulate any conditions for the theorem to hold. The mathematical formulation of the theorem is due to the St. Petersburg School of probability, from 1870 until 1910, with Chebyshev, Markov and Liapounov.

Mathematical Formulation

Let be independent random variables that are identically distributed, with mean  and finite variance. Let

then from the law of large numbers  tend to 0 as  tends to infinity. The central limit theorem establishes that the distribution of  tends to a centered normal distribution when  goes to infinity. More specificaly,

We can also write

or .

A limiting result as an approximation

This central limit thereom is used to approximate distributions derived from summing, or averaging, identical random variables.

Consider for instance a course where 7 students out of 8 pass. What is the probability that (at least) 4 failed in a class of 25 students. Let be the dichotomous variable that describe failure : 1 if the student failed and 0 if he passed. That random variable has a Bernoulli distribution with parameter=1/8, with mean1/8 and variance 7/64. Consequently, if students’ grades are independent, the sum  follows a binomial distribution, with mean  and variance , which can be approximated, by the central limit theorem, by a normal distribution with mean  and variance . Here, =3.125 while =2.734. To compute  either either the binomial distribution, or the Gaussian approximation. In the first case, the probability is 80.47 %,

In the second case, use a continuity correction, and compute the probability that  is less than 4+1/2. From the central limit theorem

The probability that a standard Gaussian variable is less than this quantity is

which can be compared with 80.47% obtained without the approximation, see Figure 1. Note that this approximation was obtained by De Moivre, in 1713, and is usually known as « Bernoulli’s law of large numbers ».

Figure 1: Gaussian approximation of the binomial distribution.

Asymtptotic Confidence Intervals

The intuition is that a confidence interval is an interval in which one may be confident that a parameterof interest lies. For instance, that some quantity is measured , but the measurement is subject to a normally distributed error, with known variance . If  has a distribution, we know that

Equivalently, we could write

or 

Thus, if   is measured to be , then the 95 % confidence interval for  is .

In the context of Bernoulli trials (described above), the asyymptotic 95 % confidence interval for p is

A popular rule of thumb can be derived when ~50%. In that context is close to 1.96 (or 2), and a 95 % approximated confidence interval is then

see Figure 2. If that confidence interval provides a good approximation for the 95 % confidence interval when ~50 %, it is an over-estimation when p is either much smaller, or much larger.

Figure 2: law of large numbers on the left, with the convergence of towards  as  increases, and central limit theorem, on the right, with the convergence of  towards a Gaussian distribution. The red area is the 95% confidence region.

The Delta Method and Method of Moments

This method is used to approximate a general transformation of a parameter that is known to be asymptotically normal,

then

Consider now a parametric model, … independent, with identical distribution (which can be a Weibull distribution to model a duration, a Pareto distribution to model the income or the wealth, etc). The method of moments is a method of estimating parameters based on equating population and sample values of certain moments of the distribution. For instance, if , then the estimator of the unknown parameter is given by equation or equivalently . From the central limit theorem

and applying the delta-method with , then

where a numerical approximation for the variance can be derived.This method has a long history, and has been intensively studied. Furthermore, this asymptotic normality can be used to compute a confidence interval, and also to derive an asymptotic testing procedure.

An Asymptotic Testing Procedure

Based on that asymptotic normality, it is possible to derive a simple testing procedure. Consider a test of the hypothesis  against , usually called a “significant” test for parameter  (or significance of an explanatory variance in the context of regression model). Under the assumption that  is valid, then

for some variance , that can be computed using the delta method. The p-value associated with that test is

where  is the observed empirical estimator of the parameter and  is a standard normal variable. Thus, the p-value can easily be computed using quantiles of the standard normal distribution. Here, the p-value is above 5% if

Weaker Forms of the Central Limit Theorem

As stated by Laplace, the Central Limit Theorem relies on strong assumption. Hopefully, most of them can be relaxed. In a first variant of the theorem, random variables have to be independent, but not necessarily identically distributed. If random variables  have averages  and , then  and  in the Central Limit Theorem should be replaced by averages of  and ‘s, with an additional technical assumption related to the existence of some higher moments (the so called Lyapounov condition).

For a second variant of the theorem, random variables can be dependent, as in ergodic Markov chain, or in autoregressive time series. In that context, if  is a stationnary time series, with mean \mu, then define

and with that limit, the central limit theorem hold

even if the variance term has here a different interpretation.

Finally, a third variant that can be mentioned is the one obtained be Paul Lévy about asymptotic properties of the empirical average, when the variance is not finite (actually, even when the first moment in not finite). In that case, the limiting distribution is no longer Gaussian.

References

Laplace, P.S. de (1810). Mémoire sur les approximations des formules qui sont fonctions de très grands nombres et sur leur application aux probabilités. Mémoires de l’Académie Royale des Sciences de Paris, 10.

Le Cam, L. (1986), The Central Limit Theorem around 1935, Statistical Science 1(1): 78-96

Polyà, G. (1920). Ueber den zentralen Grenzwertsatz der Wahrscheinlichkeitsrechnung und das Momentproblem. Mathematische Zeitschrift 8, 171–181.

Reverse Engineering with Correlated Features

In econometric modeling, I usually have a problem with correlated features. A few weeks ago, I was discussing feature selection when features are correlated. This week, I was wondering about reverse engineering when features might be correlated (not to say very correlated). The way I see reverse engineering is the following

  1. someone has some dataset, and based on that dataset, a model was fitted. But we cannot see how it works….
  2. we can generate “fake data”, feed the model with those data, and get predictions
  3. based on those predictions, we wish we can fit a model that should be close to the the ‘true’ model used
  4. one way to measure how good our model is is to compare predictions on the initial data with our model with the original dataset (or the initial ‘true’ values if we use generated datasets).

Continue reading Reverse Engineering with Correlated Features

Clustering French Cities (based on Temperatures)

In order to illustrate hierarchical clustering techniques and k-means, I did borrow François Husson‘s dataset, with monthly average temperature in several French cities.

> temp=read.table(
+ "http://freakonometrics.free.fr/FR_temp.txt",
+ header=TRUE,dec=",")

We have 15 cities, with monthly observations

> X=temp[,1:12]
> boxplot(X)

Since the variance seems to be rather stable, we will not ‘normalize’ the variables here,

> apply(X,2,sd)
    Janv     Fevr     Mars     Avri 
2.007296 1.868409 1.529083 1.414820 
     Mai     Juin     juil     Aout 
1.504596 1.793507 2.128939 2.011988 
    Sept     Octo     Nove     Dece 
1.848114 1.829988 1.803753 1.958449

In order to get a hierarchical cluster analysis, use for instance

> h <- hclust(dist(X), method = "ward")
> plot(h, labels = rownames(X), sub = "")

An alternative is to use

> library(FactoMineR)
> h2=HCPC(X)
> plot(h2)

Here, we visualise observations with a principal components analysis. We have here also an automatic selection of the number of classes, here 3. We can get the description of the groups using

> h2$desc.ind

or directly

> cah=hclust(dist(X))
> groups.3 <- cutree(cah,3)

We can also visualise those classes by ourselves,

> acp=PCA(X,scale.unit=FALSE)
> plot(acp$ind$coord[,1:2],col="white")
> text(acp$ind$coord[,1],acp$ind$coord[,2],
+ rownames(acp$ind$coord),col=groups.3)

It is possible to plot the centroïds of those clusters

> PT=aggregate(acp$ind$coord,list(groups.3),mean)
> points(PT$Dim.1,PT$Dim.2,pch=19)

If we add Voroid sets around those centroïds, here we do not see them (actually, we see the point – in the middle – that is exactly at the intersection of the three regions),

> library(tripack)
> V <- voronoi.mosaic(PT$Dim.1,PT$Dim.2)
> plot(V,add=TRUE)

To visualize those regions, use

> p=function(x,y){
+   which.min((PT$Dim.1-x)^2+(PT$Dim.2-y)^2)
+ }
> vx=seq(-10,12,length=251)
> vy=seq(-6,8,length=251)
> z=outer(vx,vy,Vectorize(p))
> image(vx,vy,z,col=c(rgb(1,0,0,.2),
+ rgb(0,1,0,.2),rgb(0,0,1,.2)))
> CL=c("red","black","blue")
> text(acp$ind$coord[,1],acp$ind$coord[,2],
+ rownames(acp$ind$coord),col=CL[groups.3])

Actually, those three groups (and those three regions) are also the ones we obtain using a k-mean algorithm,

> km=kmeans(acp$ind$coord[,1:2],3)
> km
K-means clustering 
with 3 clusters of sizes 3, 7, 5

(etc). But actually, since again we have some spatial data, it is possible to visualize them on a map

> library(maps)
> map("france")
> points(temp$Long,temp$Lati,col=groups.3,pch=19)

or, to visualize the regions, use e.g.

> library(car)
> for(i in 1:3) 
+ dataEllipse(temp$Long[groups.3==i],
+ temp$Lati[groups.3==i], levels=.7,add=TRUE,
+ col=i+1,fill=TRUE)

Those three regions actually make sense, geographically speaking.

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/…

Clusters of Texts

Another popular application of classification techniques is on texmining (see e.g. an old post on French president speaches). Consider the following example,  inspired by Nobert Ryciak’s post, with 12 wikipedia pages, on various topics,

> library(tm)
> library(stringi)
> library(proxy)
> titles = c("Boosting_(machine_learning)",
+            "Random_forest",
+            "K-nearest_neighbors_algorithm",
+            "Logistic_regression",
+            "Boston_Bruins",
+            "Los_Angeles_Lakers",
+            "Game_of_Thrones",
+            "House_of_Cards_(U.S._TV_series)",
+            "True Detective (TV series)",
+            "Picasso",
+            "Henri_Matisse",
+            "Jackson_Pollock")
> articles = character(length(titles))
> wiki = "http://en.wikipedia.org/wiki/"
> for (i in 1:length(titles)) {
+   articles[i] = stri_flatten(readLines(stri_paste(wiki, titles[i])), col = " ")
+ }

Here, we store all the contents of the pages in a corpus (from the text mining package).

> docs = Corpus(VectorSource(articles))

This is what we have in that corpus

> a = stri_flatten(readLines(stri_paste(wiki, titles[1])), col = " ")
> a = Corpus(VectorSource(a))
> a[[1]]

Thoughts on Hypothesis Boosting</i></a>, Unpublished manuscript (Machine Learning class project, December 1988)</span></li> <li id="cite_note-4"><span class="mw-cite-backlink"><b><a href="#cite_ref-4">^</a></b></span> <span class="reference-text"><cite class="citation journal"><a href="/wiki/Michael_Kearns" title="Michael Kearns">Michael Kearns</a>; <a href="/wiki/Leslie_Valiant" title="Leslie Valiant">Leslie Valiant</a> (1989). <a rel="nofollow" class="external text" href="http://dl.acm.org/citation.cfm?id=73049">"Crytographic limitations on learning Boolean formulae and finite automata"</a>. <i>Symposium on T

This is because we read an html page.

> a = tm_map(a, function(x) 
> a = tm_map(a, function(x) stri_replace_all_fixed(x, "\t", " "))
> a = tm_map(a, PlainTextDocument)
> a = tm_map(a, stripWhitespace)
> a = tm_map(a, removeWords, stopwords("english"))
> a = tm_map(a, removePunctuation)
> a = tm_map(a, tolower)
> a 

can  set  weak learners create  single strong learner  a weak learner  defined    classifier    slightly correlated   true classification  can label examples better  random guessing in contrast  strong learner   classifier   arbitrarily wellcorrelated   true classification robert 

Now we have the text of the wikipedia document. What we did was

  • replace all “” elements with a space. We do it because there are not a part of text document but in general a html code.
  • replace all “/t” with a space.
  • convert previous result (returned type was “string”) to “PlainTextDocument”, so that we can apply the other functions from tm package, which require this type of argument.
  • remove extra whitespaces from the documents.
  • remove punctuation marks.
  • remove from the documents words which we find redundant for text mining (e.g. pronouns, conjunctions). We set this words as stopwords(“english”) which is a built-in list for English language (this argument is passed to the function removeWords.
  • transform characters to lower case.

Now we can do it on the entire corpus

> docs2 = tm_map(docs, function(x) stri_replace_all_regex(x, "<.+?>", " "))
> docs3 = tm_map(docs2, function(x) stri_replace_all_fixed(x, "\t", " "))
> docs4 = tm_map(docs3, PlainTextDocument)
> docs5 = tm_map(docs4, stripWhitespace)
> docs6 = tm_map(docs5, removeWords, stopwords("english"))
> docs7 = tm_map(docs6, removePunctuation)
> docs8 = tm_map(docs7, tolower)

Now, we simply count words in each page,

> dtm <- DocumentTermMatrix(docs8)
> dtm2 <- as.matrix(dtm)
> dim(dtm2)
[1] 12 13683
> frequency <- colSums(dtm2)
> frequency <- sort(frequency, decreasing=TRUE)
> mots=frequency[frequency>20]
> s=dtm2[1,which(colnames(dtm2) %in% names(mots))]
> for(i in 2:nrow(dtm2)) s=cbind(s,dtm2[i,which(colnames(dtm2) %in% names(mots))])
> colnames(s)=titles

 

Once we have that dataset, we can use a PCA to visualise the ‘variables’ i.e. the pages

> library(FactoMineR)
> PCA(s)

We can also use non-supervised classification to group pages. But first, let us normalize the dataset

> s0=s/apply(s,1,sd)

Then, we can run a cluster dendrogram, using the Ward distance

> h <- hclust(dist(t(s0)), method = "ward")
> plot(h, labels = titles, sub = "")

Groups are consistent with intuition: painters are in the same cluster, as well as TV series, sports teams, and statistical techniques.

Clusters of (French) Regions

For the data scienec course of tomorrow, I just wanted to post some functions to illustrate cluster analysis. Consider the dataset of the French 2012 elections

> elections2012=read.table(
"http://freakonometrics.free.fr/elections_2012_T1.csv",sep=";",dec=",",header=TRUE)
> voix=which(substr(names(
+ elections2012),1,11)=="X..Voix.Exp")
> elections2012=elections2012[1:96,]
> X=as.matrix(elections2012[,voix])
> colnames(X)=c("JOLY","LE PEN","SARKOZY","MÉLENCHON","POUTOU","ARTHAUD","CHEMINADE","BAYROU","DUPONT-AIGNAN","HOLLANDE")
> rownames(X)=elections2012[,1]

The hierarchical cluster analysis is obtained using

> cah=hclust(dist(X))
> plot(cah,cex=.6)

To get five groups, we have to prune the tree

> rect.hclust(cah,k=5)
> groups.5 <- cutree(cah,5)

We have to zoom-in to visualize the French regions,

It is also possible to use

> library(dendroextras)
> plot(colour_clusters(cah,k=5))

And again, if we zoom-in, we get

The interpretation of the clusters can be obtained using

> aggregate(X,list(groups.5),mean)
  Group.1     JOLY   LE PEN  SARKOZY
1       1 2.185000 18.00042 28.74042
2       2 1.943824 23.22324 25.78029
3       3 2.240667 15.34267 23.45933
4       4 2.620000 21.90600 34.32200
5       5 3.140000  9.05000 33.80000

It is also possible to visualize those clusters on a map, using

> library(RColorBrewer)
> CL=brewer.pal(8,"Set3")
> carte_classe <- function(groupes){
+ library(stringr)
+ elections2012$dep <- elections2012[,2]
+ elections2012$dep <- tolower(elections2012$dep)
+ elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
+ library(maps)
+ france<-map(database="france")
+ france$dep <- france$names
+ france$dep <- tolower(france$dep)
+ france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
+ corresp_noms <- elections2012[, c(1,2, ncol(elections2012))]
+ corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
+ col2001<-groupes+1
+ names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms[,1])]
+ color <- col2001[match(france$dep, names(col2001))]
+ map(database="france", fill=TRUE, col=CL[color])
+ }
> carte_classe(cutree(cah,5))

or, if we simply want 4 clusters

> carte_classe(cutree(cah,4))

 

Chaire de Recherche, actinfo

Je n’avais pas trop eu l’occasion d’en parler, mais avec Romuald Elie, professeur à l’université Paris-Est, on porte – depuis la fin 2015 – une chaire de recherche financée par Covéa (réunissant notamment les marques GMF, MAAF et MMA), sur les usages de l’information en assurance, sur la compréhension des mécanismes de marché d’un marché de plus en plus segmenté et des contradictions avec l’hypothèse de mutualisation des risques. Cette chaire, valorisation et nouveau usages actuariel de l’information – le petit nom est actinfo – est soutenue par le GENES, l’université Paris Est et l’université de Rennes 1. Un blog dédié permettra de disséminer la recherche qui sera menée dans le cadre de cette chaire, actinfo.hypotheses.org. A suivre….

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.