使用另一个时间序列的增长率插入时间序列(替换 NA)

Interpolate time series (replace NAs) using growth rate of another time series

假设我有这样一个数据集:

trt <- data.table(group = rep(c("a","b"), each = 5), 
                  val1= c(60,62,NA,NA,71, NA, 21,22,NA,25),
                  val2 = c(1,1,1,NA,2, 1,1,NA,NA,2),
                  reflev = rep(c(1.01, 1.03, 1.061, 1.104,1.159), 2))
trt[ , ref:= round(reflev/shift(reflev), 2), by = group]


> trt
    group val1 val2 reflev  ref
 1:     a   60    1  1.010   NA
 2:     a   62    1  1.030 1.02
 3:     a   NA    1  1.061 1.03
 4:     a   NA   NA  1.104 1.04
 5:     a   71    2  1.159 1.05
 6:     b   NA    1  1.010   NA
 7:     b   21    1  1.030 1.02
 8:     b   22   NA  1.061 1.03
 9:     b   NA   NA  1.104 1.04
10:     b   25    2  1.159 1.05

在每个组中,我想通过乘以先前的可用值(例如 shift(val1)lag(val1)) 与来自 ref 列的值。如果在非 NA 值之后的序列中出现多个 NA,则应使用先前插值作为起点对所有 NA 进行插值。

所以,这是我设想的计算方式:

    group val1          val2         reflev  ref
 1:     a   60            1           1.010   NA
 2:     a   62            1           1.030 1.02
 3:     a   62*1.03       1           1.061 1.03
 4:     a   62*1.03*1.04  1*1.04      1.104 1.04
 5:     a   71            2           1.159 1.05
 6:     b   NA            1           1.010   NA
 7:     b   21            1           1.030 1.02
 8:     b   22           1*1.03       1.061 1.03
 9:     b   22*1.04      1*1.03*1.04  1.104 1.04
10:     b   25            2           1.159 1.05

有什么想法吗?我能想到的一切都很脏,并且会涉及两个循环,一个用于组,一个用于所需的列。

这是我的 'quick an dirty' 解决方案。很高兴听到改进的可能性:

interpolate <- function(data, int.column){

  for(row in 2:nrow(data)){
    if(is.na(data[row,get(int.column)]) & !is.na(data[row-1, get(int.column)])){
      data[[int.column]][row] <- data[row,ref]*data[[int.column]][row-1]
      }
  }
  return(data[ , get(int.column)])
}

我对 int.column(要插入的列的名称)使用了如此奇怪的调用,因为我没能从适当的环境中调用 int.column,否则总是出错) .

然后 trt[ , val1:= interpolate(.SD,"val1"), by = group] 用于单列插值或

columns.to.int <- c("val1", "val2")
trt[ , (columns.to.int):= lapply(columns.to.int, function(x)interpolate(.SD,x)), by = group]

对于多个。

有没有更好的方法?

这是另一个选项:

cols <- paste0("val", 1L:2L)
trt[, paste0("prev", cols) := lapply(.SD, nafill, type="locf"), group, .SDcols=cols]

trt[, outval1 := fifelse(is.na(val1), prevval1 * cumprod(ref), val1), .(group, rleid(is.na(val1)))]

trt[, outval2 := fifelse(is.na(val2), prevval2 * cumprod(ref), val2), .(group, rleid(is.na(val2)))]

编辑多个 val 列。也许是这样的:

cols <- paste0("V", 1L:30L)
for (x in cols) {
    trt[, c("prev", "ri") := {
            v <- get(x)
            .(nafill(v, "locf"), rleid(is.na(v)))
        }, group]
    trt[, paste0("out", x) := {
            v <- get(x)
            fifelse(is.na(v), prev * cumprod(ref), v)
        }, .(group, ri)]
}

或使用 melt 会更快:

mDT <- melt(trt[, rn := .I], measure.vars=patterns("^V"))
mDT[, pv := nafill(value, "locf"), group]
mDT[, nv := fifelse(is.na(value), pv * cumprod(ref), value),
    .(group, variable, rleid(is.na(value)))]
dcast(mDT, rn + group + reflev + ref ~ variable, value.var="nv")

示例数据:

library(data.table)
set.seed(0L)
nc <- 30L
nr <- 3e3L
trt <- data.table(group = rep(1:(nr/5L), each=5L), 
    reflev = 1+runif(nr)/10,
    as.data.table(matrix(sample(c(NA,10,20,30), nc*nr, TRUE), ncol=nc)))
trt[ , ref:= round(reflev/shift(reflev), 2), by = group]