如果我理解你,这里是一个使用data.table
包的解决方案。我找到了两个选项(但第一个具有更好的性能)
将原始数据框转换为data.table
对象:
dt <- data.table(df) # Create a data table from the data frame
setorder(dt, patient_id, visit_date) # Sort by patient_id, then by visit_date
定义周阈值参数:
weekNum = 20L # Considering a threshold of: 20-weeks.
选项 1:直接计算visit_date
列中的周数
我们定义了以下函数,对每个组进行计算:
visitFreq <- function(x) {
n <- length(x)
result <- numeric(n)
if (n > 1) {
for (i in 1:n) {
# For each row of the column by patient_id
ref <- x[i] # reference date
x.prev <- x[x < ref] # select previous dates
if (length(x.prev) > 0) {
x.prev <- sapply(x.prev, function(y) {
ifelse(difftime(ref, y, units = "weeks") <= weekNum, 1, 0)
})
result[i] <- sum(x.prev)
}
}
}
return(result)
}
对于每一个x[i]
,它都会找到之前访问的次数,然后计算之前的日期是否在定义的阈值内。然后只剩下在阈值内统计之前的访问次数。
一旦我们知道如何进行计算,我们只需要将这个函数应用于visit_date
每个列patient_id
:
dt[, visits := visitFreq(visit_date), by = patient_id]
注意:visitFreq
必须考虑一个向量函数来定义该函数,该函数接收一个数组visit_date
并且应该返回一个相同维度的数组。
选项 2:创建一个人工变量来收集给定患者的所有就诊日期。
现在我们需要创建一个函数来计算周数:
calc <- function(vec, x) {
vec.prev <- vec[vec < x] # Select all dates before x
n <- 0
if (length(vec.prev) > 0) {
vec.prev <- sapply(vec.prev, function(y) {
ifelse(difftime(x, y, units = "weeks") <= weekNum, 1, 0)
})
n <- sum(vec.prev)
}
return(n)
}
在哪里:
我们仅按 date 之前的日期进行过滤x
。现在我们sapply
对 的每个元素应用该函数,以使用周数为单位vec
计算y
(的每个元素vec
)与参考日期之间的时间差。x
结果将是1
任何差异日期小于weekNum
或为零。然后,从参考日期开始少于特定周数的先前访问次数将仅计算1
我们获得的所有内容。
data.table
现在我们在这样的对象中使用这个函数:
dt[, visits := .(list(visit_date)), by = patient_id]
[, visits := mapply(calc, visits, visit_date)][order(patient_id)][]
让我们稍微解释一下:
- 我们创建一个
visits
变量,它是给定的所有日期的列表patient_id
(因为by
子句)。
如果我们执行第一个表达式,它将产生如下内容:
> dt[, visits := .(list(visit_date)), by = patient_id][]
visit_id patient_id visit_date visits
1: 1 1 2016-12-02 2016-12-02,2016-12-15,2016-12-30,2017-02-15
2: 4 1 2016-12-15 2016-12-02,2016-12-15,2016-12-30,2017-02-15
3: 3 1 2016-12-30 2016-12-02,2016-12-15,2016-12-30,2017-02-15
4: 7 1 2017-02-15 2016-12-02,2016-12-15,2016-12-30,2017-02-15
5: 2 2 2016-12-02 2016-12-02,2017-02-01
6: 6 2 2017-02-01 2016-12-02,2017-02-01
7: 5 3 2016-12-30 2016-12-30
8: 8 4 2017-02-10 2017-02-10
9: 9 5 2017-01-15 2017-01-15
10: 10 6 2017-03-01 2017-03-01
>
- 第二条语句(第二
[]
块)只是重新分配先前创建的变量进行计算visits
,但现在计算参考日期的次数或先前访问。我们需要该mapply
函数来进行向量计算,在每次调用cal
函数时,我们都有输入参数:(dt[i]$visits
一个列表)和对应的dt[i]$visit_date[i]
. mapply
只是遍历所有i
调用函数的元素calc
。
结果
最后,结果将是:
> dt
visit_id patient_id visit_date visits
1: 1 1 2016-12-02 0
2: 4 1 2016-12-15 1
3: 3 1 2016-12-30 2
4: 7 1 2017-02-15 3
5: 2 2 2016-12-02 0
6: 6 2 2017-02-01 1
7: 5 3 2016-12-30 0
8: 8 4 2017-02-10 0
9: 9 5 2017-01-15 0
10: 10 6 2017-03-01 0
>
我想这就是你想要的。
注意:这可能是一种即时计算的方法,但我不知道如何。也许其他人可以提出一种在语法上更简洁的方法。
表现
我想知道哪个选项具有更好的性能(我期望 OPC1),让我们检查一下:
library(microbenchmark)
op <- microbenchmark(
OP1 = copy(dt)[, visits := visitFreq(visit_date), by = patient_id],
OP2 = copy(dt)[, visits := .(list(visit_date)), by = patient_id][, visits := mapply(calc, visits, visit_date)],
times=100L)
print(op)
它产生以下输出:
Unit: milliseconds
expr min lq mean median uq max neval cld
OP1 3.467451 3.552916 4.165517 3.642150 4.200413 7.96348 100 a
OP2 4.732729 4.832695 5.799648 5.063985 6.073467 13.17264 100 b
>
因此,第一个选项具有最佳性能。
编辑(添加了由@DavidArenburg 提出的解决方案)
让我们将连接解决方案作为第三个选项,但增加重复输入向量的输入参数的大小,例如:
nSample <- 100
patient_id <- rep(c(1, 2, 1, 1, 3, 2, 1, 4, 5, 6), nSample)
visit_id <- 1:nSample
visit_date <- rep(as.Date(c('2016-12-02', '2016-12-02', '2016-12-30',
'2016-12-15', '2016-12-30', '2017-02-01',
'2017-02-15', '2017-02-10', '2017-01-15', '2017-03-01')), nSample)
df <- data.frame(visit_id, patient_id, visit_date)
opc3 <- function(df) {
df[, visit_date20 := visit_date - 20 * 7] # Create a 20 weeks boundry
## Count previous visits within the range
df[df, .(visits = .N),
on = .(patient_id, visit_date < visit_date, visit_date > visit_date20),
by = .EACHI]
}
dt <- data.table(df)
dt3 <- copy(dt)[, visit_date := as.IDate(visit_date)] # Convert visit_date to a proper Date class
library(microbenchmark)
op <- microbenchmark(
OP1 = copy(dt)[, visits := visitFreq(visit_date), by = patient_id],
OP2 = copy(dt)[, visits := .(list(visit_date)), by = patient_id][, visits := mapply(calc, visits, visit_date)],
OP3 = opc3(copy(dt3)),
times = 10L)
print(op)
我得到以下结果:
Unit: milliseconds
expr min lq mean median uq max neval cld
OP1 6315.73724 6485.111937 10744.808669 11789.230998 15062.957734 15691.445961 10 b
OP2 6266.80130 6431.330087 11074.441187 11773.459887 13928.861934 15335.733525 10 b
OP3 2.38427 2.845334 5.157246 5.383949 6.711482 8.596792 10 a
>
当行数增加时,@DavidArenburg 解决方案可以更好地扩展。