雙因子變異數分析(Two-way ANOVA)

1
2
3

setwd("C:/R/biostatistics")
data <- read.csv("data_twoway_anova.csv", header=T)
data$temp <- as.factor(data$temp)

1
2
3
4

#boxplot
boxplot(gain ~ supp + temp, col=c("red", "orange"), ylab="Body weight gain", xlab="supplement and temperature", data=data)
libraray(ggplot2)
ggplot(data, aes(x=supp, y=gain, col=temp)) + geom_boxplot()

1
2
3
4
5

# class by supp and temp, and then count means and standard error
library(reshape); library(agricolae)
mystats <- function(x) {c(n=length(x), mean=mean(x, na.rm=T), sd=sd(x))}
dfm <- melt(data, measure.vars="gain", id.vars=c("supp", "temp"))
cast(dfm, supp + temp + variable ~ ., mystates)

1
2
3
4
5
6
7

   supp	temp	variable	n	mean	sd
1	A	20	gain	10	7.98	2.746634
2	A	25	gain	10	16.77	2.515309	
3	A	30	gain	10	26.14	4.797731	
4	B	20	gain	10	13.23	4.459709
5	B	25	gain	10	22.70	3.910953
6	B	30	gain	10	26.06	2.655058

1
2
3
4
5
6

# main effect ANOVA with no interaction
model_aov <- aov(data=data, gain ~ supp + temp)
summary(model_aov)
model_aov$cofficients  
# if we don't care supp and temp, mean of gain is 12.455 
# if suppOJ = 0. and then suppVC = -3.7