title: Two Factor ANOVA Errorbar plot in R
body: <p>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.  </p>
  
  <p>So the question is... is there an easier/faster/simpler way to do all of this?</p>
  
  <pre><code>#   LIZARD LENGTH DATA
  island.1 &lt;- c(0.2, 5.9, 6.1, 6.5)
  island.2 &lt;- c(5.6, 14.8, 15.5, 16.4)
  island.3 &lt;- c(0.8, 3.9, 4.3, 4.9)
  sex.codes &lt;- c("Male", "Female", "Male", "Female")
  
  #   PUTTING DATA TOGETHER IN A DATA FRAME
  df.1 &lt;- data.frame(island.1, island.2, island.3, sex.codes)
  
  #   MELTING THE DATA FRAME INTO LONG FORM
  library(reshape)
  df.2 &lt;- melt(df.1)
  
  #   MEAN BY CELL
  mean.island1.male &lt;- with(df.2, mean(value[variable == "island.1" &amp; sex.codes == "Male"]))
  mean.island1.female &lt;- with(df.2, mean(value[variable == "island.1" &amp; sex.codes == "Female"]))
  mean.island2.male &lt;- with(df.2, mean(value[variable == "island.2" &amp; sex.codes == "Male"]))
  mean.island2.female &lt;- with(df.2, mean(value[variable == "island.2" &amp; sex.codes == "Female"]))
  mean.island3.male &lt;- with(df.2, mean(value[variable == "island.3" &amp; sex.codes == "Male"]))
  mean.island3.female &lt;- with(df.2, mean(value[variable == "island.3" &amp; sex.codes == "Female"]))
  
  #   ADDING CELL MEANS TO DATA FRAME
  df.2$means[df.2$variable == "island.1" &amp; df.2$sex.codes == "Male"] &lt;- mean.island1.male
  df.2$means[df.2$variable == "island.1" &amp; df.2$sex.codes == "Female"] &lt;- mean.island1.female
  df.2$means[df.2$variable == "island.2" &amp; df.2$sex.codes == "Male"] &lt;- mean.island2.male
  df.2$means[df.2$variable == "island.2" &amp; df.2$sex.codes == "Female"] &lt;- mean.island2.female
  df.2$means[df.2$variable == "island.3" &amp; df.2$sex.codes == "Male"] &lt;- mean.island3.male
  df.2$means[df.2$variable == "island.3" &amp; df.2$sex.codes == "Female"] &lt;- mean.island3.female
  
  #   LINEAR MODEL
  lizard.model &lt;- lm(value ~ variable*sex.codes, data=df.2)
  
  #   CALCULATING RESIDUALS BY HAND:
  df.2$residuals.1 &lt;- df.2$value - df.2$means
  
  #   CONFIRMING RESIDUALS FROM LINEAR MODEL:
  df.2$residuals.2 &lt;- residuals(lizard.model)
  
  #   TWO FACTOR MAIN EFFECT ANOVA
  lizard.anova &lt;- anova(lizard.model)        
  
  #   INTERACTION PLOT
  interaction.plot(df.2$variable, df.2$sex.codes, df.2$value)
  
  #   SAMPLE SIZE IN EACH CELL
  n &lt;- length(df.2$value[df.2$variable == "island.1" &amp; df.2$sex.codes == "Male"])
  # &gt; n
  # [1] 2
  
  #   NOTE: JUST FOR CLARITY, PRETEND n=10
  n &lt;- 10
  
  #   CALCULATING STANDARD ERROR
  island.se &lt;- sqrt(lizard.anova$M[4]/n)
  
  #   HALF CONFIDENCE INTERVAL
  island.ci.half &lt;- qt(0.95, lizard.anova$D[4]) * island.se
  
  #   MAKING SUMMARY DATA FRAME
  summary.df &lt;- 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)        
          )
  
  # &gt; 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()
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/sdgj2.png" alt="ggplot2 errobar plot of Two-way Main Effects Anova"></p>