R语言学习系列27-方差分析

一、方差分析原理

1. 方差分析概述

方差分析可用来研究多个分组的均值有无差异，其中分组是按影响因素的不同水平值组合进行划分的。

方差分析是对总变异进行分析。看总变异是由哪些部分组成的，这些部分间的关系如何。

方差分析，是用来检验两个或两个以上均值间差别显著性（影响观察结果的因素：原因变量（列变量）的个数大于2，或分组变量（行变量）的个数大于1）。一元时常用F检验（也称一元方差分析），多元时用多元方差分析（最常用Wilks’∧检验）。

方差分析可用于：

（1）完全随机设计（单因素）、随机区组设计（双因素）、析因设计、拉丁方设计和正交设计等资料；

（2）可对两因素间交互作用差异进行显著性检验； （3）进行方差齐性检验。

要比较几组均值时，理论上抽得的几个样本，都假定来自正态总体，且有一个相同的方差，仅仅均值可以不相同。还需假定每一个观察值都由若干部分累加而成，也即总的效果可分成若干部分，而每一部分都有一个特定的含义，称之谓效应的可加性。所谓的方差是离均差平方和除以自由度，在方差分析中常简称为均方（Mean Square）。

2. 基本思想

基本思想是，将所有测量值上的总变异按照其变异的来源分解为多个部份，然后进行比较，评价由某种因素所引起的变异是否具有统计学意义。

根据效应的可加性，将总的离均差平方和分解成若干部分，每一部分都与某一种效应相对应，总自由度也被分成相应的各个部分，各部分的离均差平方除以各自的自由度得出各部分的均方，然后列出方差分析表算出F检验值，作出统计推断。

方差分析的关键是总离均差平方和的分解，分解越细致，各部分的含义就越明确，对各种效应的作用就越了解，统计推断就越准确。

效应项与试验设计或统计分析的目的有关，一般有：主效应（包括各种因素），交互影响项（因素间的多级交互影响），协变量（来自回归的变异项），等等。

当分析和确定了各个效应项S后，根据原始观察资料可计算出各个离均差平方和SS，再根据相应的自由度df，由公式MS=SS/df，求出均方MS，最后由相应的均方，求出各个变异项的F值，F值实际上是两个均方之比值，通常情况下，分母的均方是误差项的均方。

根据F值的分子、分母均方的自由度f1和f2，在确定显著性水平为α情况下，由F(f1, f2)临界值表查得单侧Fα界限值。当Fα，不拒绝原假设H0，说明不拒绝这个效应项的效应为0的原假设，也即这个效应项是可能对总变异没有实质影响的；若F>Fα则P值≤α，拒绝原假设H0，也即这个效应项是很可能对总变异有实质影响的。

3.方差分析的实验设计

为了确定方差分析表中各个有关效应项，需要在试验设计阶段就作出安排，再根据设计要求进行试验，得出原始观察值，按原来设计方案算出方差分析表中的各项。

在试验设计阶段通常需要考虑如下4个方面： （1）研究的因变量

即试验所要观察的主要指标，一次试验时可以有多个观察指标，方差分析时也可以同时对多个因变量进行分析；

（2）因素和水平

试验的因素（factor）可以是品种、人员、方法、时间、地区等等，因素所处的状态叫水平（level）。在每一个因素下面可以分成若干水平。

（3）因素间的交互影响

多因素的试验设计，有时需要分析因素间的交互影响（interaction），2个因素间的交互影响称为一级交互影响（A×B）；3个因素间的交互影响称为二级交互影响（A×B×C）。

当交互影响项呈现统计不显著时，表明各个因素，当呈现统计显著时，就需要列出这个交互影响项的效应，以助于作出正确的统

计推断。

举例解释上述概念：要考察焦虑症的治疗疗效，一个因素是治疗方案，有2种治疗方案，即该因素有2个水平；（治疗方案称为组间因子，因为每个患者只能被分配到一个组别中，没有患者同时接受两种治疗）；再考虑一个因素治疗时间，也有两个水平：治疗5周和治疗6个月，同一患者在5周和6个月不止一次地被测量（两次），称为重复测量（治疗时间称为组内因子，因为每个患者在所有水平下都进行了测量）。

建立方差分析模型时，既要考虑两个因素治疗方案和治疗时间（主效应），又要考虑治疗方案和时间的交互影响（交互效应），此时即两因素混合模型方差分析。

当某个因素的各个水平下的因变量的均值呈现统计显著性差异时，必要时可作两两水平间的比较，称为均值间的两两比较。 二、R语言实现

方差分析对数据的要求：满足正态性（来自同一正态总体）和方差齐性（各组方差相等），在这两个条件下，若各组有差异，则只可能是来自影响因素的不同水平。

用aov()函数进行方差分析，基本格式为：

aov(formula, data=NULL, projections=FALSE, qr=TRUE,

contrasts=NULL, ...)

其中，formula为方差分析公式；

data为数据框；

projection设置是否返回预测结果；

qr设置是否返回QR分解结果； contrasts为公式中一些因子的列表。

formula公式的表示：（y为因变量，ABC为分组因子） 符号 ~ eg：y~A+B+C 用法 分隔符号，左边为响应变量，右边为解释变量 + 分隔解释变量 ： * ^ .

常见研究设计的表达式：（小写字母表示定量变量，大写字母表 示组别因子，Subject是对被试者独有的标识变量）

设计 单因素ANOVA 含单个协变量的单因素ANCOVA 双因素ANOVA 含两个协变量的双因素ANCOVA 随机化区组 单因素组内ANOVA 含单个组内因子(W)和单个组间因子(B)的重复测量ANOVA 表达式 y~A y~x+A y~A*B y~x1+x2+A*B y~B+A, B为区组因子 y~A+Error(Subject/A) y~B*W+Error(Subject/W) 表示变量的交互项 eg：y~A+B+A:B 表示所有可能交互项 eg：y~A*B*C可展开为：y~A+B+C+A:B+A:C+B:C+A:B:C 表示交互项达到次数 eg：y~(A+B+C)^2展开为：y~A+B+C+A:B+A:C+B:C 表示包含除因变量外的所有变量 eg：若一个数据框包括变量y,A、B和C，代码y~.可展开为y~A+B+C 注意：非均衡设计时或存在协变量时，效应项的顺序对结果影响较大，越基础的效应越需要放在表达式前面，首先是协变量、然后是

主效应、接着是双因素的交互项，再接着是三因素的交互项。若研究不是正交的，一定要谨慎设置效应的顺序。

有三种类型的方法可以分解y~A+B+A:B右边各效应对y所解释的方差：

类型I（序贯型）

效应根据表达式中先出现的效应做调整。A不做调整，B根据A调整，A:B交互项根据A和B调整。

类型II（分层型）

效应根据同水平或低水平的效应做调整。A根据B调整，B依据A调整，A:B交互项同时根据A和B调整。

类型III（边界型）

每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:B交互项根据A和B调整。

R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型II或类型III方法的选项，而aov()函数使用的是类型I方法。若想使结果与其他软件（如SAS和SPSS）提供的结果保持一致，可以使用Anova()函数。

三、单因素方差分析

1个因变量，1个影响因素：

总差异Yij = 平均差异μ + 因素差异αi + 随机差异εij

例1比较4种品牌的胶合板的耐磨性，各抽取5个样品，相同转速磨损相同时间测得磨损深度（mm），比较4个品牌胶合板的耐磨性有无差异？部分数据如下（ex27_ex1.Rdata）：

setwd(\"E:/办公资料/R语言/R语言学习系列/codes\") load(\"ex27_ex1.Rdata\") head(datas) wear brand 1 2.30 A 2 2.32 A 3 2.40 A 4 2.45 A 5 2.58 A 6 2.35 B attach(datas)

table(brand)#各组的样本数 brand A B C D 5 5 5 5

aggregate(wear,by=list(brand),mean)#各组均值 Group.1 x 1 A 2.410 2 B 2.404 3 C 2.046 4 D 2.572

aggregate(wear,by=list(brand),sd)#各组标准差 Group.1 x 1 A 0.11269428 2 B 0.11760102 3 C 0.11216060

4 D 0.03271085 library(car)

qqPlot(lm(wear~brand,data=datas),simulate=TRUE)#用Q-Q图检验数据的正态性

leveneTest(wear~as.factor(brand),data=datas)#方差齐性检验

Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F)

group 3 0.6987 0.56 16

fit<-aov(wear~brand,data=datas)#单因素方差分析（检验组间差异） summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

brand 3 0.7398 0.24660 24.55 3.15e-06 *** Residuals 16 0.1607 0.01005 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

说明：方差齐性检验，原假设H0：方差齐，p值=0.56>0.05, 故接受原假设，即方差齐。

单因素方差分析结果，brand是因素，Residuals是残差，各列依次为自由度、平方和、均方和、F统计量，p值=3.15e-06<0.05, 拒绝原假设，即不同品牌的磨损（均值）有显著差别。

library(gplots)

plotmeans(wear~brand,xlab=\"品牌\磨损\")#图形展示带95%置信区间的各组均值

通过前面的分析知道，不同品牌的磨损（均值）有显著差别，但并不知道哪个品牌与其它品牌有显著差别。TukeyHSD()函数提供了对各组均值差异的成对检验。

TukeyHSD(fit)

Tukey multiple comparisons of means 95% family-wise confidence level

Fit: aov(formula = wear ~ brand, data = datas) $brand

difflwrupr p adj

B-A -0.006 -0.18735345 0.1753535 0.9996826 C-A -0.3 -0.54535345 -0.18265 0.0001610 D-A 0.162 -0.01935345 0.3433535 0.0886142 C-B -0.358 -0.53935345 -0.17665 0.0001929 D-B 0.168 -0.01335345 0.3493535 0.0744337 D-C 0.526 0.344655 0.7073535 0.0000019

说明：可以看出（H0：无差异），B与A的差异非常不显著，C与A、C与B、D与C的差异非常显著。

multcomp包中的glht()函数提供了更为全面的多重均值比较方法。

library(multcomp) attach(datas)

tuk<- glht(fit, linfct = mcp(brand = \"Tukey\")) #注意datas$brand必须是因子型 summary(tuk)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: aov(formula = wear ~ brand, data = datas)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|) B - A == 0 -0.00600 0.06339 -0.095 0.9997 C - A == 0 -0.300 0.06339 -5.742 <0.001 *** D - A == 0 0.16200 0.06339 2.556 0.0886 . C - B == 0 -0.35800 0.06339 -5.8 <0.001 *** D - B == 0 0.16800 0.06339 2.650 0.0743 . D - C == 0 0.52600 0.06339 8.298 <0.001 *** ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method)

plot(cld(tuk, level = 0.05), col = \"lightgrey\")

说明：标记相同字母（标记b）的品牌ABD认为是无显著差异，在同一亚组，而品牌C（标记a）与另外三个品牌有显著差异。

另外，也可以进行多重t检验，使用函数：

pairwise.t.test(x, g, p.adjust.method=,...)

其中，x为因变量，g为因子型的分组变量；

p.adjust.method设置p值的修正方法，由于多次重复t检验会大大增加犯第一类错误的概率，为此要进行p值的修正，使用bonferroni法修正效果较好。

pairwise.t.test(wear,brand,p.adjust.method=\"bonferroni\")

Pairwise comparisons using t tests with pooled SD data: wear and brand

A B C B 1.00000 - - C 0.00018 0.00022 - D 0.12695 0.10474 2.1e-06

P value adjustment method: bonferroni

说明：原假设H0:无差异，可见A与B无差异，C与ABD有显著差异。

最后，方差分析对离群点非常敏感，检验是否有离群点：

library(car)

outlierTest(fit)

No Studentized residuals with Bonferonni p < 0.05 Largest |rstudent|:

rstudent unadjusted p-value Bonferonni p 9 2.528103 0.023182 0.463

说明：经检验无离群点。

三、两因素方差分析

1个因变量，2个影响因素：

总差异Yijk = 平均差异μ + 因素1差异αi + 因素2差异βi

+因素1,2交互作用差异γij+随机差异εijk

例2研究60只豚鼠的牙齿生长数据，按2种喂食方法：橙汁、维生素C，各喂食方法中抗坏血酸含量都有3个水平：0.5mg、1mg、2mg，

分配为6组，每组各10只，牙齿长度为因变量。做两因素方差分析。

attach(ToothGrowth) head(ToothGrowth) lensupp dose 1 4.2 VC 0.5 2 11.5 VC 0.5 3 7.3 VC 0.5 4 5.8 VC 0.5 5 6.4 VC 0.5 6 10.0 VC 0.5

table(supp, dose)#各组样本数相同，即为均衡设计 dose

supp 0.5 1 2 OJ 10 10 10 VC 10 10 10

aggregate(len, by=list(supp, dose), mean)#计算各组均值 Group.1 Group.2 x 1 OJ 0.5 13.23 2 VC 0.5 7.98 3 OJ 1.0 22.70 4 VC 1.0 16.77 5 OJ 2.0 26.06 6 VC 2.0 26.14

aggregate(len, by=list(supp, dose), sd)#计算各组标准差 Group.1 Group.2 x 1 OJ 0.5 4.459709 2 VC 0.5 2.746634 3 OJ 1.0 3.910953 4 VC 1.0 2.515309 5 OJ 2.0 2.655058 6 VC 2.0 4.797731

bartlett.test(len~supp,data=ToothGrowth)#关于因素supp的方差齐性检验

Bartlett test of homogeneity of variances data: len by supp

Bartlett's K-squared = 1.4217, df = 1, p-value = 0.2331 bartlett.test(len~dose,data=ToothGrowth)#关于因素dose的方差齐性检验

Bartlett test of homogeneity of variances data: len by dose

Bartlett's K-squared = 0.66547, df = 2, p-value = 0.717 fit<-aov(len~supp*dose,data=ToothGrowth) #做两因素方差分析，考虑全部效应

summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

supp 1 205.4 205.4 12.317 0.0004 *** dose 1 2224.3 2224.3 133.415 < 2e-16 *** supp:dose 1 88.9 88.9 5.333 0.024631 * Residuals 56 933.6 16.7 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

说明：可以看出，主效应supp和dose都非常显著（p值都远小于0.05），交互效应也显著（p值=0.0246<0.05）。若交互作用不显著，可以可以只做去掉交互效应的方差分析。

图形化展示两因素方差分析的交互效应：

par(mfrow=c(1,2))

interaction.plot(dose, supp, len, type=\"b\d\se and Supp\")

interaction.plot(supp, dose, len, type=\"b\ \"blue\"), pch = c(16, 18), main=\"Interaction between Supp and Dose\")

说明：有一个图的线有交叉，说明有交互作用。可以看出随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变长；。对于0.5 mg 和1 mg剂量，橙汁比维生素C更能促进牙齿生长；对于2 mg剂量

的抗坏血酸，两种喂食方法下牙齿长度增长相同。

也可以用HH包中的interaction2wt()函数（也适合三因素方差分析）来展示更全面的可视化结果：

library(HH)

interaction2wt(len~supp*dose)

三、重复测量方差分析

重复测量方差分析，即受试者被测量不止一次。

例3（1个组内1个组间因子的重复测量）在某浓度CO2的环境中，对寒带植物（来自魁北克）和非寒带植物的（来自密西西比）光合作用率进行比较。因变量uptake为CO2吸收量，自变量Type（组间因子）为植物类型，自变量conc（组内因子）为七种水平的CO2浓度。

attach(CO2)

head(CO2)#注意CO2是长格式的数据

Plant Type Treatmentconc uptake 1 Qn1 Quebec nonchilled 95 16.0 2 Qn1 Quebec nonchilled 175 30.4 3 Qn1 Quebec nonchilled 250 34.8 4 Qn1 Quebec nonchilled 350 37.2 5 Qn1 Quebec nonchilled 500 35.3

6 Qn1 Quebec nonchilled 675 39.2

w1b1<-subset(CO2, Treatment==\"chilled\")#先只考虑寒带植物 fit<-aov(uptake~(conc*Type)+Error(Plant/conc), data=w1b1)

summary(fit)

Error: Plant

Df Sum Sq Mean Sq F value Pr(>F)

Type 1 2667.2 2667.2 60.41 0.00148 ** Residuals 4 176.6 44.1 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Plant:conc

Df Sum Sq Mean Sq F value Pr(>F)

conc 1 888.6 888.6 215.46 0.000125 *** conc:Type 1 239.2 239.2 58.01 0.001595 ** Residuals 4 16.5 4.1 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within

Df Sum Sq Mean Sq F value Pr(>F) Residuals 30 869 28.97

说明：在0.01的显著水平下，主效应“类型”（p值=0.00148）和“浓度”（p值=0.000125）以及交叉效应“类型*浓度”（p值=0.001595）都非常显著。

attach(w1b1)

interaction.plot(conc, Type, uptake, type=\"b\ed\ant Type and Concentration\")

boxplot(uptake~Type*conc, data=w1b1, col=(c(\"gold\een\")), main=\"Chilled Quebec and Mississippi Plants\ = \"Carbon dioxide uptake rate (umol/m^2 sec)\") detach(w1b1)

可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着CO2浓度的升高，差异越来越明显。

注1：重复测量设计时，需要有长格式数据才能拟合模型，若是宽数据，需要用reshape2包中的melt()函数转化为长数据；

注2：这里是用的传统的重复测量方差分析，假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。

实际中，该假设一般不满足，可以尝试：

使用lme4包中的lmer()函数拟合线性混合模型（Bates，2005）；  使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正）；

 使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型（UCLA，2009）；

 用多元方差分析对重复测量数据进行建模（Hand，1987）。

四、多元方差分析

1. 因变量不只一个的方差分析，称为多元方差分析。

要求数据满足：多元正态性、方差—协方差矩阵同质性（即指各组的协方差矩阵相同，通常可用Box’s M检验来评估该假设）。 例4使用MASS包中的UScereal数据集，研究美国谷物中的卡路里、脂肪、糖含量是否会因为储物架位置的不同而发生变化。自变量货架位置shelf有1, 2, 3三个水平。

library(MASS) attach(UScereal)

y<-cbind(calories, fat, sugars) head(y)

calories fat sugars

[1,] 212.1212 3.030303 18.18182 [2,] 212.1212 3.030303 15.15151 [3,] 100.0000 0.000000 0.00000 [4,] 146.6667 2.666667 13.33333 [5,] 110.0000 0.000000 14.00000 [6,] 173.3333 2.666667 10.66667

aggregate(y, by=list(shelf), mean)

Group.1 calories fat sugars 1 1 119.4774 0.6621338 6.295493

2 2 129.8162 1.3413488 12.507670 3 3 180.1466 1.9449071 10.856821

cov(y)#求协方差矩阵，主对角线为方差

calories fat sugars

calories 35.24210 60.674383 180.380317 fat 60.67438 2.713399 3.995474 sugars 180.38032 3.995474 34.050018

library(mvnormtest)

mshapiro.test(t(y))#多元正态性检验，注意要对y转置 Shapiro-Wilk normality test data: Z

W = 0.6116, p-value = 7.726e-12 fit<-manova(y~shelf)#做多元方差分析 summary(fit)

DfPillaiapprox F numDf den DfPr(>F)

shelf 1 0.19594 4.955 3 61 0.00383 ** Residuals 63 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary.aov(fit) #输出每个单变量的方差分析结果

Response calories :

Df Sum Sq Mean Sq F value Pr(>F)

shelf 1 45313 45313 13.995 0.0003983 *** Residuals 63 203982 3238 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Response fat :

Df SumSq Mean Sq F value Pr(>F)

shelf 1 18.421 18.4214 7.476 0.008108 ** Residuals 63 155.236 2.41 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Response sugars :

Df SumSq Mean Sq F value Pr(>F)

shelf 1 183.34 183.34 5.787 0.01909 * Residuals 63 1995.87 31.68 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

说明：多元方差分析的p值=0.00383<0.005, 故有显著差异；每个单变量的方差分析结果也都是有显著差异的。

为了保证方差分析结果的准确性，还可以用mvoutlier包中的

aq.plot()函数，进行多元离群点检验（结果略）： library(mvoutlier) outliers<-aq.plot(y) outliers

2. 稳健的多元方差分析

若数据不满足多元正态性或方差-协方差矩阵均质性，或担心多元离群点，可以考虑用稳健或非参数版的MANOVA.

稳健单因素MANOVA可用rrcov包中的Wilks.test()函数实现。

Wilks.test(y, shelf, method=\"mcd\")#注意计算速度相当慢 Robust One-way MANOVA (Bartlett Chi2) data: x

Wilks' Lambda = 0.51073, Chi2-Value = 23.7110, DF = 4.9242, p-value = 0.0002299

sample estimates:

calories fat sugars 1 119.8210 0.7010828 5.663143 2 128.0407 1.1849576 12.537533 3 160.8604 1.6524559 10.3526

五、用回归来做ANOVA

ANOVA和回归都是广义线性模型的特例，所以，可以用lm()函数实现。

setwd(\"E:/办公资料/R语言/R语言学习系列/codes\") load(\"ex27_ex1.Rdata\")

fit.aov<-aov(wear~brand,data=datas) summary(fit.aov)

Df Sum Sq Mean Sq F value Pr(>F)

brand 3 0.7398 0.24660 24.55 3.15e-06 *** Residuals 16 0.1607 0.01005 ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.lm<-lm(wear~brand, data=datas) summary(fit.lm)

Call:

lm(formula = wear ~ brand, data = datas)

Residuals:

Min 1Q Median 3Q Max -0.1460 -0.0540 -0.0130 0.0385 0.1960

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 2.41000 0.04482 53.768< 2e-16 *** brandB -0.00600 0.06339 -0.095 0.9258 brandC -0.300 0.06339 -5.742 3.03e-05 *** brandD 0.16200 0.06339 2.556 0.0212 * ---

Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1002 on 16 degrees of freedom Multiple R-squared: 0.8215, Adjusted R-squared: 0.7881 F-statistic: 24.55 on 3 and 16 DF, p-value: 3.152e-06

说明：可以看到p值是一样的。

线性模型要求预测变量是数值型，当lm()函数遇到因子型自变量时，会用一系列与因子水平相对应的数值型对照变量来代替因子。若因子有k个水平，将会创建k-1个对照变量。

R语言中，对照变量的创建方法：

contr.helmert——第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值，以此类推；

contr.poly——基于正交多项式的对照，用于趋势分析（线性、二次、三次等）和等距水平的有序因子；

contr.sum——对照变量之和为0，也称为偏差找对，对各水平的均值与所有水平的均值进行比较；

contr.treatment——各水平对照基线水平（默认第一个水平），也

称作虚拟编码（虚拟变量、哑变量）；

contr.SAS——类似于contr.treatment，只是基线水平变成了最后一个水平，生成的系数类似于大部分SAS过程中使用的对照变量。

注1：默认将对照处理（contr.treatment）用于无序因子，正交多项式（contr.poly）用于有序因子；

注2：可以通过lm()函数的参数contrasts=\"…\"修改默认的对照方法。

contrasts(datas$brand)

#查看对照编码，默认第1水平（品牌A）作为参考基准 B C D A 0 0 0 B 1 0 0 C 0 1 0 D 0 0 1

若样本是品牌C，则变量C=1，其它变量B、D都=0，无需列出品牌A的变量值，因为它的4个变量都为0，就说明是品牌A.

回归模型的参数估计中的变量，brandB就表示品牌A与品牌B的一个对照，其余也是类似。