C’est agaçant cette manie de battre des records !

En janvier dernier, j’avais fait un billet pour noter que 2010 était une année record pour les catastrophes naturelles, mais que 2009 aussi, et 2008 aussi. Bref, tous les ans, on bat des records. Mais il n’y a pas que sur le climat que l’on bat sans cesse des records. En sport aussi.
Par exemple, dimanche, le record du monde du marathon a (encore) été battu, avec 42,195 kilomètres parcourus (ou courus, tout simplement) en 2:03:58 2:03:38 par Patrick Makau (bon, certes, avec une amélioration d’une seconde par rapport au précédant record). On devrait être content (en tous les cas Xi’an semble l’être). Sauf qu’à l’hiver 2012, je vais donner un cours sur copules et valeurs extrêmes à l’UQÀM, alias MAT8886. Et je pensais faire travailler les étudiants sur les articles de John Einmahl, sur les records au 100 mètres, ou plus généralement les records en athlétisme. Ces papiers proposent une très jolie application de la théorie des valeurs extrêmes, en particulier si on est dans le domaine d’attraction de la loi de Weibull (on peut alors chercher la borne supérieur – si on travaille sur le max – du support de la distribution, appelé ici endpoint). On apprend ainsi que selon la théorie des valeurs extrêmes, le temps minimal pour parcourir les 42,195 km du marathon devrait être de 2:04:06.

I.e. le temps pour courir un marathon devrait être borné par cette valeur. Sauf que ce record a été battu. Et à plusieurs reprises me semble-t-il.
Pour aller un peu plus loin, et comprendre ce qui se passait, j’ai été chercher les données sur wikipedia, pour le marathon de Boston, Chicago, Paris, Berlin, NYC, Stockholm, Fukuoka, Rotterdam, Amsterdam ou encore Londres. J’ai ainsi le temps du vainqueur pour les dernières années. Les données ont été copiées puis collées et donc un peu de remise en forme est ici nécessaire,

> base=read.table("http://freakonometrics.blog.free.fr/public
/data/topmarathon.csv",sep=";",header=TRUE)
> base=base[is.na(base$TIME)==FALSE,]
> base=base[(base$TIME=="")==FALSE,]
> n=nchar(as.character(base$YEAR))
> base$Y=as.numeric(substr(as.character(base$YEAR),n-3,n))
> base$Y[is.na(base$Y)]=1987
> h1=as.numeric(substr(as.character(base$TIME),1,1))
> h2=as.numeric(substr(as.character(base$TIME),1,2))
> i =which(is.na(t2))
> h=h2; h[i]=h1[i]
> m1=as.numeric(substr(as.character(base$TIME),3,4))
> m2=as.numeric(substr(as.character(base$TIME),4,5))
> m=m2; m[i]=m1[i]
> s1=as.numeric(substr(as.character(base$TIME),6,7))
> s2=as.numeric(substr(as.character(base$TIME),7,8))
> s=s2; s[i]=s1[i]
> base$T=h+m/60+s/60/60
> base=base[base$T>0,]
> base0=base[base$Y>=1982,]
> plot(base0$Y,base0$T,xlab="",ylab="")
> library(splines)
> reg=lm(T~bs(Y,6),data=base0)
> lines(1983:2010,predict(reg,newdata=
+ data.frame(Y=1983:2010)),col="red",lwd=2)

Si on regarde attentivement on observe plusieurs choses intéressantes. La première est qu’on ne court pas à la même vitesse dans toutes les villes. Par exemple à Stockholm, le temps du plus rapide est souvent bien éloigné des temps records. Mais surtout, le temps moyen du vainqueur ne cesse de baisser avec le temps,

http://freakonometrics.hypotheses.org/files/2015/12/marathon-cities_m.gif

On a alors une forme de non-stationnarité de la série qui laisse à penser qu’il conviendrait d’étudier davantage cette non-stationnarité avant d’utiliser la théorie des valeurs extrêmes. Et cette non-stationnarité et l’analyse des records n’est pas sans rappeler la discussion que l’on avait sur les catastrophes naturelles. Il serait peut-être temps de creuser davantage…

Visualisation en tarification, avec R

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

  • Quelques références en R

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

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

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

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

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

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

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

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

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

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

> graphique()

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

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

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

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

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

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

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

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

Si on regarde par classe d’age,

> graphiquecout()

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

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

ou enfin, en fonction de la zone géographique

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

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…

ACT2040, exercices de tarification

Les exercices commenceront ce vendredi, et les énoncés (ainsi que pour quelques uns des pistes de correction) sont en ligne ici. Les exercices (et les corrections) ont été tapés pour ce cours par Jean-Philippe Boucher. La premiere série d’exercices sont des exercices de remise en route sur les calculs d’espérances, et d’espérance conditionnelle. On reviendra dessus mardi prochain.