r语言 单细胞分析 基因表达差异 r语言单因素

方差分析

基本原理、一个因变量的单因子独立样本、双因子独立样本

• 单因子方差分析

library(reshape2)
table8_2<-melt(table,variable.name="品种",value.name = "产量")
mode1<-aov(table8_2$产量~table8_2$品种,data=table)
summary(mode1)

• 方差分析模型的参数估计

mode1$coefficients

• 均值图

library(gplots)
plot(table8_2$产量~table8_2$品种,data=table)

• 效应量分析

library(DescTools)
mode1<-aov(table8_2$产量~table8_2$品种)
EtaSq(mode1,anova = T)

• LSD

library(agricolae)
LSD<-LSD.test(mode1,"table8_2$品种")
#更多
library(DescTools)
PostHocTest(mode1,method = "lsd")

• HSD

TukeyHSD(mode1)
#更多
library(agricolae)
HSD<-HSD.test(mode1,"table8_2$品种")
#配对差值置信区间的比较图
plot(TukeyHSD(mode1))

• 双因子方差分析
• 主效应方差分析模型

library(reshape2)
table<-melt(table,variable.name = "品种",value.name = "产量")
model<-aov(table$产量~table$施肥方式+table$品种,data=table)
summary(model)

• 主效应参数估计

model$coefficients

• 主效应效应量分析

library(DescTools)
EtaSq(mode,anova=T)

• 交互效应方差分析模型

mode<-aov(table$产量~table$施肥方式+table$品种+table$品种:table$施肥方式,data=table)
summary(mode)

• 交互效应参数估计

mode$coefficients

• 主效应和交互效应图

library(HH)
interaction2wt(table$产量~table$施肥方式+table$品种,data=table)

• 交互效应效应量分析

mode<-aov(table$产量~table$施肥方式+table$品种+table$品种:table$施肥方式,data=table)
EtaSq(mode,anova=T)

• 正态性
• 正态性检验QQ 图

qqnorm(table$品1,xlab="qiwang",ylab="guance",datax=TRUE,main="1QQ")

• 正态性检验sw检验

shapiro.test(table$产量)

• 正态性检验ks检验

ks.test(table$产量,"pnorm",mean(table$产量),sd(table$产量))

• 方差齐性
• Bartlett方差齐性检验

bartlett.test(table$产量~table$品种,data=table)

• Levene方差齐性检验

library(car)
leveneTest(table$产量~table$品种,data=table)

---

文章目录

• 方差分析
• 8.1 方差分析的原理

• 8.1.2 什么是方差分析
• 8.1.2 误差分解
• 8.1.3 数学模型

• 8.2 单因子方差分析

• 8.2.1 效应检验

• 代码

• 8.2.2 效应量分析

• 代码

• 8.2.3 多重比较

• Fisher 的 LSD方法

• 代码

• Tukey-Kramer 的 HSD方法

• 代码

• 8.3 双因子方差分析

• 8.3.1 数学模型
• 8.3.2 主效应分析

• 1、效应检验

• 代码

• 2、效应量分析

• 代码

• 8.3.3 交互效应分析

• 1、效应检验

• 代码

• 2、效应量分析

• 8.4 方差分析的假定及其检验

• 8.4.1 正态性检验

• 1、QQ图

• 代码

• 2、SW和KS检验

• 代码

• 8.4.2 方差齐性检验

• 1、图示法

• 箱线图
• 残差图

• 2、检验法

• Bartlett

• 代码

• Levene

• 代码

• 3、独立性

---

8.1 方差分析的原理

8.1.2 什么是方差分析

分析类别自变量对数值因变量影响的一种统计方法

自变量对因变量的影响：自变量效应

通过对数据误差分析检验这种效应是否显著

因子：小麦品种（类别变量） 处理、
水平：品种1、品种2、品种3是因子的3个不同取值
实验单元：地块（接受处理的对象或实体）

8.1.2 误差分解

• SST=SSA+SSE

8.1.3 数学模型

• 处理误差不显著：每种处理的总体均值相等
• 处理误差显著：每种处理的总体均值至少有一处不相等

8.2 单因子方差分析

8.2.1 效应检验

代码

• 方差分析模型

> table<-read.csv("/Users/zhourui/Documents/example8_1.csv")
> library(reshape2)
> table8_2<-melt(table,variable.name="品种",value.name = "产量")
No id variables; using all as measure variables
#转化为长格式
> mode1<-aov(table8_2$产量~table8_2$品种,data=table)
#拟合方差分析模型
> summary(mode1)
              Df Sum Sq Mean Sq F value   Pr(>F)    
table8_2$品种  2    560  280.00   12.31 0.000158 ***
Residuals     27    614   22.74                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sum Sq：品种效应和随机效应的平方和

Df：自由度

Mean Sq：均方差

F value：检验统计量值

P：P值

将P与alpha比较：P<alpha，拒绝H0，至少有一个不等于零，品种对产量的影响效应显著

• 方差分析模型的参数估计

> mode1$coefficients
       (Intercept) table8_2$品种品种2 table8_2$品种品种3 
                84                -10                 -2

Intercept：截距

• 绘制均值图

> library(gplots)
> plot(table8_2$产量~table8_2$品种,data=table)

x轴：品种；y轴：产量

8.2.2 效应量分析

代码

> library(DescTools)
> mode1<-aov(table8_2$产量~table8_2$品种)
> EtaSq(mode1,anova = T)
                 eta.sq eta.sq.part  SS df        MS       F            p
table8_2$品种 0.4770017   0.4770017 560  2 280.00000 12.3127 0.0001583995
Residuals     0.5229983          NA 614 27  22.74074      NA           NA

8.2.3 多重比较

分析差异到底出现在哪些品种之间
谁和谁之间均值不等

Fisher 的 LSD方法

代码

• 基本信息

> library(agricolae)
> LSD<-LSD.test(mode1,"table8_2$品种")
> LSD
$statistics
   MSerror Df Mean       CV  t.value      LSD
  22.74074 27   80 5.960907 2.051831  4.375813

$parameters
        test p.ajusted        name.t ntr alpha
  Fisher-LSD      none table8_2$品种   3  0.05

$means
      table8_2$产量      std  r      LCL      UCL Min Max   Q25  Q50   Q75
品种1            84 4.546061 10 80.90583 87.09417  78  92 81.00 83.5 86.75
品种2            74 4.447221 10 70.90583 77.09417  66  81 72.00 72.5 77.00
品种3            82 5.270463 10 78.90583 85.09417  76  89 77.25 81.5 87.00

$comparison
NULL

$groups
      table8_2$产量 groups
品种1            84      a
品种3            82      a
品种2            74      b

attr(,"class")
[1] "group"

• 更多信息

> library(DescTools)
> PostHocTest(mode1,method = "lsd")

  Posthoc multiple comparisons of means : Fisher LSD 
    95% family-wise confidence level

$`table8_2$品种`
            diff     lwr.ci    upr.ci    pval    
品种2-品种1  -10 -14.375813 -5.624187   7e-05 ***
品种3-品种1   -2  -6.375813  2.375813 0.35666    
品种3-品种2    8   3.624187 12.375813 0.00085 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

有*代表有差异

Tukey-Kramer 的 HSD方法

代码

• 多重比较的HSD方法

> TukeyHSD(mode1)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = table8_2$产量 ~ table8_2$品种)

$`table8_2$品种`
            diff        lwr       upr     p adj
品种2-品种1  -10 -15.287702 -4.712298 0.0002017
品种3-品种1   -2  -7.287702  3.287702 0.6215828
品种3-品种2    8   2.712298 13.287702 0.0023770

• 绘制配对差值置信区间的比较图

> plot(TukeyHSD(mode1))

y轴：品种

• 更多信息

> library(agricolae)
> HSD<-HSD.test(mode1,"table8_2$品种")
> HSD
$statistics
   MSerror Df Mean       CV      MSD
  22.74074 27   80 5.960907 5.287702

$parameters
   test        name.t ntr StudentizedRange alpha
  Tukey table8_2$品种   3         3.506426  0.05

$means
      table8_2$产量      std  r Min Max   Q25  Q50   Q75
品种1            84 4.546061 10  78  92 81.00 83.5 86.75
品种2            74 4.447221 10  66  81 72.00 72.5 77.00
品种3            82 5.270463 10  76  89 77.25 81.5 87.00

$comparison
NULL

$groups
      table8_2$产量 groups
品种1            84      a
品种3            82      a
品种2            74      b

attr(,"class")
[1] "group"

8.3 双因子方差分析

两个类别自变量对数值因变量影响的方差分析

• 主效应、无重复双因子分析：只考虑两个因子对因变量的单独影响
• 交互效应、可重复双因子分析：还考虑两个因子的搭配

8.3.1 数学模型

8.3.2 主效应分析

1、效应检验

• 检验因子A的假设：
• 检验因子B的假设：
  
  
  
  

代码

• 方差分析模型

> table<-read.csv("/Users/zhourui/Documents/table8_4.csv")
> library(reshape2)
> table<-melt(table,variable.name = "品种",value.name = "产量")
> model<-aov(table$产量~table$施肥方式+table$品种,data=table)
> summary(model)
               Df Sum Sq Mean Sq F value   Pr(>F)    
table$施肥方式  1    480   480.0   93.13 4.42e-10 ***
table$品种      2    560   280.0   54.33 5.18e-10 ***
Residuals      26    134     5.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

P值接近于零（Pr(>F)）：两个因子对产量均有显著影响

• 主效应方差分析的参数估计

> model$coefficients
     (Intercept) table$施肥方式乙  table$品种品种2  table$品种品种3 
              80                8              -10               -2

Intercept（截距）：不考虑品种和施肥方式的影响时产量的均值
施肥方式乙：附加效应（参照物为施肥方式甲）
品种2、品种3：附加效应（参照物为品种1）

2、效应量分析

在两因子的主效应方差分析中，衡量效应量大小的统计量是偏效应

主效应量有可加性，两主效应量相加为两因子的联合效应，反应两因子联合起来对因变量的影响大小

代码

#主效应方差分析的效应量
> library(DescTools)
> EtaSq(mode,anova=T)
                  eta.sq eta.sq.part  SS df         MS        F
table$施肥方式 0.4088586   0.7817590 480  1 480.000000 93.13433
table$品种     0.4770017   0.8069164 560  2 280.000000 54.32836
Residuals      0.1141397          NA 134 26   5.153846       NA
                          p
table$施肥方式 4.422633e-10
table$品种     5.184290e-10
Residuals                NA

eta.sq：每个因子主效应量
eta.sq.part：每个因子的偏效应量

8.3.3 交互效应分析

两因子搭配对因变量的交互作用

1、效应检验

代码

• 交互效应方差分析

> mode<-aov(table$产量~table$施肥方式+table$品种+table$品种:table$施肥方式,data=table)
> summary(mode)
                          Df Sum Sq Mean Sq F value   Pr(>F)    
table$施肥方式             1  480.0   480.0   93.20 9.73e-10 ***
table$品种                 2  560.0   280.0   54.37 1.22e-09 ***
table$施肥方式:table$品种  2   10.4     5.2    1.01    0.379    
Residuals                 24  123.6     5.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

P值（Pr(>F)小于0.05：两个因子对产量的影响显著
P值（Pr(>F)大于0.05：交互效应对产量的影响并不显著

• 交互效应方差分析模型的参数估计

> mode$coefficients
                     (Intercept)                 table$施肥方式乙 
                            80.2                              7.6 
                 table$品种品种2                  table$品种品种3 
                            -9.6                             -3.0 
table$施肥方式乙:table$品种品种2 table$施肥方式乙:table$品种品种3 
                            -0.8                              2.0

• 主效应和交互效应图

> library(HH)
> interaction2wt(table$产量~table$施肥方式+table$品种,data=table)

2、效应量分析

• 交互效应方差分析的效应量

> mode<-aov(table$产量~table$施肥方式+table$品种+table$品种:table$施肥方式,data=table)
> EtaSq(mode,anova=T)
                               eta.sq eta.sq.part    SS df     MS
table$施肥方式            0.408858603  0.79522863 480.0  1 480.00
table$品种                0.477001704  0.81919251 560.0  2 280.00
table$施肥方式:table$品种 0.008858603  0.07761194  10.4  2   5.20
Residuals                 0.105281090          NA 123.6 24   5.15
                                  F            p
table$施肥方式            93.203883 9.729467e-10
table$品种                54.368932 1.220666e-09
table$施肥方式:table$品种  1.009709 3.792836e-01
Residuals                        NA           NA

8.4 方差分析的假定及其检验

8.4.1 正态性检验

要求每个处理所对应的总体都服从正态分布；

对任意一个处理，其观察值是来自正态分布总体的简单随机样本

1、QQ图

数据量较大：Q-Q图

数据量较小：数据合并，绘制正态概念图

代码

• 直接绘制Q-Q图

> par(mfrow=c(1,3),cex=0.6,mai=c(0.5,0.5,0.2,0.1))
> qqnorm(table$品种1,xlab = "qiwang",ylab = "guance",datax=TRUE,main="1QQ")
> qqline(table$品种1,datax=TRUE)
> qqnorm(table$品种2,xlab = "qiwang",ylab = "guance",datax=TRUE,main="1QQ")
> qqline(table$品种2,datax=TRUE)
> qqnorm(table$品种3,xlab = "qiwang",ylab = "guance",datax=TRUE,main="1QQ")
> qqline(table$品种3,datax=TRUE)

• 数据合并后Q-Q图

> table<-read.csv("/Users/zhourui/Documents/example8_2.csv")
> par(cex=0.8,mai=c(0.7,0.7,0.1,0.1))
> qqnorm(table$产量,xlab = "qiwang",ylab = "guance",datax=TRUE,main=" ")
> qqline(table$产量,datax=TRUE,col="red",lwd=2)
> par(fig=c(0.08,0.5,0.5,0.98),new=TRUE)
> hist(table$产量,xlab = "chanliang",ylab = " ",freq = FALSE,col="lightblue",cex.axis=0.7,cex.lab=0.7,main="")
> lines(density(table$产量),col="red",lwd=2)
> box()

2、SW和KS检验

代码

• sw

#8.2 
> shapiro.test(table$产量)

	Shapiro-Wilk normality test

data:  table$产量
W = 0.97299, p-value = 0.6237

• KW

#8.2
> ks.test(table$产量,"pnorm",mean(table$产量),sd(table$产量))

	One-sample Kolmogorov-Smirnov test

data:  table$产量
D = 0.097706, p-value = 0.9369
alternative hypothesis: two-sided

P值（p-value）均大于0.05，不拒绝H0，产量服从正态分布

SW对偏理性较为敏感，轻微偏理往往也会拒绝原假设

方差分析对正态性要求较为宽松，正态性略不满足时，分析结果影响不大

8.4.2 方差齐性检验

要求各处理的总体方差相等

1、图示法

箱线图

观察离散程度

离散程度大致相同：方差齐性假定满足

残差图

残差：实际观察值与预测值的差值

标准化残差：残差除以残差的标准差

• 方差分析的残差图和标准化残差的Q-Q图
• 8.2

> mode<-aov(table$产量~table$品种,data=table)
> par(mfrow=c(1,2),mai=c(0.5,0.5,0.2,0.1),cex=0.6,cex.main=0.7)
> plot(mode,which=c(1,2))

plot(model)默认绘制4个图，which=1和which=2代表仅输出所需的两个图

第一张图

• 方差齐性假定判定：

第一个图是残差图，图中横坐标是模型的拟合值，纵坐标是预测残差
可以看出无离群点，拟合值和残差的散点随机分布在一个水平带之内，离散程度基本一样
满足方差齐性假定

• 评价方差分析模型的拟合效果
  如果拟合的好，那么拟合值与残差值的观测点就应该在一条水平线附近波动
  拟合效果好

第二张图
各点在直线周围分布，没有固定模式，可认为正态性假定基本成立

• 8.2

> table<-read.csv("/Users/zhourui/Documents/example8_5.csv")
> mode<-aov(table$产量~table$品种+table$施肥方式+table$品种:table$施肥方式,data=table)
> par(mfrow=c(1,2),mai=c(0.5,0.5,0.2,0.1),cex=0.6,cex.main=0.7)
> plot(mode,which=1:2)

2、检验法

Bartlett

代码

• 8.5

> bartlett.test(table$产量~table$品种,data=table)

	Bartlett test of homogeneity of variances

data:  table$产量 by table$品种
Bartlett's K-squared = 0.30152, df = 2, p-value = 0.8601

> bartlett.test(table$产量~table$施肥方式,data=table)

	Bartlett test of homogeneity of variances

data:  table$产量 by table$施肥方式
Bartlett's K-squared = 0.42431, df = 1, p-value = 0.5148

P值均大于0.05，不拒绝H0，不同品种的产量和不同施肥方式的产量均满足方差齐性

Levene

代码

• 8.5

> library(car)
> leveneTest(table$产量~table$品种,data=table)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  0.9292 0.4071
      27

> leveneTest(table$产量~table$施肥方式,data=table)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.2323 0.6336
      28

P值均大于0.05，不拒绝H0，不同品种的产量和不同施肥方式的产量均满足方差齐性

在方差分析中，对方差齐性的要求较为宽松，方差略有不齐时，对分析结果要求不大

3、独立性

要求每个样本数据来自不同处理的独立样本

在实验设计前确定，不需要检验