2

我使用素食社区生态包在 R 中执行了“基于距离的冗余分析”（dbRDA）。我想在 dbRDA 结果的排序图中显示（鱼）营养组对样本之间差异（营养级鱼类组合的丰度数据）的相对贡献。即在排序图上叠加箭头和营养级组名称，其中箭头线的长度表示对差异的相对贡献。据我了解，这应该可以通过vegan::scores()函数访问，或者与dbrda.model$CCA$v对象一起存储。

但是，当使用 时，物种分数是空的 (NA) dbrda()。我知道 dbrda 函数需要在函数中定义社区矩阵才能提供物种分数。我已将其定义为这样，但仍然无法产生物种分数。令我困惑的是，当我capscale()在 vegan 包中使用相同的物种群落和环境变量数据，并在各自的函数中制定相同的公式时，会产生物种分数。dbrda素食主义者能够产生物种分数吗？这些分数与由capscale（当使用相同的数据和公式时）？我提供了我的数据示例和使用的公式。（我对实际绘制一旦获得的物种分数很有信心 - 所以将代码限制为产生物种分数。）

#Community data matrix (comm.dat): site names = row names, trophic level = column names
>head(comm.dat[1:5,1:4])

            algae/invertebrates corallivore  generalist carnivore herbivore
h_m_r_3m_18                   1           0                    3         0              
h_m_r_3m_22                   6           4                    8        26                     
h_m_r_3s_19                   0           0                    4         0                      
h_m_r_3s_21                   3           0                    7         0                      
l_pm_r_2d_7                   1           0                    5         0   

> str(comm.dat)

num [1:47, 1:8] 1 6 0 3 1 8 11 2 6 9 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:47] "h_m_r_3m_18" "h_m_r_3m_22" "h_m_r_3s_19" "h_m_r_3s_21" ...
..$ : chr [1:8] "algae/invertebrates" "corallivore" "generalist carnivore" "herbivore" ...

# environmental data (env.dat): Already standardised env data.
>head(env.dat[1:5,1:3])

      depth   water.level  time.in
-0.06017376   1.3044232   -1.7184415
-0.67902862   1.3044232   -1.7907181
-0.99619174   1.3044232   -1.7569890
-1.06581291   1.3044232   -1.7762628
 2.39203863  -0.9214933    0.1703884

# Dissimilarity distance: Modified Gower (AltGower) with logbase 10 transformation of the community data matrix

> dis.comm.mGow <- vegan::vegdist(x = decostand(comm.dat, "log", logbase = 10), method = "altGower")

# Distance based RDA model: Trophic level data logbase transformed modified Gower distance, constrained against the interaction of dept and water level (tide), and the interaction of depth and time of day sampled`

> m.dbrda <- dbrda(formula = dis.comm.mGow ~ depth*water.level + depth*time.in, data = env.dat, scaling = 2, add = "lingoes", comm = decostand(comm.dat, "log", logbase = 10), dfun = "altGower")

# Check species scores: PROBLEM: No species level scores available

> m.dbrda$CCA$v

      dbRDA1 dbRDA2 dbRDA3 dbRDA4 dbRDA5
[1,]     NA     NA     NA     NA     NA

# OR pull species scores using scores(): Also does not show species scores...

>scrs <- scores(m.dbrda,display="species"); scrs

       dbRDA1 dbRDA2
spe1     NA     NA
attr(,"const")
[1] 6.829551

# when replacing dbrda with capscale, species scores are produced, e.g.


> m.cap <- capscale(formula = dis.comm.mGow ~ depth*water.level + depth*time.in, data = env.dat, scaling = 2, add = "lingoes", comm = decostand(comm.dat, "log", logbase = 10), dfun = "altGower")

> m.cap$CCA$v[1:5,1:3]

                         CAP1        CAP2        CAP3
algae/invertebrates    0.2044097 -0.04598088 -0.37200097
corallivore            0.3832594  0.06416886 -0.27963122
generalist carnivore   0.1357668 -0.08566365 -0.06789812
herbivore              0.5745226 -0.45647341  0.73085661
invertebrate carnivore 0.1987651  0.68036211 -0.19174283
4

2 回答 2

3

该dbrda()函数确实没有获得物种分数。但是，如果我们在vegan中有以下函数，您可以添加这些分数：

#' Add Species Scores to Ordination Results
#'
#' @param object Ordination object
#' @param comm Community data
#'
#' @return Function returns the ordination object amended with species
#'     scores.
#'
#' @rdname specscores
#' @export
`specscores` <-
    function(object, comm)
{
    UseMethod("specscores")
}
#' importFrom vegan decostand
#'
#' @rdname specscores
#' @export
`specscores.dbrda` <-
    function(object, comm)
{
    comm <- scale(comm, center = TRUE, scale = FALSE)
    if (!is.null(object$pCCA) && object$pCCA$rank > 0) {
        comm <- qr.resid(object$pCCA$QR, comm)
    }
    if (!is.null(object$CCA) && object$CCA$rank > 0) {
        v <- crossprod(comm, object$CCA$u)
        v <- decostand(v, "normalize", MARGIN = 2)
        object$CCA$v <- v
        comm <- qr.resid(object$CCA$QR, comm)
    }
    if (!is.null(object$CA) && object$CA$rank > 0) {
        v <- crossprod(comm, object$CA$u)
        v <- decostand(v, "normalize", MARGIN = 2)
        object$CA$v <- v
    }
    object
}

我不知道我们是否会有这个功能（使用其他方法），但如果你在会话中获取它，你可以使用它。

于 2017-10-03T12:04:27.670 回答
1

这看起来像是文档中的错误或缺少功能，我将在 Github 上与 Jari 一起提出这个问题。

但是，请注意，dbrda() 它没有参数comm，因此传递它实际上并没有做任何事情。comm仅用于capscale。

如果这是一个文档错误，那么它的出现是因为我们在同一页面上同时记录了两者，但capscale()首先存在并且稍后才到达，并且更新了文档以适应后者，但这可能没有完美地完成，因此混淆了什么是或没有记录在案的。dbrda()capscale()dbrda()dbrda()

就目前而言，背后的代码dbrda()没有尝试计算物种分数。如果我正确地遵循代码，它可以。因此，这可能只是缺少功能。

目前，NAs 是实现所有这些模型的底层ordConstrained函数的预期结果，如果分析是在相异矩阵上时我们返回的结果，就这个函数而言，它是。

于 2017-10-03T02:45:13.943 回答