R/helpfunctions_summary.R

In JointAI: Joint Analysis and Imputation of Incomplete Data

# used in various output/summary functions (2020-06-11)
prep_MCMC <- function(object, start = NULL, end = NULL, thin = NULL,
                      subset = NULL, exclude_chains = NULL, warn = TRUE,
                      mess = TRUE, ...) {
  # general processing of the MCMC sample for various output functions;
  # creating a subset of the sample wrt iterations, chains and parameters
  # - object: an object of class JointAI
  # - start: first iteration to be used
  # - end: last iteration to be used
  # - thin: thinning to be applied
  # - subset: subset of parameters (columns of the mcmc object) to be used
  # - exclude_chains: optional vector of numbers, indexing MCMC chains to be
  #                   excluded from the output
  # - warn: logical, should warning messages be displayed?
  # - mess: logical, should messages be displayed?


  # set start, end and thin (or use values from MCMC sample)
  if (is.null(start)) {
    start <- start(object$MCMC)
  } else {
    start <- max(start, start(object$MCMC))
  }

  if (is.null(end)) {
    end <- end(object$MCMC)
  } else {
    end <- min(end, end(object$MCMC))
  }

  if (is.null(thin))
    thin <- coda::thin(object$MCMC)

  # obtain subset of parameters of the MCMC samples
  MCMC <- get_subset(object, subset, warn = warn, mess = mess)

  # exclude chains, if set by user
  chains <- seq_along(MCMC)
  if (!is.null(exclude_chains)) {
    chains <- chains[-exclude_chains]
  }

  # restrict MCMC sample to selected iterations and chains
  MCMC <- do.call(rbind,
                  window(MCMC[chains],
                         start = start,
                         end = end,
                         thin = thin))
  return(MCMC)
}



get_D_names <- function(params, varname, lvls) {
  nams <- nlapply(lvls, function(lvl) {
    params$coef[
      lvapply(params$outcome, function(x) varname %in% x) &
        grepl(paste0("^D[[:digit:]]*_[", varname, "_]*", lvl, "\\["),
              params$coef)
    ]
  })

  Filter(length, nams)
}


#' Extract the random effects variance covariance matrix
#' Returns the posterior mean of the variance-covariance matrix/matrices of
#' the random effects in a fitted JointAI object.
#'
#' @inheritParams sharedParams
#' @param outcome optional; vector of integers giving the indices of the
#'                outcomes for which the random effects variance-covariance
#'                matrix/matrices should be returned.
#' @export
rd_vcov <- function(object, outcome = NULL, start = NULL, end = NULL,
                    thin = NULL,
                    exclude_chains = NULL, mess = TRUE, warn = TRUE) {


  vars <- if (is.null(outcome)) {
    names(object$coef_list)
  } else {
    names(object$coef_list)[outcome]
  }

  rd_vcov_list <- nlapply(vars, function(varname) {
    get_Dmat(object, varname, start = start, end = end, thin = thin,
             exclude_chains = exclude_chains, mess = mess, warn = warn)
  })

  if (length(rd_vcov_list) == 1L) {
    rd_vcov_list[[1]]
  } else {
    rd_vcov_list
  }
}




# used in print.JointAI() (2020-06-10)
get_Dmat <- function(object, varname, start = NULL, end = NULL, thin = NULL,
                     exclude_chains = NULL, mess = TRUE, warn = TRUE,
                     lvls = "all") {
  # Return the posterior mean of the random effects variance matrices in
  # matrix form (one matrix per grouping level)
  # - object: object of class JointAI
  # - varname: name of the outcome of the sub-model
  # - lvls: vector of the levels for which the D matrix should be returned

  if (lvls == "all") {
    lvls <- object$Mlist$idvar
  }

    MCMC <- prep_MCMC(object, start = start, end = end, thin = thin,
                      subset = NULL, exclude_chains = exclude_chains,
                      warn = warn, mess = mess)

    Ds <- get_D_names(parameters(object), varname = varname, lvls = lvls)

    Dpos <- lapply(Ds, function(d) {
      t(sapply(strsplit(gsub("[[:print:]]+\\[|\\]", "", d), ","),
               as.numeric))
    })

    Dmat <- nlapply(names(Dpos), function(lvl) {
      nam <- get_rdvcov_names(object, varname, lvl)

      m <- matrix(nrow = max(Dpos[[lvl]][, 1]),
                  ncol = max(Dpos[[lvl]][, 2]),
                  dimnames = list(nam$nam, nam$nam))
      structure(m,
                class = "Dmat",
                structure = nam
      )
    })


    for (lvl in names(Dmat)) {
      for (k in seq_along(Ds[[lvl]])) {
        Dmat[[lvl]][Dpos[[lvl]][k, 1],
                    Dpos[[lvl]][k, 2]] <- mean(MCMC[, Ds[[lvl]][k]])
      }
      Dmat[[lvl]][is.na(Dmat[[lvl]])] <- t(Dmat[[lvl]])[is.na(Dmat[[lvl]])]
    }

  Dmat
}


#' @rdname summary.JointAI
#' @param scientific A penalty to be applied when deciding to print numeric
#'                   values in fixed or exponential notation, by default the
#'                   value obtained from `getOption("scipen")`
#' @export
print.Dmat <- function(x, digits = getOption("digits"),
                       scientific = getOption("scipen"), ...) {

  r <- rbind(c(rep("", 2), attr(x, "structure")$variable),
             c(rep("", 2), colnames(x)),
             cbind(attr(x, "structure")$variable,
                   rownames(x),
                   unname(format(x, digits = digits,
                                 scientific = scientific, ...))
             )
  )

  cat(format_Dmat(r), sep = c(rep(" ",  ncol(r) - 1), "\n"))
}


format_Dmat <- function(r) {
  spaces <- sapply(max(nchar(r)) - nchar(r), function(k) {
    if (k > 0L) {
      paste0(rep(" ", k), collapse = "")
    } else {
      ""
    }
  })
  matrix(nrow = nrow(r), ncol = ncol(r), data = paste0(spaces, r)  )

}

get_rdvcov_names <- function(object, varname, lvl) {

  pos <- nlapply(attr(object$info_list[[varname]]$rd_vcov[[lvl]], "ranef_index"),
                function(nr) eval(parse(text = nr)))

  if (length(pos) > 0L) {
    nam <- nlapply(names(pos), function(v) {
      attr(object$info_list[[v]]$hc_list$hcvars[[lvl]], "z_names")
    })
    melt_list(
      Map(function(pos, nam) {
        cbind(pos, nam)
      }, pos = pos, nam = nam), varname = "variable"
    )

  } else {
    nam <- attr(object$info_list[[varname]]$hc_list$hcvars[[lvl]], "z_names")
    if (!is.null(nam)) {
      data.frame(pos = seq_along(nam),
                 nam = nam,
                 variable = varname)
    }
  }
}

# used in print.summary.JointAI(), print.JointAI(), list_models() (2020-06-18)
print_type <- function(type, family = NULL, upper = FALSE) {
  # collection of model titles to be printed at the start of the summary of
  # each sub-model
  # - type: model type
  # - family: family in case of a (extended) exponential family model

  a <- switch(type,
              # lm = "Linear model",
              glm = switch(family,
                           gaussian = 'linear model',
                           binomial = 'binomial model',
                           Gamma = 'Gamma model',
                           poisson = 'poisson model',
                           lognorm = 'log-normal model',
                           beta = 'beta model'
              ),
              # lme = "Linear mixed model",
              glmm = switch(family,
                            gaussian = 'linear mixed model',
                            binomial = 'binomial mixed model',
                            Gamma = 'Gamma mixed model',
                            poisson = 'poisson mixed model',
                            lognorm = 'log-normal mixed model',
                            beta = 'beta mixed model'
                            ),
              # glme = 'Generalized linear mixed model',
              coxph = 'proportional hazards model',
              survreg = 'weibull survival model',
              clm = 'cumulative logit model',
              clmm = 'cumulative logit mixed model',
              mlogit = "multinomial logit model",
              mlogitmm = "multinomial logit mixed model",
              JM = "joint survival and longitudinal model"
  )

  if (upper)
    substr(a, 1, 1) <- toupper(substr(a, 1, 1))
  a
}



# used in summary() (2020-06-10)
get_intercepts <- function(stats, varname, lvls, rev = FALSE) {
  # format the output for the intercepts in ordinal models
  # - stats: matrix of posterior summaries
  # - varname: name of the ordinal outcome variable
  # - lvls: levels of the ordinal factor

  interc <- stats[grep(paste0("gamma_", varname), rownames(stats)), ,
                  drop = FALSE]
  attr(interc, "rownames_orig") <- rownames(interc)

  if (isTRUE(rev)) {
    rownames(interc) <- paste(varname, "\u2264", lvls[-length(lvls)])
  } else {
    rownames(interc) <- paste(varname, ">", lvls[-length(lvls)])
  }
  interc
}


# used in summary() (2020-06-10)
computeP <- function(x) {
  # calculate tail probability
  # - x: vector of MCMC samples for one parameter

  above <- mean(x > 0)
  below <- mean(x < 0)
  2 * min(above, below)
}