尽管提供了社区矩阵,但纯素 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()
背后的代码 没有尝试 计算物种分数;如果我正确地遵循代码,它可以。因此,这可能只是缺少功能。
目前,NA
s 是实现所有这些模型的基础 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
}
我不知道我们是否会拥有此功能(使用其他方法),但如果您在会话中获取它,则可以使用它。
我已经使用素食社区生态包在 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()
背后的代码 没有尝试 计算物种分数;如果我正确地遵循代码,它可以。因此,这可能只是缺少功能。
目前,NA
s 是实现所有这些模型的基础 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
}
我不知道我们是否会拥有此功能(使用其他方法),但如果您在会话中获取它,则可以使用它。