北京大学R语言教程(李东风)第24章：数据汇总

d.cancer <- read_csv(
  "data/cancer.csv", locale=locale(encoding="GBK"))

d.cancer |>
  summarise(
    nobs = n(),
    n = sum(!is.na(age)),
    mean.age=mean(age, na.rm=TRUE),
    sd.age=sd(age, na.rm=TRUE)) |>
  knitr::kable()

dwt <- tribble(
  ~age, ~n, ~meanheight,
  10, 20, 110,
  10, 10, 112,
  11, 25, 120
)
knitr::kable(dwt)

dwt |>
  summarise(
    nobs = sum(n), 
    meanheight = sum(meanheight*n)/sum(n)
  )

d.cancer |>
  summarise(
    mean_v0=mean(v0, na.rm=TRUE),
    sd_v0 = sd(v0, na.rm=TRUE),
    mean_1=mean(v1, na.rm=TRUE),
    sd_v1 = sd(v1, na.rm=TRUE)) |>
  knitr::kable(digits=2)

d.cancer |>
  summarise(
    across(c("v0", "v1", "age"), mean)) |>
  knitr::kable(digits=2)

d.cancer |>
  summarise(
    across(
      c("v0", "v1", "age"), 
      \(x) mean(x, na.rm=TRUE))) |>
  knitr::kable(digits=2)

d.cancer |>
  summarise(
    across(c("v0", "v1", "age"), 
           list(
             avg = \(x) mean(x, na.rm=TRUE),
             std = \(x) sd(x, na.rm=TRUE)))) |>
  knitr::kable(digits=2)

summarise(across(
  where(is.numeric), 
  list(
    变量后缀 = \(x) 统计函数名(x, <...>),
    <...>)))

d.cancer |>
  summarise(
    across(where(is.numeric),
    list(avg = \(x) mean(x, na.rm=TRUE), 
         std = \(x) sd(x, na.rm=TRUE)))) |>
  knitr::kable(digits=2)

d.cancer |>
  select(v0, v1) |>
  summarise(across(everything(),
    list(avg = \(x) mean(x, na.rm=TRUE), 
         std = \(x) sd(x, na.rm=TRUE))) ) |>
  knitr::kable(digits=2)

d.cancer |>
  select(v0, v1) |>
  summarise(across(everything(),
    list(avg = \(x) mean(x, na.rm=TRUE), 
         std = \(x) sd(x, na.rm=TRUE))) ) |>
  pivot_longer(
    everything(),
    names_sep ="_",
    names_to = c("variable", ".value")
    ) |>
  knitr::kable(digits=2)

24.2.1 修改多个变量

除了对多个变量进行汇总以外，
有时需要对多个变量进行修改。
例如，某个数据集中所有的数值型变量都用-999表示缺失值，
需要将其改为NA，
如果在mutate中逐个写出修改，
就太麻烦，
这时可以在mutate()中用across()指定要修改的列和修改方法，如：

d.999.orig <- tribble(
  ~id, ~a, ~b, ~c,
  1, 5, 10, 6,
  2, -999, 12, 8,
  3, 4, -999, -999,
  4, 7, 9, 10
)
d.999.mod <- d.999.orig |>
  mutate(across(
    a:c,
    \(x) ifelse(x == -999, NA, x)
  ))
d.999.mod

## # A tibble: 4 × 4
##      id     a     b     c
##   <dbl> <dbl> <dbl> <dbl>
## 1     1     5    10     6
## 2     2    NA    12     8
## 3     3     4    NA    NA
## 4     4     7     9    10

修改结果直接替换了原来的变量。
如果希望修改结果用新变量保存，
可以在across中用.names指定一个新变量命名模式，
如

d.999 <- d.999.orig |>
  mutate(across(
    a:c,
    \(x) ifelse(x == -999, NA, x),
    .names = "{.col}_na"
  ))
d.999

## # A tibble: 4 × 7
##      id     a     b     c  a_na  b_na  c_na
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     1     5    10     6     5    10     6
## 2     2  -999    12     8    NA    12     8
## 3     3     4  -999  -999     4    NA    NA
## 4     4     7     9    10     7     9    10

其中.names的取值是一个模式字符串，
模式中的{.col}特指原始变量名。

下面的函数expand_dates()，
输入一个数据框，
选择其中的所有日期列，
将其每一个都拆分为年、月、日：

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

expand_dates <- function(df){
  df |>
    mutate(across(
      where(is.Date),
      list(year = year, month = month, day = mday)
    ))
}

d.dates <- read_csv(
  'data/dates.csv', 
  locale=locale(encoding="GBK"),
  col_types=cols(
    `序号`=col_integer(),
    `出生日期`=col_date(format="%Y/%m/%d"),
    `发病日期`=col_date(format="%Y/%m/%d")
  ))
knitr::kable(d.dates[1:4,])

d.dates2 <- expand_dates(d.dates)
knitr::kable(d.dates2[1:4,])

24.2.2 针对多个变量进行筛选

across的类似函数也可以用在filter函数中，
同时指定多个变量进行筛选。
if_any指定是任一变量满足条件就选择该观测，
if_all指定所有变量满足该该条件才选择该观测。
如：

d.999 |>
  filter(if_any(a_na:c_na, is.na))

## # A tibble: 2 × 7
##      id     a     b     c  a_na  b_na  c_na
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     2  -999    12     8    NA    12     8
## 2     3     4  -999  -999     4    NA    NA

d.999 |>
  filter(if_all(a_na:c_na, \(x) !is.na(x)))

## # A tibble: 2 × 7
##      id     a     b     c  a_na  b_na  c_na
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     1     5    10     6     5    10     6
## 2     4     7     9    10     7     9    10

24.3.1 简单汇总

比如，
按不同性别分组计算人数与年龄平均值：

d.cancer |>
  group_by(sex) |>
  summarise(
    count = n(), 
    mean.age = mean(age, na.rm=TRUE)) |>
  knitr::kable(digits=2)

按不同性别分组计算v0和v1的人数、平均值、标准差：

d.cancer |>
  group_by(sex) |>
  summarise(across(
    c("v0", "v1"),
    list(
      count = ~n(), 
      avg = \(x) mean(x, na.rm=TRUE), 
      std = \(x) sd(x, na.rm=TRUE))) ) |>
  knitr::kable(digits=2)

上面的程序应该将v0_count和v1_count合并，
但是summarise在有多个变量和多个统计量时总会将变量两两搭配，
确实需要时可以分别统计再合并，结果程序会略繁琐一些：

bind_cols(
  d.cancer |>
    group_by(sex) |>
    summarise(count = n() ),
  d.cancer |>
    group_by(sex) |>
    summarise(
      across(c("v0", "v1"),
      list(
        avg = \(x) mean(x, na.rm=TRUE), 
        std = \(x) sd(x, na.rm=TRUE)))) |>
    select(-sex)
  ) |>
  knitr::kable(digits=2)

24.3.2 满足条件汇总

为了查询数值型变量取值满足某种条件的个数和比例，
可以将该条件用sum()和mean()函数统计，
比如，
对男女两组分别计算年龄在60岁以上的人数和比例：

d.cancer |>
  group_by(sex) |>
  summarise(
    nold = sum(age >= 60, na.rm=TRUE), 
    pold = nold / n()) |>
  knitr::kable(digits=4)

24.3.3 分组筛选

用group_by()分组后除了可以分组汇总，
还可以分组筛选：

d.cancer |>
  group_by(sex) |>
  filter(rank(desc(v0)) <= 2) |>
  arrange(sex, desc(v0)) |>
  knitr::kable()

以上程序按性别分组后，
在每组中找出疗前体积排名在前两名的。

24.3.4 分组计算新变量

在分组后也可以根据每组的统计量用mutate()定义新变量。
比如，每组的v0, v1除以本组的最大值（结果略）：

d.cancer |>
  group_by(sex) |>
  mutate(
    v0s = v0 / max(v0),
    v1s = v1 / max(v1))

24.3.5 按前后变化分组

有时我们的数据中行次序是有意义的，
可能需要根据某个变量的变化来分组。
如：

d.cid <- tribble(
  ~state, ~y,
  10, 1.5,
  10, 1.6,
  10, 1.8,
  20, 2.3,
  20, 2.2,
  10, 1.1,
  10, 1.2
)
knitr::kable(d.cid)

按state的变化生成一个分组变量：

d.cid |>
  mutate(
    grp = consecutive_id(state)) |>
  knitr::kable()

可以直接按临时生成的分组变量分组汇总，如：

d.cid |>
  group_by(grp = consecutive_id(state)) |>
  slice_head(n=1) |>
  knitr::kable()

24.4.1 简单例子

下面的程序对d.cancer数据框分性别与病理类型分别统计人数：

d.cancer |>
  group_by(sex, type) |>
  summarise(freq=n()) |>
  knitr::kable()

## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.

如果仅需要计算交叉分类频数，
不需要用group_by()，
可以用dplyr的count()函数，如：

d.cancer |>
  count(sex, type) |>
  knitr::kable()

下面的程序数出NHANES数据框中ID(受访者编码)与SurveyYr(考察年份)每一对组合的出现次数,
筛选出二次及以上者，并降序排列，仅显示前10行结果：

NHANES |>
  count(ID, SurveyYr) |>
  filter(n >= 2 ) |>
  arrange(desc(n)) |>
  slice_head(n = 10) |>
  knitr::kable()

24.4.2 分组信息残留问题

用group_by()交叉分组汇总后的结果不是普通的tibble，
总是带有外层分组信息，
最内层的分组信息不再有效。
如果有更多的分组变量，
每次用summarise()汇总都会去掉最内层的分组信息，
但保留其它层的分组信息。
不注意这种规定在后续的使用中可能会产生问题，
为此，
可以在summarise()函数中加选项.groups="drop"要求汇总后不保留任何分组信息，
或者汇总后再用ungroup()函数取消所有分组。

例如，
希望在用group_by()按照性别和病理类别交叉分类计算频数后求所有病人的总人数，
用了如下程序：

d.cancer |>
  group_by(sex, type) |>
  summarise(freq=n()) |>
  summarise(ntotal=sum(freq)) |>
  knitr::kable()

## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.

可以看出并没有能够计算总人数，
而是按原来交叉分类的外层分类即性别分组计算了总人数。
这是因为交叉分组计算频数后的结果仍按照外层分类变量sex分组，
所以summarise(ntotal=sum(freq))也是分两组进行的。
在第一个summarise()中加.groups = "drop"可以不保留分组信息，如：

d.cancer |>
  group_by(sex, type) |>
  summarise(freq=n(), .groups="drop") |>
  summarise(ntotal=sum(freq)) |>
  knitr::kable()

得到了需要的结果。

为了避免在汇总后仍保留分组信息，
在summarise()中加选项.groups="drop"，
这样可以在汇总结果中去掉残留的分组信息。
默认的选项是.groups="drop_last"，
即仅去掉最内层分组信息。
还可以要求.groups="keep"，
这时汇总结果中包含所有的分组信息，
包括最内层分组。

用ungroup()去掉所有分组信息的例子如下：

d.cancer |>
  group_by(sex, type) |>
  summarise(freq=n()) |>
  ungroup() |>
  summarise(ntotal=sum(freq)) |>
  knitr::kable()

## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.

group_summ_1 <- function(df, group_var, ana_var){
  df |>
    group_by({{ group_var }}) |>
    summarise(
      mean = mean({{ana_var}}, na.rm=TRUE),
      max = max({{ana_var}}, na.rm=TRUE))
}
group_summ_1(d.class, sex, height)

group_summ_1(d.class, "sex", "height")
## # A tibble: 1 × 3
##   `"sex"`  mean max   
##   <chr>   <dbl> <chr> 
## 1 sex        NA height

group_summ_2 <- function(df, group_var, ana_var){
  df |>
    group_by(.data[[group_var]]) |>
    summarise(
      mean = mean(.data[[ana_var]], na.rm=TRUE),
      max = max(.data[[ana_var]], na.rm=TRUE))
}
group_summ_2(d.class, "sex", "height")

group_summ_3 <- function(df, group_var, ana_var){
  df |>
    group_by(pick({{ group_var }})) |>
    summarise(
      mean = mean({{ ana_var }}, na.rm=TRUE),
      max = max({{ ana_var }}, na.rm=TRUE))
}
group_summ_3(d.class, c(sex, age), height)

group_summ_4 <- function(df, group_var, ana_vars){
  df |>
    group_by({{ group_var }}) |>
    summarise(across(
      {{ ana_vars }},
      list(avg = ~ mean(., na.rm=TRUE), 
           max = ~ max(., na.rm=TRUE))))
}
group_summ_4(d.class, sex, c(height, weight))

group_summ_5 <- function(df, group_var, ana_vars){
  df |>
    group_by({{ group_var }}) |>
    summarise(across(
      all_of(ana_vars),
      list(avg = ~ mean(., na.rm=TRUE), 
           max = ~ max(., na.rm=TRUE))))
}
ana_vars <- c("height", "weight")
group_summ_5(d.class, sex, ana_vars)

24.7.1 summary()函数

对数值型向量x，用summary(x)可以获得变量的平均值、中位数、
最小值、最大值、四分之一和四分之三分位数。 如

summary(d.cancer[["v0"]])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.58   43.77   93.40  110.08  157.18  330.24

summary(d.cancer[["v1"]])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.30   12.73   30.62   44.69   72.94  128.34

v0是放疗前的肿瘤体积，
v1是放疗后的体积，
可以看出放疗后体积减小了很多。

可以用盒形图表现类似的信息，如

boxplot(d.cancer[,c("v0", "v1")], main="肿瘤体积")

对一个数据框d， 用summary(d)可以获得每个连续型变量的基本统计量，
和每个离散取值变量的频率。如

##        id             age            sex                type          
##  Min.   : 1.00   Min.   :49.00   Length:34          Length:34         
##  1st Qu.: 9.25   1st Qu.:55.00   Class :character   Class :character  
##  Median :17.50   Median :67.00   Mode  :character   Mode  :character  
##  Mean   :17.50   Mean   :64.13                                        
##  3rd Qu.:25.75   3rd Qu.:70.00                                        
##  Max.   :34.00   Max.   :79.00                                        
##                  NA's   :11                                           
##        v0               v1        
##  Min.   : 12.58   Min.   :  2.30  
##  1st Qu.: 43.77   1st Qu.: 12.73  
##  Median : 93.40   Median : 30.62  
##  Mean   :110.08   Mean   : 44.69  
##  3rd Qu.:157.18   3rd Qu.: 72.94  
##  Max.   :330.24   Max.   :128.34  
## 

对数据框d，用str(d)可以获得各个变量的类型和取值样例。 如

## spec_tbl_df [34 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ id  : num [1:34] 1 2 3 4 5 6 7 8 9 10 ...
##  $ age : num [1:34] 70 70 69 68 67 75 52 71 68 79 ...
##  $ sex : chr [1:34] "F" "F" "F" "M" ...
##  $ type: chr [1:34] "腺癌" "腺癌" "腺癌" "腺癌" ...
##  $ v0  : num [1:34] 26.5 135.5 209.7 61 237.8 ...
##  $ v1  : num [1:34] 2.91 35.08 74.44 34.97 128.34 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   id = col_double(),
##   ..   age = col_double(),
##   ..   sex = col_character(),
##   ..   type = col_character(),
##   ..   v0 = col_double(),
##   ..   v1 = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

用head(d)可以返回数据框（或向量、矩阵）的前几行，
用tail(d)可以返回数据框的后几行。

24.7.2 连续型变量概括函数

对连续取值的变量x，
可以用mean, std, var, sum, prod,
min, max等函数获取基本统计量。
加na.rm=TRUE选项可以仅对非缺失值计算。

sort(x)返回排序后的结果。
rev(x)把x所有元素次序颠倒后返回。
quantile(x, c(0.05, 0.95))可以求x的样本分位数。
rank(x)对x求秩得分（即名次，但从最小到最大排列）。

24.7.3 分类变量概括

分类变量一般输入为因子。
对因子或其它向量x，
table(x)返回x的每个不同值的频率（出现次数），
结果为一个类（class）为table的一维数组。
每个元素有对应的元素名，为x的各水平值。
如

res <- table(d.cancer[["sex"]]); res

## 
##  F  M 
## 13 21

##  F 
## 13

对单个分类变量, table结果是一个有元素名的向量。
用as.data.frame()函数把table的结果转为数据框:

##   Var1 Freq
## 1    F   13
## 2    M   21

用prop.table()将频数转换成百分比：

## 
##         F         M 
## 0.3823529 0.6176471

table作的单变量频数表可以用barplot表现为图形，如:

barplot(res, main="性别分布")

对两个分类变量x1和x2，
其每个组合的出现次数可以用table(x1,x2)函数统计， 结果叫做列联表。
如

res2 <- with(d.cancer, table(sex, type)); res2

##    type
## sex 鳞癌 腺癌
##   F    4    9
##   M   18    3

结果是一个类为table的二维数组（矩阵），
每行以第一个变量x1的各水平值为行名，
每列以第二个变量x2的各水平值为列名。
这里用了with()函数引入一个数据框，
后续的参数中的表达式可以直接使用数据框的变量。

对两个分类变量, table结果是一个矩阵。
用as.data.frame函数把table的结果转为数据框:

##   sex type Freq
## 1   F 鳞癌    4
## 2   M 鳞癌   18
## 3   F 腺癌    9
## 4   M 腺癌    3

列联表的结果可以用条形图表示。如

barplot(res2, legend=TRUE)

或

barplot(res2, legend=TRUE, beside=TRUE)

对于table()的结果列联表，
可以用addmargins()函数增加行和与列和。 如

##      type
## sex   鳞癌 腺癌 Sum
##   F      4    9  13
##   M     18    3  21
##   Sum   22   12  34

用margin.table()可以计算列联表行或列的和并返回，如

## sex
##  F  M 
## 13 21

## type
## 鳞癌 腺癌 
##   22   12

用prop.table(r)把一个列联表r转换成百分比表。 如

##    type
## sex       鳞癌       腺癌
##   F 0.11764706 0.26470588
##   M 0.52941176 0.08823529

用prop.table(res,1)把列联表res转换成行百分比表。
用prop.table(res,2)把列联表res转换成列百分比表。 如

##    type
## sex      鳞癌      腺癌
##   F 0.3076923 0.6923077
##   M 0.8571429 0.1428571

##    type
## sex      鳞癌      腺癌
##   F 0.1818182 0.7500000
##   M 0.8181818 0.2500000

在有多个分类变量时，
用as.data.frame(table(x1, x2, x3))
形成多个分类变量交叉分类的频数统计数据框。

dplyr包的count()功能与table()类似。

24.7.4 数据框概括

用colMeans()对数据框或矩阵的每列计算均值，
用colSums()对数据框或矩阵的每列计算总和。
用rowMeans()和rowSums()对矩阵的每行计算均值或总和。

数据框与矩阵有区别，
某些适用于矩阵的计算对数据框不适用，
例如矩阵乘法。
用as.matrix()把数据框的数值子集转换成矩阵。

对矩阵，用apply(x, 1, FUN)对矩阵x的每一行使用函数FUN计算结果，
用apply(x, 2, FUN)对矩阵x的每一列使用函数FUN计算结果。

如果apply(x,1,FUN)中的FUN对每个行变量得到多个m结果，
结果将是一个矩阵，行数为m，列数等于nrow(x)。
如果apply(x,2,FUN)中的FUN对每个列变量得到多个m结果，
结果将是一个矩阵，行数为m，列数等于ncol(x)。
例如：

apply(as.matrix(iris[,1:4]), 2,
        function(x)
        c(n=sum(!is.na(x)),
          mean=mean(x, na.rm=TRUE),
          sd=sd(x, na.rm=TRUE)))

##      Sepal.Length Sepal.Width Petal.Length Petal.Width
## n     150.0000000 150.0000000   150.000000 150.0000000
## mean    5.8433333   3.0573333     3.758000   1.1993333
## sd      0.8280661   0.4358663     1.765298   0.7622377

24.8.1 用tapply()分组汇总

用tapply()函数进行分组概括, 格式为：

其中X是一个向量， INDEX是一个分类变量， FUN是概括统计函数。

比如，下面的程序分性别组计算疗前体积的均值：

with(
  d.cancer, 
  tapply(v0, sex, mean))

##        F        M 
## 113.2354 108.1214

24.8.2 用aggregate()分组概括数据框

aggregate函数对输入的数据框用指定的分组变量（或交叉分组）
分组进行概括统计。
例如，下面的程序按性别分组计算年龄、疗前体积、疗后体积的平均值:

aggregate(
  d.cancer[, c("age", "v0", "v1")],
  list(sex=d.cancer[["sex"]]), 
  mean, na.rm=TRUE)

##   sex      age       v0       v1
## 1   F 66.14286 113.2354 42.65538
## 2   M 63.25000 108.1214 45.95524

aggregate()第一个参数是数据框，
第二个参数是列表，列表元素是用来分组或交叉分组的变量，
第三个参数是概括用的函数，
概括用的函数的选项可以在后面给出。

可以交叉分组后概括，如

with(
  d.cancer,
  aggregate(
    cbind(v0, v1), 
    list(sex=sex, type=type), 
    mean))

##   sex type        v0       v1
## 1   F 鳞癌 126.99250 45.54750
## 2   M 鳞癌 113.55722 49.65556
## 3   F 腺癌 107.12111 41.37000
## 4   M 腺癌  75.50667 23.75333

24.8.3 用split()函数分组后概括

split函数可以把数据框的各行按照一个或几个分组变量分为子集的列表，
然后可以用sapply()或vapply()对每组进行概括。
如

sp <- split(d.cancer[,c("v0","v1")], d.cancer[["sex"]])
sapply(sp, colMeans)

##            F         M
## v0 113.23538 108.12143
## v1  42.65538  45.95524

返回矩阵，行为变量v0, v1，列为不同性别，
值为相应的变量在各性别组的平均值。
当sapply()对列表每项的操作返回一个向量时，
总是列表每项的输出保存为结果的一列。
colMeans函数计算分组后数据框子集每列的平均值。

if("dplyr" %in% .packages()) detach("package:dplyr")
library(plyr)

ddply(d.cancer, "sex", 
      function(d) c(mean.v0 = mean(d[["v0"]], na.rm=TRUE)))

ddply(d.cancer, "sex", 
      function(d) colMeans(d[,c("v0", "v1")]))

f1 <- function(dsub){
  tab <- tibble(
    "变量"=c("v0", "v1"),
    "均值"=c(mean(dsub[,"v0"], na.rm=TRUE), 
           mean(dsub[,"v1"], na.rm=TRUE)),
    "标准差"=c(sd(dsub[,"v0"], na.rm=TRUE), 
            sd(dsub[,"v1"], na.rm=TRUE)))
  tab
}
ddply(d.cancer, "sex", f1)

f2 <- function(d, variables){
  d1 <- d[,variables,drop=FALSE]
  nnotmiss <- vapply(d1, function(x) sum(!is.na(x)), 1L)
  xm <- vapply(d1, mean, 0.0, na.rm=TRUE)
  xsd <- vapply(d1, sd, 1.0, na.rm=TRUE)
  xmed <- vapply(d1, median, 0.0, na.rm=TRUE)
  data.frame(variable=variables,
             n=nnotmiss,
             mean=xm,
             sd=xsd,
             median=xmed)
}
ddply(d.cancer, "sex", f2, variables = c("v0", "v1"))

ddply(d.cancer, c("sex", "type"), f2, variables=c("v0", "v1"))

f3 <- function(d){
  ff <- function(x){
    c(n=sum(!is.na(x)),
      each(mean, sd, median)(x, na.rm=TRUE))
  }
  ldply(Filter(is.numeric, d), ff)
}
ddply(d.cancer, "sex", f3)

f4 <- function(x){
  c(n=sum(!is.na(x)),
    each(mean, sd, median)(x, na.rm=TRUE))
}
ddply(d.cancer, "sex", 
      function(d) 
        ldply(Filter(is.numeric, d), f4))