假设检验：方差分析

knitr::opts_chunk$set(echo = TRUE)

安装包

library(dplyr)
library(ggplot2)
library(foreign) # load .sav dataset in spss
library(knitr)
# library(car) #  leveneTest() [in car package] will be used
opts_template$set(fig.height = 4, fig.width = 4)

加载数据

为研究三种不同饲料对生猪体重增加（wyh）的影响，将生猪随机分成三组各喂养不同的饲料（sl），得到喂养后的体重。由于喂养前的体重不作为可控变量但却对结果会产生一定影响的假设，数据集中也包含了喂养前的体重（wyq），作为自身身体条件的指标。

my_data <- read.spss("data/生猪与饲料.sav") %>% 
  as.data.frame()
# Show a random sample
set.seed(1234)
sample_n(my_data, 5)

EDA

先看一下数据，

print(\'修正前:\')
summary(my_data)
my_data <- my_data %>% 
  mutate(SL = as.factor(SL))
print(\'修正后:\')
summary(my_data)

方差分析

假设条件的验证

协变量有作用？

为分析猪喂养前的体重是否能够作为covariate协变量，我们绘制前后散点图，看协变量是否对目标变量有一定线性关系。在方差分析中，这里的原假设时协变量对观测变量的线性影响不显著，也就是说这个x不会影响其它x对y产生的作用。

# append linear coeffient to raw data
result = data.frame(SL = 0, intercept=0, slope=0)
for(i in c(1:3)) {
  res <- lm(WYH ~ WYQ, my_data %>% filter(SL == i))
  tmp <- data.frame(SL = i,
                    intercept = as.numeric(res$coefficients[1]),
                    slope = as.numeric(res$coefficients[2]))
  result <- rbind(result, tmp)
}
result <- result[-1,]
result$SL <- as.factor(result$SL)

my_data <- my_data %>%
  left_join(result)

# plot
ggplot(my_data, aes(WYQ, WYH))+
  geom_point(aes(colour = SL), size = 3)+
  geom_smooth(aes(WYQ, WYH))+
  geom_abline(aes(intercept = intercept, slope = slope, colour = SL)) +
  labs(x = "喂养前体重", y = "喂养后体重增加")+
  scale_fill_discrete(name = "饲料")+
    theme_classic()

ggplot(my_data, aes(WYQ, WYH))+
  geom_point(aes(colour = SL), size = 3)+
  labs(x = "喂养前体重", y = "喂养后体重增加")+
  scale_fill_discrete(name = "饲料")+
  theme_classic()

ggplot(my_data, aes(SL, WYH))+
  geom_boxplot(aes(fill = SL))+
  labs(x = "饲料", y = "喂养后体重增加")+
  scale_fill_discrete(name = "饲料")
  # facet_wrap(~ SL)

可见三种饲料的间，生猪喂养前后的体重增加量均称为相似的线性关系。因此，可作为协变量。

# MANOVA test
res.man <- aov(WYH ~ SL + WYQ, data = my_data)
summary(res.man)

方差检验

为了说明数据，将生猪的单因素方差分析结果一同显示，进行比较。

# ANOVA test
res.aov <- aov(WYH ~ SL, data = my_data)
summary(res.aov)

可见不同饲料对体重增加结果是有显著影响的，那么具体是哪些组之间差异最明显呢？

两两比较

# Tukey multiple pairwise-comparisons
TukeyHSD(res.aov)

多重比较

# Multiple comparisons using multcomp package
library(multcomp)
summary(glht(res.aov, linfct = mcp(SL = "Tukey")))

正太？

# 1. Normality
plot(res.aov, 2)

# 2. Extract the residuals
aov_residuals <- residuals(object = res.aov )
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )

各组方差齐性？

# 1. Homogeneity of variances
plot(res.aov, 1)

# 2. Levene’s test, which is less sensitive to departures from normal distribution.
# leveneTest(WYH ~ SL, data = my_data)

选取文献

统计分析与SPSS应用（第三版）薛薇编撰
数据来源官网地址