r - R data.table：嵌套回归结果之间的差异

Question

我正在比较两种替代策略，以使用data.tablefor 包估计数据子集的线性回归模型R。这两种策略产生相同的系数，因此它们看起来是等价的。这个外表是骗人的。我的问题是：

为什么存储在lm模型中的数据不同？

library(data.table)

dat = data.table(mtcars)

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

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

乍一看，这两种方法似乎产生了相同的结果：

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

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

但是，如果我尝试从（扩展的）model.frame 中提取数据，我会得到不同的结果：

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

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

这是一个简单的最小工作示例。我真正的用例是，在sandwich::vcovCL为我的模型应用计算的聚类标准误差时，我得到了完全不同的结果。

编辑：

我接受@TimTeaFan 的回答（出色的侦探工作！），但在这里为未来的读者添加了一些有用的信息。

正如@achim-zeileis 在其他地方指出的那样，我们可以在全球环境中复制类似的行为：

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)

这似乎不是一个data.table特定的问题。一般来说，我们在重新评估模型中的数据时必须小心。

score 6 · Accepted Answer

我没有完整的答案，但我能够在一定程度上查明问题。

当我们比较两个模型的输出时，我们可以看到结果是相等的，除了调用不同（这是有道理的，因为它们实际上是不同的）：

# 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

所以看起来初始调用在这两种方法下都可以正常工作，一旦我们尝试访问底层数据，问题就会出现。

如果我们看一下如何expand.model.frame获取它的数据，我们可以看到它调用eval(model$call$data, envir)whereenvir被定义为environment(formula(model))与对象的公式关联的环境lm。

如果我们查看每个模型的关联环境中的数据并将其与我们期望它保存的数据进行比较，我们可以看到第二种方法产生了我们期望的数据，而.SD在调用中使用的第一种方法产生了一些不同的数据。

我仍然不清楚为什么会发生什么以及发生了什么，但我们现在知道问题出在对.SD. 我首先想到，这可能是由命名 a 引起的data.table .SD，但是在使用了数据为 a 的模型之后，data.table这.SD似乎不是问题所在。

# 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

添加在

我试着深入挖掘一下，看看发生了什么。首先我打电话debug(as.formula)，然后在每次迭代中查看以下对象：

object
ls(environment(object))

我们可以看到，在“策略 2”中，每个公式都与不同的环境相关联，当查看环境时，我们看到它包含一个对象x，当检查时 ( environment(object)$x) 包含预期的mtcars数据。

然而，在“策略 1”中，我们可以观察到每次调用都将相同的环境与正在创建的公式as.formula相关联。此外，在检查环境时，我们可以看到它填充了子集数据的单个向量（例如，等）以及一些函数（例如，，mtcarsamcarbcyl.POSIXtCfastmeanstrptimeETC。）。这可能是事情出错的地方。我怀疑当将相同的环境与两个不同的公式（模型）相关联时，第一个模型的基础数据在计算第二个模型时会“更新”。这也应该是模型输出本身正确的原因。在计算第一个模型之前，数据仍然是正确的。它被第二个模型覆盖，因此也是正确的。但是当之后访问底层数据时，事情变得一团糟。

旁注
我很好奇我们是否可以在使用 tidyverse 时观察到类似的问题和差异，expand.model.frame答案是“是”。在这里，新的rowwise表示法会引发错误，而方法group_map和map方法一样有效：

# 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")

r - R data.table：嵌套回归结果之间的差异

编辑：

1 回答 1

Related

Reference