冗余分析_ 字节宝

上一次给大家介绍了如何用R语言进行主成分分析，今天介绍的主角也是PCA的好朋友噢，掌声欢迎我们的第二位小伙伴——冗余分析（RDA）。

1 冗余分析简介

冗余分析（Redundancy Analysis,RDA），是一种回归分析结合主成分分析的排序方法。RDA建模[1]的大致思想是先将响应变量矩阵与解释变量之间进行多元线性回归，再对得到的拟合值进行主成分分析。

2 计算步骤

数据预处理：如果响应变量或者解释变量具有不同的测量单位，可以进行标准化处理。

符号说明：

{X}

：标准化后的解释变量矩阵，

{X_i}

为第

{i}

个解释变量。

{Y}

：标准化后的响应变量矩阵，

{Y_j}

为第

{j}

个响应变量。

step 1：将

{Y}

中的每个响应变量分别与

{X}

进行多元回归，获得对应的响应变量的拟合值向量

hat{y_i}

和残差向量

{y_{res}}

，

hat{y_i}

构成拟合值矩阵

hat{Y}

;

step 2：对

{Y}

进行PCA分析，将得到典型特征根向量矩阵

{U}

；

step 3：计算样方得分

{YU}

和样方约束

hat{Y}

{U}

;

step 4：对

{y_{res}}

进行PCA分析；

3 R语言实战

R语言中为我们提供了可直接用来进行简单冗余分析的函数，通过下载相应的程序包就可以使用，本文中使用的是easyCODA程序包，这个程序包也包含了许多实用的函数，可以在Rstudio中help一下去学习具体函数用法。

在R语言的帮助页面里，使用的是fish数据集对RDA() 进行说明。

代码语言：javascript复制

install.packages("easyCODA")
library(easyCODA)
data(fish)
sex     <- fish[,1]
habitat <- fish[,2]
mass    <- fish[,3]
# 将第4个至第29个描述鱼的形态变量的数据存放在fishm矩阵中
fishm   <- as.matrix(fish[,4:29])
# 对fishm矩阵中的每一行数据进行中心化处理
fishm   <- fishm / apply(fishm, 1, sum)
# fishm包含了每种鱼的26个形态指标经过中心化处理后的数据【【
> head(fishm)
           Hw         Bg         Bd        Bcw         Jw         Jl
FP 0.02682445 0.03157327 0.03388351 0.01142285 0.01244962 0.03632210
ML 0.02662625 0.03318040 0.03413622 0.01146977 0.01174286 0.03632094
ML 0.02651337 0.03352411 0.03658335 0.01402149 0.01376655 0.03518120
ML 0.02787768 0.03439454 0.03729092 0.01243031 0.01351645 0.03572204
FL 0.02585009 0.03201063 0.03430573 0.01232107 0.01377061 0.03744640
ML 0.02920211 0.03570563 0.03927620 0.01326209 0.01479233 0.03672579
           Bp        Bac        Bch         Fc        Fdw        Faw
FP 0.10696409 0.04271376 0.01918782 0.04385605 0.03321611 0.02692712
ML 0.10321427 0.04488230 0.02087771 0.03950243 0.03442296 0.02597083
ML 0.10395024 0.04359409 0.02099400 0.03783253 0.03358784 0.02267658
ML 0.10307499 0.04222683 0.02185562 0.04045280 0.03054476 0.02525886
FL 0.09964366 0.03913752 0.02132029 0.03964486 0.03558616 0.02610376
ML 0.09907038 0.04157156 0.02341269 0.03902116 0.03390760 0.02566980
            Bc         Fp        Fpl        Fal        Fdl         Hh
FP 0.011089150 0.04697487 0.03357548 0.03019996 0.03587289 0.02768437
ML 0.010391065 0.04392648 0.03128243 0.02810093 0.03581572 0.03199246
ML 0.011280927 0.04625817 0.03193076 0.02749487 0.03552536 0.03157385
ML 0.009799426 0.04374744 0.03276532 0.02994135 0.03470831 0.03064131
FL 0.009434076 0.04620402 0.03326690 0.03310986 0.03422118 0.02961889
ML 0.010571418 0.04537166 0.02966118 0.02820745 0.03596067 0.02939339
           Hg        Ba         Jm        Hal        Hpl         Ed
FP 0.04327849 0.1319788 0.02869831 0.02868547 0.01333522 0.01259081
ML 0.04338030 0.1356846 0.02631220 0.03694904 0.01097821 0.01339505
ML 0.04448637 0.1357153 0.02619469 0.03473506 0.01074556 0.01274681
ML 0.04655934 0.1321836 0.02589848 0.03468417 0.01175448 0.01235790
FL 0.04595035 0.1352781 0.02833847 0.03339977 0.01303376 0.01326327
ML 0.04782004 0.1364720 0.02554228 0.03013300 0.01124727 0.01369566
           Hs         Hl
FP 0.06667608 0.06401930
ML 0.06619695 0.06324758
ML 0.06577354 0.06331341
ML 0.06667713 0.06363592
FL 0.06497554 0.06276499
ML 0.06363254 0.06067407
# 对质量进行对数变换
logmass <- log(mass)
# 计算性别和栖息地的组合代表数
sexhab  <- 2*(sex-1) habitat
sexhab.names <- c("FL","FP","ML","MP")
# 对sexhab进行命名
rownames(fishm) <- sexhab.names[sexhab]
# 将sexhab转换为01哑变量矩阵
sexhab.Z <- DUMMY(sexhab, catnames=sexhab.names)
vars     <- cbind(logmass, sexhab.Z)
# 对fishm数据进行对数中心化变化
require(ca)
fish.RDA <- RDA(CLR(fishm), cov=vars)
# 绘制冗余分析排序图
PLOT.RDA(fish.RDA, map="contribution", rescale=0.05, indcat=2:5, colrows=rainbow(4, start=0.1, end=0.8)[sexhab], cexs=c(0.8,0.8,1))

使用函数RDA()得到的分析结果包含了更多的信息:ChiDist是列联表的卡方检验结果；Inertia是特征根；Dim. 1、Dim. 2、Dim. 3、Dim. 4是提取的四个约束轴。

还可以可通过names()查看冗余分析输出的对象列表。

代码语言：javascript复制

> names(fish.RDA)
 [1] "sv"         "nd"         "rownames"   "rowmass"    "rowdist"   
 [6] "rowinertia" "rowcoord"   "rowsup"     "colnames"   "colmass"   
[11] "coldist"    "colinertia" "colcoord"   "colsup"     "N"         
[16] "call"       "rowpcoord"  "colpcoord"  "covcoord"   "covnames"  
[21] "cov"

在RDA的排序图中，线条长度代表该影响因子所占的比重，各因子之间的相关性可以根据夹角来判断。例如，Hh与Fdl存在正相关关系，Fdl与Hg近似正交可视为不存在相关关系，而Hg与Fdw箭头相反则存在负相关关系。