Writing a function for summary statistics in R

Question

I'm having a problem I can't figure out... Basically I want to generate mean, SD, and N per group for a number of variables. My data looks like this:

dataSet <- data.frame(study_id=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4),
                      Timepoint=c(1,6,12,18,1,6,12,18,1,6,12,18,1,6,12,18),
                      Secretor=c(0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1),
                      Gene1=runif(16, min=0, max=100),
                      Gene2=runif(16, min=0, max=100),
                      Gene3=runif(16, min=0, max=100),
                      Gene4=runif(16, min=0, max=100))

Then I group it...

library(tidyverse)

grouped_dataSet <- dataSet %>%
  group_by(Secretor, Timepoint)

When I run the following line of code, I get what I want:

summarise(grouped_dataSet, mean = mean(Gene1, na.rm=T), sd = sd(Gene1, na.rm=T), n = n())

Output:

# A tibble: 8 x 5
# Groups:   Secretor [2]
  Secretor Timepoint  mean    sd     n
     <dbl>     <dbl> <dbl> <dbl> <int>
1        0         1  21.8 18.6      2
2        0         6  34.8 33.2      2
3        0        12  43.1  4.34     2
4        0        18  72.6 38.0      2
5        1         1  13.3 15.3      2
6        1         6  41.2 22.8      2
7        1        12  44.9 25.7      2
8        1        18  37.0  8.49     2

However, when I write this same line of code as a function (which I'm intending to then map onto many columns using tidyverse's purrr package), it doesn't work, instead returning "NA" for everything except the n column:

summary_function <- function(x) {
  summary <- summarise(grouped_dataSet, mean = mean(x, na.rm=T), sd = sd(x, na.rm=T), n = n())
  return(summary)
}

summary_function("Gene1")

Output:

# A tibble: 8 x 5
# Groups:   Secretor [2]
  Secretor Timepoint  mean    sd     n
     <dbl>     <dbl> <dbl> <dbl> <int>
1        0         1    NA    NA     2
2        0         6    NA    NA     2
3        0        12    NA    NA     2
4        0        18    NA    NA     2
5        1         1    NA    NA     2
6        1         6    NA    NA     2
7        1        12    NA    NA     2
8        1        18    NA    NA     2

This is the warning I get:

In var(if (is.vector(x) || is.factor(x)) x else as.double(x),  ... :
  NAs introduced by coercion

Could anyone please provide advice as to why it works as a line of code, but not as a function?

Many thanks in advance.

score 2 · Accepted Answer · answered Jun 25 '20 at 21:55

@akrun's suggestion for how to immediately solve your question is right on.

An alternative is to use the nesting functionality of tidyr by returning a single element list which contains a data.frame of your results.

summary_function <- function(x) {
  summary <- list(tibble(mean = mean(x, na.rm=T), sd = sd(x, na.rm=T), n = length(x[!is.na(x)])))
  return(summary)
}

Then you can use across to do the same function to multiple columns:

dataSet %>%
  group_by(Secretor, Timepoint) %>% 
  summarize(across(Gene1:Gene4, summary_function))
# A tibble: 8 x 6
# Groups:   Secretor [2]
#  Secretor Timepoint Gene1            Gene2            Gene3            Gene4           
#     <dbl>     <dbl> <list>           <list>           <list>           <list>          
#1        0         1 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#2        0         6 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#3        0        12 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#4        0        18 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#5        1         1 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#6        1         6 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#7        1        12 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>
#8        1        18 <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]> <tibble [1 × 3]>

Now we can unnest those same columns using unnest with names_sep =:

dataSet %>%
  group_by(Secretor, Timepoint) %>% 
  summarize(across(Gene1:Gene4, summary_function)) %>%
  unnest(Gene1:Gene4, names_sep = "_")
# A tibble: 8 x 14
# Groups:   Secretor [2]
#  Secretor Timepoint Gene1_mean Gene1_sd Gene1_n Gene2_mean Gene2_sd Gene2_n Gene3_mean Gene3_sd Gene3_n
#     <dbl>     <dbl>      <dbl>    <dbl>   <int>      <dbl>    <dbl>   <int>      <dbl>    <dbl>   <int>
#1        0         1      71.2     28.6        2       62.3     27.0       2       28.4    33.3        2
#2        0         6       5.40     7.43       2       58.6     29.1       2       37.0    33.9        2
#3        0        12      91.8     11.4        2       53.9     31.0       2       33.2    46.0        2
#4        0        18      51.5     65.0        2       65.3     40.2       2       63.8    32.7        2
#5        1         1      30.8     18.0        2       50.0     19.9       2       22.8     6.71       2
#6        1         6      63.9     49.2        2       59.9     41.8       2       30.9    39.5        2
#7        1        12      85.3      6.74       2       51.0     41.1       2       28.5    22.9        2
#8        1        18      41.7     44.8        2       80.2     24.0       2       64.7    17.4        2
## … with 3 more variables: Gene4_mean <dbl>, Gene4_sd <dbl>, Gene4_n <int>

This is a recent addition to tidyr and dplyr (version >=1.0.0), but can come handy.

Ah, OK. Nice to see a different solution. – Allan Cameron Jun 25 '20 at 21:59 — Allan Cameron, Jun 25 '20 at 21:59

score 1 · Answer 2 · answered Jun 25 '20 at 21:20

We can use ensym so that we can pass either quoted or unquoted and it can be evaluated (!!)

summary_function <- function(x) {
   x <- ensym(x)
    summarise(grouped_dataSet, 
        mean = mean(!! x, na.rm=T), sd = sd(!!x, na.rm=T), n = n())

  }

summary_function("Gene1")
# A tibble: 8 x 5
# Groups:   Secretor [2]
#  Secretor Timepoint  mean    sd     n
#     <dbl>     <dbl> <dbl> <dbl> <int>
#1        0         1 69.4   2.25     2
#2        0         6  9.67 13.6      2
#3        0        12 39.5  10.6      2
#4        0        18 17.4  19.2      2
#5        1         1 41.0  54.0      2
#6        1         6 58.5   7.57     2
#7        1        12 75.5   1.42     2
#8        1        18 80.5  24.7      2


summary_function(Gene1)
# A tibble: 8 x 5
# Groups:   Secretor [2]
#  Secretor Timepoint  mean    sd     n
#     <dbl>     <dbl> <dbl> <dbl> <int>
#1        0         1 69.4   2.25     2
#2        0         6  9.67 13.6      2
#3        0        12 39.5  10.6      2
#4        0        18 17.4  19.2      2
#5        1         1 41.0  54.0      2
#6        1         6 58.5   7.57     2
#7        1        12 75.5   1.42     2
#8        1        18 80.5  24.7      2

Also, for reusability in different datasets, it may be better to have additional argument that takes the dataset object

Writing a function for summary statistics in R

2 Answers2