#r #data.table
#r #данные.таблица
Вопрос:
Я сравниваю две альтернативные стратегии для оценки моделей линейной регрессии на подмножествах данных с использованием data.table пакета for 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 специфическая проблема. И в целом, мы должны быть осторожны при переоценке данных из модели.
Ответ №1:
У меня нет полного ответа, но я смог в какой-то степени точно определить проблему.
Когда мы сравниваем выходные данные двух моделей, мы видим, что результат одинаков, за исключением вызовов, которые отличаются (что имеет смысл, поскольку они на самом деле разные):
# 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) where envir определяется как 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 связывает одну и ту же среду с создаваемой формулой. Кроме того, при проверке среды мы можем видеть, что она заполнена отдельными векторами подмножеств mtcars данных (например am , carb , cyl и т.д.), А также некоторыми функциями (например .POSIXt , Cfastmean , strptime и т.д.). Вероятно, именно здесь все идет наперекосяк. Я подозреваю, что при сопоставлении одной и той же среды с двумя разными формулами (моделями) базовые данные первых моделей «обновляются» при вычислении второй модели. Это также должно быть причиной, по которой сам вывод модели является правильным. На момент вычисления первой модели данные остаются правильными. Она перезаписывается второй моделью, которая, следовательно, тоже верна. Но при последующем доступе к базовым данным все становится беспорядочным.
Дополнительное примечание
Мне было любопытно, можем ли мы наблюдать аналогичные проблемы и различия в 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")
Комментарии:
1. Спасибо за эту детективную работу! При использовании
identicalвместоall.equal, вы видите, чтоcall,terms,modelотличаются междуmod1иmod2где среда вызывает различия в последних двух. Кратчайший способ извлечь это, вероятноenvironment(terms(mod1$models[[1]])), etc. Как вы правильно указали дляmod1env одинаковы для двух моделей, поскольку данные, по-видимому, обрабатываются последовательно в одной и той же среде. Чтоexpand.model.frame()касается этой среды для переоценки данных, то здесь все идет не так. Также не знаю, почему именно это происходит.