12.3　lmPerm包的置换检验

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版，能够进行置换检验，而非正态理论检验。

lmp()和aovp()函数的参数与lm()和aov()函数类似，只额外添加了perm =参数。perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样，直至估计的标准差在估计的p值0.1之下，判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意，若观测数大于10，perm = "Exact"将自动默认转为perm = "Prob"，因为精确检验只适用于小样本问题。

为深入了解函数的工作原理，我们将对简单回归、多项式回归、多元回归、单因素方差分析、单因素协方差分析和双因素因子设计进行置换检验。

12.3.1　简单回归和多项式回归

第8章中，我们使用线性回归研究了15名女性的身高和体重间的关系。用lmp()代替lm()，可获得代码清单12-2中置换检验的结果。

代码清单12-2　简单线性回归的置换检验

  1. > library(lmPerm)
  2. > set.seed(1234)
  3. > fit <- lmp(weight~height, data=women, perm="Prob")
  4. [1] "Settings:  unique SS : numeric variables centered"
  5. > summary(fit)
  7. Call:
  8. lmp(formula = weight ~ height, data = women, perm = "Prob")
  10. Residuals:
  11.    Min     1Q Median      3Q   Max
  12. -1.733 -1.133 -0.383  0.742  3.117
  14. Coefficients:
  15.        Estimate Iter Pr(Prob)
  16. height     3.45 5000   <2e-16 ***
  17. ---
  18. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  20. Residual standard error: 1.5 on 13 degrees of freedom
  21. Multiple R-Squared: 0.991,      Adjusted R-squared: 0.99
  22. F-statistic: 1.43e+03 on 1 and 13 DF, p-value: 1.09e-14

要拟合二次方程，可使用代码清单12-3中的代码。

代码清单12-3　多项式回归的置换检验

  1. > library(lmPerm)
  2. > set.seed(1234)
  3. > fit <- lmp(weight~height + I(height^2), data=women, perm="Prob")
  4. [1] "Settings:  unique SS : numeric variables centered"
  5. > summary(fit)
  7. Call:
  8. lmp(formula = weight ~ height + I(height^2), data = women, perm = "Prob")
  10. Residuals:
  11.     Min      1Q  Median      3Q     Max
  12. -0.5094 -0.2961 -0.0094  0.2862  0.5971
  14. Coefficients:
  15.             Estimate Iter Pr(Prob)
  16. height       -7.3483 5000   <2e-16 ***
  17. I(height^2)   0.0831 5000   <2e-16 ***
  18. ---
  19. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  21. Residual standard error: 0.38 on 12 degrees of freedom
  22. Multiple R-Squared: 0.999,      Adjusted R-squared: 0.999
  23. F-statistic: 1.14e+04 on 2 and 12 DF, p-value: <2e-16

可以看到，用置换检验来检验这些回归是非常容易的，修改一点代码即可。输出结果也与lm()函数非常相似。值得注意的是，增添的Iter栏列出了要达到判停准则所需的迭代次数。

12.3.2　多元回归

在第8章，多元回归被用来通过美国50个州的人口数、文盲率、收入水平和结霜天数预测犯罪率。将lmp()函数应用到此问题，结果参见代码清单12-4。

代码清单12-4　多元回归的置换检验

  1. > library(lmPerm)
  2. > set.seed(1234)
  3. > states <- as.data.frame(state.x77)
  4. > fit <- lmp(Murder~Population + Illiteracy+Income+Frost,
  5.              data=states, perm="Prob")
  6. [1] "Settings:  unique SS : numeric variables centered"
  7. > summary(fit)
  9. Call:
  10. lmp(formula = Murder ~ Population + Illiteracy + Income + Frost,
  11.     data = states, perm = "Prob")
  13. Residuals:
  14.      Min       1Q   Median        3Q     Max
  15. -4.79597 -1.64946 -0.08112  1.48150  7.62104
  17. Coefficients:
  18.             Estimate Iter Pr(Prob)
  19. Population 2.237e-04   51   1.0000
  20. Illiteracy 4.143e+00 5000   0.0004 ***
  21. Income     6.442e-05   51   1.0000
  22. Frost      5.813e-04   51   0.8627
  23. ---
  24. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '. ' 0.1 ' ' 1
  26. Residual standard error: 2.535 on 45 degrees of freedom
  27. Multiple R-Squared: 0.567,      Adjusted R-squared: 0.5285
  28. F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

回顾第8章，正态理论中Population和Illiteracy均显著（P < 0.05）。而该置换检验中，Population不再显著。当两种方法所得结果不一致时，你需要更加谨慎地审视数据，这很可能是因为违反了正态性假设或者存在离群点。

12.3.3　单因素方差分析和协方差分析

第9章中任意一种方差分析设计都可进行置换检验。首先，让我们看看9.1节中的单因素方差分析问题——各种疗法对降低胆固醇的影响。代码和结果见代码清单12-5。

代码清单12-5　单因素方差分析的置换检验

  1. > library(lmPerm)
  2. > library(multcomp)
  3. > set.seed(1234)
  4. > fit <- aovp(response~trt, data=cholesterol, perm="Prob")
  5. [1] "Settings:  nique SS "
  6. > summary(fit)
  7. Component 1 :
  8.             Df R Sum Sq R Mean Sq Iter  Pr(Prob)
  9. trt          4  1351.37    337.84 5000 < 2.2e-16 ***
  10. Residuals   45   468.75     10.42
  11. ---
  12. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '. ' 0.1 ' ' 1

结果表明各疗法的效果不全相同。

协方差分析的置换检验取自第9章的问题——当控制妊娠期时间相同时，观测四种药物剂量对鼠崽体重的影响。代码和结果参见代码清单12-6。

代码清单12-6　单因素协方差分析的置换检验

  1. > library(lmPerm)
  2. > set.seed(1234)
  3. > fit <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob")
  4. [1] "Settings:  unique SS : numeric variables centered"
  5. > summary(fit)
  6. Component 1 :
  7.             Df R Sum Sq R Mean Sq Iter Pr(Prob)
  8. gesttime     1   161.49   161.493 5000   0.0006 ***
  9. dose         3   137.12    45.708 5000   0.0392 *
  10. Residuals   69  1151.27    16.685
  11. ---
  12. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

依据p值可知，当控制妊娠期时间相同时，四种药物剂量对鼠崽的体重影响不相同。

12.3.4　双因素方差分析

本节最后，我们对析因实验设计进行置换检验。以第9章中的维生素C对豚鼠牙齿生长的影响为例，该实验两个可操作的因子是剂量（三水平）和喂食方式（两水平）。10只豚鼠分别被分配到每种处理组合中，形成一个3 × 2的析因实验设计。置换检验结果参见代码清单12-7。

代码清单12-7　双因素方差分析的置换检验

  1. > library(lmPerm)
  2. > set.seed(1234)
  3. > fit <- aovp(len~supp*dose, data=ToothGrowth, perm="Prob")
  4. [1] "Settings:  unique SS : numeric variables centered"
  5. > summary(fit)
  6. Component 1 :
  7.             Df R Sum Sq R Mean Sq Iter Pr(Prob)
  8. supp         1   205.35    205.35 5000  < 2e-16 ***
  9. dose         1  2224.30   2224.30 5000  < 2e-16 ***
  10. supp:dose    1    88.92     88.92 2032  0.04724 *
  11. Residuals   56   933.63     16.67
  12. ---
  13. Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

在0.05的显著性水平下，三种效应都不等于0，在0.01的水平下，只有主效应显著。

值得注意的是，当将aovp()应用到方差分析设计中时，它默认使用唯一平方和法（SAS也称为类型III平方和）。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和（SAS是类型I平方和）。每种效应依据模型中先出现的效应做相应调整。对于平衡设计，两种方法结果相同，但是对于每个单元格观测数不同的不平衡设计，两种方法结果则不同。不平衡性越大，结果分歧越大。若在aovp()函数中设定seqs = TRUE，可以生成你想要的序贯平方和。关于类型I和类型III的平方和的更多细节，请参考9.2节。