Fusion de catégories pour un modèle GLM
library(outilscoursglm)1 Fusion guidée par des tests statistiques
1.1 Récupération des données
On récupère les données, construites à partir d’un jeu de données
d’assurance automobile retraitées de CASdatasets, déjà utilisées dans de
précédents bloc-notes. On va s’intéresser ici à la variable explicative
DriverAge, qui a fait l’objet d’une catégorisation a
priori par tranches de 5 ans.
load(url("https://irma.math.unistra.fr/~jberard/freMTPLfreq_categ_age_cat_5"))
str(donnees) ## 'data.frame': 413169 obs. of 8 variables:
## $ ClaimNb : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Exposure : num 0.09 0.84 0.52 0.45 0.15 0.75 0.81 0.05 0.76 0.34 ...
## $ Power : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 1 1 1 3 ...
## $ CarAge : Factor w/ 2 levels "[0,15]","(15,Inf]": 1 1 1 1 1 1 1 1 1 1 ...
## $ DriverAge: Factor w/ 17 levels "[18,23]","(23,28]",..: 6 6 4 4 5 5 2 2 1 6 ...
## $ Brand : Factor w/ 2 levels "other","J": 2 2 2 2 2 2 2 2 1 2 ...
## $ Gas : Factor w/ 2 levels "Diesel","Regular": 1 1 2 2 1 1 2 2 2 2 ...
## $ Density : Factor w/ 5 levels "[0,40]","(40,200]",..: 2 2 4 4 2 2 4 4 5 5 ...1.2 Ajustement d’un modèle GLM log-Poisson
On ajuste un modèle GLM log-Poisson dont la réponse est la fréquence de sinistres, l’exposition étant gérée via un terme d’offset.
ajustement <- glm(ClaimNb ~
DriverAge +
CarAge +
Density +
Brand +
Power +
Gas +
offset(log(Exposure)),
family = poisson(link = "log"),
data = donnees)
summary(ajustement)##
## Call:
## glm(formula = ClaimNb ~ DriverAge + CarAge + Density + Brand +
## Power + Gas + offset(log(Exposure)), family = poisson(link = "log"),
## data = donnees)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.04246 0.04274 -47.790 < 2e-16 ***
## DriverAge(23,28] -0.67329 0.04057 -16.598 < 2e-16 ***
## DriverAge(28,33] -0.94742 0.03888 -24.370 < 2e-16 ***
## DriverAge(33,38] -1.02156 0.03853 -26.511 < 2e-16 ***
## DriverAge(38,43] -0.97440 0.03810 -25.577 < 2e-16 ***
## DriverAge(43,48] -0.86788 0.03742 -23.193 < 2e-16 ***
## DriverAge(48,53] -0.89575 0.03709 -24.147 < 2e-16 ***
## DriverAge(53,58] -1.04153 0.03997 -26.055 < 2e-16 ***
## DriverAge(58,63] -1.05345 0.04570 -23.049 < 2e-16 ***
## DriverAge(63,68] -1.12680 0.04981 -22.622 < 2e-16 ***
## DriverAge(68,73] -1.06917 0.05173 -20.668 < 2e-16 ***
## DriverAge(73,78] -1.06340 0.05859 -18.150 < 2e-16 ***
## DriverAge(78,83] -1.01245 0.07241 -13.982 < 2e-16 ***
## DriverAge(83,88] -0.95630 0.13162 -7.266 3.71e-13 ***
## DriverAge(88,93] -0.71956 0.19852 -3.625 0.000289 ***
## DriverAge(93,98] -1.09003 0.57812 -1.885 0.059367 .
## DriverAge(98,103] -0.87677 0.57815 -1.517 0.129390
## CarAge(15,Inf] -0.24372 0.03084 -7.902 2.75e-15 ***
## Density(40,200] 0.18180 0.02677 6.790 1.12e-11 ***
## Density(200,500] 0.31521 0.02969 10.617 < 2e-16 ***
## Density(500,4.5e+03] 0.52547 0.02609 20.144 < 2e-16 ***
## Density(4.5e+03,Inf] 0.65616 0.03584 18.309 < 2e-16 ***
## BrandJ -0.17610 0.02480 -7.102 1.23e-12 ***
## Power2 0.09167 0.02344 3.910 9.23e-05 ***
## Power3 0.24161 0.03014 8.015 1.10e-15 ***
## GasRegular -0.17135 0.01675 -10.232 < 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: 105613 on 413168 degrees of freedom
## Residual deviance: 103787 on 413143 degrees of freedom
## AIC: 135096
##
## Number of Fisher Scoring iterations: 6trace_terme_univar(ajustement = ajustement,
nom_variable = "DriverAge",
nom_terme = "DriverAge",
nom_expo = "Exposure",
trace_int_conf = TRUE,
trace_expo = TRUE)## Warning in trace_terme_univar(ajustement = ajustement, nom_variable = "DriverAge", : Utilisation des données de ajustement$data
En examinant le tracé des coefficients, on constate que certaines différences (notamment pour les âges élevés, où le volume de données est faible, mais pas uniquement) semblent faibles en comparaison des incertitudes reflétées par les intervalles de confiance individuels.
1.3 Fusion via le test de Wald sur les coefficients
On va envisager ici de fusionner certaines catégories, à l’aide d’une approche (très rudimentaire !) reposant sur des tests de Wald successifs. L’algorithme utilisé est “glouton” : initialement, toutes les catégories sont distinctes. Puis, étant donné le regroupement de catégories courant, on examine toutes les fusions possibles entre paires de catégories, et l’on retient la fusion pour laquelle la p-valeur du test de Wald est la plus élevée, tant que celle-ci reste au-dessus du seuil fixé. Chaque test effectué porte sur l’intégralité des fusions de catégories envisagées (pas seulement sur la paire en cours de test), à partir de la matrice de variance-covariance estimée dans le modèle initial (le modèle n’est pas ré-ajusté au fur et à mesure des fusions.)
On utilise pour cela la fonction fus_categ_Wald, en
spécifiant le caractère ordinal de la variable étudiée (seules des
fusions entre catégories adjacentes seront envisagées), et en
choisissant un seuil de p-valeur à \(0.05\).
corresp_recat <- fus_categ_Wald(ajustement = ajustement,
noms_variables = "DriverAge",
ordonne = "DriverAge",
seuil = 0.05)## Etape 1
## Etape 2
## Etape 3
## Etape 4
## Etape 5
## Etape 6
## Etape 7
## Etape 8
## Etape 9
## Etape 10
## Etape 11
## Etape 12
## Etape 13
## [1] "p_val"On examine maintenant le résultat obtenu. On constate qu’un nombre substantiel de fusions est proposé.
corresp_recat## $corresp
## $corresp$DriverAge
## $corresp$DriverAge$`[18,23]`
## [1] "[18,23]"
##
## $corresp$DriverAge$`(23,28]`
## [1] "(23,28]"
##
## $corresp$DriverAge$`(28,33]-(33,38]-(38,43]`
## [1] "(28,33]" "(33,38]" "(38,43]"
##
## $corresp$DriverAge$`(43,48]-(48,53]`
## [1] "(43,48]" "(48,53]"
##
## $corresp$DriverAge$`(53,58]-(58,63]-(63,68]-(68,73]-(73,78]-(78,83]-(83,88]-(88,93]-(93,98]-(98,103]`
## [1] "(53,58]" "(58,63]" "(63,68]" "(68,73]" "(73,78]" "(78,83]"
## [7] "(83,88]" "(88,93]" "(93,98]" "(98,103]"
##
##
##
## $cause_sortie
## [1] "p_val"On crée maintenant une fonction réalisant la fusion de catégories proposées, après un léger réarrangement des noms de catégories.
corresp_recat_DriverAge <- arrange_recat_cut(corresp_recat[["corresp"]][["DriverAge"]])[["recat"]]
corresp_recat_DriverAge ## $`[18,23]`
## [1] "[18,23]"
##
## $`(23,28]`
## [1] "(23,28]"
##
## $`(28,43]`
## [1] "(28,33]" "(33,38]" "(38,43]"
##
## $`(43,53]`
## [1] "(43,48]" "(48,53]"
##
## $`(53,103]`
## [1] "(53,58]" "(58,63]" "(63,68]" "(68,73]" "(73,78]" "(78,83]"
## [7] "(83,88]" "(88,93]" "(93,98]" "(98,103]"f_recat_DriverAge <- cree_fonction_recat(corresp_recat_DriverAge)On crée maintenant la variable recatégorisée correspondante.
donnees[["DriverAge_recat"]] <- f_recat_DriverAge(donnees[["DriverAge"]])Et l’on réajuste le modèle avec cette nouvelle variable en lieu et place de la précédente.
ajustement_2 <- glm(ClaimNb ~
DriverAge_recat +
CarAge +
Density +
Brand +
Power +
Gas +
offset(log(Exposure)),
family = poisson(link = "log"),
data = donnees)
summary(ajustement_2)##
## Call:
## glm(formula = ClaimNb ~ DriverAge_recat + CarAge + Density +
## Brand + Power + Gas + offset(log(Exposure)), family = poisson(link = "log"),
## data = donnees)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.04237 0.04274 -47.789 < 2e-16 ***
## DriverAge_recat(23,28] -0.67303 0.04056 -16.592 < 2e-16 ***
## DriverAge_recat(28,43] -0.98156 0.03320 -29.567 < 2e-16 ***
## DriverAge_recat(43,53] -0.88183 0.03390 -26.010 < 2e-16 ***
## DriverAge_recat(53,103] -1.05818 0.03384 -31.268 < 2e-16 ***
## CarAge(15,Inf] -0.24304 0.03083 -7.883 3.20e-15 ***
## Density(40,200] 0.18253 0.02676 6.820 9.09e-12 ***
## Density(200,500] 0.31523 0.02968 10.620 < 2e-16 ***
## Density(500,4.5e+03] 0.52633 0.02606 20.195 < 2e-16 ***
## Density(4.5e+03,Inf] 0.65683 0.03581 18.341 < 2e-16 ***
## BrandJ -0.17741 0.02477 -7.161 8.01e-13 ***
## Power2 0.09012 0.02344 3.845 0.00012 ***
## Power3 0.23949 0.03010 7.956 1.78e-15 ***
## GasRegular -0.17029 0.01665 -10.228 < 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: 105613 on 413168 degrees of freedom
## Residual deviance: 103800 on 413155 degrees of freedom
## AIC: 135085
##
## Number of Fisher Scoring iterations: 6On constate une amélioration du critère AIC (diminution d’environ \(11\) entre le modèle initial et le modèle après recatégorisation).
Le tracé des coefficients permet de visualiser les effets attribués à
la variable DriverAge recatégorisée.
trace_terme_univar(ajustement = ajustement_2,
nom_variable = "DriverAge_recat",
nom_terme = "DriverAge_recat",
nom_expo = "Exposure",
trace_int_conf = TRUE,
trace_expo = TRUE)## Warning in trace_terme_univar(ajustement = ajustement_2, nom_variable = "DriverAge_recat", : Utilisation des données de ajustement$data
Une manière alternative de procéder consiste à laisser la variable recatégorisée sous la forme explicite d’une fonction de la précédente dans la formule.
ajustement_2_bis <- glm(ClaimNb ~ f_recat_DriverAge(DriverAge) +
CarAge +
Density +
Brand +
Power +
Gas +
offset(log(Exposure)),
family = poisson(link = "log"),
data = donnees)
summary(ajustement_2_bis)##
## Call:
## glm(formula = ClaimNb ~ f_recat_DriverAge(DriverAge) + CarAge +
## Density + Brand + Power + Gas + offset(log(Exposure)), family = poisson(link = "log"),
## data = donnees)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.04237 0.04274 -47.789 < 2e-16 ***
## f_recat_DriverAge(DriverAge)(23,28] -0.67303 0.04056 -16.592 < 2e-16 ***
## f_recat_DriverAge(DriverAge)(28,43] -0.98156 0.03320 -29.567 < 2e-16 ***
## f_recat_DriverAge(DriverAge)(43,53] -0.88183 0.03390 -26.010 < 2e-16 ***
## f_recat_DriverAge(DriverAge)(53,103] -1.05818 0.03384 -31.268 < 2e-16 ***
## CarAge(15,Inf] -0.24304 0.03083 -7.883 3.20e-15 ***
## Density(40,200] 0.18253 0.02676 6.820 9.09e-12 ***
## Density(200,500] 0.31523 0.02968 10.620 < 2e-16 ***
## Density(500,4.5e+03] 0.52633 0.02606 20.195 < 2e-16 ***
## Density(4.5e+03,Inf] 0.65683 0.03581 18.341 < 2e-16 ***
## BrandJ -0.17741 0.02477 -7.161 8.01e-13 ***
## Power2 0.09012 0.02344 3.845 0.00012 ***
## Power3 0.23949 0.03010 7.956 1.78e-15 ***
## GasRegular -0.17029 0.01665 -10.228 < 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: 105613 on 413168 degrees of freedom
## Residual deviance: 103800 on 413155 degrees of freedom
## AIC: 135085
##
## Number of Fisher Scoring iterations: 6Le modèle obtenu est équivalent au précédent, mais on peut tracer les effets en se référant toujours à la variable précédente.
trace_terme_univar(ajustement = ajustement_2_bis,
nom_variable = "DriverAge",
nom_terme = "f_recat_DriverAge(DriverAge)",
nom_expo = "Exposure",
trace_int_conf = TRUE,
trace_expo = TRUE)## Warning in trace_terme_univar(ajustement = ajustement_2_bis, nom_variable = "DriverAge", : Utilisation des données de ajustement$data
2 Fusion par maximum de vraisemblance pénalisé LASSO
2.1 Récupération des données
On travaille de nouveau avec le jeu de données des loyer à Munich.
load(url("https://irma.math.unistra.fr/~jberard/rent_Munich")) # Données retraitées de catdata
str(donnees)## 'data.frame': 2053 obs. of 13 variables:
## $ rent : num 10.9 11.01 8.38 8.52 6.98 ...
## $ size : int 68 65 63 65 100 81 55 79 52 77 ...
## $ rooms : Factor w/ 6 levels "1","2","3","4",..: 2 2 3 3 4 4 2 3 1 3 ...
## $ year : Factor w/ 10 levels "1910","1920",..: 1 9 1 8 9 8 2 2 5 4 ...
## $ area : Factor w/ 25 levels "1","2","3","4",..: 2 2 2 16 16 16 6 6 6 6 ...
## $ good : Factor w/ 2 levels "0","1": 2 2 2 1 2 1 1 1 1 1 ...
## $ best : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ warm : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ central : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ tiles : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ bathextra: Factor w/ 2 levels "0","1": 1 1 1 2 2 1 2 1 1 1 ...
## $ kitchen : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 1 1 ...
## $ standing : Factor w/ 3 levels "0","1","2": 2 2 2 1 2 1 1 1 1 1 ...Afin de travailler uniquement avec des variables catégorielles, on
catégorise (en déciles) la variable surface size.
# Catégorisation en déciles
donnees[["size"]] <- cut(x = donnees[["size"]],
breaks = quantile(donnees[["size"]], seq(0, 1, length.out = 11)),
include.lowest = TRUE)Afin de pouvoir comparer les résultats, on commence par ajuster un modèle GLM avec les variables telles qu’elles se présentent initialement.
ajustement <- glm(rent ~ size +
rooms +
year +
area +
warm +
central +
tiles +
bathextra +
kitchen +
standing,
data = donnees)
summary(ajustement)##
## Call:
## glm(formula = rent ~ size + rooms + year + area + warm + central +
## tiles + bathextra + kitchen + standing, data = donnees)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.15369 0.36306 30.721 < 2e-16 ***
## size(39,50] -1.54642 0.23965 -6.453 1.37e-10 ***
## size(50,55] -1.93748 0.27880 -6.949 4.95e-12 ***
## size(55,62] -1.95185 0.27493 -7.099 1.73e-12 ***
## size(62,67] -2.23771 0.29492 -7.587 4.97e-14 ***
## size(67,73] -2.17995 0.30036 -7.258 5.61e-13 ***
## size(73,80] -2.17080 0.30736 -7.063 2.25e-12 ***
## size(80,87] -2.12216 0.32697 -6.490 1.08e-10 ***
## size(87,100] -2.33471 0.33065 -7.061 2.27e-12 ***
## size(100,185] -2.62680 0.35537 -7.392 2.12e-13 ***
## rooms2 0.01668 0.22877 0.073 0.941872
## rooms3 -0.13957 0.27230 -0.513 0.608313
## rooms4 -0.45972 0.31821 -1.445 0.148703
## rooms5 -0.45110 0.43388 -1.040 0.298617
## rooms6 -0.73676 0.63072 -1.168 0.242901
## year1920 -1.22295 0.26704 -4.580 4.94e-06 ***
## year1930 -0.93486 0.54604 -1.712 0.087039 .
## year1940 -0.93334 0.23238 -4.016 6.13e-05 ***
## year1950 -0.32890 0.16523 -1.991 0.046662 *
## year1960 0.12480 0.16200 0.770 0.441166
## year1970 0.38932 0.17898 2.175 0.029732 *
## year1980 1.15170 0.19881 5.793 8.02e-09 ***
## year1990 1.65941 0.19769 8.394 < 2e-16 ***
## year2000 1.71908 0.43650 3.938 8.49e-05 ***
## area2 -0.53763 0.34506 -1.558 0.119375
## area3 -0.20928 0.35544 -0.589 0.556073
## area4 -0.59637 0.35254 -1.692 0.090867 .
## area5 -0.60158 0.34982 -1.720 0.085650 .
## area6 -1.00296 0.40095 -2.501 0.012448 *
## area7 -1.62264 0.40179 -4.039 5.58e-05 ***
## area8 -1.09134 0.40740 -2.679 0.007450 **
## area9 -0.86820 0.34276 -2.533 0.011386 *
## area10 -1.14872 0.41779 -2.750 0.006022 **
## area11 -1.69309 0.40374 -4.194 2.87e-05 ***
## area12 -0.57155 0.38223 -1.495 0.134994
## area13 -0.81200 0.37648 -2.157 0.031137 *
## area14 -1.78812 0.41033 -4.358 1.38e-05 ***
## area15 -1.28227 0.44232 -2.899 0.003784 **
## area16 -1.88694 0.37085 -5.088 3.96e-07 ***
## area17 -1.28408 0.40134 -3.199 0.001398 **
## area18 -0.57046 0.38690 -1.474 0.140520
## area19 -1.39464 0.37318 -3.737 0.000191 ***
## area20 -1.21542 0.42413 -2.866 0.004205 **
## area21 -1.29668 0.41007 -3.162 0.001590 **
## area22 -1.87391 0.52049 -3.600 0.000326 ***
## area23 -1.59991 0.62316 -2.567 0.010318 *
## area24 -1.83534 0.49202 -3.730 0.000197 ***
## area25 -1.31770 0.36737 -3.587 0.000343 ***
## warm1 -2.01921 0.28315 -7.131 1.38e-12 ***
## central1 -1.33502 0.19420 -6.874 8.29e-12 ***
## tiles1 -0.54654 0.11667 -4.684 3.00e-06 ***
## bathextra1 0.49929 0.16129 3.096 0.001991 **
## kitchen1 1.19925 0.17783 6.744 2.02e-11 ***
## standing1 0.40122 0.11274 3.559 0.000381 ***
## standing2 1.45645 0.31425 4.635 3.81e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 3.912073)
##
## Null deviance: 12486.0 on 2052 degrees of freedom
## Residual deviance: 7816.3 on 1998 degrees of freedom
## AIC: 8682.8
##
## Number of Fisher Scoring iterations: 2On va utiliser le paquet smurf pour l’inférence avec
pénalisation LASSO.
library(smurf)On commence par spécifier la formule caractérisant l’ajustement souhaité :
pénalité de type “fused lasso” aux variables catégorielles ordinales (pénalisation de la différence entre coefficients associés à des modalités consécutives)
pénalité de type “lasso” aux variables catégorielles binaires (pénalisation de l’unique coefficient)
pénalité de type “generalized fused lasso” à l’unique variable purement catégorielle
area(pénalisation de toutes les différences entre coefficients).
formule <- rent ~
p(size, pen = "flasso") +
p(rooms,pen = "flasso") +
p(year, pen = "flasso") +
p(area, pen = "gflasso") +
p(warm, pen = "flasso") +
p(central, pen = "flasso") +
p(tiles, pen = "flasso") +
p(bathextra, pen = "flasso") +
p(kitchen, pen = "flasso") +
p(standing, pen = "flasso")On effectue ensuite l’ajustement, en spécifiant un certain nombre de
paramètres : la combinaison de lambda="is.aic" et
lambda.reest=TRUE fait que la valeur du paramètre de
pénalisation lambda va être effectuée en optimisant la
valeur de l’AIC du modèle GLM ré-estimé à partir des fusions de
catégories produites lors de l’ajustement pénalisé. L’utilisation de
“glm.stand” standardise les poids appliqués pour les termes de pénalité
à partir d’une estimation initiale des coefficients du GLM
(non-pénalisé).
ajustement_lasso <- glmsmurf(formula = formule,
family = gaussian(),
data = donnees,
pen.weights = "glm.stand",
lambda = "is.aic",
control = list(lambda.reest = TRUE))summary(ajustement_lasso)##
## Call: glmsmurf(formula = formule, family = gaussian(), data = donnees,
## lambda = "is.aic", pen.weights = "glm.stand", control = list(lambda.reest = TRUE))
##
## Deviance residuals of estimated model:
## Min. 1st Qu. Median 3rd Qu. Max.
## -6.73645 -1.26900 -0.00965 1.25494 8.90310
##
## Deviance residuals of re-estimated model:
## Min. 1st Qu. Median 3rd Qu. Max.
## -6.5682 -1.2127 -0.0181 1.2592 8.5690
##
##
## Coefficients:
## Estimated Re-estimated
## Intercept 10.6600 11.0720
## size(39,50] -1.6732 -1.5464
## size(50,55] -2.0172 -1.9632
## size(55,62] -2.0172 -1.9632
## size(62,67] -2.2777 -2.3030
## size(67,73] -2.2777 -2.3030
## size(73,80] -2.2777 -2.3030
## size(80,87] -2.2777 -2.3030
## size(87,100] -2.3073 -2.4928
## size(100,185] -2.3616 -2.7943
## rooms2 * *
## rooms3 * *
## rooms4 -0.2537 -0.3278
## rooms5 -0.2537 -0.3278
## rooms6 -0.2537 -0.3278
## year1920 -0.9019 -1.0810
## year1930 -0.9019 -1.0810
## year1940 -0.9019 -1.0810
## year1950 -0.1901 -0.3711
## year1960 0.1003 0.0728
## year1970 0.1317 0.3469
## year1980 1.2073 1.1245
## year1990 1.5781 1.6513
## year2000 1.5781 1.6513
## area2 -0.3013 -0.4287
## area3 * *
## area4 -0.3013 -0.4287
## area5 -0.3013 -0.4287
## area6 -0.6784 -1.1076
## area7 -0.9050 -1.6387
## area8 -0.6784 -1.1076
## area9 -0.4107 -0.7080
## area10 -0.6784 -1.1076
## area11 -0.9050 -1.6387
## area12 -0.3013 -0.4287
## area13 -0.4107 -0.7080
## area14 -0.9050 -1.6387
## area15 -0.6784 -1.1076
## area16 -0.9050 -1.6387
## area17 -0.6784 -1.1076
## area18 -0.3013 -0.4287
## area19 -0.6784 -1.1076
## area20 -0.6784 -1.1076
## area21 -0.6784 -1.1076
## area22 -0.9050 -1.6387
## area23 -0.9050 -1.6387
## area24 -0.9050 -1.6387
## area25 -0.6784 -1.1076
## warm1 -1.9562 -2.0509
## central1 -1.2423 -1.3296
## tiles1 -0.4160 -0.5433
## bathextra1 0.0196 0.4945
## kitchen1 1.0584 1.1853
## standing1 0.4266 0.3762
## standing2 1.1314 1.4360
## ---
## '*' indicates a zero coefficient
##
##
## Null deviance: 12486.047 on 2052 degrees of freedom
## -------------------
## Estimated model:
##
## Residual deviance: 8026.104 on 2029 degrees of freedom
## AIC: 8675.216; BIC: 8810.265; GCV score: 4.002
## Penalized minus log-likelihood: 2.143388
## -------------------
## Re-estimated model:
##
## Residual deviance: 7839.303 on 2029 degrees of freedom
## AIC: 8626.869; BIC: 8761.919; GCV score: 3.909
## Objective function: 2.188856
## -------------------
## lambda: 0.02072845; lambda1: 0; lambda2: 0
## Number of iterations: 199
## Final step size: 0.8019531
## Convergence: Succesful convergenceOn constate la réduction très nette de la valeur du critère AIC par rapport au modèle non-fusionné, avec une valeur de \(8634\) environ, contre \(8682.8\) pour le modèle non-fusionné.
plot_reest(ajustement_lasso)
On peut manipuler “manuellement” les re-catégorisations ainsi
obtenues. On utilise ci-après la fonction (expérimentale, caveat emptor
!) categ_smurf pour extraire les regroupements
produits.
recat_lasso <- categ_smurf(ajustement_lasso)
recat_lasso## $size
## $size$`[17,39]`
## [1] "[17,39]"
##
## $size$`(39,50]`
## [1] "(39,50]"
##
## $size$`(50,55]-(55,62]`
## [1] "(50,55]" "(55,62]"
##
## $size$`(62,67]-(67,73]-(73,80]-(80,87]`
## [1] "(62,67]" "(67,73]" "(73,80]" "(80,87]"
##
## $size$`(87,100]`
## [1] "(87,100]"
##
## $size$`(100,185]`
## [1] "(100,185]"
##
##
## $rooms
## $rooms$`1-2-3`
## [1] "1" "2" "3"
##
## $rooms$`4-5-6`
## [1] "4" "5" "6"
##
##
## $year
## $year$`1910`
## [1] "1910"
##
## $year$`1920-1930-1940`
## [1] "1920" "1930" "1940"
##
## $year$`1950`
## [1] "1950"
##
## $year$`1960`
## [1] "1960"
##
## $year$`1970`
## [1] "1970"
##
## $year$`1980`
## [1] "1980"
##
## $year$`1990-2000`
## [1] "1990" "2000"
##
##
## $area
## $area$`1-3`
## [1] "1" "3"
##
## $area$`2-4-5-12-18`
## [1] "2" "4" "5" "12" "18"
##
## $area$`6-8-10-15-17-19-20-21-25`
## [1] "6" "8" "10" "15" "17" "19" "20" "21" "25"
##
## $area$`7-11-14-16-22-23-24`
## [1] "7" "11" "14" "16" "22" "23" "24"
##
## $area$`9-13`
## [1] "9" "13"
##
##
## $warm
## $warm$`0`
## [1] "0"
##
## $warm$`1`
## [1] "1"
##
##
## $central
## $central$`0`
## [1] "0"
##
## $central$`1`
## [1] "1"
##
##
## $tiles
## $tiles$`0`
## [1] "0"
##
## $tiles$`1`
## [1] "1"
##
##
## $bathextra
## $bathextra$`0`
## [1] "0"
##
## $bathextra$`1`
## [1] "1"
##
##
## $kitchen
## $kitchen$`0`
## [1] "0"
##
## $kitchen$`1`
## [1] "1"
##
##
## $standing
## $standing$`0`
## [1] "0"
##
## $standing$`1`
## [1] "1"
##
## $standing$`2`
## [1] "2"On peut ensuite créer un jeu de données recatégorisé.
On crée ci-dessous une liste de fonctions de re-catégorisation (une par variable).
f_recat <- list()
for(boucle in (1:length(recat_lasso))) {
f_recat[[names(recat_lasso)[boucle]]] <- cree_fonction_recat(recat_lasso[[boucle]])
}On crée ensuite (automatiquement, plutôt qu’à la main) la formule correspondante.
explicatives <- lapply(X = names(recat_lasso),
FUN = function(x){paste('f_recat[[', '"', x, '"]](', x, ')', sep = "")})
droite <- do.call(what = paste, args = c(explicatives, sep = "+"))
gauche <- "rent ~"
formule_glm <- as.formula(paste(gauche, droite, sep = ""))
formule_glm## rent ~ f_recat[["size"]](size) + f_recat[["rooms"]](rooms) +
## f_recat[["year"]](year) + f_recat[["area"]](area) + f_recat[["warm"]](warm) +
## f_recat[["central"]](central) + f_recat[["tiles"]](tiles) +
## f_recat[["bathextra"]](bathextra) + f_recat[["kitchen"]](kitchen) +
## f_recat[["standing"]](standing)ajustement_glm <- glm(formule_glm, data = donnees)
summary(ajustement_glm)##
## Call:
## glm(formula = formule_glm, data = donnees)
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 11.07196 0.24726
## f_recat[["size"]](size)(39,50] -1.54640 0.18850
## f_recat[["size"]](size)(50,55]-(55,62] -1.96317 0.17040
## f_recat[["size"]](size)(62,67]-(67,73]-(73,80]-(80,87] -2.30298 0.15600
## f_recat[["size"]](size)(87,100] -2.49278 0.21546
## f_recat[["size"]](size)(100,185] -2.79427 0.24411
## f_recat[["rooms"]](rooms)4-5-6 -0.32780 0.16222
## f_recat[["year"]](year)1920-1930-1940 -1.08099 0.17886
## f_recat[["year"]](year)1950 -0.37108 0.15921
## f_recat[["year"]](year)1960 0.07279 0.14933
## f_recat[["year"]](year)1970 0.34690 0.16742
## f_recat[["year"]](year)1980 1.12449 0.18951
## f_recat[["year"]](year)1990-2000 1.65129 0.18035
## f_recat[["area"]](area)2-4-5-12-18 -0.42875 0.17424
## f_recat[["area"]](area)6-8-10-15-17-19-20-21-25 -1.10762 0.18494
## f_recat[["area"]](area)7-11-14-16-22-23-24 -1.63866 0.20392
## f_recat[["area"]](area)9-13 -0.70796 0.19874
## f_recat[["warm"]](warm)1 -2.05092 0.27800
## f_recat[["central"]](central)1 -1.32961 0.19159
## f_recat[["tiles"]](tiles)1 -0.54330 0.11536
## f_recat[["bathextra"]](bathextra)1 0.49455 0.15865
## f_recat[["kitchen"]](kitchen)1 1.18533 0.17478
## f_recat[["standing"]](standing)1 0.37624 0.10282
## f_recat[["standing"]](standing)2 1.43597 0.30753
## t value Pr(>|t|)
## (Intercept) 44.779 < 2e-16 ***
## f_recat[["size"]](size)(39,50] -8.204 4.08e-16 ***
## f_recat[["size"]](size)(50,55]-(55,62] -11.521 < 2e-16 ***
## f_recat[["size"]](size)(62,67]-(67,73]-(73,80]-(80,87] -14.763 < 2e-16 ***
## f_recat[["size"]](size)(87,100] -11.569 < 2e-16 ***
## f_recat[["size"]](size)(100,185] -11.447 < 2e-16 ***
## f_recat[["rooms"]](rooms)4-5-6 -2.021 0.043436 *
## f_recat[["year"]](year)1920-1930-1940 -6.044 1.79e-09 ***
## f_recat[["year"]](year)1950 -2.331 0.019865 *
## f_recat[["year"]](year)1960 0.487 0.626000
## f_recat[["year"]](year)1970 2.072 0.038389 *
## f_recat[["year"]](year)1980 5.934 3.48e-09 ***
## f_recat[["year"]](year)1990-2000 9.156 < 2e-16 ***
## f_recat[["area"]](area)2-4-5-12-18 -2.461 0.013950 *
## f_recat[["area"]](area)6-8-10-15-17-19-20-21-25 -5.989 2.49e-09 ***
## f_recat[["area"]](area)7-11-14-16-22-23-24 -8.036 1.56e-15 ***
## f_recat[["area"]](area)9-13 -3.562 0.000376 ***
## f_recat[["warm"]](warm)1 -7.377 2.34e-13 ***
## f_recat[["central"]](central)1 -6.940 5.26e-12 ***
## f_recat[["tiles"]](tiles)1 -4.710 2.65e-06 ***
## f_recat[["bathextra"]](bathextra)1 3.117 0.001851 **
## f_recat[["kitchen"]](kitchen)1 6.782 1.55e-11 ***
## f_recat[["standing"]](standing)1 3.659 0.000259 ***
## f_recat[["standing"]](standing)2 4.669 3.22e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 3.863629)
##
## Null deviance: 12486.0 on 2052 degrees of freedom
## Residual deviance: 7839.3 on 2029 degrees of freedom
## AIC: 8626.9
##
## Number of Fisher Scoring iterations: 2trace_terme_univar(ajustement = ajustement_glm,
nom_variable = "size",
nom_terme = 'f_recat[["size"]](size)',
trace_int_conf = TRUE,
trace_expo = TRUE)## Warning in trace_terme_univar(ajustement = ajustement_glm, nom_variable = "size", : Utilisation des données de ajustement$data## Warning in trace_terme_univar(ajustement = ajustement_glm, nom_variable = "size", : Expositions fixées à 1
trace_terme_univar(ajustement = ajustement_glm,
nom_variable = "year",
nom_terme = 'f_recat[["year"]](year)',
trace_int_conf = TRUE,
trace_expo = TRUE)## Warning in trace_terme_univar(ajustement = ajustement_glm, nom_variable = "year", : Utilisation des données de ajustement$data## Warning in trace_terme_univar(ajustement = ajustement_glm, nom_variable = "year", : Expositions fixées à 1
Variante : pour la variable area, on pénalise seulement les différences de coefficients entre quartiers adjacents
load(url("https://irma.math.unistra.fr/~jberard/adj_mat_Munich"))
adj_mat_Munich # Matrice d'adjacence des quartiers## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## 1 0 1 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0
## 2 1 0 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## 3 1 1 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 0 0 0 0 0 0
## 6 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0
## 7 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1
## 8 0 1 1 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## 9 0 0 1 1 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1
## 10 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0
## 11 0 0 0 1 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
## 12 1 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## 13 1 0 0 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0
## 14 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0
## 15 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0
## 16 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0
## 17 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0
## 18 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0
## 19 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0
## 20 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1
## 21 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1
## 22 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0
## 23 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0
## 24 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## 25 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0On modifie la formule pour spécifier une pénalité de type
“graph-guided fused lasso” pour la variable area.
formule_2 <- rent ~
p(size, pen = "flasso") +
p(rooms, pen = "flasso") +
p(year, pen = "flasso") +
p(area, pen = "ggflasso") +
p(warm, pen = "flasso") +
p(central, pen = "flasso") +
p(tiles, pen = "flasso") +
p(bathextra, pen = "flasso") +
p(kitchen, pen = "flasso") +
p(standing, pen = "flasso")On procède maintenant à l’ajustement, en spécifiant la matrice
d’adjacence utilisée pour area.
ajustement_lasso_2 <- glmsmurf(formula = formule,
family = gaussian(),
data = donnees,
pen.weights = "glm.stand",
lambda = "is.aic",
control = list(lambda.reest = TRUE),
adj.matrix = list(area = adj_mat_Munich))summary(ajustement_lasso_2)##
## Call: glmsmurf(formula = formule, family = gaussian(), data = donnees,
## lambda = "is.aic", pen.weights = "glm.stand", adj.matrix = list(area = adj_mat_Munich),
## control = list(lambda.reest = TRUE))
##
## Deviance residuals of estimated model:
## Min. 1st Qu. Median 3rd Qu. Max.
## -6.73645 -1.26900 -0.00965 1.25494 8.90310
##
## Deviance residuals of re-estimated model:
## Min. 1st Qu. Median 3rd Qu. Max.
## -6.5682 -1.2127 -0.0181 1.2592 8.5690
##
##
## Coefficients:
## Estimated Re-estimated
## Intercept 10.6600 11.0720
## size(39,50] -1.6732 -1.5464
## size(50,55] -2.0172 -1.9632
## size(55,62] -2.0172 -1.9632
## size(62,67] -2.2777 -2.3030
## size(67,73] -2.2777 -2.3030
## size(73,80] -2.2777 -2.3030
## size(80,87] -2.2777 -2.3030
## size(87,100] -2.3073 -2.4928
## size(100,185] -2.3616 -2.7943
## rooms2 * *
## rooms3 * *
## rooms4 -0.2537 -0.3278
## rooms5 -0.2537 -0.3278
## rooms6 -0.2537 -0.3278
## year1920 -0.9019 -1.0810
## year1930 -0.9019 -1.0810
## year1940 -0.9019 -1.0810
## year1950 -0.1901 -0.3711
## year1960 0.1003 0.0728
## year1970 0.1317 0.3469
## year1980 1.2073 1.1245
## year1990 1.5781 1.6513
## year2000 1.5781 1.6513
## area2 -0.3013 -0.4287
## area3 * *
## area4 -0.3013 -0.4287
## area5 -0.3013 -0.4287
## area6 -0.6784 -1.1076
## area7 -0.9050 -1.6387
## area8 -0.6784 -1.1076
## area9 -0.4107 -0.7080
## area10 -0.6784 -1.1076
## area11 -0.9050 -1.6387
## area12 -0.3013 -0.4287
## area13 -0.4107 -0.7080
## area14 -0.9050 -1.6387
## area15 -0.6784 -1.1076
## area16 -0.9050 -1.6387
## area17 -0.6784 -1.1076
## area18 -0.3013 -0.4287
## area19 -0.6784 -1.1076
## area20 -0.6784 -1.1076
## area21 -0.6784 -1.1076
## area22 -0.9050 -1.6387
## area23 -0.9050 -1.6387
## area24 -0.9050 -1.6387
## area25 -0.6784 -1.1076
## warm1 -1.9562 -2.0509
## central1 -1.2423 -1.3296
## tiles1 -0.4160 -0.5433
## bathextra1 0.0196 0.4945
## kitchen1 1.0584 1.1853
## standing1 0.4266 0.3762
## standing2 1.1314 1.4360
## ---
## '*' indicates a zero coefficient
##
##
## Null deviance: 12486.047 on 2052 degrees of freedom
## -------------------
## Estimated model:
##
## Residual deviance: 8026.104 on 2029 degrees of freedom
## AIC: 8675.216; BIC: 8810.265; GCV score: 4.002
## Penalized minus log-likelihood: 2.143388
## -------------------
## Re-estimated model:
##
## Residual deviance: 7839.303 on 2029 degrees of freedom
## AIC: 8626.869; BIC: 8761.919; GCV score: 3.909
## Objective function: 2.188856
## -------------------
## lambda: 0.02072845; lambda1: 0; lambda2: 0
## Number of iterations: 199
## Final step size: 0.8019531
## Convergence: Succesful convergence3 Fusion par ajustement itératif d’arbres
On reprend les données de l’exemple précédent, en transformant les variables ordinales en facteurs ordonnés, de manière à ce qu’elles soient reconnues comme telles.
On travaille de nouveau avec le jeu de données des loyers à Munich.
load(url("https://irma.math.unistra.fr/~jberard/rent_Munich")) # Données retraitées de catdata
donnees$size<-cut(donnees$size, breaks=quantile(donnees$size, seq(0,1,length.out=11)),include.lowest=TRUE)
donnees$size<-factor(donnees$size, ordered=TRUE)
donnees$year<-factor(donnees$year, ordered=TRUE)
donnees$rooms<-factor(donnees$rooms, ordered=TRUE)
donnees$standing<-factor(donnees$standing, ordered=TRUE)On charge d’abord le paquet structree.
library(structree)## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmList## Loading required package: penalized## Loading required package: survival## Welcome to penalized. For extended examples, see vignette("penalized").On spécifie ensuite le fait que les variables possédant plus d’une
modalité seront traitées de manière arborescente (avec les termes
tr).
formule <- rent ~
tr(size) +
tr(rooms) +
tr(year) +
tr(area) +
warm +
central +
tiles +
bathextra +
kitchen +
tr(standing)On effectue ensuite l’ajustement. On spécifie le fait que le critère d’arrêt du partitionnement est donné par la valeur du critère AIC du modèle GLM estimé.
ajustement_arbre <- structree(formula = formule, family = gaussian(), data = donnees, stop_criterion = "AIC")## Split 1
## Split 2
## Split 3
## Split 4
## Split 5
## Split 6
## Split 7
## Split 8
## Split 9
## Split 10
## Split 11
## Split 12
## Split 13
## Split 14
## Split 15
## Split 16
## Split 17
## Split 18
## Split 19
## Split 20
## Split 21
## Split 22
## Split 23
## Split 24
## Split 25
## Split 26
## Split 27
## Split 28
## Split 29
## Split 30
## Split 31
## Split 32
## Split 33
## Split 34
## Split 35
## Split 36
## Split 37
## Split 38
## Split 39
## Split 40
## Split 41
## Split 42
## Split 43
## Split 44
## Split 45
## Split 46
## Split 47
## Split 48
## Split 49Le tracé ci-après permet de visualiser les arbres associés au partitionnement de chaque variable.
plot(ajustement_arbre, select = 1)
plot(ajustement_arbre, select = 2)
plot(ajustement_arbre, select = 3)
plot(ajustement_arbre, select = 4)
plot(ajustement_arbre, select = 5)
Le tracé ci-après permet de visualiser les groupements retenus par la méthode.
plot(ajustement_arbre, select = 1, result = TRUE)
plot(ajustement_arbre, select = 2, result = TRUE)
plot(ajustement_arbre, select = 3, result = TRUE)
plot(ajustement_arbre, select = 4, result = TRUE)
plot(ajustement_arbre, select = 5, result = TRUE)
On peut également visualiser l’évolution des coefficients estimés au fur et à mesure du partitionnement.
plot(ajustement_arbre, select = 1, paths = TRUE)
plot(ajustement_arbre, select = 2, paths = TRUE)
plot(ajustement_arbre, select = 3, paths = TRUE)
plot(ajustement_arbre, select = 4, paths = TRUE)
plot(ajustement_arbre, select = 5, paths = TRUE)
Et encore la valeur du critère retenu (ici AIC) au cours des étapes successives de partitionnement.
plot(ajustement_arbre[["tune_values"]], ylab = "AIC du GLM", xlab = "nombre de bifurcations")
On examine rapidement Le modèle GLM ajusté in fine. On constate la réduction très nette de la valeur du critère AIC par rapport au modèle non-fusionné, avec une valeur de \(8626\) environ, contre \(8682.8\) pour le modèle non-fusionné.
summary(ajustement_arbre[["model"]])##
## Call:
## glm(formula = y ~ x6 + x7 + x8 + x9 + x10 + s10 + s36 + s45 +
## s22 + s50 + s34 + s13 + s422 + s30 + s33 + s51 + s41 + s11 +
## s37 + s18 + s412 + s35, family = family, data = dat, weights = we,
## offset = off)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4287 0.2149 43.871 < 2e-16 ***
## x61 -2.0502 0.2780 -7.374 2.39e-13 ***
## x71 -1.3328 0.1916 -6.957 4.67e-12 ***
## x81 -0.5439 0.1154 -4.715 2.58e-06 ***
## x91 0.4828 0.1583 3.050 0.002320 **
## x101 1.1747 0.1745 6.731 2.19e-11 ***
## s10 -1.5445 0.1885 -8.193 4.44e-16 ***
## s36 0.7725 0.1896 4.075 4.78e-05 ***
## s45 0.3990 0.1469 2.717 0.006640 **
## s22 -0.3979 0.1496 -2.660 0.007881 **
## s50 0.3758 0.1028 3.655 0.000264 ***
## s34 0.4383 0.1522 2.880 0.004020 **
## s13 -0.3648 0.1198 -3.046 0.002350 **
## s422 0.5245 0.1303 4.024 5.93e-05 ***
## s30 -1.0706 0.1786 -5.993 2.43e-09 ***
## s33 0.7192 0.1875 3.835 0.000129 ***
## s51 1.0437 0.3072 3.398 0.000693 ***
## s41 0.4392 0.1740 2.524 0.011672 *
## s11 -0.4158 0.1614 -2.576 0.010075 *
## s37 0.5253 0.2007 2.617 0.008926 **
## s18 -0.4027 0.1833 -2.197 0.028117 *
## s412 0.2728 0.1460 1.869 0.061810 .
## s35 0.2721 0.1480 1.838 0.066168 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 3.864108)
##
## Null deviance: 12486.0 on 2052 degrees of freedom
## Residual deviance: 7844.1 on 2030 degrees of freedom
## AIC: 8626.1
##
## Number of Fisher Scoring iterations: 2