vegan dbrda species scores are empty despite community matrix provided

Question

I have performed "Distance based redundancy analysis" (dbRDA) in R using the vegan community ecology package. I would like to show the relative contribution of (fish) trophic groups to the dissimilarity between samples (abundance data of trophic level fish assemblages) in an ordination plot of the dbRDA results. I.e. Overlay arrows and trophic level group names onto the ordination plot, where the length of the arrow line indicates the relative contribution to dissimilarity. This should be accessible via the vegan::scores() function, or stored with the dbrda.model$CCA$v object, as I understand.

However, species scores are empty (NA) when using dbrda(). I understand that the dbrda function requires the community matrix to be defined within the function in order to provide species-scores. I have defined it as such, but still unable to produce the species scores. What puzzles me is that when I use capscale() in the vegan package, with the same species-community and environmental variable data, and formulated the same within the respective functions, species scores are produced. Is dbrda in vegan able to produce species scores? How are these scores different to that produced by capscale (when the same data and formula are used)? I provide an example of my data, and the formula used. (I am fairly confident in actually plotting the species-scores once obtained - so limit the code to producing the species-scores.)

#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

Answer 1

The dbrda() function really does not get the species scores. However, if we had the following function in vegan, you could add those scores:

#' 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
}

I don't know if we are going to have this function (with other methods), but you can use it if you source it in your session.

Answer 2

This looks like either a bug in documentation or missing feature and I'll raise this with Jari on Github.

However, note that dbrda() doesn't have argument comm and hence passing this doesn't actually do anything. comm is for capscale only.

If this is a documentation bug, then it arises because we document both capscale() and dbrda() on the same page, but capscale() existed first and dbrda() only arrived later and the documentation was updated to accommodate the latter, but this may not have have been done perfectly, hence confusion about what dbrda() is or is not documented to do.

As it stands, the code behind dbrda() makes no attempt to compute species scores; it could, if I follow the code correctly. As such, this may just be missing functionality.

For the moment, the NAs are the expected result from the underlying ordConstrained function that implements all these models, and is what we return if the analysis was on a dissimilarity matrix, which as far as this function is concerned, it was.