# Classification from scratch, logistic with splines 2/8

Today, second post of our series on classification from scratch, following the brief introduction on the logistic regression.

## Piecewise linear splines

To illustrate what’s going on, let us start with a “simple” regression (with only one explanatory variable). The underlying idea is natura non facit saltus, for “nature does not make jumps”, i.e. process governing equations for natural things are continuous. That seems to be a rather strong assumption, because we can assume that there is a fixed threshold to explain death. For instance, if patients die (for sure) if the “stroke index” exceeds a threshold, we might expect some discontinuity. Exceept that if that threshold is an heterogeneous (non-observable continuous) variable, then we get back to the continuity assumption.

The most simple model we can think of to extend the linear model we’ve seen in the previous post is to consider a piecewise linear function, with two parts : small values of $x$, and larger values of $x$. The most convenient way to do so is to use the positive part function $(x-s)_+$ which is the difference between $x$ and $s$ if that difference is positive, and $0$ otherwise. For instance $$\beta_1 x+\beta_2(x-s)_+$$ is the following piecewise linear function, continuous, with a “rupture” at knot $s$.

Observe also the following interpretation: for small values of $x$, there is a linear increase, with slope $\beta_1$, and for lager values of $x$, there is a linear decrease, with slope $\beta_1+\beta_2$. Hence, $\beta_2$ is interpreted as a change of the slope.

And of course, it is possible to consider more than one knot. The function to get the positive value is the following

pos = function(x,s) (x-s)*(x&gt;=s)

then we can use it direcly in our regression model

reg = glm(PRONO~INSYS+pos(INSYS,15)+ pos(INSYS,25),data=myocarde,family=binomial)

The output of the regression is here

summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.1109 3.2783 -0.034 0.9730 INSYS -0.1751 0.2526 -0.693 0.4883 pos(INSYS, 15) 0.7900 0.3745 2.109 0.0349 * pos(INSYS, 25) -0.5797 0.2903 -1.997 0.0458 *

Hence, the original slope, for very small values is not significant, but then, above 15, it become significantly positive. And above 25, there is a significant change again. We can plot it to see what’s going on

u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,type="l") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Using bs() linear splines

Using the GAM function, things are slightly different. We will use here so called b-splines,

library(splines)

We can define spline functions with support $(5,55)$ and with knots $\{15,25\}$

clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02") x = seq(0,60,by=.25) B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1) matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

as we can see, the functions defined here are different from the one before, but we still have (piecewise) linear functions on each segment $(5,15)$, $(15,25)$ and $(25,55)$. But linear combinations of those functions (the two sets of functions) will generate the same space. Said differently, if the interpretation of the output will be different, predictions should be the same

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=1), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.9863 2.0555 -0.480 0.6314 bs(INSYS,..)1 -1.7507 2.5262 -0.693 0.4883 bs(INSYS,..)2 4.3989 2.0619 2.133 0.0329 * bs(INSYS,..)3 5.4572 5.4146 1.008 0.3135

Observe that there are three coefficients, as before, but again, the interpretation is here more complicated…

v=predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Nevertheless, the prediction is the same… and that’s nice.

Let us go one step further… Can we have also the continuity of the derivative ? Yes, and that’s easy actually, considering parabolic functions. Instead of using a decomposition on $x,(x-s_1)_+$ and $(x-s_2)_+$ consider now a decomposition on $x,x^{\color{red}{2}},(x-s_1)^{\color{red}{2}}_+$ and $(x-s_2)^{\color{red}{2}}_+$.

 pos2 = function(x,s) (x-s)^2*(x&gt;=s) reg = glm(PRONO~poly(INSYS,2)+pos2(INSYS,15)+pos2(INSYS,25), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 29.9842 15.2368 1.968 0.0491 * poly(INSYS, 2)1 408.7851 202.4194 2.019 0.0434 * poly(INSYS, 2)2 199.1628 101.5892 1.960 0.0499 * pos2(INSYS, 15) -0.2281 0.1264 -1.805 0.0712 . pos2(INSYS, 25) 0.0439 0.0805 0.545 0.5855

As expected, there are here five coefficients: the intercept and two for the part on the left (three parameters for the parabolic function), and then two additional terms for the part in the center – here $(15,25)$ – and for the part on the right. Of course, for each portion, there is only one degree of freedom since we have a parabolic function (three coefficients) but two constraints (continuity, and continuity of the first order derivative).

On a graph, we get the following

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Of course, we can do the same with our R function. But as before, the basis of function is expressed here differently

 x = seq(0,60,by=.25) B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=2) matplot(x,B,type="l",xlab="INSYS",col=clr6)

If we run R code, we get

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=2),data=myocarde, family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 7.186 5.261 1.366 0.1720 bs(INSYS, ..)1 -14.656 7.923 -1.850 0.0643 . bs(INSYS, ..)2 -5.692 4.638 -1.227 0.2198 bs(INSYS, ..)3 -2.454 8.780 -0.279 0.7799 bs(INSYS, ..)4 6.429 41.675 0.154 0.8774

But that’s not really a big deal since the prediction is exactly the same

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Cubic splines

Last, but not least, we can reach the cubic splines. With our previous notions, we would consider a decomposition on (guess what) $x,x^2,x^{\color{red}{3}},(x-s_1)^{\color{red}{3}}_+,(x-s_2)^{\color{red}{3}}_+$, to get this time continuity, as well as continuity of the first two derivatives (and to get a very smooth function, since even variations will be smooth). If we use the bs function, the basis is the followin

B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=3) matplot(x,B,type="l",lwd=2,col=clr6,lty=1,ylim=c(-.2,1.2)) abline(v=c(5,15,25,55),lty=2)

and the prediction will now be

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=3), data=myocarde,family=binomial) u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Two last things before concluding (for today), the location of the knots, and the extension to additive models.

## Location of knots

In many applications, we do not want to specify the location of the knots. We just want – say – three (intermediary) knots. This can be done using

reg = glm(PRONO~1+bs(INSYS,degree=1,df=4),data=myocarde,family=binomial)

We can actually get the locations of the knots by looking at

attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 1L, knots = c(15.8, 21.4, 27.15), Boundary.knots = c(8.7, 54), intercept = FALSE) which provides us with the location of the boundary knots (the minumun and the maximum from from our sample) but also the three intermediary knots. Observe that actually, those five values are just (empirical) quantiles quantile(myocarde$INSYS,(0:4)/4) 0% 25% 50% 75% 100% 8.70 15.80 21.40 27.15 54.00

If we plot the prediction, we get

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2) If we get back on what was computed before the logit transformation, we clealy see ruptures are the different quantiles B = bs(x,degree=1,df=4) B = cbind(1,B) y = B%*%coefficients(reg) plot(x,y,type="l",col="red",lwd=2) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2)

Note that if we do specify anything about knots (number or location), we get no knots…

reg = glm(PRONO~1+bs(INSYS,degree=2),data=myocarde,family=binomial) attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 2L, knots = numeric(0), Boundary.knots = c(8.7,54), intercept = FALSE) and if we look at the prediction u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) actually, it is the same as a quadratic regression (as expected actually) reg = glm(PRONO~1+poly(INSYS,degree=2),data=myocarde,family=binomial) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) ## Additive models Consider now the second dataset, with two variables. Consider here a model like $$\mathbb{P}[Y|X_1=x_1,X_2=x_2]=\frac{\exp[\eta(x_1,x_2)]}{1+\exp[\eta(x_1,x_2)]}$$ where $$\exp[\eta(x_1,x_2)]=\beta_0+\color{red}{s_1(x_1)}+\color{blue}{s_2(x_2)}$$ $$\color{red}{s_1(x_1)}=\beta_{1,0}x_1+\beta_{1,1}(x_1-s_{11})_++\beta_{1,2}(x_1-s_{12})_+$$ and $$\color{blue}{s_2(x_2)}=\beta_{2,0}x_2+\beta_{2,1}(x_2-s_{21})_++\beta_{2,2}(x_2-s_{22})_+$$ It might seem a little bit restrictive, but that’s actually the idea of additive models. reg = glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Now, if think about is, we’ve been able to get a “perfect” model, so, somehow, it seems no longer continuous…

persp(u,u,v,theta=20,phi=40,col="green"

Of course, it is… it is piecewise linear, with hyperplane, some being almost vertical.

And one can also consider piecewise quadratic functions

reg = glm(y~bs(x1,degree=2,df=3)+bs(x2,degree=2,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) Funny thing, we now have two “perfect” models, with different areas for the white and the black dots… Don’t ask me how to choose on that one. In R, it is possible to use the mgcv package to run a gam regression. It is used for generalized additive models, but here, we have only one variable, so it is difficult to see the “additive” part, actually. And to be more specific, mgcv is using penalized quasi-likelihood from the nlme package (but we’ll get back on penalized routines later on). But maybe I should also mention another smoothing tool before, kernels (and maybe also $k$-nearest neighbors). To be continued # Classification from scratch, logistic regression 1/8 Let us start today our series on classification from scratch The logistic regression is based on the assumption that given covariates $\mathbf{x}$, $Y$ has a Bernoulli distribution,$$Y|\mathbf{X}=\mathbf{x}\sim\mathcal{B}(p_{\mathbf{x}}),~~~~p_\mathbf{x}=\frac{\exp[\mathbf{x}^T\mathbf{\beta}]}{1+\exp[\mathbf{x}^T\mathbf{\beta}]}$$The goal is to estimate parameter $\mathbf{\beta}$. Recall that the heuristics for the use of that function for the probability is that$$\log[\text{odds}(Y=1)]=\log\frac{\mathbb{P}[Y=1]}{\mathbb{P}[Y=0]}=\mathbf{x}^T\mathbf{\beta}$$ ## Maximimum of the (log)-likelihood function The log-likelihood is here$$\log\mathcal{L} = \sum_{i=1}^n y_i\log p_i+(1-y_i)\log (1-p_i)$$ where $p_{i}=(1+\exp[-\mathbf{x}_i^T\mathbf{\beta}])^{-1}$. Numerical techniques are based on (numerical) gradient descent to compute the maximum of the likelihood function. The (negative) log-likelihood is the following function y = myocarde$PRONO X = cbind(1,as.matrix(myocarde[,1:7])) negLogLik = function(beta){ -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta))) }

We use the minus sign since standard optimization routines compute minima, not maxima. Now, to find the minimum of that function, we need a starting point to initiate the algorithm

beta_init = lm(PRONO~.,data=myocarde)$coefficients Why not start with the parameter of the OLS. Somehow, we might think that at least, sign should be ok for instance. Anyway, we need a starting point, and let us use that one. logistic_opt = optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) Here, we obtain  logistic_opt$par (Intercept) FRCAR INCAR INSYS 1.656926397 0.045234029 -2.119441743 0.204023835 PRDIA PAPUL PVENT REPUL -0.102420095 0.165823647 -0.081047525 -0.005992238

Let us verify here that this output is valid. For instance, what if we change the value of the starting point (randomly)

simu = function(i){ logistic_opt_i = optim(par = rnorm(8,0,3)*beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) logistic_opt_i$par[2:3] } v_beta = t(Vectorize(simu)(1:1000)) plot(v_beta) par(mfrow=c(1,2)) hist(v_beta[,1],xlab=names(myocarde)[1]) hist(v_beta[,2],xlab=names(myocarde)[2]) Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine library(optimx) logit = function(mX, vBeta) { exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) } logLikelihoodLogitStable = function(vBeta, mX, vY) { -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta))) + (1-vY)*(-log(1 + exp(mX %*% vBeta)))) } likelihoodScore = function(vBeta, mX, vY) { return(t(mX) %*% (logit(mX, vBeta) - vY) ) } optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable, method = 'L-BFGS-B', gr = likelihoodScore, mX = X, vY = y, hessian=TRUE) The optimum is here attr(optimLogitLBFGS, "details")[[2]] [,1] 0.066680272 FRCAR 0.003080542 INCAR 0.079031364 INSYS -0.001586194 PRDIA 0.040500697 PAPUL -0.041870705 PVENT -0.014162756 REPUL 0.195632244 Let’s be honest here, I do not feel confortable with those techniques. So, what happened here ? Here, the technique we use is based on the following idea,$$\mathbf{\beta}_{new}=\mathbf{\beta}_{old} -\left(\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}\right)^{-1}\cdot \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}$$The problem is that my computer does not know this first and second derivatives. So it will compute them using approximation techniques. Actually, it is possible to use functions dedicated to such computation library(numDeriv) library(MASS) logit = function(x){1/(1+exp(-x))} logLik = function(beta, X, y){ -sum(y*log(logit(X%*%beta)) + (1-y)*log(1-logit(X%*%beta))) } optim_second = function(beta, num_iter){ LL = vector() for(i in 1:num_iter){ grad = (t(X)%*%(logit(X%*%beta) - y)) H = hessian(logLik, beta, method = "complex", X = X, y = y) beta = beta - ginv(H)%*%grad LL[i] = logLik(beta, X, y) } result = list(beta, H) return(result) } With our OLS starting point, we obtain opt0 = optim_second(beta_init,500) opt0[[1]] [,1] [1,] 0.951074420 [2,] 0.018860280 [3,] 0.275428978 [4,] 0.144803636 [5,] -0.058535606 [6,] 0.001182178 [7,] -0.108651776 [8,] -0.002940315 But if we try with another starting point opt1 = optim_second(beta_init*runif(8),500) opt1[[1]] [,1] [1,] 0.052894794 [2,] 0.024718435 [3,] 0.167953661 [4,] 0.171662947 [5,] -0.057458066 [6,] -0.011361034 [7,] -0.107532114 [8,] -0.002679064 Clearly, some coefficients are rather close. But other aren’t. From my point of viezw, that is a major problem (keep in mind that we do not deal here with massive data ! There are only 7 explanatory variables, and only 71 observations). Why not try to be clever, and use the analytical values of those derivatives ? Even if some people claim the oppositive, sometimes, it can actually be usefull to do the maths, instead of considering only numerical values. ## Newton (or Fisher) Algorithm If you open any Econometrics textbooks (one can also try to derive it), you will get $$\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})$$ while$$\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}$$ Y=myocarde$PRONO X=cbind(1,as.matrix(myocarde[,1:7])) colnames(X)=c("Inter",names(myocarde[,1:7])) beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} Observe that here, I use only ten iterations of the algorithm !  beta[,8:10] [,1] [,2] [,3] XInter -10.187641685 -10.187641696 -10.187641696 XFRCAR 0.138178119 0.138178119 0.138178119 XINCAR -5.862429035 -5.862429037 -5.862429037 XINSYS 0.717084018 0.717084018 0.717084018 XPRDIA -0.073668171 -0.073668171 -0.073668171 XPAPUL 0.016756506 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187 The thing is that is seems to converge extremely fast. And it is rather robust ! Look at what we get if we change our starting point beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} beta[,8:10] [,1] [,2] [,3] XInter -10.187641586 -10.187641696 -10.187641696 XFRCAR 0.138178118 0.138178119 0.138178119 XINCAR -5.862429017 -5.862429037 -5.862429037 XINSYS 0.717084013 0.717084018 0.717084018 XPRDIA -0.073668172 -0.073668171 -0.073668171 XPAPUL 0.016756508 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187

Nice, isn’t it? Looks like we got our winner, don’t we? And one can use the inverse of the Hessian matrix to get standard deviations.

## Weighted Least-Squares

Let us go one step further. We’ve seen that we want to compute something like$$\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}$$(if we do substitute matrices in the analytical expressions) where $\mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]$. But actually, that’s simply a standard least-square problem$$\mathbf{\beta}_{new} = \text{argmin}\left\lbrace(\mathbf{z}-\mathbf{X}\mathbf{\beta})^T\mathbf{\Delta}_{old}^{-1}(\mathbf{z}-\mathbf{X}\mathbf{\beta})\right\rbrace$$The only problem here is that weights $\mathbf{\Delta}_{old}$ are functions of unknown $\mathbf{\beta}_{old}$. But actually, if we keep iterating, we should be able to solve it : given the $\mathbf{\beta}$ we got the weights, and with the weights, we can use weighted OLS to get an updated $\mathbf{\beta}$. That’s the idea of iteratively reweighted least squares.

The algorithm will be

df = myocarde beta_init = lm(PRONO~.,data=df)$coefficients X = cbind(1,as.matrix(myocarde[,1:7])) beta = beta_init for(s in 1:1000){ p = exp(X %*% beta) / (1+exp(X %*% beta)) omega = diag(nrow(df)) diag(omega) = (p*(1-p)) df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p) beta = lm(Z~.,data=df[,-8], weights=diag(omega))$coefficients }

and the output is here

 beta (Intercept) FRCAR INCAR INSYS PRDIA -10.187641696 0.138178119 -5.862429037 0.717084018 -0.073668171 PAPUL PVENT REPUL 0.016756506 -0.106776012 -0.003154187

which is almost what we’ve obtained before. Nice isn’t it ? Actually, here we also have standard deviations of estimators

summary( lm(Z~.,data=df[,-8], weights=diag(omega)))   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) -10.187642 10.668138 -0.955 0.343 FRCAR 0.138178 0.102340 1.350 0.182 INCAR -5.862429 6.052560 -0.969 0.336 INSYS 0.717084 0.503527 1.424 0.159 PRDIA -0.073668 0.261549 -0.282 0.779 PAPUL 0.016757 0.306666 0.055 0.957 PVENT -0.106776 0.099145 -1.077 0.286 REPUL -0.003154 0.004386 -0.719 0.475

## The standard glm function

Of course, it is possible to use an R built-in function to get our estimate

summary(glm(PRONO~.,data=myocarde,family=binomial(link = "logit")))   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -10.187642 11.895227 -0.856 0.392 FRCAR 0.138178 0.114112 1.211 0.226 INCAR -5.862429 6.748785 -0.869 0.385 INSYS 0.717084 0.561445 1.277 0.202 PRDIA -0.073668 0.291636 -0.253 0.801 PAPUL 0.016757 0.341942 0.049 0.961 PVENT -0.106776 0.110550 -0.966 0.334 REPUL -0.003154 0.004891 -0.645 0.519

## Application and visualisation

Let us visualize the prediction obtained from the logistic regression, on our second dataset

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(x1=x,x2=y,y=as.factor(z)) reg = glm(y~x1+x2,data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Here level curves – or iso-probabilities – are linear, so the space is divided in two (0 and 1, survival and death, white and black) by a straight line (or an hyperplane in higher dimension). Furthermore, since we have a linear model, if we change the cutoff (the threshold used to create the two classes), we obtain another straight line (or hyperplane) parallel to the first one.

Next time, we will introduce splines to smooth those continuous covariates… to be continued.

# Classification from scratch, overview 0/8

Before my course on « big data and economics » at the university of Barcelona in July, I wanted to upload a series of posts on classification techniques, to get an insight on machine learning tools.

According to some common idea, machine learning algorithms are black boxes. I wanted to get back on that saying. First of all, isn’t it the case also for regression models, like generalized additive models (with splines) ? Do you really know what the algorithm is doing ? Even the logistic regression. In textbooks, we can easily find math formulas. But what is really done when I run it, in R ?

When I started working on academia, someone told me something like « if you really want to understand a theory, teach it ». And that has been my moto for more than 15 years. I wanted to add a second part to that statement: « if you really want to understand an algorithm, recode it ». So let’s try this… My ambition is to recode (more or less) most of the standard algorithms used in predictive modeling, from scratch, in R. What I plan to mention, within the next two weeks, will be

I will use two datasets to illustrate. The first one is inspired by the cover of « Foundations of Machine Learning » by Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar. At least, with this dataset, it will be possible to plot predictions (since there are only two – continuous – features)

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(x1=x,x2=y,y=as.factor(z)) plot(x,y,pch=c(1,19)[1+z])

Here is some code to get a visualization of the prediction (here the probability to be a black point)

rmatrix_model = function(model){ u = seq(0,1,length=101) p = function(x,y) predict(model,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) return(v)} nice_graph=function(v){ u = seq(0,1,length=101) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10[c(1,10)],breaks=c(0,5,10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) } reg = glm(y~x1+x2,data=df,family=binomial) nice_graph(rmatrix_model(reg))

Note that colors are defined here as

clr10= c("#ffffff","#f7fcfd","#e5f5f9","#ccece6","#99d8c9","#66c2a4","#41ae76","#238b45","#006d2c","#00441b")

or with some nonlinear model

The second one is a dataset I got from Gilbert Saporta, about heart attacks and decease (our binary variable).

myocarde = read.table("http://freakonometrics.free.fr/myocarde.csv",head=TRUE, sep=";") myocarde$PRONO = (myocarde$PRONO=="SURVIE")*1 y = myocarde$PRONO X = as.matrix(cbind(1,myocarde[,1:7])) So far, I do not plan to talk (too much) on the choice of tunning parameters (and cross-validation), on comparing models, etc. The goal here is simply to understand what’s going on when we call either glm, glmnet, gam, random forest, svm, xgboost, or any function to get a predict model. # 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[1]) segments(x[I1],y[I1],xbar1,ybar1,col=CL2c[2]) points(xbar0,ybar0,pch=19,cex=1.5,col=CL2c[1]) points(xbar1,ybar1,pch=19,cex=1.5,col=CL2c[2])} and then, simply cluster_viz(I) But that was possible only because $n$ is not to large (since the total number of scenarios – with only 2 clusters – is $2^n$, or $2^{n-1}$ if we changes zeros in ones). # Les Arbres de Classification J’animerai une formation lundi 28 de 14:00 à 16:00 au local N-6320 de l’UQAM sur le thème introduction aux arbres de classification. Cette formation est organisée dans le cadre des séminaires en méthodes d’analyses quantitatives et qualitatives qui se tiennent régulièrement depuis un peu plus d’un mois. animé par le collectif pour le développement et les applications en mesure et évaluation (Cdame). Les slides sont disponibles en pdf (il y a quelques animations, qui ne passent qu’avec Acrobat) La base utilisée tout au long des exposés est la suivante > MYOCARDE=read.table("http://freakonometrics.free.fr/saporta.csv",head=TRUE,sep=";") # Bar des Sciences: Débat sur le Big Data Le Cœur des Sciences, à Université du Québec à Montréal, organise le 13 février prochain, à 18h, un débat grand public sur le Big Data, dans le cadre d’un bar des sciences, auquel je devrais participer, avec Vincent Gautrais (a.k.a. @gautrais), Yves-Alexandre de Montjoye (a.k.a. @yvesalexandre) et Jean-Hughes Roy (a.k.a. @jeanhuguesroy). Les mauvaises langues diront que je n’aurais pas pu refuser d’intervenir dans un bar (et elles n’auront probablement pas tort). Cela dit, le public ne semble plus autorisé : dans les heures qui ont suivi l’annonce, l’événement a été complet ! Je vais en profiter, aujourd’hui, pour livrer quelques éléments de réflexion… les commentaires sont ouverts, mais la suite du débat se fera au Cœur des Sciences (j’essayerais de faire un résumé sur le blog). Car c’est la première fois que je fais une intervention sur le sujet (de manière aussi explicite). Je ne peux m’empêcher d’avoir une pensée pour les mots de Dan Ariely, sur le Big Data, Pour avoir discuté avec de nombreux professionnels, dans l’industrie, qui manipulent des données quotidiennement, je les ai vu passer par les cinq grandes étapes décrites par Zubin Dowalty, 1. le déni, “There’s nothing in that big data that we don’t already know. 2. la colère, “There’s nothing in that big data that we don’t already know!” (avec un point d’exclamation cette fois) 3. la négociation, “If we could just get the budget to expand the project… 4. la dépression, “These data sets are just overwhelming. There’s no way we can do this.” 5. l’acceptation, “This isn’t going to happen overnight. We need to be strategic about how and when to undertake big data analysis. C’est devenu un cliché que l’utilisation de l’informatique exploser, de manière exponentielle. Aussi bien au niveau du nombre d’utilisateurs, de la vitesse de calcul, des capacités de stockage que du volume de données récoltées via internet Si on admet qu’il y a matière à réflexion, il va falloir dire davantage de quoi parle, car si on fonctionne par association d’idées, le terme Big Data peut évoquer beaucoup de choses. En tant que statisticien, il me fait penser à des problèmes assez profonds de modélisation, de high frequency-data (sans parler de HFT), de différence entre causalité et corrélation, ainsi que des aspects importants de statistique computationnelle (difficulté de faire tourner des algorithmes sur des très très grosses matrices, sans parler des mathématiques plus complexes qui vont avec). En tant que spécialiste de l’analyse de données, ce terme m’évoque deux autres termes connexes, que sont l’open data et le hacking de données. En tant qu’économiste, je pense également au business associé aux logiciels d’analyse de données, aux espaces publicitaires de Google, aux suggestions que me fait Amazon dès que je commande un livre, à l’utilisation des informations en ligne (par exemple les amis sur Facebook) dans les scores de crédit, etc. Mais l’accent qui a été mis ici (de part la seconde partie du titre), c’est l’aspect Big Brother. Donc un aspect important sera probablement autour des applications mobiles, du marketing, d’internet, et des algorithmes qui semblent nous contrôler. Bref, on va parler de data science. Sean Owen notait (toujours avec ironie) On pourrait donc croire qu’il n’y a rien de nouveau sous le soleil…? • le Big Data avec les yeux d’un statisticien Pour reprendre les mots d’Alan Mitchell, dans big data, big dead end, Big data is all about statistics: divining patterns and trends from large data sets. Statistics are incredibly powerful and useful for the way they challenge the assumptions and inferences naturally made by human minds – many of them faulty. As I said, that’s great. Avant de rentrer dans le vif du sujet, deux petites réflexions pour débuter. Autour du mot Big pour commencer. En mathématiques financières, beaucoup de modèles de valorisation ont été développés en temps continu, avec des processus de prix $(S_t)_{t\in\mathbb{R}^+}$. Délicat de réconcilier avec l’économétrie financière, où les prix sont observés à des dates $\{0\leq t_1\leq t_2\leq \cdots\leq t_n\leq\cdots\}$. En statistique mathématique et en économétrie, la majorité des théorèmes sont asymptotiques, i.e. valides lors que la taille de l’échantillon $\{x_1,\cdots,x_n\}$ tend vers l’infini. Aussi, traditionnellement, avoir des Small Data était un soucis, car il semble que les modèles ont été pensé, précisément, sous l’hypothèse qu’un grand nombre d’observations serait disponible. Narinder Singh, écrivait en décembre dernier Big Data is misnamed in our (academic) world, because data sets have always been big. What is different is that we now have the technology to simply run every scenario. Mais ça n’est pas si simple que ça, loin de là. Le terme Data tout d’abord s’est élargie. Les données étaient auparavant des niveaux de production industrielle, ou des prix de matière première. Classiquement, ces jeux de données était présenté sous forme matricielle, un rectangle de $n$ lignes (les observations) et $k$ colonnes (les variables). Historiquement, on avait un nombre important (par forcément Big) de lignes, un peu moins de variables pour les expliquer. Un certain nombre de problèmes se sont posé en génétique, quand $k$ s’est mis à croitre plus rapidement que $n$. Et les données on commencé à devenir plus compliquées. On pourra penser aux données de l’assurance maladie, avec des consultations médicales, des consommations de médicaments (consécutives à une visite chez le médecin), puis la visite d’un spécialiste, etc. Des transcriptions d’appels téléphonique à un centre d’appels, ou des collections de livres dans différentes langues éventuellement. Des images, ou des vidéos. Voire de l’imagerie médicale en trois dimensions (en plus du temps). Des achats de consommateurs sur un site de vente en ligne (des grosses matrices creuses, car non informées). Bref, le mot Data, il est devenu plus complexe qu’il ne l’était auparavant (nous reviendrons un peu sur les outils mathématiques pour les appréhender). 1. Le Big Data, renoncer à chercher des cause et identifier des corrélations (dont j’avais déjà parlé rapidement dans un billet suite à ma lecture des livres de Nate Silver, et de Kenneth Cukier et Viktor Mayer-Schönberger). Pour résumer rapidement, Kenneth Cukier et Viktor Mayer-Schönberger notaient dans mise en données du monde, le déluge numérique, Pareil usage suppose trois changements majeurs dans notre approche. Le premier consiste à recueillir et à utiliser le plus grand nombre possible d’informations plutôt que d’opérer un tri sélectif comme le font les statisticiens depuis plus d’un siècle. Le deuxième implique une certaine tolérance à l’égard du désordre : mouliner des données innombrables, mais de qualité inégale, s’avère souvent plus efficace qu’exploiter un petit échantillon impeccablement pertinent. Enfin, le troisième changement implique que, dans de nombreux cas, il faudra renoncer à identifier les causes et se contenter de corrélations. Au lieu de chercher à comprendre précisément pourquoi une machine ne fonctionne plus, les chercheurs peuvent collecter et analyser des quantités massives d’informations relatives à cet événement et à tout ce qui lui est associé afin de repérer des régularités et d’établir dans quelles circonstances la machine risque de retomber en panne. Ils peuvent trouver une réponse au « comment », non au « pourquoi » ; et, bien souvent, cela suffit. […] Ce changement d’approche à l’égard des données numériques — exhaustives et non plus échantillonnées, désordonnées et non plus méthodiques — explique le glissement de la causalité vers la corrélation. On s’intéresse moins aux raisons profondes qui président à la marche du monde qu’aux associations susceptibles de relier entre eux des phénomènes disparates. L’objectif n’est plus de comprendre les choses, mais d’obtenir une efficacité maximale. On y retrouve ici les propos de Chris Anderson, tenu en 2008 dans The End of Theory: The Data Deluge Makes the Scientific Method Obsolete, In short, the more we learn about biology, the further we find ourselves from a model that can explain it. There is now a better way. Petabytes allow us to say: “Correlation is enough.” We can stop looking for models. We can analyze the data without hypotheses about what it might show. We can throw the numbers into the biggest computing clusters the world has ever seen and let statistical algorithms find patterns where science cannot. C’est un changement profond de mentalité. En 1963, Karl Popper (bien que traditionnellement vu comme un philosophe critique de l’inductivisme) affirmait qu’il était inutile de recueillir des données en espérant faire ressortir des similitudes, permettant ensuite de faire émerger une nouvelle théorie : the belief that we can start with pure observations alone, without anything in the nature of a theory, is absurd; as may be illustrated by the story of the man who dedicated his life to natural science, wrote down everything he could observe, and bequeathed his priceless collection of observations to the Royal Society to be used as inductive evidence. This story should show us that though beetles may profitably be collected, observations may not. Malgré tout, le big data propose un changement de paradigme sur la méthode scientifique, peut être plus que sur la modélisation statistique. 2. Petit problème épistémologique Au delà du débat corrélation-causalité, la statistique pose des problèmes fondamentaux d’épistémologie, pour les philosophes des sciences. Pour reprendre l’idée de Karl Popper, dans The Logic of Scientific Discovery, reprise dans Gilles (1971) for although probability statements play such a vitally important role in empirical science, they turn out to be impervious to strict falcification Un exemple est donné par Katja de Vries: le postulat “tous les cygnes sont blancs” a été invalidé dès lors qu’un cygne noir a été observé pour la première fois. C’est ainsi que fonctionne la science. Mais si on reprend l’exemple suivant (pris chez un autre grand philosophe) on note qu’il est difficile (“impossible” dirait Karl Popper) d’invalider une hypothèse statistique. Pour pouvoir affirmer avec certitude que la probabilité qu’un pièce tombe sur “face” est 50%, il faut une infinité de tirages de cette pièce: only an infinite sequence of events (…) could contracit a probability estimate C’est d’ailleurs l’argumentaire que l’on retrouve abordé avec maints exemples dans Merchants of Doubt, de Naomi Oreskes et Erik Conway (je reviendrais sur ce livre dans les semaines qui viennent). C’est probablement à cause de cela que beaucoup de prévisionnistes (je préfère le terme “prédicateurs“) sont vus comme des charlatans. Si j’affirme qu’il y a 3 chances sur 10000 pour qu’un incendie se déclare au pavillon du Coeur des Science le 13 février, il sera difficile de valider ou d’invalider mon estimation. Quoi qu’il se passe, je pourrais toujours dire “je vous l’avais dit” (soit qu’un incendie était possible, soit qu’il était extrêment improbable). On va alors quiter le monde où un évènement serait “impossible” à celui où il serait “hautement improbable“. Ce qu’un scientifique appelerait “peu vraisemblable” correspondrait à ce que Karl Popper appelerait “practical falsified“, It is fairly clear that this ‘practical falsification’ can be obtained only through a methodological decision to regard highly improbable events as ruled out – as prohibited. But with what right can they be so regarded? Where are we to draw the line? Where does this ‘high improbability’ begin? En pratique, la réponse a été apportée par l’utilisation d’une p-value, une sorte de ‘statistical significant falsification‘, pour reprendre les termes de Karl Popper. Nous allons revenir sur ce point dans quelques paragraphes. 3. Du modèle paramétrique à la statistique non-paramétrique Dans the end of theory: the data deluge makes the scientific method obsolete, Chris Anderson prétendait que Data without a model is just noise. But faced with massive data, this approach to science — hypothesize, model, test — is becoming obsolete. C’est là encore probablement un peu plus complexe qu’il n’y parait. La commission Cowles, qui a posé les bases de la théorie économétrique (tout en fondant l’Econometric Society, et la prestigieuse revue Econometrica) partait du postulat qu’un modèle économétrique devait être le reflet d’un modèle économique. C’est ce qu’on retrouve dans les SEM (Structural Equation Model) comme le modèle de Klein, proposé en 1950, Techniquement, on parlera de modèle paramétrique, puisqu’à partir de jeux de données, les seuls inconnues sont les paramètres $\gamma$ et $\beta$. La transition vers les modèles non-paramétriques, dans lesquels les données l’emportent sur un éventuel modèle, ne s’est pas fait sans douleur. Leo Breiman l’analyse dans statistical modeling: the two cultures, paru en 2001. Pour comprendre la distinction, on peut reprendre l’exemple classique du modèle linéaire. Considérons les 50 observations suivantes, avec des couples de variables $(X,Y)$. Dans le modèle linéaire, on suppose que $Y_i=\beta_0+\beta_1X_i+\varepsilon_i$ mais on peut aussi considérer un modèle quadratique $Y_i=\beta_0+\beta_1X_i+\beta_2X_i^2+\varepsilon_i$ (ou toute autre relation fonctionnelle ayant une forme paramétrique). Si on suppose que $\varepsilon_i$ est un bruit, imprévisible, on peut espérer que ce bruit soit centré, au sens où – en moyenne – ce bruit doit être nul. Pour ajuster notre modèle, on va alors chercher à minimiser la somme des carrés de ces erreurs, on on cherche $(\widehat{\beta}_0,\widehat{\beta}_1)=\text{argmin}\left\{\sum_{i=1}^n(Y_i-[\beta_0+\beta_1X_i])^2\right\}$ On peut alors montrer que cette méthode permet d’obtenir la valeur moyenne de $Y$, à $X$ donné, $\widehat{Y}(x)=\widehat{\beta}_0+\widehat{\beta}_1x$ Mais si on veut connaître la moyenne de $Y$ sachant $X=x$, pourquoi par juste prendre une moyenne, calculée sur les $X_i$ proches de la valeur x qui nous intéresse ? C’est ce qu’on va appeler un estimateur non-paramétrique de la “droite” de régression $\widetilde{Y}(x)=\frac{\displaystyle{\sum_{i=1}^n Y_i\cdot \boldsymbol{1}(\vert X_i-x\vert\leq h)}}{\displaystyle{\sum_{i=1}^n \boldsymbol{1}(\vert X_i-x\vert\leq h)}}$ On ne fait pas ici de modèle, et la prévision sera imposée par nos données. Pour aller plus loin sur les outils statistiques, on pourra relire Hal Varian. 4. La p-value à l’ère du Big Data Dans too big to fail: large samples and the p-value problem, Mingfeng Lin, Henry Lucas, Jr. et Galit Shmueli notent le point suivant A key issue with applying small-sample statistical inference to large samples is that even minuscule effects can become statistically significant. The increased power leads to a dangerous pitfall as well as to a huge opportunity. The issue is one that statisticians have long been aware of: “the p-value problem.” Chatfield (1995, p. 70) comments, ‘The question is not whether differences are ‘significant’ (they nearly always are in large samples), but whether they are interesting. Forget statistical significance, what is the practical significance of the results?’ Sur l’utilisation de la p-value, ainsi que son histoire, je renvois à quelques billets. Le problème de la taille de l’échantillon est plus technique. Je remets ça à un billet dans les mois qui viennent… 5. Le Big-Data pour remplacer des algorithmes qui ne peuvent être améliorés Une des difficulté est de modéliser les prises de décisions d’êtres humains. Dans des comportements d’achat (nous allons y revenir un peu plus tard), mais aussi dans le cadre de jeux de stratégie. Les jeux sont plus simples que la vraie vie, car les règles y sont claires et non modifiables. En particulier le jeu d’échec. Comme le notent Viktor Meyer-Schönberger et Kenneth Cukier dans Big data, a revolution that will transform how we live, work and think For example, chess algorithms have changed only slightly in the past few decades, since the rules of chess are fully known and tightly constrained. The reason computer chess programs play far better today than in the past is in part that they are playing their endgame better. And they’re doing that simply because the systems have been fed more data. In fact, endgames when six or fewer pieces are left on the chessboard have been completely analyzed and all possible moves (N=all) have been represented in a massive table that when uncompressed fills more than a terabyte of data. This enables chess computers to play the endgame flawlessly. No human will ever be able to outplay the system. Cet exemple est également repris par Nate Silver dans the Signal and the Noise, Chess might be thought of as analogous to prediction. The players must process information—the position of the thirty-two pieces on the board and their possible moves. They use this information to devise strategies to place their opponent in checkmate. These strategies in essence represent different hypotheses about how to win the game. Whoever succeeds in that task had the better hypothesis. Chess is deterministic—there is no real element of luck involved. But the same is theoretically true of the weather (…).Our knowledge of both systems is subject to considerable imperfections. In weather, much of the problem is that our knowledge of the initial conditions is incomplete. Even though we have a very good idea of the rules by which the weather system behaves, we have incomplete information about the position of all the molecules that form clouds and rainstorms and hurricanes. Hence, the best we can do is to make probabilistic forecasts. In chess, we have both complete knowledge of the governing rules and perfect information—there are a finite number of chess pieces, and they’re right there in plain sight. But the game is still very difficult for us. Chess speaks to the constraints on our information-processing capabilities—and it might tell us something about the best strategies for making decisions despite them. The need for prediction arises not necessarily because the world itself is uncertain, but because understanding it fully is beyond our capacity. Both computer programs and human chess masters therefore rely on making simplifications to forecast the outcome of the game. We can think of these simplifications as “models,” but heuristics is the preferred term in the study of computer programming and human decision making. It comes from the same Greek root word from which we derive eureka.10 A heuristic approach to problem solving consists of employing rules of thumb when a deterministic solution to a problem is beyond our practical capacities Et cet exemple du jeu d’échec se retrouve sur la traduction automatique, entre autres. The degree to which more data trumps better algorithms has been powerfully demonstrated in the area of natural language processing: the way computers learn how to parse words as we use them in everyday speech. Around 2000, Microsoft researchers Michele Banko and Eric Brill were looking for a method to improve the grammar checker that is part of the company’s Word program. They weren’t sure whether it would be more useful to put their effort into improving existing algorithms, finding new techniques, or adding more sophisticated features. Before going down any of these paths, they decided to see what happened when they fed a lot more data into the existing methods. Most machine-learning algorithms relied on corpuses of text that totaled a million words or less. Banko and Brill took four common algorithms and fed in up to three orders of magnitude more data: 10 million words, then 100 million, and finally a billion words. The results were astounding. As more data went in, the performance of all four types of algorithms improved dramatically. In fact, a simple algorithm that was the worst performer with half a million words performed better than the others when it crunched a billion words. Its accuracy rate went from 75 percent to above 95 percent. Inversely, the algorithm that worked best with a little data performed the least well with larger amounts, though like the others it improved a lot, going from around 86 percent to about 94 percent accuracy. “These results suggest that we may want to reconsider the tradeoff between spending time and money on algorithm development versus spending it on corpus development,” Banko and Brill wrote in one of their research papers on the topic. Aussi, sur http://translate.google.com/, on peut avoir une traduction du portugais vers le français, sans que l’ordinateur ne “connaisse” aucune des deux langues. Il pourrait tout aussi bien traduire vers de l’elfique ou du Klingon, comme le sugérait Chris Anderson dans the end of theory: the data deluge makes the scientific method obsolete. 6. Les outils mathématiques du Big-Data Un problème classique, vu comme un problème trop complexe car aboutissant à des calculs impossibles à mener en un temps raisonnable (disons avant la fin du système solaire dans quelques milliards d’années) est le problème du voyageur de commerce. En 1962, Procter & Gamble ont offert un prix de 10,000$ à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis.

Comme le calcule William Cook, dans In Pursuit of the Traveling Salesman, avec 33 villes, il y aurait

131,565,418,466,846,765,083,609,006,080,000,000

chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Une des ébauches de solution est celle proposée par Euler, dans son problème des sept ponts de Königsberg, dans solutio problematis ad geometriam situs pertinentis.

Cette approche à l’aide de graphs a donné naissance à la notion de Topological Data Analysis (trois tomes sur le sujet ont été publiés récemment chez Springer).

Topological data analysis is a way of getting structured data out of unstructured data so that machine-learning algorithms can act more directly on it

Parmi les techniques utilisées, on va retrouver les méthodes d’échantillonnage. En 2007, Jim Gray (lauréat du prix Turing) affirmait (cette idée a été reprise dans The Fourth Paradigm, paru en 2009)

The world of science has changed, and there is no question about this. The new model is for the data to be captured by instruments or generated by simulations before being processed by software and for the resulting information or knowledge to be stored in computers. Scientists only get to look at their data fairly late in this pipeline. The techniques and technologies for such data-intensive science are so different that it is worth distinguishing data-intensive science from computational science as a new, fourth paradigm for scientific exploration

Sur la visualisation des réseaux, on pourra relire à ce sujet cartographier les réseaux de Martin Grandjean, qui montre un cas pratique sur des réseaux entre parlementaires et activités de lobbying, en Suisse

avec en prime une visualisation élégante

Les ordinateurs ont transformé les mathématiques, plus personnes ne remettra ce point en doute. Les mathématiciens ont longtemps cherché des solutions analytiques aux équations obtenues; et si l’on sait qu’il sera impossible d’avoir une forme analytique élégante, on sait se contenter, depuis plusieurs décennies, de solutions numériques (on pourra relire science in the age of computer simuulations d’Eric Winsberg). Dans computational complexity: a modern approach, Sanjeev Arora et Boaz Barak présentent de manière presque compréhensible le concept de preuves interactives (dans le chapitre 8), introduites en 1985 par Goldwasser, Micali et Rackoff d’une part, Babai d’autre part.

As an example for a probabilistic interactive proof system, consider the following scenario: Marla claims to Arthur that she can distinguish between the taste of Coke (Coca-Cola) and Pepsi. To verify this statement, Marla and Arthur repeat the following experiment 50 times: Marla turns her back to Arthur, as he places Coke in one unmarked cup and Pepsi in another, choosing randomly whether Coke will be in the cup on the left or on the right. Then Marla tastes both cups and states which one contained which drinks. While, regardless of her tasting abilities, Marla can answer correctly with probability 1/2 by a random guess, if she manages to answer correctly for all the 50 repetitions, Arthur can indeed be convinced that she can tell apart Pepsi and Coke.

On retrouve en quelque sorte qu’avec suffisamment de données, on peut tirer des conclusions robustes.

Dans un article passionnant (datant de 2000), intitulé the curses and blessings of dimensionality, David Donoho nous rappelle que la grande dimension n’est pas qu’un fléau (“curve of dimensionsality”) comme l’évoquait Richard Bellman en 1957,

all [problems due to high dimension] may be subsumed under the header the curse of dimensionality. Since it is a curse (…) there is no need to feel discouraged about the possibility of obtaining significant results despite it

En fait, la grande dimension, ce n’est pas si simple, et les difficultés ne vont pas forcément croître avec la dimension. Considérons par exemple le volume d’une sphère de rayon 1, en dimension $d$. On peut montrer que le volume est

$V(d)=\frac{\pi^{d/2}}{\Gamma(1+d/2)}$

Si on regarde cette fonction, effectivement, passer de la dimension 3 a la dimension 4 puis a la dimension 5 faire grossir la taille de notre boule. Mais passée la dimension 10, de manière surprenante, la dimension diminue… Amusant non ?

En plus on n’est pas n’importe où ! En fait, on a très peu de chances d’être au centre: on peu montrer assez simplement que si on tire des points au hasard dans la sphère unité, on a en fait de très très fortes chances d’être au bord de la sphère. Et le calcul est simple: la probabilité d’être a une distance supérieure a $\alpha\in[0,1]$ (et donc d’être vraiment proche du bord) si $\alpha\sim1$ est

$\mathbb{P}(\|X\|\in[\alpha,1]) = 1-\alpha^d \rightarrow1$

lorsque $d\rightarrow \infty$ (quelle que soit la valeur de $\alpha$). Bref, les grands espaces, c’est grand, mais c’est rempli de vide !

Il y a plus de 5 ans, lors des journées UseR! à Rennes, Trevor Hastie avait fait un exposé brillant qui m’avait ouvert les yeux sur des problèmes que j’ignorais : les très grosses matrices sont souvent très vide, en évoquant l’exemple de la prévision pour Netflix ou Amazon (ce dernier exemple est repris dans big data, a revolution that will transform how we live, work and think)

Dans les années 60, John Tukey s’est battu pour que l’analyse de données (Data Analysis) devienne une discipline reconnue, distincte de la statistique mathématique. En 1962, il prononçait un discours au congrés annuel de l’IMS (Institute of Mathematical Statistics) intitulé the future of data analysis, qui servira de base à son ouvrage paru 15 ans plus tard, Explanatory Data Analysis. Dans ce discours, il quitte le monde des preuves mathématiques pour explorer les données. Comme le notera Howard Wainer,

he legitimized that, because he wasn’t doing it because he wasn’t good at math; he was doing it because it was the right thing to do

Comme le notait très justement John Tukey dans the technical tools of statistics (via http://cm.bell-labs.com/…)

Today, software and hardware together provide far more powerful factories than most statisticians realize, factories that many of today’s most able young people find exciting and worth learning about on their own. Their interest can help us greatly, if statistics starts to make much more nearly adequate use of the computer. However, if we fail to expand our uses, their interest in computers can cost us many of our best recruits, and set us back many years.

Car le Big Data, c’est résolument autre chose (par rapport à ce qui se faisait traditionnellement en statistique). Tout d’abord, les simulations se font autour d’un modèle mathématique: on est face à un espace immense, mais il y a un modèle. Alors que le Big Data repose davantage sur un problème d’exploration face à une forêt vierge. Dans le premier cas, on continue à progresser de manière déductive, alors que le Big Data force à avancer par induction. Pour reprendre les termes de Cosma Shalizi dans data mining

Data mining, more stuffily ‘knowledge discovery in databases’, is the art of finding and extracting useful patterns in very large collections of data. It’s not quite the same as machine learning, because, while it certainly uses ML techniques, the aim is to directly guide action (praxis!), rather than to develop a technology and theory of induction. In some ways, in fact, it’s closer to what statistics calls ‘exploratory data analysis’, though with certain advantages and limitations that come from having really big data to explore.

7. “With great power comes great responsibility

Tout le monde le dit, le Big Data redonne le pouvoir aux statisticiens. Pour reprendre l’analyse de Mike Loukides dans what is data science,

According to Mike Driscoll (@medriscoll), statistics is the ‘grammar of data science.’ It is crucial to ‘making data speak coherently.’. But it takes statistics to know whether this difference is significant, or just a random fluctuation.

C’est ce que Nate Silver appelle the Signal and the Noise (pour reprendre le titre de son livre).

Data science isn’t just about the existence of data, or making guesses about what that data might mean; it’s about testing hypotheses and making sure that the conclusions you’re drawing from the data are valid.

Cela dit raconter une histoire (et donc extraire les informations importantes) à partir d’observations en (très) grande dimension est toujours un exercice délicat. Et souvent subjectif.

8. Les applications importantes du Big-Data

Il existe des milliers d’applications concrètes, dans la vie de tous les jours (ou presque). Un exemple important peut être l’imagerie. Il s’agit d’un problème de Big-Data au même sens que les jeux d’échecs. Par exemple, pour reprendre un graphique tiré de Cireşan et al. (2013), sur la détection de mitoses dans des cancers du sein (par des réseaux de neurones), on cherche des tâches sombres atypiques. Avec une fréquence, et une distribution que l’on retrouve chez des personnes malades, et pas chez des personnes saines,

En tant que statisticien, le big data, c’est un porte-ouverte fabuleuse. Pour reprendre le diagramme de Venn proposé par Vincent Granville,

C’est plus ni moins ce qu’affirmait Mike Loukides dans what is data science,

What differentiates data science from statistics is that data science is a holistic approach. We’re increasingly finding data in the wild, and data scientists are involved with gathering data, massaging it into a tractable form, making it tell its story, and presenting that story to others.

Et je ne peux que confirmer ! Lorsqu’auparavant, dans une soirée, je disais que j’étais mathématicien, ou statisticien, j’étais constamment vu comme un nerd (et ce n’était probablement pas infondé). Maintenant que je suis data scientist, expert en big data, je me fais inviter par des jolies filles dans des bars en tant que VIP (des bars des sciences, certes, mais c’est un début !).

• le Big Data avec les yeux d’un économiste

Quand on pense Big Data, on pense aussi au commerce en ligne, et aux applications interactives (en tous les cas qui doivent fournir des réponses beaucoup plus vite qu’une analyse statistique poussée, menée sur plusieurs semaines, voir plusieurs mois).Pour reprendre l’analyse de Dave Rich et Jeanne Harris dans why predictive analytics is a game-changer, how companies use real-time data to plan for the future,

In simple terms analytics means using quantitative methods to derive insights from data, and then drawing on those insights to shape business decisions and, ultimately, improve business performance. Thus predictive analytics is emerging as a game-changer. Instead of looking backward to analyze ‘what happened?’ predictive analytics help executives answer ‘What’s next?’ and ‘What should we do about it?’

1. Enjeux financiers du Big Data

En février 2010, The Economist, dans l’article data, data everywhere, évaluait que la vente de logiciels spécialisés dans les notions de “data management” et autres “analytics” représentait un enjeux financier colossal, avec des montants de l’ordre de 100 milliards de dollars, avec un taux de croissance de l’ordre de 10% par an,

The business of information management—helping organisations to make sense of their proliferating data—is growing by leaps and bounds. In recent years Oracle, IBM, Microsoft and SAP between them have spent more than $15 billion on buying software firms specialising in data management and analytics. This industry is estimated to be worth more than$100 billion and growing at almost 10% a year, roughly twice as fast as the software business as a whole.

Et le travail des data engineers va au delà (en quelque sorte) du travail traditionnel du statisticien. Le titre complet de l’article what is data science? de Mike Loukides (publié en juin 2010) était the future belongs to the companies and people that turn data into products. Pour une analyse plus poussées de liens entre l’économie et le Big Data, on pourra relire the data revolution and economic analysis par Liran Einav et Jonathan Levin, ou seven big data trends for 2014 pour des aspects plus business.

2. Marketing et publicité

Quand on pense Big Data, on pense aussi aux publicités en ligne, ciblées en fonction de mon profil (des mots clés que j’ai pu taper, et des sites que j’ai pu visiter auparavant). Dans pêcher le client dans une baignoire, Ariane Krol et Jacques Nantel écrivent

Bienvenue dans le nouveau monde du marketing personnalisé. Un monde qui veut votre bien… et qui fera tout pour l’obtenir. Au début des années 1980, décliner, à coups de sondages auprès des consommateurs, ses stratégies en fonction des « segments » visés — ménagères de plus de 50 ans, professions libérales de moins de 35 ans ayant des revenus supérieurs à 210 000 francs et jouant au tennis au moins deux fois par mois, etc. — était le fin du fin. Les spécialistes du secteur pêchaient en quelque sorte au filet dérivant, après que leur sonar avait signalé un banc de poissons de la bonne espèce. Aujourd’hui, ce ne sont plus les jeunes de moins de 35 ans ou tout autre segment qui les intéressent : c’est vous. Plus de filet, plus de sonar : la pêche se fait dans une baignoire.

En 1998, le moteur de recherche Google se contentait de suggérer des sites (comme tout moteur de recherche)

Aujourd’hui, Google va beaucoup plus loin. En 2013, sur 14 Mds$de chiffre d’affaire de la compagnie (au premier trimestre), 85% provenait de la vente d’espace publicitaire (8.640 Mds$ pour Google AdWords, et 3.262 Mds$pour Google AdSense). En tant qu’économiste (théorique), ce marché pose des problèmes fascinants, de recherche d’équilibre dans un jeu non coopératif, répété, en temps continu (on pourra relire à ce sujet computer scientists optimize innovative ad auction, de Sara Robinson, ou toward an integrated framework for automateddevelopment and optimization of online advertising campaigns, de Stamatina Thomaidou, Michaelis Vazirgiannis et Kyriakos Liakopoulos). Mais cela pose aussi des difficultés aux entreprises qui souhaitent acheter de la publicité lorsqu’un mot clé est tapé. Comme le rappelle Stephen Baker dans une petite anecdote dans l’introduction de the Numerati, on peut trouver des corrélations intéressantes en regardant les données, en particulier pour essayer de comprendre quand les gens louent de voitures de location (afin de lisser une bannière publicitaire): il évoque ainsi une corrélation avec les billets d’avion, les couronnes funéraires, et les films romantiques. ‘What is it about romantic-movie lovers?’ Morgan asks, as we sit in his New York office on a darkening summer afternoon. The advertising entrepreneur is flush with details about our ramblings online. He can trace the patterns of our migrations, as if we were swallows or humpback whales, while we move from site to site. Recently he’s become intrigued by the people who click most often on an ad for car rentals. Among them, the largest group had paid a visit to online obituary listings. That makes sense, he says, over the patter of rain against the windows. ‘Someone dies, so you fly to the funeral and rent a car.’ But it’s the second-largest group that has Morgan scratching his head. Romantic-movie lovers. For some reason Morgan can’t fathom, loads of them seem drawn to a banner ad for Alamo Rent A Car. (…) I ask him about the correlation he told me about earlier, the one between romantic-movie fans and Alamo Rent A Car. It takes a moment for him to recall it. ‘Oh yeah. They were off the charts.’ Did his researchers, I ask, ever come up with an explanation for it. He nods. “It had to do with weekends. It was Alamo ads promoting ‘escapes’ that attracted the attention of these web surfers, he says. The romantic-movie fans booked leisure rentals, largely for weekend getaways. Perhaps they wanted to act out the kind of scenes that drew them to the cinema. Banners for weekday rentals apparently left them cold. Il est important d’avoir accès à des gros volumes de données afin d’acheter des espaces publicitaires, sur les bons mots clés. Mais en plus, il faut intégrer la composante dynamique : acheter un espace publicitaire, c’est un premier pas. Mais il peut être intéressant de réapparaitre régulièrement, après, en tenant compte de la courbe d’oubli (dans l’idée des travaux d’Hermann Ebbinghaus, par exemple) Il ne faut pas s’étonner des voir réapparaître régulièrement des espaces publicitaires sur les crèmes solaires sur plusieurs sites visités, si on a demandé le prix du billet Montréal-Miami chez un vendeur en ligne.. • la perte de contrôle, Big Data et Big Brother Revenons maintenant au second terme du titre de la soirée, Big Brother. 1. Statistique et collecte de données Avant de commencer, juste une petite réflexion en passant. Comme le note Alain Desrosières, pour définir le mot statistiqueon remonte traditionnellement à Abriß der neuen Staatswissenschaft der vornehmen Europäischen Reiche und Republiken, publié en 1749 par Gottfried Achenwall: dans cet ouvrage, le mot Statistik (ou Staatwissenschaft ou Kameralwissenschaft) est souvent associé à des données démographiques, de taux de fécondité, ou de nombre de conscrits. On notera que ces deux quantités sont associées à la puissance et la richesse d’un état, au XVIIIème siècle en tous cas. En ce sens, je rejoindrais Daniel Solove notait que la bonne comparaison n’était peut-être pas 1984 d’Orwell, mais davantage le Procès de Kafka (le film d’Orson Welles est d’ailleurs en visionnage libre sur http://openculture.com/). Dans le Procès, il n’y a pas explicitement de surveillance, mais juste une bureaucratie qui stocke des données, toutes sortes de données. Le héros se sent dépossédé quand on l’accuse, et impuissant. Comme quand on se fait refuser un retrait dans une machine ATM au Mexique, parce qu’un algorithme a jugé qu’il était improbable qu’un jeudi, à 14:17, vous soyez sur une plage de la Riviera Maya. Comme on l’a rappelé au tout début, cette statistique purement descriptive (on se contentait de collecter des données) a évolué vers une statistique inférentielle par la suite. Comme le rappellent Maurice Kendall en 1942 et Victor Hilts en 1978, la Royal Statistical Society, fondée en 1834 avait pour devise Aliis Exterendum, expliquant clairement qu’il n’était dans dans les attributions d’un statisticien d’interpréter les données (on comprend qu’en 1857, les membres de cette éminente société savante l’aient changé). Cette transition vers l’inférence statistique (ce qui est aujourd’hui appelé analytics par certains) a révolutionné la science, comme le notent Ronald Nelson, Mats Pettersson et Örjan Carlborg dans a century after Fisher: time for a new paradigm in quantitative geneticsOn pourra aussi relire Maurice Kendall, parlant de l’expansion de la statistique, They have already overrun every branch of science with a rapidity of conquest rivalled only by Attila, Mohammed, and the Colorado beetle Mais mon point ici est que collecter des données n’est pas récent, et pas uniquement lié au développement des ordinateurs, et de l’internet, loin de là. Sur le sectre de Big Brother, on pourra relire à ce sujet l’article paru dans Newsweek en juillet 1970, 2. De la granularité des statistiques La Statistik pose un problème profond de granularité: à quel niveau de détail peut-on avoir accès à des données? i) il est possible d’avoir accès à des données très fines, mais par des techniques d’enquêtes (consommation de riz par mois, temps passé devant la télévision, etc), et donc non-exhaustives ii) de l’autre, quand des données exhaustives sont connues, il est délicat de les communiquer (on pensera au patrimoine financier, et tout autre information connues par les centres des impôts et des taxes). La réponse classique est alors d’agréger les données par région géographique ou par classe d’âge. Un contre-exemple (suffisamment rare pour être repris) est peut-être le cas des données du Census américain (repris par Bill Rankin) qui a permis d’obtenir des informations sur le revenu et des informations raciales à des degrés de granularité très très fin iii) entre les deux, reposant sur le mouvement de libération de données (open-data), ces organismes sont désormais tenu de divulguer des informations, mais souvent, les données sont agrégées, afin d’éviter tout soucis de confidentialité (les textes légaux sont plus anciens, comme le Freedom of Information Act, datant de 1966). Mais l’agrégation pose des soucis d’interprétation, comme le rappellent toutes les études sur la notion d’ecological fallacy. 3. Les données individuelles On peut raconter beaucoup de choses à partir de données individuelles. En octobre dernier, dans I challenged hackers to investigate me and what they found out is chilling, Adam Penenberg tente ainsi de reproduire la filature de Sophie Calle mais uniquement à l’aide des traces numériques. Dans un autre style, en mars 2012, Stephen Wolfram s’est amusé à retracer 10 ans de sa vie à partir du détail de ses communications téléphoniques, des courriels envoyés et reçus, etc. Dans the personal analytics of my life il revient sur le rythme de ses journées, A quoi cette information peut-elle servir? Je me souviens en mai dernier, je passais une visite médicale pour l’immigration, et le médecin m’a demandé “et à part ça, vous dormez bien?”. Il ne s’agit pas d’une vraie question piège, juste d’une question de routine. Je suppose. Et pourtant, ça m’a déstabilisé. Car ça dépend. Je dors moins que mes enfants. Moins que je devrais, probablement. J’ai des périodes d’insomnies, souvent quand les enfants font des cauchemars. Mais comment un médecin peut-il travaillé avec des données aussi pourries? C’est un peu comme quand mon fils m’explique qu’il a mal au ventre. Ça fait pas très mal, mais un peu. Et par moments. Pour savoir si je dors bien, il pourrait me mettre des électrodes, et me suivre pendant 2 mois. Si j’avais un téléphone ‘intelligent‘, je pourrais aussi avoir un traqueur qui enregistre des informations pendant mon sommeil Cela dit, je peux aussi regarder mes logs de connexion, via mes envois de courriels, ou via Twitter (j’avais réellement fait cette étude il y a quelques mois). 4. Le terrorisme, la NSA et les données individuelles A en lire certains médias, un moment clé dans l’histoire du Big Data aura été les révélations d’Edward Snowden, qui a rendu public, pendant l’été 2013 (via le Guardian et le Washington Post), un certain nombre d’information considérées comme secrètes de la NSA, concernant par exemple la captation des métadonnées d’appels téléphoniques aux États-Unis, ainsi que les systèmes d’écoute sur internet (le programme de surveillance PRISM, mais aussi le programme de surveillance Tempora du gouvernement britannique). Cela dit, ces révélations ont surtout permis de confirmer, par des documents officiels, certaines informations connues de beaucoup de monde. En 2008, dans the Numerati, Stephen Baker s’étonnait déjà du fait qu’à la NSA, la statistique avait remplacé la cryptographie. Depuis 2001, il ne s’agit plus de décrypter les messages coder, mais d’identifier des comportements atypiques, qui pourraient être associées à des activités terroristes. What sorts of data would fuel the hunt for terrorists? Practically anything the government could get its hands on. In the years following 9/11/1 , the government spent more than$1 billion to merge its enormous databases, including those of the FBI and the CIA. This would give data miners a single unified resource. But that wasn’t all. They would also trawl, oceans of consumer and demographic details, airline records and hotel receipts, along with videos, photos, and millions of hours of international phone and Internet traffic harvested by the NSA. This trove matched anything that the Web giants Yahoo and Google were grappling with. In May 2006, news surfaced that the NSA was secretly extending its nets even fur­ther. USA Today reported that major phone companies had delivered hundreds of billions of phone records to the govern­ment. These provided details on who was calling whom, from where, for how long, and whether the call was forwarded. Were the NSA staff also listening in on the calls and reading the e-mails? There was no telling.

On essaye ici de quitter l’homme moyen d’Adolphe Quetelet pour détecter les déviations par rapport à une norme.

For the re­searchers to pick out these oudiers, they must first figure out what’s ‘normal’. Picture our society on a big piece of poster­ board. At first glance it looks entirely blue, monochromatic. But step closer, and you’ll see tiny dots and strings of red. That background of blue represents boring, law-abiding (for the most part) us. Our only function on this display is to bring into relief the bits of red. Those are the suspected terrorists. (…) Are there times, I ask, when you just have too much data? When it gets in the way and confuses things? He seems taken aback by this line of questioning. ‘More data is always better,’ he says.

5. Big Data est mon ami…

Dans l’histoire de la NSA, le Big Data fait peur car on ne sait pas trop qui a accès à toutes ces informations. Et en l’occurrence, le problème visé (détecter un terroriste) est un problème plus complexe que détecter des tendances de consommation. Si j’achète un disque sur amazon, j’ai moins tendance à essayer d’effacer les traces que doivent probablement le faire des terroristes s’ils commandent en ligne des grosses quantités de nitrate d’ammonium (mais je suppose que ça doit dépendre de quel disque j’achète en ligne, on a tous nos petites faiblesses et nos petites lâchetés…). La NSA (et tous les hommes en noirs, de manière générale) font peur. Mais il existe des cas où les utilisateurs de Big Data sont (relativement) bien identifiés, et le but est lui aussi clair.

Par exemple, Google avait montré (dans une étude sur la grippe) que l’analyse des mots clés tapés dans le moteur de recherche permettait de détecter des épidémies avant que les médecins ne lancent l’alerte.

A titre individuel, je n’aimerais pas que l’on rende public toutes les recherches que je peux faire en ligne. Mais collectivement, on imagine très bien, sur cet exemple que l’agrégation d’autant d’information permet de détecter d’éventuelles tendance, malgré tout le bruit. Un autre exemple connu est celui de l’utilisation des données de Twitter pour détecter une épidémie de choléra suite au tremblement de terre en Haiti, deux semaines avant sa reconnaissance officielle. Plus proche de nous, dans how New York’s Fire Department uses data mining, Elizabeth Dwoskin. nous explique comment le NYFD utilise le Big Data pour détecter des débuts d’incendies.

Le crowdsourcing est un des gros intérêt associé au Big Data. Pour reprendre un exemple emprunté à Dick Kasperowski (et évoqué au début de ce billet) si le crowdsourcing avait existé au XVIème siècle (l’expression black swan, du poète Decimus Iunius Iuvenalis, dit Juvenal, “rara avis in terris nigroque simillima cygno” était utilisée à Londres pour décrire un évènement impossible), le débat sur les blacks swans n’aurait pas existé, car on peut espérer que quelqu’un se soit manifester pour dire que les cygnes noirs existent (avant que Willem de Vlamingh ne les découvre, en 1697) . De manière générale le crowdsourcing permet d’éviter les biais cognitifs.

On pourrait continuer à l’infini je pense… et je veux en garder un peu pour la semaine prochaine. Et si ça ne suffit pas (je pense aux aspects de visualisation des données), je reviendrais une nouvelle fois au Cœur des Sciences !

# ROC curves and classification

To get back to a question asked after the last course (still on non-life insurance), I will spend some time to discuss ROC curve construction, and interpretation. Consider the dataset we’ve been using last week,

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)

The first step is to get a model. For instance, a logistic regression, where some factors were merged together,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db\$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)

From this model, we can predict a probability, not a $\{0,1\}$ variable,

> S=predict(reg,type="response")

Let $\widehat{S}$ denote this variable (actually, we can use the score, or the predicted probability, it will not change the construction of our ROC curve). What if we really want to predict a $\{0,1\}$ variable. As we usually do in decision theory. The idea is to consider a threshold , so that

• if , then  will be $1$, or “positive” (using a standard terminology)
• si , then  will be $0$, or “negative

Then we derive a contingency table, or a confusion matrix

 observed value predicted value “positive“ “négative“ “positive“ TP FP “négative“ FN TN

where TP are the so-called true positive, TN  the true negative, FP are the false positive (or type I error) and FN are the false negative (type II errors). We can get that contingency table for a given threshold

> roc.curve=function(s,print=FALSE){
+ Ps=(S>s)*1
+ FP=sum((Ps==1)*(Y==0))/sum(Y==0)
+ TP=sum((Ps==1)*(Y==1))/sum(Y==1)
+ if(print==TRUE){
+ print(table(Observed=Y,Predicted=Ps))
+ }
+ vect=c(FP,TP)
+ names(vect)=c("FPR","TPR")
+ return(vect)
+ }
> threshold = 0.5
> roc.curve(threshold,print=TRUE)
Predicted
Observed   0   1
0   5 231
1  19 745
FPR       TPR
0.9788136 0.9751309

Here, we also compute the false positive rates, and the true positive rates,

• TPR = TP / P = TP / (TP + FN) also called sentivity, defined as the rate of true positive: probability to be predicted positve, given that someone is positive (true positive rate)
• FPR = FP / N = FP / (FP + TN) is the rate of false positive: probability to be predicted positve, given that someone is negative (false positive rate)

The ROC curve is then obtained using severall values for the threshold. For convenience, define

> ROC.curve=Vectorize(roc.curve)

First, we can plot $(\widehat{S}_i,Y_i)$ (a standard predicted versus observed graph), and visualize true and false positive and negative, using simple colors

> I=(((S>threshold)&(Y==0))|((S<=threshold)&(Y==1)))
> plot(S,Y,col=c("red","blue")[I+1],pch=19,cex=.7,,xlab="",ylab="")
> abline(v=threshold,col="gray")

And for the ROC curve, simply use

> M.ROC=ROC.curve(seq(0,1,by=.01))
> plot(M.ROC[1,],M.ROC[2,],col="grey",lwd=2,type="l")

This is the ROC curve. Now, to see why it can be interesting, we need a second model. Consider for instance a classification tree

> library(tree)
> ctr <- tree(Y~X1+X2+X3bis,data=db)
> plot(ctr)
> text(ctr)

To plot the ROC curve, we just need to use the prediction obtained using this second model,

> S=predict(ctr)

All the code described above can be used. Again, we can plot $(\widehat{S}_i,Y_i)$ (observe that we have 5 possible values for $\widehat{S}_i$, which makes sense since we do have 5 leaves on our tree). Then, we can plot the ROC curve,

An interesting idea can be to plot the two ROC curves on the same graph, in order to compare the two models

> plot(M.ROC[1,],M.ROC[2,],type="l")
> lines(M.ROC.tree[1,],M.ROC.tree[2,],type="l",col="grey",lwd=2)

The most difficult part is to get a proper interpretation. The tree is not predicting well in the lower part of the curve. This concerns people with a very high predicted probability. If our interest is more on those with a probability lower than 90%, then, we have to admit that the tree is doing a good job, since the ROC curve is always higher, comparer with the logistic regression.

# Think academic journals look the same ? Well, some do…

We have seen yesterday that finding an optimal strategy to publish is not that simple. And actually, it can be even more difficult in the case the journal rejects the paper (not because it is not correct, but because “it does not fit” with the standards, the quality of the journal, the audience, the editor’s mood, or whatever). The author has basically two choices,

• forget about the article and move to something else (e.g. start a blog where he/she will be the author and the editor)
• pretend that the article is worth publishing and then try to find another journal with similar interests

But this last choice is not that easy, since sometimes the author think that this journal was indeed the one that should publish it (e.g. all the articles on the subject have been published in that journal).
So I was wondering if there were clusters of journals, i.e. journals that publish almost the same kind of articles (so that next time one of my paper is rejected by the editor, I just go to for some journal in the same cluster).
So what I did is extremely simple: I looked at articles titles and looked for correlations between words frequency (I could have done that in key words, but I am not a big fan of those key words). I looked at 35 journals (that are somehow related to my areas of interest) and looked at titles of all articles published over the last 20 years. Then I kept the top 1000 of words, and I removed standard short words (“a“, “the“, “is“, etc). Actually, my top words looks like

"models" "model" "data" "estimation" "analysis" "time"
"processes" "risk" "random" "stochastic" "regression"
"market" "approach" "optimal" "based" "information"
"evidence" "linear" "games" "bayesian" "theory" "effects"
"distribution" "multivariate" "tests" "markets" "markov"
"equilibrium" "dynamic" "process" "distributions"
"application" "stock" "likelihood"

Then, I ran a principal component analysis on my dataset (containing 960 variables – here words – and 35 observations – here journal names).

library("FactoMineR")
res.pca = PCA(MATRICE, scale.unit=TRUE, ncp=5,
graph=FALSE)
plot.PCA(res.pca, axes=c(1, 2), choix="ind")

The projection of the journals on the first two axis looks like that

Here, we can clearly observe some clusters : on the up-left Journal of Finance and Journal of Banking and Finance (say financial journals) on the top-right Biometrika, Biometrics, Computational Statistics and Data Analysis and Journal of Econometrics (JASA is not far away, i.e. applied statistics journal). And below, on the right, Stochastic Processes and their Applications, Annals of Applied Probability, Journal of Applied Probability, Annals of Probability, Proceedings of AMS and Topology and Applications (ie more theoretical journal).
Note that the projection is rather robust: if I consider my first 200 words, the graph is the same

In order to go further in the interpretation, we can also plot variables, i.e. words from titles,

where we cannot distinguish anything. So if I just look at my top 30, here they are,

On top left we see market(s), risk or information; on top right analysis, effects, models or tests; while below we see Markov or process(es). And we can observe interesting facts: in finance in statistics, we talk about dynamics while in theoretical (mathematical) journal it is about processes.
But the goal was to find cluster, i.e. classes of journals that publish papers with similar titles.

DISTANCE = dist(MATRICE)
cah = hclust(DISTANCE)
plot(cah)

Here we have

If some classes a rather natural (Journal of Applied Proba. and Advances in Applied Proba.or Economic Theory, Journal of Economic Theory and Journal of Mathematical Economics) some strong correlation are not simple to understand, (e.g. Insurance: Mathematics and Economics and Management Science or Annals of Statistics and the Journal of Multivariate Analysis).
Again, it might be possible to spend hours on the graphs, but if I want – someday – to submit something to one of those journals, I guess I have to stop here, and move to something else…