Tag Archives: regression

De l’interprétation d’un effet nonlinéaire en régression

Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

  • Analyser un effet nonlinéaire

Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,

Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-1.png

notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-2.png

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-3.png

Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.

> m=mean(base$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")

On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)

Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.

 

  • Changer de modèle (et les comparer)

On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,

> REG0 = glm(nbsin~age,base,
+  family=poisson(link="log"),
+  offset=expo)
> REG1a = glm(nbsin~age,base,
+  family=binomial(link = "logit"),
+  offset=expo)
> REG1b = glm(nbsin~age,base,
+  family=binomial(link = "probit"),
+  offset=expo)
> REG1c = glm(nbsin~age,base, 
+  family=gaussian(link = "identity"),
+  offset=expo)

Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici

> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83

Mais on peut également regarder le critère BIC,

> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577

Autrement dit, en utilisant le règle “the smaller, the better” on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique1, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,

Bref, ce dernier modèle n’est pas convainquant… Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle… Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles…..

Propriétés des estimateurs dans une régression

Comme les rappels (je devrais plutôt dire “remise à niveau“) ont été particulièrement rapides, je vais prendre un peu de temps pour revenir sur quelques éléments du cours d’économétrie (d’autant plus que j’ai reçu quelques questions par mail).

  • sur le calcul des estimateurs

Je vais reprendre les questions que j’ai reçu, ça sera plus simple, “je comprend pas très bien quelle méthode utilise l’opération lm(Y~X…..) sur R pour trouver les coef” ou encore “je regresse un modèle simple, je sors le summary, et j’obtiens un betaREG 3 étoiles. Toujours sur R, je compare ce betaREG avec le betaCalc = (X’X)-1 * X’Y, il est différent ! Comment expliquer cela ? Est-ce possible ?“.
Commençons par la première question: la fonction lm() utilise la méthode dite des moindres carrés, ce qui revient à calculer

en adoptant une écriture matricielle, i.e. la solution est alors

Cet estimateur est celui qui minimise la somme des carrés des erreurs. Et on peut vérifier numériquement que cet estimateur est bien celui calculé par R,

>data(cars)
>X=cbind(rep(1,nrow(cars)),cars$speed); Y=cars$dist
> X[1:5,]
      [,1] [,2]
 [1,]    1    4
 [2,]    1    4
 [3,]    1    7
 [4,]    1    7
 [5,]    1    8
> solve(t(X)%*%X)
            [,1]         [,2]
[1,]  0.19310949 -0.011240876
[2,] -0.01124088  0.000729927
> t(X)%*%Y
      [,1]
[1,]  2149
[2,] 38482
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
           [,1]
[1,] -17.579095
[2,]   3.932409

ce qui correspond très précisément à ce que calcule la fonction de R

> lm(dist~speed,cars)
Coefficients:
(Intercept)        speed 
    -17.579        3.932

Et on peut aussi vérifier qu’il minimise bien la somme des carrés des erreurs,

> b0=seq(-30,10,by=3);
> b1=seq(1,7,by=.5)
> SC=matrix(NA,length(b0),length(b1))
> for(i in 1:length(b0)){
> for(j in 1:length(b1)){
 + SC[i,j]=sum((cars$dist-(b0[i]+b1[j]*cars$speed))^2)
}}
> contour(b0,b1,SC)

et je peux certifier que la valeur obtenue est bien le minimum de cette fonction (en fait c’est un résultat théorique d’optimisation).
Ca se généralise bien entendu en dimension plus grande. Alors attention, cette formule permet de calculer un estimateur du vecteur des paramètres, mais la formule est bien sûre fausse si on raisonne composante par composante.

> X=cars$speed; Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
         [,1]
[1,] 2.909132

Mais ça c’est un résultat (trivial) d’algèbre linéaire. En effet, le fait d’avoir

ne signifie pas que l’on ait ce genre de formule composante par composante, i.e.

Dit de manière un peu formelle, le produit matriciel n’est pas un produit terme à terme. La traduction économétrique de cette idée est qu’une régression multiple n’est pas une succession de régressions simples. En effet, dans l’équation précédante, le dernier terme correspond à la régression sans la constante, comme on peut le voir ci-dessous

> lm(dist~0+speed,cars)
Coefficients:
speed  
2.909

et le premier est tout simplement la moyenne empirique,

> lm(dist~1,cars)
Coefficients:
(Intercept)  
      42.98  
> mean(cars$dist)
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
      [,1]
[1,] 42.98
  • sur les propriétés des estimateurs

Une question était “doit-on vérifier les propriétés de l’estimateur (sans biais, efficace etc…) du summary de la Reg ? pour chaque régression ?
Bon, la réponse est non car c’est impossible ! Pour répéter tranquillement ce que j’avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d’hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l’estimateur par moindre carré est BLUE (best linear unbiased estimator).
Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d’être sans biais, par exemple.
On ne
peut pas
vérifier qu’un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l’absence de biais.
Mais comme on va le voir par la suite, le biais ce n’est pas forcément gênant, ce qui est important, c’est surtout la convergence. Comme le notait Clive Granger  “if you can’t even get a consistant estimator, you shouldn’t be in this business“.

  • et sur des points plus théoriques

sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un “rapport avec les plim X’E/n et plim (X’X)-1/ n“.
Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c’est à dire non corrélé avec la variable explicative, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-1.png

En effet, sinon

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-2.png

autrement dit l’estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce “problème” est d’utiliser un instrument (ou une variable instrumentale).
Mais surtout, sur cet exemple, on voit qu’on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement.
Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur).

  • un estimateur sans biais est généralement convergent (et de manière consistante)
  • un estimateur biaisé peut ne pas converger

Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l’ai expliqué auparavant, la distribution de l’estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d’aujourd’hui). La figure ci-dessous correspond au cas sans biais, et convergent

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-sans-biais-convergent_m.jpg

La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l’estimateur sera sans biais.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-convergent_m.jpg

Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l’estimateur sera toujours biaisé,
https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-non-convergent_m.jpg

Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-3.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-4.png

(c’est l’hypothèse d’homoscédasticité et d’indépendance des résidus). Bon, il faut aussi l’absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité.
On peut s’intéresser aux propriétés asymptotiques de l’estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-5.png converge en probabilité vers https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-6.png, parfois noté

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-7.png

si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-8.png

Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-9.png

si l’espérance des variables est finie, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-10.png

Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-11.png

On parlera alors d’estimateur convergent. On peut aussi montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-12.png

Mais formellement,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-13.png

si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-14.png

Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d’économétrie 1 qui commencera très bientôt… et surtout d’économétrie 2 qui abordera les aspects dynamiques.