按组插值估算缺失答案

【问题标题】：impute missing with interpolation by groups按组插值估算缺失
【发布时间】：2020-06-19 19:33:19
【问题描述】：

我正在尝试通过多组插值来估算缺失值 NA。

我只是子集一个简单的例子：

Year    ST   CC   ID     MP     PS 
2002    15   3     3     NA    1.5
2003    15   3     3     NA    1.5
2004    15   3     3    193    1.5
2005    15   3     3    193    1.5
2006    15   3     3    348    1.5
2007    15   3     3    388    1.5
2008    15   3     3    388    1.5
1999    53   33    1     NA    3.4 
2000    53   33    1     NA    3.4        
2002    53   33    1     NA    2.9           
2003    53   33    1     NA    2.6           
2004    53   33    1     NA    2.6         
2005    53   33    1    170    3.8         
2006    53   33    1    170    3.0           
2007    53   33    1    330    4.2         
2008    53   33    1    330    5.0

我使用了na.approx()，但弄错了。似乎我的数据在每组的第一次观察中缺少价值。

setDT(dt)[, MP_interpolate := na.approx(MP, na.rm = T), .(Year, ST, CC, ID)]

setDT(dt)[, MP_interpolate := if(length(na.omit(MP))<2) MP else na.approx(MP, na.rm=TRUE), .(Year, ST, CC, ID)]

我也试过包imputeTS，但没用。

两者都不正确。这是否意味着使用插值来估算缺失不是一个好主意？

我不确定哪种插补方式更好（不想使用mean 或median）。我正在考虑使用PS 趋势来估算MP。（只是我的想法，不是问题）

【问题讨论】：

标签： r data.table data-manipulation missing-data imputets

【解决方案1】：

您可以在ave 中尝试imputeTS::na_kalman。这也可以推断出您可能想要什么。

library(imputeTS)

dt$MP.imp <- with(dt, ave(MP, ST, CC, ID, FUN=na_kalman))
#    Year ST CC ID  MP  PS   MP.imp
# 1  2001 15  3  3  NA 1.5 193.0000
# 2  2002 15  3  3  NA 1.5 193.0000
# 3  2003 15  3  3  NA 1.5 193.0000
# 4  2004 15  3  3 193 1.5 193.0000
# 5  2005 15  3  3 193 1.5 193.0000
# 6  2006 15  3  3 348 1.5 348.0000
# 7  2007 15  3  3 388 1.5 388.0000
# 8  2008 15  3  3 388 1.5 388.0000
# 9  1999 53 33  1  NA 3.4 170.2034
# 10 2000 53 33  1  NA 3.4 166.3867
# 11 2002 53 33  1  NA 2.9 164.4496
# 12 2003 53 33  1  NA 2.6 165.0018
# 13 2004 53 33  1  NA 2.6 168.6527
# 14 2005 53 33  1 170 3.8 170.0000
# 15 2006 53 33  1 170 3.0 170.0000
# 16 2007 53 33  1 330 4.2 330.0000
# 17 2008 53 33  1 330 5.0 330.0000

您可能希望在 ave 中包含 year，否则在您的示例数据中将无法正常工作。

【讨论】：

您能解释一下na_kalman 的工作原理吗？在ID=1 中，MP.imp 减少和增加。看起来很酷。
@PeterChen 它基于Kalman 1960 提出的空间状态模型。卡尔曼滤波器应用基于水平、趋势和季节性分量的平滑算法来内推/外推缺失，并且还处理非平稳数据。因此它不同于线性插值。
谢谢。我一开始还以为是根据另一列来推算的，比如根据PS的趋势推算MP...
@PeterChen 实际上它只是以单变量方式插入分组的时间序列。对于多元插值，您可以做出预测，例如在loess regression。或使用多重插补，如mice（参数）或missForest（非参数）。但这导致了一个相当统计的讨论，更适合Cross Validated。
loess regression 有很大帮助。由于专栏中有很多NA，看来我需要研究多元插值。谢谢

【解决方案2】：

如果您在组的开头和结尾有NA 值，则应使用na.rm = FALSE 来获得与输入大小相同的输出。使用na.rm = TRUE，将删除前导和尾随NAs。

此外，您不应该按year 分组，因为那样只会给您一个组中的一个观察结果。

library(data.table)

setDT(dt)
dt[, MP_interpolate := as.integer(zoo::na.approx(MP, na.rm = FALSE)), .(ST, CC, ID)]

数据

dt <- structure(list(Year = c(2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 
2008L, 1999L, 2000L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 
2008L), ST = c(15L, 15L, 15L, 15L, 15L, 15L, 15L, 53L, 53L, 53L, 
53L, 53L, 53L, 53L, 53L, 53L), CC = c(3L, 3L, 3L, 3L, 3L, 3L, 
3L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L), ID = c(3L, 
3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), 
MP = c(NA, NA, 193L, 193L, 348L, 388L, 388L, NA, NA, NA, 
NA, NA, 170L, 170L, 330L, 330L), PS = c(1.5, 1.5, 1.5, 1.5, 
1.5, 1.5, 1.5, 3.4, 3.4, 2.9, 2.6, 2.6, 3.8, 3, 4.2, 5)), 
class = "data.frame", row.names = c(NA, -16L))

【讨论】：

我在新栏目中有NA 或TRUE。
它对我有用。我还在那里包装了as.integer，所以实际上你不应该得到逻辑值。你能检查一下我在帖子中分享的数据吗？
收到警告消息：Warning messages: 1: In [.data.table(setDT(dt), , :=(MP_interpolate, as.integer(na.approx(MP, : Group 1 column 'MP_interpolate': 193 (type 'integer') at RHS position 3 taken as TRUE when assigning to type 'logical' 2: In [.data.table(setDT(dt), , :=(MP_interpolate, as.integer(na.approx(MP, : Group 2 column 'MP_interpolate': 170 (type 'integer') at RHS position 6 taken as TRUE when assigning to type 'logical'
您对我分享的数据也有同样的理解吗？你的packageVersion('data.table') 是什么？我有‘1.12.2’。
packageVersion('data.table'): ‘1.12.8’ 和你的样本，我在每个组的第一次观察中仍然得到NA...$ MP_interpolate: int NA NA 193 193 348 388 388 NA NA NA ...