R：基于另一个矩阵和运行 rle 的子集矩阵

Question

样本数据：

      year <- rep(1980:2015, each = 365) 
      doy <- rep(1:365, times = 36)

      set.seed(125) 
      val <- sample(0:1, size = 365*36,replace = TRUE) 
      mat <- as.matrix(cbind(year,doy,val))

mat 是包含年份、doy 和值为 1 或 0 的列的数据。

我还有另一个矩阵mat1。

      set.seed(123) 
      mat1 <- apply(matrix(sample(c(230:365), replace = TRUE, size = 4L * 36L), nrow = 36L), 2L, sort)
      mat1 <- t(apply(mat1, 1, function(x) x[order(x)]))
      head(mat1)

            [,1] [,2] [,3] [,4]
      [1,]  230  231  233  236
      [2,]  235  238  242  242
      [3,]  236  242  243  246
      [4,]  243  245  247  249
      [5,]  247  248  249  250
      [6,]  249  250  253  263

我想使用mat1 每年从mat 进行子集化。例如，mat 中的 1980 年将是 分为三组：

      group 1 from 230 till 231 (1st and second column of row 1 from `mat1`)
      group 2 from 232 till 233 (second column + 1 to third column of row 1 from `mat1`)
      group 3 from 234 till 236 (third column + 1 to fourth column of row 1 from `mat1`)

这将为我提供 1980 年的三个向量。我想对每个向量进行 rle 以找到最长的 连续出现 1. 类似

       group1.rle <- rle(group1)
       group2.rle <- rle(group2)
       group3.rle <- rle(group3)

       max(group1.rle$lengths[group1.rle$values == 1])
       max(group2.rle$lengths[group2.rle$values == 1])
       max(group3.rle$lengths[group3.rle$values == 1])

然后多年来一直重复此操作mat。

谢谢。

Answer 1

这可能有点过头了，但我在 tidyverse 中迷失了（谁不知道那种感觉？;)）

包

# You might as well go with library(tidyverse)
library(dplyr)
library(purrr)
library(tidyr) 

代码和说明

# Preparation
mat <- as.data.frame(mat)
colnames(mat1) <- c("D1", "D2", "D3", "D4")
mat1 <- cbind(year = 1980:2015, mat1)
mat1 <- as.data.frame(mat1)

我将列命名为mat1，添加列year，并将mat 和mat1 都转换为漂亮的tibbles，所以我在它们上使用left_join。这对于在正确的年份获得正确的指标列很重要。

mat_new <- mat %>% 
  left_join(mat1, by = "year") %>% 
  mutate(group1 = (doy >= D1 & doy <=D2), 
         group2 = (doy >= D2 + 1 & doy <=D3),  
         group3 = (doy >= D3  + 1 & doy <=D4)) 

mat_new

# A tibble: 13,140 x 10
#     year   doy   val    D1    D2    D3    D4 group1 group2 group3
#    <int> <int> <int> <int> <int> <int> <int>  <lgl>  <lgl>  <lgl>
#  1  1980     1     1   230   231   233   236  FALSE  FALSE  FALSE
#  2  1980     2     0   230   231   233   236  FALSE  FALSE  FALSE
#  3  1980     3     0   230   231   233   236  FALSE  FALSE  FALSE
#  4  1980     4     0   230   231   233   236  FALSE  FALSE  FALSE
#  5  1980     5     1   230   231   233   236  FALSE  FALSE  FALSE
#  6  1980     6     1   230   231   233   236  FALSE  FALSE  FALSE
#  7  1980     7     1   230   231   233   236  FALSE  FALSE  FALSE
#  8  1980     8     0   230   231   233   236  FALSE  FALSE  FALSE
#  9  1980     9     1   230   231   233   236  FALSE  FALSE  FALSE
# 10  1980    10     1   230   231   233   236  FALSE  FALSE  FALSE
# ... with 13,130 more rows

所以首先加入两个tibble 然后添加组列doy 是否在子集中。

mat_new <- mat_new %>% 
  gather(group, indicator, group1:group3) %>% 
  nest(doy, val, indicator) 

mat_new

# A tibble: 108 x 7
# year    D1    D2    D3    D4  group               data
# <int> <int> <int> <int> <int>  <chr>             <list>
#  1  1980   230   231   233   236 group1 <tibble [365 x 3]>
#  2  1981   235   238   242   242 group1 <tibble [365 x 3]>
#  3  1982   236   242   243   246 group1 <tibble [365 x 3]>
#  4  1983   243   245   247   249 group1 <tibble [365 x 3]>
#  5  1984   247   248   249   250 group1 <tibble [365 x 3]>
#  6  1985   249   250   253   263 group1 <tibble [365 x 3]>
#  7  1986   250   250   255   269 group1 <tibble [365 x 3]>
#  8  1987   255   258   259   269 group1 <tibble [365 x 3]>
#  9  1988   259   259   263   274 group1 <tibble [365 x 3]>
# 10  1989   261   270   273   285 group1 <tibble [365 x 3]>
# ... with 98 more rows

将组列收集在一起，然后nest 数据。现在每一行都包含一年 - 组组合，并且数据列存储此特定组合的 doy、val 和 indicator 列。这样可以更轻松地在下一步中为全年 - 组组合计算 rle。

mat_new <- mat_new %>% 
 mutate(group.rle = map(data, ~ .x %>% filter(indicator) %>% pull(val) %>% rle), 
        max.group.rle = map_dbl(group.rle, ~max(.x$lengths[.x$values == 1])))

mat_new

# A tibble: 108 x 9
#     year    D1    D2    D3    D4  group               data group.rle max.group.rle
#    <int> <int> <int> <int> <int>  <chr>             <list>    <list>         <dbl>
#  1  1980   230   231   233   236 group1 <tibble [365 x 3]> <S3: rle>             1
#  2  1981   235   238   242   242 group1 <tibble [365 x 3]> <S3: rle>             2
#  3  1982   236   242   243   246 group1 <tibble [365 x 3]> <S3: rle>             1
#  4  1983   243   245   247   249 group1 <tibble [365 x 3]> <S3: rle>             1
#  5  1984   247   248   249   250 group1 <tibble [365 x 3]> <S3: rle>          -Inf
#  6  1985   249   250   253   263 group1 <tibble [365 x 3]> <S3: rle>             1
#  7  1986   250   250   255   269 group1 <tibble [365 x 3]> <S3: rle>             1
#  8  1987   255   258   259   269 group1 <tibble [365 x 3]> <S3: rle>             2
#  9  1988   259   259   263   274 group1 <tibble [365 x 3]> <S3: rle>          -Inf
# 10  1989   261   270   273   285 group1 <tibble [365 x 3]> <S3: rle>             2
# ... with 98 more rows

通过两次调用map，我们可以获得每个组的最大权限。在第一次调用中，data 列中的每个tibble 都由存储在indicator 列中的值过滤，然后提取val 列（使用pull），最后应用rle到这些值。

在第二个map 调用中，存储在group.rle 列中的rle 将根据您的条件（仅值1 的长度）进行过滤，并计算max。因为这返回一个长度为 1 的数字向量，所以我使用 map_dbl 直接存储结果。

注意，此调用将产生警告，因为并非所有组都包含值 1，因此在过滤后最大值没有非缺失参数。

mat_new  %>% 
  select(year, group, max.group.rle) %>% 
  spread(group, max.group.rle)

# A tibble: 36 x 4
#     year group1 group2 group3
# *  <int>  <dbl>  <dbl>  <dbl>
# 1   1980      1   -Inf      1
# 2   1981      2      1   -Inf
# 3   1982      1   -Inf   -Inf
# 4   1983      1   -Inf      1
# 5   1984   -Inf      1   -Inf
# 6   1985      1   -Inf      2
# 7   1986      1      3      1
# 8   1987      2   -Inf      1
# 9   1988   -Inf   -Inf      2
# 10  1989      2      1      3
# # ... with 26 more rows  

为了更好地查看结果，我只选择列year、group、max.group.rle，然后使用spread，将组分散到单独的列中。现在我们每年有一条线路提供相关信息。