R/coef.R

In florale/multilevelcoda: Estimate Bayesian Multilevel Models for Compositional Data

#' Calculate "coefficients" based on substitutions for each compositional part
#'
#' @param object An object of class \code{brmcoda}.
#' @param level A character string specifying the level of the coefficients to be calculated.
#'   Either \dQuote{between} or \dQuote{within}.
#' @param h A numeric value specifying the step size for the substitution.
#' @return A data table of results.
#' @importFrom data.table as.data.table copy
#' @importFrom stats fitted model.frame
#' @importFrom testthat expect_equal
#' @importFrom brms posterior_summary
#' @importFrom nlme fixef
#' @export
#' @examples
#' \donttest{
#' if(requireNamespace("cmdstanr")){
#' sbp2 <- sbp
#' sbp2[1, ] <- c(-1, 1, -1, -1, -1)
#' sbp2[2, ] <- c( 1, 0, -1, -1, -1)
#' sbp2[3, ] <- c( 0, 0,  1, -1, -1)
#' sbp2[4, ] <- c( 0, 0,  0,  1, -1)
#' 
#' cilr <- complr(data = mcompd, sbp = sbp2, 
#'   parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
#'   idvar = "ID")
#' 
#' m1 <- brmcoda(complr = cilr,
#'               formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 +
#'                                  wilr1 + wilr2 + wilr3 + wilr4 + (1 | ID),
#'               chain = 4, iter = 1000, cores = 4L,
#'               backend = "cmdstanr")
#' substition_coef(m1, level = "between", h = 10)
#' substition_coef(m1, level = "within", h = 10)
#' rm(sbp2, cilr, m1) ## cleanup
#' }
#' }
substition_coef <- function(object, level = c("between", "within"), h = 10) {
  level <- match.arg(level)
  expect_s3_class(object, "brmcoda")

  if (object$model$family$family == "gaussian" && object$model$family$link == "identity") {
    linear <- TRUE
  } else {
    linear <- FALSE
  }

  parts <- object$complr$parts
  x <- object$complr$data[, ..parts]

  if (isFALSE(linear)) {
    y0 <- fitted(object$model,
      newdata = model.frame(object),
      re_formula = NA,
      summary = FALSE
    )
  }

  out <- vector("list", length(parts))

  for (k in seq_along(parts)) {
    w <- (-x) / rowSums(x[, -..k])
    w <- as.data.table(w)
    w[, (parts[k]) := +1]

    x2 <- x + (h * w)

    expect_equal(rowSums(x2), rowSums(x), tolerance = 1e-3)

    rm(w) ## cleanup

    d2 <- copy(object$complr$data)
    d2[, (parts) := x2]

    rm(x2) ## cleanup

    cilr2 <- complr(
      data = d2,
      sbp = object$complr$sbp,
      parts = parts,
      idvar = object$complr$idvar
    )

    if  (isFALSE(linear)) {
    switch(level,
      within = {
        bilr2 <- object$complr$between_logratio
        wilr2 <- cilr2$logratio - object$complr$between_logratio
        names(wilr2) <- names(object$complr$within_logratio)
      },
      between = {
        # expect_equal(cilr$within_logratio, cilr2$within_logratio, tolerance = 1e-3)
        bilr2 <- cilr2$between_logratio
        wilr2 <- object$complr$within_logratio
      }
    )

    rm(cilr2) ## cleanup

    d2 <- cbind(d2, bilr2, wilr2)

    rm(bilr2, wilr2)

    y2 <- fitted(object$model,
      newdata = d2,
      re_formula = NA,
      summary = FALSE
    )

    rm(d2) ## cleanup

    out[[k]] <- rowMeans((y2 - y0) )

    rm(y2)
    } else if (isTRUE(linear)) {
      switch(level,
        within = {
        wilr2 <- cilr2$logratio - object$complr$between_logratio
        out[[k]] <- fixef(object$model, summary = FALSE)[, colnames(object$complr$within_logratio)] %*% 
          colMeans(wilr2 - object$complr$within_logratio)
      },
      between = {
        out[[k]] <- fixef(object$model, summary = FALSE)[, colnames(object$complr$between_logratio)] %*% 
          colMeans(cilr2$between_logratio - object$complr$between_logratio)
      }
    )
    }
  }

  if (isTRUE(linear)) rm(x) else rm(x, y0) ## cleanup

  finalout <- cbind(
    Component = parts,
    as.data.table(do.call(rbind, lapply(out, posterior_summary)))
  )

  return(finalout)
}