最近赖江山老师发布了一个R包: 原创R包：rdaenvpart（层次分割获取RDA和CCA单解释变量的贡献） http://wap.sciencenet.cn/home.php?mod=space&uid=267448&do=blog&id=1199315
将Hierarchical Partitioning与rda(cca)相结合，可得到单个环境因子的解释度。 看了一下源代码，其核心的Hierarchical Partitioning用的是我之前介绍过的hier.part包。
以下对rdaenvpart的过程进行一下探索
代码语言：javascript复制 1#主函数rdaenvpart
2#其中包含了另外四个子函数
3rdaenvpart=function (Y, X,type="RDA", pieplot = "tv")
4{
5#环境因子个数
6 Env.num <- dim(X)[2]
7#大于13个环境因子会报错
8 if (Env.num > 13)
9 stop("Number of variables must be < 13 for current implementation",call. = FALSE)
10 else {
11#allR2是一个子函数，其中又包含了2个子函数
12 gfs <- allR2(Y, X,type)
13#partition.rda是一个子函数，得到每个因子的解释度。
14 HP <- partition.rda(gfs, Env.num, var.names = names(data.frame(X)))
15#tv是展示每个解释变量占总变化量（total variation）的比例的饼图
16 if (pieplot=="tv") {
17 lbls<- c("unexplained",row.names(HP$I.perc))
18 pct <- round(c(1-sum(HP$IJ$I),HP$IJ$I)*100,1)
19 lbls <- paste(lbls, pct) # add percents to labels
20 lbls <- paste(lbls,"%",sep="") # ad % to labels
21#画图
22 pie(pct,labels = lbls,
23 main="individual % on total variation" )
24
25 }
26#tev是展示每个解释变量占总被解释变化量（total explained variation）的比例的饼图
27 if (pieplot=="tev") {
28 lbls<- row.names(HP$I.perc)
29 pct <- round(HP$I.perc$I,1)
30 lbls <- paste(lbls, pct) # add percents to labels
31 lbls <- paste(lbls,"%",sep="") # ad % to labels
32#画图
33 pie(pct,labels = lbls,
34 main="individual % on total explained variation" )
35
36 }
37#计算R2及adjusted R2
38 if(type=="RDA")
39 rsq <- RsquareAdj(rda(Y~., data = X))
40 if(type=="CCA")
41 rsq <- RsquareAdj(cca(Y~., data = X))
42 list(R2=rsq$r.squared,gfs = gfs, IJ.R2 = HP$IJ, I.perc = HP$I.perc,adj.R2=rsq$adj.r.squared,IJ.adjR2=HP$I.perc*0.01*rsq$adj.r.squared)
43 }
44}allR2,这里面还嵌套了两个函数，先往下拉看他们是干啥的。之后就知道allR2是得到环境因子的各种组合与OTU做cca或rda后得到的R2，并作为最后的gfs(拟合的优异度)输出。
代码语言：javascript复制 1allR2=function (Y,X,type)
2{
3 Env.num <- dim(X)[2]
4 n <- (2^Env.num) - 1
5 combs <- hpmatrix(Env.num)
6 gfs <- 0
7
8 for (i in 1:n) {
9 current.comb <- as.vector(combs[i, ][combs[i, ] > 0])
10 combn <- paste(names(data.frame(X)[current.comb]),
11 "", collapse = "")
12 new.line <- R2current(Y, current.comb, X,type)
13 gfs <- c(gfs, new.line)
14 }
15 gfs
16}hpmatrix，从环境因子中依次挑出来2，3，4，…n个因子，得到每种条件下可能的组合情况及组合数量。最后在汇总为所有环境因子组合的可能，并输出到allR2中的combs。
代码语言：javascript复制1hpmatrix=function (n)
2{ require(gtools)
3 x <- cbind(1:n, matrix(0, n, n - 1))
4 for (i in 2:n) {
5 nc <- dim(combinations(n, i, 1:n))[1]
6 x <- rbind(x, cbind(combinations(n, i, 1:n), matrix(0,nc, n - i)))
7 }
8x
9}R2current，每种组合的条件下与OTU计算cca或rda，得到R2，并输出到allR2的new.line。
代码语言：javascript复制1R2current=function(Y, current.comb, X,type)
2{require(vegan)
3 data <- data.frame(X[, current.comb])
4 if(type=="RDA")
5 gf <- RsquareAdj(rda(Y~., data = data))$r.squared
6 if(type=="CCA")
7 gf <- RsquareAdj(cca(Y~., data = data))$r.squared
8 gf
9}到这为止主函数前12行完成。14行是另一个函数： partition.rda,用到了hier.part进行Hierarchical Partitioning。
代码语言：javascript复制 1partition.rda=function (gfs, pcan, var.names = NULL)
2{require(hier.part)
3 if (pcan > 12)
4 stop("Number of variables must be < 13 for current implementation",call. = FALSE)
5 else if (pcan > 9)
6 warning("rdaenvpart produces a rounding error if number of variables >9nSee documentation.",
7 call. = FALSE)
8 {
9 n <- 2^pcan
10 theta <- gfs[1]
11 fin <- gfs[2:n]
12
13 len <- length(fin)
14 IJ <- vector("numeric", pcan * 2)
15 storage.mode(pcan) <- "integer"
16 storage.mode(len) <- "integer"
17 storage.mode(theta) <- "double"
18 storage.mode(fin) <- "double"
19 storage.mode(IJ) <- "double"
20 IJ <- .C("hierpart", pcan, len, theta, fin, IJ = IJ,
21 PACKAGE = "hier.part")$IJ
22 IJ <- array(IJ, dim = c(pcan, 2))
23 IJ <- data.frame(t(data.frame(t(IJ), row.names = c("I",
24 "J"))), row.names = var.names)
25 IJ.perc <- IJ * 100/sum(IJ)
26 I <- data.frame(I = IJ[, 1], row.names = var.names)
27 I.perc <- I * 100/sum(I)
28 IJ <- cbind(IJ, Total = IJ$I IJ$J)
29 list(gfs = gfs, IJ = IJ, I.perc = I.perc)
30 }
31}效果图，感觉画图这块还可以美化一下
另外还看到了两篇赖老师关于rda的文章： 关于RDA中每个环境因子解释率的说明 http://wap.sciencenet.cn/home.php?mod=space&uid=267448&do=blog&id=1186999
关于冗余分析（RDA）中环境因子共同解释部分出现负值的说明 http://wap.sciencenet.cn/home.php?mod=space&uid=267448&do=blog&id=1187530
END
