处理 growthcurver 中的缺失值

Question

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

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

这是使用一次性函数 SummarizeGrowth() 的代码（当我手动删除 NAs 时）。我应该注意到，当一个人对 analyze/visualize 仅有少量观察时，此函数很有用，但理想情况下，我会使用函数 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

Answer 1

这是解决 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])

Answer 2

老鼠的问题，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 假设有规律的间隔时间序列——如果你能让动物园工作，你可以获得更好的结果，因为它也考虑了不规则的间隔。

处理 growthcurver 中的缺失值

Handling missing values in growthcurver

r

apply

missing-data

imputation