Category Archives: Actuarial science

Quelques ordres de grandeurs sur les accidentés de la route

Hier soir, Benoît m’a posé une colle par courriel, et je dois avouer que j’ai horreur de ça… Il me demandait un ordre de grandeur descoûts moyens en France (du point de vue de l’assureur) d’un blessé et d’un mort pour les sinistres automobiles“. La colle m’a embarrassé parce que déteste ne pas avoir la moindre idée des ordres de grandeurs sur des sujets qui pourtant m’intéressent. En fait, je connais les coûts extrêmes, les plus gros accidents, mais pas les coûts moyens… Moralité je me suis mis au travail, et comme toujours sur ces sujets, Audrey m’a beaucoup aidé. Alors avant de commencer, il faut savoir que quand on parle de gravité d’un accident corporel (dans les précédant billets, ici ou , je ne parlais que des nombres) il faut parler de taux d’AIPP, atteinte permanente à l’intégrité physique et/ou psychique (anciennement appelées taux d’IPP, i.e. taux d’incapacité permanente partielle). Alors c’est très technique, c’est fixé par des docteurs, des caisses d’assurance maladie, des juges, bref, il y a tellement de paramètres que c’est dur de s’y retrouver. Sauf qu’il y a des espèces de grilles… Elles sont évoquées ici par exemple,

  • Perte, paralysie d’un doigt ou perte de l’odorat ou d’une partie de la voix =  1 à 5 %
  • Perte de l’ ouie d’une oreille ou d’un orteil genou remplacé par une prothèse =  5 à 10 %
  • Perte de tous les orteils, perte de toutes les dents, perte de la voix =  10 à 15 %

(bref, un taux de 10% c’est pas très drôle.. mais ce n’est pas ce qui est retenu pour définir un accident grave)

  • Perte du pouce droit (pour un droitier) ou perte de la vision d’un œil =  20 à 25 %
  • Perte de la main droite (pour un droitier) ou perte du bras gauche =  40 à 50 %

Là on va franchir le cap des 50% qui correspondent à une classification particulière, à savoir les graves,

  • Perte du bras droit (pour un droitier), perte de la jambe ou perte totale de la vue ou de l’ouïe =  50 à 60 %
  • Para ou tétraplégie et  traumatisme crâniens =  60 à 100 %

Voilà pour les ordres de grandeurs sur la gravité d’un accidenté.
En 2008, il y a eu de l’ordre de 7 millions de sinistres (avec une indemnisation de l’assureur). Près de la moitié concernent des bris de glaces, et environ 250 000 correspondaient à un accident corporel. Si on creuse un peu, parmi les accidentés, dans 75% des accidents, il s’agit de blessés légers, sans AIPP. 24% sont accidentés avec un AIPP, avec la répartition suivante: 19% ont un taux d’AIPP inférieur à 5%, et 2% dépassent 20%. Enfin, 1% décèdent. Voilà pour la répartition des nombres, regardons maintenant les coûts.
Ces 3% des sinistres (en nombres) représentent 30% du coût total des sinistres, soit 4,5 milliards d’euros (vu à la fin de la première année, selon la FFSA). Et ces milliards d’euros sont répartis de la manière suivante

  • 11% pour les blessés léger
  • 36% pour les AIPP compris entre 1% et 20%
  • 42% pour les AIPP supérieurs à 20%
  • 11% pour les décès

Si on regarde l’indemnité moyenne en 2008 (pour en venir à l’objet de mon billet)

  • les blessés avec AIPP ont été indemnisé à hauteur de 45,680 euros (dont 14,165 euros au titre des dépenses de santé, actuelles et projetées)
  • les ayants droits, en cas de décès, ont été indemnisé à hauteur de 89,800 euros,

Si l’on regarde les blessés légers (sans AIPP), ils ont touché en moyenne 1,500 euros (la médiane se situant plutôt vers 750 euros).
Pour les blessés avec AIPP

  • taux entre 1% et 5%, 9,195 euros
  • taux entre 6% et 9%, 30,828 euros
  • taux entre 10% et 14%, 53,295 euros,
  • taux entre 15% et 19%, 86,950 euros,
  • taux entre 20% et 29%, 130,585 euros,
  • taux entre 20% et 29%, 130,585 euros,
  • taux supérieur à 30%, 644,100 euros,
  • taux supérieur à 50%, 988,000 euros,
  • taux supérieur à 80%, 2,018,000 euros,

Si l’on regarde plus en détail les blessés graves (taux d’AIPP supérieur à 50%) l’indemnité total était proche du million d’euros, répartis de la manière suivante

  • 21% au titre du déficit fonctionnel permanent
  • 6% au titre de préjudices personnels,
  • 25% au titres des dépenses de santé (actuelles et projetées)
  • 5% au titre de pertes de gains professionnels
  • 35% au titre d’assistance par une tierce personne

(le reste étant d’autres frais). De plus, l’indemnité décroît avec l’âge

  • moins de 30 ans, 1,521,000 euros
  • entre 30 et 60 ans, 1,287,000 euros,
  • plus de 60 ans, 580,000 euros,

Comme le note la FFSA, l’évolution de l’indemnité moyenne croît régulièrement depuis 10 ans, avec une croissance annuelle de l’ordre de 9%.

Catastrophes naturelles et population mondiale

Pour reprendre une question posée au téléphone par François, et par Eric dans un commentaire (ici), il est possible de comparer l’évolution du nombre de catastrophes et de l’évolution de la population (mondiale). Pour les données de population on peut en trouver sur wikipedia (ici), sur le site des nations unies () après 1950, ou encore sur le site du census.gov (). Ensuite on peut moyenner les doublons.

POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=data.frame(years=1900:2010))

Revenons à nos catastrophes “naturelles“. Par exemple les tremblements de terre. La première piste (évoquée ) était de lisser cette série pour l’utiliser en baseline (multiplicatif). On retrouve ici une série qui oscille autour de 1.

earthquake=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet.csv",
sep=",",header=TRUE)
flood=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_2_.csv",
sep=",",header=TRUE)
storm=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_3_.csv",
sep=",",header=TRUE)
library(splines)
n=nrow(earthquake)
earthquake=earthquake[-((n-3):n),]
DATE=1900:2010
EARTHQUAKE=FLOOD=STORM=rep(0,length(DATE))
EARTHQUAKE[as.numeric(as.character(earthquake$Year))-1899]=
as.numeric(as.character(earthquake$Number))
plot(DATE,EARTHQUAKE,type="b")
reg=lm(EARTHQUAKE~bs(DATE,6))
lines(DATE,predict(reg),col="red")
plot(DATE,EARTHQUAKE/predict(reg),type="b")
abline(h=1,lty=2,col="blue")

Si on utilise cette fois la population mondiale comme baseline, on obtient

POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=
     data.frame(years=1900:2010))/1000000
plot(DATE,EARTHQUAKE/baseline,type="b")

La série de points présente une jolie rupture. Afin de trouver la date optimale de rupture, utilisons un test basé sur le critère d’Akaike.

Y=EARTHQUAKE/baseline
X=DATE
DUMMY=1910:2000
VAIC=rep(NA,length(DUMMY))
for(x in DUMMY){
IX=(X>x)
reg=lm(Y~1+IX)
VAIC[x-1909]=AIC(reg)
}
plot(DUMMY,VAIC)
k=which.min(VAIC)
DUMMY[k]

Manifestement, il y a une rupture en 1975 (que je ne sais trop comment interpréter).

Si on regarde une autre série, comme celle des inondations.

n=nrow(flood)
flood=flood[-((n-3):n),]
DATE=1900:2010
FLOOD[as.numeric(as.character(flood$Year))-1899]=
as.numeric(as.character(flood$Number))
plot(DATE,FLOOD,type="b")
reg=lm(FLOOD~bs(DATE,6))
lines(DATE,predict(reg),col="red")
population=tapply(POP$POP,POP$YEAR,
mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
baseline=predict(reg,newdata=data.frame(years=1900:2010))
plot(DATE,FLOOD/baseline,type="b")

On retrouve ici clairement le fait que la croissance de la population ne peut expliquer la croissance du nombre d’inondations.

Increasing number of catastrophes and human factor

Peter asked me (here) to discuss the occurrence of natural catastrophes. I have looked at the EM-DAT dataset (here), and more specifically three kinds of risks: earthquakes, floods and storms (in any country over the past century). For earthquakes, we have the following figures,

the blue points are raw data (I mean, public ones, those that can be downloaded on the website), and I added two smoothed patterns, in red (based on splines regressions). This is no reason to assume that there could be an increase of the number of earthquakes. The trend is due to the fact that we look at earthquakes considered as catastrophic, i.e. that caused deaths and severe damages. And as pointed out in Peter’s comment, it means that “the increase in ‘disasters’ might have a lot to do with the development of the population”.
If we look now at storms, we have

i.e. also an increasing trend. But here, we can consider that the trend can be explained by two factors: an inflation due to the value at risk (e.g. increase of concentration of population in some areas as for earthquakes) but maybe also a change that can be related to climate change.  If we remove the pattern du to the increase of value at risk on the raw data (i.e. we divide raw data by the trend observed on earthquakes, the red line), we have

Here it looks like the number of catastrophic storms is increasing (removing the human factor). And it might be possible to relate that trend to climate change. If we look finally at major flood events, we also observe an increase of the pattern,

and we can divide by the trend observed on earthquakes, to obtain

Looking at that graph, it is possible to propose another interpretation about the increase: people might have a shorter memory about rivers than the one they have on earthquakes; and that they think they might influence climate risk, and flood risk, (but not physical ones). I mean that people know where earthquake risk is, and did not really concentrate in those areas during the past century. On the other hand, it is possible to assume that people imagine they can change river design, e.g. building dykes, and though they can control the risk. But obviously they couldn’t.

Evolution of the number of natural catastrophes

Recently (here and there) I mentioned Reinsurers comments on the increase on the number of catastrophes. But actually, there are other sources of information about catastrophes. For instance the http://www.cred.be/, Centre for Research on the Epidemiology of Disasters. As the mention it below, there has been an increase in the number of natural catastrophes. For a disaster to be entered into the so-called EM-database at least one of the following criteria must be fulfilled:

• Ten (10) or more people reported killed.
• Hundred (100) or more people reported affected.
• Declaration of a state of emergency.
• Call for international assistance.

Note that it is possible to download data (here), so I guess I will work further with that database. To be continued…

Encore une année record…

Comme l’a noté la blogosphère ces dernières heures, il semble que des records aient (encore) été battus en terme de catastrophes naturelles: “catastrophes naturelles : une année 2010 exceptionnelle” (ici) ou “2010, année record des catastrophes naturelles” (). Maintenant on peut remonter un peu dans le temps aussi: “catastrophes naturelles: 2008 bat des records” (ici) ou “2008 a battu des records de catastrophes naturelles” (), “en 2007, les catastrophes naturelles ont atteint un nombre record” (ici), “catastrophes naturelles : 2005 bat des records en matière d’intensité et de coûts” (ici) ou “les catastrophes ont entraîné une facture record en 2005” (), “déjà record, le bilan 2004 des catastrophes naturelles s’alourdit” (ici), “année 2000 : record de catastrophes naturelles” (ici) et “catastrophes naturelles : encore une année record” (en 2000, ) etc, etc.

Bon, en tant que statisticien, les rares résultats que je connaisse sur les records sont des théorèmes en log (i.e. en variation lente). Si http://freakonometrics.hypotheses.org/files/2015/12/record01.png est une suite de variables i.i.d. et que l’on définie les durées entre records,

http://freakonometrics.hypotheses.org/files/2015/12/records02.png

http://freakonometrics.hypotheses.org/files/2015/12/record03.png, avec http://freakonometrics.hypotheses.org/files/2015/12/record04.png. Alfred Rényi a montré que

http://freakonometrics.hypotheses.org/files/2015/12/record08.png

alors que Marcel Neuts a montré que

http://freakonometrics.hypotheses.org/files/2015/12/record09.png

En particulier, on note que

http://freakonometrics.hypotheses.org/files/2015/12/record11.png

Autrement dit, plus le temps avance, plus les records doivent devenir rares.
Par exemple, regardons la simulation suivante. Considérons un processus Poissonnien homogène, http://freakonometrics.hypotheses.org/files/2015/12/catnat01.png, avec les dates des catastrophes, et http://freakonometrics.hypotheses.org/files/2015/12/catnat02.png leur coût.

On note

http://freakonometrics.hypotheses.org/files/2015/12/catnat04.png

le nombre annuel de catastrophes (éventuellement, on ne considère que les catastrophes dont le cout dépasse un seuil prédéfini, comme cela semble généralement être le cas).

On peut définir le processus de records,

http://freakonometrics.hypotheses.org/files/2015/12/catnat07.png

qui ressemble à ça

Bon, maintenant, si on a des records tous les ans (ou presque) c’est probablement que l’hypothèse i.i.d. n’était pas valide.
Supposons par exemple qu’il y ait de l’inflation, i.e. les coût sont http://freakonometrics.hypotheses.org/files/2015/12/catnat03.png.

Dans ce cas, le processus de nombre d’événements annuels est

et celui des records

Si on répète l’expérience un grand nombre de fois, le nombre moyen d’événements par année me cesse de croître (en rouge ci-dessous)

et le processus de records moyen (avec le cas i.i.d. en bleu),

Autrement dit, si l’inflation est trop importante, il est possible d’avoir régulièrement de nouveaux records.
Une autre possibilité est de supposer qu’il n’y a pas d’inflation, mais que le nombre de sinistres ne cesse de croître (ici un processus de Poisson à intensité linéaire), c’est à dire une fréquence plus importante des événements climatiques extrêmes

Le nombre annuel de catastrophes est alors

soit, si l’on considère le processus de records

Là encore, si l’on répète l’expérience un grand nombre de fois, le nombre moyen de catastrophes ne cesse de croître

alors que le processus de records a la moyenne suivante

Bref, il est possible d’avoir effectivement des records presque tous les ans. Cela signifie probablement que les sinistres coûtent de plus en plus cher (hausse de la valeur à risque), mais aussi probablement que les catastrophes sont de plus en plus fréquentes. On retrouve finalement ce que je disais il y a quelques mois maintenant (ici).

Is it that stupid to make extremely long term forecast when studying mortality ?

I received recently a comment by FCA (here) who raised an important question, about forecast in dynamic mortality models. (S)he mentioned that from his(her) point of view, the econometric models I considered were “good to predict for the next, say, 3 or 4 years. Not for the next 50 years…”. Which was the message I tried to stress last year in a conference about retirement in France (here). But from a quantitativepoint of view, how inconsistent were forecasts made 35 years ago, or 60 years ago ?

Consider here the Lee Carter model, obtained on the periods 1816-1950 (in black below), 1816-1975 (in red) and 1816-2000 (in blue), unfortunately, it is difficult to compare http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s since we have identifiability problems here. Nevertheless, we if consider affine transformation so that  http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s are equal in 1900 and 1950 (say), we obtain

On that graph, we considered an ETS (AAN) forecast. If we do not consider the entire series for forecasting, but only observations following WWI (1945), we obtain

For sketches of the R code,

T=1980
base0=data.frame(D,E,A,Y,a=as.factor(A),
y=as.factor(Y))
base=base0[base0$Y<=T,]
LC2=gnm(D~a+Mult(a,y),offset=log(E),family=
poisson,data=base)
A=LC2$coefficients[1]+LC2$coefficients[2:110]
B=LC2$coefficients[111:220]
K0=LC2$coefficients[221:length(LC2$coefficients)]
Y=as.numeric(K0)
K1=c(K0,forecast(ets(Y,model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=240)$mean)
MU=matrix(NA,length(A),length(K1))
MU1=MU2=MU
for(i in 1:length(A)){
for(j in 1:length(K1)){
MU1[i,j]=exp(A[i]+B[i]*K1[j])
MU2[i,j]=exp(A[i]+B[i]*K2[j])
}}
x=40
s=seq(0,109-x-1)
t=2000
Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-base1$Year[1]-1])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-base1$Year[1]-1])))
r=.035
m=70
h=seq(0,39)
V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h]
V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h]
M=cbind(V1,V2)
apply(M,2,sum)

Actually, it is not that bad…. even if it is only a qualitative intuition. Again, I am not a demographer, and my interest is more on actuarial science… so if we look at the estimation of annuities (still the same insurance contract, as here) for some insured of age 40 in 2000, we get the following graph (where forecasts http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s were obtained on the complete series, i.e. from 1816 until the year we consider),

(here it means that in 1900, I had to forecast mortality for someone of age 40 in 2000… so we had to forecast mortality with a 150 year horizon). Obviously, even if we are able to forecast improvement of mortality rates, it is not enough since it looks like, each year, improvement are alway higher than what what expected. Note that if we run it twice (since there might be problem with initial values in the econometric procedure) we obtain something similar,

So, the output is consistent. And if we change the way we predict future values, e.g. on focusing only on the past 50 years, i.e.

K1=c(K0,forecast(ets(Y[(length(Y)-50):length(Y)],
model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y[(length(Y)-50):length(Y)],
allowdrift=TRUE),h=240)$mean)

we obtain the following graph for the annuity associated to an insurance contract sold in 2000,

so that relative changes compared with 1980 are (in %)

Hence, over a bit more than 25 years, we underestimated annuities of 25%. We if start to take into account possible investments, it is not so bad, I think….  don’t you think ?

 

Finding roots of functions in actuarial science

The following simple code can be used to find roots of functions (based on the secant algorithm),

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
	x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
	if (abs(fun(x2)) < tolerence)
		return(x2)
	x0 <- x1
	x1 <- x2
}}

It can be interesting in actuarial science, e.g. to find the actuarial rate so that to present values are equal. For instance, consider the following capital, given only if the insured is still alive (this example was initially considered here). We would like to find the rate so that the probable discounted value is 600,

> Lx=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ header=TRUE,sep=";")
> capital=c(100,100,125,125,150,150)
> n=length(capital)
> x=0.035
> X=45
> f=function(x){
+ capital.act=capital*(1/(1+x))^(1:n)
+ PROBA=Lx[((Lx[,1]>X)*(Lx[,1]<=(X+n)))==1,2]/Lx[(Lx[,1]==X)==1,2]
+ return(sum(capital.act*PROBA))}
>
> f1=function(x){f(x)-600}
> secant(f1,0,0.1)
[1] 0.06022313
> f(0.06022313)
[1] 600

*

Millenium bridge, endogeneity and risk management

In less than 48 hours, last week two friends mentioned the Millennium Bridge as an illustration of a risk management concept. There are several documents with that example, here (for the initial idea of using the Millennium Bridge to illustrate issues in risk management) here or there, e.g.

When we mention resonance effects on bridges, we usually thing of the Tacoma Narrows Bridge (where strong winds set the bridge oscillating) or the Basse-Chaîne Bridge (in France, which collapsed on April 16, 1850, when 478 French soldiers marched across it in lockstep). In the first case, there is nothing we can do about it, but for the second one, this is why soldiers are required to break step on bridges.

But for the Millennium bridge, a ‘positive feedback‘ phenomenon (known as Synchronous Lateral Excitation in physics) has been observed: the natural sway motion of people walking caused small sideways oscillations in the bridge, which in turn caused people on the bridge to sway in step, increasing the amplitude of the oscillations and continually reinforcing the effect. That has been described in a nice paper in 2005 (here). In the initial paper by Jon Danielsson and Hyun Song Shin, they note that “what is the probability that a thousand people walking at random will end up walking exactly in step? It is tempting to say “close to zero”, or “negligible”. After all, if each person’s step is an independent event, then the probability of everyone walking in step would be the product of many small numbers – giving us a probability close to zero. Presumably, this is the reason why Arup – the bridge engineers – did not take this into account. However, this is exactly where endogenous risk comes in. What we must take into account is the way that people react to their environment. Pedestrians on the bridge react to how the bridge is moving. When the bridge moves under your feet, it is a natural reaction for people to adjust their stance to regain balance. But here is the catch. When the bridge moves, everyone adjusts his or her stance at the same time. This synchronized movement pushes the bridge that the people are standing on, and makes the bridge move even more. This, in turn, makes the people adjust their stance more drastically, and so on. In other words, the wobble of the bridge feeds on itself. When the bridge wobbles, everyone adjusts their stance, which sets off an even worse wobble, which makes the people adjust even more, and so on. So, the wobble will continue and get stronger even though the initial shock (say, a gust of wind) has long passed. It is an example of a force that is generated and amplified within the system. It is an endogenous response. It is very different from a shock that comes from a storm or an earthquake which are exogenous to the system.

And to go further, they point out that this event is rather similar to what is observed in financial markets (here) by quoting The Economist from October 12th 2000 “So-called value-at-risk models (VaR) blend science and art. They estimate how much a portfolio could lose in a single bad day. If that amount gets too large, the VAR model signals that the bank should sell. The trouble is that lots of banks have similar investments and similar VAR models. In periods when markets everywhere decline, the models can tell everybody to sell the same things at the same time, making market conditions much worse. In effect, they can, and often do, create a vicious feedback loop.

Course on risk measures (in French)

The course on risk measure, in Luminy, starts at 16.00 on Monday (here). The slides can be found here,

Note that additional references can be downloaded on the internet, e.g. the short course on risk measures by Freddy Delbaen (here) or the article from the Encyclopedia of quantitative finance, by Hans Föllmer and Alexander Schied (there). See also here for the paper by Jean Marc Tallon, Johanna Etner and Meglena Jeleva, on decision theory under uncertainty.

Lecture notes on risk and insurance

I just finished some lectures notes on risk and insurance. The notes, that can be downloaded [pdf], are in French, and will be used at the JES (Journées d’Etudes Statistiques), organised at the CIRM (mentioned here). Previous notes on risk measures [pdf] and copulas [pdf]. Again, all comments are welcome…

La tarification avec SAS

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble  que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

DATA contrats;
SET lib.contrats;
lnexpo = log(expo);
run;

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1_Param
           Type3=Genmod1_Var
           Modelfit=Genmod1_InfoModele; 
MODEL nbsin = ageconducteur /
                  dist = poisson   
                  link = log   
                  offset = lnexpo 
                  type3;
RUN; QUIT;

La sortie SAS a alors l’allure suivante

                                  The GENMOD Procedure
                    Critère pour évaluer la qualité de l'ajustement
              Critère                   DF          Valeur       Valeur/DF
              Deviance                63E3      26872.5334          0.4237
              Scaled Deviance         63E3      26872.5334          0.4237
              Pearson Chi-Square      63E3      73275.5362          1.1553
              Scaled Pearson X2       63E3      73275.5362          1.1553
              Log Likelihood                   -18474.2667

       Algorithm converged.
                       Analyse des résultats estimés de paramètres

                              Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

                         Statistiques LR pour Analyse de Type 3
                      Source           DF      Khi 2    Pr > Khi 2
                      ageconducteur     1     148.72        <.0001

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

PROC GAM DATA = base;
MODEL nbsin = spline(ageconducteur) / dist = Poisson;
OUTPUT OUT=gam PREDICTED; 
RUN; QUIT;
PROC SORT DATA = gam NODUPKEY; BY age_cond; RUN; QUIT;

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;
RUN;

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
PROC MEANS DATA = base NOPRINT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo
y_uni_nbsin;
RUN; QUIT;

ensuite, on fait un modèle GLM, et  un modèle GAM, et on récupère les sorties

PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
      gam;
BY age_cond;
RUN;

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1bis_Param
           Type3=Genmod1bis_Var
           Modelfit=Genmod1bis_InfoModele; 
MODEL nbsin = ageconducteur /
                 dist = poisson  
                 link = log      
                 offset = lnexpo 
                 type3           
                 scale = deviance;
RUN; QUIT;

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

                      Analyse des résultats estimés de paramètres

                                   Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

%MACRO CALCUL_AIC_BIC(infomodel=, param=);
    DATA _null_;
    SET &infomodel.;
    IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
    IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
    RUN;
    DATA _null_;
    SET &param.  end=fin;
    RETAIN nb_df 0;
    nb_df = nb_df + df;
    IF fin THEN CALL SYMPUT("k", nb_df);
    RUN;
    DATA genmod_aic_bic;
    SET &param.;
    FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
    Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
    IF Parameter = "Intercept";
    KEEP Loglike K N AIC_CALC BIC_CALC;
    RUN;
    DATA genmod_aic_bic;
    SET genmod_aic_bic;
    Loglike = &loglike.;
    K = &k.;
    N = %eval(&n_etoile. + &k.);
    AIC_CALC = 2 * Loglike + 2 * K;
    BIC_CALC = 2 * Loglike + K * log(N);
    RUN;
    PROC PRINT DATA = genmod_aic_bic;
    RUN; QUIT;
%MEND CALCUL_AIC_BIC;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);

Nombre de morts sur la route en baisse… et ? (partie 1)

Tout le monde semble se réjouir de la baisse du nombre de morts sur la route depuis quelques années (longtemps attribué à la présence de radars automatiques, mais qu’on pourrait aussi associer à l’instauration du permis à point). Et c’est tant mieux. Mais sans vouloir jouer les cyniques, et en tous les cas c’est une statistique qui intéresse probablement davantage les assureurs, il peut sembler intéressant de s’intéresser aussi aux accidents corporels graves….

image volée chez Andy Warhol

Le principal soucis méthodologique est qu’au 1er janvier 2005, la définition d’accidents graves ou encore de décès a changé. Avant 2005, un tué, était une victime décédée sur le coup ou dans les six jours qui suivent l’accident. A partir de 2005, les 6 jours passent à 30 jours. De même, un blessé grave était un blessé dont l’état nécessitait plus de six jours d’hospitalisation. Après 2005, cette classe disparait, et seuls restent les blessés qui ont passé plus de 24 heures à l’hôpital…. Bref, on est un peu bloqué pour vraiment analyser les tendances sur une longue période…

  • Nombre de morts sur la route

Le graphique suivant montre le nombre quotidien de morts sur les routes, passant d’une vingtaine en 2002 à une douzaine (en 2008).

On notera qu’il n’y a pas eu de réelle rupture en 2005, malgré la nouvelle définition de décès.

  • Nombre de blessés graves

En revanche, les statistiques sur les blessés (graves) ont connu une réelle rupture en 2005.

On peut essayer de corriger afin de rendre les séries comparables (au moins en terme d’ordre de grandeur),

Au delà de la tendance (et de la rupture de la forte baisse observée en 2002-2004), on notera que le cycle annuel est plus faible. De la même manière que les cycles dans les naissances semblent avoir diminué (comme je l’évoquais ici), on notera que le cycle dans les accidents graves (mais aussi dans une moindre mesure les décès) a fortement diminué…

L’équipe de France de foot devrait jouer seulement en semaine (ça sauverait des vies)

Il y a quelques semaines, le service de centralisation des statistiques de l’assurance accidents LAA, par l’intermédiaire de Stefan Scholz Odermatt, a publié une étude sur la survenance d’accident automobile les jours de match de foot.

Selon l’étude (mentionnée ici ou ) le nombre d’accidents les jours où il y a match est supérieur aux jours sans match,

Ayant réussi à récupérer une base presque exhaustive d’accidents corporels en France, je peux faire la même étude, entre janvier 2002 et décembre 2007. Les conclusions sont moins flagrantes que sur le cas suisse, en particulier les jours de semaine, où manifestement, personne ne s’intéresse aux matchs (où en tous les cas pas assez pour être moins vigilent au volant).

En revanche, le week end, l’effet devient plus prononcé

avec davantage d’accidents dans la journée où il y a un match de l’équipe de France.

Moralité, l’équipe de France de foot devrait jouer seulement en semaine, à chaque match, cela éviterait une vingtaine d’accidents de la route à chaque fois (et je ne parle que d’accidents corporels, ayant causé des blessures et ayant donné lieu à un constat de police) ! Comme le disait un copain à qui je racontais ça, “l’équipe de France devrait jouer, ça serait un début“…. mais c’est un autre sujet.