用最接近的值替换 R 中的 NA答案

【问题标题】：Replacing NAs in R with nearest value用最接近的值替换 R 中的 NA
【发布时间】：2012-04-22 01:47:04
【问题描述】：

我在 zoo 包中寻找类似于 na.locf() 的东西，但不是总是使用 previous 非NA 值我想使用 最近的非NA值。一些示例数据：

dat <- c(1, 3, NA, NA, 5, 7)

将NA 替换为na.locf（3 结转）：

library(zoo)
na.locf(dat)
# 1 3 3 3 5 7

和na.locf 与fromLast 设置为TRUE（5 向后携带）：

na.locf(dat, fromLast = TRUE)
# 1 3 5 5 5 7

但我希望使用最近非NA 值。在我的示例中，这意味着应该将 3 前移到第一个 NA，并将 5 后移到第二个 NA：

1 3 3 5 5 7

我编写了一个解决方案，但想确保我没有重新发明轮子。是不是已经有什么东西飘过来了？

仅供参考，我当前的代码如下。也许如果不出意外，有人可以建议如何提高效率。我觉得我缺少一个明显的改进方法：

  na.pos <- which(is.na(dat))
  if (length(na.pos) == length(dat)) {
    return(dat)
  }
  non.na.pos <- setdiff(seq_along(dat), na.pos)
  nearest.non.na.pos <- sapply(na.pos, function(x) {
    return(which.min(abs(non.na.pos - x)))
  })
  dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]

回答以下smci的问题：

不，任何条目都可以是 NA
如果都是 NA，则保持原样
没有。我当前的解决方案默认为左侧最接近的值，但这没关系
这些行通常有几十万个元素，因此理论上上限为几十万。实际上，这里和那里最多只有几个，通常是一个。

更新所以事实证明我们正在朝着完全不同的方向前进，但这仍然是一个有趣的讨论。谢谢大家！

【问题讨论】：

你看过na.locf的可选参数了吗？ fromLast 看起来可以做你想做的事。
它没有，因为它只是在相反的方向上取先前的值。它不会找到最近的非NA值
'介意发布您的解决方案吗？我很想看看你有什么。
刚刚做了，意识到如果不出意外，这可能会变成如何让我做得更好
我们可以遍历 rle(which(is.na(dat)))。并不是说这是最有效的，但它是一种改进。另请参阅 "How can I count runs in R?"，它需要调整 rle.na() 来处理 NA。

标签： r na missing-data

【解决方案1】：

代码如下。最初的问题没有完全明确，我曾要求这些澄清：

是否保证至少第一个和/或最后一个条目是非 NA？ [否]
如果一行中的所有条目都是 NA，该怎么办？ [保持原样]
您是否关心关系是如何拆分的，即如何处理1 3 NA NA NA 5 7 中的中间 NA？ [不关心/离开]
在最长的连续 NA 跨度上是否有上限 (S)？（如果 S 很小，我正在考虑递归解决方案。或者如果 S 很大并且行数和列数很大，则使用 ifelse 的数据框解决方案。）[最坏情况 S可能病态很大，因此不应使用递归]

geoffjentry，你的解决方案是你的瓶颈将是 nearest.non.na.pos 的串行计算和串行分配 dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]] 对于长度为 G 的大间隙，我们真正需要计算的是第一个（G/2，向上取整）项从左填充，其余从右填充。（我可以使用 ifelse 发布答案，但看起来很相似。）您的标准是运行时、big-O 效率、临时内存使用还是代码易读性？

Coupla 可能的调整：

只需要计算一次N <- length(dat)
普通情况下速度提升：if (length(na.pos) == 0) 跳过行，因为它没有 NAs
if (length(na.pos) == length(dat)-1) 只有一个非 NA 条目的（罕见）情况，因此我们用它填充整行

大纲解决方案：

遗憾的是 na.locf 不适用于整个数据帧，您必须按行使用 sapply：

na.fill_from_nn <- function(x) {
  row.na <- is.na(x)
  fillFromLeft <- na.locf(x, na.rm=FALSE) 
  fillFromRight <- na.locf(x, fromLast=TRUE, na.rm=FALSE)

  disagree <- rle(fillFromLeft!=fillFromRight)
  for (loc in (disagree)) { ...  resolve conflicts, row-wise }
}

sapply(dat, na.fill_from_nn)

或者，既然您说连续的 NA 很少见，请执行一个快速而愚蠢的 ifelse 从左侧填充孤立的 NA。这将在数据帧方面进行操作 => 使常见情况快速。然后使用逐行 for 循环处理所有其他情况。（这会影响长时间NA中中间元素的抢七，但你说你不在乎。）

【讨论】：

@joran，如果 cols 的数量很高并且 S 也很高，最有效的答案将会有很大的不同。它还会根据没有 NA 的行的比例而有所不同。
我希望能找到跑得更快的东西。我最初的假设是我缺少类似于 na.locf() 的东西，它比我拥有的要快得多。
@geoffjentry：添加常见情况的速度增强：if (length(na.pos) == 0) 我建议（这种情况有多普遍？）
我实际上在我的真实代码中有它，我把它拿出来，因为我认为有人会指出它没有多大作用。如果 na.pos 的长度为 0，则 sapply() 将简单地迭代任何内容。

【解决方案2】：

我想不出一个明显简单的解决方案，但是，在查看了建议（特别是smci 使用rle 的建议）后，我想出了一个似乎更有效的复杂函数。

这是代码，下面我会解释：

# Your function
your.func = function(dat) {
  na.pos <- which(is.na(dat))
  if (length(na.pos) == length(dat)) {
    return(dat)
  }
  non.na.pos <- setdiff(seq_along(dat), na.pos)
  nearest.non.na.pos <- sapply(na.pos, function(x) which.min(abs(non.na.pos - x)))
  dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]
  dat
}

# My function
my.func = function(dat) {
    nas=is.na(dat)
    if (!any(!nas)) return (dat)
    t=rle(nas)
    f=sapply(t$lengths[t$values],seq)
    a=unlist(f)
    b=unlist(lapply(f,rev))
    x=which(nas)
    l=length(dat)
    dat[nas]=ifelse(a>b,dat[ ifelse((x+b)>l,x-a,x+b) ],dat[ifelse((x-a)<1,x+b,x-a)])
    dat
}


# Test
n = 100000
test.vec = 1:n
set.seed(1)
test.vec[sample(test.vec,n/4)]=NA

system.time(t1<-my.func(test.vec))
system.time(t2<-your.func(test.vec)) # 10 times speed improvement on my machine

# Verify
any(t1!=t2)

我的功能依赖于rle。我正在阅读上面的 cmets，但在我看来 rle 对 NA 工作得很好。用一个小例子来解释是最容易的。

如果我从一个向量开始：

dat=c(1,2,3,4,NA,NA,NA,8,NA,10,11,12,NA,NA,NA,NA,NA,18)

然后我得到所有 NA 的位置：

x=c(5,6,7,8,13,14,15,16,17)

然后，对于 NA 的每次“运行”，我都会创建一个从 1 到运行长度的序列：

a=c(1,2,3,1,1,2,3,4,5)

然后我再做一次，但我颠倒了顺序：

b=c(3,2,1,1,5,4,3,2,1)

现在，我可以比较向量 a 和 b：如果 ab 则向前看并获取 x+b 处的值。剩下的只是处理所有 NA 或 NA 在向量的末尾或开头运行时的极端情况。

可能有更好、更简单的解决方案，但我希望这可以帮助您入门。

【讨论】：

@smci：没问题。是否可以格式化用户名以便链接？
就像任何其他超链接：open-square-bracket 名称 close-square-bracket (URL) 或点击右侧的“帮助”。

【解决方案3】：

这是我的尝试。我从不喜欢在 R 中看到 for 循环，但在稀疏 NA 向量的情况下，它看起来实际上会更有效（下面的性能指标）。代码要点如下。

  #get the index of all NA values
  nas <- which(is.na(dat))

  #get the Boolean map of which are NAs, used later to determine which values can be used as a replacement, and which are just filled-in NA values
  namask <- is.na(dat)

  #calculate the maximum size of a run of NAs
  length <- getLengthNAs(dat);

  #the furthest away an NA value could be is half of the length of the maximum NA run
  windowSize <- ceiling(length/2)

  #loop through all NAs
  for (thisIndex in nas){
    #extract the neighborhood of this NA
    neighborhood <- dat[(thisIndex-windowSize):(thisIndex+windowSize)]
    #any already-filled-in values which were NA can be replaced with NAs
    neighborhood[namask[(thisIndex-windowSize):(thisIndex+windowSize)]] <- NA

    #the center of this neighborhood
    center <- windowSize + 1

    #compute the difference within this neighborhood to find the nearest non-NA value
    delta <- center - which(!is.na(neighborhood))

    #find the closest replacement
    replacement <- delta[abs(delta) == min(abs(delta))]
    #in case length > 1, just pick the first
    replacement <- replacement[1]

    #replace with the nearest non-NA value.
    dat[thisIndex] <- dat[(thisIndex - (replacement))]
  }

我喜欢您提出的代码，但我注意到我们正在计算矩阵中每个 NA 值与每个其他非 NA 索引之间的增量。我认为这是最大的性能猪。相反，我只是提取每个 NA 周围的最小尺寸邻域或窗口，并在该窗口内找到最近的非 NA 值。

因此，性能与 NA 的数量和窗口大小成线性关系——其中窗口大小（上限）是 NA 最大运行长度的一半。要计算 NA 的最大运行长度，可以使用以下函数：

getLengthNAs <- function(dat){
  nas <- which(is.na(dat))
  spacing <- diff(nas)
  length <- 1;
  while (any(spacing == 1)){        
    length <- length + 1;
    spacing <- diff(which(spacing == 1))
  }
    length
}

性能比较

#create a test vector with 10% NAs and length 50,000.
dat <- as.integer(runif(50000, min=0, max=10))
dat[dat==0] <- NA

#the a() function is the code posted in the question
a <- function(dat){
  na.pos <- which(is.na(dat))
    if (length(na.pos) == length(dat)) {
        return(dat)
    }
    non.na.pos <- setdiff(seq_along(dat), na.pos)
    nearest.non.na.pos <- sapply(na.pos, function(x) {
        return(which.min(abs(non.na.pos - x)))
    })
    dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]
    dat
}

#my code
b <- function(dat){
    #the same code posted above, but with some additional helper code to sanitize the input
    if(is.null(dat)){
      return(NULL);
    }

    if (all(is.na(dat))){
      stop("Can't impute NAs if there are no non-NA values.")
    }

    if (!any(is.na(dat))){
      return(dat);
    }

    #starts with an NA (or multiple), handle these
    if (is.na(dat[1])){
      firstNonNA <- which(!is.na(dat))[1]
      dat[1:(firstNonNA-1)] <- dat[firstNonNA]
    }

    #ends with an NA (or multiple), handle these
    if (is.na(dat[length(dat)])){
      lastNonNA <- which(!is.na(dat))
      lastNonNA <- lastNonNA[length(lastNonNA)]
      dat[(lastNonNA+1):length(dat)] <- dat[lastNonNA]
    }

    #get the index of all NA values
    nas <- which(is.na(dat))

    #get the Boolean map of which are NAs, used later to determine which values can be used as a replacement, and which are just filled-in NA values
    namask <- is.na(dat)

    #calculate the maximum size of a run of NAs
    length <- getLengthNAs(dat);

    #the furthest away an NA value could be is half of the length of the maximum NA run
    #if there's a run at the beginning or end, then the nearest non-NA value could possibly be `length` away, so we need to keep the window large for that case.
    windowSize <- ceiling(length/2)

    #loop through all NAs
    for (thisIndex in nas){
      #extract the neighborhood of this NA
      neighborhood <- dat[(thisIndex-windowSize):(thisIndex+windowSize)]
      #any already-filled-in values which were NA can be replaced with NAs
      neighborhood[namask[(thisIndex-windowSize):(thisIndex+windowSize)]] <- NA

      #the center of this neighborhood
      center <- windowSize + 1

      #compute the difference within this neighborhood to find the nearest non-NA value
      delta <- center - which(!is.na(neighborhood))

      #find the closest replacement
      replacement <- delta[abs(delta) == min(abs(delta))]
      #in case length > 1, just pick the first
      replacement <- replacement[1]

      #replace with the nearest non-NA value.
      dat[thisIndex] <- dat[(thisIndex - (replacement))]
    }
    dat
}

#nograpes' answer on this question
c <- function(dat){
  nas=is.na(dat)
  if (!any(!nas)) return (dat)
  t=rle(nas)
  f=sapply(t$lengths[t$values],seq)
  a=unlist(f)
  b=unlist(lapply(f,rev))
  x=which(nas)
  l=length(dat)
  dat[nas]=ifelse(a>b,dat[ ifelse((x+b)>l,x-a,x+b) ],dat[ifelse((x-a)<1,x+b,x-a)])
  dat
}

#run 10 times each to get average performance.
sum <- 0; for (i in 1:10){ sum <- sum + system.time(a(dat))["elapsed"];}; cat ("A: ", sum/10)
A:  5.059
sum <- 0; for (i in 1:10){ sum <- sum + system.time(b(dat))["elapsed"];}; cat ("B: ", sum/10)
B:  0.126
sum <- 0; for (i in 1:10){ sum <- sum + system.time(c(dat))["elapsed"];}; cat ("C: ", sum/10)
C:  0.287

所以它看起来像这段代码（至少在这些条件下），与问题中发布的原始代码相比提供了大约 40 倍的加速，并且比下面@nograpes 的答案提供了 2.2 倍的加速（尽管我想象一个 rle在某些情况下，解决方案肯定会更快——包括更富含 NA 的向量）。

【讨论】：

毕竟，我想我还是把它打包吧。 GitHub 上提供了完整代码以及用于验证正确行为的 RUnit 测试：github.com/trestletech/R-Utils

【解决方案4】：

这是一个非常快的。它使用findInterval 来查找原始数据中每个NA 应考虑的两个位置：

f1 <- function(dat) {
  N <- length(dat)
  na.pos <- which(is.na(dat))
  if (length(na.pos) %in% c(0, N)) {
    return(dat)
  }
  non.na.pos <- which(!is.na(dat))
  intervals  <- findInterval(na.pos, non.na.pos,
                             all.inside = TRUE)
  left.pos   <- non.na.pos[pmax(1, intervals)]
  right.pos  <- non.na.pos[pmin(N, intervals+1)]
  left.dist  <- na.pos - left.pos
  right.dist <- right.pos - na.pos

  dat[na.pos] <- ifelse(left.dist <= right.dist,
                        dat[left.pos], dat[right.pos])
  return(dat)
}

我在这里测试它：

# sample data, suggested by @JeffAllen
dat <- as.integer(runif(50000, min=0, max=10))
dat[dat==0] <- NA

# computation times
system.time(r0 <- f0(dat))    # your function
# user  system elapsed 
# 5.52    0.00    5.52
system.time(r1 <- f1(dat))    # this function
# user  system elapsed 
# 0.01    0.00    0.03
identical(r0, r1)
# [1] TRUE

【讨论】：

哇。我认为我们有一个赢家。非常好。
我相信你的第 4 行应该是 if (length(na.pos) == 0)
谢谢@Morten。实际上，这两种极端情况都需要提前退出。它是固定的。

【解决方案5】：

速度大约比所选答案慢 3-4 倍。不过我的很简单。这也是一个罕见的while循环。

f2 <- function(x){

  # check if all are NA to skip loop
  if(!all(is.na(x))){

    # replace NA's until they are gone
    while(anyNA(x)){

      # replace from the left
      x[is.na(x)] <- c(NA,x[1:(length(x)-1)])[is.na(x)]

      # replace from the right
      x[is.na(x)] <- c(x[-1],NA)[is.na(x)]
    }
  }

  # return original or fixed x
  x
}

【讨论】：

【解决方案6】：

我喜欢所有严谨的解决方案。虽然不是直接问什么，但我发现这篇文章正在寻找一种用插值填充 NA 值的解决方案。在查看这篇文章后，我发现 na.fill 在 zoo 对象（向量、因子或矩阵）上：

z <- c(1,2,3,4,5,6,NA,NA,NA,2,3,4,5,6,NA,NA,4,6,7,NA)
z1 <- zoo::na.fill(z, "extend")

注意 NA 值之间的平滑过渡

round(z1, 0)
#>  [1] 1 2 3 4 5 6 5 4 3 2 3 4 5 6 5 5 4 6 7 7

也许这会有所帮助

【讨论】：

我稍微编辑了这个答案，因为这个解决方案不需要转换成动物园对象