在R中纳入重复样本缺失因子水平以计算均值与标准误

现有数据框df，包含站点（Site）、重复样本（Replica）、生境（Habitat）及对应占比（Porcentage），示例数据如下：

df <- data.frame(Site = c("A","A","A","A","A","B","B","B","B","B","B"),
                 Replica = c(1,1,1,2,2,1,1,1,2,2,2),
                 Habitat = c("X","Y","Z","X","Y","X","Y","M","X","M","Z"),
                 Porcentage = c(46,38,16,40,60,20,60,20,35,55,10))

每个站点内单个重复样本的所有生境占比之和为100。需要按站点计算各生境的平均占比（Mean_Porcentage）及均值标准误（SE），但直接分组汇总的结果未将缺失的生境按0计算，导致均值和标准误不准确：

library(dplyr)

df %>% 
  group_by(Site,Habitat) %>% 
  summarise(Mean_Porcentage = mean(Porcentage),
            SE = sd(Porcentage)/sqrt(length(Porcentage)))

输出结果：

# A tibble: 7 × 4
# Groups:   Site [2]
  Site  Habitat Mean_Porcentage    SE
  <chr> <chr>             <dbl> <dbl>
1 A     X                  43     3  
2 A     Y                  49    11  
3 A     Z                  16    NA  
4 B     M                  37.5  17.5
5 B     X                  27.5   7.5
6 B     Y                  60    NA  
7 B     Z                  10    NA 

例如站点A的生境Z在Replica=2中缺失，应按0计算，均值应为8、SE为8；站点A中未出现的生境M也需显示均值为0。期望结果如下：

Site Habitat Mean_Porcentage     SE
1    A       M             0.0      0
2    A       X            43.0      3
3    A       Y            49.0     11
4    A       Z             8.0      8
5    B       M            37.5   17.5
6    B       X            27.5    7.5
7    B       Y            30.0     30
8    B       Z             5.0      5

解决方案 #

核心思路是先构建每个站点-重复样本下所有生境的完整组合，将缺失的生境占比补0后再进行统计：

library(dplyr)
library(tidyr)

# 获取所有出现过的生境类型
all_habitats <- unique(df$Habitat)

# 构建完整的站点-重复样本-生境组合，补全缺失值为0后计算统计量
df_result <- df %>%
  # 按站点分组，生成每个站点下所有重复样本与所有生境的笛卡尔积
  group_by(Site) %>%
  expand(Replica, Habitat = all_habitats) %>%
  # 左连接原数据，将缺失的占比替换为0
  left_join(df, by = c("Site", "Replica", "Habitat")) %>%
  replace_na(list(Porcentage = 0)) %>%
  # 按站点和生境分组计算均值和标准误
  group_by(Site, Habitat) %>%
  summarise(
    Mean_Porcentage = mean(Porcentage),
    # 处理所有值为0的情况，避免sd返回NA，此时SE为0
    SE = ifelse(all(Porcentage == 0), 0, sd(Porcentage)/sqrt(n()))
  ) %>%
  ungroup()

# 查看结果
df_result

执行后即可得到符合预期的输出。

内容的提问来源于stack exchange，提问作者Dekike