R中的两因子ANOVA误差线图 [英] Two Factor ANOVA Errorbar plot in R

问题描述

我们正在为生物学生教授一个统计类，并试图将R用作计算和数据可视化平台。尽可能多的，我们想避免使用额外的软件包，并在R中做任何非常奇特的事情;课程的重点是统计数据，而不是编程。尽管如此，我们还没有找到一种在R中为双因素ANOVA设计生成误差线图的非常好的方法。我们使用ggplot2软件包来绘制图表，虽然它有一个生成95％CI错误条的内置stat_summary方法，但这些方法的计算方式可能并不总是正确的。下面，我手工检查ANOVA的代码，并手工计算95％CIs（根据总剩余方差估算标准误，而不仅仅是组内方差ggplot的汇总方法使用的标准误）。最后，实际上有一个情节。

所以问题是......有没有更容易/更快速/更简单的方法来做到这一切？

< （0.2,5.9,6.1,6.5）
岛.2 <-C（5.6,14.8） （0.8,3.9,4.3,4.9）
sex.codes< -c（男，女，男，女性）

＃将数据放在数据框中
df.1 < - data.frame（island.1，island.2，island.3，sex.codes）

＃将数据帧融入到较长的格式中
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& - （df.2，mean（value [variable ==island.2& sex.codes ==Male]））
mean.island2.female& - （df.2，mean（值[variable ==island.2& sex.codes ==Fema （df.2，mean（value [variable ==island.3& sex.codes ==Male]））
mean.island3.female < - with（df.2，mean（value [variable ==island.3& sex.codes ==Female ]））

＃将数据单元添加到数据帧
df.2 $表示[df.2 $ variable ==island.1& df.2 $ sex.codes ==Male]< - mean.island1.male
df.2 $表示[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 $表示[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 $表示[df.2 $ variable ==island.3& df.2 $ sex.codes ==Female]< - mean.island3.female

＃LINEAR MODEL
lizard.model< - lm（value_variable * sex。代码，数据= df.2）

＃用手计算剩余物：
df.2 $ residuals.1 < - df.2 $ value - df.2 $表示

＃确认来自线性模型的残差：
df.2 $ residuals.2< - 残差（lizard.model）

＃双因素主效应方差分析
lizard.anova< - anova（lizard.model）

＃INTERACTION PLOT
interaction.plot（df.2 $ variable，df.2 $ sex.codes，df.2 $ value ）

＃每个单元格的样本大小
n < - 长度（df.2 $ value [df.2 $ variable ==island.1& df.2 $ sex。 codes ==Male]）
＃> n
＃[1] 2

＃注意：只适用于清晰度，PRETEND n = 10
n < - 10

＃计算标准错误
island.se < - sqrt（lizard.anova $ M [4] / n）

＃半身信心区间
island.ci.half <-qt（0.95， lizard.anova $ D [4]）* island.se

＃制作汇总数据帧
summary.df< - data.frame（
Means = c（mean。 island1.male，
mean.island1.female，
mean.island2.male，
mean.island2.female，
mean.island3.male，
mean。
位置= c（island1，
island1，
island2，
island2，
island3，
island3），
性别= c（男性，b $ b女性，
男性，
女性，
男性，
女性），
CI.half = rep（island.ci.half，6）
）

＃> summary.df
＃表示位置性别CI.half
＃1 3.15 island1男性2.165215
＃2 6.20 island1女性2.165215
＃3 10.55 island2男性2.165215
＃4 15.60 island2女2.165215
＃5 2.55 island3男2.165215
＃6 4.40 island3女2.165215

＃生成ERRORBAR PLOT
库（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（）

解决方案

使用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））
  

We're teaching a stats class for biology students and trying to use R as the computing and data visualization platform. As much as possible, we'd like to avoid using extra packages and doing anything terribly "fancy" in R; the focus of the course is on the statistics, not the programming. Nevertheless, we haven't found a very good way of generating an errorbar plot in R for a two factor ANOVA design. We're using the ggplot2 package to make the plot, and while it does have a built-in stat_summary method of generating 95% CI errorbars, the way these are calculated may not always be the right way . Below, I go through the code for the ANOVA by hand and calculate the 95% CIs by hand also (with standard error estimated from the total residual variance, not just the within-group variance ggplot's summary method would use). At the end, there's actually a plot.

So the question is... is there an easier/faster/simpler way to do all of this?

#   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()

解决方案

Here is another attempt using the sciplot package. Alternative ways to compute the confidence intervals can be passed in parameter 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))

这篇关于R中的两因子ANOVA误差线图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！