为什么在分组的data.table内使用update()函数会丢失线性模型数据?

12

好的,这是一个奇怪的问题。我怀疑这是 data.table 内部的一个 bug,但如果有人能解释为什么会发生这种情况 - update 到底在做什么,那就太有用了。

我在使用 data.table 中的 list(list()) 技巧来存储拟合模型。当您创建一系列针对不同分组的 lm 对象,然后更新这些模型时,所有模型的模型数据都变成了最后一组的模型数据。这似乎像是某个引用挂在了某个地方,而应该制作副本,但我找不到在 lmupdate 之外复现此问题的方法。

具体示例:

从鸢尾花数据开始,首先使三个物种具有不同的样本大小,然后对每个物种拟合一个 lm 模型,最后更新这些模型:

set.seed(3)
DT = data.table(iris)
DT = DT[rnorm(150) < 0.9]
fit = DT[, list(list(lm(Sepal.Length ~ Sepal.Width + Petal.Length))),
          by = Species]
fit2 = fit[, list(list(update(V1[[1]], ~.-Sepal.Length))), by = Species]

原始数据表中每个物种的数量不同。
DT[,.N, by = Species]
#       Species  N
# 1:     setosa 41
# 2: versicolor 39
# 3:  virginica 42

而第一个匹配确认了这一点:

fit[, nobs(V1[[1]]), by = Species]
#       Species V1
# 1:     setosa 41
# 2: versicolor 39
# 3:  virginica 42

但更新后的第二个适配器对所有型号都显示为42。
fit2[, nobs(V1[[1]]), by = Species]
#       Species V1
# 1:     setosa 42
# 2: versicolor 42
# 3:  virginica 42

我们还可以查看模型属性,该属性包含用于拟合的数据,并且可以看到所有模型确实使用了最终分组的数据。问题是这是如何发生的?
head(fit$V1[[1]]$model)
#   Sepal.Length Sepal.Width Petal.Length
# 1          5.1         3.5          1.4
# 2          4.9         3.0          1.4
# 3          4.7         3.2          1.3
# 4          4.6         3.1          1.5
# 5          5.0         3.6          1.4
# 6          5.4         3.9          1.7
head(fit$V1[[3]]$model)
#   Sepal.Length Sepal.Width Petal.Length
# 1          6.3         3.3          6.0
# 2          5.8         2.7          5.1
# 3          6.3         2.9          5.6
# 4          7.6         3.0          6.6
# 5          4.9         2.5          4.5
# 6          7.3         2.9          6.3
head(fit2$V1[[1]]$model)
#   Sepal.Length Sepal.Width Petal.Length
# 1          6.3         3.3          6.0
# 2          5.8         2.7          5.1
# 3          6.3         2.9          5.6
# 4          7.6         3.0          6.6
# 5          4.9         2.5          4.5
# 6          7.3         2.9          6.3
head(fit2$V1[[3]]$model)
#   Sepal.Length Sepal.Width Petal.Length
# 1          6.3         3.3          6.0
# 2          5.8         2.7          5.1
# 3          6.3         2.9          5.6
# 4          7.6         3.0          6.6
# 5          4.9         2.5          4.5
# 6          7.3         2.9          6.3

1
这种行为可以用以下更简单的方式再现:m1 <- fit$V1[[2]]; m2 <- update(m1, .~.); m1m2 会变得不同。或许这会有助于解决问题。 - Aaron left Stack Overflow
谢谢,我试了很久想把这个简化一点。 - Corvus
3
你读过https://dev59.com/LGvXa4cB1Zd3GeqPK6D1吗?它可能有关联。 - hadley
@hadley 谢谢 - 它肯定看起来有关联。如果更新正在查找“全局”范围内的对象,那么当它在 data.table 中时它会在哪里查找呢?值得注意的是它并没有捕获整个 data.table 列 - 只有最终分组的列。 - Corvus
1
我已经仔细查看了,但我也感到困惑。现在已将其标记为bug#2590,并会回头研究它。顺便说一句,你发现得很好,问题布局也很棒。 - Matt Dowle
@MatthewDowle -- 我在某种程度上详细说明了一些发现。不确定是否有data.table的解决方案。 - mnel
1个回答

6

这不是一个答案,但对于一个评论来说太长了。

每个生成的模型中,术语组件的 .Environment 是相同的。

e1 <- attr(fit[['V1']][[1]]$terms, '.Environment')
e2 <- attr(fit[['V1']][[2]]$terms, '.Environment')
e3 <- attr(fit[['V1']][[3]]$terms, '.Environment')
identical(e1,e2)
## TRUE
identical(e2, e3)
## TRUE

看起来 data.table 在每次对分组进行 j 的评估时都使用相同的 一小块内存(这是我的非技术术语,意思是高效的)。然而,当调用 update 时,它会使用此内存来重新拟合模型。这将包含上一个组的值。

因此,如果你在这方面弄虚作假,它就能正常工作。

fit = DT[, { xx <-list2env(copy(.SD))

             mymodel <-lm(Sepal.Length ~ Sepal.Width + Petal.Length)
             attr(mymodel$terms, '.Environment') <- xx
             list(list(mymodel))}, by= 'Species']





lfit2 <- fit[, list(list(update(V1[[1]], ~.-Sepal.Width))), by = Species]
lfit2[,lapply(V1,nobs)]
V1 V2 V3
1: 41 39 42
# using your exact diagnostic coding.
lfit2[,nobs(V1[[1]]),by = Species]
      Species V1
1:     setosa 41
2: versicolor 39
3:  virginica 42

不是长期解决方案,但至少是一个解决方法。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接