Search code examples
rfunctionstatisticssummary

Writing a function for summary statistics in R


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.


Solution

  • @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.