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 Values defining 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{diagonal = -0.5})
#' @param connectance Numeric frequency of inter-species interactions.
#' i.e. proportion of non-zero off-diagonal terms.
#' Should be in the interval 0 <= connectance <= 1.
#' (default: \code{connectance = 0.2})
#' @param scale_off_diagonal Numeric: scale of the off-diagonal elements
#' compared to the diagonal.
#' (default: \code{scale_off_diagonal = 0.1})
#' @param mutualism Numeric: relative proportion of interactions terms
#' consistent with mutualism (positive <-> positive)
#' (default: \code{mutualism = 1})
#' @param commensalism Numeric: relative proportion of interactions terms
#' consistent with commensalism (positive <-> neutral)
#' (default: \code{commensalism = 1})
#' @param parasitism Numeric: relative proportion of interactions terms
#' consistent with parasitism (positive <-> negative)
#' (default: \code{parasitism = 1})
#' @param amensalism Numeric: relative proportion of interactions terms
#' consistent with amensalism (neutral <-> negative)
#' (default: \code{amensalism = 1})
#' @param competition Numeric: relative proportion of interactions terms
#' consistent with competition (negative <-> negative)
#' (default: \code{competition = 1})
#' @param interactions Numeric: 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{interactions = NULL})
#' @param symmetric Logical: 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{symmetric=FALSE})
#' @param list_A 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{list_A = 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)
}