Chose promise, chose due: mardi prochain, le 27 septembre, le cours d’Actuariat IARD se déroulera en salle informatique (PK-S1525) à l’horaire habituel. Sinon quelques notes de cours sur la tarification a priori sont en ligne ici.

• Quelques références en R

Avant de commencer à programmer en R, quelques références qui peuvent être utiles, en français pour commencer,

• “R pour les débutants” d’Emmanuel Paradis, (PDF)
• “Introduction à la programmation en S” par Vincent Goulet, (PDF)
• “Statistique de l’Assurance” par Arthur Charpentier (PDF) qui sont mes notes de cours de l’an passé, mais qui insistent sur l’utilisation de R, pas sur la programmation en R, qui est supposée un peu connue.

Sinon en anglais, les références sont un peu plus nombreuses,

• “R for Beginners” d’Emmanuel Paradis (PDF),
• “An Introduction to R” par Longhow Lam (PDF)
• “The R language — a short companion” par Marc Vandemeulebroecke (PDF),
• “The R Guide” par Jason Owen (PDF),
• “Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
• “Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
• “Statistics with R and S-Plus” d’Hugo Quené (PDF)
• “Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
• “Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Le code pour importer les données est le suivant,

> sinistre <- read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
+ header=TRUE,sep=";")
> sinistres=sinistre[sinistre$garantie=="1RC",]
> contrat <- read.table("http://freakonometrics.free.fr/contractACT2040.txt",
+ header=TRUE,sep=";")
> T=table(sinistres$nocontrat)
> T1=as.numeric(names(T))
> T2=as.numeric(T)
> nombre1 = data.frame(nocontrat=T1,nbre=T2)
> I = contrat$nocontrat%in%T1
> T1= contrat$nocontrat[I==FALSE]
> nombre2 = data.frame(nocontrat=T1,nbre=0)
> nombre=rbind(nombre1,nombre2)
> basenb = merge(contrat,nombre)
> head(basenb)
> basesin=merge(sinistres,contrat)
> basesin=basesin[basesin$cout>0,]

• Faire des graphiques en R

Dans A practicioner’s guide to Generalized Linear Models (en ligne ici), on peut voir des graphiques comme celui ci-dessous, avec la fréquence de sinistre en fonction de l’age (ou plutôt de classes d’ages). Sur le graphique ci-dessous, on a également une distinction entre les hommes (en bleu) et les femmes (en rouge),

sur ce graphique, les fréquences sont exprimées en logarithme du multiplicateur (i.e. en variation par rapport à la moyenne du portefeuille).
On peut utiliser le code suivant pour générer (automatiquement) des graphiques similaires,

> graphique=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=as.factor(basenb[,nom])}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyenne=variance=n=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyenne[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variance[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyenne[k])^2,E[X==FREQ[k]])
+ n[k]       =sum(E[X==FREQ[k]])
+}
+ w=barplot(n,names.arg=FREQ,col="light green",axes=FALSE,
+ xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ axis(2)
+ par(new=TRUE)
+ IC1=moyenne+1.96/sqrt(n)*sqrt(variance)
+ IC2=moyenne-1.96/sqrt(n)*sqrt(variance)
+ moyenneglobale=sum(Y)/sum(E)
+ 
+ if(type==1){
+ plot(mid,moyenne,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyenne,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobale,lty=2,col="red")}
+
+ if(type==2){
+ plot(mid,log(moyenne/moyenneglobale),ylim=
+ range(c(log(IC1/moyenneglobale),log(IC2/moyenneglobale))),
+ type="b",col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,log(IC1/moyenneglobale),mid,
+ log(IC2/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC1/moyenneglobale),mid+.1,
+ log(IC1/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC2/moyenneglobale),mid+.1,
+ log(IC2/moyenneglobale),col="blue")
+ points(mid,log(moyenne/moyenneglobale),pch=19,col="red")
+ axis(4)
+ abline(h=0,lty=2,col="red")}
+
+ mtext("Exposition", 2, line=2, cex=1.2,col="light green")
+ if(type==1){mtext("Fréquence annualisée", 
+     4, line=-2, cex=1.2,col="red")}
+ if(type==2){mtext("Fréquence annualisée (log multiplicateur)", 
+    4, line=-2, cex=1.2,col="red")}
+ }

Par exemple en utilisant un découpage arbitraire par classe d’age (comme cela est fait par défaut dans la fonction),

> graphique()

Mais on pourrait aussi utiliser un découpage assurant d’avoir des classes plus grandes, par exemple en utilisant les quantiles,

> Q=quantile(basenb[,"ageconducteur"],(0:10)/10)
> Q[1]=Q[1]-1
> graphique(nom="ageconducteur",niveau=Q,continu=TRUE)

Ces deux graphiques permettent de visualiser la fréquence empirique de sinistre, par classe d’age, sans modèle paramétrique. Des intervalles de confiance sont également représentés (basés sur une hypothèse de normalité). Notons que l’on peut aussi faire une représentation relative à la fréquence moyenne (en log des multiplicateurs),

> graphique(type=2)

ou encore, si on cherche à analyser la fréquence en fonction de la zone géographique d’habitation

> graphique(nom="zone",continu=FALSE,type=2)

Une autre utilisation peut être faite sur la sévérité des sinistres (coût moyen) et la fréquence, par exemple un assureur,

Il est possible de modifier un peu la fonction pour ajouter au graphique la sévérité des sinistres, e.g.

> graphiquecout=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=basenb[,nom]}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyennen=variancen=nn=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennen[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variancen[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyennen[k])^2,E[X==FREQ[k]])
+ nn[k]       =sum(E[X==FREQ[k]])
+ }
+ moyenneglobalen=sum(Y)/sum(E)
+ 
+ if(continu==TRUE){X=cut(basesin[,nom],niveau)}
+ if(continu==FALSE){X=basesin[,nom]}
+ Y=basesin$cout
+ FREQ=levels(X)
+ moyennes=variances=ns=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennes[k] =mean(Y[X==FREQ[k]])
+ variances[k]=var(Y[X==FREQ[k]])
+ ns[k]=length(Y[X==FREQ[k]])
+ }
+ moyenneglobales=mean(Y)
+ 
+ w=barplot(nn,names.arg=FREQ,col="light green",
+ axes=FALSE,xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ 
+  par(new=TRUE)
+ IC1=moyennen+1.96/sqrt(nn)*sqrt(variancen)
+ IC2=moyennen-1.96/sqrt(nn)*sqrt(variancen)
+ plot(mid,moyennen,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyennen,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobalen,lty=2,col="red")
+ 
+ par(new=TRUE)
+ IC1=moyennes+1.96/sqrt(ns)*sqrt(variances)
+ IC2=moyennes-1.96/sqrt(ns)*sqrt(variances)
+ plot(mid,moyennes,ylim=range(c(IC1,IC2)),type="b",
+ col="blue",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="blue")
+ segments(mid-.1,IC1,mid+.1,IC1,col="blue")
+ segments(mid-.1,IC2,mid+.1,IC2,col="blue")
+ points(mid,moyennes,pch=19,col="blue")
+ axis(2)
+ abline(h=moyenneglobales,lty=2,col="blue")
+ 
+ mtext("Cout moyen", 2, line=2, cex=1.2,col="blue")
+ mtext("Fréquence annualisée", 4, line=-2, cex=1.2,col="red")
+ }

Si on regarde par classe d’age,

> graphiquecout()

ou encore, en fonction de l’age du véhicule,

> Q=quantile(basenb[,"agevehicule"],(0:10)/10)
> Q[1]=Q[1]-1
> graphiquecout(nom="agevehicule",niveau=Q,continu=TRUE)

ou enfin, en fonction de la zone géographique

> graphiquecout(nom="zone",continu=FALSE)