如何在R的ggplot绘图中用标准误(SE)替代标准差(SD)

问题背景 #

我在R Studio中使用ggplot和plyr库编写代码，原本可以计算均值、标准差(SD)并绘制带SD误差棒的柱状图，但尝试将误差棒替换为标准误(SE)时，出现报错：

Error in se(x[[col]], na.rm = TRUE) : could not find function "se"

不清楚如何为数据框计算SE并应用到绘图中，相关代码如下：

原均值&SD计算函数 #

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}

原柱状图绘制代码 #

ggplot(datossummary, aes(x=Treatment, y=Phosphatase, fill=Treatment)) +
  geom_bar(stat="identity",color="black", width=0.4) +  
  scale_fill_manual(values=c("lightpink2", "lightpink2", "lightblue2", "lightblue2", "lightsalmon", "lightsalmon", "lightcyan4", "lightcyan4", "mediumpurple1", "mediumpurple1")) + 
  theme_light() +
  ggtitle("Phosphatase") +
  xlab("Treatment") +
  ylab("CE") +
  geom_errorbar(aes(ymin=Phosphatase-sd, ymax=Phosphatase+sd), width=.2, position=position_dodge(.9))

尝试替换SE的报错代码 #

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      se = se(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
} 

解决方案 #

R中没有内置的se()函数，标准误的计算公式为：SE = SD / √n（n为去除NA后的有效样本量），按以下步骤修改代码即可：

1. 自定义标准误计算函数 #

先定义一个计算SE的函数，你可以把它放在全局环境，或者嵌入到数据汇总函数中：

se <- function(x, na.rm = TRUE) {
  # 先计算去除NA后的标准差，再除以有效样本量的平方根
  sd(x, na.rm = na.rm) / sqrt(length(na.omit(x)))
}

2. 修改数据汇总函数 #

将原函数中的SD计算替换为SE计算，确保调用自定义的se()函数：

data_summary <- function(data, varname, groupnames){
  require(plyr)
  # 若未在全局定义se函数，可在此嵌入定义
  se <- function(x, na.rm = TRUE) {
    sd(x, na.rm = na.rm) / sqrt(length(na.omit(x)))
  }
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      se = se(x[[col]], na.rm=TRUE))
  }
  data_sum <- ddply(data, groupnames, .fun=summary_func, varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}

3. 修改绘图代码 #

将误差棒的计算逻辑从sd改为se：

# 先生成带SE的汇总数据（替换成你的原始数据框名称）
datossummary <- data_summary(your_raw_data, "Phosphatase", "Treatment")

# 绘制带SE误差棒的柱状图
ggplot(datossummary, aes(x=Treatment, y=Phosphatase, fill=Treatment)) +
  geom_bar(stat="identity", color="black", width=0.4) +
  scale_fill_manual(values=c("lightpink2", "lightpink2", "lightblue2", "lightblue2", "lightsalmon", "lightsalmon", "lightcyan4", "lightcyan4", "mediumpurple1", "mediumpurple1")) +
  theme_light() +
  ggtitle("Phosphatase") +
  xlab("Treatment") +
  ylab("CE") +
  geom_errorbar(aes(ymin=Phosphatase - se, ymax=Phosphatase + se), width=.2, position=position_dodge(.9))

内容的提问来源于stack exchange，提问作者Celia Garcí­a Dí­az