R语言实战(9)——方差分析

引言:上一章我们学习了如何通过量化的预测变量来预测量化的响应变量的回归模型。本期我们将一起学习如何对类别型预测变量建立合适的统计模型进行分析。

当包含的因子是解释变量时,关注的重点通常会从预测转向组别差异的分析,这种分析法称作方差分析(ANOVA)

9.1 ANOVA模型拟合

lm()函数也能分析ANOVA模型,但基本都使用aov()函数,两个函数的结果是等同的。

1. aov()函数

aov()函数的语法为 :

aov(formula, data=dataframe)

表9-1 常见研究设计的表达式

设计 表达式
单因素ANOVA y ~ A
含单个协变量的单因素 ANCOVA y ~ x + A
双因素 ANOVA y ~ A* B
含两个协变量的双因素 ANCOVA y ~ x1 + x2 + A * B
随机化区组 y ~ B + A(B是区组因子)
单因素组内 ANOVA y ~ A + Error(Subject/A)
含单个组内因子(W)和单个组间因子(B)的重复测量 ANOVA y ~ B * W + Error(Subject/W)

表中小写字母表示定量变量,大写字母表示组别因子,Subject是对被试者独有的标识变量。

2. 表达式中各项的顺序

当出现:(a)因子不止一个,并且是非平衡设计;(b)存在协变量

以上任意一种情况时,等式右边的变量都与其他每个变量相关,无法清晰地划分它们对因变量的影响。
对于双因素方差分析,若不同处理方式中的观测数不同,那么模型y ~ A * B与模型y ~ B * A的结果不同。

R默认类型I(序贯型)方法计算ANOVA效应。R中的ANOVA表的结果将评价:

  • A对y的影响;
  • 控制A时,B对y的影响;
  • 控制A和B的主效应时,A与B的交互效应

一般来说,越基础性的效应越需要放在表达式前面。

  1. 首先是协变量,然后是主效应,接着是双因素的交互项,再接着是三因素的交互项,以此类推。
  2. 对于主效应,越基础性的变量越应放在表达式前面,因此性别要放在处理方式之前。
  3. 基本准则:若研究设计不是正交的(也就是说,因子和/或协变量相关),一定要谨慎设置效应的顺序。

9.2 单因素方差分析

单因素方差分析(one-way ANOVA)比较分类因子定义的两个或多个组别中的因变量均值。

例:multcomp包中的cholesterol数据集
50个患者均接受降低胆固醇药物治疗(trt)五种疗法中的一种疗法。五种疗法为:相同药物分别20mg一天一次(1time)、10mg一天两次(2times)和5mg一天四次(4times)治疗及(drugD和drugE)候选药物治疗。

分析哪种药物疗法降低胆固醇(响应变量)最多:

library(multcomp)attach(cholesterol)#列出各组样本大小table(trt)          trt1time 2times 4times  drugD  drugE10     10     10     10     10#每10个患者接受其中一个药物疗法#计算各组均值      aggregate(response, by=list(trt), FUN=mean)   Group.1        x1   1time  5.781972  2times  9.224973  4times 12.374784   drugD 15.361175   drugE 20.94752#均值显示drugE降低胆固醇最多,而1time降低胆固醇最少#计算各组标准差aggregate(response, by=list(trt), FUN=sd)     Group.1        x1   1time 2.8781132  2times 3.4830543  4times 2.9231194   drugD 3.4546365   drugE 3.345003#各组的标准差相对恒定,在2.88到3.48间浮动#检验组间差异(ANOVA)fit <- aov(response ~ trt)   summary(fit)            Df Sum Sq Mean Sq F value    Pr(>F)trt          4 1351.4   337.8   32.43  9.82e-13 ***Residuals   45  468.8    10.4---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANOVA对治疗方式(trt)的F检验非常显著(p<0.0001),说明五种疗法的效果不同。

gplots包中的plotmeans()绘制带有置信区间的组均值图形:

library(gplots)plotmeans(response ~ trt, xlab="Treatment", ylab="Response",   #绘制各组均值及其置信区间的图形          main="Mean Plotnwith 95% CI")detach(cholesterol)

图9-1 五种降低胆固醇药物疗法的均值,含95%的置信区间

1. 多重比较

ANOVA对各疗法的F检验表明五种药物疗法效果不同,多重比较可以知道哪种疗法与其他疗法不同。

TukeyHSD()函数提供了对各组均值差异的成对检验:

TukeyHSD(fit)  Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = response ~ trt)$trt                  diff        lwr       upr     p adj2times-1time   3.44300 -0.6582817  7.544282 0.13809494times-1time   6.59281  2.4915283 10.694092 0.0003542drugD-1time    9.57920  5.4779183 13.680482 0.0000003drugE-1time   15.16555 11.0642683 19.266832 0.00000004times-2times  3.14981 -0.9514717  7.251092 0.2050382drugD-2times   6.13620  2.0349183 10.237482 0.0009611drugE-2times  11.72255  7.6212683 15.823832 0.0000000drugD-4times   2.98639 -1.1148917  7.087672 0.2512446drugE-4times   8.57274  4.4714583 12.674022 0.0000037drugE-drugD    5.58635  1.4850683  9.687632 0.0030633

1time和2times的均值差异不显著(p=0.138);1time和4times间的差异非常显著 (p<0.001)。

par(las=2)           #旋转轴标签par(mar=c(5,8,4,2))  #增大左边界的面积,可使标签摆放更美观plot(TukeyHSD(fit))

图9-2 Tukey HSD均值成对比较图

multcomp包中的glht()函数提供了既适用于线性模型,也适用于广义线性模型的多重均值比较方法。
重现Tukey HSD检验,并用一个不同的图形对结果进行展示:

library(multcomp)par(las=1)par(mar=c(5,4,6,2))  #增大顶部边界面积,适合字母阵列摆放tuk <- glht(fitlinfct=mcp(trt="Tukey"))plot(cld(tuklevel=.05),col="lightgrey")#cld()函数中的level选项设置显著水平(0.05,即本例中的95%的置信区间)

图9-3 multcomp包中的TukeyHSD检验

有相同字母的组(用箱线图表示)说明均值差异不显著。

2. 评估检验的假设条件

单因素方差分析中,假设因变量服从正态分布,各组方差相等。
可以使用Q-Q图来检验正态性假设

library(car)qqPlot(lm(response ~ trt, data=cholesterol),       simulate=TRUE, main="Q-Q Plot", labels=FALSE)

图9-4 正态性检验

注意qqPlot()要求用lm()拟合。

方差齐性检验:

可以通过如下代码来做Bartlett检验:

bartlett.test(response ~ trt, data=cholesterol)   Bartlett test of homogeneity of variancesdata:  response by trtBartlett's K-squared = 0.57975, df = 4, p-value = 0.9653

Bartlett检验表明五组的方差并没有显著不同(p=0.97)。

方差齐性分析对离群点非常敏感。
可利用car包中的outlierTest()函数来检测离群点

library(car)outlierTest(fit)No Studentized residuals with Bonferroni p < 0.05Largest |rstudent|:rstudent unadjusted p-value Bonferroni p19 2.251149           0.029422           NA

没有证据说明胆固醇数据中含有离群点(当p>1时将产生NA)。
根据Q-Q图、Bartlett检验和离群点检验,该数据可以用ANOVA模型拟合得很好。

9.3 单因素协方差分析

单因素协方差分析(ANCOVA)包含一个或多个定量的协变量。

例:multcomp包中的litter数据集
怀孕小鼠被分为四个小组,每个小组接受不同剂量(0、5、50或500)的药物处理产下幼崽的体重均值为因变量,怀孕时间为协变量。

library(multcomp)data(litter, package="multcomp")attach(litter)table(dose)dose0   5  50 50020  19  18  17      #每种剂量下所产的幼崽数并不相同aggregate(weight, by=list(dose), FUN=mean)Group.1        x1       0 32.308502       5 29.308423      50 29.866114     500 29.64647  #获得各组均值,可以发现未用药组幼崽体重均值最高(32.3)fit <- aov(weight ~ gesttime + dose)summary(fit)Df Sum Sq Mean Sq F value  Pr(>F)gesttime     1  134.3  134.30   8.049 0.00597 **dose         3  137.1   45.71   2.739 0.04988 *Residuals   69 1151.3   16.69---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANCOVA的F检验表明:(a)怀孕时间与幼崽出生体重相关;(b)控制怀孕时间,药物剂量与出生体重相关。

可使用effects包中的effects()函数获取调整的组均值,即去除协变量效应后的组均值:

library(effects)effect("dose", fit)dose effectdose       0        5       50      50032.35367 28.87672 30.56614 29.33460

同样可使用multcomp包来对所有均值进行成对比较,还可以用来检验用户自定义的均值假设:

library(multcomp)contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))   #c(3, -1, -1, -1)设定第一组和其他三组的均值进行比较summary(glht(fit, linfct=mcp(dose=contrast)))Simultaneous Tests for General Linear HypothesesMultiple Comparisons of Means: User-defined ContrastsFit: aov(formula = weight ~ gesttime + dose)Linear Hypotheses:Estimate Std. Error t value Pr(>|t|)no drug vs. drug == 0    8.284      3.209   2.581    0.012 *---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Adjusted p values reported -- single-step method)

假设检验的t统计量(2.581) 在p<0.05水平下显著,可以得出未用药组比其他用药条件下的出生体重高的结论。

1. 评估检验的假设条件

  • ANCOVA与ANOVA相同,都需要正态性和同方差性假设,可用上文提到的Q-Q图、Bartlett检验和离群点检验。
  • ANCOVA还假定回归斜率相同

例:multcomp包中的litter数据集
假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。ANCOVA模型包含怀孕时间×剂量的交互项时,可对回归斜率的同质性进行检验,交互效应若显著,则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。

检验回归斜率的同质性:

library(multcomp)fit2 <- aov(weight ~ gesttime*dose, data=litter)summary(fit2)Df Sum Sq Mean Sq F value  Pr(>F)gesttime       1  134.3  134.30   8.289 0.00537 **dose           3  137.1   45.71   2.821 0.04556 *gesttime:dose  3   81.9   27.29   1.684 0.17889Residuals     66 1069.4   16.20---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到交互效应不显著,支持了斜率相等的假设。
若假设不成立,可以尝试变换协变量或因变量,或使用能对每个斜率独立解释的模型,或使用不需要假设回归斜率同质性的非参数 ANCOVA方法。

2. 结果可视化

HH包中的ancova()函数可以绘制因变量、协变量和因子之间的关系图。

library(HH)ancova(weight ~ gesttime + dose, data=litter)Analysis of Variance Table  Response: weight           Df  Sum Sq Mean Sq F value   Pr(>F)     gesttime   1  134.30 134.304  8.0493 0.005971 **  dose       3  137.12  45.708  2.7394 0.049883 *   Residuals 69 1151.27  16.685                      ---  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

图9-5 四种药物处理组的怀孕时间和出生体重的关系图

用怀孕时间来预测出生体重的回归线相互平行,只是截距项不同。

若用ancova(weight ~ gesttime*dose),生成的图形将允许斜率和截距项依据组别而发生变化,可以可视化违背回归斜率同质性的实例。

ancova(weight ~ gesttime * dose, data=litter)

图9-6 四种药物处理组的怀孕时间和出生体重的关系图(斜率和截距项可改变)

9.4 双因素方差分析

在双因素方差分析中,受试者被分配到两因子的交叉类别组中。

例:基础安装中的ToothGrowth数据集
随机分配60只豚鼠,分别采用两种喂食方法(橙汁或维生素C),各喂食方法中抗坏血酸含量有三种水平(0.5mg、1mg或2mg),每种处理方式组合都被分配10只豚鼠。

attach(ToothGrowth)table(supp, dose)dosesupp 0.5  1  2OJ  10 10 10VC  10 10 10#获得各单元的均值aggregate(len, by=list(supp, dose), FUN=mean)    Group.1 Group.2     x1      OJ     0.5 13.232      VC     0.5  7.983      OJ     1.0 22.704      VC     1.0 16.775      OJ     2.0 26.066      VC     2.0 26.14    #获得各单元的标准差aggregate(len, by=list(supp, dose), FUN=sd)     Group.1 Group.2        x1      OJ     0.5 4.4597092      VC     0.5 2.7466343      OJ     1.0 3.9109534      VC     1.0 2.5153095      OJ     2.0 2.6550586      VC     2.0 4.797731#dose变量被转换为因子变量dose <- factor(dose)   #aov()函数会将它当做一个分组变量,而不是一个数值型协变量fit <- aov(len ~ supp*dose)       summary(fit)Df Sum Sq Mean Sq F value   Pr(>F)supp         1  205.4   205.4  15.572 0.000231 ***dose         2 2426.4  1213.2  92.000  < 2e-16 ***supp:dose    2  108.3    54.2   4.107 0.021860 *Residuals   54  712.1    13.2---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

用summary()函数得到方差分析表,可以看到主效应(supp和dose)和交互效应都非常显著。

结果可视化

用interaction.plot()函数展示双因素方差分析的交互效应:

attach(ToothGrowth)interaction.plot(dose, supp, len, type="b",                 col=c("red","blue"), pch=c(16, 18),                 main = "Interaction between Dose and Supplement Type")detach(ToothGrowth)

图9-7 喂食方法和剂量对牙齿生长的交互作用

用interaction.plot()函数绘制了牙齿长度的均值。

用gplots包中的plotmeans()函数展示交互效应:

library(gplots)attach(ToothGrowth)plotmeans(len ~ interaction(supp, dose, sep=" "),          connect=list(c(1,3,5),c(2,4,6)),          col=c("red", "darkgreen"),          main = "Interaction Plot with 95% CIs",          xlab="Treatment and Dose Combination")detach(ToothGrowth)

图9-8 喂食方法和剂量对牙齿生长的交互作用

用plotmeans()函数绘制的95%的置信区间的牙齿长度均值

用HH包中的interaction2wt()函数来可视化结果:

library(HH)attach(ToothGrowth)interaction2wt(len ~ supp * dose)     detach(ToothGrowth)

图9-9 ToothGrowth数据集的主效应和交互效应

推荐HH包中的interaction2wt()函数,能展示任意复杂度设计(双因素方差分析、三因素方差分析等)的主效应(箱线图)和交互效应。

9.5 重复测量方差分析

重复测量方差分析(单因素组内方差分析),即受试者被测量不止一次。

例:基础安装包中的CO2数据集
包含了北方和南方牧草类植物Echinochloa crus-galli的寒冷容忍度研究结果。比较某浓度二氧化碳的环境中寒带植物与非寒带植物的光合作用率。自变量是植物类型Type(组间因子)和七种水平的二氧化碳浓度conc(组内因子)。

含一个组间因子和一个组内因子的重复测量方差分析:

CO2$conc <- factor(CO2$conc)w1b1 <- subset(CO2, Treatment=='chilled')fit <- aov(uptake ~ conc*Type + Error(Plant/(conc)), w1b1)summary(fit)Error: PlantDf 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 ‘ ’ 1Error: Plant:conc     Df Sum Sq Mean Sq F value   Pr(>F)conc       6 1472.4  245.40   52.52 1.26e-12 ***conc:Type  6  428.8   71.47   15.30 3.75e-7 ***Residuals 24  112.1    4.67---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

方差分析表表明在0.01的水平下,主效应类型和浓度以及交叉效应类型×浓度都非常显著。

通过interaction.plot()函数展示了交互效应:

par(las=2)par(mar=c(10,4,4,2))with(w1b1, interaction.plot(conc,Type,uptake,type="b", col=c("red","blue"), pch=c(16,18),main="Interaction Plot for Plant Type and Concentration"))

图9-10 interaction.plot()函数绘制CO2浓度和植物类型对CO2吸收的交互影响

可以使用boxplot()函数对相同的数据画图,展示交互效应其他不同的侧面:

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

图9-11 boxplot()函数绘制CO2浓度和植物类型对CO2吸收的交互效应

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

9.6 多元方差分析

当因变量(结果变量)不止一个时,可用多元方差分析(MANOVA)对它们同时进行分析。

例:MASS包中的UScereal数据集。
研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化。

单因素多元方差分析:

library(MASS)attach(UScereal)#将shelf变量转换为因子变量,从而使它在后续分析中能作为分组变量shelf <- factor(shelf)   #cbind()函数将三个因变量合并成一个矩阵y <- cbind(calories, fat, sugars)  #获取货架的各个均值aggregate(y, by=list(shelf), FUN=mean)     Group.1 calories       fat    sugars1       1 119.4774 0.6621338  6.2954932       2 129.8162 1.3413488 12.5076703       3 180.1466 1.9449071 10.856821#cov()函数输出各谷物间的方差和协方差cov(y)         calories       fat     sugarscalories 3895.24210 60.674383 180.380317fat        60.67438  2.713399   3.995474sugars    180.38032  3.995474  34.050018#manova()函数能对组间差异进行多元检验fit <- manova(y ~ shelf)   summary(fit)Df Pillai approx F num Df den Df    Pr(>F)shelf      2 0.4021   5.1167      6    122 0.0001015 ***Residuals 62---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

F值显著,说明三个组的营养成分测量值不同。

#summary.aov()函数对每一个变量做单因素方差分析summary.aov(fit)                              Response calories :Df Sum Sq Mean Sq F value    Pr(>F)shelf        2  50435 25217.6  7.8623 0.0009054 ***Residuals   62 198860  3207.4---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response fat :Df Sum Sq Mean Sq F value  Pr(>F)shelf        2  18.44  9.2199  3.6828 0.03081 *Residuals   62 155.22  2.5035---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response sugars :Df  Sum Sq Mean Sq F value   Pr(>F)shelf        2  381.33 190.667  6.5752 0.002572 **Residuals   62 1797.87  28.998---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

三组中每种营养成分的测量值都是不同的。

1. 评估假设检验

单因素多元方差分析有两个前提假设:

  • 多元正态性,指因变量组合成的向量服从一个多元正态分布,可用Q-Q图来检验。
  • 方差-协方差矩阵同质性,R中还没有好的方法检验。

检验多元正态性:

center <- colMeans(y)n <- nrow(y)p <- ncol(y)cov <- cov(y)d <- mahalanobis(y,center,cov)coord <- qqplot(qchisq(ppoints(n),df=p),                d, main="Q-Q Plot Assessing Multivariate Normality",                ylab="Mahalanobis D2")abline(a=0,b=1)identify(coord$x, coord$y, abels=row.names(UScereal))      #identify()函数交互性地对图中的点进行鉴别

图9-12 检验多元正态性的Q-Q图

观测点“Wheaties Honey Gold”和“Wheaties”异常,数据集似乎违反了多元正态性,可以删除这两个点再重新分析。

可以使用mvoutlier包中的ap.plot()函数来检验多元离群点:

library(mvoutlier)outliers <- aq.plot(y)Projection to the first and second robust principal components.  Proportion of total variation (explained variance): 0.9789888outliers$outliers    [1FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE[16FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE[31]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE[46FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE[61FALSE FALSE FALSE FALSE FALSE

图9-13 检验多元离群点

2. 稳健多元方差分析

如果多元正态性或者方差-协方差均值假设都不满足,或者担心多元离群点,那么可以考虑用稳健或非参数版本的MANOVA检验。

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

library(rrcov)Wilks.test(y, shelf, method="mcd")    Robust One-way MANOVA (Bartlett Chi2)data:  xWilks' Lambda = 0.511, Chi2-Value = 23.96, DF = 4.98, p-value = 0.0002167sample estimates:  calories    fat  sugars1      120  0.701    5.662      128  1.185   12.543      161  1.652   10.35

稳健检验对离群点和违反MANOVA假设的情况不敏感,验证了存储在货架顶部、中部和底部的谷物营养成分含量不同。

9.7 用回归来做ANOVA

例:multcomp包中的cholesterol数据集
比较五种降低胆固醇药物疗法(trt)的影响:

library(multcomp)levels(cholesterol$trt)[1] "1time"  "2times" "4times" "drugD"  "drugE"fit.aov <- aov(response ~ trt, data=cholesterol)       summary(fit.aov)  Df   Sum Sq  Mean Sq  F value     Pr(>F)trt          4  1351.37   337.84   32.433  9.819e-13 ***Residuals   45   468.75    10.42---  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

用lm()函数拟合同样的模型:

fit.lm <- lm(response ~ trt, data=cholesterol)summary(fit.lm)Call:  lm(formula = response ~ trt, data = cholesterol)Residuals:    Min      1Q  Median      3Q     Max-6.5418 -1.9672 -0.0016  1.8901  6.6008Coefficients:            Estimate Std. Error t value Pr(>|t|)(Intercept)    5.782      1.021   5.665 9.78e-07 ***trt2times      3.443      1.443   2.385   0.0213 *    #变量trt2times表示水平1time和2times的一个对照trt4times      6.593      1.443   4.568 3.82e-05 ***  #trt4times是1time和4times的一个对照trtdrugD       9.579      1.443   6.637 3.53e-08 ***trtdrugE      15.166      1.443  10.507 1.08e-13 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 3.227 on 45 degrees of freedomMultiple R-squared:  0.7425,  Adjusted R-squared:  0.7196F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13

从输出的概率值来看,各药物条件与第一组(1time)显著不同。

因为线性模型要求预测变量是数值型,lm()函数会用一系列与因子水平相对应的数值型对照变量来代替因子。默认情况下,对照处理用于无序因子,正交多项式用于有序因子。

表9-5 R提供的五种创建对照变量的内置方法

对照变量创建方法 描述
contr.helmert 第二个水平对照第一个水平,第三个水平对照前两个的均值,第四个水平对照前三个的均值,以此类推
contr.poly 基于正交多项式的对照,用于趋势分析(线性、二次、三次等)和等距水平的有序因子
contr.sum 对照变量之和限制为0。也称作偏差找对,对各水平的均值与所有水平的均值进行比较
contr.treatment 各水平对照基线水平(默认第一个水平)。也称作虚拟编码
contr.SAS 类似于contr.treatment,只是基线水平变成了最后一个水平。生成的系数类似于大部分SAS过程中使用的对照变量

以对照(treatment contrast)为例,因子的第一个水平变成了参考组,随后的变量都以它为标准。
可以通过contrasts()函数查看它的编码过程:

contrasts(cholesterol$trt)2times 4times drugD drugE1time       0      0     0     02times      1      0     0     04times      0      1     0     0drugD       0      0     1     0drugE       0      0     0     1

若患者处于drugD条件下,变量drugD等于1,其他变量2times、4times和drugE都等于0;无需列出第一组的变量值,因为其他四个变量都为0,这已经说明患者处于1time条件。

通过设定contrasts选项,可以修改lm()中默认的对照方法:

fit.lm <- lm(response ~ trtdata=cholesterol, contrasts"contr.helmert")

还能通过options()函数修改R会话中的默认对照方法:

options(contrasts = c("contr.SAS""contr.helmert"))

设定无序因子的默认对比方法为contr.SAS,有序因子的默认对比方法为contr.helmert。

再查看它的编码过程

contrasts(cholesterol$trt)         1time 2times 4times drugD  1time      1      0      0     02times     0      1      0     04times     0      0      1     0drugD      0      0      0     1drugE      0      0      0     0

9.8 小结

本章我们通过组内和组间设计的示例介绍了方差分析模型以及假设检验。学习完回归和方差分析这两种常用的统计方法后,下一章我们将学习对研究设计非常重要的功效分析。敬请期待~

生物信息学

利用Word和Excel批量获取Oncomine原始数据

2019-12-17 21:19:41

生物信息学

​R语言实战(10)——功效分析

2019-12-17 21:44:16

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索