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

2020-10-23 16:27:27 浏览数 (2)

在前面

本期大猫课堂将继续上期的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)

0 人点赞