我有一个脚本，在一周内的若干天里，当温度超过某个阈值时，会生成指示变量（例如，连续两天或以上超过温度的第90百分位数）。

temp <- temp %>%
  mutate(
    warm=if_else(month < 3 | month > 8, 1, 0),
    per_warm=case_when(
      warm == 0 ~ 0,
      warm == 1 & month < 3 ~ year - 1991,
      warm == 1 & month >= 9 ~ year - 1990
    )
  )

## 2 consecutive days over the 90th percentile
temp <- temp %>%
  group_by(per_warm, comuna) %>%
  mutate(
    p90=round(as.numeric(quantile(tmax, 90/100, na.rm=TRUE)), 1)
    ) %>%
  ungroup()

temp <- temp %>%
  arrange(comuna, year, month, day)

## Heat day: tmax > percentile threshold
temp <- temp %>%
  mutate(
    hd90=if_else(tmax > p90, 1, 0),
    cd90=if_else(tmax > p90, 1, 0)
  )

## Consecutive days
temp <- temp %>%
  arrange(comuna, year, month, day) %>%
  group_by(comuna) %>%
  mutate(
    cd90=fcumsum(x=hd90, g=fcumsum(!hd90)),
    cd90=if_else(warm == 0, 0, cd90),
    hd90=if_else(warm == 0, 0, hd90)
  ) %>%
  ungroup()

## Heat week
## In the past 7 days, did the pregnant woman have at least 3
## consecutive days over p90?
## Is the delivery day at least the third consecutive day over p90?
## Important: does this 7-day window include the present day?

## Number of consecutive days over threshold
last <- 2

temp <- temp %>%
  group_by(comuna) %>%
  mutate(hk90=if_else(cd90 >= last, 1, 0)) %>%
  ungroup()

temp <- temp %>%
  group_by(comuna) %>%
  mutate(hk90=if_else(lag(cd90, 1, default=0) >= last, 1, hk90)) %>%
  ungroup()

temp <- temp %>%
  group_by(comuna) %>%
  mutate(hk90=if_else(lag(cd90, 2, default=0) >= last, 1, hk90)) %>%
  ungroup()

temp <- temp %>%
  group_by(comuna) %>%
  mutate(hk90=if_else(lag(cd90, 3, default=0) >= last, 1, hk90)) %>%
  ungroup()

temp <- temp %>%
  group_by(comuna) %>%
  mutate(hk90=if_else(lag(cd90, 4, default=0) >= last, 1, hk90)) %>%
  ungroup()

## Same for p95 and p98
#...

现在的问题在于，当我尝试修改参数时会出现问题（把 tmax 改为 tmean，或者把 2 天改为 3，或者把每周的计数天数改为1天）。

我可以把脚本复制粘贴多次，改动其中的具体部分，然后运行，以获得所有组合，总共18种：2个温度、计数天数为1到3天、阈值对应3个百分位数。我的想法是得到18个变量，如此排列：“tmax_p90_l1 tmax_p90_l2 tmax_p90_l3 .... tmean_p98_l1 tmean_p98_l2 tmean_p98_l3”。

但我完全相信，这可以用一种“优雅”和“高效”的方式来完成。

有什么建议？

（我编辑了帖子，使之更容易理解）

解决方案

这是一个基于 data.table 的函数，能够实现这一点，按照评论中建议的运行长度编码来实现。

library(data.table)

f_dt <- \(temp, last=2, tx='tmax', qs=c(0.90, 0.95, 0.98), round_p=FALSE) {
  qs <- setNames(qs, qs*100)
  qn <- names(qs)
  p_cols <- paste0("p", qn)
  hd_cols <- paste0("hd", qn)
  cd_cols <- paste0("cd", qn)
  hk_cols <- paste0("hk", qn)
  base_cols <- names(temp)
  heat_week <- \(cd, last) {
    z <- cd >= last
    as.integer(Reduce(`|`, shift(z, 0:4, fill=FALSE)))
  }
  dt <- as.data.table(temp)
  dt[, warm:=as.integer(month < 3 | month > 8)]
  dt[, per_warm:=0L]
  dt[warm == 1L & month < 3, per_warm:=year - 1991L]
  dt[warm == 1L & month >= 9, per_warm:=year - 1990L]
  setorder(dt, comuna, year, month, day)
  if (round_p) {  ## rounding here (not ideal; kept for compatibility)
    dt[, (p_cols) := as.list(
      round(
        unname(stats::quantile(get(tx), qs, na.rm=TRUE)),
        1
      )
    ), by=.(per_warm, comuna)]
  } else {
    dt[, (p_cols) := as.list(
      unname(stats::quantile(get(tx), qs, na.rm=TRUE))
    ), by=.(per_warm, comuna)]
  }
  for (j in seq_along(qn)) {
    p_col <- p_cols[j]
    hd_col <- hd_cols[j]
    cd_col <- cd_cols[j]
    hk_col <- hk_cols[j]
    dt[, (hd_col):=as.integer(get(tx) > get(p_col))]
    dt[, (cd_col):={
      x <- get(hd_col)
      rid <- rleid(x)
      cd <- sequence(tabulate(rid))
      cd[x == 0L] <- 0L
      cd
    }, by=comuna]
    dt[warm == 0L, c(hd_col, cd_col):=0L]
    dt[, (hk_col):=heat_week(get(cd_col), last), by=comuna]
  }
  setcolorder(dt, c(
    base_cols, "warm", "per_warm",
    p_cols, hd_cols, cd_cols, hk_cols
  ))
  dt
}

用法

> f_dt(temp) |> head(3)  ## using defaults
   comuna  year month   day     tmax    tmean  warm per_warm      p90      p95      p98  hd90  hd95  hd98  cd90  cd95  cd98  hk90  hk95  hk98
   <char> <int> <int> <int>    <num>    <num> <int>    <int>    <num>    <num>    <num> <int> <int> <int> <int> <int> <int> <int> <int> <int>
1:      a  1990     1     1 27.37383 17.15873     1       -1 31.57008 33.46768 35.40864     0     0     0     0     0     0     0     0     0
2:      a  1990     1     1 18.29882 21.79981     1       -1 31.57008 33.46768 35.40864     0     0     0     0     0     0     0     0     0
3:      a  1990     1     1 30.20806 19.31374     1       -1 31.57008 33.46768 35.40864     0     0     0     0     0     0     0     0     0
> f_dt(temp, last=3, tx='tmean', qs=c(0.75, 0.90)) |> head(3)  ## customize
   comuna  year month   day     tmax    tmean  warm per_warm      p75      p90  hd75  hd90  cd75  cd90  hk75  hk90
   <char> <int> <int> <int>    <num>    <num> <int>    <int>    <num>    <num> <int> <int> <int> <int> <int> <int>
1:      a  1990     1     1 27.37383 17.15873     1       -1 21.39404 24.57328     0     0     0     0     0     0
2:      a  1990     1     1 18.29882 21.79981     1       -1 21.39404 24.57328     1     0     1     0     0     0
3:      a  1990     1     1 30.20806 19.31374     1       -1 21.39404 24.57328     0     0     0     0     0     0

检查

> temp_OP <- f_OP(temp)  ## your code
> temp_dt <- f_dt(temp, round_p=TRUE)
> 
> stopifnot(all.equal(
+   as.data.frame(temp_dt),
+   as.data.frame(temp_OP[, sub("\\.", "", names(temp_dt))])
+ ))
> 
> microbenchmark::microbenchmark(
+   f_OP(temp),
+   f_dt(temp),
+   times=10L
+ )
Unit: milliseconds
       expr      min       lq     mean   median       uq      max neval cld
 f_OP(temp) 400.6320 423.1116 521.9866 514.2539 621.2008 660.8288    10  a 
 f_dt(temp) 168.6245 178.8905 184.7490 184.2289 194.0279 195.9787    10   b

它的执行速度也略快一些。

数据：

set.seed(42)

n <- 1e6

temp <- data.frame(
  comuna=sample(letters[1:5], n, replace=TRUE),
  year=sample(1990:1995, n, replace=TRUE),
  month=sample(1:12, n, replace=TRUE),
  day=sample(1:28, n, replace=TRUE),
  tmax=rnorm(n, mean=25, sd=5),
  tmean=rnorm(n, mean=18, sd=5)
)