上一次给大家介绍了如何用R语言进行主成分分析,今天介绍的主角也是PCA的好朋友噢,掌声欢迎我们的第二位小伙伴——冗余分析(RDA)。
1 冗余分析 简介
冗余分析(Redundancy Analysis,RDA),是一种回归分析结合主成分分析的排序方法。RDA建模[1]的大致思想是先将响应变量矩阵与解释变量之间进行多元线性回归,再对得到的拟合值进行主成分分析。
2 计算步骤
数据预处理:如果响应变量或者解释变量具有不同的测量单位,可以进行标准化处理。
符号说明:
:标准化后的解释变量矩阵, 为第 个解释变量。
:标准化后的响应变量矩阵, 为第 个响应变量。
step 1:将 中的每个响应变量分别与 进行多元回归,获得对应的响应变量的拟合值向量 和残差向量 , 构成拟合值矩阵 ;
step 2:对 进行PCA分析,将得到典型特征根向量矩阵 ;
step 3:计算样方得分 和样方约束 ;
step 4:对 进行PCA分析;
3 R语言实战
R语言中为我们提供了可直接用来进行简单冗余分析的函数,通过下载相应的程序包就可以使用,本文中使用的是easyCODA
程序包,这个程序包也包含了许多实用的函数,可以在Rstudio中help一下去学习具体函数用法。
在R语言的帮助页面里,使用的是fish数据集对RDA()
进行说明。
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转换为0\1哑变量矩阵 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()查看冗余分析输出的对象列表。
> 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箭头相反则存在负相关关系。
4 结语
冗余分析在生物统计中应用较多,概念比较难懂,本文中也只是对RDA做出了一个简短的解释,想进行更深入的学习可以参考下述资料:
R语言实现冗余分析完整代码[2]
数量统计学生态笔记||冗余分析[3]