# Growing one Tree

Consider the following toy dataset, with some spam/ham information, and two words, “viagra” and “lottery”.

> load(spam.RData)
Y viagra lottery
27 spam      0       1
37  ham      0       1
57 spam      0       0
89  ham      0       0
20 spam      1       0
86  ham      0       0

For the first node, compute Gini index for the two variables,

> gini=function(variable){
+ T=table(db$Y,db[,variable]) + nx=apply(T,2,sum) + ProbCond=T/matrix(rep(nx,each=2),2,2) + ProbCond + Gini=-ProbCond*(1-ProbCond) + sum(matrix(rep(nx,each=2),2,2)/sum(nx)*Gini)} > gini("viagra")  -0.44 > gini("lottery")  -0.487 Here Gini index is maximal for “viagra”, so that will be the first node. # Visualising a Classification in High Dimension So far, when discussing classification, we’ve been playing on my toy-dataset (actually, I should no claim it’s mine, it is inspired by the one used in the introduction of Boosting, by Robert Schapire and Yoav Freund). But in ral life, there are more observations, and more explanatory variables.With more than two explanatory variables, it starts to be more complicated to visualise. For instance, consider MYOCARDE=read.table( "http://freakonometrics.free.fr/saporta.csv", head=TRUE,sep=";") where we have observations from people in E.R., for infarctus, and we want to understand who did survive, to get a predictive model. But before running some classifier, let us visualise our data. Since we have seven explanatory variables and our class (survival or death), we can go for a PCA. library(FactoMineR) # ACP (sur les var continues) X=MYOCARDE[,1:7] acp=PCA(X) To add the death/survival variable, treat it as numerical 0/1 variable (at least to get a direction) MYOCARDE2=MYOCARDE MYOCARDE2$PRONO=(MYOCARDE2$PRONO=="SURVIE")*1 acp=PCA(MYOCARDE2,quanti.sup=8,graph=TRUE) The nice thing is that we see here where variables are colinear with that one. It is also possible to visualise individuals, and classes, too acp=PCA(MYOCARDE,quali.sup=8,graph=TRUE) plot(acp, habillage = 8,col.hab=c("red","blue")) # Supervised Classification, beyond the logistic In our data-science class, after discussing limitations of the logistic regression, e.g. the fact that the decision boundary line was a straight line, we’ve mentioned possible natural extensions. Let us consider our (now) standard dataset  clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1)) clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2)) x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z <- c(1,1,1,1,1,0,0,1,0,0) df <- data.frame(x,y,z) plot(x,y,pch=19,cex=2,col=clr1[z+1]) One can consider a quadratic function of the covariates (instead of a linear one)  reg=glm(z~x+y+I(x^2)+I(y^2)+I(x*y), data=df,family=binomial) summary(reg) pred_1 <- function(x,y){ predict(reg,newdata=data.frame(x=x, y=y),type="response")>.5 } x_grid<-seq(0,1,length=101) y_grid<-seq(0,1,length=101) z_grid <- outer(x_grid,y_grid,pred_1) image(x_grid,y_grid,z_grid,col=clr2) points(x,y,pch=19,cex=2,col=clr1[z+1]) # Supervised Classification, discriminant analysis Another popular technique for classification (or at least, which used to be popular) is the (linear) discriminant analysis, introduced by Ronald Fisher in 1936. Consider the same dataset as in our previous post > clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1)) > x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) > y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) > z <- c(1,1,1,1,1,0,0,1,0,0) > df <- data.frame(x,y,z) > plot(x,y,pch=19,cex=2,col=clr1[z+1]) The main interest of that technique is not the output, but more the fact that we can make here simple (and explicit) computations. Especially to get a better understanding of theoretical concepts on classification. # Supervised Classification, Logistic and Multinomial We will start, in our Data Science course, to discuss classification techniques (in the context of supervised models). Consider the following case, with 10 points, and two classes (red and blue) > clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1)) > clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2)) > x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) > y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) > z <- c(1,1,1,1,1,0,0,1,0,0) > df <- data.frame(x,y,z) > plot(x,y,pch=19,cex=2,col=clr1[z+1]) To get a prediction, i.e. a partition of the space in two parts, consider some logistic regression > reg=glm(z~x+y,data=df,family=binomial) > summary(reg) Call: glm(formula = z ~ x + y, family = binomial, data = df) Deviance Residuals: Min 1Q Median 3Q Max -1.6593 -0.4400 0.2564 0.5830 1.5374 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.706 1.999 -0.854 0.393 x -5.489 5.360 -1.024 0.306 y 8.568 5.515 1.554 0.120 (Dispersion parameter for binomial family taken to be 1) Null deviance: 13.4602 on 9 degrees of freedom Residual deviance: 8.1445 on 7 degrees of freedom AIC: 14.144 Number of Fisher Scoring iterations: 5 Given some point, the predicted class is obtained using > pred_1 <- function(x,y){ + predict(reg,newdata=data.frame(x=x, + y=y),type="response")>.5 + } (here, the predicted class is simply the one that is the most likely). To visualize it use > x_grid<-seq(0,1,length=101) > y_grid<-seq(0,1,length=101) > z_grid <- outer(x_grid,y_grid,pred_1) > image(x_grid,y_grid,z_grid,col=clr2) > points(x,y,pch=19,cex=2,col=clr1[z+1]) Since the logistic regression is a (generalized) linear model, the line that separate the two regions is a straight line. # Visualizing Clusters Consider the following dataset, with (only) ten points x=c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y=c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) plot(x,y,pch=19,cex=2) We want to get – say – two clusters. Or more specifically, two sets of observations, each of them sharing some similarities. Since the number of observations is rather small, it is actually possible to get an exhaustive list of all partitions, and to minimize some criteria, such as the within variance. Given a vector with clusters, we compute the within variance using within_var = function(I){ I0=which(I==0) I1=which(I==1) xbar0=mean(x[I0]) xbar1=mean(x[I1]) ybar0=mean(y[I0]) ybar1=mean(y[I1]) w=sum(I0)*sum( (x[I0]-xbar0)^2+(y[I0]-ybar0)^2 )+ sum(I1)*sum( (x[I1]-xbar1)^2+(y[I1]-ybar1)^2 ) return(c(I,w)) } Then, to compute all possible partitions, use base2=function(z,n=10){ Base.b=rep(0,n) ndigits=(floor(logb(z, base=2))+1) for(i in 1:ndigits){ Base.b[ n-i+1]=(z%%2) z=(z%/%2)} return(Base.b)} L=function(x) within_var(base2(x)) S=sapply(1:(2^10),L) The cluster indices at the mimimum is here I=S[1:n,which.min(S[n+1,])]  To visualize those clusters, use cluster_viz = function(indices){ library(RColorBrewer) CL2palette=rev(brewer.pal(n = 9, name = "RdYlBu")) CL2f=CL2palette[c(1,9)] plot(x,y,pch=19,xlab="",ylab="",xlim=0:1,ylim=0:1,cex=2,col=CL2f[1+I]) CL2c=CL2palette[c(3,7)] I0=which(indices==0) I1=which(indices==1) xbar0=mean(x[I0]) xbar1=mean(x[I1]) ybar0=mean(y[I0]) ybar1=mean(y[I1]) segments(x[I0],y[I0],xbar0,ybar0,col=CL2c) segments(x[I1],y[I1],xbar1,ybar1,col=CL2c) points(xbar0,ybar0,pch=19,cex=1.5,col=CL2c) points(xbar1,ybar1,pch=19,cex=1.5,col=CL2c)} and then, simply cluster_viz(I) But that was possible only because is not to large (since the total number of scenarios – with only 2 clusters – is , or if we changes zeros in ones). # k-means clustering and Voronoi sets In the context of -means, we want to partition the space of our observations into classes. each observation belongs to the cluster with the nearest mean. Here “nearest” is in the sense of some norm, usually the (Euclidean) norm. Consider the case where we have 2 classes. The means being respectively the 2 black dots. If we partition based on the nearest mean, with the (Euclidean) norm we get the graph on the left, and with the (Manhattan) norm, the one on the right, Points in the red region are closer to the mean in the upper part, while points in the blue region are closer to the mean in the lower part. Here, we will always use the standard (Euclidean) norm. Note that the graph above is related to Voronoi diagrams (or Voronoy, from Вороний in Ukrainian, or Вороно́й in Russian) with 2 points, the 2 means. # Analyse des Données et Cartes Vendredi, on continuera en cours la classification (non supervisée) et en particulier les méthodes hiérarchiques, et les arbres. On va utiliser la base suivante, avec les résultats des élections présidentielles (premier tour) de 2012, elections2012 <- read.csv("http://komodo.regardscitoyens.org/public/presidentielles2012/resultats/resultats_departements_final_T1.csv", sep=";",header=TRUE,dec=",") La base qu’on garde contient juste les voix exprimées, en pourcentage voix <- which(substr(names(elections2012),1,12)=="X..Voix.Exp.") X <- as.matrix(elections2012[,voix]) colnames(X) <- substr(names(elections2012)[voix],13,nchar(names(elections2012)[voix])) rownames(X) <- elections2012[,3] Parfois, il convient de normaliser les données. Ici, sur les département, ça n’est pas nécessaire. Par contre, sur les candidats, ça le serait (comme pour l’ACP en fait). Pour visualiser les distances (entre départements), on peut utiliser heatmap(X) On peut faire une CAH, sur notre matrice de distance entre rangs cah <- hclust(dist(X)) plot(cah,cex=.6) et si on souhaite garder 5 groupes (par exemple) on utilise rect.hclust(cah,k=5) groups.5 <- cutree(cah,5) On peut aussi visualiser les classes avec la fonction library(dendroextras) plot(colour_clusters(cah,k=5)) Pour comprendre qui se trouve dans nos groupes, on peut utiliser aggregate(X,list(groups.5),mean) qui va nous renvoyer les “votes moyens” pour chaque candidat, dans chaque groupe. Maintenant, pour en finir avec la visualisation, on peut utiliser le code suivant (un peu long, certes) qui permettra de visualiser sur un carte les différents groupes carte_classe <- function(groupes){ library(stringr) elections2012$dep <- elections2012$Libellé.du.département 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("Libellé.du.département", "dep")] corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud" col2001<-groupes+1 names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms$Libellé.du.département)] color <- col2001[match(france$dep, names(col2001))]
map(database="france", fill=TRUE, col=color)
}

carte_classe(groups.5) Avec ces fonctions, on devrait pouvoir tester différentes méthodes pour constituer des groupes.

# Analyse des Données

Ce vendredi, on fera un TD en analyse des données. En plus des bases en ligne sur un ancien billet, je rajoute deux autres bases.

library(xts)
library(YieldCurve)
data(FedYieldCurve)
maturity <- c(3/12,6/12,1,2,3,5,7,10)
plot(maturity,FedYieldCurve,type="b")

load(url("http://freakonometrics.free.fr/titanic.rdata"))

# Analyse des Correspondances, suite

Hier soir, on voyait comment faire une analyse des correspondances à partir du tableau de contingence. Mais on peut avoir le problème autrement, à partir des individus. Certes, ces derniers ne sont pas observés (vraiment), mais peu importe.

> data(HairEyeColor)
> N = HairEyeColor[,,"Male"] + HairEyeColor[,,"Female"]
> Hair = rep(rep(rownames(N),ncol(N)),
+ as.vector(N))
> Eye =  rep(colnames(N),apply(N,2,sum))
> df = data.frame(Ind=1:sum(N),Hair,Eye)
> tail(df)
Ind  Hair   Eye
587 587 Blond Green
588 588 Blond Green
589 589 Blond Green
590 590 Blond Green
591 591 Blond Green
592 592 Blond Green

# Analyse des Correspondances

Lors du dernier cours d’analyse des données, on était parti sur l’analyse (simple) des correspondance, à partir d’un tableau de contingence, pour deux variables qualitatives,

On définit alors les effets marginaux,

pour les lignes, et pour les colonnes,

# Analyse des Données, plan de cours Ce semestre, je vais donner le cours d’Analyse des Données du Master Statistique & Économétrie. Le plan de cours sera, en gros,

1. Introduction à l’analyse des données non-supervisée
2. Réduction de dimension et Analyse en Composantes Principales
3. Analyse Factorielle et Analyse des Correspondances
4. Analyse Discriminante et k-Means
5. Classification Hiérarchique

Parmi les références utiles,