R 中的双因素方差分析误差条图

Question

我们正在为生物学学生教授统计课程，并尝试使用 R 作为计算和数据可视化平台。 我们希望尽可能避免使用额外的包并在 R 中做任何非常“花哨”的事情；该课程的重点是统计，而不是编程。 尽管如此，我们还没有找到在 R 中为二因素方差分析设计生成误差条图的非常好的方法。 我们使用 ggplot2 包来制作绘图，虽然它确实有一个内置的 stat_summary 方法来生成 95% CI 误差条，但这些计算方式可能并不总是正确的。 下面，我手动检查方差分析的代码，并手动计算 95% CI（根据总残差方差估计标准误差，而不仅仅是 ggplot 的汇总方法使用的组内方差）。 到最后，其实还有一个剧情。

所以问题是...有没有一种更容易/更快/更简单的方法来完成这一切？

#   LIZARD LENGTH DATA
island.1 <- c(0.2, 5.9, 6.1, 6.5)
island.2 <- c(5.6, 14.8, 15.5, 16.4)
island.3 <- c(0.8, 3.9, 4.3, 4.9)
sex.codes <- c("Male", "Female", "Male", "Female")

#   PUTTING DATA TOGETHER IN A DATA FRAME
df.1 <- data.frame(island.1, island.2, island.3, sex.codes)

#   MELTING THE DATA FRAME INTO LONG FORM
library(reshape)
df.2 <- melt(df.1)

#   MEAN BY CELL
mean.island1.male <- with(df.2, mean(value[variable == "island.1" & sex.codes == "Male"]))
mean.island1.female <- with(df.2, mean(value[variable == "island.1" & sex.codes == "Female"]))
mean.island2.male <- with(df.2, mean(value[variable == "island.2" & sex.codes == "Male"]))
mean.island2.female <- with(df.2, mean(value[variable == "island.2" & sex.codes == "Female"]))
mean.island3.male <- with(df.2, mean(value[variable == "island.3" & sex.codes == "Male"]))
mean.island3.female <- with(df.2, mean(value[variable == "island.3" & sex.codes == "Female"]))

#   ADDING CELL MEANS TO DATA FRAME
df.2$means[df.2$variable == "island.1" & df.2$sex.codes == "Male"] <- mean.island1.male
df.2$means[df.2$variable == "island.1" & df.2$sex.codes == "Female"] <- mean.island1.female
df.2$means[df.2$variable == "island.2" & df.2$sex.codes == "Male"] <- mean.island2.male
df.2$means[df.2$variable == "island.2" & df.2$sex.codes == "Female"] <- mean.island2.female
df.2$means[df.2$variable == "island.3" & df.2$sex.codes == "Male"] <- mean.island3.male
df.2$means[df.2$variable == "island.3" & df.2$sex.codes == "Female"] <- mean.island3.female

#   LINEAR MODEL
lizard.model <- lm(value ~ variable*sex.codes, data=df.2)

#   CALCULATING RESIDUALS BY HAND:
df.2$residuals.1 <- df.2$value - df.2$means

#   CONFIRMING RESIDUALS FROM LINEAR MODEL:
df.2$residuals.2 <- residuals(lizard.model)

#   TWO FACTOR MAIN EFFECT ANOVA
lizard.anova <- anova(lizard.model)        

#   INTERACTION PLOT
interaction.plot(df.2$variable, df.2$sex.codes, df.2$value)

#   SAMPLE SIZE IN EACH CELL
n <- length(df.2$value[df.2$variable == "island.1" & df.2$sex.codes == "Male"])
# > n
# [1] 2

#   NOTE: JUST FOR CLARITY, PRETEND n=10
n <- 10

#   CALCULATING STANDARD ERROR
island.se <- sqrt(lizard.anova$M[4]/n)

#   HALF CONFIDENCE INTERVAL
island.ci.half <- qt(0.95, lizard.anova$D[4]) * island.se

#   MAKING SUMMARY DATA FRAME
summary.df <- data.frame(
        Means = c(mean.island1.male,
                mean.island1.female,
                mean.island2.male,
                mean.island2.female,
                mean.island3.male,
                mean.island3.female),
        Location = c("island1",
                "island1",
                "island2",
                "island2",
                "island3",
                "island3"),
        Sex = c("male",
                "female",
                "male",
                "female",
                "male",
                "female"),      
        CI.half = rep(island.ci.half, 6)        
        )

# > summary.df
# Means Location    Sex  CI.half
# 1  3.15  island1   male 2.165215
# 2  6.20  island1 female 2.165215
# 3 10.55  island2   male 2.165215
# 4 15.60  island2 female 2.165215
# 5  2.55  island3   male 2.165215
# 6  4.40  island3 female 2.165215

#   GENERATING THE ERRORBAR PLOT
library(ggplot2)

qplot(data=summary.df,
        y=Means,
        x=Location,
        group=Sex,
        ymin=Means-CI.half,
        ymax=Means+CI.half,
        geom=c("point", "errorbar", "line"),
        color=Sex,
        shape=Sex,
        width=0.25) + theme_bw()

Answer 1

这是使用 sciplot 包的另一次尝试。计算置信区间的替代方法可以在参数 ci.fun 中传递。

lineplot.CI(variable,value, group =sex.codes , data = df.2, cex = 1.5,
            xlab = "Location", ylab = "means", cex.lab = 1.2, x.leg = 1,
            col = c("blue","red"), pch = c(16,16))

Answer 2

我不得不承认我对你的代码感到非常困惑。不要将此视为个人批评，但我强烈建议您学习您的学生尽可能多地使用 R 的力量。他们只能从中受益，而我的经验是，如果我不向他们灌输一行又一行的混乱代码，他们会更快地理解正在发生的事情。

首先，你不必手工计算平均值。只要这样做：

df.2$mean <- with(df.2,ave(value,sex.codes,variable,FUN=mean))

另请参阅

?ave

。这比您的示例中混乱的代码更清晰。如果你有lizard.model，你就可以使用

fitted(lizard.model)

并将这些值与平均值进行比较。

那么我强烈不同意你的观点。您计算的并不是您预测的标准误差。要正确执行此操作，请使用

predict()

函数

outcome <- predict(lizard.model,se.fit=TRUE)
df.2$CI.half <- outcome$se / 2

要获得预测平均值的置信区间，如果您希望学生正确理解这一点，则必须使用正确的公式。请看一下来自 Faraway 的使用 R 的令人难以置信的伟大实用回归和方差分析的第 3.5 节。它包含大量代码示例，其中所有内容都是以方便而简洁的方式手动计算的。它将为您和您的学生服务。我从中学到了很多东西，并且在向学生解释这些事情时经常使用它作为指导。

现在要获取摘要数据框，您有几个选项，但这个选项有效并且很容易理解。

summary.df <- unique(df.2[,-c(3,5,6)])
names(summary.df) <- c('Sex','Location','Means','CI.half')

现在您可以直接运行您的绘图代码。

或者，如果您想要值的预测误差，您可以使用以下内容：

lizard.predict <- predict(lizard.model,interval='prediction')
df.2$lower <- lizard.predict[,2]
df.2$upper <- lizard.predict[,3]

summary.df <- unique(df.2[,-3])
names(summary.df)[1:3] <- c('Sex','Location','Means')


qplot(data=summary.df,
        y=Means,
        x=Location,
        group=Sex,
        ymin=lower,
        ymax=upper,
        geom=c("point", "errorbar", "line"),
        color=Sex,
        shape=Sex,
        width=0.25) + theme_bw()

PS：如果我在这里或那里听起来很刺耳，那不是故意的。英语不是我的母语，我仍然不熟悉这种语言的微妙之处。

Answer 3

[潜在的无耻促销]您应该考虑 HandyStuff 包中的compareCats 和 rxnNorm 函数，可在 www.github.com/bryanhanson/HandyStuff 警告：我不确定它是否可以与 R 2.14 无缝配合。 特别是， rxnNorm 看起来像您正在尝试生成的图，而且它在总结统计数据和装饰图方面为您提供了多种选项。 但是，这需要让您的学生安装一个单独的软件包，所以也许您会排除它（但它允许学生专注于呈现和分析数据）。 此处包含 ?rxnNorm 示例的绘图。

使用 rxnNorm，您可以选择多种计算 CI 的方法，由参数“方法”控制。 以下是实际功能（来自 ChemoSpec 包）。

> seX <- function (x)  sd(x, na.rm = TRUE)/sqrt(length(na.omit(x)))
> <environment: namespace:ChemoSpec>
> 
> seXy <- function (x)  {
>     m <- mean(na.omit(x))
>     se <- seX(x)
>     u <- m + se
>     l <- m - se
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>
> 
> 
> seXy95 <- function (x)  {
>     m <- mean(na.omit(x))
>     se <- seX(x)
>     u <- m + 1.96 * se
>     l <- m - 1.96 * se
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>
> 
> 
> seXyIqr <- function (x)  {
>     i <- fivenum(x)
>     c(y = i[3], ymin = i[2], ymax = i[4]) } <environment: namespace:ChemoSpec>
> 
> seXyMad <- function (x)  {
>     m <- median(na.omit(x))
>     d <- mad(na.omit(x))
>     u <- m + d
>     l <- m - d
>     c(y = m, ymin = l, ymax = u) } <environment: namespace:ChemoSpec>

Answer 4

使用该库可以得到一个简单的均值图以及 95% 置信区间superb（请注意，我是该库的维护者）。

在您的代码中，您提供了 data.frame

df.1

，所以我从那里开始

library(superb)

superb( crange(island.1, island.3) ~ sex.codes, df.1, 
WSFactors="island(3)", plotStyle="line")

使用

crange()

，您可以提供一系列列来进行说明（使用

cbind()

逐一命名列）。

您可以使用任何附加的 ggplot 指令来个性化绘图。例如，与

superb( crange(island.1, island.3) ~ sex.codes, df.1, 
  WSFactors="island(3)", plotStyle="line") + 
ylab("Means") + theme_bw()

你得到