尽管提供了社区矩阵,但纯素 dbrda 物种得分为空

vegan dbrda species scores are empty despite community matrix provided

我已经使用素食社区生态包在 R 中执行了 "Distance based redundancy analysis" (dbRDA)。我想在 dbRDA 结果的排序图中显示(鱼类)营养组对样本(营养级鱼类组合的丰度数据)之间差异的相对贡献。 IE。将箭头和营养级别组名称叠加到排序图上,其中箭头线的长度表示对差异的相对贡献。据我所知,这应该可以通过 vegan::scores() 函数访问,或者与 dbrda.model$CCA$v 对象一起存储。

但是,使用 dbrda() 时,物种分数为空 (NA)。我知道 dbrda 函数需要在函数内定义社区矩阵才能提供物种分数。我已经这样定义了它,但仍然无法产生物种分数。令我困惑的是,当我在 vegan 包中使用 capscale() 时,使用相同的物种群落和环境变量数据,并在各自的函数中制定相同的公式,就会产生物种分数。 dbrda in vegan 是否能够产生物种分数?这些分数与 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

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

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

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

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

目前,NAs 是实现所有这些模型的基础 ordConstrained 函数的预期结果,如果分析是基于相异矩阵,就此函数而言,它是。

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
}

我不知道我们是否会拥有此功能(使用其他方法),但如果您在会话中获取它,则可以使用它。