# Fairness and discrimination, PhD Course, #3 Machine Learning, losses and distances

For the third course, we will get back a little bit on machine learning (slides are still online on the github repository). The starting point will be loss functions and risk.

#### Loss functions and risk

A general definition for a loss is that it is positive, and null when we consider $\ell(y,y)$. As we will discuss further, it is neither a distance, nor a dissimilarity measure

Then, define the empirical risk (and the associated empirical risk minimization principle, as coined in Vapnik (1991))

Given a loss $\ell$ and some probabilistic space, define the optimal decision, also called Bayes decision rule

And instead of the risk of a model, define the excess of risk.

A classical loss for a classifier is $\ell_{0/1}$,

In that case, Bayes decision rule, is$$m^\star(\boldsymbol{x}) = \boldsymbol{1}(\mu(\boldsymbol{x})>1/2) =\begin{cases}1 \text{ if }\mu(\boldsymbol{x})>1/2\\0 \text{ if }\mu(\boldsymbol{x})\leq1/2\\\end{cases}$$where (of course), one needs to know $\mu$, otherwise, we can consider some plug-in estimator based on $\widehat\mu$. For continuous variables $y$, consider the quadratic loss $\ell_2$,

In that case, Bayes decision rule (the optimal model) is the conditional expectation

Observe that we can also define the quantile loss (or the expectile loss)

Observe that this loss is not symmetric..

#### From loss functions to distances

Let us discuss a bit more the fact that losses are not distances. As mentioned, it is neither necessarily symmetric nor seperable,

But furthermore, it has no reason to satisfy the triangle inequality. Actually, if $d$ is the distance, it is very likely that $d^2$ is not (since exponentiating is not a subadditive transformation)

Another related concept could be the concept of similarity, or dissimilarity.

Another one is the concept of divergence, that we will use much more. For instance, Bregman divergence is

which safisfies desirables properties.

Interestingly it is possible to define “projections” even if we have neither an orthogonal projection (since there is no orthogonal concept since there is no inner-product), nor a distance. But still

One can use a nice algorithm to estimate that quantity, if the convex set can we expressed simply

When considering “distances” between distributions, instead of $y$‘s, among other interesting properties in statistics, we can mention the one of unbiased gradients,

and Müller (1997) defined integral probability metrics

#### Standard “distances” between distributions

The first one will be Hellinger distance

that can lead so simple expressions for standard parametric distributions, such as Beta distributions,

or (multivariate) Gaussian ones

We can also mention Pearson divergence

More interesting (and popular in probability theory), total variation

There are several ways to express that distance.

If instead of general sets $\mathcal{A}$, we can consider half lines, (-\infty,t][\latex], and we obtain Kolmogorov distance (or Kolmogorov-Smirnov)

Another important one in statistics is Kullback–Leibler divergence

For instance, with Gaussian vectors

Observe that the measure is actually a dissimilarity measure

If we want a symmetric version, we can consider Jeffreys divergence

Finally, we will mention f-divergence

We will discuss a little bit more those "distances" (yes, I usually use that term, abusively) and next week, we will present the most interesting distance, that will be Wasserstein's.

# Puissance et tests statistiques (simples)

Pour comprendre la notion de puissance dans un test statistique, reprenons Statistical power analysis for the behavioral sciences de Jacob Cohen. Comme il le rappelle Fisher (dans The Design of Experiments) ne souhaite pas “prouver” qu’une hypothèse (appelée H_0) est valide, mais il espère, en réalité, rejeter l’hypothèse en question (c’est ce que rappellent d’ailleurs Hubbard et Bayarri en 2003). Autrement dit, supposons que l’on cherche à faire un test médical, et à détecter si une personne est malade (ou pas). Si une personne est généralement malade dès lors qu’un taux de quelque chose est faible (noté $\theta$), on va chercher à tester $H_0:\theta_\star>\theta_0$. Dans un test statistique, on espère rejeter l’hypothèse nulle $H_0$.

The power of a statistical test of a null hypothesis is the probability that it will lead to the rejection of the null hypothesis.

# Independence and correlation

A short post to get back on a property I gave briefly in the MAT8595 class in January, and again in the MAT8181 class this week (to illustrate the distinction between weak and strong white noises). Recall that (real-valued) random variables $X$ and $Y$ are independent if for all $A, B\subset\mathbb{R}$, $\mathbb{P}((X,Y)\in A\times B)=\mathbb{P}(X\in A)\cdot \mathbb{P}(Y\in B)$ Another characterization, for integrable variable is that $\mathbb{E}(f(X)\cdot g(Y)) = \mathbb{E}(f(X))\cdot \mathbb{E}( g(Y))$ for all $f,g:\mathbb{R}\mapsto \mathbb{R}$, which can be written, if variables are square integrable $\text{corr}(f(X), g(Y))=0$ The idea to prove this characterization is to observe that if $X$ and $Y$ are independent $\mathbb{P}((X,Y)\in A\times B)=\mathbb{P}(X\in A)\cdot \mathbb{P}(Y\in B)$ can be written $\mathbb{E}(\boldsymbol{1}_A(X)\cdot \boldsymbol{1}_B(Y)) = \mathbb{E}(\boldsymbol{1}_A(X))\cdot \mathbb{E}( \boldsymbol{1}_B(Y))$ Using a standard argument in integration theory, equality $\mathbb{E}(f(X)\cdot g(Y)) = \mathbb{E}(f(X))\cdot \mathbb{E}( g(Y))$ is valid for step functions (not only indicators), and then to positive measurable functions, and finally to integrable functions. Proving this result is not that difficult. Observe that Rényi (1959) – inspired by Gebelein (1947) – followed by Sarmanov (1958) introduced the concept of maximal correlation, that can be related to this result, $\text{corr}^{\star}(X,Y)=\max\{\text{corr}(f(X),g(Y))\}$ where the maximum is taken over all functions $f$ and $g$ such that the correlation exist. Actually, it is possible to consider only transformations such that $\mathbb{E}(f(X))=0$ and $\text{Var}(f(X))=1$ (and similarly for $g$, the idea is that we simple center and scale, which does not impact the correlation.Thus, $X$ and $Y$ are independent if and only if $\text{corr}^{\star}(X,Y)=0$ Algorithm to estimate that coefficient are interesting. The problem can be written, equivalently $(f^\star,g^\star)=\text{argmin}\{\mathbb{E}[f(X)-g(Y)]^2\}$ And if the minimization is considered over $g$, assuming that $f$ is fixed, then the optimal transformation is $g(y)=\mathbb{E}[f(X)\vert Y=y]$ And similarly for $f$. So using an iterative algorithm, it is possible to get $f^\star$ and $g^\star$ (see Breiman and Friedman (1985) for more details). Actually, those functions appear in nonlinear canonical analysis. As mentioned in Lancaster (1957), for a Gaussian random vector $(X,Y)$ $\text{corr}^{\star}(X,Y)=\vert \text{corr}(X,Y)\vert$ and in that case  $f^\star$ and $g^\star$ are affine functions. This can be related to Hermite’s polynomial and to the expansion of the bivariate Gaussian density. I still hope that someone will go further for the project in the MAT8181 course.

# 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
Y
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 Y 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 Ys 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) 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 Y 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(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> 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$$ et $$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] va suivre une loi de Student (connue depuis les travaux de William Gosset), avec [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 $\mathbb{P}(T>t)$$T\sim\mathcal{S}td(d)$ 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 ! # Overdispersed Poisson et bootstrap Pour le dernier cours sur les méthodes de provisionnement, on s’est arrête aux méthodes par simulation. Reprenons là où on en était resté au dernier billet où on avait vu qu’en faisant une régression de Poisson sur les incréments, on obtenait exactement le même montant que la méthode Chain Ladder, > Y [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA > y=as.vector(as.matrix(Y)) > base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) > reg2=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson) > summary(reg2) Call: glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson, data = base) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 8.05697 0.01551 519.426 < 2e-16 *** as.factor(ai)2001 0.06440 0.02090 3.081 0.00206 ** as.factor(ai)2002 0.20242 0.02025 9.995 < 2e-16 *** as.factor(ai)2003 0.31175 0.01980 15.744 < 2e-16 *** as.factor(ai)2004 0.44407 0.01933 22.971 < 2e-16 *** as.factor(ai)2005 0.50271 0.02079 24.179 < 2e-16 *** as.factor(bj)1 -0.96513 0.01359 -70.994 < 2e-16 *** as.factor(bj)2 -4.14853 0.06613 -62.729 < 2e-16 *** as.factor(bj)3 -5.10499 0.12632 -40.413 < 2e-16 *** as.factor(bj)4 -5.94962 0.24279 -24.505 < 2e-16 *** as.factor(bj)5 -5.01244 0.21877 -22.912 < 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: 46695.269 on 20 degrees of freedom Residual deviance: 30.214 on 10 degrees of freedom (15 observations deleted due to missingness) AIC: 209.52 Number of Fisher Scoring iterations: 4 > base$py2=predict(reg2,newdata=base,type="response")
> round(matrix(base$py2,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.7 1202.1 49.8 19.1 8.2 21.0 [2,] 3365.6 1282.1 53.1 20.4 8.8 22.4 [3,] 3863.7 1471.8 61.0 23.4 10.1 25.7 [4,] 4310.1 1641.9 68.0 26.1 11.2 28.7 [5,] 4919.9 1874.1 77.7 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.4 31.6 13.6 34.7 > sum(base$py2[is.na(base$y)]) [1] 2426.985 Le plus intéressant est peut être de noter que la loi de Poisson présente probablement trop peu de variance… > reg2b=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson) > summary(reg2) Call: glm(formula = y ~ as.factor(ai) + as.factor(bj), family = quasipoisson, data = base) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.05697 0.02769 290.995 < 2e-16 *** as.factor(ai)2001 0.06440 0.03731 1.726 0.115054 as.factor(ai)2002 0.20242 0.03615 5.599 0.000228 *** as.factor(ai)2003 0.31175 0.03535 8.820 4.96e-06 *** as.factor(ai)2004 0.44407 0.03451 12.869 1.51e-07 *** as.factor(ai)2005 0.50271 0.03711 13.546 9.28e-08 *** as.factor(bj)1 -0.96513 0.02427 -39.772 2.41e-12 *** as.factor(bj)2 -4.14853 0.11805 -35.142 8.26e-12 *** as.factor(bj)3 -5.10499 0.22548 -22.641 6.36e-10 *** as.factor(bj)4 -5.94962 0.43338 -13.728 8.17e-08 *** as.factor(bj)5 -5.01244 0.39050 -12.836 1.55e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 3.18623) Null deviance: 46695.269 on 20 degrees of freedom Residual deviance: 30.214 on 10 degrees of freedom (15 observations deleted due to missingness) AIC: NA Number of Fisher Scoring iterations: 4 Mais on en reparlera dans un instant. Ensuite, on avait commencé à regarder erreurs commises, sur la partie supérieure du triangle. Classiquement, par construction, les résidus de Pearson sont de la forme $\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\text{Var}(Y_i)}}$ On avait vu dans le cours de tarification que la variance au dénominateur pouvait être remplacé par le prévision, puisque dans un modèle de Poisson, l’espérance et la variance sont identiques. Donc on considérait $\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$ > base$erreur=(base$y-base$py2)/sqrt(base$py2) > round(matrix(base$erreur,n,n),1)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]  0.9 -1.1 -1.5 -0.5 -0.4    0
[2,]  0.0  0.3 -2.2  0.8  0.4   NA
[3,]  0.1  0.1 -1.0 -0.3   NA   NA
[4,] -1.1  0.9  4.2   NA   NA   NA
[5,]  0.1 -0.2   NA   NA   NA   NA
[6,]  0.0   NA   NA   NA   NA   NA

Le soucis est que si $\widehat{Y}_i$ est – asymptotiquement – un bon estimateur pour $\text{Var}(Y_i)$, ce n’est pas le cas à distance finie, car on a alors un estimateur biaisé pour la variance, et donc la variance des résidus n’a que peu de chances d’être de variance unitaire. Aussi, il convient de corriger l’estimateur de la variance, et on pose alors

$\varepsilon_i=\sqrt{\frac{n}{n-k}}\cdot\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$

qui sont alors les résidus de Pearson tel qu’on doit les utiliser.

> E=base$erreur[is.na(base$y)==FALSE]*sqrt(21/(21-11))
> E
[1]  1.374976e+00  3.485024e-02  1.693203e-01 -1.569329e+00  1.887862e-01
[6] -1.459787e-13 -1.634646e+00  4.018940e-01  8.216186e-02  1.292578e+00
[11] -3.058764e-01 -2.221573e+00 -3.207593e+00 -1.484151e+00  6.140566e+00
[16] -7.100321e-01  1.149049e+00 -4.307387e-01 -6.196386e-01  6.000048e-01
[21] -8.987734e-15
> boxplot(E,horizontal=TRUE)

En rééchantillonnant dans ces résidus, on peut alors générer un pseudo triangle. Pour des raisons de simplicités, on va générer un peu rectangle, et se restreindre à la partie supérieure,

> Eb=sample(E,size=36,replace=TRUE)
> Yb=base$py2+Eb*sqrt(base$py2)
> Ybna=Yb
> Ybna[is.na(base$y)]=NA > Tb=matrix(Ybna,n,n) > round(matrix(Tb,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3115.8 1145.4 58.9 46.0 6.4 26.9 [2,] 3179.5 1323.2 54.5 21.3 12.2 NA [3,] 4245.4 1448.1 61.0 7.9 NA NA [4,] 4312.4 1581.7 68.7 NA NA NA [5,] 4948.1 1923.9 NA NA NA NA [6,] 4985.3 NA NA NA NA NA Cette fois, on a un nouveau triangle ! on va alors pouvoir faire plusieurs choses, 1. compléter le triangle para la méthode Chain Ladder, c’est à dire calculer les montants moyens que l’on pense payer dans les années futures 2. générer des scénarios de paiements pour les années futurs, en générant des paiements suivant des lois de Poisson (centrées sur les montants moyens que l’on vient de calculer) 3. générer des scénarios de paiements avec des lois présentant plus de variance que la loi de Poisson. Idéalement, on voudrait simuler des lois qusi-Poisson, mais ce ne sont pas de vraies lois. Par contre, on peut se rappeler que dans ce cas, la loi Gamma devrait donner une bonne approximation. Pour ce dernier point, on va utiliser le code suivant, pour générer des quasi lois, > rqpois = function(n, lambda, phi, roundvalue = TRUE) { + b = phi + a = lambda/phi + r = rgamma(n, shape = a, scale = b) + if(roundvalue){r=round(r)} + return(r) + } Je renvois aux diverses notes de cours pour plus de détails sur la justification, ou à un vieux billet. On va alors faire une petite fonction, qui va soit somme les paiements moyens futurs, soit sommer des générations de scénarios de paiements, à partir d’un triangle, > CL=function(Tri){ + y=as.vector(as.matrix(Tri)) + base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) + reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson) + py2=predict(reg,newdata=base,type="response") + pys=rpois(36,py2) + pysq=rqpois(36,py2,phi=3.18623) + return(list( + cl=sum(py2[is.na(base$y)]),
+ sc=sum(pys[is.na(base$y)]), + scq=sum(pysq[is.na(base$y)])))
+ }

Reste alors à générer des paquets de triangles. Toutefois, il est possible de générer des triangles avec des incréments négatifs. Pour faire simple, on mettra des valeurs nulles quand on a un paiement négatif. L’impact sur les quantiles sera alors (a priori) négligeable.

> for(s in 1:1000){
+ Eb=sample(E,size=36,replace=TRUE)*sqrt(21/(21-11))
+ Yb=base$py2+Eb*sqrt(base$py2)
+ Yb=pmax(Yb,0)
+ scY=rpois(36,Yb)
+ Ybna=Yb
+ Ybna[is.na(base$y)]=NA + Tb=matrix(Ybna,6,6) + C=CL(Tb) + VCL[s]=C$cl
+ VR[s]=C$sc + VRq[s]=C$scq
+ }

Si on regarde la distribution du best estimate, on obtient

> hist(VCL,proba=TRUE,col="light blue",border="white",ylim=c(0,0.003))
> D=density(VCL)
> lines(D)
> I=which(D$x<=quantile(VCL,.05)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
> I=which(D$x>=quantile(VCL,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

Mais on peut aussi visualiser des scénarios basés sur des lois de Poisson (équidispersé) ou des scénarios de lois quasiPoisson (surdispersées), ci-dessous

Dans ce dernier cas, on peut en déduire le quantile à 99% des paiements à venir.

> quantile(VRq,.99)
99%
2855.01

Il faut donc augmenter le montant de provisions de l’ordre 15% pour s’assurer que la compagnie pourra satisfaire ses engagements dans 99% des cas,

> quantile(VRq,.99)-2426.985
99%
428.025

# Quels résidus et quelle loi simuler ?

Suite à une question sur les résidus, je vais essayer de prendre deux minutes pour reprendre un point que je n’avais pas abordé en cours, faute de temps. Les résidus les plus intéressants en provisionnement (et dans les GLM) sont les résidus de Pearson,

Mais il existe également des résidus dits ajustés, permettant de prendre en compte le fait que le nombre de paramètres est ici relativement grand, ramené au nombre d’observations. On pose alors

A première vue, utiliser l’un ou l’autre devrait donner la même chose si l’on génère des pseudo-triangles, puisque dans un cas, on utiliserait

et dans l’autre

Bref, multiplier par une constante pour ensuite diviser par la même constante, ça revient au même. Sauf qu’il peut sembler légitime de poser

(on reviendra là dessus tout à l’heure). En attendant, voilà nos deux résidus

> (E=residuals(regp,"pearson"))
1             2             3             4
9.488238e-01  2.404895e-02  1.168421e-01 -1.082940e+00
5             6             7             8
1.302749e-01 -1.007348e-13 -1.128013e+00  2.773332e-01
9            10            11            13
5.669707e-02  8.919633e-01 -2.110748e-01 -1.533031e+00
14            15            16            19
-2.213449e+00 -1.024162e+00  4.237393e+00 -4.899687e-01
20            21            25            26
7.929194e-01 -2.972380e-01 -4.275912e-01  4.140426e-01
31
-6.202125e-15
> n=sum(is.na(Y)==FALSE)
> k=ncol(PAID)+nrow(PAID)-1
> (R=residuals(regp,"pearson")*sqrt(n/(n-k)))
1             2             3             4
1.374976e+00  3.485024e-02  1.693203e-01 -1.569329e+00
5             6             7             8
1.887862e-01 -1.459787e-13 -1.634646e+00  4.018940e-01
9            10            11            13
8.216186e-02  1.292578e+00 -3.058764e-01 -2.221573e+00
14            15            16            19
-3.207593e+00 -1.484151e+00  6.140566e+00 -7.100321e-01
20            21            25            26
1.149049e+00 -4.307387e-01 -6.196386e-01  6.000048e-01
31
-8.987734e-15

L’autre point sur lequel je voulais revenir sur les simulations, est que la loi de Poisson n’est peut-être pas adaptée pour simuler des scénarios de paiements futurs. En effet, si on fait une régression quasipoisson, on voit que le paramètre de surdispersion n’est pas négligeable,

> regqp=glm(Y~as.factor(D)+as.factor(A),
> summary(regqp)

Call:
glm(formula = Y ~ as.factor(D) + as.factor(A),
data = base)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3426  -0.4996   0.0000   0.2770   3.9355

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)       8.05697    0.02769 290.995  < 2e-16 ***
as.factor(D)2    -0.96513    0.02427 -39.772 2.41e-12 ***
as.factor(D)3    -4.14853    0.11805 -35.142 8.26e-12 ***
as.factor(D)4    -5.10499    0.22548 -22.641 6.36e-10 ***
as.factor(D)5    -5.94962    0.43338 -13.728 8.17e-08 ***
as.factor(D)6    -5.01244    0.39050 -12.836 1.55e-07 ***
as.factor(A)2002  0.06440    0.03731   1.726 0.115054
as.factor(A)2003  0.20242    0.03615   5.599 0.000228 ***
as.factor(A)2004  0.31175    0.03535   8.820 4.96e-06 ***
as.factor(A)2005  0.44407    0.03451  12.869 1.51e-07 ***
as.factor(A)2006  0.50271    0.03711  13.546 9.28e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Dispersion parameter for quasipoisson family taken to be 3.18623

Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

On peut alors envisager de simuler non pas des lois de Poisson, mais des lois quasipoissons, en reprenant un vieux code,

> rqpois = function(n, lambda, phi, roundvalue = TRUE) {
+ b = phi
+  a = lambda/phi
+ r = rgamma(n, shape = a, scale = b)
+  if(roundvalue){r=round(r)}
+ return(r)
+ }

A partir de ces deux remarques, on peut reprendre le code qui permettait de générer des scénarios de paiements pour les années futures,

> Yp=predict(regp,type="response",newdata=base)
> Rs=rep(NA,20000)
> for(s in 1:20000){
+ serreur=sample(erreur,
+ size=36,replace=TRUE)
+ E=matrix(serreur,6,6)
+ sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6))
+ sbase=data.frame(sY=as.vector(sY),D,A)
+ sbase$sY[is.na(Y)==TRUE]=NA + sreg=glm(sY~as.factor(D)+as.factor(A), + data=sbase,family=poisson(link="log")) + sYp=predict(sreg,type="response", + newdata=sbase) + sYpscenario=rqpois(36,sYp,phi=3.18623) + Rs[s]=sum(sYpscenario[is.na(Y)==TRUE]) + } si on regarde la densité de la distribution des paiements futurs, on obtient, plot(density(Rs)) avec en bleu la distribution obtenue sur les résidus de Pearson brut, et en simulant une loi de Poisson, et en rouge la distribution des provisions avec les deux modifications. On observe clairement que les quantiles (par exemple) ont fortement changé. On peut le voir encore plus précisément sur un box-plot Si on revient cinq minutes sur ce qu’on vient de faire. Dans le premier rappelons que l’on a besoin d’utiliser, pour prédire les paiements futurs, que (avec la filtration correspondant à la partie supérieure du triangle disponible), i.e. Comme on utilise un modèle de Poisson, on suppose aussi que i.e. c’est à dire que Il faut avoir des résidus centrés et de variance unitaire. A première vue, c’est le cas pour nos résidus de Pearson (à peu près), > mean(residuals(regp,"pearson")) [1] -0.02462518 > sd(residuals(regp,"pearson")) [1] 1.261934 sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur (classique en statistique quand on possède observations) > sqrt(sum((residuals(regp,"pearson")- + mean(residuals(regp,"pearson")))^2)/ + (length(residuals(regp,"pearson"))-1)) [1] 1.261934 Sauf qu’ici on ne perd pas un degré de liberté (remplacer l’espérance des observations par leur moyenne empirique): en régression, il faut corriger par le nombre de variables explicatives. Et donc, pour avoir un estimateur sans biais de la variance de nos résidus, il faut corriger par un facteur . D’où l’utilisation des résidus dits ajustés afin d’avoir des résidus de variance (vraiment) unitaire. Pour le second point, on utilise le fait qu’en pratique, la dispersion des paiements est plus grande que ce que donnerait un modèle de Poisson. Et il convient d’en tenir compte dans nos simulations. En l’occurrence, on veut que , mais aussi que. On a vu que le paramètre de surdispersion ne pouvait être supposé unitaire. On a alors le choix, entre simuler une loi Gamma de paramètres et ou bien simuler une loi binomiale négative, de moyenne , et de paramètre de surdispersion telle que la variance s’écrive Aussi, on prend alors que Avec ces méthodes, on obtient quelque chose… qui est programmé sous R, > BootChainLadder(PAID,10000,"od.pois") BootChainLadder(Triangle = PAID, R = 10000, process.distr = "od.pois") Latest Mean Ultimate Mean IBNR SD IBNR IBNR 75% IBNR 95% 1 4,456 4,456 0.0 0.0 0 0 2 4,730 4,752 22.2 12.2 28 44 3 5,420 5,456 35.5 15.2 43 64 4 6,020 6,086 65.7 19.8 78 102 5 6,794 6,946 151.9 28.7 170 202 6 5,217 7,364 2,147.4 111.0 2,218 2,339 Totals Latest: 32,637 Mean Ultimate: 35,060 Mean IBNR: 2,423 SD IBNR: 132 Total IBNR 75%: 2,506 Total IBNR 95%: 2,653 > quantile(Rs,c(.75,.95)) 75% 95% 2509 2653 et qui donne les mêmes résultats que ce qu’on vient de reprogrammer… # Cochrane, Pearson et le test du chi-deux En cours, nous avons poursuivi sur la loi multinomiale, et le test du chi-deux. Je voulais mettre un petit billet pour récapituler les différents points, et montrer une application numérique (nous reviendrons en détails mercredi sur des applications des outils vus jusqu’à présent). • Inférence avec la loi multinomiale On suppose qu’une variable peut prendre modalités, notées , avec . On posera en notant que appartient au simplexe de au sens où On a vu que l’estimateur du maximum de vraisemblance s’obtenait en faisant un peu d’optimisation sous contrainte, et que (en reprenant les notations du cours). On avait montré que et on a vu (ce qui peut se retrouver en introduisant la variable binomiale ). Mais plus généralement, on finira les calculs permettant d’établir que, pour (ce qui permet d’obtenir la matrice de variance covariance de ). En utilisant le théorème central limite on peut montrer que Sous une forme multivariée, on écrira avec et pour , . On notera que la somme de la ième colonne de est . Aussi, n’est pas inversible (c’est le fait que notre paramètre appartient au simplexe). Pour s’en sortir, la première idée est d’utiliser un peu d’algèbre linéaire. Une matrice est une matrice de projection si elle est idempotente, i.e. . Ses valeurs propres sont alors 0 ou 1, et si est le nombre de fois où 1 est valeur propre, et si , alors . Posons . Alors Or compte tenu de la forme de , qui est une matrice de projection dont la trace est (qui est aussi le nombre de fois où 1 est valeur propre). Donc Le test du chi-deux sera basé sur le fait qu’asymmptotiquement Une autre idée consiste à construire variables aléatoires qui seront indépendantes. Mais on peut plutôt regarder les applications de ce test, en particulier comme test d’indépendance. Pour information, Frank Yates a proposé un correction “pour continuité“, ici. La statistique considérée est alors • Retour sur le théorème de Cochrane Soit de dimension . Posons , pour , où on notera le rang de , en supposant que les sont positives semidefinies, alors on a équivalence entre • • pour , • les sont des variables indépendantes. Les s’interprètent comme des longueurs (euclidienne) de projections d’un vecteur Gaussien sur des sous-espaces orthogonaux (de dimension respective ). Si ces vecteurs sont indépendants, et suivent des lois du chi-deux à degrés de libertés, avec , alors les sous-espaces sont orthogonaux, et supplémentaires. On peut y voir une espèce d’extension du théorème de Pythagore, en remplaçant la notion de vecteurs orthogonaux par des variables indépendantes suivant des lois du chi-deux, et la somme des carrés des longueurs devient la somme des degrés de liberté. • Application comme test d’indépendance Considérons deux variables pouvant prendre toutes les deux deux modalités (disons deux lois binomiales). On est alors face a une loi multinomiale à 4 modalités • avec probabilité • avec probabilité • avec probabilité • avec probabilité Un test d’indépendance revient à tester si la loi multinomiale peut s’écrire pour des vecteurs et . tels que et . On a alors contraintes sur les paramètres. Ces deux contraintes additionnelles font que la statistique de test s’écrit qui va suivre asymptotiquement une loi i.e. d’après le théorème de Cochrane. • Peine de mort pour les condamnées pour meurtre en Floride 1976-1987 en fonction de la “race” du meurtrier et de la victime, • meurtrier de “race blanche” et victime de “race blanche“: 53 condamnés à mort, 414 non condamnés à mort • meurtrier de “race blanche” et victime de “race noire“: 0 condamné à mort, 16 non condamnés à mort • meurtrier de “race noire“et victime de “race blanche“: 11 condamnés à mort, 37 non condamnés à mort • meurtrier de “race noire“et victime de “race noire“: 4 condamnés à mort, 139 non condamnés à mort On peut alors faire des tests d’indépendance, entre la “race” du meurtrier et le verdict par exemple. MEURTRIER=matrix(c(53+0,11+4,414+16,139+37),2,2) VICTIME =matrix(c(53+11,0+4,414+37,139+16),2,2) n=sum(MEURTRIER) (PROBMEUTR=MEURTRIER/n) [,1] [,2] [1,] 0.07863501 0.6379822 [2,] 0.02225519 0.2611276 SL=rowSums(PROBMEUTR) SC=colSums(PROBMEUTR) (MEUTRINDEP=outer(SL, SC, "*")) [,1] [,2] [1,] 0.07229966 0.6443176 [2,] 0.02859055 0.2547922 (Q=n*sum((PROBMEUTR - MEUTRINDEP)^2/MEUTRINDEP)) [1] 1.468519 (Qcorrec=n*sum((abs(PROBMEUTR - MEUTRINDEP)-.5/n)^2/MEUTRINDEP)) [1] 1.144741 pchisq(Qcorrec, (2-1)*(2-1), lower.tail = FALSE) [1] 0.2846528 qchisq(.95, (2-1)*(2-1)) [1] 3.841459 chisq.test(MEURTRIER) Pearson's Chi-squared test with Yates' continuity correction data: MEURTRIER X-squared = 1.1447, df = 1, p-value = 0.2847 On rejette donc l’hypothèse d’indépendance. # Comprendre les tableaux croisés, 1 Je voulais faire quelques billets rapides sur les tableaux croisés, suite à des demandes consécutives à mon billet sur le sexe des bébés (ici). Pour commencer, je vais m’intéresser au cas où les variables supposées indépendantes (d’autres billets seront l’occasion de tester d’autres types d’hypothèses). Histoire de trouver des données un peu originales, on peut traîner sur internet. Etant tombé (presque) par hasard sur le fichier pdf ici, j’ai voulu m’intéresser à l’origine des enseignants-chercheurs recrutés lors de la campagne 2008 (comme le dit le titre du document en ligne sur le site du ministère). Prenons deux variables catégorielles géographiques, par exemple • , ville (au sens large) où l’on a soutenu sa thèse • , ville (au sens large) où l’on a obtenu un poste de maître de conférence Compte tenu des règles de non-recrutement local, on devrait pouvoir supposer que ces variables sont indépendantes. En fait, je ne raisonne pas ici par université, mais par ville, ce qui correspond à un recrutement local au sens large (passer de Paris 6 à Paris 7 sera ici considéré comme “local“). L’étude a été faite à partir des données disponibles ici et . Petite précaution d’usage: ce sont les données de la campagne 2008, avant que les comités de sélection ne soient mis en place. Il conviendrait de refaire l’étude une fois les nouvelles données publiques. On a deux bases, une en droit-économie-gestion, et une en sciences. J’insiste sur le fait que l’on a des données sur une année seulement, donc avec peu d’observations… Mais comme l’indique le titre, je vais parler des outils visualisation des tableaux croisés, pas faire un billet sur le recrutement local. Voilà pour les précautions d’usage… Pour visualiser ici, j’ai fait des regroupements, car avec 25 villes, ça ne tient pas dans mon billet (mais par la suite, l’analyse sera plus fine). Les données de base sont les suivantes, les lignes étant l’endroit où un doctorant a soutenu sa thèse (en économie ou droit), et en colonnes, l’endroit où a été recruté un maître de conférence,  Rennes Province (autre) Région parisienne total Rennes 5 8 0 13 Province (autre) 3 124 23 150 Région parisienne 1 13 49 63 total 9 145 72 226 Classiquement, on peut retraduire ces chiffres (qui sont des comptages) sous la forme de probabilités,  Rennes Province (autre) Région parisienne total Rennes 2,21% 3,54% 0,00% 5,75% Province (autre) 1,32% 54,86% 10,18% 66,37% Région parisienne 0.44% 5,75% 21,68% 27,88% total 3,98% 64,15% 31,86% 100,00% Les probabilités qui apparaissent sur les totaux sont les probabilités marginales. Aussi, au bout de la première ligne on trouve que 5,75% des doctorants qui ont eu un poste en 2008 avaient soutenu leur thèse à Rennes (ou dans la région, i.e. les universités de Rennes 1, de Rennes 2 et de Brest). De manière similaire, en bas de la première colonne, 3,98% des maîtres de conférences recrutés en 2008 ont été recrutés à Rennes. Deux autres tableaux classiques sont alors les profils ligne et profils colonnes, c’est à dire les probabilités conditionnelles,  Rennes Province (autre) Région parisienne total Rennes 38,46% 61,54% 0,00% 100,00% Province (autre) 2,00% 82,67% 15,33% 100,00% Région parisienne 1,58% 20,64% 77,78% 100,00% total 3,98% 64,15% 31,86% 100,00% par ligne, autrement dit, pour un étudiant qui a soutenu sa thèse à Rennes (et qui a eu un poste de Maître de Conférences), il avait 38,46% de chance de rester à Rennes (ou plutôt en Bretagne), et 61,54% de chance d’aller dans une autre ville de province. Côté colonnes,  Rennes Province (autre) Région parisienne total Rennes 55,56% 5,52% 0,00% 5,75% Province (autre) 33,33% 85,52% 31,94% 66,37% Région parisienne 11,11% 8,96% 68,06% 27,88% total 100,00% 100,00% 100,00% 100,00% Ici, on lit que pour un Maître de Conférences recruté sur la région parisienne, il y avait 68,06% de chances pour qu’il ait fait en région parisienne. Bref, à partir de ces tableaux, on peut essayer de tester l’hypothèse d’indépendance, et si cette dernière n’est pas vérifiée, on va essayer de comprendre pourquoi. Mais avant, rappelons que l’hypothèse d’indépendance se traduit par le fait que La valeur de gauche est la probabilité jointe (que l’on retrouve au coeur du second tableau), alors que les deux autres sont les probabilités marginales, qui apparaissent en bas, et à droite (dans les totaux). Bref, on pourrait construire le tableau des probabilités que l’on devrait avoir sous hypothèse d’indépendance,  Rennes Province (autre) Région parisienne total Rennes 0,22% 3,69% 1,83% 5,75% Province (autre) 2,64% 42,58% 21,14% 66,37% Région parisienne 1,11% 17,88% 8,88% 27,88% total 3,98% 64,15% 31,86% 100,00% que l’on peut retraduire en terme d’effectif espéré, toujours sous l’hypothèse d’indépendance. Pour ceux qui connaissent les copules, on change juste la structure de dépendance ici, pas les lois marginales (ou les effectifs marginaux) qui restent inchangés.  Rennes Province (autre) Région parisienne total Rennes 0,5 8,3 4,2 13 Province (autre) 6,0 96,2 47,8 150 Région parisienne 2,5 40,5 20 63 total 9 145 72 226 L’outils de base pour comparer ces données est la distance du chi-deux. La statistique du chi-deux est ou encore Karl Pearson avait été un des premiers à étudier ces grandeurs, sous des hypothèses de normalité (des sommes de lois normales indépendantes au carré se ramenant à des lois du chi-deux). On notera qu’une approche alternative est de supposer que les nombres (ou les effectifs) suivent une loi de Poisson. On peut alors se demander si Si on suppose avoir une loi de Poisson, on peut s’intéresser aux résidus de Pearson, autrement dit, on retrouve des grandeurs très proches de la contribution au chi-deux (en prenant le carré). On aussi le signe qui nous indique s’il y a trop peu de monde, ou au contraire par assez…. Si je travaille sur la base totale, soit environ 25 “villes“, on obtient > cont1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/contingence-universite_1.csv", + header=TRUE,sep=";") > cont2=cont1 > for(i in 1:ncol(cont2)){ + for(j in 1:nrow(cont2)){ + if(is.na(cont2[j,i]==TRUE)){cont2[j,i]=0} } } > X=as.matrix(cont2[1:nrow(cont2),2:ncol(cont2)]) > rownames(X)=as.character(cont2[,1]) > colnames(X)=names(cont2[2:ncol(cont2)]) > X=X[,-2] > X=X[-2,] > X=X[,-9] > X=X[-9,] > Z=X[-nrow(X),] > Z=Z[,-ncol(Z)] > Z=Z[-nrow(Z),] > chisq.test(Z) Pearson's Chi-squared test data: Z X-squared = 1342.546, df = 576, p-value < 2.2e-16 La base est assez sale, j’ai fait un copier/coller rapide du fichier pdf (et je vire les Antilles et la Réunion car après je voudrais faire un peu de visualisation sur une carte). Autrement dit, l’hypothèse d’indépendance ne semble pas validée par le test du chi-deux. Essayons maintenant de comprendre un peu mieux les interactions dans notre carré… Si on regarde les problèmes de recrutement local, il faut regarder ce qui se passe sur la diagonale de la matrice des résidus, > CT=chisq.test(Z) Message d'avis : In chisq.test(Z) : l'approximation du Chi-2 est peut-être incorrecte > Dia=diag(CT$residuals)
> names(Dia)=rownames(Z)
> Dia
AIX-MARSEILLE        BESANCON        BORDEAUX            CAEN        CLERMONT
6.43901593      3.26116568      8.37261202     10.53607373      5.12692873
CORSE           DIJON        GRENOBLE           LILLE         LIMOGES
8.56426340      4.10931024      9.18869705      7.59456681      8.56426340
LYON     MONTPELLIER           NANCY          NANTES            NICE
3.16119625      4.18366071      8.84765951      2.67087742      7.86973591
ORLEANS        POITIERS           REIMS          RENNES           ROUEN
1.79480842      5.81144228     -0.06651901      6.22963060      6.57435350
STRASBOURG        TOULOUSE           PARIS         CRETEIL      VERSAILLES
5.64850464      4.12289099      5.46092001     -0.22061860      4.01351350

Visuellement, sur une carte de France, on obtient la distribution ci-dessous, pour les postes en économie et droit,

et ci-dessous pour les postes en sciences et pharmacie,

Si on regarde le tableau des résidus on obtient les valeurs suivantes sur la ligne correspondant à Rennes, i.e. les endroits où vont les étudiants ayant soutenu leur thèse à Rennes,

avec en bleu les valeurs positives, et en rouge les valeurs négatives. A droite sont représentées les lignes, i.e. les endroits d’où viennent les Maîtres de Conférences recrutés à Rennes.
Au lieu de regarder ville par ville, on peut aussi faire une petite animation, soit sur les endroits où sont partis les doctorants,

ou les endroits d’où viennent les Maîtres de Conférence recrutés,

On peut également s’il existe des interactions, ou des échanges d’une ville à l’autre. Dans ce cas, les résidus croisées  et  doivent être positifs.

Bref, au delà du “recrutement local” on voit également apparaître de fortes interactions régionales (Rennes-Nantes, Toulouse-Bordeaux, Orléans-Paris).

Avant de conclure, il existe une autre technique permettant de mieux visualiser les interactions, c’est l’analyse factorielle des correspondances (je renvois icipour des slides). Si on regarde les 2 premiers axes principaux, en sciences (à gauche) et en droit et économie (à droite), on obtient les projections suivantes des modalités, avec en bleu les villes où les étudiants ont soutenu leur thèse, enrouge les villes où les postes de Maître de Conférences ont été obtenus,

et sur les les 2 axes suivants

Bref, comme souvent en statistique, afin de mieux comprendre un ajustement de modèle (ou comme ici lors la réalisation d’un test d’hypothèse) la visualisation des résidus est la technique de base (c’était l’idée d’ailleurs du commentaire que j’avais fait ).
Voilà pour la première partie, sur les modèles où l’on suppose l’indépendance. La prochaine fois, on verra d’autres formes d’hypothèses.

# Khi de la poule ou de l’oeuf…

Samedi soir, Djalil et Florent se demandaient ce que Karl Pearson avait inventé en premier, la loi du khi (ou du chi), ou la loi du khi-deux (pensant que la loi du chi-deux était relativement naturelle, et que la loi du khi découlait de la loi du khi-deux). La loi du khi admet pour densité

et est aussi parfois définie à partir de la fonction Gamma incomplète, i.e. sa fonction de répartition est

avec

et classiquement,

Pour la loi du khi-deux, sa densité est

dont la fonction de répartition est

Bref, on peut définir ces deux lois complètement indépendamment de leurs applications: si on considère des variables indépendantes , alors

suit une loi du khi-deux, alors que

suit une loi du khi. Bref, intuitivement, on définit d’abord la loi du chi-deux avant d’introduire la loi du khi, si l’on part de l’interprétation. Et si l’on reprend le texte fondamental de Fisher, ici, datant de 1922, on ne retrouve que le khi-deux

Mais dans “The lady drinking tea“, David Salsburg suggère que Karl Pearson avait introduit le khi avant le khi-deux… Si l’on y regarde de plus près, Karl Pearson a proposé le “test du chi-deux” en 1900 (ici), avec l’idée que l’exposant qui apparait dans la loi normale pouvait être noté, comme on le retrouve ci-dessous

On voit aussi interprété dans ce texte le khi, et pas seulement le khi-carré.

Effectivement, si l’on regarde ici par exemple, la loi du chi-deux est reliée à l’étude du volume d’un ellipsoïde. Mais ce n’est pas la première apparition de cette variation dite du khi-deux. Si on y regarde de plus près, par exemple ici, ou , on notera que Karl Pearson avait défini, dès 1896 aussi bien la loi du khi que la loi du khi-deux, dans le cadre de ce qui s’appelle les distribution de Pearson de type III,

Autrement dit, Karl Pearson a commencé par définir les deux lois simultanément, via les fonction gamma incomplète, et autres équations d’analyse fonctionnelle, et c’est seulement plus tard que l’utilisation de la loi du khi-deux l’a emporté, de part son interprétation relativement élégante (somme de carrés de lois normales centrées réduites) et le test d’indépendance (ou d’adéquation) du khi-deux.
Mais en 1931, Karl Pearson renvoie à un papier plus ancien où l’on trouve pour la première fois la loi du khi-deux, suggérant d’ailleurs d’appeler cette loi la loi de Helmert, en mémoire à Friedrich Robert Helmert, qui avait également beaucoup travailler sur les moindres carrés (dont il est naturel d’arriver à des sommes de carrés de lois normales, ici), qui aurait trouvé cette loi en 1875. Mais on peut également remonter aux travaux d’Ernst Karl Abbe en 1863 qui avait également mentionné la loi.

Mais c’est de Karl Pearson que l’on a hérité la notation et le nom du khi-deux.