检测并纠正纵向数据中的单位误差

Question

我有一个非常大的纵向数据集，其中有数千人的权重。使用的重量单位（磅/千克）与重量一起输入。数据集中存在很大比例的错误，其中输入的权重单位与所选单位不同。显示随时间变化的体重曲线（转换为千克）时，单个时间点的单位误差明显突出，因为这些测量值通常与个人体重曲线偏差 2.2 倍（1 千克 = 2.20462 磅），如下图所示.

我意识到这个问题没有完美的解决方案，并且考虑到数据集的大小，不可能手动检查或纠正个别错误。所以我正在寻找一种有效的算法来检测最明显的错误并纠正它们。

我当前的方法并不完美，但包含一个主题循环，其中：

1. 我检查该科目的最大和最小权重条目的比率是否大于1.5，如果是，我转到2），如果不是，我转到下一个科目
2. 如果受试者只有 2 个观测值且比率大于 3.5，则通过将最大观测值乘以 1/2.20462、将最小值乘以 2.20462 来计算校正权重，如果比率小于 3.5，则将校正权重设置为不适用
3. 如果受试者有超过 2 个观察值：观察值按与中位数的绝对偏差降序排列。迭代排序的观测值，通过将权重乘以 1/2.20462 或 2.20462 来计算受试者观测值的标准差是否减小，如果是这样，则将校正后的权重设置为权重乘以相关因子，并根据关于下一步修正后的权重

我有两个问题：

1. R 中是否已经实现了针对此类问题的标准方法？
2. 我的算法或类似的算法可以在 tidyverse 中紧凑地实现吗？

这是一个示例数据集：

weight_dat <- structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
                                    2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                    3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
                                    4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), 
                             years_bl = c(0, 1.46, 2.34, 3.06, 5.28, 6.36, 7.01, 
                                          9.37, 10.08, 10.31, 0, 4.45, 4.76, 
                                          7.07, 10.3, 11.33, 0, 1.13, 1.86, 
                                          2.26, 5.67, 7.37, 7.5, 11.94, 12.31, 
                                          0, 0.7, 2.54, 4.26, 4.55, 5.79, 7.93, 
                                          9.89, 11.85, 0, 1.04, 4.01, 5.98, 
                                          8.16, 9.57, 13.68), 
                             weight_kg = c(78.7, 78.8, 75.8, 79.9, 82.8, 77.1, 
                                           81.2, 80.8, 82.9, 85.3, 74.1, 75.7, 
                                           78, 80.3, 90.5, 85, 90.2, 196.5, 
                                           85.7, 79.8, 83, 83, 80, 74.8, 71, 
                                           39.5, 81.5, 80, 79, 79.1, 77.8, 70.7, 
                                           33.3, 65.4, 72.1, 75.6, 79, 67.4, 
                                           69.4, 64.5, 67.3)), 
                        class = 'data.frame')

Answer 1

您基本上是在寻找异常值。有多种方法可以做到这一点，但它们都取决于您如何定义异常值，即给定值必须与其余值相差多远才能被视为异常值。

正如您所知，自动化方法始终取决于误差的大小，以及误差的大小与（正确测量的）值如何变化的关系。如果您确定引入的误差远大于随时间的预期变化，我建议使用两个常用选项：

使用 z 分数

z 分数，指示每个值与平均值的距离，并以标准差来衡量。您可以通过减去平均值并将结果除以标准差来计算它们。由于每个人都有自己的平均值和标准差，因此您需要为每个 ID 单独计算它们：

library(dplyr)
library(ggplot2)

weight_zs <-
  weight_dat %>%
  mutate_at("id", factor) %>%
  group_by(id) %>%
  mutate(z = (weight_kg - mean(weight_kg)) / sd(weight_kg))

z 分数为 2 将表示给定值比该个体的平均值高 2 个标准差。然后，您需要找到构成离群值的适当阈值：1.96 是通常的选择，因为 -1.96SD 和 1.96SD 之间的区间包含正态分布变量值的 95%（我不会更进一步）这个，否则我们就会进入交叉验证领域）。在下面的代码中，我使用

abs(z) > 2

来检测两侧的异常值，但如果您绝对确定所有错误都是正数，则可以只使用

z > 2

。

ggplot(weight_zs, aes(x = years_bl, y = weight_kg, colour = id)) +
geom_line() +
geom_point(aes(size = abs(z) > 2))

在您的情况下，2 似乎太高（有一些出乎意料的低值没有被检测为异常值）。您可以测试其他值，也许是在数据的子集上，但恐怕总是需要反复试验。

中位数 + IQR

另一个常见的定义是在

ggplot2

的箱线图中使用的定义， 将离群值阈值设置为中位数 +/- 1.5 * 四分位数范围：

weight_out <-
  weight_dat %>%
  mutate_at("id", factor) %>%
  group_by(id) %>%
  mutate(outlier = (weight_kg > (median(weight_kg) + 1.5 * IQR(weight_kg))) |
                   (weight_kg < (median(weight_kg) - 1.5 * IQR(weight_kg))))

weight_out %>%
  ggplot(aes(x = years_bl, y = weight_kg, colour = id)) +
  geom_line() +
  geom_point(aes(size = outlier))

这似乎效果更好一些。

其他注意事项

我相信如果你想实现自己的算法，你可以轻松地修改上面的代码。

case_when

函数可能有助于第一步。另请注意，我上面展示的方法和您的算法都假设平均值（或中位数）随着时间的推移保持不变。例如，如果您正在测量预计体重会随着时间的推移而增长的个体，那么当个体小于或大于其体重时，您更有可能（自动）在时间序列的开始和结束时找到异常值。平均/中位权重，使用整个时间序列计算得出。这可以通过为每个个体拟合回归模型并通过检查其残差或库克距离来找到异常值之类的方法来解决。