【发布时间】：2015-07-07 06:09:57
【问题描述】：
我正在使用“功能指数用户指南”中的脚本来计算功能多样性指数。但是，该脚本充满了错误和问题。我现在设法让脚本运行，但它产生了不正确的答案。
它使用如下三个矩阵：
丰富：
S1 S2 S3 S4 S5 S6 S7 S8 S9
Palm 6 3 0 1 0 16 0 2 3
Forest 2 0 2 1 2 1 3 0 2
分钟：
min
S1 25.3038
S2 19.5750
S3 60.5880
S4 16.2864
S5 46.1040
S6 10.9056
S7 8.7570
S8 2.1289
S9 4.1730
最大值：
max
S1 44.7344
S2 22.6966
S3 75.1817
S4 17.8176
S5 50.7472
S6 33.3660
S7 14.3341
S8 3.3947
S9 10.2510
我认为它应该根据Palm 和Forest 两个站点产生两个不同的数字。相反，它会产生两个完全相同的数字。我查看了它，发现问题接缝来自这些代码行：
Siteminsp <- abundance*minsp
Siteminsp <- minsp[apply(Siteminsp,1,function (x) sum(x,na.rm=T))!=0, ,drop = FALSE]
Sitemaxsp <- abundance*maxsp
Sitemaxsp <- maxsp[apply(Sitemaxsp,1,function (x) sum(x,na.rm=T))!=0, ,drop = FALSE]
而不是abundance*minsp 生成一个新矩阵，其中每个物种的最小值乘以每个站点该物种的丰度。该函数是将第一个物种丰度乘以第一个物种最小值，然后将另一个站点的第一个物种丰度乘以第二个最小值，然后将第二个物种丰度乘以第三个值，依此类推。
我是否错过了这里应该发生的事情？
有人可以解释一下出了什么问题，因为我完全没有想法。如果有帮助，这是整个脚本：
rownames(Abundance1) <- Abundance1[,1]
Abundance <- Abundance1[,-1]
Abundance <- Abundance[,order(colnames(Abundance))]
rownames(min1) <- min1[,1]
minsp <- min1[,-1, drop= FALSE]
minsp <- minsp[order(rownames(min1)), , drop = FALSE]
rownames(max1) <- max1[,1]
maxsp <- max1[,-1, drop=FALSE]
maxsp <- maxsp[order(rownames(max1)), , drop = FALSE]
globalFRI <- function(minsp,maxsp) {
deltaS <- list()
for (j in 1:ncol(minsp))
{
xx <- cbind(minsp[,j], maxsp[,j])
xx <- xx[apply(xx,1,function(z) sum(is.finite(z))==2),]
xx <- xx[order(xx[,1]),]
z <- c(0, nrow(xx))
i <- 1
b <- xx[1, 2]
while ( i < nrow(xx) )
{
if (b < xx[i+1,1]) z <- c(z, i)
b <- ifelse(b >= xx[i+1,2], b, xx[i+1,2])
i <- i+1
if (i==nrow(xx)) break
}
group <- factor(rep(1:(length(z)-1), diff(sort(z))))
deltaS[[j]] <- tapply(xx[,2], group, max) - tapply(xx[,1], group, min)
}
globalFRIs <- sapply(deltaS, sum)}
globalFRIAll <- globalFRI(minsp,maxsp)
Abundance2 <- Abundance
Abundance2[Abundance2 != 0] <- 1
abundance <- data.matrix(Abundance)
IndexFRIs <- function(Abundance,minsp,maxsp,globalFRI) {
Siteminsp <- abundance*minsp
Siteminsp <- minsp[apply(Siteminsp,1,function (x) sum(x,na.rm=T))!=0, ,drop = FALSE]
Sitemaxsp <- abundance*maxsp
Sitemaxsp <- maxsp[apply(Sitemaxsp,1,function (x) sum(x,na.rm=T))!=0, ,drop = FALSE]
deltaS <- list()
for (j in 1:ncol(Siteminsp))
{
xx <- cbind(Siteminsp[,j], Sitemaxsp[,j])
xx <- xx[apply(xx,1,function(z) sum(is.finite(z))==2),]
xx <- xx[order(xx[,1]),]
z <- c(0, nrow(xx))
i <- 1
b <- xx[1, 2]
while ( i < nrow(xx) )
{
if (b < xx[i+1,1]) z <- c(z, i)
b <- ifelse(b >= xx[i+1,2], b, xx[i+1,2])
i <- i+1
if (i==nrow(xx)) break
}
group <- factor(rep(1:(length(z)-1), diff(sort(z))))
deltaS[[j]] <- tapply(xx[,2], group, max) - tapply(xx[,1], group, min)
}
RI <- sapply(deltaS, sum)
FRI1 <- RI/globalFRI
FRIs <- mean(FRI1)
}
FRIs <- apply(Abundance2,1,IndexFRIs,minsp,maxsp,globalFRIAll)
输入（丰度）：
structure(c(6L, 2L, 3L, 0L, 0L, 2L, 1L, 1L, 0L, 2L, 16L, 1L,
0L, 3L, 2L, 0L, 3L, 2L), .Dim = c(2L, 9L), .Dimnames = list(c("Palm",
"Forest"), c("S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8",
"S9")))
输入（分钟）：
structure(list(min = c(25.3038, 19.575, 60.588, 16.2864, 46.104,
10.9056, 8.757, 2.1289, 4.173)), .Names = "min", row.names = c("S1",
"S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9"), class = "data.frame")
输入（最大值）：
structure(list(max = c(44.7344, 22.6966, 75.1817, 17.8176, 50.7472,
33.366, 14.3341, 3.3947, 10.251)), .Names = "max", row.names = c("S1",
"S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9"), class = "data.frame")
我认为丰度*minsp 的输出应该是这样的：
Palm Forest
S1 151.8228 50.6076
S2 58.725 0
S3 0 121.176
S4 16.2864 16.2864
S5 92.208 0
S6 174.489 10.9056
S7 0 26.271
S8 4.2578 0
S9 8.346 12.519
【问题讨论】：
-
R 中的矩阵-矩阵乘法是使用
%*%运算符执行的。 -
为了使该示例可重现，请发布
dput(abundance)、dput(minsp)和dput(maxsp)的输出。然后，发布您想要的输出。 -
@nicola Ive 添加了您要求的 dput 输出。