# Forecasts with ARIMA Models

In our time series class this morning, I was discussing forecasts with ARIMA Models. Consider some simple stationnary AR(1) simulated time series

> n=95
> set.seed(1)
> E=rnorm(n)
> X=rep(0,n)
> phi=.85
> for(t in 2:n) X[t]=phi*X[t-1]+E[t]
> plot(X,type="l")

If we fit an AR(1) model,

> model=arima(X,order=c(1,0,0),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*P$se[20],lty=2,col="red")
> abline(h=-2*P$se[20],lty=2,col="red") we observe the exponential decay of the forecast towards 0, and the increasing confidence interval (where the variance increases, from the variance of the white noise to the variance of the stationnary time series). Plain lines are conditional forecast (given our latest observation, since the AR(1) is a first order Markov process), and dotted lines are unconditional. Let us store some values, to use them as benchmark > s=P$se[20]
> y=P$pred If we fit a MA(1) model > model=arima(X,order=c(0,0,1), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

after two lags, the forecast is null, and the (conditional) variance remains constant. But if we consider a moving average process with a longer order,

> model=arima(X,order=c(0,0,14),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*s,lty=2,col="red") > abline(h=-2*s,lty=2,col="red") > lines(y,col="grey") we get an output that can be compared with the AR(1) processes. Which makes sense since our AR(1) process can also be seen as a MA(∞), with infinite order. But if we think that our time series is not stationary, an we fit an integrated model > model=arima(X,order=c(0,1,0), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

we observe the (standard) martingale property: the forecast is flat, and the confidence interval keeps increasing, and actually, the variance increases towards infinity (at a linear rate). So one should be very careful when differentiation a time series… it will have a huge impact on the forecasts….

# Where People Live

There was an interesting map on reddit this morning, with a visualisation of latitude and longituge of where people live, on Earth. So I tried to reproduce it. To compute the density, I used a kernel based approch

> library(maps)
> data("world.cities")
> X=world.cities[,c("lat","pop")]
> liss=function(x,h){
+   w=dnorm(x-X[,"lat"],0,h)
+   sum(X[,"pop"]*w)
+ }
> vx=seq(-80,80)
> vy=Vectorize(function(x) liss(x,1))(vx)
> vy=vy/max(vy)
> plot(world.cities$lon,world.cities$lat,)
> for(i in 1:length(vx))
+ abline(h=vx[i],col=rgb(1,0,0,vy[i]),lwd=2.7)

For the other axis, we use a miror technique, to ensure that -180 is close the +180

> Y=world.cities[,c("long","pop")]
> Ya=Y; Ya[,1]=Y[,1]-360
> Yb=Y; Yb[,1]=Y[,1]+360
> Y=rbind(Y,Ya,Yb)
> liss=function(y,h){
+   w=dnorm(y-Y[,"long"],0,h)
+   sum(Y[,"pop"]*w)
+ }
> vx=seq(-180,180)
> vy=Vectorize(function(x) liss(x,1))(vx)
> vy=vy/max(vy)
> plot(world.cities$lon,world.cities$lat,pch=19)
> for(i in 1:length(vx))
+ abline(v=vx[i],col=rgb(1,0,0,vy[i]),lwd=2.7)

Now we can add the two, on the same graph

# R Crash Course, Data Science for Actuaries, Year 2

This Monday, we will start the second year of the Actuary: Data Science (ADS) program, supported by the (French) Institute of Actuaries. I will be there on monday morning for the opening, and we will start the R & Datamining course. The slides are now online,

In order to get nice slides, I have been using slidify. The R code is available online, as well as some pdf version of the slides.

# 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)

# 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="")

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

# 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).

# 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. # 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))

# Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates$X$$Y|X$ has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that $Y$ has a similar distribution… so there is no reason to test for a Gamma model for $Y$ 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 $(X,Y)$ has a joint Gaussien distribution, then both marginals are Gaussian, but also $Y|X$. So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for $Y$. The econometric interpretation is that with a standard Gaussian linear model, if$X$ is normally distributed, not only the conditional distribution $Y|X$ is Gaussian but also the non-conditional distribution of $Y$.

> 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 $Y$ is also Gaussian

> library(nortest)

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 $X$ is continuous. What if we consider a finite mixture here, i.e.$X$ takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for $Y$. But in practice, would we really reject the Gaussian assumption, for $Y$? 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)
>  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") # 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") # 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$X_1$, except that$X_1$ 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 $2^{26}$ 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 $p_i=\mathbb{E}[Y \vert X=\{x_i\}]$ > cond_prob=aggregate(df$Y,by=list(df$X2),mean) • Then sort those values, $p_{1:k}\leq p_{2:k} \leq \cdots \leq p_{k:k}$, • Based on that ordering, consider $\{x_{1:k},x_{2:k},\cdots,x_{k:k}\}$ > Group_Letters=cond_prob[order(cond_prob$x),2]

• Then consider (only) $k$ possible partitions,

$\{x_{1:k},x_{2:k},\cdots,x_{i:k}\}$ against $\{x_{i+1:k},x_{i+2:k},\cdots,x_{k:k}\}$

> 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?

# Inter-relationships in a matrix

Last week, I wanted to displaying inter-relationships between data in a matrix. My friend Fleur, from AXA, mentioned an interesting possible application, in car accidents. In car against car accidents, it might be interesting to see which parts of the cars were involved. On https://www.data.gouv.fr/fr/, we can find such a dataset, with a lot of information of car accident involving bodily injuries (in France, a police report is necessary, and all of them are reported in a big dataset… actually several dataset, with information of people involved, cars, locations, etc). For 2014 claims, the dataset is

> base = read.csv("https://www.data.gouv.fr/s/resources/base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-153355/vehicules_2014.csv")


Let us keep only claims involving two vehicules,

> T=table(base$Num_Acc) > idx=names(T)[which(T==2)] For 2014, we have 32,222 claims. > length(idx) [1] 32222 In this dataset, we have information about where cars were hit, plus ‘9’ for multiple hot (in rollover accidents) and ‘0’ should be missing information. > nom=c("NA","Front","Front R",'Front L',"Back","Back R","Back L","Side R","Side L","Multiple")  Now, we simply have to go through our dataset, and get the matrix. My first idea was to get a symmetric one, > B=base[base$Num_Acc %in% idx,]
> B=B[order(B$Num_Acc),] > M=matrix(0,10,10) > for(i in seq(1,nrow(B),by=2)){ + a=B$choc[i]+1
+   b=B$choc[i+1]+1 + M[a,b]=M[a,b]+1 + M[b,a]=M[b,a]+1 + } > rownames(M)=nom > colnames(M)=nom The problem, when we ask for a symmetric chord diagram, is that we cannot have Front – Front claims (since values on the diagonal are removed) > library(circlize) > chordDiagramFromMatrix(M,symmetric=TRUE) So let’s pretend that there could be some possible distinction in the dataset, between the first and the second row. Like the first one is the ‘responsible’ driver. Or like, for insurer, the first one is your insured. Just to avoid this symmetry problem > M=matrix(0,10,10) > for(i in seq(1,nrow(B),by=2)){ + a=B$choc[i]+1
+   b=B\$choc[i+1]+1
+ M[a,b]=M[a,b]+1
+ }
> rownames(M)=paste("A",nom,sep=" ")
> colnames(M)=paste("B",nom,sep=" ")

If we visualize the chord diagram, this time it is more complex to analyze,

> chordDiagram(M)

Below we have the first row (say our driver, letter A) and on top, the second row (say the other driver, letter B),

In bodily injury claims, we observe a large proportion of Front – Front claims, as well as Front – Back. And as expected Back-Back are not that common….

A few month ago, I did mention a graph, of some so-called Lorenz curves to compare regression models, see e.g. Progressive’s slides (thanks Guillaume for the reference)

The idea is simple. Consider some model for the pure premium (in insurance, it is the quantity that we like to model), i.e. the conditional expected valeur

$m(\boldsymbol{x})=\mathbb{E}[Y\vert \boldsymbol{X}=\boldsymbol{x}]$

On some dataset, we have our predictions,$\widehat{y}_i=m(\boldsymbol{x}_i)$ as well as observed quantities, $y_i$. The curve are obtained simply :

• sort the observations so that

$\widehat{y}_1\geq \widehat{y}_2 \geq \cdots \geq \widehat{y}_i \geq \cdots\geq \widehat{y}_n$

• based on that ordering (from high risks to low risks, based on our predictions), we plot Lorenz curve

$\left(\frac{i}{n},\frac{\displaystyle\sum_{k=1}^i y_k}{\displaystyle\sum_{k=1}^n y_k}\right)$