Tag Archives: R

On the conjugate function

In the MAT7381 course (graduate course on regression models), we will talk about optimization, and a classical tool is the so-called conjugate. Given a function f:\mathbb{R}^p\to\mathbb{R} its conjugate is function f^{\star}:\mathbb{R}^p\to\mathbb{R} such that f^{\star}(\boldsymbol{y})=\max_{\boldsymbol{x}}\lbrace\boldsymbol{x}^\top\boldsymbol{y}-f(\boldsymbol{x})\rbraceso, long story short, f^{\star}(\boldsymbol{y}) is the maximum gap between the linear function \boldsymbol{x}^\top\boldsymbol{y} and f(\boldsymbol{x}).

Just to visualize, consider a simple parabolic function (in dimension 1) f(x)=x^2/2, then f^{\star}(\color{blue}{2}) is the maximum gap between the line x\mapsto\color{blue}{2}x and function f(x).

x = seq(-100,100,length=6001)
f = function(x) x^2/2
vf = Vectorize(f)(x)
fstar = function(y) max(y*x-vf)
vfstar = Vectorize(fstar)(x)

We can see it on the figure below.

viz = function(x0=1,YL=NA){
idx=which(abs(x)<=3) par(mfrow=c(1,2)) plot(x[idx],vf[idx],type="l",xlab="",ylab="",col="blue",lwd=2) abline(h=0,col="grey") abline(v=0,col="grey") idx2=which(x0*x>=vf)
if(is.na(YL)) YL=range(vfstar[idx])



In that case, we can actually compute f^{\star}, since f^{\star}(y)=\max_{x}\lbrace xy-f(x)\rbrace=\max_{x}\lbrace xy-x^2/2\rbraceThe first order condition is here x^{\star}=y and thusf^{\star}(y)=\max_{x}\lbrace xy-x^2/2\rbrace=\lbrace x^{\star}y-(x^{\star})^2/2\rbrace=\lbrace y^2-y^2/2\rbrace=y^2/2And actually, that can be related to two results. The first one is to observe that f(\boldsymbol{x})=\|\boldsymbol{x}\|_2^2/2 and in that case f^{\star}(\boldsymbol{y})=\|\boldsymbol{y}\|_2^2/2 from the following general result : if f(\boldsymbol{x})=\|\boldsymbol{x}\|_p^p/p with p>1, where \|\cdot\|_p denotes the standard \ell_p norm, then f^{\star}(\boldsymbol{y})=\|\boldsymbol{y}\|_q^q/q where\frac{1}{p}+\frac{1}{q}=1The second one is the conjugate of a quadratic function. More specifically if f(\boldsymbol{x})=\boldsymbol{x}^{\top}\boldsymbol{Q}\boldsymbol{x}/2 for some definite positive matrix \boldsymbol{Q}f^{\star}(\boldsymbol{y})=\boldsymbol{y}^{\top}\boldsymbol{Q}^{-1}\boldsymbol{y}/2. In our case, it was a univariate problem with \boldsymbol{Q}=1.

For the conjugate of the \ell_p norm, we can use the following code to visualize it

p = 3
f = function(x) abs(x)^p/p
vf = Vectorize(f)(x)
fstar = function(y) max(y*x-vf)
vfstar = Vectorize(fstar)(x)


p = 1.1
f = function(x) abs(x)^p/p
vf = Vectorize(f)(x)
fstar = function(y) max(y*x-vf)
vfstar = Vectorize(fstar)(x)
viz(1, YL=c(0,10))

Actually, in that case, we almost visualize that if f(x)=|x| then\displaystyle{f^{\star}\left(y\right)={\begin{cases}0,&\left|y\right|\leq 1\\\infty ,&\left|y\right|>1.\end{cases}}}

To conclude, another popular case, f(x)=\exp(x) then{\displaystyle f^{\star}\left(y\right)={\begin{cases}y\log(y)-y,&y>0\\0,&y=0\\\infty ,&y<0.\end{cases}}}We can visualize that case below

f = function(x) exp(x)
vf = Vectorize(f)(x)
fstar = function(y) max(y*x-vf)
vfstar = Vectorize(fstar)(x)

Combining automatically factor levels with trees

Last year, in a post, I discussed how to merge levels of factor variables, using combinatorial techniques (it was for my STT5100 cours, and trees are not in the syllabus), with an extension on trees at the end of the post.

consider the following (simulated dataset)

'data.frame':	200 obs. of  3 variables:
 $ y : num  1.345 1.863 1.946 2.481 0.765 ...
 $ x1: num  0.266 0.372 0.573 0.908 0.202 ...
 $ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ...
 A  B  C  D  E  F  G  H  I  J 
11 12 23 34 23 36 12 32  3 14

Just by looking at the data (see the previous post), we could easily get the feeling that 10 levels was too much.

Following my post, Przemyslaw sent a comment suggesting to use


It is indeed a nice package (unless you have really really big datasets with a lot of categories in your factor variables – as I experienced recently), and you can get great graphs

MF = mergeFactors(response = b$y, 
             factor = b$x2, 
             family = "gaussian")

Here is suggests to create three categories. Recall that with student t-tests (changing the reference), we got

Another interesting package, by Piro Polo, is


To use it, we simply call the following function, and we transform automatically our dataset : the continuous variables remain unchanged, and (possibly) categories of categorical variables are merged

b.bins = tree.bins(data=b, y=y)
Classes ‘data.table’ and 'data.frame':	200 obs. of  3 variables:
 $ y : num  1.345 1.863 1.946 2.481 0.765 ...
 $ x1: num  0.266 0.372 0.573 0.908 0.202 ...
 $ x2: chr  "Group.4" "Group.4" "Group.4" "Group.4" ...
 - attr(*, ".internal.selfref")= 

Group.1 Group.2 Group.3 Group.4 
     23      35      26     116

here in four groups. To get the correspondance, use

tree.bins(data=b, y=y, return = "lkup.list")
   x2 Categories
1   E    Group.1
2   G    Group.2
3   C    Group.2
4   B    Group.3
5   J    Group.3
6   I    Group.4
7   A    Group.4
8   H    Group.4
9   F    Group.4
10  D    Group.4

(we have a list with one element, one dataframe, since there is only one factor variable). Cool, isn’t it ? I miss Przemyslaw’s plot, but this is rather quick, and efficient..


On leverage

Last week, in our STT5100 (applied linear models) class, I’ve introduce the hat matrix, and the notion of leverage. In a classical regression model, \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta} (in a matrix form), the ordinary least square estimator of parameter \boldsymbol{\beta} is \widehat{\boldsymbol{\beta}}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y}The prediction can then be written\widehat{\boldsymbol{y}}=\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\underbrace{\color{blue}{\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}}_{\color{blue}{\boldsymbol{H}}}\boldsymbol{y}where \color{blue}{\boldsymbol{H}} is called the hat matrix.

The matrix is idempotent, i.e. \boldsymbol{H}^2={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\textcolor{grey}{\boldsymbol{X}^\top{\boldsymbol{X}}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}}\boldsymbol{X}^\top}={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}=\boldsymbol{H}so it can be interpreted as a projection matrix. Furthermore, since\boldsymbol{H}\boldsymbol{X}=\boldsymbol{X} (just do the maths), the projection is on a subspace that contains all the linear combinations of columns of \boldsymbol{X}. One can also observe that \mathbb{I}-\boldsymbol{H} is also a projection matrix. And we can write\boldsymbol{y}=\underbrace{\boldsymbol{H}\boldsymbol{y}}_{\widehat{\boldsymbol{y}}}+\underbrace{(\mathbb{I}-\boldsymbol{H})\boldsymbol{y}}_{\widehat{\boldsymbol{\varepsilon}}}where \widehat{\boldsymbol{y}} is the orthogonal projection of \boldsymbol{y} on the (linear) space of linear combinations of columns of \boldsymbol{X}, and \widehat{\boldsymbol{y}}\perp\widehat{\boldsymbol{\varepsilon}}, which gives the classical interpretation of residuals, being unpredictible (at least with a linear model using variables \boldsymbol{X}).

Let’s move a bit faster now (we’ve seen many other properties last week), and consider elements on the diagonal of matrix \boldsymbol{H}. Recall that we have

so entry \boldsymbol{H}_{i,i} is a measure of the influence of entry \boldsymbol{y}_i on its prediction latex]\widehat{\boldsymbol{y}}_i[/latex].

We have seen that\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}(\boldsymbol{H})=\text{trace}(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top)which can be written\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}\boldsymbol{X}^\top(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1})=\text{trace}(\mathbb{I})=pwhere classically p=k+1, where k is the number of explanatory variables. Further, since \boldsymbol{H} is idempotent, we can write (from \boldsymbol{H}=\boldsymbol{H}^2) that\boldsymbol{H}_{i,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}\boldsymbol{H}_{j,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}^2One the one hand, since the second term is positive \boldsymbol{H}_{i,i}\geq\boldsymbol{H}_{i,i}^2, i.e. 1\geq\boldsymbol{H}_{i,i}. And since both terms are positive, then \boldsymbol{H}_{i,i}\in[0,1]. And there was a question in the course on the sharpeness of the bounds.

Using Anscombe’s dataset, we’ve seen that it was possible to get a leverage of 1. Using something rather similar

df = data.frame(x = c(rep(1,10),6), y = c(1:10,8))

we obtain

model = lm(y~x,data=df)
H = lm.influence(model)$hat

The very last observation, the one one the right, is here extremely influencial : if we remove it, the model is completely different ! And here, we reach the upper bound, \boldsymbol{H}_{11,11}=1. Observe that all other points are equally influencial, and because on the constraint on the trace of the matrix, \boldsymbol{H}_{i,i}=1/10 when i\in\{1,2,\cdots,10\}.

Now, what about the lower bound ? In order to have some sort of “non-influencial” observations, consider the two following case.

  • the case where one observation (below the first one) is such that \widehat{\boldsymbol{y}}_{i}=\boldsymbol{y}_{i} (perfect prediction)
  • the case where one observation (below the tenth one) is such that \boldsymbol{x}_{i}=\overline{\boldsymbol{x}} and \boldsymbol{y}_{i}=\overline{\boldsymbol{y}} (from the first order condition – or normal equation), the fitted regression line always go through point (\overline{\boldsymbol{x}},\overline{\boldsymbol{y}})

Let us move two observations from our dataset,

[1] 1.8
df = data.frame(x = c(4,rep(1,8),6,1.8),
y = c(predict(model,newdata=data.frame(x=4)),

We now have

If we compute the leverages, we obtain

model = lm(y~x,data=df)
H = lm.influence(model)$hat

so, for the first observation, its leverage actually increased (the blue part), and for the tenth one, we have the lowest influence, but it is not zero. Is it possible to reach zero ?

Here, observe that for the tenth observation, \boldsymbol{H}_{i,i}=1/n. And actually, that’s the best we can do… We can prove that, in the case of a simple regression (as above)\boldsymbol{H}_{i,i}=\frac{1}{n}+\frac{(x_i-\overline{x})^2}{n\text{Var}(x)}which is minimum when x_i=\overline{x}, and then \boldsymbol{H}_{i,i}=1/n, otherwise \boldsymbol{H}_{i,i}>1/n. And this property is also valid in a multiple regression (as soon as an intercept is included in the regression – which should always be the case). To prove that result, let \tilde{\boldsymbol{X}} denote the matrix of centered variables \boldsymbol{X}, then we can prove that \boldsymbol{H}_{i,i}=\frac{1}{n}+\big[\tilde{\boldsymbol{X}}(\tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}})^{-1}\tilde{\boldsymbol{X}}^\top\big]_{i,i}(which is basically a matrix version of the previous equation).

I can maybe add another comment on Anscombe’s data. We’ve seen that on the right that we did reach 1. But I did not prove it. One way to prove it is actually to focus on the remaining n-1 points, on the left. Those have all the same x values. We can prove that if \boldsymbol{X}_{i_1}=\boldsymbol{X}_{i_2}, then \boldsymbol{H}_{i_1,i_2}=\boldsymbol{X}_{i_1}^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_{i_2}=\boldsymbol{H}_{i_1,i_1}hence, using the relationship obtained since the hat matrix is idempotent\boldsymbol{H}_{i_1,i_1}=2\boldsymbol{H}_{i_1,i_1}^2+\sum_{j\notin\{i_1,i_2\}}\boldsymbol{H}_{i_1,j}^2thus, we now have\boldsymbol{H}_{i_1,i_1}\big(1-2\boldsymbol{H}_{i_1,i_1}\big)>0i.e. \boldsymbol{H}_{i_1,i_1}\in[0,1/2], where the upper bound becomes 1/(n-1) “duplicates”. So for n-1 \boldsymbol{H}_{i,i}‘s, we have values below 1/(n-1), the last one should be below 1 and the sum has to be k=2 . So we have the value of the n \boldsymbol{H}_{i,i}‘s.


Insurance data science : use and value of unusual data #1

Next week, with , I will be at the Summer School of the Swiss Association of Actuaries, in Lausanne, with Jean-Philippe Boucher (UQAM) and Ewen Gallic (AMSE).

I will give an introductionary talk on Monday morning, and the slides are now available

There will be some hands-on applications, on R. I will share some codes in the slides.

Optimal transport on large networks

With Alfred Galichon and Lucas Vernet, we recently uploaded a paper entitled optimal transport on large networks on arxiv.

This article presents a set of tools for the modeling of a spatial allocation problem in a large geographic market and gives examples of applications. In our settings, the market is described by a network that maps the cost of travel between each pair of adjacent locations. Two types of agents are located at the nodes of this network. The buyers choose the most competitive sellers depending on their prices and the cost to reach them. Their utility is assumed additive in both these quantities. Each seller, taking as given other sellers prices, sets her own price to have a demand equal to the one we observed. We give a linear programming formulation for the equilibrium conditions. After formally introducing our model we apply it on two examples: prices offered by petrol stations and quality of services provided by maternity wards (only the later is described here for privacy issues). These examples illustrate the applicability of our model to aggregate demand, rank prices and estimate cost structure over the network. We insist on the possibility of applications to large scale data sets using modern linear programming solvers such as Gurobi.

Demand for gas in gas stations in Britanny, and demand for maternity in France (with border correction)

In addition to this paper we released a R toolbox to implement our results and an online tutorial, optimalnetwork.github.io.

On my way to Manizales (Colombia)

Next week, I will be in Manizales, Colombia, for the Third International Congress on Actuarial Science and Quantitative Finance. I will be giving a lecture on Wednesday with Jed Fress and Emilianos Valdez.

I will give my course on Algorithms for Predictive Modeling on Thursday morning (after Jed and Emil’s lectures). Unfortunately, my computer locked itself last week, and I could not unlock it (could not IT team at the university, who have the internal EFI password). So I will not be able to work further on the slides, so it will be based on the version as-at now (clearly in progress).


Pareto Models for Top Incomes

With Emmanuel Flachaire, we uploaded on hal a paper on Pareto Models for Top Incomes,

Top incomes are often related to Pareto distribution. To date, economists have mostly used Pareto Type I distribution to model the upper tail of income and wealth distribution. It is a parametric distribution, with an attractive property, that can be easily linked to economic theory. In this paper, we first show that modelling top incomes with Pareto Type I distribution can lead to severe over-estimation of inequality, even with millions of observations. Then, we show that the Generalized Pareto distribution and, even more, the Extended Pareto distribution, are much less sensitive to the choice of the threshold. Thus, they provide more reliable results. We discuss different types of bias that could be encountered in empirical studies and, we provide some guidance for practice. To illustrate, two applications are investigated, on the distribution of income in South Africa in 2012 and on the distribution of wealth in the United States in 2013.

This paper was presented at and UCSB and in several workshops this spring, and this Summer, Emmanuel will present it at ECINEQ.

Note that a R package is also available on github, TopIncomes.

Estimates on training vs. validation samples

Before moving to cross-validation, it was natural to say “I will burn 50% (say) of my data to train a model, and then use the remaining to fit the model”. For instance, we can use training data for variable selection (e.g. using some stepwise procedure in a logistic regression), and then, once variable have been selected, fit the model on the remaining set of observations. A natural question is usually “does it really matter ?”.

In order to visualize this problem, consider my (simple) dataset


Let us generate 100 training samples (where we keep about 50% of the observations). On each of them, we use a stepwise procedure, and we keep the estimates of the remaining variables (and their standard deviation actually)

for(i in 1:100){
idx = which(sample(0:1,size=n, replace=TRUE)==1)

Then, for the 7 covariates (and the constant) we can look at the value of the coefficient in the model fitted on the training sample, and the value on the model fitted on the validation sample (of course, only when they were remaining)

for(j in 1:8){

For instance, with the intercept, we have the following


where horizontal segments are confidence intervals of the parameter on the model fitted on the training sample, the vertical on the validation sample. The green part means some sort of consistency, while the red one means that actually, the coefficient was negative with one model, positive with the other one. Which is odd (but in that case, observe that coefficients are rarely significant).

We can also visualize the joint distribution of the two estimators,

for(j in 1:8){
idx = which(!is.na(M1[,j]))
Z = cbind(M1[idx,j],M2[idx,j])
H = Hpi(x=Z)
fhat = kde(x=Z, H=H)

which are here, almost on the diagonal,

meaning that the intercept on the two samples is (more or less) the same. We can then look at other parameters (which is actually more interesting).

On that variable, it seems that it is significant on the training dataset (somehow, it is consistent with the fact that it is remaining in the model after the stepwise procedure) but not on the validation sample (or hardly significant).

Others are much more consistent (with some possible outliers)



On the next one, we have again significance on the training sample, but not on the validation sample,



and probably more interesting

where the two are very consistent.

Bailey (1963) and Poisson regression on two factors

Consider the following dataset, from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy, with 563 observations, and nine variables : eight covariates, and the variable of interest, the number of extramarital affairs, over a year,

base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE)
'data.frame':	563 obs. of  9 variables:
 $ SEX         : int  1 0 0 1 1 0 0 1 0 1 ...
 $ AGE         : num  37 27 32 57 22 32 22 57 32 22 ...
 $ YEARMARRIAGE: num  10 4 15 15 0.75 1.5 0.75 15 15 1.5 ...
 $ CHILDREN    : int  0 0 1 1 0 0 0 1 1 0 ...
 $ RELIGIOUS   : int  3 4 1 5 2 2 2 2 4 4 ...
 $ EDUCATION   : int  18 14 12 18 17 17 12 14 16 14 ...
 $ OCCUPATION  : int  7 6 1 6 6 5 1 4 1 4 ...
 $ SATISFACTION: int  4 4 4 5 3 5 3 4 2 5 ...
 $ Y           : int  0 0 0 0 0 0 0 0 0 0 ...

Let us focus on two categorical covariates, related to the importance of religion, and the occupation

              expo = 1)
religion  1  2  3  4  5  6  7
       1  4  1  8  4 16  9  0
       2 23  3 11 17 56 36  6
       3 29  1 10 12 39 25  2
       4 38  7 12 21 59 44  2
       5 13  1  3 10 19 19  3
religion  1  2  3  4  5  6  7
       1  4  1 13  3 13  7  0
       2  1  1 13 10 25 43 10
       3 15  0 12 11 34 35  1
       4 24  1  3 15 11  9 10
       5  6  0  0  6 11  7  0

The two tables above are the exposure (number of observations) and the number of extramarital affairs, here as contingency tables. Without any covariate, one can assume that N\sim\mathcal{P}(\lambda\cdot E), where \lambda would be

[1] 0.6305506

The idea with the margin method is to assume that N_{i,j}=E_{i,j}\cdot\lambda_{i,j} where \lambda_{i,j}=A_i\cdot B_j. Bailey (1963) added two series of constraints : per row, \sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j for any i and similarly, for any j \sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_jFrom the first series of constraints, write A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j} and use the second series to write B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}Because we need A_i‘s to compute B_j‘s, and conversely, it is natural to consider some iterative procedure to solve it. Observe that we do not have unicity…

Consider here some starting values for A_i‘s and B_j‘s

[1] 1 1 1 1 1
[1] 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506

The predicted number of extramarital affairs would be \hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j

E * A%*%t(B)
religion          1          2          3          4          5          6          7
       1  2.5222025  0.6305506  5.0444050  2.5222025 10.0888099  5.6749556  0.0000000
       2 14.5026643  1.8916519  6.9360568 10.7193606 35.3108348 22.6998224  3.7833037
       3 18.2859680  0.6305506  6.3055062  7.5666075 24.5914742 15.7637655  1.2611012
       4 23.9609236  4.4138544  7.5666075 13.2415631 37.2024867 27.7442274  1.2611012
       5  8.1971581  0.6305506  1.8916519  6.3055062 11.9804618 11.9804618  1.8916519
[1] 26.48313
[1] 95.84369
        1         2         3         4         5 
 26.48313  95.84369  74.40497 115.39076  42.87744 
[1] 107
[1] 13
  1   2   3   4   5   6   7 
107  13  44  64 189 133  13

From expressions above, observe that one can very easily write expressions of A_i‘s and B_j‘s as functions of B_j‘s and A_i‘s respectively


Let it iterate one thousand times

for(i in 1:1000){

We obtain here

        1         2         3         4         5 
1.5404346 1.0447195 1.4825650 0.6553159 0.6634763 
        1         2         3         4         5         6         7 
0.4685515 0.2629769 0.8454435 0.7245310 0.4889697 0.7770553 1.6753750 
E * A%*%t(B)
religion          1          2          3          4          5          6          7
       1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
       2 11.2586111  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
       3 20.1450811  0.3898804 12.5342484 12.8899708 28.2722423 28.8008726  4.9677044
       4 11.6678702  1.2063307  6.6483904  9.9707299 18.9053460 22.4055332  2.1957997
       5  4.0413463  0.1744790  1.6827951  4.8070914  6.1639760  9.7955975  3.3347148

That is our prediction, per category, of the number of affairs. Observe that here, sums per row are equal to observed numbers,

  1   2   3   4   5 
 41 103 108  73  30 
apply(E * A%*%t(B),1,sum)
  1   2   3   4   5 
 41 103 108  73  30

as well as sums per colums

  1   2   3   4   5   6   7 
 50   3  41  45  94 101  21 
apply(E * A%*%t(B),2,sum)
  1   2   3   4   5   6   7 
 50   3  41  45  94 101  21

Now, why should I mention that here, in the section on the Poisson regression in our course ? Because actually, this is exactly what we get if we run a Poisson regression on those two covariates

            Estimate Std. Error z value Pr(&gt;|z|)    
(Intercept) -0.32604    0.21325  -1.529 0.126285    
religion2   -0.38832    0.18791  -2.066 0.038783 *  
religion3   -0.03829    0.18585  -0.206 0.836771    
religion4   -0.85470    0.19757  -4.326 1.52e-05 ***
religion5   -0.84233    0.24416  -3.450 0.000561 ***
occupation2 -0.57758    0.59549  -0.970 0.332083    
occupation3  0.59022    0.21349   2.765 0.005699 ** 
occupation4  0.43588    0.20603   2.116 0.034381 *  
occupation5  0.04265    0.17590   0.242 0.808399    
occupation6  0.50587    0.17360   2.914 0.003569 ** 
occupation7  1.27415    0.26298   4.845 1.27e-06 ***
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

First of all, observe that the total sum of predictions equals the total sum of observations

yp = predict(reg,type="response")
[1] 355
[1] 355

But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey’s technique

df$religion          1          2          3          4          5          6          7
          1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
          2 11.2586112  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
          3 20.1450813  0.3898804 12.5342484 12.8899708 28.2722424 28.8008726  4.9677044
          4 11.6678703  1.2063307  6.6483904  9.9707300 18.9053460 22.4055332  2.1957997
          5  4.0413464  0.1744790  1.6827951  4.8070914  6.1639761  9.7955975  3.3347148
E * A%*%t(B)
religion          1          2          3          4          5          6          7
       1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
       2 11.2586111  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
       3 20.1450811  0.3898804 12.5342484 12.8899708 28.2722423 28.8008726  4.9677044
       4 11.6678702  1.2063307  6.6483904  9.9707299 18.9053460 22.4055332  2.1957997
       5  4.0413463  0.1744790  1.6827951  4.8070914  6.1639760  9.7955975  3.3347148

To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for A_i‘s

          religion2 religion3 religion4 religion5 
1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 
        1         2         3         4         5 
1.0000000 0.6781979 0.9624329 0.4254098 0.4307072

but also for B_j‘s

            occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 
  1.0000000   0.5612551   1.8043769   1.5463210   1.0435773   1.6584203   3.5756477 
        1         2         3         4         5         6         7 
1.0000000 0.5612551 1.8043770 1.5463210 1.0435773 1.6584203 3.5756478

This will have major implications in non-life insurance models (for claims reserving).

October, grant proposal season

In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad.

Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes

download.file(url,destfile = "GSC.html")

We’re going to need a small function, to remove the $ and other symbols that pollute the data (and prevent them from being treated as numbers)

Correction = function(x) as.numeric(gsub('[$,]', '', x))

We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before)

grants= function(gsc){
    download.file(url,destfile = "GSC.html")
    X=as.character(tables[[1]]$"Awarded Amount")

The average amounts of individual grants can be compared,


In mathematics, the average grant amount is $24400. If we normalize by this quantity, we obtain


In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths…

We can also look at the median values (rather than the averages)


Here again, it is in mathematics that it is the weakest….


in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry.

Cumulative functions can also be ploted,


with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics…

Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know…

« Dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales »

Le statisticien et économiste Alfred Sauvy est resté dans les mémoires pour avoir inventé en 1952 le terme “tiers-monde”. Mais on lui a aussi attribué la paternité de la phrase suivante « dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales ». J’ai du l’entendre alors que j’étais étudiant, et depuis, elle me suit partout.

J’y repensais l’autre jour, quand Mathieu Gallard mentionnait sur Twitter le graphique suivant (correspondant a la popularité du président de la république, en France, dans les 18 mois qui suivent l’élection).

Le tweet disait (entre autres) “ est à ce stade de son quinquennat légèrement moins populaire que “. Pour rappel, ces courbes “de popularité” sont construites par un sondage, avec a chaque fois environ 1000 personnes interrogées (“sur la taille de l’échantillon on est toujours entre 950 et 970 interviews” me disait Mathieu). Bon, tous ceux qui ont des souvenirs de cours de stats se souviennent qu’avec 1000 personnes interrogées, 2 points de différence, c’est rarement significatif. Mais plus globalement, compte tenu de la marge d’erreur, je me suis demande pourquoi les courbes n’étaient pas lissées ? Ça éviterait les discussions stériles pour une variation de 2 points par exemple..

Si on reprend les données brutes (merci Mathieu), on a ici

lines(rate[,5],type="b",col="dark green")
text(18.15,rate[17,5],"JC",col="dark green")

ce qui donne le même que dans le tweet (même si je n’ose pas interpoler linéairement les valeurs manquantes – il y en a deux dans mon fichier)

Prenons la courbe la plus récente, celle d’Emmanuel Macron (avec en plus les valeurs manquantes pour pimenter un peu) et rajoutons les intervalles de confiance ponctuels.


Personnellement, j’aurais bien voulu (1) lisser tout ça, (2) rajouter quelque chose qui s’apparente a des bandes de confiance. Mais avec des erreurs de mesure (c’est comme ça qu’on peut interpréter le fait que les points viennent d’un sondage), je ne sais pas trop quoi faire. J’ai tenté la méthode suivante : le tire au hasard des points dans l’intervalle de confiance, puis je lisse sur ce nouveau jeu de points. Et je répète mille fois

for(s in 1:1000){
  if(s&lt;100) lines(seq(0,18,by=.25),yp,col="light blue")
lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.95)),col="red",lty=2)
lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.05)),col="red",lty=2)

On voit que notre courbe lissée est réaliste, voire même les pseudo-bandes de confiance autour. Pour obtenir ces trois courbes, on peut utiliser la fonction suivante

for(s in 1:1000){
upr=apply(Y,1,function(x) quantile(x,.975)),
lwr=apply(Y,1,function(x) quantile(x,.025)))

Sur les quatre colonnes de notre tableau, ça donne

lines(Y$x,Y$pred,col="dark green",lwd=2)
text(18.65,Y$pred[73],"JC",col="dark green")

Pourquoi les instituts de sondages, qui produisent les courbes de popularité, ne montrent pas ce genre de courbes ? Elles sont – a mon avis – aussi justes que celles qu’ils fournissent, au centième près, jouant sur une précision que l’incertitude ne devrait pas autoriser…

Analyse des résultats au baccalauréat des séries générales

Pour continuer sur les manipulation de données publiques, je voulais m’inspirer du projet de Cédric, de la formation Data Science pour l’Actuariat sur les résultats au baccalauréat. Les données nécessaires à cette étude sont disponibles sur plusieurs sites,

Il ne s’agit aucunement d’une analyse poussée des résultats, juste un peu de visualisation, sans aucune autre prétention ! Ah oui, même si on ne va pas faire de carte (je les trouve peu lisibles) on va quand même utiliser les données spatiales : les établissements scolaires sont géolocalises, et on peut obtenir des informations locales, sur le taux de chômage, ou le revenu médian. Et faire des graphiques.

Ce préambule passé, on peut commencer.


On va commencer par récupérer par établissement, les résultats au bac.

url_resultat_etab = "https://data.education.gouv.fr/explore/dataset/fr-en-indicateurs-de-resultat-des-lycees-denseignement-general-et-technologique/download/?format=csv&amp;timezone=Europe/Berlin&amp;use_labels_for_header=true"
download.file(url_resultat_etab,destfile = paste0(librairie,"import_resultat_etab.csv"), method="curl")
df_resultat_etab = read.csv("import_resultat_etab.csv",header=TRUE, sep= ";", encoding="UTF-8")

Comme bien souvent avec les données des administrations françaises, on a souvent des soucis de typographie. Pour simplifier, on va supprimer les accents, et uniformiser un peu les noms

MiseEnForme_Colonnes = function(text) {
  text &lt;- gsub("è", "e", text)  
  text &lt;- gsub("é", "e", text)         
  text &lt;- gsub("_", ".", text)
  text &lt;- gsub("serie.", "", text)
  text &lt;- gsub("Effectif.Presents.", "Effectif.", text)
  text &lt;- gsub("Taux.","Tx.",text)
  text &lt;- gsub("Brut.de.Reussite.", "Admis.Etab.", text)
  text &lt;- gsub("Reussite.Attendu.", "Admis.", text)
  text &lt;- gsub("brut", "Etab", text)
  text &lt;- gsub("attendu", "Academie", text)
  text &lt;- gsub("toutes.", "TOTAL", text)
  text &lt;- gsub("Total.", "TOTAL", text)
  text &lt;- gsub("..Etablissement", ".Etab", text)
  text &lt;- gsub("Pourcentage", "Tx", text)
for(i in 1:ncol(df_resultat_etab)){
  colnames(df_resultat_etab)[i] &lt;- MiseEnForme_Colonnes(names(df_resultat_etab)[i])

On va ensuite supprimer les départements et régions d’outre-mer,

df_resultat_etab = df_resultat_etab[-which(toupper(df_resultat_etab$Departement) %in% c("GUADELOUPE","MAYOTTE","MARTINIQUE","REUNION","GUYANE")),]

récupérer les noms des colonnes

Colonnes = colnames(df_resultat_etab)

et comme on s’intéresse aux premières variables

Colonnes_Generiques = Colonnes[1:8]

on les recupere, pour construire quelques statistiques pour colonnes relatives aux séries L, ES et S

Colonnes_Series = Colonnes[grepl("([a-zA-Z]*?.)*\\.S$|([a-zA-Z]*?.)*\\.ES$|([a-zA-Z]*?.)*\\.L$|([a-zA-Z]*?.)*\\.TOTAL$",Colonnes)]

Et on finit avec les autres

Colonnes_Autres = Colonnes[grepl("(Tx.Bacheliers.*)|(Tx.acces.*)|(Effectif.de.*)|(libelle.region)|(code.region)|(element)",Colonnes)] 
df_resultat_etab = cbind(df_resultat_etab[Colonnes_Generiques],df_resultat_etab[Colonnes_Series],df_resultat_etab[Colonnes_Autres])

On peut aussi localiser les établissements

url_carto_etab &lt;- "https://www.data.gouv.fr/s/resources/adresse-et-geolocalisation-des-etablissements-denseignement-du-premier-et-second-degres/20160526-143453/DEPP-etab-1D2D.csv"
df_carto_etab = read.csv2("import_carto_etab.csv",header=TRUE

On récupère ici la géolocalisation de 66556 établissements nationaux ! On peut croiser avec des données socio-économiques des communes

nom_base_emploi = "base-cc-emploi-pop-act-2014"
url_baseemploi_popactive = paste0("https://www.insee.fr/fr/statistiques/fichier/2862207/",nom_base_emploi,".zip")
unzip(paste0(nom_base_emploi,".zip"),overwrite = TRUE) 
df_base_emploi_source = read_excel(paste0(nom_base_emploi,".xls"),sheet="COM_2014",skip=5)

On va exclure les territoires d’outre-mer ici

df_base_emploi_source &lt;- df_base_emploi_source[-which(df_base_emploi_source$DEP %in% c("971","972","973","974","975")),]
df_base_emploi_colonnes = c("CODGEO","P14_POP1564","P14_H1564","P14_F1564","P14_ACT1564","P14_ACTOCC1564","P14_CHOM1564","P14_INACT1564", "P14_ETUD1564", "P14_RETR1564", "P14_AINACT1564", "P14_HCHOM1524", "P14_FCHOM1524", "C14_ACT1564","C14_ACT1564_CS1","C14_ACT1564_CS2","C14_ACT1564_CS3","C14_ACT1564_CS4","C14_ACT1564_CS5","C14_ACT1564_CS6","P14_POP15P")
df_base_emploi = df_base_emploi_source[,names(df_base_emploi_source) %in% df_base_emploi_colonnes]

et corriger les soucis classiques de la Corse,

MiseEnForme_CodeGeo = function(text) {
  text &lt;- gsub("2A", "20", text)  
  text &lt;- gsub("2B", "20", text)  
df_base_emploi$CODGEO = MiseEnForme_CodeGeo(df_base_emploi$CODGEO)

On peut aussi utiliser des données de revenus, par communes

nom_base_revenus = "indic-struct-distrib-revenu-2014-COMMUNES"
url_baserevenus = paste0("https://www.insee.fr/fr/statistiques/fichier/3126151/",nom_base_revenus,".zip")
unzip(paste0(nom_base_revenus,".zip"),overwrite = TRUE)
df_base_revenus = read_excel("FILO_DISP_COM.xls",sheet="ENSEMBLE",skip=5)[,c(1,4,7)]
df_base_revenus$CODGEO = MiseEnForme_CodeGeo(df_base_revenus$CODGEO)

On recupere des donnees spatiales relatives aux communes

url_geoloc_communes = "http://www.nosdonnees.fr/wiki/images/b/b5/EUCircos_Regions_departements_circonscriptions_communes_gps.csv.gz"
df_geoloc_communes = read.csv2(gzfile("geoloc_communes.csv.gz"),header=TRUE, stringsAsFactors = FALSE,encoding="UTF-8")

et comme toujours, un peu de corrections s’imposent

df_geoloc_communes = df_geoloc_communes[-which(df_geoloc_communes$numéro_département %in% c("971","972","973","974","975")),]
df_geoloc_communes = df_geoloc_communes[,names(df_geoloc_communes) %in% c("code_insee","latitude","longitude","codes_postaux")]
df_geoloc_communes_nb &lt;- nrow(df_geoloc_communes)

On va ensuite creer une fonction de remplacement des valeurs manquantes, et de correction des séparateurs décimaux

MiseEnForme_CoordonneesGeo = function(valeur){
pretraitement = ifelse(as.character(valeur)=="-","0",as.character(valeur))
traitement = as.numeric(ifelse(pretraitement==".","0",gsub(pattern=",",replacement=".",pretraitement)))
df_geoloc_communes$latitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$latitude)
df_geoloc_communes$longitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$longitude)

On passe ensuite a l’élimination des lignes en double

df_geoloc_communes = unique(df_geoloc_communes)

On va ensuite changer les noms des colonnes pour harmoniser avec les autres bases

names(df_geoloc_communes) = c("Codes_Postaux","CODGEO","coordonnee_y","coordonnee_x")

On peut ensuite rechercher les lignes en double sur les codes insee

liste_CODGEO2 = aggregate(x=df_geoloc_communes$Codes_Postaux,by=list(df_geoloc_communes$CODGEO),FUN="length")
list_geoloc_communes_CODGEO2 = liste_CODGEO2[liste_CODGEO2$x&gt;1,1]
df_geoloc_communes_CODGEO2 = df_geoloc_communes[df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,1:2]

Ici, un correction manuelle s’impose pour 4 configurations : les données propres à Lyon, Paris et Marseille ne sont pas géolocalisées, les données propres à la ville de Laguépie sont géolocalisées en doubles

df_geoloc_communes_propre = df_geoloc_communes[!df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,]
df_geoloc_communes_corrige = data.frame(Codes_Postaux=c("13001","69001","75001","82250"),                        CODGEO=c("13055","69123","75056","82088"),
df_geoloc_communes = rbind(df_geoloc_communes_propre,df_geoloc_communes_corrige)

On peut enfin fusionner les bases

df_etab = merge(df_resultat_etab,df_carto_etab,by="Cod.Etab")

Certains établissements ne peuvent être géolocalisés pour certaines années

df_etab_total_nongeolocalises &lt;- df_resultat_etab[!df_resultat_etab$Cod.Etab %in% df_carto_etab$Cod.Etab,]

Comme l’étude ne porte que sur les seuls lycées d’enseignement général et technologique, le dataframe est réduit aux observations relatives d’une part aux lycées, d’autre part aux établissements d’enseignement polyvalent, général ou général et technologique.

df_etab = df_etab[grep("LYCÉE",toupper(df_etab$nature_uai_libe)),]
df_etab = df_etab[grep("GÉNÉRAL|POLYVALENT",toupper(df_etab$nature_uai_libe)),]
df_etab_nongeolocalises = df_etab[df_etab$Cod.Etab %in% df_etab_total_nongeolocalises$Cod.Etab,]
df_etab_geolocalise = df_etab[!is.na(df_etab$coordonnee_x),]
df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$coordonnee_y),]

Enfin, on va convertir les code géographiques (ici reconnus comme facteur) en chaines de caractères (de 5 caractères) pour pouvoir fusionner les tables

ConvertCODGEO = function(code) {
  if(is.character(code)) {
    code_character = ifelse(nchar(code)&lt;5, paste0("0",code), code)
  else if(is.factor(code)){
    code_character = ifelse(code&lt;10000, paste0("0",as.numeric(as.character(code))), as.numeric(as.character(code)))
  else if(is.numeric(code)){
    code_character = ifelse(code&lt;10000, paste0("0",code), as.character(code))
df_etab_geolocalise$Code.commune = ConvertCODGEO(df_etab_geolocalise$Code.commune) 
df_etab_geolocalise$Secteur.Public.Prive = sapply(df_etab_geolocalise$Secteur.Public.Prive,function(nature) {ifelse(nature=="PU","Lycées Publics","Lycées Privés")})

On conservation alors les établissements dont la commune n’est pas manquante

df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$Code.commune),]

Pour finir, on va creer une base, pour ensuite faire une graphique

tbl_etab_nature_res_source = df_etab_geolocalise[,c(3,8,9,10,11,13,14,15)]
for(i in c(6,7,8)){
  temp = tbl_etab_nature_res_source[!is.na(tbl_etab_nature_res_source[i]),c(1,2,i-3,i)]
  temp$Serie = ifelse(i==6,"L",ifelse(i==7,"ES","S"))
  names(temp)[2:4] = c("Nature","Effectif","Tx.Admis")
    tbl_etab_nature_result = temp
    tbl_etab_nature_result = rbind(tbl_etab_nature_result,temp)
graph = ggplot(tbl_etab_nature_result,aes(x=Effectif,y=Tx.Admis,colour=factor(Annee))) 
graph = graph + geom_point(alpha=0.45)
graph = graph + facet_grid(Serie~Nature)
graph = graph + xlab("Effectifs de l'établissement en terminale (par série)") + ylab("Taux d'admission (%)") 
graph = graph + scale_color_discrete(name="Année des\nrésultats")
graph = graph + theme(legend.title = element_text(size=9,face="bold"),
       legend.text = element_text(size=9),
       strip.background = element_rect(colour="black", fill="gray95"),
       panel.border = element_rect(linetype = "solid"),
       panel.grid.major = element_line(colour = "gray75",linetype = "dashed"),
       panel.grid.minor = element_line(colour = "gray95",linetype = "dashed"),
       axis.title.x = element_text(size=9, face="bold"),
       axis.text.x  = element_text(size=8),
       axis.title.y = element_text(size=9, face="bold"),
       axis.text.y  = element_text(size=8))

On a ici l’evolution des resultats en fonction de la taille des etablissements.

df_communes_CorrNaN = df_communes_Corr[which(!df_communes_Corr$TxChomage == "NaN" &amp; !df_communes_Corr$TxCadres == "NaN" &amp; !df_communes_Corr$TxOuvriers == "NaN" &amp; !df_communes_Corr$NbPopulation == "NaN" &amp; !df_communes_Corr$TxSenior == "NaN" &amp; !df_communes_Corr$RevenusMedians == "NaN"),]
df_communes_sp = SpatialPointsDataFrame(coords = df_communes_CorrNaN[, c("coordonnee_x", "coordonnee_y")], data = df_communes_CorrNaN) 
Grille              = as.data.frame(makegrid(df_communes_sp, nsig=2, cellsize = 0.1))
names(Grille)       = c("X", "Y")
coordinates(Grille) = c("X", "Y")
gridded(Grille)     = TRUE  
fullgrid(Grille)    = TRUE  
proj4string(Grille) = proj4string(df_communes_sp)

On peut ensuite faire du krigeage, histoire de lisser un peu nos donnees de chomage et de revenu

df_communes_sp.TxChomage = krige(TxChomage ~ 1, df_communes_sp, Grille, nmax=1)
df_communes_sp.RevenusMedians = krige(RevenusMedians ~ 1, df_communes_sp, Grille, nmax=1)
sp_lycee_WGS84@data$TxChomage      = extract(R.TxChomage,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxCadres       = extract(R.TxCadres,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxOuvriers     = extract(R.TxOuvriers,sp_lycee_WGS84)
sp_lycee_WGS84@data$NbPopulation   = extract(R.NbPopulation,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxSenior       = extract(R.TxSenior,sp_lycee_WGS84)
sp_lycee_WGS84@data$RevenusMedians = extract(R.RevenusMedians,sp_lycee_WGS84)

On peut enfin conclure, en faisant une fonction generique de visualisation

Creation_Graphique = function(df, AnneeObs_Ouv, AnneeObs_Clo, Effectifs, Abscisses, Ordonnees, TitreAbs, TitreOrd, CouleurGraph, CouleurLiss, Serie) {
  df_temp = df[which(df$Annee&gt;=AnneeObs_Ouv &amp; df$Annee&lt;=AnneeObs_Clo),]
  df_temp = df_temp[which(!is.na(df_temp[,Effectifs]) &amp; !is.na(df_temp[,Abscisses]) &amp; !is.na(df_temp[,Ordonnees])),]
  df_temp = df_temp[!df_temp[,Effectifs]==0,]
  df_temp = df_temp[,c(Abscisses,Ordonnees)]
  graphique = ggplot(df_temp,aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees])) 
  graphique = graphique + geom_point(data = df_temp, aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees]),size=1, color=CouleurGraph,alpha=0.25) 
  graphique = graphique + geom_density2d(aes(colour=..level..),show.legend=F) + scale_colour_gradient(low="gray55",high="gray25") 
  graphique = graphique + scale_y_continuous(breaks= seq(80,100,by=2), limits = c(80,100))
  graphique = graphique + xlab(TitreAbs) + ylab(TitreOrd) 
  graphique = graphique + ggtitle(Serie) 
  graphique = graphique + theme(plot.title   = element_text(size=13,color=CouleurLiss, face="bold", hjust=0),
       axis.title.x = element_text(size=8, face="bold"),
       axis.text.x  = element_text(size=8),
       axis.title.y = element_text(size=8, face="bold"),
       axis.text.y  = element_text(size=8),
       panel.border = element_rect(linetype = "solid"),
       panel.grid.major = element_line(colour = "gray55",linetype = "dashed"),
       panel.grid.minor = element_line(colour = "gray75",linetype = "dashed")) 
graphique = graphique + stat_smooth(method = "loess",fill=CouleurLiss,color=CouleurLiss)
Production_Graphique_VI_1 = function(df, Titre_General, Axe_Abscisses, Titre_Abscisses, Annee_Observee_Ouv, Annee_Observee_Clo){
  Graph_S = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.S", Axe_Abscisses, "Tx.Admis.Etab.S", Titre_Abscisses, "Taux d'admission (%)", "dodgerblue3","dodgerblue4","Série S")
  Graph_ES = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.ES", Axe_Abscisses, "Tx.Admis.Etab.ES", Titre_Abscisses, "Taux d'admission (%)","darkorange2","darkorange3","Série ES")
  Graph_L = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.L", Axe_Abscisses, "Tx.Admis.Etab.L", Titre_Abscisses, "Taux d'admission (%)", "chartreuse4","darkgreen","Série L")
  Graph_TS = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.Etab", Axe_Abscisses, "Tx.Admis.Etab", Titre_Abscisses, "Taux d'admission (%)", "indianred1","red4","Toutes séries")
  p = plot_grid(Graph_S, Graph_ES, Graph_L, Graph_TS, ncol = 2, nrow = 2,align = 'hv',
  scale = c(0.95, 0.95, 0.95, 0.95),vjust = 0.9, hjust=-0.5)
  titre &lt;- ggdraw() + draw_label(Titre_General,fontface="bold", size=10)
  plot_grid(titre, p, ncol = 1, rel_heights=c(.25,5))

On note ici

df_lycee &lt;- sp_lycee_WGS84@data

et on peut faire un premier graphique, avec le taux de chomage

Production_Graphique_VI_1(df              = df_lycee,
                          Titre_General   = "Taux d'admission par série en fonction du taux de chômage \n dans la population active - Tous lycées confondus",
                          Axe_Abscisses   = "TxChomage",
                          Titre_Abscisses = "Taux de chômage dans la population active (%)",
                          Annee_Observee_Ouv     = "2013",
                          Annee_Observee_Clo     = "2015")

On ne va pas enfoncer les portes ouvertes de l’inference ecologique, en affirmant des choses aussi stupides que “on a moins de chances d’avoir le bac quand on est au chômage”. Mais on peut noter que dans les zones avec un fort taux de chômage, les résultats au bac sont moins bons.

On peut ensuite regarder en fonction du revenu de la commune du lycée

Production_Graphique_VI_1(df                     = df_lycee,
                          Titre_General          = "Taux d'admission par série en fonction du niveau des revenus disponibles médians - Tous lycées confondus",
                          Axe_Abscisses          = "RevenusMedians",
                          Titre_Abscisses        = "Quantile du niveau des revenus disponibles médians (%)",
                          Annee_Observee_Ouv     = "2013",
                          Annee_Observee_Clo     = "2015")

Fascinant, non ? Mais c’est clairement juste une première approche… il faudrait aller plus loin ensuite !

Scraper, ou pas ?

Ce matin, je mettais en ligne un billet scraper la base d’incendies de forêts expliquant comment scraper la base Promethee, en remplaçant automatique une formulaire. Dans la soirée, sur Twitter me faisait remarquer que c’était intéressant, mais peut être inutile (sur cet exemple particulier en tous cas).

Sur Firefox, il existe un “Moniteur Réseau“, qui s’ouvre en cliquant sur le clavier CtrlMaj + E ( CommandOption + E quand on est est sur un Mac). Ça ouvre un espace en bas de la fenêtre

Si maintenant je lance a la main ma requête, on voit ce qui est fait

Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite

Dans l’onglet ‘Reponse’, on voit l’information sur le fichier json créé,

on a plus d’information bien entendu, si on agrandit

et surtout si on va dans l’onglet ‘En-tetes’

on a l’URL complet de la requête ! On peut ensuite copier le lien, et l’ouvrir

On a ainsi a accès a presque toute la table. Quand on regarde l’URL, seuls 20 incendies sont renvoyés, a cause du &nbrLigne=20& mais on peut tenter de le changer. propose de mettre &nbrLigne=10000& et il affirmait que ça marchait. J’avoue avoir essayé, sans succès.

Mais je retiens le tuyau en tous cas. Parfois, on peut faire simple !