R Tricks: 如何处理Gaps & Islands问题？

写

在前面

本期大猫课堂将继续上期的R Tricks系列。在这一期中，大猫将向大家介绍“Gaps & Islands Problem”。这是在处理时间序列或者基因组数据中常见的一项任务。虽然常见，但要高效解决可不容易哦！

PS：大猫发现好多人给大猫留了言，但是因为超过48小时以后就不能回复大家了。所以如果小伙伴们有问题，可以再试着给大猫留言哦，大猫看到一定第一时间回复哈。

提

出问题

话说有个擅长使用SQL的小伙伴在StackOverflow上提出了这样一个问题，他说，Gaps & Islands问题在SQL中能很容易解决，那么在R中也能高效解决吗？

慢着——什么是Gaps & Islands问题？这个小伙伴举了个栗子来说明，让我们看一下。假如我们有如下数据集：

这是一个记录时间的数据集。每一行都有ID、起始时间（stime）、结束时间（etime）。我们可以发现，第1至4行的时间是有重叠的，其中最早的起始时间是（2014-01-15 08:00:00），最晚的结束时间是（2014-01-15 11:00:00）。而第5与第6行的时间也有重叠。Gaps & Islands问题就是说：能不能把时间有重叠的行合并起来？也即最终的结果数据集应该是这个样子的：

这个任务就像要把相连的岛屿链接起来，而把不相连的岛屿分割开，这也就是Gaps & Islands名称的由来。实际上，大猫听说在处理基因数据的时候也常遇见这个问题，但是大猫自己没有接触，欢迎有经验的小伙伴分享经历哦。

那么，如何处理Gaps & Islands问题呢？

原问题大家可以访问以下链接：http://stackoverflow.com/questions/30629894/how-to-solve-gaps-and-island-problems-in-r-and-performance-vs-sql

（生成样例数据集的代码附在见文末）

解

题思路

在解决本问题的过程中我们需要用到data.table包！

”

我们的思路很简单，分成四步：

▶ 将数据集按照ID与起始时间（stime）进行排序

▶ 找到结束时间（etime）的累计最大值

▶ 一旦完成以上两步，那么重叠的行即为当前结束时间（etime)累计最大值仍旧大于下一行的观测

▶ 标记出这些行，合并。

解

题步骤

首先，我们将原数据集按照ID以及起始时间（stime）排序：

▶ setorder(dat, ID, stime)

其次，也是最关键的一步，我们需要建立一个新变量etime.max。顾名思义，它记录了每个ID中结束时间的累计最大值：

▶ dat[, etime.max := as.POSIXct(cummax(as.numeric(etime)), origin = '1970-01-01'), by = ID]

结果是这个样子的：

上一行代码中，使用的关键函数是累计最大值函数cummax。此外，由于cummax不能直接处理日期格式，所以需要先将日期转化为数字进行比较，完了再转换回日期。

接下来，我们需要新建一个grp分组变量，它用于将一个个“islands”区分开来——即如果当前行的stime小于etime.max，那么grp的数字不变（意味着观测之间有重叠）；但如果stime比etime.max要大，那么grp则+1，代表现在出现了一个gap，我们进入了一个“新的islands”。最终效果如下：

从上图中我们可以看到，1-4行的grp值都为0，说明属于同一组；而5-6行的grp值为1，说明属于新的一组。为了实现这一步，我们的代码是：

▶ dat[, grp := cumsum(c(FALSE, stime[2:.N] > etime.max[1:(.N - 1)]))[1:.N], by = ID]

其中，stime[2, .N]表示截取stime向量的第2个元素至最后一个元素，etime.max[1, (.N - 1)]表示截取etime.max向量的第1个元素至倒数第二个元素。cumsum(stime[2:.N] > etime.max[1:(.N - 1)])表示如果当前行的stime比上一行的etime.max的值要大，那么返回TRUE，同时grp+1（我们用cumsum函数完成grp的累加）；而如果比上一行小，那么返回FALSE，同时grp不累加。

关于如何巧用cumsum函数，大猫在上一期的《R Tricks：如何巧为分组观测编号》中也有详细讲解哦

最后，我们只要把每个grp组中起始时间（stime）的最小值和结束时间（etime）的最大值找出来就行啦：

▶ dat[, .(stime = min(stime), etime = max(etime)), by = .(ID, grp)][, grp := NULL][]

结果如下：

本

期总结

本期大猫带领大家学习了如何处理Gaps & Islands问题——也即如何合并时间上有重叠的观测。我们灵活使用了cummax与cumsum函数，他们在处理分组数据的时候尤其有用。关于如何巧用cumsum函数，大猫在上一期的《R Tricks：如何巧为分组观测编号》中也有详细讲解哦。

我是大猫，咱们下期见！

附：样例数据集生成代码

▶ dat <- structure(

list(ID = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L),

stime = structure(c(as.POSIXct("2014-01-15 08:00:00"),

as.POSIXct("2014-01-15 10:00:00"),

as.POSIXct("2014-01-15 08:30:00"),

as.POSIXct("2014-01-15 09:00:00"),

as.POSIXct("2014-01-15 11:30:00"),

as.POSIXct("2014-01-15 12:00:00"),

as.POSIXct("2014-01-15 07:30:00"),

as.POSIXct("2014-01-15 08:00:00"),

as.POSIXct("2014-01-15 08:30:00"),

as.POSIXct("2014-01-15 09:00:00"),

as.POSIXct("2014-01-15 09:00:00"),

as.POSIXct("2014-01-15 09:30:00"),

as.POSIXct("2014-01-15 10:00:00")),

class = c("POSIXct", "POSIXt"), tzone = ""),

etime = structure(c(as.POSIXct("2014-01-15 09:30:00"),

as.POSIXct("2014-01-15 11:00:00"),

as.POSIXct("2014-01-15 10:00:00"),

as.POSIXct("2014-01-15 09:30:00"),

as.POSIXct("2014-01-15 12:30:00"),

as.POSIXct("2014-01-15 13:00:00"),

as.POSIXct("2014-01-15 08:30:00"),

as.POSIXct("2014-01-15 09:00:00"),

as.POSIXct("2014-01-15 09:30:00"),

as.POSIXct("2014-01-15 10:00:00"),

as.POSIXct("2014-01-15 10:00:00"),

as.POSIXct("2014-01-15 10:30:00"),

as.POSIXct("2014-01-15 11:00:00")

),

class = c("POSIXct", "POSIXt"), tzone = "")

),

.Names = c("ID", "stime", "etime"),

sorted = c("ID", "stime", "etime"),

class = c("data.table", "data.frame"),

row.names = c(NA, -13L)

)

setDT(dat)