04方差分析

1.One way ANOVA

Fisher，第三代弟子，很厉害。实验设计，大部分是方差分析的内容。

用分类变量来解释连续变量的一个情况。

t检验比较两个，同时方差分析。

方差分析可以分析多个的因素。

下面这张图是指的自由度

第三种是最好，组间差异最大，最

看一看几个分布是否存在差异。

这个事实上是组间差异除以error MS

F分布的情况，x表示的是值。

# F distribution
df.1 <- 10
df.2 <- 1
f.d <- df (seq (0.1, 15, by=0.1), df.1, df.2)
plot (f.d, xlab='X', ylab='p', type='l',
xlim=c(0,15), ylim=c(0,1))
for (df.2 in 1:10){
f.d <- df (seq (0.1, 15, by=0.1), 10, df.2)
lines(f.d, type='l', col=rainbow(10)[df.2])
legend(df.2*1.3-1, 1, paste('df=', df.2, sep=''),
text.col = rainbow(10)[df.2], box.lty=0, cex=1)
}

卡方分布具有相加性，多个卡方分布，想加，自由度。

各类分布的区别

F、Z、T和chi分布的区别：F分布是两个数的平方。Z分布就是相当与一个正态分布。T分布则是一个平方。chi分布像是一个单变量的平方的分布。

上文的参考文献： (https://blog.csdn.net/anshuai_aw1/article/details/82735201)

(https://m.ximalaya.com/ask/q13347376)

方差分析的内容

# One way ANOVA
# Rodent weight at 7 sites
site1 <- c(9.4, 8.7, 13.3, 13.6, 15, 15.2, 17.7, 18.6, 22.2)
site2 <- c(16.8, 30.8, 33.6, 40.5, 48.9)
site3 <- c(27.0, 28.9, 32, 32.7, 35.5, 45.6)
site4 <- c(21.0, 23.4, 27.5, 27.5, 30.5, 31.9, 32.5, 33.8, 33.8)
site5 <- c(24.3, 29.7, 19.9)
site6 <- c(17.7, 19.7, 21.5, 27.9, 34.8, 40.2)
site7 <- c(16.5, 20.7, 23.5, 26.4, 26.7, 29.5, 29.8, 31.9, 35.5)
rodent.survey <- data.frame(weight=c(site1,site2,site3,site4,site5,site6,site7),
site=factor(c(rep("1",9),rep("2",5),rep("3",6),rep("4",9),rep("5",3),
rep("6",6),rep("7",9))))
options(digits=3) # default value = 7
tapply(rodent.survey$weight, rodent.survey$site, mean)
tapply(rodent.survey$weight, rodent.survey$site, var)
tapply(rodent.survey$weight, rodent.survey$site, sum)
boxplot(weight~site, data=rodent.survey, xlab='Sites', ylab='Weight')

输入数据和检查数据

tapply(rodent.survey$weight, rodent.survey$site, mean)

tapply(rodent.survey$weight, rodent.survey$site, var)

tapply(rodent.survey$weight, rodent.survey$site, sum)

# One Way Anova (Completely Randomized Design)
fit <- aov (weight~site, data=rodent.survey)
summary(fit)
# fit=lm(weight~site, data=rodent.survey); anova(fit)
# report the means and the number of subjects/cell
print(model.tables(fit,"means"), digits=3)

response 代表的是y值

residuals是不等于均值的值（事实上它就是残差）

# ANOVA fit
fit[[1]] # coefficients
fit[[2]] # residuals 残差
fit[[5]] # predicted (fitted.values)
all <- cbind (rodent.survey,
predicted=fit[[5]], residual=fit[[2]])
all

par(mfrow=c(2, 2))
plot(fit)

如何判断值是否是正态的，后面会再说。

shapiro检验可以用来检验正态性。

t检验 随机性 连续性 正态性

方差分析需要 方差齐性，独立性，随机性。

检查齐性——Bartlett‘s test

例子

Bartlett和Fisher检验是检查方差的齐性。

前者要求数据的正态性，而后者对数据的正态性的要求并不高。

参考：http://www.datasoldier.net/archives/350

# Bartlett Test of Homogeneity of Variances (parametric)
bartlett.test (weight~site, data=rodent.survey)
# Figner-Killeen Test of Homogeneity of Variances (non-parametric)
fligner.test (weight~site, data=rodent.survey) 

方差分析，

harmonic 1/（1/n+1/m）

mse 残差的平方，处以自由度40-1-6

lwr 和upr是啥意思，是lower和upper的意思，p adj是指调整过的p值。

# Tukey's multiple comparisons tests (Honestly Significantly Different (HSD) test)
TukeyHSD(fit)$site

使用一般线性模型进行数据分析的通用方法

df自由度 ms均方差

# R code
fit <- aov (Width ~ Host, data=tick)
TukeyHSD(fit)$Host

组平均偏离总平均的host

!

datafilename <- "http://personality-project.org/R/datasets/R.appendix1.data" #tell where the data come from
data.ex1 <- read.table(datafilename, header = T) #read the data into a table
aov.ex1 <- aov(Alertness ~ Dosage, data = data.ex1) #do the analysis of variance
print(model.tables(aov.ex1,"means"), digits = 3) #report the means and the number of subjects/cell
boxplot(Alertness ~ Dosage, data = data.ex1) #graphical summary
summary(aov.ex1) #show the results

2.Randomized block designs

用block来解释了一部分。负面信息的话，减少了自由度。这个有没有必要，只要取决于

#Randomized Block Design
#Carbon dioxygen density at 8 incubators and 4 treatments
CO2 <- data.frame(ID=1:32, group=NA, treat=NA, density=NA)
n <- 0
for(i in 1:8){
for(j in c('A','B','C','D')){
n <- n+1
CO2$group[n] = i
CO2$treat[n] = j
}}
CO2$group <- factor(CO2$group)
CO2$treat <- factor(CO2$treat)
CO2$density <- c(5.27,5.27,5.94,5.53,5.27,5.22,4.88,4.96,5.88,5.83,
5.38,5.53,5.44, 5.38,5.27,5.32,5.66, 5.44,5.38,4.88,6.22,
6.22,5.61,5.92,5.83,5.72,5.38,4.88,5.27,5.11,5.12,4.44)
fit1 <- aov(density ~ treat, data = CO2) # one way ANOVA
fit2 <- aov(density ~ group + treat, data = CO2) # Randomized Block Design
summary(fit1)
summary(fit2)
par(mfrow=c(2,2))
plot(fit2)

head(CO2)
library(reshape2)
Tab = melt(CO2, id=c("group", "treat"), na.rm=TRUE)
Mat = acast(Tab, group ~ treat ~ variable)
dim(Mat); Mat[ , , 2]

Factorial

3.Two way ANOVA

interaction：a, b, c a change, b and c’s guanxi change.

有没有相互作用，可以用R来画

# Two-way Interaction Plot
attach(mtcars)
gears <- factor(gear)
cyl <- factor(cyl)
interaction.plot(cyl, gear, mpg,
type="b", col=c(1:3),
leg.bty="o", leg.bg="beige",
lwd=2, pch=c(18,24,22),
xlab="Number of Cylinders",
ylab="Mean Miles Per Gallon",
main="Interaction Plot")

方差可以可减性

显著的情况，需要再做分析，

没有显著的情况下，不考虑相互作用

leveneTest是检方差齐不齐的情况

# two explanatory variables
head(mtcars)
str(mtcars)
mtcars$cyl = as.factor(mtcars$cyl)
mtcars$gear = as.factor(mtcars$gear)
library(car)
leveneTest(mpg ~ cyl * gear, data = mtcars)

双因素分析代码

# Two Way ANOVA
weight.gain <- data.frame(ID=1:60, amount=NA, food=NA, gain=NA)
n <- 0
for(i in c('high','low')){
for(j in c('beef','cereal','port')){
for(k in 1:10){
n <- n+1
weight.gain$amount[n] = i
weight.gain$food[n] = j
}
}
}
weight.gain$gain <- c(73,102,118,104,81,107,100,87,117,111,
98,74,56,111,95,88,82,77,86,92, 94,79,96,98,102,102,108,91,120,105,
90,76,90,64,86,51,72,90,95,78, 107,95,97,80,98,74,74,67,89,58,
49,82,73,86,81,97,106,70,61,82)
fit <- aov(gain ~ amount + food + amount:food, data = weight.gain)
fit <- aov(gain ~ amount * food, data = weight.gain) # same thing 这个是包括的线性和交互项。
summary(fit)
par(mfrow=c(2,2)); plot(fit)

# One Way Anova (Completely Randomized Design)
attach(mtcars); head(mtcars)
fit1 <- aov(mpg ~ cyl, data=mtcars)
B=cyl; A=gear; x=wt
# Randomized Block Design (B is the blocking factor)
fit2 <- aov(mpg ~ A + B, data=mtcars)
# Two Way Factorial Design
fit3 <- aov(mpg ~ A + B + A:B, data=mtcars)
fit4 <- aov(mpg ~ A*B, data=mtcars) # same thing
# Analysis of Covariance
fit5 <- aov(mpg ~ A + x, data=mtcars)
summary(fit1) # display Type I ANOVA table
drop1(fit1,~.,test="F") # type III SS and F Tests

数据的取值和顺序没有关系。

（1）方差分析的假定

（2）检查数据的独立性

如何验证数据的独立性（不能衡量齐独立性）

par(mfrow=c(1,2))
D = rnorm(30)
plot(D[1:29],D[2:30])
D.order = sort(D)
plot(D.order[1:29],D.order[2:30])

par(mfrow=c(2,2))
plot(fit)

3.Two way ANOVA without replication

这里是没有重复的双因素，交互作用就是差别。交互作用

fit1 <- aov(W ~ plot + type, data = mydata)

不光在上面的情况，同时有一个因素只有两个值。

下面加上了配对，上面没有表现配对的情况。多加一个变量

多了个2，是更保守的情况。nh是调和平均数

第一个出现是对照。本质是比较组间差异与组内差异的比较，当组内差大得多，我们称让显著。
两个是t分析，两列多则是f分析，但是又不完全是这个，f分析配对。
随机区组和交互作用。
没法量化交互作用的话，样本量不够，我们可以配对，线上的偏移是说明的，就是他们的相互作用。