Tag Archives: Poisson

Régression sur des variables catégorielles

Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression

> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, 
+ family=poisson(link="log"))
> summary(reg0)

Call:
glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), 
    data = base, offset = log(exposition))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5717  -0.3968  -0.2996  -0.1547  12.6722  

Coefficients:
      Estimate Std. Error z value Pr(>|z|)    
zoneB -2.54187    0.06287  -40.43   <2e-16 ***
zoneA -2.54912    0.05285  -48.23   <2e-16 ***
zoneC -2.38525    0.03753  -63.56   <2e-16 ***
zoneD -2.13454    0.03878  -55.05   <2e-16 ***
zoneE -2.00204    0.03965  -50.49   <2e-16 ***
zoneF -2.06932    0.11547  -17.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 50966  on 50000  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(reg0,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
        1         2         3         4         5 
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,

> base$zone=relevel(base$zone,"B")
> regB=glm(nbre~zone,offset=log(exposition),data=base,
+ family=poisson(link="log"))
> summary(regB)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54187    0.06287 -40.431  < 2e-16 ***
zoneA       -0.00725    0.08213  -0.088 0.929661
zoneC        0.15662    0.07322   2.139 0.032432 *
zoneD        0.40733    0.07387   5.514 3.50e-08 ***
zoneE        0.53983    0.07433   7.263 3.80e-13 ***
zoneF        0.47255    0.13148   3.594 0.000325 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(regB,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
1         2         3         4         5
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,

> base$zone=relevel(base$zone,"A")
> reg=glm(nbre~zone,offset=log(exposition),
> data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54912    0.05285 -48.232  < 2e-16 ***
zoneB        0.00725    0.08213   0.088 0.929661
zoneC        0.16387    0.06482   2.528 0.011471 *
zoneD        0.41458    0.06555   6.324 2.54e-10 ***
zoneE        0.54708    0.06607   8.280  < 2e-16 ***
zoneF        0.47980    0.12699   3.778 0.000158 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.

> base$zonesimple=base$zone
> base$zonesimple[base$zone%in%c("A","B")]="A"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54612    0.04046 -62.937  < 2e-16 ***
zonesimpleC  0.16087    0.05518   2.915  0.00355 **
zonesimpleD  0.41158    0.05604   7.345 2.06e-13 ***
zonesimpleE  0.54408    0.05665   9.605  < 2e-16 ***
zonesimpleF  0.47681    0.12235   3.897 9.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence

> base$zonesimple=relevel(base$zonesimple,"C")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38525    0.03753 -63.557  < 2e-16 ***
zonesimpleA -0.16087    0.05518  -2.915  0.00355 **
zonesimpleD  0.25071    0.05396   4.646 3.39e-06 ***
zonesimpleE  0.38321    0.05460   7.019 2.24e-12 ***
zonesimpleF  0.31593    0.12142   2.602  0.00927 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.06932    0.11547 -17.921  < 2e-16 ***
zonesimpleC -0.31593    0.12142  -2.602  0.00927 **
zonesimpleA -0.47681    0.12235  -3.897 9.74e-05 ***
zonesimpleD -0.06522    0.12181  -0.535  0.59232
zonesimpleE  0.06727    0.12209   0.551  0.58161
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles

> base$zonesimple[base$zone%in%c("D","E","F")]="F"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5660  -0.3959  -0.3004  -0.1547  12.5929

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.07182    0.02696 -76.853  < 2e-16 ***
zonesimpleC -0.31344    0.04621  -6.783 1.18e-11 ***
zonesimpleA -0.47431    0.04861  -9.757  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15698  on 49997  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> library(car)
> linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0"))
Linear hypothesis test

Hypothesis:
zonesimpleD = 0
zonesimpleE = 0

Model 1: restricted model
Model 2: nbre ~ zonesimple

Res.Df Df  Chisq Pr(>Chisq)
1  49997
2  49995  2 5.7073    0.05763 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)

Les doctorants, le sexe et le directeur de thèse

Pourquoi certains chercheurs publient beaucoup, et d’autres moins ? ou de l’hétérogénéité dans une loi de Poisson (oui, le but principal du billet est de parler de la régression de Poisson).

Il y a une vingtaine d’années, Scott Long publiait un article intitulé ‘The Origins of Sex Differences in Science‘. Mais surtout, les données sont téléchargeables ici.

> library(foreign)
> base=read.dta("http://www.ats.ucla.edu/stat/stata/examples/long/couart2.dta")
> head(base)
  art   fem     mar kid5  phd ment
1   0   Men Married    0 2.52    7
2   0 Women  Single    0 2.05    6
3   0 Women  Single    0 3.75    6
4   0   Men Married    1 1.18    3
5   0 Women  Single    0 3.75   26
6   0 Women Married    2 3.59    2

On dispose de plusieurs informations, recueillies auprès de 950 docteurs en biochimie. On a ainsi

  • le nombre de publications au cours des 3 dernières années de leur doctorat,
  • le sexe du docteur,
  • son statut marital (marié(e) ou célibataire), a priori à la fin du doctorat,
  • le nombre d’enfant(s) du docteur agé(s) de moins de 5 ans, là encore à la fin du doctorat
  • le nombre de publications du directeur de thèse (appelé ici mentor) au cours des 3 dernières années

(on a aussi le prestige du département, mais je trouve cette variable trop subjective). Comme on fait un comptage de publications, on pourrait espérer voir une loi de Poisson.

> n=nrow(base)
> lambda=mean(base$art)
> compte=table(base$art)
> N=as.numeric(names(compte))
> empirique=as.numeric(compte[as.character(0:19)])/n
> theorique=dpois(0:19,lambda)
> cbind(0:19,empirique,theorique)
        empirique theorique
 [1,]    0 0.3005  0.1839
 [2,]    1 0.2688    0.3114
 [3,]    2 0.1945    0.2636
 [4,]    3 0.0918    0.1487
 [5,]    4 0.0732  0.0629
 [6,]    5 0.0295  0.0213
 [7,]    6 0.0185  0.0060
 [8,]    7 0.0131  0.0014
 [9,]    8 0.0010  0.0003
[10,]    9 0.0021  0.0000
[11,]   10 0.0010  0.0000
[12,]   11 0.0010  0.0000
[13,]   12 0.0021  0.0000
[14,]   13     NA  0.0000
[15,]   14     NA  0.0000
[16,]   15     NA  0.0000
[17,]   16 0.0010  0.0000
[18,]   17     NA  0.0000
[19,]   18     NA  0.0000
[20,]   19 0.0010  0.0000
> plot(0:19,theorique,col="red",type="p")
> lines(0:19,empirique,type="p",col="blue ")

Mais non… graphiquement, on est relativement éloignés de la loi de Poisson

avec en rouge la fréquence que l’on aurait eu avec une loi de Poisson (de paramètre la moyenne obtenue sur l’ensemble de la population), et en bleu, la fréquence empirique.
Un test rapide permettra de convaincre les sceptiques,

> library(vcd)
> summary(goodfit(base$art,type="poisson"))

	 Goodness-of-fit test for poisson distribution

                      X^2 df     P(> X^2)
Likelihood Ratio 296.3715 13 1.381259e-55

Manifestement, il y a de l’hétérogénéité dans la population des doctorants. Il y a ceux qui publient beaucoup, ceux qui publient peu et ceux qui ne publient pas.

  • Le sexe et le nombre de publications

Peut-être qu’en discriminant par le sexe, on aurait, par sous-populations, des lois de Poisson. Si on étudie la distribution du nombre de publications pour les hommes

versus la distribution du nombre de publications pour les femmes

ça ne semble pas suffire pour expliquer toute l’hétérogénéité. Même s’il semble y avoir une différence ne serait-ce que sur le nombre moyen de publications par sexe,

> tapply(base$art,base$fem,mean)
     Men    Women 
1.882591 1.470309

Là encore, on peut faire un test rapide de comparaison de moyennes

> t.test(base$art~base$fem)

	Welch Two Sample t-test

data:  base$art by base$fem 
t = 3.3298, df = 885.945, p-value = 0.0009049
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.1692780 0.6552866 
sample estimates:
  mean in group Men mean in group Women 
           1.882591            1.470309

Visiblement les deux sous-populations sont différentes, mais le critère du sexe ne suffit pas à expliquer l’hétérogénéité.
Maintenant, notons que le critère retenu est le nombre de publications. La stratégie des personnes de sexe féminin peut être de viser moins de publications, mais dans de meilleures revues. En
tous les cas, le sexe pourrait expliquer en partie (et en partie seulement) l’hétérogénéité observée.

  • Le statut marital et le nombre de publications

On peut aussi comparer le statut marital, en particulier marié(e)

versus célibataire

Là encore, on peut comparer les moyennes, et faire un test d’égalité. Mais cette fois, l’influence est beaucoup plus faible que le sexe.

> t.test(base$art~base$mar)

	Welch Two Sample t-test

data:  base$art by base$mar 
t = -1.1869, df = 710.307, p-value = 0.2356
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.4033993  0.0994165 
sample estimates:
 mean in group Single mean in group Married 
             1.592233              1.744224

la différence entre les moyennes n’étant pas, ici, significative.

  • Les enfants et le nombre de publications

Une variable probablement plus intéressante que le statut marital pourrait être le nombre d’enfants (moi qui finissait la rédaction de ma thèse alors que ma deuxième apprenait à marcher On a la distribution suivante pour le nombre de publications pour un(e) docteur(e) ayant au moins un enfant,

versus un(e) docteur(e) sans enfant,

Mais la encore, le test de comparaison de moyenne n’est pas vraiment concluant,

> base$kid=base$kid5>0
> t.test(base$art~base$kid)

	Welch Two Sample t-test

data:  base$art by base$kid 
t = 0.6136, df = 646.359, p-value = 0.5397
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1803154  0.3442384 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.721202            1.639241

Certes, les docteur(e)s sans enfant publient plus, mais ce n’est pas non plus flagrant. On pourrait se dire que ce sont surtout les femmes qui ont eu des enfants qui sont pénalisées (la encore, car l’objectif est d’expliquer une quantité de publications), mais tout d’abord autant les hommes que les femmes répercutent le fait d’avoir eu, ou non, un enfant; mais surtout la variable de sexe est beaucoup plus discriminante que le fait d’avoir eu des enfants,

> baseH=base[base$fem=="Men",]
> baseF=base[base$fem=="Women",]
> t.test(baseH$art~baseH$kid)

	Welch Two Sample t-test

data:  baseH$art by baseH$kid 
t = 1.3831, df = 491.989, p-value = 0.1673
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1136025  0.6538074 
sample estimates:
mean in group FALSE  mean in group TRUE 
           2.011628            1.741525 

> t.test(baseF$art~baseF$kid)

	Welch Two Sample t-test

data:  baseF$art by baseF$kid 
t = 0.9525, df = 137.906, p-value = 0.3425
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.1764049  0.5043375 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.501466            1.337500
  • Impact du directeur de thèse sur le nombre de publications

Une fois analysées les variables personnelles, si on regardait l’impact du directeur de thèse sur le travail de l’étudiant ? Qui sait, les variables pourraient être corrélées (je ne parle pas ici de relation causale). Par exemple, on peut regarder la distribution du nombre de publications pour un(e) docteur(e) ayant eu un directeur de thèse qui publiait peu (moins que la moyenne ici)

versus un directeur de thèse qui publiait beaucoup (plus que la moyenne),

Cette fois, on a une nette différence

> base$publis=base$ment>=9
> t.test(base$art~base$publis)

	Welch Two Sample t-test

data:  base$art by base$publis 
t = -6.7568, df = 497.017, p-value = 3.96e-11
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -1.2506429 -0.6871632 
sample estimates:
mean in group FALSE  mean in group TRUE 
           1.341338            2.310241

Un directeur de thèse qui publie semble avoir des étudiant(e)s qui publient aussi (ou réciproquement).

  • Régression de Poisson et nombre de publications

La population était initialement très hétérogène. Trop hétérogène pour une loi de Poisson, avec une variance valant plus du double de la moyenne empirique,

> var(base$art)/mean(base$art)
[1] 2.191358

> summary(glm(art~1,data=base,family=quasipoisson))

Call:
glm(formula = art ~ 1, family = quasipoisson, data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8401  -1.8401  -0.5770   0.2294   7.5677  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.52644    0.03761   14.00   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 2.191389)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1817.4  on 914  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Si on régresse sur les différentes variables on retrouve les interprétations que nous avions eues,

> summary(glm(art~mar+fem+I(kid5>0)+ment,data=base,family=quasipoisson))

Call:
glm(formula = art ~ mar + fem + I(kid5 > 0) + ment, family = quasipoisson, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.5080  -1.5615  -0.3626   0.5614   5.4494  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      0.338732   0.081482   4.157 3.53e-05 ***
marMarried       0.149527   0.085110   1.757  0.07928 .  
femWomen        -0.218643   0.074151  -2.949  0.00327 ** 
I(kid5 > 0)TRUE -0.250286   0.085834  -2.916  0.00363 ** 
ment             0.026000   0.002656   9.788  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 1.837376)

    Null deviance: 1817.4  on 914  degrees of freedom
Residual deviance: 1641.1  on 910  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

On peut d’ailleurs faire une régression GAM avec une estimation par splines de l’effet du nombre de publications du directeur de thèse. Ce qui permet de légitimer le fait que la variable intervient de manière linéaire

On pourrait noter la concavité de la courbe qui semble sous-entendre qu’un directeur qui publie énormément peut toujours s’occuper de ses étudiants – et les faire publier plus qu’un directeur qui publie moins, mais en première approximation, la courbe peut être supposée linéaire.
On peut aussi regarder rapidement qui ne publie pas pendant sa thèse,

> summary(glm((art==0)~mar+fem+I(kid5>0)+ment,data=base,family=binomial))

Call:
glm(formula = (art == 0) ~ mar + fem + I(kid5 > 0) + ment, family = binomial, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1609  -0.9043  -0.7019   1.2697   2.5221  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -0.28843    0.18591  -1.551   0.1208    
marMarried      -0.33565    0.18800  -1.785   0.0742 .  
femWomen         0.23697    0.15829   1.497   0.1344    
I(kid5 > 0)TRUE  0.42957    0.19032   2.257   0.0240 *  
ment            -0.08141    0.01261  -6.455 1.08e-10 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1118.7  on 914  degrees of freedom
Residual deviance: 1052.0  on 910  degrees of freedom
AIC: 1062.0

Number of Fisher Scoring iterations: 5

Certain(e)s ont eu un enfant, et semblent avoir été pénalisé(e)s, mais si un doctorant n’a pas publié, c’est presque uniquement à cause de son directeur (qui publie peu).
Pour reprendre l’interprétation sur l’analyse du nombre de publications, par exemple, le fait d’être marié, ou pas, a peu d’influence sur le nombre de publications. Par contre les femmes – toutes choses étant égales par ailleurs – publient 20% de moins que les hommes,

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[3])
 femWomen 
0.8036087

Par contre le fait d’avoir – ou pas – des enfants semble ici beaucoup plus fort, car les docteur(e)s qui ont eu au moins un enfant publient 25% de moins que ceux qui n’en ont pas eu.

> exp(glm(art~mar+fem+I(kid5>0)+ment,data=base,
+ family=quasipoisson)$coefficients[4])
I(kid5 > 0)TRUE 
      0.7785778

Aussi, une femme ayant eu un enfant aura publié 37% de moins qu’un homme qui n’a pas eu d’enfants. Toutes choses étant égales par ailleurs.

  • Toutes choses étant égales par ailleurs ?

C’est ici qu’est l’arnaque, car rien ne garantie que toutes choses puissent être égales par ailleurs. Par exemple, une femme avec enfant (ou pas) peut elle avoir un directeur de thèse qui publie beaucoup ? N’y a-t-il pas un peu d’endogénéité cachée ?
Si on regarde le nombre de publications du directeur de thèse versus le sexe de l’étudiant,

> t.test(base$ment~base$fem)

	Welch Two Sample t-test

data:  base$ment by base$fem 
t = 2.7101, df = 901.2, p-value = 0.006854
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.4587075 2.8673525 
sample estimates:
  mean in group Men mean in group Women 
           9.532389            7.869359

on note que les directeurs de thèse des hommes publient plus que les directeurs de thèse des femmes. Autrement dit, si les femmes publient moins, c’est peut être parce que les directeurs de thèses qui pourraient les aider à publier ne les prennent pas en thèse, non ? Car il n’y a pas de lien entre le fait d’avoir ou pas des enfants à la fin de sa thèse, et le nombre de publications de son directeur,

> t.test(base$ment~base$kid)

	Welch Two Sample t-test

data:  base$ment by base$kid 
t = -1.5626, df = 499.06, p-value = 0.1188
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.5485856  0.2905184 
sample estimates:
mean in group FALSE  mean in group TRUE 
           8.377295            9.506329

Bref, c’est manifestement le directeur de these qui influence le plus le nombre de publication d’un doctorant. Mais il serait intéressant d’avoir les memes statistiques 5 ans apres la soutenance…

STT2700, estimation, tests et coupes du monde

Mercredi, pour le dernier cours, nous allons revenir sur l’estimation, les tests, et plus généralement sur la modélisation statistique. Pour cela, j’avais pensé travailler sur les nombres de buts marqués, par match, lors de différentes coupes du monde de soccer (1982, 1998 et 2010). Je ne mets pas l’intégralité du code aujourd’hui, l’idée est pour l’instant de mettre en ligne des données qui serviront à répondre aux questions qui seront posées mercredi. Le code (accompagné – éventuellement – d’explications théoriques) sera posté par la suite.

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
soccer1998=read.table("http://freakonometrics.free.fr/soccer1998")
S98=(soccer1998$V1+soccer1998$V2)
soccer2010=read.table("http://freakonometrics.free.fr/soccer2010")
S10=(soccer2010$V1+soccer2010$V3)

Les boxplot associés à ces trois échantillons sont les suivants,

On va se poser des questions autour de ces données, par exemple voir s’il est vraisemblance (ou pas) que le nombre moyen de but dans un match (avant prolongation, s’il y en a eu). On peut commencé par essayer de se demander quel modèle utiliser. Classiquement, la loi de Poisson est la plus utilisée (en plus, c’est la seule loi qui est autorisée lorsqu’on publie un billet le 1er avril). Les histogrammes sont les suivants

boxplot(S82,S98,S10,col=c("red","yellow","blue"),
label=("1982","1998","2010"))
hist(S82,breaks=0:11,col="red")
hist(S98,breaks=0:11,col="yellow")
hist(S10,breaks=0:11,col="blue")

Si on compare les fonctions de répartition empiriques à celles de lois de Poisson ajustées par maximum de vraisemblance, on obtient, pour la coupe du monde de 1982

et pour celle de 2010,

Visuellement, l’ajustement semble relativement bon, surtout en 2010. On peut aussi faire un test du chi-deux,

> library(vcd)
> (GF=goodfit(S10,type="poisson"))

Observed and fitted values for poisson distribution
with parameters estimated by ML

 count observed     fitted
     0        7  6.6409703
     1       17 15.0459484
     2       13 17.0442384
     3       14 12.8719508
     4        7  7.2907534
     5        5  3.3036226
     6        0  1.2474617
     7        1  0.4037543

> summary(GF)

	 Goodness-of-fit test for poisson distribution

                      X^2 df  P(> X^2)
Likelihood Ratio 5.586765  5 0.3485255

On voit que l’on accepte l’ajustement par une loi de Poisson. Pour ceux qui veulent une visualisation, sur la figure ci-dessous, on a la densité d’une loi du chi-deux. Le premier trait vertical est la valeur observée, et l’aire jaune est alors la p-value (qui excède largement 5%). En rouge on a 5%, donc le second trait vertical est la borne de la région critique associé au test pour une erreur de première espèce valant 5%,

On peut ensuite faire plein de tests.  On suppose que . On va pouvoir tester

http://freakonometrics.free.fr/test-soccer-04.gif

contre une hypothèse alternative

http://freakonometrics.free.fr/test-<br /><br /> soccer-06.gif

Comme on a une hypothèse sur la loi des observations qui semble robuste, on peut utiliser un test de type rapport de vraisemblance.
On peut aussi faire un test de la forme

http://freakonometrics.free.fr/test-soccer-09.gif

contre

http://freakonometrics.free.fr/test-soccer-10.gif

(histoire de tester des hypothèses simples – qui ont une interprétation). Sinon, comme ce qui nous intéresse, c’est de savoir si on a plus de trois buts par matchs, on peut définir la variable binomiale http://freakonometrics.free.fr/test-soccer-03.gif, en notant que

http://freakonometrics.free.fr/test-soccer-02.gif

est une proportion – donc facile à tester – qui nous intéresse ici compte tenu du problème que l’on cherchera à résoudre. On pourra alors tester, par exemple

http://freakonometrics.free.fr/test-soccer-08.gif

contre

http://freakonometrics.free.fr/test-soccer-07.gif

Ces derniers tests sont alors facile à mettre en œuvre,

> Z=(S10>=3)*1
> prop.test(sum(Z),length(Z),p=1/2,alternative="less")

	1-sample proportions test with continuity correction

data:  sum(Z) out of length(Z), null probability 1/2 
X-squared = 1.2656, df = 1, p-value = 0.1303
alternative hypothesis: true p is less than 0.5 
95 percent confidence interval:
 0.0000000 0.5322764 
sample estimates:
       p 
0.421875

On peut aussi faire des tests de moyenne sur

http://freakonometrics.free.fr/test-soccer-11.gif

Un test de l’hypothèse

http://freakonometrics.free.fr/test-soccer-04.gif

contre une hypothèse alternative

http://freakonometrics.free.fr/test-soccer-06.gif

s’écrit alors
> t.test(S10,mu=3,alternative ="less")

	One Sample t-test

data:  S10 
t = -3.7763, df = 63, p-value = 0.0001775
alternative hypothesis: true mean is less than 3 
95 percent confidence interval:
     -Inf 2.590273 
sample estimates:
mean of x 
 2.265625
Mais on l’aura l’occasion de revoir tous les points du cours, y compris aller peut être un peu plus loin, par exemple sur la comparaison de moyenne entre échantillons,
> t.test(S82,S98,var.equal=FALSE)

	Welch Two Sample t-test

data:  S82 and S98 
t = 0.427, df = 85.266, p-value = 0.6704
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.5503669  0.8514658 
sample estimates:
mean of x mean of y 
 2.807692  2.657143

Fin des débats sur la valeur de π

La rumeur avait agité pas mal de monde en août dernier lors du congrès international de mathématiques, à Hyderabad (même si à l’époque, ce sont surtout les médailles Fields qui avaient retenu toute l’attention, en France en tous les cas), mais finalement le communiqué de l’Union Internationale de Mathématiques (IMU) est finalement tombé hier soir: à compter du 1er juillet, π sera officiellement égal à 4.
Pour ceux qui ont suivi les débats (je n’ai eu que des bruits de couloir), Microsoft avait augmenté la pression ces derniers mois, quand le cap des 5000 milliards de décimales avait été franchi (ici). Comme l’avait dit Bill Gates avant le congrès à Hyderabad, bientôt la moitié de ma mémoire d’un processeur sera dédiée à stocker les décimales de π.

Et comme il l’avait souligné “la recherche sur les décimales de π allant plus vite que la recherche des améliorations de Windows, nous allons rapidement faire face à un choc informatique sans précédant”. Il avait comparé la situation au bug de l’an 2000 (en ajoutant que – comme il y a 11 ans – la transition se fera, selon lui, en douceur).
Dans le communiqué de l’IMU, il est mentionné que “π gardera son interprétation originelle « περίμετρος »” (périmètre en grec), et en particulier, l’IMU utilise la définition géométrique suivante:

http://freakonometrics.blog.free.fr/public/perso2/animation-pi.gif

Mais en quoi cela peut avoir un intérêt pour mon blog (hormis faire mon frimeur pour montrer que j’ai compris une démonstration géométrique). Tout simplement car π joue un rôle central en statistique et en probabilité (même si on a souvent tendance à l’oublier) au travers de la loi normale ! Pareil en risk management, y compris pour va valorisation des options, mais aussi le capital réglementaire calculé dans les accords de Bâle. Les banquiers se sont réjouit de l’annonce de l’IMU hier soir, car pour ceux qui l’auraient oublié, la loi normale intervient partout en finance . En particulier dans les calculs de quantiles. Rappelons que la densité s’écrit

https://freakonometrics.hypotheses.org/files/2015/12/dens-gauss-pi.gif (pour les utilisateurs de R, la version 2.12.3. sera lancée exceptionnellement plus tôt afin d’intégrer cette mise à jour – et je crois que MSExcel a prévu un addins qui sera bientôt en ligne, comme pour le passage à l’an 2000).
Si on regarde “à la main” ce que vont devenir les probabilités de dépassement de seuils, on obtient

> 1-pnorm(2)
[1] 0.02275013
> integrate(f=function(x){exp(-x^2/2)/sqrt(2*pi)},2,+Inf)
0.02275013 with absolute error < 1.5e-05
> integrate(f=function(x){exp(-x^2/2)/sqrt(2*4)},2,+Inf)
0.02016178 with absolute error < 1.3e-05

autrement dit, la probabilité de dépasser 2 va passer de 2.27% à 2.02%. Pour les quantiles à 99.5% (utilisées par les institutions financières), on a

> qnorm4(.995)
[1] 2.543701
> qnorm(.995)
[1] 2.575829

Autrement dit la baisse est finalement relativement faible. On peut alors s’attendre à une (très) légère baisse des fonds propres des banques et institutions financières (même si pour l’instant, les instances de Bâle ne se sont pas encore prononcées).

 

Test d’ajustement et lâcher de bombes

Avant les applications demain en cours, un petit billet (presque d’actualité) sur une application évoquée vendredi dernier sur le test du chi-deux: l’ajustement de lois.

Le problème est le suivant: un pays se fait bombarder. Et les dirigeants doivent se demander si certaines cibles sont visées (auquel cas il peut être légitime de les déplacer) ou au contraire si les tirs sont aléatoires. Quand je dis que des cibles sont visées, j’entends aussi par là que les pilotes savent viser. Car les pilotes ont (je pense, ou j’espère) un carnet de route à suivre, avec des cibles à viser

Ce que l’on va se demander c’est plutôt si les pilotes savent – ou arrive à – viser. Bref, le problème peut se modéliser en supposant que le lancer de bombes se fait (globalement) selon un processus de Poisson. Localement, si les tirs sont aléatoires, on devrait observer des tirages de lois de Poisson. C’est tout du moins la théorie que R. D. Clarke (alors actuaire chez Prudential Assurance Company) a utilisé pendant la seconde guerre mondiale, sur les bombardements à Londres (en ligne ici). Cette histoire (vraie) est reprise par Thomas Pynchon dans Gravity’s Rainbow (en ligne ici)

(l’exemple est même repris dans Feller). Bref, le problème n’est pas de savoir quel serait le paramètre de la loi de Poisson, mais si la loi de Poisson est adaptée, ou pas. C’est ce qu’on appelle un problème d’ajustement de loi.

  • Test du chi-deux et ajustement de lois

Jusqu’à présent, on avait supposé que les observations suivaient une certaine loi, e.g. une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/chi2-16.gif, et on cherchait à tester une hypothèse de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-13.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-14.gif

Ici on va cherche à utiliser un test sur une loi multinomiale, de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-01.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-02.gif

L’hypothèse nulle est ici une égalité vectorielle,

http://freakonometrics.hypotheses.org/files/2015/12/chi2-03.gif

ou encore

http://freakonometrics.hypotheses.org/files/2015/12/chi2-56.gif

Dans le cas d’un test d’ajustement de lois de Poisson, si on suppose que les observations suivent une loi http://freakonometrics.hypotheses.org/files/2015/12/chi2-09.gif, on va utiliser un test sur une loi multinomiale, avec

http://freakonometrics.hypotheses.org/files/2015/12/chi2-08.gif

Le soucis est que, puisque l’on se limite à un vecteur de taille finie, le vecteur ne sera pas dans le simplexe (cf ici). Donc classiquement, pour la dernière valeur, on retient une hypothèse de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-10.gif

Le test est alors basé sur la statistique de Pearson

http://freakonometrics.hypotheses.org/files/2015/12/chi2-11.gif

qui va suivre (si effectivement les observations suivent une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/chi2-09.gif) une loi http://freakonometrics.hypotheses.org/files/2015/12/chi2-12.gif.
Rappelons que cette statistique peut aussi s’écrire

http://freakonometrics.hypotheses.org/files/2015/12/chi2-21.gif

 

  • Application pour étudier la précision d’un lancer de bombes

Pendant la second guerre mondiale, R.D. Clarke étudia les impacts de bombes V1 et V2 tombées dans une zone de 144 km2 dans le sud de Londres (l’article original, publié après guerre en 1946, est en ligne ici). Il divisa cette zone en 576 zones de 0,25 km2 et compta le nombre d’impact dans chacune des zones. Il obtint plus précisément la distribution suivante

nbre impacts par zone 0 1 2 3 4 5 et plus
fréquence (nbre zones) 229 211 93 35 7 1

(en fait, on sait que le “5
et plus” correspond à 7, car sait qu’il y a eu 537 bombes sur 576 zones). Avant de se lancer tête baissée, réfléchissons un peu au type de loi que l’on pourrait utiliser. Pour cela, notons http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif le nombre de points tombés dans un ensemble http://freakonometrics.hypotheses.org/files/2015/12/asubb.gif (ou http://freakonometrics.hypotheses.org/files/2015/12/calA.gif désigne la région globale). Si on suppose qu’un nombre aléatoire http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif de points sont lancés aléatoirement dans http://freakonometrics.hypotheses.org/files/2015/12/calA.gif, et que http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif suit une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/lambda.gif, alors http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif suit une loi de Poisson de paramètre

http://freakonometrics.hypotheses.org/files/2015/12/ratio-aire.gif

Bref, si on regarde plusieurs plusieurs régions http://freakonometrics.hypotheses.org/files/2015/12/calB.gif (de même taille, éventuellement – afin de garder toutes les observations – formant une partition de http://freakonometrics.hypotheses.org/files/2015/12/calA.gif) et si on observe une loi de Poisson, c’est que dans la région http://freakonometrics.hypotheses.org/files/2015/12/calA.gif, les tirs sont faits au hasard.

> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> y=0:7
> (lambda=sum(n*y)/sum(n))
[1] 0.9322917
> prob=dpois(y,lambda)
> freq.theo=sum(n)*prob
> freq.emp =n
> cbind(y,trunc(freq.theo),freq.emp)
     y     freq.emp
[1,] 0 226      229
[2,] 1 211      211
[3,] 2  98       93
[4,] 3  30       35
[5,] 4   7        7
[6,] 5   1        0
[7,] 6   0        0
[8,] 7   0        1

On peut commencer par regarder la log-vraisemblance

> logvrais=function(L){sum(log(dpois(y,L))*n)}
> param=seq(.5,1.5,by=.025)
> LV=sapply(param,logvrais)
> plot(param,LV,type="b",col="blue")
http://freakonometrics.hypotheses.org/files/2015/12/logvrais1-bombes.png

La statistique du chi-deux, elle, ressemble à ça

> chi2=function(L)sum(n)*sum(((n/sum(n)-dpois(y,L))^2)/dpois(y,L))
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")xxx.
http://freakonometrics.hypotheses.org/files/2015/12/chi2-bombes.png

Bon, le soucis c’est qu’en tronquant le vecteur (on suppose que le nombre maximum d’impact est 7), la somme des probabilités ne fait pas exactement un. On peut regrouper dans une même classe les fréquences élevées (comme le fait Clarke dans le papier initial d’ailleurs).

> (n=c(229,211,93,35,7,1))
[1] 229 211  93  35   7   1
> y=0:4
> prob=c(dpois(y,lambda),1-ppois(4,lambda))
> freq.theo=sum(n)*prob
> freq.emp =n
> (Q=sum(1
[1] 1.169155
> 1-pchisq(Q,length(y)-1)
[1] 0.8831505

On retrouve les quantités évoquées dans l’article de Clarke. La  statistique de test vaut 1.16 et la p-value associée est de l’ordre de 88%. On va donc accepter l’hypothèse de loi de Poisson,

> chi2=function(L){
+ sum(n)*sum(2)^2)/
+ c(dpois(y,L),1-ppois(4,L)) )}
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")

La valeur de la statistique du chi-deux en fonction du paramètre de la loi de Poisson est représentée ci-dessous,

http://freakonometrics.hypotheses.org/files/2015/12/chi2-bombes-v2.png

et le trait horizontal est la valeur seuil de la région critique (en dessous, on accepte l’ajustement d’une loi de Poisson),

> abline(h=qchisq(.95,length(y)-1),lty=2)

Il existe plusieurs fonctions sous R permettant de faire des choses semblables,

> library(vcd)
> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> nsim=c(rep(y[0],n[0]),rep(y[1],n[1]),
       rep(y[2],n[2]),rep(y[3],n[3]),
       rep(y[4],n[4]),rep(y[5],n[5]),
       rep(y[6],n[6]),rep(y[7],n[7]),
       rep(y[8],n[8]))
> gf=goodfit(nsim,type="poisson",method="ML")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Likelihood Ratio 4.007784 3 0.2606249
 
> gf=goodfit(nsim,type="poisson",method="MinChisq")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Pearson 1.275499 3 0.7349592

Bref, quelle que soit la méthode utilisée, on notera que l’on accepte toujours l’hypothèse d’une loi de Poisson. Autrement dit les bombes étaient envoyés un peu au hasard dans Londres….

Et si on y réfléchit un peu, d’un point de vue de la théorie des jeux, c’est effectivement une stratégie optimale…

  1. freq.theo-freq.emp)^2)/freq.theo []
  2. n/sum(n)-c(dpois(y,L),1-ppois(4,L []

Talk at Desjardins General Insurance

This afternoon, I will give a talk at the seminar of the R&D department at Desjardins General Insurance, on correlation in claims reserving. A lot of interesting papers have been published recently on that topic. On multivariate Chain Ladder, some interesting articles have been published, e.g. the one by Carsten Prohl and Klaus Schmidt (here) or the one by Michael Merz and Mario Wuthrich (there).

But I think another interesting perspective (so far, not in claims reserving, but one should find some time to look at it) should be about multivariate regression (multivariate GLM’s), e.g.

All that will be mentioned in the talk. Slides can be downloaded here,

The dataset used in the example can be obtained with the code below

> P.corp=read.table("http://freakonometrics.blog.free.fr/public/data/auto-corporel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.corp=as.matrix(P.corp) > n=nrow(P.corp) > P.mat =read.table("http://freakonometrics.blog.free.fr/public/data/auto-materiel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.mat=as.matrix(P.mat) > P.mat=P.mat[1:n,1:n] >  P.mat = P.mat[2:10,1:9] >  P.corp= P.corp[2:10,1:9] > n=9 > P.tot = P.mat + P.corp

Estimateurs des moments

Dans le cours, nous avions vu que si http://freakonometrics.hypotheses.org/files/2015/12/moment-est.gif, un estimateur naturel pour le paramètre était http://freakonometrics.hypotheses.org/files/2015/12/moment-est2.gif. On appelle cette méthode la méthode des moments. Mais il est possible, que, plus généralement http://freakonometrics.hypotheses.org/files/2015/12/mmomot3.gif, et on pourrait considérer comme estimateur naturel

http://freakonometrics.hypotheses.org/files/2015/12/kvmoment4.gif

On avait noté en cours que le premier soucis était un soucis d’inversibilité. Pour utiliser cette méthode, il faut une relation bijective entre le moment et le paramètre. On avait vu avec la loi binomiale – que ce n’est pas le cas pour la variance, par exemple.

Le second soucis est un problème numérique: il faut inverser cette fonction (ce qui peut ne pas être trivial) mais avec une méthode de type sécante (programméeici), ça pourrait marcher. Par exemple pour la loi de Poisson, on sait que

http://freakonometrics.hypotheses.org/files/2015/12/poisson-mu-1.gif
http://freakonometrics.hypotheses.org/files/2015/12/poisson-m-2.gif
http://freakonometrics.hypotheses.org/files/2017/10/poisson-m-3.gif
http://freakonometrics.hypotheses.org/files/2015/12/poisson-m4.gif

(…etc). On peut numériquement construire des estimateurs de la manière suivante

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
}}
X=rpois(50,10)
m1=mean(X)
f=function(x){x*(1+x)-mean(X^2)}
secant(f,0,40)
m2=secant(f,0,50)
f=function(x){x*(1+3*x+x^2)-mean(X^3)}
secant(f,0,50)
m3=secant(f,0,50)
f=function(x){x*(1+7*x+6*x^2+x^3)-mean(X^4)}
secant(f,0,100)
m4=secant(f,0,20)

Le dernier soucis est un problème plus statistique: il semble que l’estimateur basé sur le premier moment soit meilleur que les autres. Sur le dessin ci-dessous, on a simulé 10,000 fois des lois de Poisson de paramètre 10 (qui est la vraie valeur du paramètre, ici représentée verticalement en pointillés), et on représente en rouge l’estimateur construit sur le premier moment, avec la croix au centre qui représente la valeur moyenne, et le segment qui est un intervalle de confiance à 80% de notre estimateur (i.e. les quantiles à 10% et 90% obtenus sur les 10,000 simulations),

En bleu, on a l’estimateur construit à l’aide du moment d’ordre deux, en vert du moment d’ordre trois, et en mauve du moment d’ordre quatre. Bref, plus on prend des moments d’ordre élevé, plus le biais semble important, et plus la variance de notre estimateur augmente.

Moralité, l’estimateur de la méthode des moments (et l’unique) sera construit sur les moments les plus faibles possibles…

Open data might be a false good opportunity…

I am always surprised to see many people on Twitter tweeting about #opendata, e.g. @data4all, @usdatagov, @datapublicatwit, @ProPublica or @open3 among so many others… Initially, I was also very enthousiastic, but I have to admit thatopen data are rarely raw data. Which is what I am usually looking for, as a statistician…
Consider the following example: I was wondering (Valentine’s day is approaching)when will a man born in 1975 (say) get married – if he ever gets married ?More technically, I was looking for a distribution of the age of first marriage (given the year of birth), including the proportion of men that will never get married, for that specific cohort.

The only data I found on the internet is the following, on statistics.gov.uk/

Note that we can also focus on women (e.g. here). Is it possible to use that opendata to get an estimation of the distribution of first marriage for some specific cohort ? (and to answer the question I asked). Here, we have two dimensions: on line http://freakonometrics.free.fr/blog/latex/marriage01.gif, the year (of the marriage), and on column http://freakonometrics.free.fr/blog/latex/marriage02.gif, the age of the man when he gets married. Assume that those were rawdata, i.e. that we have the number of marriages of men of age http://freakonometrics.free.fr/blog/latex/marriage02.gif during the year http://freakonometrics.free.fr/blog/latex/marriage01.gif.

We are interested at a longitudinal lecture of the table, i.e. consider some man born year http://freakonometrics.free.fr/blog/latex/marriage03.gif, we want to estimate (or predict) the age he will get married, if he gets married. With raw data, we can do it… The first step is to build up triangles (to have a cohort vs. age lecture of the data), and then to consider a model, e.g.

http://freakonometrics.free.fr/blog/latex/marriage04.gif

where http://freakonometrics.free.fr/blog/latex/marriage05.gif is a year effect, and http://freakonometrics.free.fr/blog/latex/marriage06.gif is a cohort effect.

base=read.table("http://freakonometrics.free.fr/mariage-age-uk.csv",
sep=";",header=TRUE)
m=base[1:16,]
m=m[,3:10]
m=as.matrix(m)
triangle=matrix(NA,nrow(m),ncol(m))
n=ncol(m)
for(i in 1:16){
triangle[i,]=diag(m[i-1+(1:n),])
}
triangle[nrow(m),1]=m[nrow(m),1]
 
triangle
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]   12  104  222  247  198  132   51   34
 [2,]    8   89  228  257  202  102   75   49
 [3,]    4   80  209  247  168  129   92   50
 [4,]    4   73  196  236  181  140   88   45
 [5,]    3   78  242  206  161  114   68   47
 [6,]   11  150  223  199  157  105   73   39
 [7,]   12  117  194  183  136   96   61   36
 [8,]   11  118  202  175  122   92   62   40
 [9,]   15  147  218  162  127   98   72   48
[10,]   20  185  204  171  138  112   82   NA
[11,]   31  197  240  209  172  138   NA   NA
[12,]   34  196  233  202  169   NA   NA   NA
[13,]   35  166  210  199   NA   NA   NA   NA
[14,]   26  139  210   NA   NA   NA   NA   NA
[15,]   18  104   NA   NA   NA   NA   NA   NA
[16,]   10   NA   NA   NA   NA   NA   NA   NA
 
Y=as.vector(triangle)
YEARS=seq(1918,1993,by=5)
AGES=seq(22,57,by=5)
X1=rep(YEARS,length(AGES))
X2=rep(AGES,each=length(YEARS))
reg=glm(Y~as.factor(X1)+as.factor(X2),family="poisson")
summary(reg)
 
Call:
glm(formula = Y ~ as.factor(X1) + as.factor(X2), family = "poisson")
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.4502  -1.1611  -0.0603   1.0471   4.6214  
 
Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        2.8300461  0.0712160  39.739  < 2e-16 ***
as.factor(X1)1923  0.0099503  0.0446105   0.223 0.823497    
as.factor(X1)1928 -0.0212236  0.0449605  -0.472 0.636891    
as.factor(X1)1933 -0.0377019  0.0451489  -0.835 0.403686    
as.factor(X1)1938 -0.0844692  0.0456962  -1.848 0.064531 .  
as.factor(X1)1943 -0.0439519  0.0452209  -0.972 0.331082    
as.factor(X1)1948 -0.1803236  0.0468786  -3.847 0.000120 ***
as.factor(X1)1953 -0.1960149  0.0470802  -4.163 3.14e-05 ***
as.factor(X1)1958 -0.1199103  0.0461237  -2.600 0.009329 ** 
as.factor(X1)1963 -0.0446620  0.0458508  -0.974 0.330020    
as.factor(X1)1968  0.1192561  0.0450437   2.648 0.008107 ** 
as.factor(X1)1973  0.0985671  0.0472460   2.086 0.036956 *  
as.factor(X1)1978  0.0356199  0.0520094   0.685 0.493423    
as.factor(X1)1983  0.0004365  0.0617191   0.007 0.994357    
as.factor(X1)1988 -0.2191428  0.0981189  -2.233 0.025520 *  
as.factor(X1)1993 -0.5274610  0.3241477  -1.627 0.103689    
as.factor(X2)27    2.0748202  0.0679193  30.548  < 2e-16 ***
as.factor(X2)32    2.5768802  0.0667480  38.606  < 2e-16 ***
as.factor(X2)37    2.5350787  0.0671736  37.739  < 2e-16 ***
as.factor(X2)42    2.2883203  0.0683441  33.482  < 2e-16 ***
as.factor(X2)47    1.9601540  0.0704276  27.832  < 2e-16 ***
as.factor(X2)52    1.5216903  0.0745623  20.408  < 2e-16 ***
as.factor(X2)57    1.0060665  0.0822708  12.229  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
(Dispersion parameter for poisson family taken to be 1)
 
    Null deviance: 5299.30  on 99  degrees of freedom
Residual deviance:  375.53  on 77  degrees of freedom
  (28 observations deleted due to missingness)
AIC: 1052.1
 
Number of Fisher Scoring iterations: 5

Here, we have been able to derive http://freakonometrics.free.fr/blog/latex/marriage12.gif and http://freakonometrics.free.fr/blog/latex/marriage13.gif, where now http://freakonometrics.free.fr/blog/latex/marriage14.gifdenotes the cohort.
We can now predict the number of marriages per year, and per cohort

http://freakonometrics.free.fr/blog/latex/marriage15.gif

Here, given the cohort http://freakonometrics.free.fr/blog/latex/marriage03.gif, the shape of http://freakonometrics.free.fr/blog/latex/marriage16.gif is the following

Yp=predict(reg,type="response")
tYp=matrix(Yp,nrow(m),ncol(m))
tYp[16,]
tYp[16,]
[1]  10.00000 222.94525 209.32773 159.87855 115.06971  42.59102
[7]  18.70168 148.92360
The errors (Pearson error) look like that
Ep=residuals(reg,type="pearson")
 

(where the darker the blue, the smaller the residuals, and the darker the red, the higher the residuals). Obviously, we are missing something here, like a diagonal effect. But this is not the main problem here…

I guess that study here is not valid. The problem is that we deal with open data, and numbers of marriages are not given here: what is given is a he proportion of marriage of men of age http://freakonometrics.free.fr/blog/latex/marriage02.gif during the year http://freakonometrics.free.fr/blog/latex/marriage01.gif, with a yearly normalization. There is a constraint on lines, i.e. we observe

http://freakonometrics.free.fr/blog/latex/marriage08.gif

so that

http://freakonometrics.free.fr/blog/latex/marriage09.gif

This is mentioned in the title

It is still possible to consider a Poisson regression on the http://freakonometrics.free.fr/blog/latex/marriage10.gif, but unfortunately, I do not think any interpretation is valid (unless demography did not change last century). For instance, the following sum

http://freakonometrics.free.fr/blog/latex/marriage17.gif

looks like that

apply(tYp,1,sum)
 [1] 919.948 838.762 846.301 816.552 943.559 930.280 857.871 896.113
 [9] 905.086 948.087 895.862 853.738 826.003 816.192 813.974 927.437

i.e. if we look at the graph

But I do not think we can interpret that sum as the probability (if we divide by 1,000) that a man in that cohort gets married…. And more basically, I cannot do anything with that dataset…

So open data might be interesting. The problem is that most of the time, the data are somehow normalized (or aggregated). And then, it becomes difficult to use them…

So I will have to work further to be able to write something (mathematically valid) on marriage strategy before Valentine’s day…. to be continued.

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);

Méthodes numériques et convolution, suite

Pour poursuivre un ancien billet, je voulais revenir plus particulièrement sur l’algorithme de Panjer, et sur la précision de l’approximation faite en discrétisant la loi des coûts individuels. Les références sont ici

Je reprends ici un partie de l’article de wikipedia sur le sujet. Comme toujours, on souhaite approcher la loi de

 S = \sum_{i=1}^N X_i\,

On suppose ici que les coûts individuels sont discrets, et on pose

f_k = P[X_i = hk].\,

On suppose de plus que les nombres sont dans une classe un peu particulière, vérifiant

 P[N=k] = p_k=

En fait, cette condition n’est pas trop forte, la plupart des lois classiques étant de ce type.

La formule permettant d’obtenir

g_k =P[S = hk] \,.

est alors obtenue par récurrence, en posant

g_0=p_0\cdot \exp(f_0 b)\text{ if }a = 0,\,

g_0=\frac{p_0}{(1-f_0a)^{1+b/a}}\text{ for }a \ne 0,\,

et ensuite

g_k=\frac{1}{1-f_0a}\sum_{j=1}^k \left( a+\frac{b\cdot j}{k} \right) \cdot f_j \cdot g_{k-j}.\,

L’avantage c’est que tout cela a été programmé sous R… Pour la discrétisation, je reprends ici actuar qui propose plusieurs techniques. En particulier on distinguera entre une borne supérieure et une borne inférieure, une version arrondir et une version sans biais.

> fx <- discretize(pgamma(x, 2, 1), method = "upper", from = 0, to = 17, step = 0.5)
> fx <- discretize(pgamma(x, 2, 1), method = "lower", from = 0, to = 17, step = 0.5)
> fx <- discretize(pgamma(x, 2, 1), method = "rouding", from = 0, to = 17, step = 0.5)
> fx <- discretize(pgamma(x, 2, 1),method = "unbiased", lev = levgamma(x, 1), from = 0, to = 17, step = 0.5)

On peut alors utiliser l’algorithme de Panjer pour calculer la loi composée basée sur ces discrétisations. Au niveau de la programmation, le code est le suivant 

>  fx <- discretize(pgamma(x, 2, 1), from = 0, to = 22,
+  step = 0.5, method = "unbiased", lev = levgamma(x,2, 1))
>  Fs <- aggregateDist("recursive", model.freq = "poisson",
+  model.sev = fx, lambda = 10, x.scale = 0.5)

Graphiquement, on obtient la distribution suivante pour la loi composée (fonction de répartition et fonction quantile).

Simulation d’un processus de Poisson

Le prochain cours reviendra sur la solvabilité des compagnies d’assurance, et en particulier sur la modélisation et le calcul de la probabilité de ruine. En particulier, nous passerons pas mal de temps sur le “processus de Poisson“.

Nous verrons en cours les notions de base sur les processus de Poisson$ \{N_t\}$. En particulier, nous verrons qu’ils se caractérisent soit par le processus d’arrivées $ \{T_i\}$, correspondant aux dates d’arrivée des sinistres, ou bien par$ W_i=, correspondants aux durées intersinistres. En effet, on peut obtenir des expressions de la forme $ N_t=.

Dans sa forme la plus simple, on s’intéresse aux processus homogènes. $ \{N_t: t \ge 0 \}$ sera un processus de Poisson homogène d’intensité $ \lambda > 0$ si les durées entre arrivées, $ W_i$ i.i.d. de loi exponentielle de paramètre $ \lambda $. L’algorithme naturel pour simuler les dates d’arrivées de sinistre $ T_1,T_2,\ldots,T_n$ est alors

  1. poser $ T_0 =
  2. pour $ i = on génère$ E$ suivant une loi exponentielle de paramètre $ \lambda $, puis on pose $ T_i =

Les choses vont se compliquer si on suppose que le processus de Poisson n’est plus homogène. Cette hypothèse est peut être davantage réaliste pour l’assurance. Par exemple, les assureurs savent qu’il y a plus d’accidents automobiles les jours de brouillard, ou de pluie, ou de verglas. On peut donc supposer qu’il existe un cycle lié au climat.

On suppose ici que le processus de Poisson est d’intensité $ \lambda(t)$. La première idée pour simuler un tel processus est de supposer qu’il existe $ \overline{\lambda}$ tel que $ \lambda(t)\leq\overline{\lambda}$ pour tout $ t$. On note $ T_1^* ,T_2^* , T_3^* ,\ldots$ des arrivées successives d’un processus de Poisson de paramètre $ \overline{\lambda}$. On met alors en place une stratégie acceptation/rejet, telle que l’on garde le $ i$ème sinistre $ T_i^*$ avec probabilité $ \lambda(T_i^* )/\overline{\lambda}$. La suite $ T_1, T_2, \ldots$ ainsi construite suit un processus de Poisson nonhomogène d’intensité $ \lambda(t)$. D’où l’algorithme

  1. On pose $ T_0 = et $ T^* = 0$
  2. on génère une variable exponentielle$ E$ d’intensité $ \overline{\lambda}$.
  3. On pose alors$ T^* = T^* + E$.
  4. On génère une variable $ U$ uniforme sur $ (0,1)$: si $ U>\lambda(T^*)/\overline{\lambda}$ on rejette (et on retourne en 2.), et sinon on pose $ T_i = .

Une autre méthode consiste à noter que l’incrément

$ N_t-N_s$, $ 0 < s < t$, suit une loi de Poisson d’intensité$ \widetilde{\lambda} = \int_s^t \lambda(u)du$. Aussi, la fonction de répartition $ F_s$ du temps d’attente $ W_s$ est qui vaut

On peut utiliser la méthode suivante pour simuler le processus de Poisson

  1. Poser $ T_0 =
  2. Générer $ U$ suivant une loi uniforme sur $ (0,1)$
  3. Poser $ T_i =

Pour aller plus loin, je renvoie à des bouquins sur la simulation de processus, ou pour les applications assurantielles au Statistical Tools for Finance and Insurance, chez Springer.