On commence par charger les données utilisées. Il s’agit du même jeu de données que celui utilisé dans le premier bloc-notes, auquel on renvoie pour davantage d’explications sur les données.
load(url("https://irma.math.unistra.fr/~jberard/donnees_swautoins_ext_sim"))Les données ainsi récupérées se trouvent dans une variable appelée
donnees.
str(donnees)## 'data.frame': 239444 obs. of 7 variables:
## $ Kilometres: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ Zone : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Bonus : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Make : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Claims : int 0 1 3 3 6 1 1 1 2 2 ...
## $ Insured : num 6.36 7.37 6.81 6.56 13.53 ...
## $ Payment : num 0 5475 21866 15230 31354 ...On va ici s’intéresser à la variable Payment, qui
représente le montant total (en couronnes suédoises de 1977) de
sinistres associé à chaque observation (ici, une observation est
représentée par une ligne du jeu de données, et correspond à une ou
plusieurs polices d’assurance).
Plus précisément, nous allons ajuster un modèle log-Gamma pour le
coût individuel d’un sinistre, en fonction des variables
explicatives Kilometres, Bonus,
Zone et Make.
On commence par extraire du jeu de données les lignes correspondant
effectivement à des polices sinistrées, c’est-à-dire comportant un
nombre de sinistres non-nul. L’extraction est effectuée au moyen de la
commande subset, les lignes sélectionnées étant
caractérisées par une valeur strictement positive de la variable
Claims. On constate que le jeu de données ainsi obtenu est
fortement réduit par rapport au jeu de données intial, puisque l’on ne
conserve que 83980 lignes sur les 239444 initialement présentes.
donnees <- subset(donnees, subset = (Claims > 0))
str(donnees)## 'data.frame': 83980 obs. of 7 variables:
## $ Kilometres: Factor w/ 5 levels "1","2","3","4",..: 2 3 4 5 1 2 3 4 5 2 ...
## $ Zone : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 2 3 ...
## $ Bonus : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Make : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Claims : int 1 3 3 6 1 1 1 2 2 3 ...
## $ Insured : num 7.37 6.81 6.56 13.53 8.2 ...
## $ Payment : num 5475 21866 15230 31354 7508 ...Chaque ligne du jeu de données ainsi sélectionné correspond à un
nombre de sinistres individuels variable d’une ligne à l’autre, dont la
variable Payment indique le coût total, c’est-à-dire la
somme des coûts des sinistres individuels correspondants. Notre but
étant de modéliser le coût d’un sinistre individuel, on travaillera donc
avec les coûts de sinistres normalisés par ce nombre, en utilisant le
nombre de sinistres comme un poids dans l’ajustement du modèle GLM
correspondant. On ajoute à cette fin une variable nommée
Payments_n dans le jeu de données, qui contenant les coûts
ainsi normalisés.
donnees[["Payment_n"]] <- donnees[["Payment"]] / donnees[["Claims"]]
str(donnees)## 'data.frame': 83980 obs. of 8 variables:
## $ Kilometres: Factor w/ 5 levels "1","2","3","4",..: 2 3 4 5 1 2 3 4 5 2 ...
## $ Zone : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 2 3 ...
## $ Bonus : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Make : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Claims : int 1 3 3 6 1 1 1 2 2 3 ...
## $ Insured : num 7.37 6.81 6.56 13.53 8.2 ...
## $ Payment : num 5475 21866 15230 31354 7508 ...
## $ Payment_n : num 5475 7289 5077 5226 7508 ...On procède maintenant à l’ajustement du modèle GLM, en prenant soin
de choisir les coûts normalisés (Payment_n) comme variable
réponse, et les nombres de sinistres (Claims) comme
poids.
ajustement <- glm(formula = Payment_n ~ Kilometres + Zone + Bonus + Make,
family = Gamma(link = "log"),
weights = Claims,
data = donnees)
summary(ajustement)##
## Call:
## glm(formula = Payment_n ~ Kilometres + Zone + Bonus + Make, family = Gamma(link = "log"),
## data = donnees, weights = Claims)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.487525 0.007525 1127.849 < 2e-16 ***
## Kilometres2 0.004559 0.004124 1.106 0.26891
## Kilometres3 0.001646 0.004738 0.347 0.72830
## Kilometres4 0.003231 0.006589 0.490 0.62392
## Kilometres5 -0.004860 0.006979 -0.696 0.48617
## Zone2 0.025130 0.005209 4.824 1.41e-06 ***
## Zone3 0.044600 0.005309 8.401 < 2e-16 ***
## Zone4 0.135446 0.004748 28.528 < 2e-16 ***
## Zone5 0.056600 0.007942 7.127 1.03e-12 ***
## Zone6 0.151309 0.006526 23.187 < 2e-16 ***
## Zone7 0.045258 0.022049 2.053 0.04011 *
## Bonus2 0.008902 0.006631 1.343 0.17939
## Bonus3 -0.001263 0.007398 -0.171 0.86446
## Bonus4 0.014629 0.007959 1.838 0.06604 .
## Bonus5 0.001818 0.007710 0.236 0.81353
## Bonus6 -0.002833 0.006374 -0.444 0.65674
## Bonus7 0.003935 0.004758 0.827 0.40818
## Make2 -0.030249 0.011623 -2.603 0.00926 **
## Make3 0.079690 0.013809 5.771 7.92e-09 ***
## Make4 -0.211314 0.012994 -16.263 < 2e-16 ***
## Make5 -0.088844 0.011153 -7.966 1.66e-15 ***
## Make6 -0.057513 0.009525 -6.038 1.56e-09 ***
## Make7-9 -0.140698 0.005414 -25.989 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2993475)
##
## Null deviance: 26953 on 83979 degrees of freedom
## Residual deviance: 26217 on 83957 degrees of freedom
## AIC: 2039224
##
## Number of Fisher Scoring iterations: 4On constate que le paramètre de dispersion \(\phi\) estimé est égal à 0.2993475.
On peut extraire cette valeur estimée à l’aide de la commande ci-dessous.
summary(ajustement)[["dispersion"]]## [1] 0.2993475On examine maintenant de plus près l’estimation du paramètre de dispersion \(\phi\).
L’estimation précédente est obtenue à l’aide de l’estimateur classique \[\hat{\phi} = \frac{1}{N-(d+1)}\sum_{i=1}^N w_i \cdot \frac{(Y_i -\hat{\mu}_i)^2}{v(\hat{\mu}_i)}.\]
On peut retrouver cette valeur “à la main”, en extrayant directement
les résidus de Pearson \(\frac{(Y_i
-\hat{\mu}_i)}{\sqrt{v(\hat{\mu}_i)/w_i}}\) du modèle ajusté à
l’aide de la fonction residuals, et en calculant la somme
de leurs carrés.
stat_Pearson <- sum(residuals(ajustement, type = "pearson")^2)
stat_Pearson # Statistique de Pearson## [1] 25132.32On accède ensuite à la quantité \(N-(d+1)\), appelée “nombre de degrés de
liberté” du modèle ajusté, en faisant appel à la fonction
df.residual.
deg_lib <- df.residual(ajustement)
deg_lib # Nombre de degrés de liberté## [1] 83957On peut ainsi recalculer l’estimation du paramètre de dispersion, et constater que l’on obtient la même valeur que celle constatée plus haut.
phi_est <- stat_Pearson / deg_lib
phi_est # Valeur estimée de phi## [1] 0.2993475On peut aller un cran plus loin dans l’extraction “manuelle” des quantités précédentes, comme illustré maintenant. L’intérêt de ce qui suit est de montrer comment accéder individuellement aux éléments intervenant dans ce type de calculs.
On commence par extraire la fonction de variance \(v\) associée à la loi du modèle que l’on a ajusté. (Dans notre cas précis, nous savons que, pour la loi Gamma, la fonction de variance est \(v(\mu)=\mu^2\)), mais il est plus commode d’obtenir celle-ci automatiquement.)
fonc_variance <- family(ajustement)[["variance"]]On vérifie que la fonction ainsi extraite est bien celle attendue.
fonc_variance## function (mu)
## mu^2
## <bytecode: 0x7f8205f7e0c0>
## <environment: 0x7f8205f83048>On extrait maintenant les valeurs modélisées pour l’espérance des réponses normalisées (puisque c’est sur des réponses normalisées que le modèle a été ajusté) sur les données utilisées pour l’ajustement.
esp_mod_n <- fitted(ajustement)On récupère également les poids de normalisation utilisés (ici, les nombres de sinistres).
poids <- donnees[["Claims"]]Enfin, on récupère les valeurs des réponses normalisées dans les données d’ajustement.
reponse_n <- donnees[["Payment_n"]]On peut maintenant calculer la statistique de Pearson en appliquant élément par élément la formule qui la définit.
sum(poids * (reponse_n - esp_mod_n)^2 / fonc_variance(esp_mod_n))## [1] 25132.32Si l’on y tient, on peut encore extraire “à la main” le nombre de coefficients (\(d+1\)) et le nombre de lignes (\(N\)) du jeu de données utilisé pour l’ajustement.
N <- nrow(donnees)
d_plus_un <- length(coef(ajustement))
N-d_plus_un## [1] 83957Une autre approche possible pour estimer le paramètre de dispersion
est d’estimer celui-ci par maximum de vraisemblance, ce qui est par
exemple possible pour un GLM gamma en utilisant la fonction
gamma.shape du paquet MASS, qui fournit une
estimation du paramètre de forme \(k\),
appelé ici alpha (attention : on a \(\phi=1/k\)), avec non seulement une
estimation ponctuelle mais également un écart-type estimé permettant, si
on le souhaite, de construire des intervalles de confiance en utilisant
l’approximation normale. La composante alpha de l’objet
renvoyé par la fonction donne la valeur estimée de \(k\), tandis que la composante
SE donne la valeur estimée de l’écart-type de l’estimateur
de \(k\) fournissant cette valeur
estimée.
library(MASS)
forme_est <- gamma.shape(ajustement)
str(forme_est)## List of 2
## $ alpha: num 3.34
## $ SE : num 0.0156
## - attr(*, "class")= chr "gamma.shape"phi_est_bis <- 1 / forme_est[["alpha"]]
phi_est_bis # Valeur estimée de phi## [1] 0.2994553On constate que la valeur ainsi estimée pour \(\phi\) est voisine de celle obtenue avec l’approche basée sur la statistique de Pearson.
Une troisième approche non-conseillée utilise la déviance, en estimant \(\phi\) par le rapport \(D/(N-d+1)\). Ce calcul est effectué ci-après pour mémoire.
phi_est_ter <- deviance(ajustement)/df.residual(ajustement)
phi_est_ter # Valeur estimée de phi## [1] 0.3122628Dans le cas présent, la valeur obtenue n’est pas très éloignée de celle obtenue par les deux autres méthodes.
La déviance associée au modèle ajusté (appelée “residual deviance”)
peut être extraite facilement à l’aide de la fonction
deviance.
D <- deviance(ajustement)
D## [1] 26216.65On rappelle que celle-ci est définie par la formule \(D=\sum_{i=1}^N w_i \cdot D\left(y_i, \mu_i^{\mbox{est.}} \right)\), où \(D\) est la fonction de déviance associée à la loi utilisée par le modèle. On recalcule maintenant cette valeur “à la main”.
On commence par extraire la fonction de déviance dont on a besoin. (On sait qu’il s’agit de la fonction de déviance associée à la loi Gamma, soit \(D(y,\mu)=-2 \left[\log(y/\mu) - \frac{y-\mu}{\mu} \right]\), mais on l’extrait ici de manière systématique.)
fonc_deviance <- family(ajustement)[["dev.resids"]]
fonc_deviance## function (y, mu, wt)
## -2 * wt * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu)
## <bytecode: 0x7f8205f7d950>
## <environment: 0x7f8205f83048>La fonction ainsi extraite inclut comme argument non-seulement la
valeur de la réponse (y) et de l’espérance modélisée
(mu), mais également le poids (wt), de manière
à renvoyer directement la valeur \(w_i \cdot
D\left(y_i, \mu_i^{\mbox{est.}} \right)\) lorsque y
et mu sont les valeurs normalisées de la réponse et de
l’espérance modélisée, et wt le poids utilisé. On peut donc
re-calculer la déviance en sommant comme ci-dessous.
sum(fonc_deviance(y = reponse_n, mu = esp_mod_n, wt = poids))## [1] 26216.65On obtient bien une valeur identique à celle obtenue précédemment.
La déviance \(D_0\) (“null deviance”) est associée à un modèle ajusté en utilisant uniquement un terme constant, c’est-à-dire un modèle dont l’espérance modélisée est toujours égale à la moyenne globale des réponses observées, indépendamment des valeurs des variables explicatives.
Il n’y a pas de fonction dédiée pour extraire celle-ci, mais on peut
obtenir sa valeur en accédant au champ nommé null.deviance
de l’objet renvoyé par la fonction glm.
D_0 <- ajustement[["null.deviance"]]
D_0## [1] 26953.47(Remarque : dans ce qui précède, nous avons utilisé diverses
fonctions d’extraction : coef, family,
df.resid, deviance, fitted, etc.
pour accéder à des valeurs associées au modèle ajusté. Dans chacun de
ces cas, il était également possible d’accéder à ces valeurs à l’aide du
champ correspondant de l’objet renvoyé par la fonction glm
: ajustement[["coefficients"]],
ajustement[["family"]],
ajustement[["df.residual"]],
ajustement[["deviance"]],
ajustement[["fitted.values"]], etc. On peut néanmoins
considérer qu’il est préférable de passer par des fonctions d’extraction
dédiées, lorsqu’elles existent, plutôt que de manipuler la structure
interne des objets.)
Pour retrouver la valeur de \(D_0\), on calcule la moyenne globale des réponses.
rep_moyenne <- sum(reponse_n * poids) / sum(poids)On calcule ensuite la somme correspondante, et l’on vérifie que la valeur obtenue est bien identique à celle obtenue précédemment.
sum(fonc_deviance(y = reponse_n, mu = rep_moyenne, wt = poids))## [1] 26953.47On peut maintenant calculer le pseudo-\(R^2\) donné par la formule \(1-D_0/D\).
r_D_sq <- (1 - D / D_0)
r_D_sq## [1] 0.02733705La valeur obtenue (0.0273371) ne semble pas très élevée si on la compare à des références couramment mentionnées dans des contextes de modélisation différents. Cependant, lorsque le phénomène étudié présente un caractère fortement aléatoire à l’échelle d’une unité d’observation (ici, une ligne du jeu de données), il n’est pas surprenant que le modèle utilisé ne parvienne à rendre compte que d’une faible part de la variabilité de la réponse. En conservant le même modèle, mais en utilisant des unités d’observation associées à des expositions plus importantes, les valeurs obtenues peuvent changer considérablement.
On illustre cela ci-après, en regroupant les lignes du jeu de données possédant exactement la même combinaison de valeurs des variables explicatives.
donnees_groupe <- aggregate(x = cbind(Claims, Payment, Insured) ~ Kilometres + Zone + Bonus + Make,
data = donnees,
FUN = sum)
str(donnees_groupe)## 'data.frame': 1428 obs. of 7 variables:
## $ Kilometres: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ Zone : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Bonus : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Make : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Claims : num 73 94 46 10 12 69 109 50 8 7 ...
## $ Payment : num 336602 478366 201823 56473 67723 ...
## $ Insured : num 357.4 514.9 178.2 36.3 39.4 ...donnees_groupe[["Payment_n"]] <- donnees_groupe[["Payment"]] / donnees_groupe[["Claims"]]ajustement_groupe <- glm(formula = Payment_n ~ Kilometres + Zone + Bonus + Make,
family = Gamma(link = "log"),
weights = Claims,
data = donnees_groupe)
summary(ajustement_groupe)##
## Call:
## glm(formula = Payment_n ~ Kilometres + Zone + Bonus + Make, family = Gamma(link = "log"),
## data = donnees_groupe, weights = Claims)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.487525 0.007674 1105.978 < 2e-16 ***
## Kilometres2 0.004560 0.004206 1.084 0.2785
## Kilometres3 0.001646 0.004831 0.341 0.7334
## Kilometres4 0.003231 0.006719 0.481 0.6307
## Kilometres5 -0.004860 0.007117 -0.683 0.4948
## Zone2 0.025130 0.005312 4.730 2.47e-06 ***
## Zone3 0.044600 0.005414 8.238 3.96e-16 ***
## Zone4 0.135446 0.004842 27.975 < 2e-16 ***
## Zone5 0.056600 0.008099 6.989 4.27e-12 ***
## Zone6 0.151309 0.006655 22.737 < 2e-16 ***
## Zone7 0.045258 0.022485 2.013 0.0443 *
## Bonus2 0.008902 0.006762 1.317 0.1882
## Bonus3 -0.001263 0.007545 -0.167 0.8671
## Bonus4 0.014629 0.008116 1.802 0.0717 .
## Bonus5 0.001818 0.007862 0.231 0.8171
## Bonus6 -0.002833 0.006500 -0.436 0.6630
## Bonus7 0.003935 0.004852 0.811 0.4175
## Make2 -0.030249 0.011853 -2.552 0.0108 *
## Make3 0.079689 0.014082 5.659 1.85e-08 ***
## Make4 -0.211315 0.013251 -15.947 < 2e-16 ***
## Make5 -0.088844 0.011373 -7.812 1.10e-14 ***
## Make6 -0.057514 0.009713 -5.921 4.01e-09 ***
## Make7-9 -0.140698 0.005521 -25.485 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3113038)
##
## Null deviance: 1176.14 on 1427 degrees of freedom
## Residual deviance: 439.31 on 1405 degrees of freedom
## AIC: 1599345
##
## Number of Fisher Scoring iterations: 4r_D_sq_groupe <- 1 - ajustement_groupe[["deviance"]] / ajustement_groupe[["null.deviance"]]
r_D_sq_groupe## [1] 0.6264814