Category Archives: Course

Variables Catégorielles et Modèle Logistique

Petit complément, suite aux coquilles qu’il y avait dans les slides, sur une propriété de la régression logistique quand on régresse sur des variables catégorielles. On était sur la base des avocats

> avocat <- read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") > avocat$CLMSEX <- factor(avocat$CLMSEX, labels=c("M","F")) > avocat$MARITAL <- factor(avocat$MARITAL, labels=c("M","C","V","D")) > avocat=avocat[!is.na(avocat$CLMSEX),]
> attach(avocat)
> sum((ATTORNEY==2)&(CLMSEX=="F"),na.rm=TRUE)/sum(CLMSEX=="F",na.rm=TRUE)
[1] 0.5256065
> sum((ATTORNEY==2)&(CLMSEX=="M"),na.rm=TRUE)/sum(CLMSEX=="M",na.rm=TRUE)
[1] 0.4453925

Autrement dit, la proportion d’hommes qui se sont fait représentés par un avocat est de 44.539%, et la proportion de femmes 52.56%. On peut visualiser le tableau croisé ci-dessous

> tab=xtabs(~ATTORNEY+CLMSEX,data=avocat)
> require(vcd)
> mosaic(tab, shade=TRUE, legend=TRUE)

Quand on fait une régression linéaire (Gaussienne), on retrouve ces probabilités dans les valeurs de coefficients: pour les femmes, on a 44.539%, et pour les hommes, on a la différence avec les femmes (cette dernière modalité étant la modalité de référence ici)

> reglm = lm((ATTORNEY==2) ~ CLMSEX, data=avocat)
> summary(reglm)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.44539 0.02060 21.620 < 2e-16 ***
CLMSEXF     0.08021 0.02756 2.911  0.00367 **

C’est ce que l’on a ci-dessous

> sum(coefficients(reglm))
[1] 0.5256065
> coefficients(reglm)[1]
(Intercept)
0.4453925

On a un résultat similaire avec une régression logistique

> reglogit = glm((ATTORNEY==2) ~ CLMSEX, data=avocat,family=binomial)
> summary(reglogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.21930 0.08312 -2.639 0.00833 **
CLMSEXF      0.32182 0.11097  2.900 0.00373 **

même si les coefficients n’ont pas la même interprétation. En transformant ces coefficients, on retrouve très exactement les mêmes valeurs

> exp(sum(reglogit$coefficients[1:2])) /(1+exp(sum(reglogit$coefficients[1:2])))
[1] 0.5256065
> exp(reglogit$coefficients[1]) /(1+exp(reglogit$coefficients[1]))
(Intercept)
0.4453925

Désolé pour la typo.

Actuariat de l’Assurance Non-Vie #3

Mardi, seconde partie du cours d’actuariat, avec les modèles de classification le matin, mais l’après midi, on devrait commencer les modèles de comptage. Les slides sont en ligne.

Le matin, je dois intervenir 15 minutes à l’IHP vers 9 heures, dans un colloque sur Artificial Intelligence for Fintech and Insurtech: si le RER traîne un peu trop, entre Luxembourg et Lozère, j’aurais peut-être 5 minutes de retard….

Articles for the Probability and Statistics Project

Here are some articles for the project for the graduate crash course on probability and statistics,

For those willing to work on datasets, consider success per school in a national exam, over time (here is the file)

base=read.table("http://freakonometrics.free.fr/brevet_rennes.csv", header=TRUE,sep=";",dec=",")

Some datasets will be uploaded soon

Métro: centralité et robustesse

Demain matin, nous aurons un TP pour le cours sur les réseaux, les flux et les transports. En particulier, en nous inspirant des travaux de Sybil Derrible, nous allons commencer par étudier la centralité dans les différents systèmes de métro, mais aussi la robustesse. Les matrices d’adjacence d’une trentaine de métros dans le monde sont en ligne dans un fichier xls. Histoire de gagner un peu de temps, le code pour créer une matrice d’adjacence peut être le suivant

loc="/data/Metro_Networks_Adjacency.xls"
library(xlsx)
E=read.xlsx(loc,"StPetersburg")
n=nrow(E)
nom=as.character(E[3:(n-2),1])
Adj=E[3:(n-2),(4:ncol(E)-1)]
Adj[is.na(Adj)]=0
Adj=as.matrix(Adj)
colnames(Adj)=rownames(Adj)=nom

On est ensuite prêt à manipuler le réseau,

library(igraph)
iflo=graph_from_adjacency_matrix(Adj,mode = "undirected")
plot(iflo)

On va utiliser les notions vues en cours, sur la centralité, mais surtout, on travaillera sur Quantifying the robustness of metro networks, inspiré de The complexity and robustness of metro networks. Plusieurs fonctions utiles sont déjà programmées dans R, comme l’assortativité.

 

 

Base de Données, Devoir d’Actuariat

Pour les deux premiers devoirs du cours d’actuariat de l’assurance non-vie, on utilisera les mêmes données, dont le descriptif complet est en ligne ici. La première base d’apprentissage (sa version corrigée) est aussi en ligne

> base=read.csv2("http://freakonometrics.free.fr/PG3/PG_2017_SURV_YEAR0.csv",sep=",")
> dim(base)
[1] 100000 32
> table(base$survenance)

FALSE TRUE
85636 14364

On reviendra en cours mardi prochain sur les attendus, mais l’idée sera de modélisation la variable survenance. Pour se second projet, on passera à la modélisation de la prime pure.

Actuariat de l’Assurance Non-Vie #2

Mardi, on commencera le cours d’actuariat, avec un cours d’introduction le matin, mais l’après midi, on attaquera plus sérieusement, avec une mise en bouche, avec les modèles de classification (de type {survenance} vs {non survenance} d’accident). Là encore, les slides sont en ligne.

Je parlerais aussi du premier devoir, qui portera sur les modèles que nous (re)verrons ici.

 

Actuariat de l’Assurance Non-Vie #1

Reprise du cours d’assurance non-vie à l’ENSAE la semaine prochaine. Les slides sont en ligne dans le fichier ENSAE-2017.7z, pour la globalité des slides, sinon je mettrais les liens au fur et à mesure (en corrigeant les typos, voire en rajoutant des sections, si besoin), avec cette semaine les slides pour l’introduction du cours. Et pour illustrer, j’ai emprunté quelques images au très beau livre Karambolage d’Arnold Odermatt.

Networks with R

In order to practice with network data with R, we have been playing with the Padgett (1994) Florentine’s wedding dataset (discussed in the lecture). The dataset is available from

> library(network)
> data(flo)
> nflo=network(flo,directed=FALSE)
> plot(nflo, displaylabels = TRUE,
+ boxed.labels =
+ FALSE)

The next step was to move from the network package to igraph. Since we have the adjacency matrix, we can use it

> library(igraph)
> iflo=graph_from_adjacency_matrix(flo,
+ mode = "undirected")
> plot(iflo)

The good thing is that a lot of functions are available, for instance we can get shortest paths, between two specific nodes. And we can give appropriate colors to the nodes that we’ll cross

> AP=all_shortest_paths(iflo,
+ from="Peruzzi",
+ to="Ginori")
> L=AP$res[[1]]
> V(iflo)$color="yellow"
> V(iflo)$color[L[2:4]]="light blue"
> V(iflo)$color[L[c(1,5)]]="blue"
> plot(iflo)

We can also visualize edges, but I found it slightly more complicated (to extract edges from the output)

> liens=c(paste(as.character(L)[1:4],
+ "--",
+ as.character(L)[2:5],sep=""),
+ paste(as.character(L)[2:5],
+ "--",
+ as.character(L)[1:4],sep=""))
> df=as.data.frame(ends(iflo,E(iflo)))
> names(df)=c("src","target")
> lstn=sort(unique(c(as.character(df[,1]),as.character(df[,2]),"Pucci")))
> Eliens=paste(as.numeric(factor(df[,1],levels=lstn)),"--",
+ as.numeric(factor(df[,2],levels=lstn)),sep="")
> EU=unlist(lapply(Eliens,function(x) x%in%liens))
> E(iflo)$color=c("grey","black")[1+EU]
> plot(iflo)

But it works. It is also possible to use some D3js visualization

> library( networkD3 )
> simpleNetwork (df)

Then the next question was to add a vertice to the network. The most simple way to do it is probability through the adjacency matrix

> flo2=flo
> flo2["Pucci","Bischeri"]=1
> flo2["Bischeri","Pucci"]=1
> nflo2=network(flo2,directed=FALSE)
> plot(nflo2, displaylabels = TRUE,
+ boxed.labels =
+ FALSE)

Then, we’ve been playing with centrality measures.

> plot(iflo,vertex.size=betweenness(iflo))

The goal was to see how related they were. Here, for all of them, “Medici” is the central node. But what about the others?

> B=betweenness(iflo)
> C=closeness(iflo)
> D=degree(iflo)
> E=eigen_centrality(iflo)$vector
> base=data.frame(betw=B,close=C,deg=D,eig=E)
> cor(base)
betw close deg eig
betw 1.0000000 0.5763487 0.8333763 0.6737162
close 0.5763487 1.0000000 0.7572778 0.7989789
deg 0.8333763 0.7572778 1.0000000 0.9404647
eig 0.6737162 0.7989789 0.9404647 1.0000000

Those measures are quite correlated. It is also possible to use a hierarchical graph to visualize how close those centrality measures can be

> H=hclust(dist(t(base)),
+ method="ward")
> plot(H)

Instead of looking at values of centrality measures, it is possible to looks are ranks

> rbase=base
> for(i in 1:4) rbase[,i]=rank(base[,i])
> H=hclust(dist(t(rbase)),
+ method="ward")
> plot(H)

Here the eigenvector measure is very close to the degree of vertices.

Finally, it is possible to seek clusters (in the context of coalition here, in case a war should start between those families)

> kc <- fastgreedy.community ( iflo )

Here we have 3 classes (+1 for the node that is disconnected from the other families)

> V(iflo)$color=c("yellow","orange",
+ "light blue")[membership ( kc )]
> plot(iflo)

> plot(kc,iflo)