Concurrent ordination in ggplot

[mike-kratz]

Hi there!

Thank you for creating this awesome package! I generated a concurrent ordination model for microbial community data with four predictor variables. The ordiplot.gllvm function is handy, but I would like to display my results in ggplot2 since it allows me to quickly change visualization parameters. Unfortunately some of the scaling seems off when I ordinate using extracted site LV scores and predictor LV scores. I am not sure if I need to multiply the scores by a constant in order for the ggplot2 to match the ordiplot.gllvm plot. Here is the code I used to generate my ggplot:

############################Constrained Ordination Dataframes############################

Extract Predictor scores

Env.con = as.data.frame(fit.gllvm.ASV.con$params$LvXcoef)

Extract site scores

LVs = as.data.frame(getLV(fit.gllvm.ASV.con, type = "conditional"))

Bind LV site scores to metadata

LVs = cbind(LVs, sim.meta)

####################################GGPLOT######################################## ggplot(LVs, aes(x = CLV1, y = CLV2))+ geom_point(aes(color = EC, #Colored by Electrical conductivity, proxy for salinity shape = Date), size = 8, data = LVs)+ geom_text( #Create labels for the predictor arrows data=Env.con, aes(x= CLV1, y= CLV2, label=row.names(Env.con), angle = (180/pi) atan(Env.con[ ,2]/Env.con[ ,1]), #Code to angle the labels parallel to predictor arrows hjust = (1 - 2 sign(Env.con[ ,1])) / 2.75), size= 10, colour="black", fontface = "bold")+ geom_hline(yintercept=0, #Creates dotted lines to pass through plot origin linetype=3, size=1) + geom_vline(xintercept=0, linetype=3, linewidth=1)+ geom_segment(data = Env.con, #Create arrows for predictors aes(x = 0, y = 0, xend = CLV1, yend = CLV2), arrow = arrow(angle=22.5, length = unit(0.35,"cm"), type = "closed"), linetype=1, size=3, colour = "black")+ mytheme+ scale_color_viridis_c()

Plots from both formats

'ordiplot.gllvm'

ggplot2

Thank you,

Mike

[BertvanderVeen]

Thanks Mike.

The ordiplot function does quite a lot of scaling, and some rotation. First a rotation matrix is determined based on the latent variables, and extra scaling is added based on both the latent variables and the species loadings.

The arrows are additionally scaled by the scale of the predictor variables.

There is no functionality in gllvm currently that allows to extract the rotation & scaling matrix, so you will have to mine the code to get it (I can help, but I do not have time this week).

[mike-kratz]

@BertvanderVeen Ok that makes sense, thank you for your response!

[RodolfoPelinson]

Hello @mike-kratz and @BertvanderVeen

By mining the code, I made a function to get scaled and rotated latent variables and species loadings. I am not totally sure that it is correct, but it seems to produce similar plots as the ordiplot function. I think a future implementation of a separate function that returns scaled and rotated lvs and species loadings would be really useful. I often like to customize my plots :D

Here is the function, maybe it helps somehow.

get_scaled_lvs <- function(object, alpha = 0.5){

  choose.lvs <- getLV(object, type = "scaled")
  choose.lv.coefs <- object$params$theta

  #Compute the singular-value decomposition of a rectangular matrix for ratation
  do_svd <- svd(object$lvs)

  svd_rotmat_sites <- do_svd$v
  svd_rotmat_species <- do_svd$v

  p <- NCOL(object$y)
  n <- NROW(object$y)
  num.lv <- object$num.lv

  idx <- matrix(TRUE, ncol = num.lv,  nrow = p)

  bothnorms <- vector("numeric", ncol(choose.lv.coefs))
  for (i in 1:ncol(choose.lv.coefs)) {
    bothnorms[i] <- sqrt(sum(choose.lvs[, i]^2)) * sqrt(sum(choose.lv.coefs[idx[, i], i]^2))
  }

  #Scale sites
  scaled_cw_sites <- t(t(choose.lvs)/sqrt(colSums(choose.lvs^2)) *       (bothnorms^alpha))

  #Scale species
  scaled_cw_species <- choose.lv.coefs
  for (i in 1:ncol(scaled_cw_species)) {
    scaled_cw_species[, i] <- choose.lv.coefs[, i]/sqrt(sum(choose.lv.coefs[idx[,
                                                                                i], i]^2)) * (bothnorms[i]^(1 - alpha))
  }

  choose.lvs <- scaled_cw_sites %*% svd_rotmat_sites
  choose.lv.coefs <- scaled_cw_species %*% svd_rotmat_species

  result <- list(sites = choose.lvs, species = choose.lv.coefs)

  return(result)
}`

[BertvanderVeen]

Nice. Just a small note: if you rotate the latent variables you will also need to rotate the predictors effects. Rotating and scaling the site scores and species loadings will not be sufficient for constrained and concurrent ordinations.

[mike-kratz]

Hi @RodolfoPelinson! Thank you for sharing your code, seems like it must've taken a lot of work. @BertvanderVeen Ok I tried implementing rotated predictors based off of @RodolfoPelinson's code but this is clearly out of my wheelhouse; I should probably just wait until you are able to edit it next week. Thank you both for your help!

[BertvanderVeen]

Gettings things -exactly the same- as in ordiplot might be difficult, but probably also unnecessary. As long as the rotation is performed properly and to all entities included in the plot, things should be OK. You can try the code below (at your own peril :) ), I have put snippets of code for rotation and scaling from ordiplot together in a function for you. Note that this does not give SEs that are rotated and scaled with the same procedure (i.e., do not combine SEs from the model with what comes out of the function below).

Disclaimer: I have not tested this, so no guarantees that it works in all cases (especially not with the quadratic model), and make sure to compare the results again with ordiplot.

rotate <- function(object, alpha = 0.5, type = NULL, which.lvs = 1:2) {
    p <- ncol(object$y)
    num.lv <- object$num.lv
    num.lv.c <- object$num.lv.c
    num.RR <- object$num.RR
    quadratic <- object$quadratic

    if ((num.lv.c + num.RR) == 0 & is.null(type)) {
      type <- "residual"
    } else if (is.null(type)) {
      if (num.lv.c == 0 & num.RR > 0) {
        type <- "marginal"
      } else if (num.lv.c > 0) {
        type <- "conditional"
      }
    }

    if (type == "residual" | num.lv > 0) {
      # First create an index for the columns with unconstrained LVs
      sigma.lv <- NULL
      if (num.lv.c > 0 & type == "residual") {
        # to have similar scaling to unconstrained LVs
        sigma.lv <- object$params$sigma.lv[1:num.lv.c]
      } else if (num.lv.c > 0) {
        #here the LVs get scaled with sigma.lv not theta
        sigma.lv <- c(sigma.lv, rep(1, num.lv.c))
      }
      # constrained LVs never have a scale parameter so always get 1
      sigma.lv <- c(sigma.lv, rep(1, num.RR))

      # unconstrained LVs always get their species parameters scaled, never LVs
      if (num.lv > 0) {
        sigma.lv <-
          c(sigma.lv, object$params$sigma.lv[(num.lv.c + 1):(num.lv.c + num.lv)])
      }
      # scaling for quadratic coefficients
      if (quadratic != FALSE) {
        sigma.lv <- c(sigma.lv, sigma.lv ^ 2)
      }

      # Do the scaling
      sigma.lv <- diag(sigma.lv, length(sigma.lv))
      object$params$theta <- object$params$theta %*% sigma.lv

      #remove unconstrained LVs species loadings if type=="marginal"
      if (type == "marginal" & num.lv > 0) {
        if (quadratic != FALSE)
          object$params$theta <- object$params$theta[, -c((num.lv.c + num.RR + 1):(num.lv.c + num.RR + num.lv) +  num.RR + num.lv.c + num.lv), drop = F]
        object$params$theta <- object$params$theta[, -c((num.lv.c + num.RR + 1):(num.lv.c + num.RR + num.lv)), drop = F]
      }
      #remove constrined LVs species loadings if type=="residual"
      if (type == "residual" & num.RR > 0) {
        if (quadratic != FALSE)
          object$params$theta <-
            object$params$theta[, -c((num.lv.c + 1):(num.lv.c + num.RR) + num.RR + num.lv.c + num.lv), drop = F]
        object$params$theta <-
          object$params$theta[, -c((num.lv.c + 1):(num.lv.c + num.RR)), drop = F]
      }
    }

    lv <- getLV(object, type = type)

    do_svd <- svd(lv)

    if (type != "residual") {
      do_svd$v <- svd(getLV(object, type = type))$v
    }

    # do_svd <- svd(lv)
    # do_svd <- svd(object$lvs)
    svd_rotmat_sites <- do_svd$v
    svd_rotmat_species <- do_svd$v

    choose.lvs <- lv

    if (quadratic == FALSE) {
      choose.lv.coefs <-  object$params$theta
    } else{
      choose.lv.coefs <- optima(object, sd.errors = F)
    }

    bothnorms <- vector("numeric", ncol(choose.lv.coefs))
    for (i in 1:ncol(choose.lv.coefs)) {
      bothnorms[i] <- sqrt(sum(choose.lvs[, i] ^ 2)) * sqrt(sum(choose.lv.coefs[, i] ^ 2))
    }

    scaled_cw_sites <-  t(t(choose.lvs) / sqrt(colSums(choose.lvs ^ 2)) * (bothnorms ^ alpha))
    # scaled_cw_species <- t(t(choose.lv.coefs) / sqrt(colSums(choose.lv.coefs^2)) * (bothnorms^(1-alpha)))
    scaled_cw_species <- choose.lv.coefs
    for (i in 1:ncol(scaled_cw_species)) {
      scaled_cw_species[, i] <-
        choose.lv.coefs[, i] / sqrt(sum(choose.lv.coefs[, i] ^ 2)) * (bothnorms[i] ^ (1 - alpha))
    }

    choose.lvs <- scaled_cw_sites %*% svd_rotmat_sites
    choose.lv.coefs <- scaled_cw_species %*% svd_rotmat_species

    LVcoef  = NULL
    # Continue for env preds
    if (num.lv == 0 & (num.lv.c + num.RR) > 0 & type != "residual" | (num.lv.c + num.RR) > 0 &num.lv > 0 & type == "marginal") {
      LVcoef <- (object$params$LvXcoef %*% svd_rotmat_sites)[, which.lvs]
      LVcoef <- LVcoef / apply(object$lv.X, 2, sd)
    }

    result <- list(sites = choose.lvs, species = choose.lv.coefs, env = LVcoef)

    return(result)
  }

[mike-kratz]

The code worked perfectly, thank you so much @BertvanderVeen!!! The only things that differ between the ordiplot and ggplot are the plot origin (which I can change manually in the ggplot code if needed) and the overall length of the arrows (which likely don't matter since they convey the same message but with different scaling). Do you think you could implement this function into the package so other people can readily use ggplot2 for constrained/concurrent ordinations?

ggplot with @BertvanderVeen's function

original ordiplot