处理增长曲线中的缺失值答案

【问题标题】：Handling missing values in growthcurver处理增长曲线中的缺失值
【发布时间】：2020-03-15 05:16:57
【问题描述】：

R 包growthcurver 非常适合有效分析和可视化生物体生长，除非存在缺失值。因为我有宽格式的数据（每列都是一个变量）并且每个变量的时间是随机的，所以有大量的NAs。不幸的是，growthcurver 包不喜欢NAs，所以现在我只能选择 2 个选项。

选项 A
- 通过逻辑回归或机器学习方法对缺失数据进行插补（我不喜欢此选项，因为我尝试过 mice、Hmisc 进行回归插补但失败了，因为变量（列）多于每列中的观察值和caret 用于随机森林插补，它没有产生任何有意义的插补值）。然后，插补还会将我的数据框创建为主要是我无法证明的插补值。
选项 B
- 以某种方式调整growthcurver 函数以比当前更好地处理NAs。我尝试使用该功能四处寻找，但找不到可以放入简单 na.omit() 的位置。

这是与一次性函数 SummarizeGrowth() 一起使用的代码（当我手动删除 NAs 时）。我应该注意，当一个人只有几个观察要分析/可视化时，这个函数很有用，但理想情况下，我会使用函数SummarizeGrowthByPlate()，它是一个包派生的apply()函数，它循环遍历每一列（变量）自动生成可视化和结果。

选项 C
- 希望 SO 社区有一个快速修复！

示例数据框

        time      a      b      c      d      e      f      g
1    0.00002     NA     NA     NA     NA     NA     NA     NA
2    0.00003     NA     NA     NA     NA     NA     NA 0.0000
3   22.00000     NA     NA     NA     NA     NA     NA     NA
4   24.01000 0.1443 0.1554 0.0999 0.1110 0.0999 0.0666     NA
5   24.03000     NA     NA     NA     NA     NA     NA 0.0666
6   28.00000     NA     NA     NA     NA     NA     NA     NA
7   36.00000 0.2220 0.2775 0.2775 0.1776 0.1221 0.1221     NA
8   39.00000     NA     NA     NA     NA     NA     NA 0.2442
9   40.00000     NA     NA     NA     NA     NA     NA     NA
10  44.00000 0.3330 0.3885 0.3552 0.3108 0.2664 0.1998     NA
11  46.00000     NA     NA     NA     NA     NA     NA     NA
12  64.00000     NA     NA     NA     NA     NA     NA 0.7881
13  67.00000 0.9435 1.2210 1.1655 0.9990 1.5984 0.5217     NA
14  88.00000 1.8093 1.8093 1.8093 1.8870 1.6872 1.5096     NA
15 108.00000     NA     NA     NA     NA     NA     NA 1.6983

可重现的数据

df <- structure(list(time = c(2e-05, 3e-05, 22, 24.01, 24.03, 28, 36, 
39, 40, 44, 46, 64, 67, 88, 108), a = c(NA, NA, NA, 0.1443, NA, 
NA, 0.222, NA, NA, 0.333, NA, NA, 0.9435, 1.8093, NA), b = c(NA, 
NA, NA, 0.1554, NA, NA, 0.2775, NA, NA, 0.3885, NA, NA, 1.221, 
1.8093, NA), c = c(NA, NA, NA, 0.0999, NA, NA, 0.2775, NA, NA, 
0.3552, NA, NA, 1.1655, 1.8093, NA), d = c(NA, NA, NA, 0.111, 
NA, NA, 0.1776, NA, NA, 0.3108, NA, NA, 0.999, 1.887, NA), e = c(NA, 
NA, NA, 0.0999, NA, NA, 0.1221, NA, NA, 0.2664, NA, NA, 1.5984, 
1.6872, NA), f = c(NA, NA, NA, 0.0666, NA, NA, 0.1221, NA, NA, 
0.1998, NA, NA, 0.5217, 1.5096, NA), g = c(NA, 0, NA, NA, 0.0666, 
NA, NA, 0.2442, NA, NA, NA, 0.7881, NA, NA, 1.6983)), class = "data.frame", row.names = c(NA, 
-15L))

成功，但需要使用 SummarizeGrowth() 从单个列中手动删除 NA

library(growthcurver)

SummarizeGrowth(df$time[!is.na(df$a)], df$a[!is.na(df$a)])

Fit data to K / (1 + ((K - N0) / N0) * exp(-r * t)): 
    K   N0  r
  val:  2.121   0.004   0.085
  Residual standard error: 0.02857429 on 2 degrees of freedom

Other useful metrics:
  DT    1 / DT  auc_l   auc_e
  8.13  1.2e-01 38.16   38.77

不使用 SummarizeGrowth() 手动删除 NA 时失败

SummarizeGrowth(df$time, dfb$a)

Fit data to K / (1 + ((K - N0) / N0) * exp(-r * t)): 
    K   N0  r
  val:  0   0   0
  Residual standard error: 0 on 0 degrees of freedom

Other useful metrics:
  DT    1 / DT  auc_l   auc_e
  0 0   0   0

Note: cannot fit data

尝试使用自动 SummarizeGrowthByPlate() 时失败

SummarizeGrowthByPlate(df)

sample k n0 r t_mid t_gen auc_l auc_e sigma            note
1      a 0  0 0     0     0     0     0     0 cannot fit data
2      b 0  0 0     0     0     0     0     0 cannot fit data
3      c 0  0 0     0     0     0     0     0 cannot fit data
4      d 0  0 0     0     0     0     0     0 cannot fit data
5      e 0  0 0     0     0     0     0     0 cannot fit data
6      f 0  0 0     0     0     0     0     0 cannot fit data
7      g 0  0 0     0     0     0     0     0 cannot fit data

【问题讨论】：

growthcurver cran.r-project.org/web/packages/growthcurver/vignettes/… 的小插图

标签： r apply missing-data imputation

【解决方案1】：

老鼠的问题，Hmisc 是它们没有进行时间序列插补。他们只关注变量间的相关性。这意味着当一行完全不适用时 - 他们无法为该行计算任何内容。（从逻辑上讲，必须至少有一个回归器才能执行回归）

由于您的每个变量似乎在时间上有明显的相关性，您可以查看时间序列插补/插值。

有 imputeTS 包，它提供了很多时间序列插补算法。但是在这里很难使用它，因为它需要等间隔的时间序列（意味着每行之间的相同时间差）作为输入。要使用此包，您首先必须将时间序列转换为等距。对于这种特定情况，这似乎不是一个好主意。

据我所知 zoo 包可以对不规则间隔时间序列执行时间序列插补。所以这个包可能是你的最佳选择。我会专门尝试na.approx() - 线性插值函数。

很遗憾，我无法快速给出一个可行的示例。用法基本上是：

library(zoo)
na.approx(zooobject)

您现在唯一需要弄清楚的是如何将您的 df 转换为动物园系列（这是必需的输入）

只是为了展示它可能值得付出努力 - 这是一个使用 imputeTS 的工作示例（之前您不需要动物园对象）

library(imputeTS)
na_interpolation(df)


        time      a      b      c      d      e      f        g
1    0.00002 0.1443 0.1554 0.0999 0.1110 0.0999 0.0666 0.000000
2    0.00003 0.1443 0.1554 0.0999 0.1110 0.0999 0.0666 0.000000
3   22.00000 0.1443 0.1554 0.0999 0.1110 0.0999 0.0666 0.022200
4   24.01000 0.1443 0.1554 0.0999 0.1110 0.0999 0.0666 0.044400
5   24.03000 0.1702 0.1961 0.1591 0.1332 0.1073 0.0851 0.066600
6   28.00000 0.1961 0.2368 0.2183 0.1554 0.1147 0.1036 0.125800
7   36.00000 0.2220 0.2775 0.2775 0.1776 0.1221 0.1221 0.185000
8   39.00000 0.2590 0.3145 0.3034 0.2220 0.1702 0.1480 0.244200
9   40.00000 0.2960 0.3515 0.3293 0.2664 0.2183 0.1739 0.380175
10  44.00000 0.3330 0.3885 0.3552 0.3108 0.2664 0.1998 0.516150
11  46.00000 0.5365 0.6660 0.6253 0.5402 0.7104 0.3071 0.652125
12  64.00000 0.7400 0.9435 0.8954 0.7696 1.1544 0.4144 0.788100
13  67.00000 0.9435 1.2210 1.1655 0.9990 1.5984 0.5217 1.091500
14  88.00000 1.8093 1.8093 1.8093 1.8870 1.6872 1.5096 1.394900
15 108.00000 1.8093 1.8093 1.8093 1.8870 1.6872 1.5096 1.698300

可能这些结果已经比使用横截面数据的插补包更合理。但请记住，imputeTS 假设有规律的间隔时间序列 - 如果你可以让 zoo 工作，你可以获得更好的结果，因为它也考虑了不规则的间隔。

【讨论】：

也在这里看看如何使用来自动物园的 na.approx stackoverflow.com/questions/52795960/…

【解决方案2】：

这是解决 NA 问题的最佳方法（遵循 Angel Angelov 的 https://rpubs.com/angelov/growthcurver）。

models.all <- lapply(df[2:ncol(df)], function(x) SummarizeGrowth(df[!is.na(x), 1], x[!is.na(x)]))
df.predicted.plate <- data.frame(time = df$time)

for (i in names(df[2:ncol(df)])) 
{df.predicted.plate[[i]] <- predict(models.all[[i]]$model, newdata = list(t = df$time))}

melt1 <- melt(df, id.vars = "time", variable.name = "sample", value.name = "od")
melt2 <- melt(df.predicted.plate, id.vars = "time", variable.name = "sample", value.name = "pred.od")

df.final <- cbind(melt1, pred.od=melt2[,3])

【讨论】：