Sélection de modèle avec MuMIn et model.select

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Gilles San Martin
Messages : 211
Enregistré le : 08 Juin 2007, 17:25

Sélection de modèle avec MuMIn et model.select

Messagepar Gilles San Martin » 08 Fév 2011, 03:29

Bonjour

J'ai fait quelques comparaisons des sorties des fonctions du package MuMIn et de la fonction model.select disponible ici : viewtopic.php?t=3516

Globalement les résultats sont identiques ou similaires (pour les tableaux d'aic) mais il y a quand même l'une ou l'autre divergence de résultats (pour les "model averaged coefficients"). Ces quelques exemples peuvent aussi montrer quelques pièges dans lesquels on pourrait vite tomber...


Exemple 1 : "Cement data" de Burnham & Anderson 2002

C'est un jeu de données utilisé par ces auteurs et il permet donc de comparer les valeurs obtenues avec des valeurs externes.
Pour pouvoir comparer il faut créer les 16 modèles possibles, ne pas inclure de modèle avec juste l'intercept et utiliser la méthode "classique" pour le model averaging (par opposition à l'approche par shrinkage model averaging dans laquelle un paramètre est mis à 0 quand il est absent d'un modèle). Cette approche classique n'est pas la méthode par défaut dans MuMIn. Model.select donne les deux valeurs, av.coef2 correspondant à l'approche classique. Il semble que l'autre approche (shrinkage) soit de plus en plus privilégiée dans les cas courants malgré ce qui est présenté dans Burnham & Anderson (il me semble que cette approche ait effectivement plus d'avantages à part pour les erreurs standard)

On charge la libraire, le code de model.select, le jeu de données et on construit le modèle complet

Code : Tout sélectionner

library(MuMIn)
source("/home/gilles/stats/R/WorkingDirectory/rprojects/model.select_0.3.R")

data(Cement)
lm1 <- lm(y ~ X1 + X2 + X3 + X4, data = Cement)


Avec MuMIn, il y a une fonction pour estimer les différents modèles et une fonction pour faire le model averaging.

Code : Tout sélectionner

mumindd <- dredge(lm1)
models.list <- get.models(mumindd,subset = delta <= 10000) # extraire tous les modèles
models.list <- models.list[-(which(names(models.list)=="1"))] # éliminer le modèle avec seulement l'intercept
muminavg <- model.avg(models.list,method="NA") # obtenir les "model averaged coefficients"


model.select donne tout, on se débrouille ensuite...

Code : Tout sélectionner

res <- model.select(lm1,null.model=FALSE, dec=4)


Pour les tableaux d'aic, les résultats sont strictement identiques :
mumin :

Code : Tout sélectionner

> mumindd
Global model: lm(formula = y ~ X1 + X2 + X3 + X4, data = Cement)
---
Model selection table
    (Int)    X1      X2      X3      X4 k   R.sq Adj.R.sq     RSS    AIC   AICc  delta weight
6   52.58 1.468  0.6623                 4 0.9787   0.9744   57.90  64.31  69.31  0.000  0.566
13  71.65 1.452  0.4161         -0.2365 5 0.9823   0.9764   47.97  63.87  72.44  3.125  0.119
12  48.19 1.696  0.6569  0.2500         5 0.9823   0.9764   48.11  63.90  72.48  3.163  0.116
8  103.10 1.440                 -0.6140 4 0.9725   0.9670   74.76  67.63  72.63  3.322  0.107
14 111.70 1.052         -0.4100 -0.6428 5 0.9813   0.9750   50.84  64.62  73.19  3.879  0.081
15 203.60       -0.9234 -1.4480 -1.5570 5 0.9728   0.9638   73.81  69.47  78.04  8.727  0.007
16  62.41 1.551  0.5102  0.1019 -0.1441 6 0.9824   0.9736   47.86  65.84  79.84 10.520  0.003
11 131.30               -1.2000 -0.7246 4 0.9353   0.9223  175.70  78.74  83.74 14.430  0.000
9   72.07        0.7313 -1.0080         4 0.8470   0.8164  415.40  89.93  94.93 25.620  0.000
5  117.60                       -0.7382 3 0.6745   0.6450  883.90  97.74 100.40 31.100  0.000
3   57.42        0.7891                 3 0.6663   0.6359  906.30  98.07 100.70 31.420  0.000
10  94.16        0.3109         -0.4569 4 0.6801   0.6161  868.90  99.52 104.50 35.210  0.000
2   81.48 1.869                         3 0.5339   0.4916 1266.00 102.40 105.10 35.770  0.000
7   72.35 2.312          0.4945         4 0.5482   0.4578 1227.00 104.00 109.00 39.700  0.000
4  110.20               -1.2560         3 0.2859   0.2210 1939.00 108.00 110.60 41.310  0.000
1   95.42                               2 0.0000   0.0000 2716.00 110.30 111.50 42.220  0.000


model.select :

Code : Tout sélectionner

 > res$AICtab[,-c(6,10,11)]
            model id k  n    loglik     AICc AICc.delta AICc.w sum.w
3          X1+ X2  3 4 13 -28.15620  69.3124     0.0000 0.5657 0.566
11     X1+ X2+ X4 11 5 13 -26.93314  72.4377     3.1253 0.1186 0.684
7      X1+ X2+ X3  7 5 13 -26.95180  72.4750     3.1626 0.1164 0.801
9          X1+ X4  9 4 13 -29.81705  72.6341     3.3217 0.1075 0.908
13     X1+ X3+ X4 13 5 13 -27.30998  73.1914     3.8790 0.0813 0.989
14     X2+ X3+ X4 14 5 13 -29.73414  78.0397     8.7273 0.0072 0.997
15 X1+ X2+ X3+ X4 15 6 13 -26.91834  79.8367    10.5243 0.0029 1.000
12         X3+ X4 12 4 13 -35.37249  83.7450    14.4326 0.0004 1.000
6          X2+ X3  6 4 13 -40.96477  94.9295    25.6171 0.0000 1.000
8              X4  8 3 13 -45.87202 100.4107    31.0983 0.0000 1.000
2              X2  2 3 13 -46.03520 100.7371    31.4247 0.0000 1.000
10         X2+ X4 10 4 13 -45.76086 104.5217    35.2093 0.0000 1.000
1              X1  1 3 13 -48.20594 105.0785    35.7661 0.0000 1.000
5          X1+ X3  5 4 13 -48.00454 109.0091    39.6967 0.0000 1.000
4              X3  4 3 13 -50.97990 110.6265    41.3141 0.0000 1.000


Les résultats du model averaging sont pratiquement identiques pour les coefficients et légèrement différents pour les erreurs standard.
Ici çà ne porte clairement pas à conséquence mais je me demande quel peut être l'impact de telles différences dans d'autres cas (par exemple quand on est sur d'autres échalles : logit, log,...).
Normalement on devrait obtenir strictement la même chose. Peut-être est-ce un problème d'arrondi ou de nombre de chiffres après la virgule dans le cours des calculs ?

A titre de comparaison, Burnham & Anderson donnent les valeurs suivantes (SAS, méthode Least Squares): (p 180)
x1 : 1.4561 - se = 1.1755
x2 : 0.6110 - se = 0.1206

Il faut comparer Coefficient et Unconditional SE dans MuMIn avec av.coef2 et av.se2 dans model.select

mumin :

Code : Tout sélectionner

 Averaged model parameters:
            Coefficient Variance     SE Unconditional SE Lower CI Upper CI
(Intercept)     65.7000 4.18e+05 20.500           20.800   24.900 106.0000
X1               1.4600 1.85e-03  0.176            0.194    1.080   1.8400
X2               0.6110 1.44e-03  0.113            0.120    0.375   0.8470
X3              -0.0715 4.60e-02  0.422            0.438   -0.931   0.7880
X4              -0.4980 6.48e-03  0.231            0.239   -0.966  -0.0292


model.select :

Code : Tout sélectionner

> res$mod.av
             freq priorw      w av.coef   av.se av.coef2  av.se2
(Intercept) 1.000 1.0000 1.0000 65.7149 20.4701  65.7149 20.4701
X1          0.533 0.3681 0.9924  1.4450  0.1753   1.4561  0.1756
X2          0.533 0.3681 0.8108  0.4953  0.1539   0.6109  0.1130
X3          0.533 0.3681 0.2083 -0.0149  0.0867  -0.0715  0.4217
X4          0.533 0.3681 0.3179 -0.1582  0.1239  -0.4978  0.2305

Gilles San Martin
Messages : 211
Enregistré le : 08 Juin 2007, 17:25

Example 2

Messagepar Gilles San Martin » 08 Fév 2011, 03:33

Exemple 2 : Modèle de Poisson surdispersé avec QAICc

En résumé : les tables d'aic sont ici aussi identiques (ou en tous cas cohérentes) par contre les résultats du "model averaging" peuvent être très différents (en particulier pour les erreurs standard).
Il y a deux manières principales de faire dans mumin et une seule dans model.select, ce qui donne 3 cas de figure :

mumin 1 : on utilise un modèle estimé par quasi vraisemblance et mumin utilise automatiquement le coefficient de surdispersion estimé par le modèle ainsi que les erreurs standard "correctes" (corrigées pour la surdispersion). Par contre c'est un peu curieux d'estimer un AIC sur un modèle qui n'a pas de vraisemblance (quasi-likelihood). J'imagine qu'il utilise la vraisemblance du modèle équivalent estimé par maximum likelihood...

mumin 2 : on utilise un modèle classique (maximum de vraisemblance) et on spécifie la valeur du coefficient de surdispersion que l'on veut utiliser ("chat" - pour "c chapeau"). Ici j'ai utilisé un coefficient de surdispersion estimé sur base de la deviance divisée par les degrés de liberté résiduels comme proposé par B&A et parce que c'est cette valeur qui est utilisée dans model.select. Par contre il semble que les erreurs standard des modèles ne soient pas corrigées pour la surdisperssion ce qui donne des valeurs finales sans doute beaucoup trob basses et vraisemblablement assez loin de la réalité (bien que cette notion soit très relative...).

Model.select ne laisse pas le choix : le coefficient de surdispersion est estimé sur base de la déviance du modèle complet. Les erreurs standard sont corrigées (multipliées par la racine carrée du coefficient de surdispersion) avant le model averaging.
Donc mumin2 et model.select utilisent ici le même coefficient de surdispersion mais pas mumin1


Voici l'exemple :

Code : Tout sélectionner

 library(MASS)
data(quine)
model <- glm(Days~(Eth + Sex + Age), family=poisson, data=quine)


La surdispersion est de :

Code : Tout sélectionner

> overdisp(model)
pearsonresid     deviance
    13.88966     12.44646

> sqrt(overdisp(model))
pearsonresid     deviance
    3.726884     3.527954

La surdispersion estimée sur base des résidus de Pearson est en général très proche de la surdispersion estimée par quasi-vraisemblance (cfr mumin1)



mumin1 :

Code : Tout sélectionner

# mumin with quasilikelihood model
model <- glm(Days~(Eth + Sex + Age), family=quasipoisson, data=quine)

mumindd.quasi <- dredge(model)
models.list <- get.models(mumindd.quasi,subset = delta <= 10000)
muminavg.quasi <- model.avg(models.list,method="NA")


mumin2 :

Code : Tout sélectionner

# mumin with poisson model and overdispersion parameter based on residual deviance
model <- glm(Days~(Eth + Sex + Age), family=poisson, data=quine)

mumindd <- dredge(model, rank="QAICc", chat=overdisp(model)["deviance"])
models.list <- get.models(mumindd,subset = delta <= 10000)
muminavg <- model.avg(models.list,method="NA")


model.select :

Code : Tout sélectionner

# model select
model <- glm(Days~(Eth + Sex + Age), family=poisson, data=quine)
res <- model.select(model, srt="QAICc")



Les tables AIC qui en sortent :
mumin1 est légèrement différent de mumin2 et model.select ce qui est attendu puisqu'on utilise pas le même coefficient de surdispersion.
Curieusement les QAICc ne sont pas les mêmes entre mumin2 et model.select. Mais les delta et les poids (weights) sont strictement identiques et c'est le principal. Je suppose qu'une constante a du se perdre en route.

mumin1 :

Code : Tout sélectionner

> mumindd.quasi
Global model: glm(formula = Days ~ (Eth + Sex + Age), family = quasipoisson,
    data = quine)
---
Model selection table
  (Int) Age Eth Sex k Dev. QAICc delta  weight
5 2.939 +   +       5 1749 138.5  0.000 0.600
8 2.871 +   +   +   6 1743 140.3  1.749 0.250
3 3.056     +       2 1892 142.4  3.868 0.087
7 2.977     +   +   3 1875 143.3  4.767 0.055
2 2.698 +           4 1915 148.3  9.752 0.005
6 2.631 +       +   5 1908 150.0 11.480 0.002
1 2.801             1 2074 153.4 14.850 0.000
4 2.723         +   2 2057 154.3 15.770 0.000


mumin1 :

Code : Tout sélectionner

> mumindd
Global model: glm(formula = Days ~ (Eth + Sex + Age), family = poisson, data = quine)
---
Model selection table
  (Int) Age Eth Sex k Dev. QAICc delta  weight
5 2.939 +   +       5 1749 153.1  0.000 0.638
8 2.871 +   +   +   6 1743 154.8  1.696 0.273
3 3.056     +       2 1892 158.2  5.062 0.051
7 2.977     +   +   3 1875 158.9  5.820 0.035
2 2.698 +           4 1915 164.3 11.140 0.002
6 2.631 +       +   5 1908 165.9 12.810 0.001
1 2.801             1 2074 170.7 17.560 0.000
4 2.723         +   2 2057 171.5 18.340 0.000


model.select :

Code : Tout sélectionner

> res$AICtab[,c("model", "k","n", "loglik", "QAICc", "QAICc.delta", "QAICc.w", "sum.w")]
          model k   n    loglik   QAICc QAICc.delta QAICc.w sum.w
6      Eth+ Age 5 146 -1168.674 200.396       0.000   0.638 0.638
8 Eth+ Sex+ Age 6 146 -1165.491 202.092       1.696   0.273 0.911
2           Eth 2 146 -1240.226 205.459       5.062   0.051 0.962
4      Eth+ Sex 3 146 -1231.784 206.217       5.820   0.035 0.996
5           Age 4 146 -1251.513 211.532      11.136   0.002 0.999
7      Sex+ Age 5 146 -1248.376 213.204      12.807   0.001 1.000
1               1 146 -1331.005 217.961      17.564   0.000 1.000
3           Sex 2 146 -1322.854 218.736      18.340   0.000 1.000



Les résultats du model averaging.
Il faut comparer Coefficient et Unconditional SE dans MuMIn avec av.coef2 et av.se2 dans model.select
On voit bien que les valeurs d'erreur standard pour mumin2 sont beaucoup plus faible et elles semblent même presque identiques aux valeurs du modèle complet, comme si il n'y avait pas eu de moyenne de modèles...
Il faut prendre garde aussi que si on utilise les familles de type quasi- avec lmer, on obtiendra généralement des résultats faux. Model.select n'accepte pas de traiter des modèles lmer avec QAIC mais ce serait assez facile à changer.

mumin1

Code : Tout sélectionner

> muminavg.quasi

Averaged model parameters:
            Coefficient Variance    SE Unconditional SE Lower CI Upper CI
AgeF1            -0.251 0.003930 0.250            0.252   -0.746    0.243
AgeF2             0.342 0.002610 0.226            0.228   -0.105    0.789
AgeF3             0.282 0.003170 0.237            0.239   -0.187    0.751
EthN             -0.535 0.000604 0.157            0.158   -0.845   -0.225
(Intercept)       2.930 0.001590 0.198            0.200    2.540    3.320
SexM              0.117 0.000625 0.158            0.159   -0.196    0.429


mumin2

Code : Tout sélectionner

> muminavg

Averaged model parameters:
            Coefficient Variance     SE Unconditional SE Lower CI Upper CI
AgeF1            -0.251 2.28e-05 0.0690           0.0696   -0.388   -0.115
AgeF2             0.342 1.37e-05 0.0609           0.0614    0.222    0.462
AgeF3             0.282 1.72e-05 0.0644           0.0650    0.155    0.410
EthN             -0.534 3.24e-06 0.0424           0.0427   -0.617   -0.450
(Intercept)       2.930 2.38e-05 0.0662           0.0665    2.800    3.060
SexM              0.113 4.54e-06 0.0454           0.0458    0.023    0.202


model.select

Code : Tout sélectionner

> res$mod.av
            freq priorw     w av.coef av.se av.coef2 av.se2
(Intercept)  1.0  1.000 1.000   2.927 0.191    2.927  0.191
EthN         0.5  0.803 0.996  -0.532 0.147   -0.534  0.148
SexM         0.5  0.803 0.309   0.035 0.052    0.113  0.149
AgeF1        0.5  0.986 0.914  -0.230 0.218   -0.251  0.237
AgeF2        0.5  0.986 0.914   0.313 0.198    0.342  0.214
AgeF3        0.5  0.986 0.914   0.258 0.207    0.282  0.225

Gilles San Martin
Messages : 211
Enregistré le : 08 Juin 2007, 17:25

Example 3

Messagepar Gilles San Martin » 08 Fév 2011, 03:36

Exemple 3 : comparaison des performances

model.select est plus rapide. Mais si on reste raisonnable avec le nombre de modèles, çà reste sans importance.

Model.select utilise aussi moins de ressources = ~30Mo de RAM dans cet exemple contre ~500Mo pour mumin (sur un système 64bits). Model select ne garde par défaut que le strict nécessaire et nettoie bien derrière lui (développé au départ sur une petite configuration Windows)


Exemple de temps de calcul avec un jeu de données simulé :

Code : Tout sélectionner

## create a dataset with 20 explanatory variables

set.seed(123)

X <- matrix(nrow=200, ncol=20, data= rnorm(4000))
beta <- as.numeric(seq(5,0,length.out=20))
colnames(X) <- names(beta) <- paste("a",1:20, sep="")
Y <- X %*% beta + rnorm(200,0,5)

data <- as.data.frame(cbind (Y, X))

beta
model <- lm ( Y ~  a9 + a10 + a11 + a12 + a13 + a14 + a15 + a16 + a17 + a18 + a19 + a20, data=data)
summary(model)


Estimation du temps de calcul pour 2^12 modèles (ce n'est pas raisonnable en général...)

model.select :

Code : Tout sélectionner

> (time.model.select <- system.time(res <- model.select(model)))
utilisateur     système      écoulé
     38.570       0.000      38.453


mumin :

Code : Tout sélectionner

mumin <- function () {
   mumindd <- dredge(model)
   models.list <- get.models(mumindd,subset = delta <= 10000)
   muminavg <- model.avg(models.list,method="0")
   }
> (time.mumin <- system.time(mumin()))
utilisateur     système      écoulé
    115.530       0.540     115.585


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité