如何根据不同 R 生态系统中的另一个向量重写代码，将函数应用于行子集？

Posted 2023-02-14

【中文标题】如何根据不同 R 生态系统中的另一个向量重写代码，将函数应用于行子集？

【英文标题】：How can I re-write code that applies a function on subset of rows based on another vector in different R ecosystems?

【发布时间】：2022-01-20 18:38:38

【问题描述】：

在我的问题中，我必须根据从原始数据中提取的一组日期对单个时间序列的子集应用一个函数。 因此，我有一个 data.frame，其中包含 2005-01-01 和 2010-12-31 (test_final_ind_series) 之间的每个人的时间序列，以及理想情况下从相同的数据。

有了这些，在我的示例中，我尝试根据 sample_events 中的个人和日期来计算时间序列值 exp 子集的平均值。

我以两种不同的方式做到这一点：

1：一个简单但有效的代码，可以快速完成工作 我只是要求用户输入特定个人的数据并定义时间滞后和窗口宽度（如滚动平均值）。然后函数exp_summary 输出请求的平均值。

为了对sample_events 中的每一行重复该操作，我决定按个人 ID 嵌套单个系列，然后附上日期样本。最后，我只是运行一个循环，将该函数应用于每个单独的嵌套数据帧。

#Sample data
set.seed(111)
exp_series <- data.frame(
  id = as.character(rep(1:10000, each=2191)), 
  date = rep(seq(as.Date('2005-01-01'),
                 as.Date('2010-12-31'), by = 'day'),times=10000),
  exp = rep(rnorm(n=10000, mean=10, sd=5),times=2191)
)


sample_dates <- data.frame(
  Event_id = as.character(replicate(10000,sample(1:10000,size = 1,replace = TRUE))), 
  Event_date = sample(
    seq(as.Date('2005-01-01'),
        as.Date('2010-12-31'), by = 'day'),
    size =10000,replace = TRUE)
)



#This function, given a dataframe with dates and exposure series (df) 
#an event_date
#a lag value
#a width of the window
#Outputs the average for a user-defined time window
exp_summary<- function(df, event_date, lag=0,width=0)
    df<-as.data.table(df)
    end<-as.character(as.Date(event_date)-lag)
    start<-as.character(max(as.Date(end)-width, min(df$date)))# I need this in case the time window goes beyond the time limits (earliest date)
    return(mean(df[date %between% c(start,end)]$exp))


#Nest dataframes
exp_series_nest <- exp_series %>% 
  group_by(id) %>% 
  nest()


#Merge with sample events, including only the necessary dates
full_data<-merge(exp_series_nest,sample_dates, by.x="id", by.y="Event_id",all.x = FALSE, all.y=TRUE)


#Initialize dataframe in advance
summaries1<-setNames(data.frame(matrix(ncol = 2, nrow = nrow(full_data))), c("id", "mean"))
summaries1$id<-full_data$id

#Loop over each id, which is nasted data.frame
system.time(for (i in 1:nrow(full_data))
  summaries1$mean[i]<-exp_summary(full_data$data[[i]], full_data$Event_date[i], lag=1, width=365)
)

2：使用高度灵活的包 runner

使用相同的数据，我需要正确指定参数。我还在 Github 存储库上打开了 issue 以通过并行化加速此代码。

system.time(summaries2 <- sample_dates %>%
  group_by(Event_id) %>%
  mutate(
    mean = runner(
      x = exp_series[exp_series$id ==  Event_id[1],], 
      k = "365 days", 
      lag = "1 days",
      idx =exp_series$date[exp_series$id == Event_id[1]],
      at = Event_date,
      f = function(x) mean(x$exp),
      na_pad=FALSE
    )
  )    
)

它们给出了完全相同的结果，直到小数点后第二位，但方法 1 比方法 2 快得多，当您使用非常数据集时，您可以看到差异。

我的问题是，对于方法 1，如何在 data.table 和/或 tidyverse 生态系统中以更简洁的方式编写最后一个循环？我真的很难将嵌套列表和嵌入在同一数据框中的“普通”列一起工作。

另外，如果您有任何其他建议，我愿意听取！我在这里更多是出于好奇而不是需要，因为我的问题已经通过方法 1 解决了。

【问题讨论】：

您的示例不是很直观，sample_dates 中的某些行是重复的。我可以为方法1提供一些要点。 1.嵌套数据不是必需的，尝试用组替换它。 2.合并可能被左连接代替，例如exp_series[sample_dates, on = c(id = "Event_id")] 3. exp_summary 可以按组改写为frollmean。

感谢您的评论。 sample_dates 中的 ID 可以在我的问题中重复。其余的 cmets 并不是很有用，这意味着，我已经有预感，我可以使用不同的函数，如 frollmean 或 slider，但我无法找到一种方法将列表中的数据和其中的单个列一起工作环境。

另外，感谢第 2 点。我不知道如何在 data.table 中使用 join。

【参考方案1】：

使用data.table，您可以将exp_series 与sample_dates 中您希望的范围连接起来并计算平均值by=.EACHI：

library(data.table)

setDT(exp_series)
setDT(sample_dates)


lag <- 1
width <- 365 
# Define range
sample_dates[,':='(begin=Event_date-width-lag,end=Event_date-lag)]

# Calculate mean by .EACHI
summariesDT <- exp_series[sample_dates,.(id,mean=mean(exp))
                                      ,on=.(id=Event_id,date>=begin,date<=end),by=.EACHI][
                                      ,.(id,mean)]

请注意，这仅对 Event_id 返回与 summaries1 相同的结果，而在 sample_dates 中没有重复。

重复的情况下结果不同，例如Event_id==1002:

sample_dates[Event_id==1002]
   Event_id Event_date      begin        end
     <char>     <Date>     <Date>     <Date>
1:     1002 2010-08-17 2009-08-16 2010-08-16
2:     1002 2010-06-23 2009-06-22 2010-06-22

如果您的真实数据中没有重复，这应该不是问题。

【讨论】：

谢谢！在我的问题中，ID 可以重复，并且您的代码在这种情况下似乎仍然有效，因为它独立地为每个重复的行生成平均值。有趣的是，即使与您的代码相比，方法 1 仍然稍微快一些。如果您有足够的 RAM，请尝试将 exp_series 中的 id 数增加到 100k 并将 sample_dates 中的行数增加到 100k

我没有足够的 RAM 来使用循环测试 100k 行。 data.table 跑了 40 秒。关于相同的结果，我注意到summaries1-summariesDT 并不总是返回 0，尤其是在重复的情况下。这种情况下你觉得结果OK吗？

哦，我明白了，可能是因为“data.table”并不总是保持重复行的原始顺序吗？我检查了几个重复的 ID。