SOA Webinar on Predictive Modeling

I will give, with Qichun Xu, a joint webinar for the Reinsurance Council and the Futurism Council of the Society of Actuaries, on Perspectives of Predictive Modeling with Case Studies in a few days. The slides of my talk are now available (I do recommand to open the pdf version of the slides with Acrobat, since there are animated pictures in the slides that could not be visualized below for instance). The Society of Actuaries asked specifically for a powerpoint document, so I will use screenshots of the slides for the webinar. I do encourage to open and read the pdf file for a better quality… Sorry for the inconvenience. I will upload soon lines of codes to reproduce most of the graphs. All comments and remarks are welcome.

Halloween and candies (a ballot problem)

This year, for Halloween, a post on candies (I promise, next year I will write another post on zombies). But I don’t want to focus on the kids problems (last year, we tried to minimize their walking distance to collect as much candies as possible, with part 1 and part 2), I want to discuss my own problems. Because usually, the kids wear their costumes, and they go in the streets, they knock on the doors, while I stay at home. So I’m the one, with a bag full of candies, waiting for kids to knock on our door, and then I give them some candies (if they wear a costume). Consider the following problem. Assume that we start with  red candies, and  black ones, with . The thing is that no one like those black candies. What could be the probability that for the   kids that will get candies after knocking at my door (with  for convenience, but we will also consider the more general case where I have to many candies, , later on), the probability to get a red candy is always larger than the probability to have a black candy ? This is somehow related to the popular ballot problem, proposed (and solved) by Whitworth in 1878, but he wrote it only in the fourth edition of Choice and Chance, in 1886 (this is what the legend told us). In 1887, Joseph Bertrand proposed a similar problem, and Désiré André introduced the reflection principle to solve it. The problem is simple : consider an election between two candidates, A (who receives  votes) and B (who receives  votes). A wins the election (). If the ballots are cast one at a time, what is the probability that A will lead throughout the voting? For those who don’t remember the conclusion, the probability is here quite simple,

Observe that some geometry proofs were given, later on, by Aebly or Mirimanoff, both in 1923, as well as Howard Grossman in the 1950’s (see the discussion on…).  Actually, produced the following geometric proof (with no clear reference),

ballot box lattice diagram

We start at O, where no votes have been cast. Each vote for A moves us one point east and each vote for B moves us one point north until we arrive at E, the final count, (mn). If A is to lead throughout the contest, then our path must steer consistently east of the diagonal line OD, which represents a tie score. Any path that starts by going north, through (0,1), must cut OD on its way to E.

If any path does touch OD, let it be at C. The group of such paths can be paired off as p and q, reflections of each other in the line OD that meet at C and continue on a common track to E.

This means that the total number of paths that touch OD is twice the number of paths p that start their journey to E by going north. Now, the first segment of any path might be up to m units east or up to units north, so the proportion of paths that start by going north is n/(m + n), and twice this number is 2n/(m + n). The complementary probability — the probability of a path not touching OD — is (m –n)/(m + n).

But let’s try to solve our problem. Let  and  denote the number of black and red candies, respectively after the th kid git his (or her) candy. Yes, one at a time. Here,  and . What we want is

Using this formulation, we recognize the ballot problem. Almost. Actually, in the original ballot problem (see Bertrand (1887)), we have to compute the probability that one candidate remains strictly ahead the other one throughout the count. With a strict condition, we get the well-known probability (given previously)

Here, ties are allowed, and we can prove (easily) that

(again, there is some nice geometric interpenetration of that result). It is also possible to get numerically that value using the following function, which will generate a trajectory, and return some indicators (with or without ties)

> red_black=function(sd){
+ set.seed(sd)
+ vectcandy=sample(c(rep("R",r),rep("B",b)))
+ v1=rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B"))
+ v2=cumsum(rev(vectcandy)=="R")<= cumsum(rev(vectcandy)=="B") 
+ return(list(evol=cbind(rev(cumsum(rev(vectcandy)=="R")),
+ rev(cumsum(rev(vectcandy)=="B")),v1),list=vectcandy,test=(sum(v1)==0),
+ ballot=(sum(v2)==0),when=min(which(v1==1))))}

(here I compute the ballot-type index, where ties are not allowed, and the candy-type index). If we generate 100,000 scenarios, starting with 50 red and 25 black candies, we get

> r=50 
> b=25 
> M=sapply(1:100000,red_black) ­­
> mean(unlist(M[3,]))
[1] 0.50967

which can be compared with the theoretical value

> (r+1-b)/(r+1)
[1] 0.5098039

We can also get the distribution of the first time we have more black candies than red ones left (given that this event occur)

> Z=unlist(M[5,])
> Z=Z[Z<Inf]
> hist(Z,breaks=seq(0,80),probability=TRUE,col="light blue", 
+ border=NA,xlab="",main="")

There might be some analytically formula that can be derived, but I have to confess that I am becoming extremely lazy,

Assume now that this year, kids do not show up at my door (for some reason). Assume that  kids show up. We can see how the probability

will change, with ,

> r=50
> b=25
> impact_n = function(n){
+ red_black=function(sd,nb=n){
+ set.seed(sd)
+ vectcandy=sample(c(rep("R",r),rep("B",b)))
+ v=(rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")))[1:nb]
+ return(list(list=vectcandy,test=(sum(v)==0),when=min(which(v==1))))}
+ M=sapply(1:10000,red_black)
+ return(mean(unlist(M[2,])))}

Yes, not only I am too lazy to derive analytic formulas, I am so lazy that I do not try to optimize my code. Here, the evolution of the probability, as a function of  is

> V=Vectorize(impact_n)(25:75)
> plot(25:75,V)

Fun isn’t it? But now, I have to conclude my post, to work a little bit on my make-up : I have learnt so many thinks at the Montreal Zombie Walk a few days ago that kids willing to knock at my door will be scared to death. I guess I will keep all the candies for me this year !

More significant? so what…

Following my non-life insurance class, this morning, I had an interesting question from a student, that I will try to illustrate, and reformulate as accurately as possible. Consider a simple regression model, with one variable of interest, and one possible explanatory variable. Assume that we have two possible models, with the following output (yes, I do hide interesting parts here, but it is to get quickly to my student’s point)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.92883    0.06391  14.534   <2e-16 ***
X           -0.12499    0.06108  -2.046   0.0421 *  
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

for the first model – a GLM with some distribution, and some link function – and

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.92901    0.06270  14.817   <2e-16 ***
X           -0.09883    0.05816  -1.699   0.0909 .  
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

for the second one – with another GLM, with another distribution, but the same link function (I guess I could have changed it, but it does not really matter here). Then, I got the following statement “I would like to choose the first model because the explanatory variable is more significant, and therefore, this model should have a stronger predictive power“.

That’s a nice idea, isn’t it ? Actually, I guess this is why I love teaching, because I will never be able to think about such an idea by myself. Because when you look at that statement, somehow it could make sense. Except that from my point of view, it is not valid at all. My first thought was to recall is standard example in statistical inference : you cannot not claim that a distribution is better than another one just by looking at the parameter estimates.

> fitdistr(Y,"normal")
      mean          sd    
  0.93685011   0.90700830 
 (0.06413517) (0.04535042)
> fitdistr(Y,"exponential")

Can I claim that the Gaussian distribution is better than the exponential one because parameter estimates have smaller standard deviation ? Because somehow, this is what we did when we claimed previously that the first model was better than the second one.

Let me get back on the outputs of the two regressions, and let me explain what I did. Actually, I wanted to have a story close to the one on the Gaussian versus exponential fit. So I did generate some exponential random variable,

> set.seed(5)
> n=200
> U=runif(n); 
> Y=-log(U)

Here, we can visualize the histogram of this sample, as well as the the estimated exponential distribution

> hist(Y,proba=TRUE,col="light green",border="white",lwd=2,breaks=seq(0,5.3333333333333,by=.333333333))
> x=seq(0,6,by=.02)
> lines(x,dexp(x,1/mean(Y)),col="red",lty=2)

On top of that, let us fit a gamma distribution. Using a GLM (where the regression is here on a constant – only), just to practice because later on, we will use a gamma regression on that variable

> reg0=glm(Y~1,family=Gamma(link="identity"))
> a=reg0$coefficient
> b=summary(reg0)$dispersion
> lines(x,dgamma(x,shape=1/b,scale=a*b),col="blue")

Now, we need a covariate, to run some regressions. What I wanted is some variable slightly correlated with our previous variable. Slightly, just to make sure that our -value in the regression will be close to 5% or 10%. So here, I did generate a variable so that the pair has Clayton copula, with coefficient 0.1 (which is small, extremely small)

> a=.1
> set.seed(5)
> n=200
> U=runif(n); 
> V=(U^(-a)*(runif(n)^(-a/(1+a))-1)+1)^(-1/a)
> Y=-log(U)
> X=qnorm(V)

To visualize the copula of the variables, we can use

> cop=function(u,v){
+ (a+1)*(u*v)^(-(a+1))*
+ (u^(-a)+v^(-a)-1)^(-(2*a+1)/a) }
> x=y=seq(.05,.95,by=.05)
> z=outer(x,y,cop)
> mat=persp(x,y,z,col="green",shade=TRUE,xlim=c(0,1),ylim=c(0,1),zlim=c(0,2),theta=-30,
+ ticktype ="detailed",zlab="")

We should be not far away from the independence (actually, there is a negative – significant – correlation (Pearson’s correlation)). Now, consider two models,

  • a Gaussian model (here a standard linear model)
  • a gamma model, with a linear link function

The outputs are the following (you will recognize the outputs given previously)

> reg1=lm(Y~X)
> reg2=glm(Y~X,family=Gamma(link="identity"))
> summary(reg1)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.92883    0.06391  14.534   <2e-16 ***
X           -0.12499    0.06108  -2.046   0.0421 *  
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9021 on 198 degrees of freedom
Multiple R-squared:  0.02071,	Adjusted R-squared:  0.01576 
F-statistic: 4.187 on 1 and 198 DF,  p-value: 0.04206

> summary(reg2)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.92901    0.06270  14.817   <2e-16 ***
X           -0.09883    0.05816  -1.699   0.0909 .  
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.9086447)

    Null deviance: 229.72  on 199  degrees of freedom
Residual deviance: 226.58  on 198  degrees of freedom
AIC: 379.22

Number of Fisher Scoring iterations: 10

And here are the two predictions,

So, which model should we use? As usual, my answer will be “let’s have a look at the data” instead of looking only at tables of figures. Using some code posted a few days ago, let us visualize the two regressions. The Gaussian model is here

(for the lower part, I do not go below 0 since we do have, here, a positive variable that we would like to model) while the gamma on is here

And if we believe that the explanatory variable has no predictive power (since we can claim that the parameter is not significant in the regression), and we remove it from the regression, we get

Here, I do believe that the gamma (not to say the exponential) model is better because it is clearly more coherent with properties of the variable of interest. I trust more the confidence interval obtained above on the gamma model, than the one obtained with a Gaussian distribution. Even if the parameter in the regression is “more significant”.

Moments et variables aléatoires

Vendredi, suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur les thèmes 10, 11, 12 et 16 (tels que classifiés dans le livre de Jacques Labelle, qui sert de référence pour ce cours)

Somewhere else, part 86

Some writings worth reading,

…journalism has adopted the academic publishing model, only without the pretense of integrity. The 2008 economic crisis, combined with the transition to digital media, led to a glut of desperate writers willing to work for free—a practice that media corporations embraced and repackaged to novice journalists as “the way things have always been.” Today media outlets making healthy profits refuse to pay the freelance writers who help make them a success. Exploitative publishers tend to argue along two lines: a fake crisis (“Unfortunately, we can’t afford to pay you at this time…”) or a false promise (“Exposure will help your career.”). Academics are particularly vulnerable to media-industry exploitation. They are accustomed to writing for nothing and, in the case of adjuncts, to being treated terribly by their employers. Because academic work in professional journals is hidden behind paywalls, the prospect of reaching a wider audience can be enticing. For scholars interested in leaving academia and forging a new career, online visibility is essential. Should academics ever write for free? Maybe. Should academics write for free for a publisher that can afford to pay them? Never. [to be continued…]

In the graphic of one person’s network neighborhood (above), the cluster at the top is the individual’s co-workers. The cluster at the right is old college friends. The node (friend) in the lower left quadrant of the graphic, with links to the two dense clusters — but at a distance from those clusters — is the user’s spouse. “A spouse or romantic partner is a bridge between a person’s different social worlds,” Mr. Kleinberg explained in an interview on Sunday. Their dispersion algorithm was able to correctly identify a user’s spouse 60percent of the time, or better than a 1-in-2 chance. Since everyone in the sample had at least 50 friends, merely guessing would have at best produced a 1 in 50 chance. The algorithm also did pretty well with people who declare themselves to be “in a relationship,” correctly identifying them a third of the time — a 1 in 3 chance compared with the 1 in 50 for guesswork. Particularly intriguing is that when the algorithm fails, it looks as if the relationship is in trouble. A couple in a declared relationship and without a high dispersion on the site are 50 percent more likely to break up over the next two months than a couple with a high dispersion, the researchers found. (Their research tracked the users every two months for two years.) [to be continued…]

“I’m just not a math person.” We hear it all the time. And we’ve had enough. Because we believe that the idea of “math people” is the most self-destructive idea in America today. The truth is, you probably are a math person, and by thinking otherwise, you are possibly hamstringing your own career. Worse, you may be helping to perpetuate a pernicious myth that is harming underprivileged children—the myth of inborn genetic math ability. Is math ability genetic? Sure, to some degree. Terence Tao, UCLA’s famous virtuoso mathematician, publishes dozens of papers in top journals every year, and is sought out by researchers around the world to help with the hardest parts of their theories. Essentially none of us could ever be as good at math as Terence Tao, no matter how hard we tried or how well we were taught. But here’s the thing: We don’t have to! For high school math, inborn talent is just much less important than hard work, preparation, and self-confidence. […]  Too many Americans go through life terrified of equations and mathematical symbols. We think what many of them are afraid of is “proving” themselves to be genetically inferior by failing to instantly comprehend the equations (when, of course, in reality, even a math professor would have to read closely). So they recoil from anything that looks like math, protesting: “I’m not a math person.” And so they exclude themselves from quite a few lucrative career opportunities. We believe that this has to stop. Our view is shared by economist and writer Allison Schrager, who has written two wonderful columns in Quartz (here and here), that echo many of our views. [to be continued…]

et un peu de lecture en français, à commencer par le brillant

mais cette semaine, davantage de lecture en français,

Did I miss something?


Mercredi, on finira la modélisation des coûts individuels de sinistres en évoquant la mutualisation. Si on a le temps, on parlera aussi de réassurance. Les transparents sont en ligne.

Sinon, histoire d’illustrer les aspects pratiques de la tarification, j’utiliserais peut-être la base xls des gros sinistres en perte d’exploitation, en France, sur la période 1985-2000. Côté lectures complémentaires, je recommande la lecture de Introduction à la réassurance, publié par Swiss Re, ou ainsi que quelques documents plus techniques, comme The Pareto model in property reinsurance , Exposure rating, ou Designing property reinsurance programmes encore Introduction to reinsurance accounting. Plusieurs réassureurs (et courtiers en réassurance) publient des études techniques sur leurs sites, ou encore Sinon je renvois aux notes de cours de Peter Antal, quantitative methods in reinsurance.
Et histoire de mettre à jour mes transparents, les sinistres les plus chers, pour les compagnies d’assurance et de réassurance :… donne le classement suivant en dollars de 2010 (on pourra aussi consulter…)

  1. Hurricane Katrina (US, Bahamas, Cuba, Aug. 2005), $ 72.3 billion
  2. Tōhoku earthquake and tsunami (Japan, March 2011), $ 35 billion
  3. Hurricane Andrew (US, Bahamas, August 1992), $ 25 billion
  4. September 11 attacks (US) $ 23.1 billion
  5. Northridge earthquake (US) $ 20.6 billion
  6. Hurricane Ike (US, Haiti, Dominican Republic, Sept. 2005) $ 20.5 billion
  7. Hurricane Ivan (US, Barbados, Sept. 2004) $ 14.9 billion
  8. Hurrican Wilman (US, Mexico, Jamaica, Oct. 2005), $ 14 billion
  9. Hurricane Rita (US, Cuba, Sept. 2005) $ 11.3 billion
  10. Hurricane Charley (US, Cuba, Jamaica) $ 9.3. billion

A titre de comparaison, les chiffres d’affaires des plus gros réassureurs (prime émise en 2010) étaient, selon…

  1. Munich Reinsurance Company $ 31.3 billion
  2. Swiss Reinsurance Company Limited $ 24.7 billion
  3. Hannover Rueckversicherung AG $ 15.1 billion
  4. Berkshire Hathaway Inc. $ 14.4 billion
  5. Lloyd’s $ 13 billion
  6. SCOR S.E. $  8.8 billion
  7. Reinsurance Group of America Inc. $ 7.2 billion
  8. Allianz S.E. $ 5.7 billion
  9. PartnerRe Ltd. $ 4.9 billion
  10. Everest Re Group Ltd. $ 4.2 billion

De la significativité (statistique), suite

Suite à mon rapide billet sur l’article donnant les conclusions d’une étude des effets de la e-cigarette (et la comparaison avec le patch pour arrêter de fumer), j’avais entendu cette phrase “ça veut dire quoi ‘statistiquement non significatif’ ?“. Comme c’est une très bonne question, qui plus est relativement complexe, j’ai pris un peu de temps pour y répondre… Il ne s’agit pas de complexité mathématique, en s’entend car tous les cours de statistique abordent cette notion, mais plus d’une notion complexe à vulgariser. On va revenir un peu sur l’histoire des tests de significativité et sur les travaux de Ronald Fisher (pour comprendre l’origine des 5% utilisé comme seuil lors de l’utilisation de la p-value), mais il est possible d’aller voir directement les exemples.

  • Un peu d’histoire

Si on veut commencer par la préhistoire, les travaux les plus anciens qui posent la question de la significativité sont ceux de John Arbuthnot, en 1710. En prenant 82 années de données de natalité, à Londres, il avait observé que chacune des années, il y avait eu plus de naissances de garçons que de filles (évoqué dans Hacking (1975)).

Il s’était demandé, avec la terminologie d’il y a plus de 300 ans, si cette différence pouvait être jugée comme ‘statistiquement significative. Je pense que l’on peut utiliser le terme ‘statistiquement‘ car John Arbuthnot pose la question en termes probabilistes, très clairement,

Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un ‘test de significativité‘, là encore avec notre terminologie actuelle. Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste (que l’on verra formalisé bien plus tard), en se demandant s’il est ‘extrêmement probable‘ qu’il y ait une différence entre les deux mesures. Il avait alors fait un test de comparaison des valeurs moyennes, et noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important. A la fin du XIXème siècle, Francis Edgeworth reprendra cette idée (sur la comparaison de moyennes) lorsqu’il écrira en 1885,

où on retrouve les premières idées qui serviront de base aux travaux ultérieurs. Car pour bien comprendre l’origine de la notion de significativité en statistique, il faut relire les travaux de Ronald Fisher des années 20, mais aussi de Student (William Gosset) et Karl Pearson. Ronald Fisher s’interroge, dans son papier de 1921, sur la valeur d’un coefficient de corrélation (c’est dans ce papier qu’il introduit formelle sa transformation avec une tangente hyperbolique z = {1 \over 2}\ln\left({1+r \over 1-r}\right) = \text{arctanh}(r) pour faire un test de nullité d’un coefficient de corrélation), où il note

Aussi, dans le papier, une valeur (ici une corrélation) peut être jugée comme ‘significative‘ (c’est à dire ‘significativement non nulle‘). Encore une fois, l’idée n’est pas nouvelle, Edgeworth aurait dit ‘non accidentelle‘ mais Fisher pose ici la terminologie, en parlant de ‘significativité‘. L’année suivante, en 1922, il parle de l’idée de tester la significativité, dans le cadre d’un modele de régression (on est toujours dans le contexte de la corrélation) ou en 1924, on retrouve l’idée de tester la significativité, Mais c’est surtout dans son livre Statistical Methods for Research Workers que l’on va trouver une formalisation claire (et très illustrée) de cette idée, et surtout la connecter à la notion de p-value (on ne peut pas parler de significatité sans parler de p-value). Pour l’anecdote, dans l’ouvrage de Ronald Fisher, on y parle de p-value (oui, avec une majuscule). Cela dit, la paternité de la p-value est à attribuer davantage à Pearson, par exemple dans son papier datant de 1900, ou encore, sur un autre exemple avec en plus l’interprétation usuelle en terme de significativité (sans toutefois utiliser ce terme), puisque Pearson dit qu’on a 1 chance sur 8 d’avoir une valeur aussi improbable. Plus loin, il dit aussi qu’arriver avec 1 chance sur 70 peut être qualifier de très improbable (un long débat suivra sur le seuil à partir duquel, effectivement, un événement peut être qualifié d’improbable) Mais revenons au livre de Ronald Fisher, et plus spécifiquement au chapitre 4, qui propose précisément d’étudier la ‘significativité de la différence entre deux moyennes

On y retrouve le test de Pearson basé sur la statistique du chi-deux, comme test d’ajustement (goodness of fit) et d’indépendance. et la comparaison de la p-value avec (ce qui va correspondre aux mythiques) 5%. Il parle plus loin d’une ‘ligne conventionnelle‘ que l’on tracera à 5%, Il faut d’ailleurs un lien entre une p-value à 5% (correspondant à un événement qui survient avec 1 chance sur 20) et le fait de s’éloigner de 2 écart-types de la moyenne pour une loi normale (qui survient avec 1 chances sur 22) Tout au long du livre, on retrouve évoqué le chiffre mythique 5% qui restera à la postérité. Le plus intéressant pour apprendre comme juger la significativité d’une valeur, on peut regarder les exemples, comme l’exemple 11 l’exemple 12, l’exemple 27, ou encore l’exemple 32 (pour en prendre juste quelques uns) (on retrouve l’utilisation de 2 écart-types pour une loi normale). La discussion amorcée dans ce livre sera poursuivie dans son papier de 1926, qui commence avec Le papier est passionnant, car il pose les bases de ce que font tous les statisticiens depuis 80 ans (il n’est pas surprenant de voir cet article reproduit dans le second volume des Breakthroughs in Statistics), où le seuil de 5% est encore utilisé. Il ne cesse tout au long de l’article de justifier cette règle de ‘un sur vingt‘, tout en autorisant d’autres valeurs, comme 2%, ou 10% (Fisher semble moins orthodoxe que certains utilisateurs de la notion de significativité) Il faudra ensuite attendre les travaux de Jerzy Neyman et Egon Pearson. qui vont axiomatiser et formaliser les tests statistiques, mais ça dépasse (de loin) l’objectif de mon billet d’aujourd’hui. Donc on va repousser à plus tard la présentation de leurs apports. Pour illustrer tout ce qu’on vient de voir, considérons deux exemples simples.

  • Petit exemple, tester l’indépendance

Comme Ronald Fisher, commençons par ce qui me semble le plus simple car on n’a pas besoin de parler modèle probabiliste ici. On se demande si deux variables sont indépendantes, ou pas. Deux cas peuvent se présenter:

  • on a des variables continues, et on veut voir si elles sont non-corrélées (oui, pour commencer on va simplifier, et dire que si elles sont non-corrélées… ça va se traduire par une notion que l’on peut assimiler – rapidement – à une notion d’indépendance. Promis, d’ici 6 mois il y aura des dizaines de billets sur le sujet sur le blog puisque je redonne mon cours sur la dépendance cet hiver). Je laisse les plus motivés jouer avec les données qui traînent sur le blog pour voir si la taille et le poids sont corrélées, ou utiliser un exemple que j’avais l’habitude de présenter dans mes cours de modèles de régression, où on se demande si – chez les femmes – le tour de poitrine et le tour de taille sont corrélés.
  • on a des variables catégorielles, par exemple la couleur des yeux, et la couleur des cheveux. Et on se demande si les deux sont couleurs sont indépendantes, ou pas.

Regardons ce dernier exemple, avec un petit jeu de 25 observations.

> X=c(rep("Brun",8+5),rep("Blond",9+3))
> Y=c(rep("Noisette",8),rep("Bleu",5),rep("Noisette",3),rep("Bleu",9))

On peut résumer l’information dans un tableau de contingence,

> T=as.matrix(table(X,Y))
> T
X       Bleu Noisette
  Blond    9        3
  Brun     5        8

On va essayer de comparer ce tableau avec ce qu’on devait avoir si les deux variables étaient indépendantes (car c’est ce que l’on souhaite tester : on veut savoir si le tableau que l’on a est significativement différent de ce qu’on aurait avec des variables indépendantes). Il faut ici ressortir la définition de l’indépendance entre deux variables aléatoires, \mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B) = \mathbb{P}(X\in\mathcal A) \times\mathbb{P}(Y\in\mathcal B) pour tout \mathcal A et \mathcal B. Par exemple \mathcal A peut être avoir des cheveux blonds, et \mathcal B avoir des yeux bleus. Pour rappel, on avait en tout 25 personnes ayant ces deux caractéristiques en même temps, donc empiriquement, on avait observé une probabilité de 9/25 soit 36% (c’est le terme de gauche). Pour le terme de droite, on notera que \mathbb{P}(X\in\mathcal A), c’est la probabilité d’avoir les cheveux blonds. Empiriquement, on en avait 9+3 sur 25, soit 48%. Et \mathbb{P}(Y\in\mathcal B) c’est la probabilité d’avoir les yeux bleus. Or on avait 9+5 personnes avec ces caractéristiques, soit 56%. Si les deux variables étaient indépendantes, avec la formule précédente, la probabilité d’avoir à la fois les cheveux blonds et les yeux bleus serait 48% x 56%, soit 26.88%. A comparer aux 36% observés. En fait, au lieu de comparer des pourcentages, on peut aussi comparer des nombres: comme il y a 25 personnes en tout, 26.88% signifie que dans notre groupe, 6.72 personnes devraient avoir les cheveux blonds et les yeux bleus. On peut s’amuser à faire ces calculs pour tous les \mathcal A et \mathcal B.

> Tind=chisq.test(T)$expected
> Tind
X       Bleu Noisette
  Blond 6.72     5.28
  Brun  7.28     5.72

Visuellement, on peut mettre ça en forme,

(sur le graphique de droite, les traits noirs correspondent aux nombres sous hypothèse d’indépendance, et en couleur, les valeurs observées). Bon, maintenant, il faut réfléchir un peu… On se demande si nos deux tableaux sont proches ou pas. Ou, pour utiliser notre mot de la journée, significativement proches. Alors là, ça va devenir technique le temps d’un rapide paragraphe. Dans chaque case du tableau, on compte des gens. Sur nos n personnes, on peut s’attendre à en observer un nombre aléatoire, noté N_{\mathcal{A},\mathcal{B}}, qui suit une loi binomiale N_{\mathcal{A},\mathcal{B}}\sim\mathcal{B}(n,\mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B)) (c’est la définition de la loi binomiale). On notera que quand n est grand, la loi normale est alors proche de la loi Gaussienne. Dans chaque case, on a des lois normales (non indépendantes, car les nombres sont contraints, par exemple au total, il faut n personnes), et en prenant les carrés des lois normales (centrées et réduites), on va obtenir une loi du chi-deux (je renvoie vers un précédant billet pour les aspects techniques). Prendre le carré est intéressant, et naturel, car ça fait penser à une distance Euclidienne (usuelle). On va alors poser
Q=\sum_{\mathcal{A},\mathcal{B}} \frac{[n_{\mathcal{A},\mathcal{B}}-n_{\mathcal{A},\mathcal{B}}^\perp]^2}{n_{\mathcal{A},\mathcal{B}}^\perp}
n_{\mathcal{A},\mathcal{B}} est le nombre de personnes observées, de caractéristiques (jointes) \mathcal A et \mathcal B. et n^\perp_{\mathcal{A},\mathcal{B}} est ce qu’on aurait si les variables étaient indépendantes (calculé auparavant), soit

n^\perp_{\mathcal{A},\mathcal{B}}=n\times \frac{n_{\mathcal{A}}}{n}\times \frac{n_{\mathcal{B}}}{n}

avec les notations naturelles. Ici,

> sum((T-Tind)^2/Tind)

[1] 3.380994

Bon… on en fait quoi maintenant de ce 3.38 ? C’est là qu’on utilise le petit résultat évoqué auparavant, qui nous dit que, si n est grand, et que les variables sont effectivement indépendantes, alors Q va suivre une loi du chi-deux (en l’occurrence à 1 degré de liberté, pour tenir compte des diverses contraintes). On va alors se demander s’il est possible, ou vraisemblable, d’avoir 3.38 avec une loi du chi-deux.

> pchisq(3.380994,df=1)
[1] 0.9340477
> 1-pchisq(3.380994,df=1)
[1] 0.06595227

Aussi, une loi du chi-deux à un degré de liberté a 93.4% chances d’être inférieure à la valeur obtenue, ou encore

\mathbb{P}(Q>3.38) \approx 6.59\% avec Q\sim \chi^2(1)

C’est cette grandeur que l’on appelle p-value et que l’on va essayer de comparer au seuil désormais mythique de 5%,

  • si cette probabilité est inférieure à 5%, on va rejeter l’hypothèse d’indépendance
  • si cette probabilité est supérieure à 5%, on va retenir l’hypothèse d’indépendance (pour faire simple)
> chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 3.381, df = 1, p-value = 0.06595

Les plus attentifs devraient me faire remarquer que c’est vaseux mon histoire : tout ce que je raconte est valide si n est grand, et dire que 25 c’est grand… c’est limite ! Tout d’abord, je dirais que c’est ce qui est utilisé dans un paquet d’études en médecine (peut-être un peu plus, mais guère plus). Ensuite je dirais que, malheureusement, la plupart des résultats que l’on voit dans les cours de statistiques sont de ce genre, à savoir des résultats asymptotiques, valides seulement si n est grand (c’est pareil en économie, ou des tests asymptotiques sont parfois évoqués avec 25 années d’observations). Mais heureusement, en s’accrochant un peu, on peut s’affranchir de cette hypothèse. A condition de faire un peu de bootstrap. En fait, l’idée est très très simple. On peut utiliser des générateurs de nombres aléatoires pour générer des tableaux de contingences pour lesquels les variables sont indépendantes. Car générer des variables aléatoires indépendantes, c’est ce qui est le plus simple. En fait, ici, on va se contenter de mélanger les variables,

> s=1
> set.seed(s)
> Xs=sample(X)
> Ys=sample(Y)

En faisant une permutation, je m’assure d’avoir autant de personnes dont la couleur des cheveux est \mathcal A que dans la population initiale (je ne fais que permuter) et de personnes dont la couleur des yeux est \mathcal B. Par contre, comme je mélange indépendamment, je suis certain que les variables sont indépendantes,

> Ts=as.matrix(table(Xs,Ys))
> Ts
Xs      Bleu Noisette
  Blond    7        5
  Brun     7        6

Générer un tableau de contingence, c’est juste pour illustrer, car l’idée est de générer mille (ou plus encore) tableaux de contingence obtenus en simulant des variables indépendantes, et de voir ce que cette statistique du chi-deux donne

> Q=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ Xs=sample(X)
+ Ys=sample(Y)
+ Ts=as.matrix(table(Xs,Ys))
+ Q[s]=sum((Ts-Tind)^2/Tind)}

Si on regarde les valeurs obtenues

>  table(Q) 
0.0509906759906761  0.337162837162837   1.06560106560107   1.92411754911755 
               298                262                200                133 
  3.38099400599401   4.81185481185481    6.9971694971695   9.00037462537462 
                60                 33                  9                  3 
  11.9141275391275    14.489676989677 
                 1                  1

autrement dit, dans 29.8% des scénarios, on avait une distance (un statistique) valant 0.05, dans 26.2% des scénarios, on avait 0.33, et dans 20% des scénarios, on avait 1.06. Etc. Si on cumule on obtient

>  cumsum(table(Q))/10
0.0509906759906761  0.33162837162837   1.06560106560107    1.92411754911755 
3.38099400599401    4.81185481185481    6.9971694971695     9.00037462537462 
11.9141275391275     14.489676989677 

Pour rappel, on avait eu

>  sum((T-Tind)^2/Tind) 
[1] 3.380994

qui peut effectivement être obtenu avec des variables indépendantes. Mais seulement dans 6% des scénarios. En fait, dans 95.3% des scénarios, on a eu une valeur inférieure ou égale à celle obtenue sur notre échantillon. Et dans 4.7% des scénarios, on a strictement dépassé cette valeur. On pourrait parler de p-value obtenue par bootstrap, ou rééchantillonage. On voit que l’on flirte ici avec le seuil des 5%… Visusellement, on a

> plot(c(0,as.numeric(names(table(Q)))),c(0,
+ cumsum(table(Q))/10),type="s",xlab="",ylab="")
> abline(v=sum((T-Tind)^2/Tind),lty=2)
> lines(seq(0,15,by=.01),pchisq(seq(0,15,by=.01),1)*100,col="red")

Ici, on est bien embêté pour conclure quoi que ce soit… Ce qui est amusant, c’est qu’avec un échantillon plus grand (trois fois plus grand par exemple), et les mêmes pourcentages, il n’y a plus aucune ambiguïté,

> T=T*3
> T
X       Bleu Noisette
  Blond   27        9
  Brun    15       24
>  Tind=chisq.test(T)$expected
> Tind Y X Bleu Noisette Blond 20.16 15.84 Brun 21.84 17.16

que l’on peut encore visualiser ci-dessous,

Ici, on obtient une p-value de l’ordre de 0.1%,

>  chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 10.143, df = 1, p-value = 0.001449

On peut aussi faire du bootstrap si on trouve que 75 observations, c’est trop faible, mais la conclusion est la même,

  • Petit exemple, sur la comparaison de moyenne de deux groupes

Cette fois, on a deux populations, et on se demande si les deux groupes sont ‘comparables‘. On peut penser à comparer le nombre de personnes guéries au sein de deux populations, une qui a eu un placebo, et l’autre un médicament. Cela dit, cet exemple simple est en fait compliqué car ‘guérie’ n’est peut-être pas une notion claire. Prenons plutôt des quantités que l’on peut mesurer. Disons la taille des gens, et demandons nous si dans un groupe, les gens sont significativement plus grands que dans l’autre. Je vais tricher un peu… je vais prendre deux populations que je sais être différentes, a priori. Un groupes d’hommes, et un groupes de femmes. On va noter \{x_1,\cdots,x_{n_X}\} et \{y_1,\cdots,y_{n_Y}\} ces deux échantillons, car il va falloir formaliser un minimum.

> Davis=read.table(
+ "")
> Davis[12,c(2,3)]=Davis[12,c(3,2)] 
> Davis=Davis[order(Davis$height),]
> attach(Davis)
> set.seed(1)
> X=sample(height[sex=="F"],size=6)
> Y=sample(height[sex=="M"],size=5)

On a ici 5 personnes dans le premier groupe et 6 personnes dans le second. Et on se demande si les tailles sont significativement différentes entre les deux groupes,

> (mx=mean(X))
[1] 165.5
> (sx=sd(X))
[1] 4.679744
> (my=mean(Y))
[1] 178
> (sy=sd(Y))
[1] 7.968689

On veut savoir si l’écart entre 165,5 cm et 178 cm est significativement différent de zéro (on pourrait aussi se demander s’il est significativement positif, mais ça serait une autre histoire). On va faire une hypothèse ici : on va supposer que nos échantillons sont obtenus comme des tirages de variables aléatoires, indépendantes, x_i=X_i(\omega) et y_j=Y_j(\omega). Le modèle qu’on va supposer est que la taille des gens suit une loi normale, conditionnellement à leur sexe,

\left\{\begin{array}{l}X_i\sim\mathcal{N}(\mu_X,\sigma_X^2) \\Y_j\sim\mathcal{N}(\mu_Y,\sigma_Y^2)\end{array}\right.

Si on veut visualiser nos données, on peu utiliser

> u=seq(150,200,by=.25)
> v1=dnorm(u,mx,sx)
> v2=dnorm(u,my,sy)
> plot(u,v1,col="red",type="l",ylim=c(0,.09),axes=FALSE)
> axis(1)
> lines(u,v2,col="blue")
> points(X,rep(.075,length(X)),col="red")
> points(Y,rep(.075,length(Y)),col="blue")
> abline(v=mx,lty=2,col="red")
> abline(v=my,lty=2,col="blue")


Le plus simple aurait été de supposer que les deux lois ont la même variance, mais quand on regarde les densités sur le graphique ci-dessus, il faut rester crédible. Bref, on va utiliser le test de Welch. L’idée est d’utiliser un petit exercice de probabilité. Soient \{X_1,\cdots,X_{n_X}\} et \{Y_1,\cdots,Y_{n_Y}\} des variables aléatoires, Gaussiennes (comme supposé juste avant de faire un dessin) et indépendantes. Posons alors
\overline{X}=\frac{1}{n_X}\sum_{i=1}^{n_X} X_i
s_X=\frac{1}{n_X-1}\sum_{i=1}^{n_X} [X_i-\overline{X}]^2
(avec les quantités similaires pour le second échantillon) alors si \mu_X=\mu_Y, et si on note
s = \sqrt{{s_X^2 \over n_X} + {s_Y^2 \over n_Y}}.
la variable aléatoire
T = {\overline{X} - \overline{Y} \over s[latex]<br /> va suivre une loi de Student (connue depuis les travaux de <a href="">William Gosset</a>), avec</p> [latex display="true"]\frac{(s_X^2/n_X + s_Y^2/n_Y)^2}{(s_X^2/n_X)^2/(n_X-1) + (s_Y^2/n_Y)^2/(n_Y-1)}

degrés de liberté. On notera qu'il ne s'agit pas ici d'un résultat asymptotique. Par contre, il repose intégralement sur une hypothèse - forte - de normalité des tailles des individus. L'idée est alors la même qu'auparavant : on va calculer cette statistique sur notre échantillon,

> nx=length(X)
> ny=length(Y)
> s=sqrt(sx^2/nx+sy^2/ny)
> T=(mx-my)/s
> T
[1] -3.091371

et se demander s'il est vraisemblable d'avoir une telle valeur avec une loi de Student. Le nombre de degrés de liberté est ici

> d=(sx^2/nx+sy^2/ny)^2/((sx^2/nx)^2/(nx-1)+(sy^2/ny)^2/(ny-1))
> d
[1] 6.218682

On pourrait alors calculer

mais ça serait un test pour savoir si la différence est significativement positive, ce qui n'est pas ce qu'on avait dit que l'on cherchait : on veut ici significativement non-nulle.

\mathbb{P}(\vert T\vert> \vert t\vert)T\sim\mathcal{S}td(d)

Si on veut visualiser un peu tout ça, on peut utiliser

> u=seq(-4,4,by=.01)
> v=dt(u,df=d)
> plot(u,v,type="l")
> abline(v=T,col="red")
> u=seq(-4,T,length=100)
> polygon(c(u,rev(u)),c(dt(u,df=d),rep(0,100)),col="red",border=NA)

Pour calculer ces probabilités, on utilise

> pt(T,d)+(1-pt(-T,d))
[1] 0.02038577

et comme la loi de Student est symétrique par rapport à l’origine,

> pt(T,d)*2           
[1] 0.02038577

Aussi, la p-value est de l'ordre de 2%, ce qui est suffisamment petit pour dire que, oui, nos moyennes sont différentes. C'est d'ailleurs ce que renvoie la commande

> t.test(X,Y)

	Welch Two Sample t-test

data:  X and Y
t = -3.0914, df = 6.219, p-value = 0.02039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.31039  -2.68961
sample estimates:
mean of x mean of y 
    165.5     178.0

On avait vu tout à l'heure que l'idée de la valeur de 5% était liée au scénario qui se répète avec 1 chance sur 20 : si on pouvait générer d'autres échantillons, indépendants, alors 1 fois sur 20 on devrait atteindre une statistique de test de l'ordre de celle observée. Pour générer des échantillons vraiment indépendants, on avait utilisé du bootstrap, en rééchantillonnant dans notre population de test. Ici, on a bien plus que nos 11 observations, on en a en réalité 200. On peut se demande ce qui se serait passé si on avait constitué d'autres échantillons,

> Diff=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ X=sample(height[sex=="F"],size=6)
+ Y=sample(height[sex=="M"],size=5) 
+ Diff[s]=(mean(X)-mean(Y))/(sqrt(sd(X)^2/nx+sd(Y)^2/ny))}
> hist(Diff,col="light blue",border="white")

Sur aucun échantillon, la moyenne du premier groupe n'a dépassé la moyenne du second. Alors qu'on aurait pu penser en avoir une vingtaine (2% des mille scénarios générés). Ceci vient peut-être du fait que l'hypothèse de normalité pour les tailles peut être trop restrictive...

  • Pour aller plus loin

Le billet aura été un peu long, mais on aura eu le temps de voir en détails d'où vient cette notion de significativité statistique, et ses liens avec la p-value. On aura même eu le temps de voir sur deux exemples l'importance des hypothèses sous-jacentes, en particulier si le test est asymptotique, ou s'il repose sur des hypothèses fortes de lois. Mais ce n'était qu'une courte mise en bouche. Pour aller plus loin, je suggère la lecture de Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle, M.V. (2008) Sound and fury: McCloskey and significance testing in economics, Liberman, M. (2013). "Significance", in 1885 and today et le passionnant Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Et sur les p-value, Goodman, S. (2002) A Dirty Dozen:Twelve P-Value Misconceptions. Je devrais aussi citer Ziliak, S. & McCloskey, D. (2008). The cult of statistical significance, mais je ne l'ai pas encore lu. Par contre, Xi'an en parlait sur son blog. Et si j'étais courageux, je tenterais une lecture bayésienne de ce billet.... promis, dès que je trouve un peu de temps !

Somewhere else, part 85

Some writings worth reading, starting with

Don’t get me wrong, however, I don’t want to overly fetishize the task of coding, but in part, knowing to code makes you better able to think about and work with digital technology. Being able to code is, well, code for having a better facility with technology. And that, that my friends, is ultimately about being an effective citizen in our data driven, code-ridden, algorithm rhythm world. My answer to the Atlantic writer’s question of course is very different. I think all journalists need to code because an effective democracy depends on all citizens being able to understand and think about technology. This is not done from the sidelines. You have to get your hands dirty. You have to write some code. [to be continued…]

These studies confirm the common belief that good teachers can get bad evaluations: Teaching effectiveness, as measured by subsequent performance and career success, is negatively associated with student teaching evaluations. While one should be cautious in generalizing the conclusions because the two student populations might not be representative of students at large (or at least of Berkeley students), these are by far the best studies we know of. They are the only controlled, randomized experiments; they are from different continents and cultures; and their findings are concordant. [to be continued…]

et un peu de lecture en français,

Did I miss something?

Pricing Reinsurance Contracts

In order to illustrate the next section of the non-life insurance course, consider the following example1, inspired from…. This is the so-called “Normalized Hurricane Damages in the United States” dataset, for the period 1900-2005, from Pielke et al. (2008). The dataset is available in xls format, so we have to spend some time to import it,

> library(gdata)
> db=read.xls(
+ "",
+ sheet=1)
trying URL ''

Content type 'application/' length 119296 bytes (116 Kb)
opened URL
downloaded 116 Kb

perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
	LANGUAGE = "fr_CA:fr",
	LC_ALL = (unset),
	LANG = "fr_CA.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

The problem with excel spreadsheets is that some columns might have pre-specified format (here, losses are with a format 000,000,000 for instance)

> tail(db)
    Year Hurricane.Description State Category Base.Economic.Damage
202 2005                 Cindy    LA        1          320,000,000
203 2005                Dennis    FL        3        2,230,000,000
204 2005               Katrina LA,MS        3       81,000,000,000
205 2005               Ophelia    NC        1        1,600,000,000
206 2005                  Rita    TX        3       10,000,000,000
207 2005                 Wilma    FL        3       20,600,000,000
    Normalized.PL05 Normalized.CL05  X X.1
202     320,000,000     320,000,000 NA  NA
203   2,230,000,000   2,230,000,000 NA  NA
204  81,000,000,000  81,000,000,000 NA  NA
205   1,600,000,000   1,600,000,000 NA  NA
206  10,000,000,000  10,000,000,000 NA  NA
207  20,600,000,000  20,600,000,000 NA  NA

To get data in a format we can play with, consider the following function,

> stupidcomma = function(x){
+ x=as.character(x)
+ for(i in 1:10){x=sub(",","",as.character(x))}
+ return(as.numeric(x))}

and let’s convert those values into numbers,

> base=db[,1:4]
> base$Base.Economic.Damage=Vectorize(stupidcomma)(db$Base.Economic.Damage)
> base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05)
> base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05)

Here is the dataset we will use, from now on,

> tail(base)
    Year Hurricane.Description State Category Base.Economic.Damage
202 2005                 Cindy    LA        1             3.20e+08
203 2005                Dennis    FL        3             2.23e+09
204 2005               Katrina LA,MS        3             8.10e+10
205 2005               Ophelia    NC        1             1.60e+09
206 2005                  Rita    TX        3             1.00e+10
207 2005                 Wilma    FL        3             2.06e+10
    Normalized.PL05 Normalized.CL05
202        3.20e+08        3.20e+08
203        2.23e+09        2.23e+09
204        8.10e+10        8.10e+10
205        1.60e+09        1.60e+09
206        1.00e+10        1.00e+10
207        2.06e+10        2.06e+10

We can visualize the normalized costs of hurricanes, from 1900 till 2005, with the 207 hurricanes (here the x-axis is not time, it is simply the index of the loss)

> plot(base$Normalized.PL05/1e9,type="h",ylim=c(0,155))

As usual, there are two components when computing the pure premium of an insurance contract. The number of claims (or here hurricanes) and the individual losses of each claim. We’ve seen – above – individual losses, let us focus now on the annual frequency.

> TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+ counts=c(counts,rep(0,length(years0))))
> db[88:93,]
   years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0

On average, we experience about 2 (major) hurricanes per year,

> mean(db$counts)
[1] 1.95283

In predictive modeling (here, we wish to price a reinsurance contract for, say, 2014), we need probably to take into account some possible trend in the hurricane occurrence frequency. We can consider either a linear trend,

> reg0 <- glm(counts~years,data=db,family=poisson(link="identity"),
+ start=lm(counts~years,data=db)$coefficients)

or an exponential one,

> reg1 <- glm(counts~years,data=db,family=poisson(link="log"))

We can plot those three predictions, and get a prediction for the number of (major) hurricanes in 2014,

> plot(years,counts,type='h',ylim=c(0,6),xlim=c(1900,2020))
> cpred1=predict(reg1,newdata=data.frame(years=1890:2030),type="response")
> lines(1890:2030,cpred1,col="blue")
> cpred0=predict(reg0,newdata=data.frame(years=1890:2030),type="response")
> lines(1890:2030,cpred0,col="red")
> abline(h=mean(db$counts),col="black")
> (predictions=cbind(constant=mean(db$counts),linear=
+ cpred0[126],exponential=cpred1[126]))
    constant   linear exponential
126  1.95283 3.573999    4.379822
> points(rep((1890:2030)[126],3),prediction,col=c("black","red","blue"),pch=19)

Observe that changing the model will change the pure premium: with a flat prediction, we expect less than 2 (major) hurricanes, but with the exponential trend, we expect more than 4…

This is for the expected frequency. Now, we should find a suitable model to compute the pure premium of a reinsurance treaty, with a (high) deductible, and a limited (but large) cover. As we will seen in class next week, the appropriate model is a Pareto distribution (see Hagstrœm (1925), Huyghues-Beaufond (1991) or a survey – in French – published a few years ago).

We can use Hill’s plot to estimate the tail index,

> library(evir)
> hill(base$Normalized.PL05)

Clearly, costs of major hurricanes are heavy tailed.

Now, consider an insurance company, in the U.S., with 5% market share (just to illustrate). We will consider there \tilde Y_i= Y_i/20. The losses are given below. Consider a reinsurance treaty, with a deductible of 2 (billion) and a limited cover of 4 (billion),

For our Pareto model, consider only losses above 500 millions,

> threshold=.5
> (gpd.PL <- gpd(base$Normalized.PL05/1e9/20,threshold)$par.ests)
       xi      beta 
0.4424669 0.6705315

Keep in mind the 1 hurricane out of 8 reaches that level

> mean(base$Normalized.CL05/1e9/20>.5)
[1] 0.1256039

Given that the loss exceeds 500 millions, we can now compute the expected value of the reinsurance contact,

To compute it we can use

> E <- function(yinf,ysup,xi,beta){
+   as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta),
+   lower=yinf,upper=ysup)$value+
+   (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf))
+ }

[Nov 5th] there is a typo in the previous function, since the threshold should be used, here, as a parameter in the function, if you want to play with that function an see the impact of the threshold (see a more recent post on the same topic, but a different dataset)… but here, we do not change the threshold, so it is not a big deal.

Now, it is probably time to bring all the pieces together. We might expect a bit less than 2 (major) hurricanes per year,

> predictions[1]
[1] 1.95283

and each hurricane has 12.5% chances to cost more than 500 million for our insurance company,

> mean(base$Normalized.PL05/1e9/20>.5)
[1] 0.1256039

and given that an hurricane exceeds 500 million loss, then the expected repayment by the reinsurance company is (in millions)

> E(2,6,gpd.PL[1],gpd.PL[2])*1e3
[1] 330.9865

So the pure premium of the reinsurance contract is simply

> predictions[1]*mean(base$Normalized.PL05/1e9/20>.5)*
+ E(2,6,gpd.PL[1],gpd.PL[2])*1e3
[1] 81.18538

for a cover of 4 billion, in excess of 2.

1.This example will be found in the Reinsurance and Extremal Events chapter in the forthcoming Computational Actuarial Science with R, by Eric Gilleland and Mathieu Ribatet.

Somewhere else, part 84

Some writings worth reading

F.P. Ramsey has some claim to be the greatest philosopher of the twentieth century. In Cambridge in the 1920s, he single handedly forged a range of ideas that have since come to define the philosophical landscape. Contemporary debates about truth, meaning, knowledge, logic and the structure of scientific theories all take off from positions first defined by Ramsey. Equally importantly, he figured out the principles governing subjective probability, and so opened the way to decision theory, game theory and much work in the foundations of economics. His fertile mind could not help bubbling over into other subjects. An incidental theorem he proved in a logic paper initiated the branch of mathematics known as Ramsey theory, while two articles in the Economic Journal pioneered the mathematical analysis of taxation and saving. [to be continued…]

Students who acquire large debts putting themselves through school are unlikely to think about changing society, Chomsky suggested. “When you trap people in a system of debt . they can’t afford the time to think.” Tuition fee increases are a “disciplinary technique,” and, by the time students graduate, they are not only loaded with debt, but have also internalized the “disciplinarian culture.” This makes them efficient components of the consumer economy.

et un tout petit peu de lecture en français ces dernieres jours

Comme la radio est devenue, pour des raisons pratiques, le média privilégié de la circulation automobile, les journaux papier sont de plus en plus des objets de consommation ponctuelle, dans des situations de déconnexion, particulièrement les transports en commun. Alors que les kiosques à journaux périclitent, les gares ou les aéroports comptent parmi les derniers endroits où le commerce de l’information reste vivace. Pourtant, l’autre jour, en rentrant de voyage, je suis ressorti les mains vides de la librairie, malgré la perspective d’un long trajet en RER. Ce n’est pas la première fois que la ribambelle des Unes échoue à éveiller mon désir. Si cette offre ne me tente pas, c’est parce que mes propres outils de sélection des sources m’éloignent des récits médiatiques les plus courants, qui perdent de leur pertinence à mes yeux. J’ai donc passé mon trajet à lire et à commenter mes flux Facebook et Twitter. Une activité moins confortable que la lecture d’un magazine, compte tenu de l’étroitesse de l’écran de mon smartphone et d’une connexion 3G parfois fluctuante, mais néanmoins plus satisfaisante que la consommation d’un support d’information non interactif. La raison de cette désaffection n’est pas évidente pour tout le monde. Selon Bernard Guetta, prix Albert-Londres et éditorialiste vedette à France-Inter, «la crise de la presse occidentale est avant tout celle des grands courants de pensée européens et américains (…), la conséquence de la panne d’idées occidentale». Si l’expert en géopolitique lisait un peu moins ses organes favoris, et un peu plus internet, il serait étonné de la diversité et de la richesse de pensée qui s’y exprime, très loin de la “panne d’idées”. On pourrait plus opportunément pointer la frilosité et le conformisme du filtre médiatique, que ses contraintes économiques poussent vers un lectorat de plus en plus âgé et nanti. Construits par le jeu des affinités et par tests d’essai/erreur, les bouquets informationnels des réseaux sociaux proposent par définition un ciblage plus fin et plus adapté que n’importe quel média de masse. Mais ce n’est pas la seule raison qui les impose comme l’alternative définitive à la consommation de la presse. [a suivre…]

Did I miss something?

GLM, non-linearity and heteroscedasticity

Last week in the non-life insurance course, we’ve seen the theory of the Generalized Linear Models, emphasizing the two important components

  • the link function (which is actually the key component in predictive modeling)
  • the distribution, or the variance function

Just to illustrate, consider my favorite dataset

­lin.mod = lm(dist~speed,data=cars)

A linear model means here Y_i=\beta_0+\beta_1X_i+\varepsilon_i

where the residuals are assumed to be centered, independent, and with identical variance. If we visualize that linear regression, we usually see something like that

The idea here (in GLMs) is to assume Y\vertX=x\sim\mathcal{N}(\beta_0+\beta_1x,\sigma^2)

which will produce the same model as the one describe previously, based on some error term. That model can be visualized below,

X= cars$speed 
mat=persp(vX,vY,matrix(0,n,n),zlim=c(0,.1),theta=-30,ticktype ="detailed", box = FALSE)
y1=qnorm(.95,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
y2=qnorm(.05,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
for(j in n:1){
y=seq(min(min(Y)-15,qnorm(.05,predict(reggig,newdata=data.frame(X=vX[j]),type="response"), sdgig)),max(Y)+15,length=stp)
z=dnorm(y, mgig[j], sdgig)
polygon(C,border=NA,col="light blue",density=40)

We do have two parts here: the linear increase of the average, \mathbb{E}(Y\vert X=x)=\beta_0+\beta_1x and the constant variance of the normal distribution \text{Var}(Y\vert X=x)=\sigma^2.

On the other hand, if we assume a Poisson regression,

poisson.reg = glm(dist~speed,data=cars,family=poisson(link="log"))

we have something like

This time, two things have changed simultaneously: our model is no longer linear, it is an exponential one \mathbb{E}(Y\vert X=x)=e^{\beta_0+\beta_1x}, and the variance is also increasing with the explanatory variable \text{Var}(Y\vert X=x)=e^{\beta_0+\beta_1x}, since with a Poisson regression,
Y\vert X=x\sim\mathcal{P}(e^{\beta_0+\beta_1x})

If we adapt the previous code, we get

The problem is that we changed two things when we introduced the Poisson regression from the linear model. So let us look at what happens when we change the two components independently. First, we can change the link function, with a Gaussian model but this time a multiplicative model (with a logarithm link function)

gaussian.reg = glm(dist~speed,data=cars,family=gaussian(link="log"))

which is still, here, an homoscedasctic model, but this time non-linear. Or we can change the link function in the Poisson regression, to get a linear model, but heteroscedastic

poisson.lin = glm(dist~speed,data=cars,family=poisson(link="identity"))

So this is basically what GLMs are about….

Modélisation des coûts individuels

Cette semaine, même si le réseau de l’UQAM est down, on va continuer le cours et finir la section sur la modélisation de la surdispersion pour la fréquence de sinistres. On devrait ensuite commencer la modélisation des coûts individuels. En particulier, on passera du temps autour de deux points,

  • la distinction lognormale et gamma
  • l’écrêtement des gros sinistres

Les transparents sont en ligne. Et la base des coûts est celle évoquée au second cours.

Somewhere else, part 83

Some writings worth reading, starting with

but also


avec un peu de lecture en français

Did I miss something ?