R语言析因设计分析：线性模型中的对比

对比度可用于对线性模型中的处理进行比较。

常见的用途是使用析因设计时，除析因设计外还使用控制或检查处理。在下面的第一个示例中，有两个级别（1和2）的两个处理（D和C），然后有一个对照处理。此处使用的方法是方差的单向分析，然后使用对比来检验各种假设。

可下载资源

在下面的第二个示例中，对六种葡萄酒进行了测量，其中一些是红色，而有些是白色。我们可以比较的治疗中通过设置对比，并进行F检验红酒组。这类似于测试红酒的主要效果。

使用的软件包

如果尚未安装这些软件包，则以下命令将安装它们：

if(!require(car)){install.packages("car")}
if(!require(lsmeans){install.packages("lsmeans")}
if(!require(multcomp)){install.packages("multcomp")}

单自由度对比示例

这个假设的例子可以代表一项采用阶乘设计的实验，其中两个处理（D和C）分别处于两个级别（1 和2），并且是对照处理。

Data = read.table(textConnection(Input),header=TRUE)



Data$Treatment = factor(Data$Treatment, 
                        levels=unique(Data$Treatment))

Data

boxplot(Response ~ Treatment,
        data = Data,
        ylab="Response",
        xlab="Treatment")







###  Define linear model

model = lm(Response ~ Treatment, 
           data = Data)

library(car)

Anova(model, type="II")

summary(model)

lsmeans示例




 contrast           estimate        SE df t.ratio p.value
 D1vsD2          -0.83333333 0.1549193 10  -5.379  0.0031
 C1vsC2          -2.10000000 0.1549193 10 -13.555  <.0001
 InteractionDC    0.03333333 0.1549193 10   0.215  1.0000
 C1vsC2forD1only -1.03333333 0.1095445 10  -9.433  <.0001
 C1vsC2forD2only -1.06666667 0.1095445 10  -9.737  <.0001
 TreatsvsControl  3.96666667 0.3464102 10  11.451  <.0001
 T1vsC            0.26666667 0.1095445 10   2.434  0.3011
 T2vsC            1.30000000 0.1095445 10  11.867  <.0001
 T3vsC            0.66666667 0.1095445 10   6.086  0.0012
 T4vsC            1.73333333 0.1095445 10  15.823  <.0001

由于调整方法不同，p值与multcomp ###的p值略有不同。
   ###两个过程的调整方法，
   ### p值和其他统计信息，将是相同的。

   ###使用
   Adjust =“ none”，结果将与### aov方法相同。

multcomp示例






                     Estimate Std. Error t value Pr(>|t|)   

D1vsD2 == 0          -0.83333    0.15492  -5.379  0.00218 **

C1vsC2 == 0          -2.10000    0.15492 -13.555  < 0.001 ***

InteractionDC == 0    0.03333    0.15492   0.215  0.99938   

C1vsC2forD1only == 0 -1.03333    0.10954  -9.433  < 0.001 ***

C1vsC2forD2only == 0 -1.06667    0.10954  -9.737  < 0.001 ***

TreatsvsControl == 0  3.96667    0.34641  11.451  < 0.001 ***

T1vsC == 0            0.26667    0.10954   2.434  0.17428   

T2vsC == 0            1.30000    0.10954  11.867  < 0.001 ***

T3vsC == 0            0.66667    0.10954   6.086  < 0.001 ***

T4vsC == 0            1.73333    0.10954  15.823  < 0.001 ***

一组治疗中的全局F检验示例

该示例具有由三种红酒和三种白葡萄酒组成的处理。我们将想知道红酒组中的处理是否对响应变量有影响。这种方法之所以具有优势，是因为仍可以在红酒中进行事后比较。



boxplot(Response ~ Treatment,
        data = Data,
        ylab="Response",
        xlab="Treatment")

与lsmeans 的对比测试

问题：红酒中有功效吗？



test(Test, joint=TRUE)


df1 df2    F p.value
  2  12 24.3  0.0001

使用2个自由度进行了一次假设检验。这调查了
### 3组治疗的效果。

###结果与multcomp的结果基本相同

问题：白葡萄酒有效果吗？



test(Test, joint=TRUE)

df1 df2   F p.value
  2  12 0.3  0.7462

两行对比
使用2个自由度进行了一次假设检验。本研究调查了
###一组3种治疗方法中的效果

###结果与multcomp的结果相同

问题：红葡萄酒和白葡萄酒之间有区别吗？而且，红酒的平均分离度



contrast        estimate       SE df t.ratio p.value

 Red_vs_white          21 1.490712 12  14.087  <.0001

 Merlot_vs_Cab         -3 0.860663 12  -3.486  0.0179
 Cab_vs_Syrah          -3 0.860663 12  -3.486  0.0179
 Syrah_vs_Merlot        6 0.860663 12   6.971  0.0001

请注意，p值是
由于调整方法不同，因此与multcomp 不同。

使用Multcomp 进行对比测试

问题：红酒中有功效吗？



Global Test:
      F  DF1  DF2     Pr(>F)
1  24.3    2   12  6.029e-05

问题：白葡萄酒有效果吗？



Global Test:
    F DF1 DF2 Pr(>F)
1 0.3   2  12 0.7462

问题：红葡萄酒和白葡萄酒之间有区别吗？



   ### Adjustment options: "none", "single-step", "Shaffer", 
   ###                     "Westfall", "free", "holm", "hochberg", 
   ###                     "hommel", "bonferroni", "BH", "BY", "fdr"


Linear Hypotheses:
                     Estimate Std. Error t value Pr(>|t|)    
Red_vs_white == 0     21.0000     1.4907  14.087   <0.001 ***

Merlot_vs_Cab == 0    -3.0000     0.8607  -3.486   0.0157 *  
Cab_vs_Syrah == 0     -3.0000     0.8607  -3.486   0.0156 *  
Syrah_vs_Merlot == 0   6.0000     0.8607   6.971   <0.001 ***

(Adjusted p values reported -- single-step method)

###使用test = adjusted（“ none”），结果将与下面的aov方法相同。

aov内的对比测试

在方差分析中使用单自由度对比的另一种方法是在摘要函数中使用split选项进行aov分析。



boxplot(Response ~ Treatment,
        data = Data,
        ylab="Response",
        xlab="Treatment")


   0.018





###  Define contrasts

D1vsD2 =          c(1,  1, -1, -1,  0)
C1vsC2 =          c(1, -1,  1, -1,  0)
InteractionDC =   c(1, -1, -1,  1,  0)
TreatsvsControl = c(1,  1,  1,  1, -4)



                             Df Sum Sq Mean Sq F value   Pr(>F)    
Treatment                     4  6.189   1.547  85.963 1.06e-07 ***
  Treatment: D1vsD2           1  0.521   0.521  28.935  0.00031 ***
  Treatment: C1vsC2           1  3.307   3.307 183.750 9.21e-08 ***
  Treatment: InteractionDC    1  0.001   0.001   0.046  0.83396    
  Treatment: TreatsvsControl  1  2.360   2.360 131.120 4.53e-07 ***
Residuals                    10  0.180