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
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.