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