group_by() 和 summarise() 计算 sd 但返回 NA 作为平均值

Question

我正在尝试从按时间和某些变量分组的数据框中计算平均值和标准差。

当我用一个虚构的数据框执行以下代码时没有问题：

library(tidyverse)
library(reshape2)

experiment <- c("ex1","ex1","ex1","ex1","ex1","ex1","ex1","ex1","ex2","ex2","ex2","ex2","ex2","ex2","ex2","ex2")
time <- c("1","1","1","1","2","2","2","2","1","1","1","1","2","2","2","2")  
bottle <- c("1","2","3","4","1","2","3","4","1","2","3","4","1","2","3","4")
var1 <- runif(16, min=2, max=50)
var2 <- runif(16, min=2, max=50)
var3 <- runif(16, min=2, max=50)
df <- data.frame(experiment,time,bottle,var1,var2,var3)

df_melt <- melt(df, id = c("experiment", "time","bottle"))

df_avg <- df_melt %>% 
  filter(experiment %in% c("ex2")) %>%
  filter(bottle %in% c(1,2,3)) %>% 
  filter(variable %in% c("var1","var2")) %>%  
  group_by(time, variable) %>%
  summarise(mean = mean(value), sd= sd(value)) 

上面给了我mean和sd没问题。

但是，当我执行完全相同的代码但使用我的真实数据时，它会计算标准偏差但返回所有 NA 作为平均值。我知道在 Stack Overflow 上不鼓励共享数据文件，但我真的在这里旋转。

这是我的data file，代码如下：

data <- read.csv("path to data file",fileEncoding="UTF-8-BOM")
data_melted <- melt(data, id = c("Experiment", "Time", "Bottle"))

data_avg <- data_melted %>% 
  filter(Experiment %in% c("Vienn_U-A")) %>%
  filter(Bottle %in% c(1,2,4,5)) %>% 
  filter(variable %in% c("NH4","NO2","Urea")) %>%  
  group_by(Time, variable) %>%
  summarise(mean = mean(value), sd= sd(value)) 

这是输出的样子：

> data_avg
# A tibble: 54 x 4
# Groups:   Time [18]
    Time variable  mean     sd
   <dbl> <fct>    <dbl>  <dbl>
 1   0   NH4         NA  0.750
 2   0   Urea        NA 10.4  
 3   0   NO2         NA  0.243
 4  24.6 NH4         NA  0.539
 5  24.6 Urea        NA  8.87 
 6  24.6 NO2         NA  0.590
 7  48.2 NH4         NA  0.477
 8  48.2 Urea        NA  8.97 
 9  48.2 NO2         NA  1.21 
10  77.9 NH4         NA  0.899
# ... with 44 more rows

为什么它计算 sd，但返回 NA 表示平均值？任何指导将不胜感激！

Answer 1

受 akrun 的回答启发，我能够通过使用 .xlsx 文件并避免使用 .csv 文件来解决此问题。拥有NA 值不是问题。

这是data file in .xlsx format，下面的代码有效：

data <- read_excel("CSV for plotting.xlsx") #following akrun's suggestion 
data_melted <- melt(data, id = c("Experiment", "Time","Bottle"))
data_avg <- data_melted %>% 
  filter(Experiment %in% c("Vienn_U-A")) %>%
  filter(Bottle %in% c(1,2,4,5)) %>% 
  filter(variable %in% c("NH4","NO2","Urea")) %>%  
  group_by(Time, variable) %>%
  summarise(mean = mean(value), sd= sd(value)) 

使用正确的mean 值输出：

    Time variable   mean     sd
   <dbl> <fct>     <dbl>  <dbl>
 1   0   NH4        4.77  0.750
 2   0   Urea     281.   10.4  
 3   0   NO2       12.3   0.243
 4  24.6 NH4        9.11  0.539
 5  24.6 Urea     275.    8.87 
 6  24.6 NO2       17.8   0.590
 7  48.2 NH4       13.6   0.477
 8  48.2 Urea     268.    8.97 
 9  48.2 NO2       27.9   1.21 
10  77.9 NH4        8.75  0.899
# ... with 44 more rows

如上所述，使用 .csv 文件仍然会出现为 mean 返回 NA 的问题。

这里也是data file but in .csv format

data <- read.csv("CSV for plotting.csv",fileEncoding="UTF-8-BOM")
data_melted <- melt(data, id = c("Experiment", "Time","Bottle"))
data_avg <- data_melted %>% 
  filter(Experiment %in% c("Vienn_U-A")) %>%
  filter(Bottle %in% c(1,2,4,5)) %>% 
  filter(variable %in% c("NH4","NO2","Urea")) %>%  
  group_by(Time, variable) %>%
  summarise(mean = mean(value), sd= sd(value)) 

NA 输出 mean 但不是 sd：

 Time variable  mean     sd
   <dbl> <fct>    <dbl>  <dbl>
 1   0   NH4         NA  0.750
 2   0   Urea        NA 10.4  
 3   0   NO2         NA  0.243
 4  24.6 NH4         NA  0.539
 5  24.6 Urea        NA  8.87 
 6  24.6 NO2         NA  0.590
 7  48.2 NH4         NA  0.477
 8  48.2 Urea        NA  8.97 
 9  48.2 NO2         NA  1.21 
10  77.9 NH4         NA  0.899
# ... with 44 more rows

为什么.csv文件有这个问题我仍然很困惑，但至少有一个解决方案。

Answer 2

根据新链接，问题是在使用 read.csv 而不是在 read_excel 中读取时如何更改其中一种列类型，即列 'TN' 是 read.csv 中的 character 类，而它是 @ 987654325@

> data <- read.csv("CSV for plotting.csv", na.strings = "")
> str(data)
'data.frame':   292 obs. of  9 variables:
 $ Experiment: chr  "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" ...
 $ Time      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Bottle    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ NH4       : num  550 551 546 555 551 ...
 $ Urea      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ NO2       : num  18.6 18.3 18.6 18.1 18.9 ...
 $ NO3       : logi  NA NA NA NA NA NA ...
 $ Ln.NO2    : num  2.92 2.91 2.92 2.9 2.94 ...
 $ TN        : chr  "568.5714286" "569" "564.2857143" "573.1428571" ...

这部分是由于某些元素被读取为#REF

[81] "913.8095238" "870.1428571" "914.1190476" "946.9761905" "903.5238095" "962.6428571" "965.2619048" "982.3809524" "909.1904762" "#REF!"     

使用read_excel阅读时，同一行是NA

[85]  903.5238  962.6429  965.2619  982.3810  909.1905        NA  

这会导致melt 出现重大问题，因为指定的列是"Experiment", "Time","Bottle"，这意味着当melted 时的另一列将具有单一类型，但因为它在一个类中是character，所以整个列转换为character 类

> str(data_melted)
'data.frame':   1752 obs. of  5 variables:
 $ Experiment: chr  "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" ...
 $ Time      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Bottle    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ variable  : Factor w/ 6 levels "NH4","Urea","NO2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value     : chr  "550" "550.7142857" "545.7142857" "555" ...

即value 列是这里感兴趣的列，它现在是与 mean 预期不同的类，这对 sd 来说不是问题，因为它被转换了

> sd(c("15", "25", "35"))
[1] 10
>  sd(as.numeric(c("15", "25", "35"))) # same output
[1] 10
> mean(c("15", "25", "35"))
[1] NA

原因是根据?sd的文档

x- 数值向量或 R 对象，但不是可通过 as.double(x) 强制转换为数值的因子。

即它正在由as.double 内部转换并且函数确认它

> sd
function (x, na.rm = FALSE) 
sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x), 
    na.rm = na.rm))

melt 进行类型转换而不会引发错误。假设，如果我们使用pivot_longer（来自tidyr）会引发错误，我们就会知道问题所在

> tidyr::pivot_longer(data, cols = NH4:TN)
Error: Can't combine `NH4` <double> and `TN` <character>.
Run `rlang::last_error()` to see where the error occurred.

---

在使用read.csv 阅读的上下文中，可以选择将#REF 指定为NA 元素

> data <- read.csv("CSV for plotting.csv", na.strings = c("", "#REF!"))
> str(data)
'data.frame':   292 obs. of  9 variables:
 $ Experiment: chr  "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" ...
 $ Time      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Bottle    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ NH4       : num  550 551 546 555 551 ...
 $ Urea      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ NO2       : num  18.6 18.3 18.6 18.1 18.9 ...
 $ NO3       : logi  NA NA NA NA NA NA ...
 $ Ln.NO2    : num  2.92 2.91 2.92 2.9 2.94 ...
 $ TN        : num  569 569 564 573 570 ...

现在，melt 可以正常工作了

> data_melted <- melt(data, id = c("Experiment", "Time","Bottle"))
> str(data_melted)
'data.frame':   1752 obs. of  5 variables:
 $ Experiment: chr  "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" "Vienn_A-U" ...
 $ Time      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Bottle    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ variable  : Factor w/ 6 levels "NH4","Urea","NO2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value     : num  550 551 546 555 551 ...

mean 也是如此

> data_avg <- data_melted %>% 
+   filter(Experiment %in% c("Vienn_U-A")) %>%
+   filter(Bottle %in% c(1,2,4,5)) %>% 
+   filter(variable %in% c("NH4","NO2","Urea")) %>%  
+   group_by(Time, variable) %>%
+   summarise(mean = mean(value), sd= sd(value)) 
`summarise()` has grouped output by 'Time'. You can override using the `.groups` argument.
> data_avg
# A tibble: 54 × 4
# Groups:   Time [18]
    Time variable   mean     sd
   <dbl> <fct>     <dbl>  <dbl>
 1   0   NH4        4.77  0.750
 2   0   Urea     281.   10.4  
 3   0   NO2       12.3   0.243
 4  24.6 NH4        9.11  0.539
 5  24.6 Urea     275.    8.87 
 6  24.6 NO2       17.8   0.590
 7  48.2 NH4       13.6   0.477
 8  48.2 Urea     268.    8.97 
 9  48.2 NO2       27.9   1.21 
10  77.9 NH4        8.75  0.899
# … with 44 more rows

---

如果有NA 值，请考虑在mean 和sd 中使用na.rm = TRUE。默认情况下，这些设置为FALSE。如果base 中的mean 被屏蔽，请使用::。此外，原始链接似乎是 excel 文件。所以，用read_excel读取文件

library(dplyr)
data_avg <- data_melted %>% 
  filter(Experiment %in% c("Vienn_U-A")) %>%
  filter(Bottle %in% c(1,2,4,5)) %>% 
  filter(variable %in% c("NH4","NO2","Urea")) %>%  
  group_by(Time, variable) %>%
  summarise(mean = base::mean(value, na.rm = TRUE), sd= sd(value, na.rm = TRUE))

-输出

> data_avg
# A tibble: 54 × 4
# Groups:   Time [18]
    Time variable   mean     sd
   <dbl> <fct>     <dbl>  <dbl>
 1   0   NH4        4.77  0.750
 2   0   Urea     281.   10.4  
 3   0   NO2       12.3   0.243
 4  24.6 NH4        9.11  0.539
 5  24.6 Urea     275.    8.87 
 6  24.6 NO2       17.8   0.590
 7  48.2 NH4       13.6   0.477
 8  48.2 Urea     268.    8.97 
 9  48.2 NO2       27.9   1.21 
10  77.9 NH4        8.75  0.899
# … with 44 more rows

---

链接中的文件显示excel 数据。使用读取excel数据的包阅读可能会更好

library(readxl)
data <- read_excel("CSV for plotting.xlsx")