R/randomA.R

In microbiome/miaSim: Microbiome Data Simulation

#' Generate random interaction matrix for GLV model
#'
#' Generates a random interaction matrix for Generalized Lotka-Volterra (GLV) model.
#'
#' @template man_spe
#' @param diagonal \code{Numeric vector}. Defines the strength of 
#' self-interactions. Input can be a number (will be applied to
#' all species) or a vector of length n_species. Positive 
#' self-interaction values lead to exponential growth.
#' (Default: \code{-0.5})
#' @param connectance \code{Numeric scalar}. Specifies the frequency 
#' of inter-species interactions. i.e. proportion of non-zero 
#' off-diagonal terms. Between \code{0} and \code{1}.
#' (Default: \code{0.2})
#' @param scale_off_diagonal \code{Numeric scalar}. Indicates the 
#' scale of the off-diagonal elements compared to the diagonal.
#' (Default: \code{0.1})
#' @param mutualism \code{Numerical scalar}. Specifies the relative 
#' proportion of interactions terms consistent with mutualism 
#' (positive <-> positive). (Default: \code{1})
#' @param commensalism \code{Numeric scalar}. Indicates the relative 
#' proportion of interactions terms consistent with commensalism 
#' (positive <-> neutral). (Default: \code{1})
#' @param parasitism \code{Numeric scalar}. Indicates the relative 
#' proportion of interactions terms consistent with parasitism 
#' (positive <-> negative). (Default: \code{1})
#' @param amensalism \code{Numeric scalar}. Indicates the relative 
#' proportion of interactions terms consistent with amensalism 
#' (neutral <-> negative). (Default: \code{1})
#' @param competition \code{Numeric scalar}. Indicates the relative 
#' proportion of interactions terms consistent with competition 
#' (negative <-> negative). (Default: \code{1})
#' @param interactions \code{Numeric scalar}. Indicates the values 
#' of the n_species^2 pairwise interaction strengths. Diagonal terms 
#' will be replaced by the 'diagonal' parameter. If NULL, interactions 
#' are drawn from `runif(n_species^2, min=0, max=abs(diagonal))`.
#' Negative values are first converted to positive then the signs are defined by
#' the relative weights of the biological interactions (i.e. mutualism,
#' commensalism, parasitism, amensalism, competition).
#' (Default: \code{NULL})
#' @param symmetric \code{Logical scalar}. whether the strength of 
#' mutualistic and competitive interactions are symmetric. This is 
#' implemented by overwrite a half of the matrix, so the proportions 
#' of interactions might deviate from expectations.
#' (Default: \code{FALSE})
#' @param list_A \code{List}. A list of matrices generated by randomA. Used to support
#' different groups of interactions. If NULL (by default), no group is
#' considered. Otherwise the given list of matrices will overwrite values around
#' the diagonal. (Default: \code{NULL})
#' @examples
#'
#' dense_A <- randomA(
#'     n_species = 10,
#'     scale_off_diagonal = 1,
#'     diagonal = -1.0,
#'     connectance = 0.9
#' )
#'
#' sparse_A <- randomA(
#'     n_species = 10,
#'     diagonal = -1.0,
#'     connectance = 0.09
#' )
#'
#' user_interactions <- rbeta(n = 10^2, .5, .5)
#' user_A <- randomA(n_species = 10, interactions = user_interactions)
#'
#' competitive_A <- randomA(
#'     n_species = 10,
#'     mutualism = 0,
#'     commensalism = 0,
#'     parasitism = 0,
#'     amensalism = 0,
#'     competition = 1,
#'     connectance = 1,
#'     scale_off_diagonal = 1
#' )
#'
#' parasitism_A <- randomA(
#'     n_species = 10,
#'     mutualism = 0,
#'     commensalism = 0,
#'     parasitism = 1,
#'     amensalism = 0,
#'     competition = 0,
#'     connectance = 1,
#'     scale_off_diagonal = 1,
#'     symmetric = TRUE
#' )
#'
#' list_A <- list(dense_A, sparse_A, competitive_A, parasitism_A)
#' groupA <- randomA(n_species = 40, list_A = list_A)
#'
#' @return
#' \code{randomA} returns a matrix A with dimensions (n_species x n_species)
#'
#' @export
randomA <- function(n_species,
    names_species = NULL,
    diagonal = -0.5,
    connectance = 0.2,
    scale_off_diagonal = 0.1,
    mutualism = 1,
    commensalism = 1,
    parasitism = 1,
    amensalism = 1,
    competition = 1,
    interactions = NULL,
    symmetric = FALSE,
    list_A = NULL) {
    # input check
    if (!.isPosInt(n_species)) {
        stop("n_species must be integer.")
    }
    if (!all(vapply(
        list(
            diagonal, connectance, scale_off_diagonal,
            mutualism, commensalism, parasitism, amensalism,
            competition
        ),
        is.numeric, logical(1)
    ))) {
        stop("diagonal, connectance, scale_off_diagonal, mutualism,
                     commensalism, parasitism, amensalism and competition must
                     be numeric.")
    }

    if (connectance > 1 || connectance < 0) {
        stop("'connectance' should be in range: 0 <= connectance <= 1")
    }
    # set the default values
    if (is.null(names_species)) {
        if (!is.null(list_A)) {
            names_species <- unlist(lapply(list_A, rownames))
            # if same names found in different groups, overwrite names.
            if (length(names_species) != length(unique(names_species))) {
                list_Alengths <- unlist(lapply(list_A, nrow))
                sn <- c()
                for (i in seq_len(length(list_Alengths))) {
                    sn <- c(sn, rep(i, times = list_Alengths[i]))
                }
                names_species <- paste0("g", sn, "_", names_species)
            }
        } else {
            names_species <- paste0("sp", seq_len(n_species))
        }
    }
    if (is.null(interactions)) {
        A <- runif(n_species^2, min = 0, max = abs(diagonal))
    } else {
        A <- interactions
    }

    interaction_weights <- c(
        mutualism,
        commensalism,
        parasitism,
        amensalism,
        competition
    )

    A <- matrix(A, nrow = n_species, ncol = n_species)
    I <- .getInteractions(n_species, interaction_weights, connectance)

    if (symmetric) {
        A[lower.tri(A)] <- t(A)[lower.tri(A)]
        I[lower.tri(I)] <- t(I)[lower.tri(I)]
    }

    A <- I * abs(A) * (scale_off_diagonal * min(abs(diagonal)))
    diag(A) <- diagonal
    colnames(A) <- rownames(A) <- names_species

    if (!is.null(list_A)) {
        if (n_species != sum(unlist(lapply(list_A, nrow)))) {
            stop("'n_species' should equals to the number of species in the
                 given list")
        }
        counter <- 0
        for (repA in list_A) {
            A[(counter + 1):(counter + nrow(repA)), (counter + 1):(counter + nrow(repA))] <- repA
            counter <- counter + nrow(repA)
        }
    }

    return(A)
}