R data.table : Différence entre les résultats des régressions imbriquées

Question

Je compare deux stratégies alternatives pour estimer les modèles de régression linéaire sur des sous-ensembles de données en utilisant le paquet data.table pour R. Les deux stratégies produisent les mêmes coefficients, donc elles semblent équivalentes. Cette apparence est trompeuse. Ma question est :

Pourquoi les données stockées à l'intérieur des modèles lm sont-elles différentes ?

library(data.table)

dat = data.table(mtcars)

# stratégie 1
mod1 = dat[, .(models = .(lm(hp ~ mpg, data = .SD))), by = vs]

# stratégie 2
mod2 = dat[, .(data = .(.SD)), by = vs][
           , models := lapply(data, function(x) lm(hp ~ mpg, x))]

À première vue, les deux approches semblent produire des résultats identiques :

# stratégie 1
coef(mod1$models[[1]])
#> (Intercept)         mpg 
#>   357.97866   -10.12576

# stratégie 2
coef(mod2$models[[1]])
#> (Intercept)         mpg 
#>   357.97866   -10.12576

Cependant, si j'essaie d'extraire des données du (cadre de modèle étendu), j'obtiens des résultats différents :

# stratégie 1
expanded_frame1 = expand.model.frame(mod1$models[[1]], "am")
table(expanded_frame1$am)
#> 
#>  0  1 
#>  7 11

# stratégie 2
expanded_frame2 = expand.model.frame(mod2$models[[1]], "am")
table(expanded_frame2$am)
#> 
#>  0  1 
#> 12  6

Il s'agit d'un exemple trivial de travail minimal. Mon cas d'utilisation réel est que j'ai obtenu des résultats radicalement différents lors de l'application de sandwich::vcovCL pour calculer des erreurs standard clusterisées pour mes modèles.

Éditer :

Je valide la réponse de @TimTeaFan (excellent travail de détective !) mais j'ajoute un peu d'information utile ici pour les lecteurs futurs.

Comme l'a souligné @achim-zeileis ailleurs, nous pouvons reproduire un comportement similaire dans l'environnement global :

d <- subset(mtcars, subst = vs == 0)
m0 <- lm(hp ~ mpg, data = d)
d <- mtcars[0, ]
expand.model.frame(m0, "am")

[1] hp  mpg am
<0 rows> (or 0-length row.names)

Cela ne semble pas être un problème spécifique à data.table. En général, nous devons être prudents lors de la réévaluation des données à partir d'un modèle.

Answer 1

Je n'ai pas de réponse complète, mais j'ai pu localiser le problème dans une certaine mesure.

Lorsque nous comparons les résultats des deux modèles, nous pouvons constater que le résultat est identique, à l'exception des appels, qui sont différents (ce qui est logique, puisqu'ils sont effectivement différents) :

# compare models
purrr::map2(mod1$models[[1]], mod2$models[[1]], all.equal)
#> $coefficients
#> [1] TRUE
#> 
#> $residuals
#> [1] TRUE
#> 
#> $effects
#> [1] TRUE
#> 
#> $rank
#> [1] TRUE
#> 
#> $fitted.values
#> [1] TRUE
#> 
#> $assign
#> [1] TRUE
#> 
#> $qr
#> [1] TRUE
#> 
#> $df.residual
#> [1] TRUE
#> 
#> $xlevels
#> [1] TRUE
#> 
#> $call
#> [1] "target, current do not match when deparsed"
#> 
#> $terms
#> [1] TRUE
#> 
#> $model
#> [1] TRUE

Il semble donc que l'appel initial fonctionne correctement avec les deux approches, le problème se pose dès que nous essayons d'accéder aux données sous-jacentes.

Si nous regardons comment expand.model.frame reçoit ses données, on peut voir qu'il appelle eval(model$call$data, envir) donde envir est défini comme suit environment(formula(model)) l'environnement associé à la formule du lm objet.

Si nous examinons les données de l'environnement associé à chaque modèle et que nous les comparons avec les données que nous attendons, nous pouvons constater que la deuxième approche produit les données que nous attendons, tandis que la première approche utilisant des .SD dans l'appel donne des données différentes.

Je ne sais toujours pas pourquoi ni ce qui se passe, mais nous savons maintenant que le problème se situe au niveau de l'appel à l'aide. .SD . J'ai d'abord pensé que cela pouvait être causé par le fait de nommer un data.table .SD mais après avoir joué avec des modèles où les données sont un data.table appelé .SD cela ne semble pas être le problème.

# data of model 2 (identical to subsetted mtcars)
environment(formula(mod2$models[[1]]))$x[order(mpg),]
#>      mpg cyl  disp  hp drat    wt  qsec am gear carb
#>  1: 10.4   8 472.0 205 2.93 5.250 17.98  0    3    4
#>  2: 10.4   8 460.0 215 3.00 5.424 17.82  0    3    4
#>  3: 13.3   8 350.0 245 3.73 3.840 15.41  0    3    4
#>  4: 14.3   8 360.0 245 3.21 3.570 15.84  0    3    4
#>  5: 14.7   8 440.0 230 3.23 5.345 17.42  0    3    4
#>  6: 15.0   8 301.0 335 3.54 3.570 14.60  1    5    8
#>  7: 15.2   8 275.8 180 3.07 3.780 18.00  0    3    3
#>  8: 15.2   8 304.0 150 3.15 3.435 17.30  0    3    2
#>  9: 15.5   8 318.0 150 2.76 3.520 16.87  0    3    2
#> 10: 15.8   8 351.0 264 4.22 3.170 14.50  1    5    4
#> 11: 16.4   8 275.8 180 3.07 4.070 17.40  0    3    3
#> 12: 17.3   8 275.8 180 3.07 3.730 17.60  0    3    3
#> 13: 18.7   8 360.0 175 3.15 3.440 17.02  0    3    2
#> 14: 19.2   8 400.0 175 3.08 3.845 17.05  0    3    2
#> 15: 19.7   6 145.0 175 3.62 2.770 15.50  1    5    6
#> 16: 21.0   6 160.0 110 3.90 2.620 16.46  1    4    4
#> 17: 21.0   6 160.0 110 3.90 2.875 17.02  1    4    4
#> 18: 26.0   4 120.3  91 4.43 2.140 16.70  1    5    2

# subset and order mtcars data
mtcars_vs0 <- subset(mtcars, vs == 0)
mtcars_vs0[order(mtcars_vs0$mpg), ]
#>                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
#> Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
#> Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
#> Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
#> Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
#> Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
#> Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
#> Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
#> AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
#> Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
#> Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
#> Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
#> Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
#> Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
#> Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
#> Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
#> Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
#> Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
#> Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2

# data of model 1 (not identical to mtcars)
environment(formula(mod1$models[[1]]))$.SD[order(mpg),]
#>      mpg cyl  disp  hp drat    wt  qsec am gear carb
#>  1: 15.0   8 301.0 335 3.54 3.570 14.60  1    5    8
#>  2: 15.8   8 351.0 264 4.22 3.170 14.50  1    5    4
#>  3: 17.8   6 167.6 123 3.92 3.440 18.90  0    4    4
#>  4: 18.1   6 225.0 105 2.76 3.460 20.22  0    3    1
#>  5: 19.2   6 167.6 123 3.92 3.440 18.30  0    4    4
#>  6: 19.7   6 145.0 175 3.62 2.770 15.50  1    5    6
#>  7: 21.4   6 258.0 110 3.08 3.215 19.44  0    3    1
#>  8: 21.4   4 121.0 109 4.11 2.780 18.60  1    4    2
#>  9: 21.5   4 120.1  97 3.70 2.465 20.01  0    3    1
#> 10: 22.8   4 108.0  93 3.85 2.320 18.61  1    4    1
#> 11: 22.8   4 140.8  95 3.92 3.150 22.90  0    4    2
#> 12: 24.4   4 146.7  62 3.69 3.190 20.00  0    4    2
#> 13: 26.0   4 120.3  91 4.43 2.140 16.70  1    5    2
#> 14: 27.3   4  79.0  66 4.08 1.935 18.90  1    4    1
#> 15: 30.4   4  75.7  52 4.93 1.615 18.52  1    4    2
#> 16: 30.4   4  95.1 113 3.77 1.513 16.90  1    5    2
#> 17: 32.4   4  78.7  66 4.08 2.200 19.47  1    4    1
#> 18: 33.9   4  71.1  65 4.22 1.835 19.90  1    4    1

Ajouter

J'ai essayé de creuser un peu plus pour voir ce qui se passe. J'ai d'abord appelé debug(as.formula) et a ensuite examiné les objets suivants à chaque itération :

object
ls(environment(object))

Nous pouvons voir que dans la "stratégie 2" chaque formule est associée à un environnement différent, et en regardant l'environnement nous voyons qu'il contient un objet x qui, lorsqu'il est inspecté ( environment(object)$x ) contient la valeur attendue mtcars données.

Dans la "stratégie 1" cependant, nous pouvons observer que chaque appel à as.formula associe le même environnement avec la formule en cours de création. De plus, en inspectant l'environnement, nous pouvons voir qu'il est peuplé des vecteurs uniques du sous-ensemble mtcars données (par exemple am , carb , cyl etc.) ainsi que certaines fonctions (par ex. .POSIXt , Cfastmean , strptime etc.). C'est probablement là que les choses se gâtent. Je soupçonne qu'en associant le même environnement à deux formules (modèles) différentes, les données sous-jacentes du premier modèle sont "mises à jour" lorsque le second modèle est calculé. Cela devrait également être la raison pour laquelle la sortie du modèle lui-même est correcte. Au moment où le premier modèle est calculé, les données sont encore correctes. Elles sont écrasées par le deuxième modèle, qui est donc lui aussi correct. Mais lorsqu'on accède ensuite aux données sous-jacentes, les choses se gâtent.

Note complémentaire
J'étais curieux de savoir si nous pouvions observer des problèmes et des différences similaires dans le tidyverse en utilisant expand.model.frame et la réponse est "oui". Ici, le nouveau rowwise entraîne une erreur, tandis que la notation group_map ainsi que le map travail d'approche :

# dplyr approaches:

# group_map: works
mod3 <- mtcars %>%
  group_by(vs) %>%
  group_map(~ lm(hp ~ mpg, data = .x))

expand.model.frame(mod3[[1]], "am")

# mutate / rowwise: does not work
mod4 <- mtcars %>%
  nest_by(vs) %>%
  mutate(models = list(lm(hp ~ mpg, data = data)))

expand.model.frame(mod4$models[[1]], "am")

# mutate / map: works
mod5 <- mtcars %>%
  tidyr::nest(data = !vs) %>%
  mutate(models = purrr::map(data, ~ lm(hp ~ mpg, data = .x)))

expand.model.frame(mod5$models[[1]], "am")