R/mars_main.r

In mars: Meta Analysis and Research Synthesis

#' Master mars function
#'
#' The primary function used for input and estimation. The function takes the
#' data inputs and routes the estimation and structure type based on data structure.
#' The function can handle univariate, multivariate, longitudinal, and multilevel meta-analytic models.
#'
#' @param data Data used for analysis
#' @param studyID Character string representing the study ID
#' @param effectID Character string representing the effect size ID
#' @param sample_size Character string representing the sample size of the studies.
#' @param effectsize_type Type of effect size being analyzed
#' @param formula The formula used for specifying the fixed and random structure.
#' Used for univariate and multilevel structures.
#' @param variable_names Vector of character strings representing the attributes with
#'   correlations. The attributes that are correlated should be separated by an
#'   underscore.
#' @param effectsize_name Character string representing the name of the effect size
#'   column in the data.
#' @param estimation_method Type of estimation used, either "REML" or "MLE", REML is the default
#' @param variance Character string representing the name of the variance of the effect size in the data.
#' @param varcov_type Type of variance covariance matrix computed. Default is
#' 'cor_weighted' for correlations or 'smd_outcome' for standardized mean differences.
#' @param weights User specified matrix of weights for analysis.
#' @param structure Between studies covariance structure, default is "UN" or unstructured. See details for more specifics.
#' @param intercept Whether a model intercept should be specified, default is FALSE
#'  meaning no intercept. See details for more information.
#' @param missing Whether missing data should be removed, or kept. Default is
#' removing.
#' @param optim_method Optimization method that is passed to the optim function.
#'  Default is 'L-BFGS-B'.
#' @param robustID A character vector specifying the cluster group to use for computing the robust standard errors.
#' @param multivariate_covs A one-sided formula to specify the covariates used in a multivariate analysis.
#' @param tol Tolerance of the optimization, default is 1E-10.
#' @param ... Not currently used.
#'
#' @importFrom Matrix bdiag sparseMatrix sparseVector
#' @importFrom corpcor sm2vec
#' @importFrom matrixcalc elimination.matrix svd.inverse vec
#' @importFrom utils combn
#'
#' @returns Returns a list of class mars; The returned object contains elements from the
#' estimation.
#'
#' @export
mars <- function(data, studyID, effectID, sample_size, effectsize_type = NULL,
                 formula = NULL,
                 variable_names = NULL, effectsize_name = NULL,
                 estimation_method = 'REML', variance = NULL,
                 varcov_type, weights = NULL,
                 structure = 'UN', intercept = FALSE,
                 missing = 'remove', optim_method = "L-BFGS-B",
                 robustID = NULL, multivariate_covs = NULL,
                 tol = 1E-10, ...) {

  if(!is.null(formula)) {
    if(structure == 'univariate') {
      est <- mars_univariate(data = data, formula = formula,
                             studyID = studyID,
                             variance = variance, varcov_type = varcov_type,
                             weights = weights,
                             structure = structure, estimation_method = estimation_method,
                             missing = missing, optim_method = optim_method,
                             robustID = robustID,
                             tol = tol)
    }
    if(structure == 'multilevel') {
      est <- mars_multilevel(data = data, formula = formula,
                             studyID = studyID,
                             variance = variance, varcov_type = varcov_type,
                             weights = weights,
                             structure = structure, estimation_method = estimation_method,
                             missing = missing, optim_method = optim_method,
                             robustID = robustID,
                             tol = tol)
    }
  }

  if(!is.null(effectsize_type)) {
    if(effectsize_type == 'cor') {

      est <- mars_cor(data = data, studyID = studyID, effectID = effectID,
                      sample_size = sample_size, variable_names = variable_names,
                      estimation_method = estimation_method, varcov_type = varcov_type,
                      structure = structure, intercept = intercept,
                      missing = missing, optim_method = optim_method,
                      robustID = robustID, multivariate_covs = multivariate_covs,
                      tol = tol)
    }
    if(effectsize_type == 'smd') {
      est <- mars_smd(data = data, studyID = studyID, effectID = effectID,
                      sample_size = sample_size, effectsize_name = effectsize_name,
                      estimation_method = estimation_method, varcov_type = varcov_type,
                      structure = structure, intercept = intercept,
                      missing = missing, optim_method = optim_method,
                      robustID = robustID, multivariate_covs = multivariate_covs,
                      tol = tol)
    }
  }

  class(est) <- 'mars'

  return(est)
}

mars_cor <- function(data, studyID, effectID, sample_size, variable_names,
                     estimation_method = 'REML', varcov_type,
                     structure = "UN", intercept = FALSE,
                     missing = 'remove', optim_method = "L-BFGS-B",
                     robustID = NULL, multivariate_covs = NULL,
                     tol = 1E-10, ...) {
  data_input <- prep_data(data, studyID = studyID, samp_size_name = sample_size,
                          variable_names = variable_names)

  q <- length(unique(data_input[['numID']]))

    mod <- estimation(effect_name = 'ri', data = data_input,
                      studyID = 'study',
                      effectID = 'numID',
                      estimation_method = estimation_method,
                      varcov_type = varcov_type,
                      structure = structure,
                      intercept = intercept,
                      N = sample_size, missing = missing,
                      optim_method = optim_method,
                      robustID = robustID, multivariate_covs = multivariate_covs,
                      tol = tol)

    c(mod,
         list(variable_names = variable_names,
         sample_size = sum(data[[sample_size]])))

}

mars_smd <- function(data, studyID, effectID, sample_size,
                     effectsize_name,
                     estimation_method = 'REML', varcov_type,
                     structure = 'UN', intercept = FALSE,
                     missing = 'remove', optim_method = "L-BFGS-B",
                     robustID = NULL, multivariate_covs = NULL,
                     tol = 1E-10, ...) {

  mod <- estimation(effect_name = effectsize_name, data = data,
             studyID = studyID,
             effectID = effectID,
             estimation_method = estimation_method,
             varcov_type = varcov_type,
             structure = structure,
             intercept = intercept,
             N = sample_size, missing = missing,
             optim_method = optim_method,
             robustID = robustID, multivariate_covs = multivariate_covs,
             tol = tol)

  return(mod)

}

mars_univariate <- function(data, formula, studyID,
                            variance, varcov_type,
                            weights,
                            structure, estimation_method,
                            missing = 'remove', optim_method = "L-BFGS-B",
                            robustID = NULL,
                            tol = 1E-10, ...) {

  mod <- estimation(formula = formula,
                    data = data, studyID = studyID,
                    variance = variance, varcov_type = varcov_type,
                    weights = weights,
                    structure = structure,
                    estimation_method = estimation_method,
                    missing = missing, optim_method = optim_method,
                    robustID = robustID,
                    tol = tol)

  return(mod)
}

mars_multilevel <- function(data, formula, studyID,
                            variance, varcov_type,
                            weights,
                            structure, estimation_method,
                            missing = 'remove', optim_method = "L-BFGS-B",
                            robustID = NULL,
                            tol = 1E-10, ...) {

  mod <- estimation(formula = formula,
             data = data, studyID = studyID,
             variance = variance, varcov_type = varcov_type,
             weights = weights,
             structure = structure,
             estimation_method = estimation_method,
             missing = missing, optim_method = optim_method,
             robustID = robustID,
             tol = tol)

  return(mod)
}

#' @importFrom stats cov2cor
#' @importFrom stats aggregate median
#' @export
summary.mars <- function(object, fit_measures = TRUE,
                         digits = max(3, getOption("digits") - 3),
                         ...) {

  if(object[['estimator']] != 'FE') {
    if(is.null(object[['random_names']])) {
      if(object[['structure']] == 'univariate') {
        correlations <- data.frame(term = 'intercept',
                                   var = diag(object[['est_values']][['Tau']]),
                                   SD = sqrt(diag(object[['est_values']][['Tau']]))
        )
      } else {
        correlations <- cov2cor(object[['est_values']][['Tau']])

        names_dim <- paste0(object[['effect_size_name']], "_", gsub("^factor\\(.+\\)", "",
                                                                    colnames(object[['design_matrix']][1:object[['q_r']], 1:object[['q_r']]])))

        correlations[upper.tri(correlations, diag = TRUE)] <- object[['est_values']][['Tau']][upper.tri(object[['est_values']][['Tau']], diag = TRUE)]
        colnames(correlations) <- rownames(correlations) <- names_dim
      }
    } else {
      correlations <- data.frame(term = object[['random_names']],
                                 var = diag(object[['est_values']][['Tau']]),
                                 SD = sqrt(diag(object[['est_values']][['Tau']]))
      )
    }
  } else {
    correlations <- NULL
  }


  if(is.null(object[['beta_r']])) {
    names_dim <- paste0(object[['effect_size_name']], "_", gsub("^factor\\(.+\\)", "",
                                                                colnames(object[['design_matrix']])))
    beta <- cbind.data.frame(names_dim,
                                   object[['beta_fe']])
    fixed_effect <- TRUE
  } else {
    names_dim <- paste0(object[['effect_size_name']], "_", gsub("^factor\\(.+\\)", "",
                                                                colnames(object[['design_matrix']])))
    beta <- cbind.data.frame(names_dim,
                          object[['beta_r']])
    fixed_effect <- FALSE
  }

  names(beta) <- c("attribute", 'estimate')
  beta <- cbind(beta, SE = sqrt(diag(object[['varcov_beta']])))
  beta <- cbind(beta, z_test = beta[['estimate']] / beta[['SE']])
  beta <- cbind(beta, p_value = pnorm(abs(beta[['z_test']]), lower.tail = FALSE) * 2,
                lower = beta[['estimate']] + qnorm(.025, lower.tail = TRUE) * beta[['SE']],
                upper = beta[['estimate']] + qnorm(.975, lower.tail = TRUE) * beta[['SE']])

  if(!is.null(object[['robust_se']])) {
    beta_robust <- data.frame(beta[, c("attribute", 'estimate')])
    beta_robust <- cbind(beta_robust, SE = sqrt(diag(matrix(object[['robust_se']][['vb_1']]))))
    beta_robust <- cbind(beta_robust, z_test = beta_robust[['estimate']] / beta_robust[['SE']])
    beta_robust <- cbind(beta_robust, p_value = pnorm(abs(beta_robust[['z_test']]), lower.tail = FALSE) * 2,
                  lower = beta_robust[['estimate']] + qnorm(.025, lower.tail = TRUE) * beta_robust[['SE']],
                  upper = beta_robust[['estimate']] + qnorm(.975, lower.tail = TRUE) * beta_robust[['SE']])
  } else{
    beta_robust = NULL
  }
  if(is.null(object[['beta_r']])) {
    dimensions <- list(
      number_random = 0,
      number_fixed = length(object[['beta_fe']]),
      number_effectsize = nrow(object[['design_matrix']]),
      dim_random = unlist(object[['dim_random']])
    )
  } else {
    dimensions <- list(
      number_random = nrow(correlations),
      number_fixed = length(object[['est_values']][['mu']]),
      number_effectsize = nrow(object[['design_matrix']]),
      dim_random = unlist(object[['dim_random']])
    )
  }

  if(!is.null(object[['random_names']])) {

    #min_length_r <- object[['dim_random']][which.min(object[['dim_random']])]
    length_r <- aggregate(object[['random_data']],
                          by = list(object[['random_data']][, names(object[['random_data']][which.min(object[['dim_random']])])]),
                          length)[, 1:2]
    names(length_r) <- c(names(object[['dim_random']][which.min(object[['dim_random']])]),
                         'length')

    random_stats <- list(
      min_length = min(length_r[['length']]),
      med_length = median(length_r[['length']]),
      max_length = max(length_r[['length']])
    )
  } else {
    random_stats <- NULL
  }

  object[['QE']][['p_r']] <- ifelse(object[['QE']][['p']] < .0001,
                                  'p < 0.0001', paste0('p = ', round(object[['QE']][['p']], 4)))

  if(length(object[['I2_within']]) > 1) {
    if(object[['structure']] == 'multilevel') {
      object[['I2_within']] <- data.frame(names = names(dimensions[['dim_random']]),
                                          values = object[['I2_within']])
    } else {
      names_dim <- paste0(object[['effect_size_name']], "_", gsub("^factor\\(.+\\)", "",
                                                                  colnames(object[['design_matrix']][1:object[['q_r']], 1:object[['q_r']]])))
      object[['I2_within']] <- data.frame(names = names_dim,
                                          values = object[['I2_within']])
    }
  }
  if(object[['estimator']] == 'FE') {
    i2_fe <- object[['QE']][['value']] / object[['QE']][['df']]
    i2_fe_part2 <- as.numeric(((i2_fe - 1) / i2_fe) * 100)
  } else {
    i2_fe_part2 <- NULL
  }
  if(length(object[['I2_jackson']]) > 1) {
    if(object[['structure']] == 'multilevel') {
      object[['I2_jackson']] <- data.frame(names = names(dimensions[['dim_random']]),
                                          values = object[['I2_jackson']])
    } else {
      names_dim <- paste0(object[['effect_size_name']], "_", gsub("^factor\\(.+\\)", "",
                                                                  colnames(object[['design_matrix']])))
      object[['I2_jackson']] <- data.frame(names = names_dim,
                                          values = object[['I2_jackson']])
    }
  }



  # Compute number of studies within the higher level

  structure(list(formula = object[['formula']],
                 structure = object[['structure']],
                 tau2_cor = correlations,
                 random_stats = random_stats,
                 fixed_effect = fixed_effect,
                 beta = beta,
                 beta_robust = beta_robust,
                 #model = object[['model']],
                 fit_measures = object[['fit_stats']],
                 QE = object[['QE']],
                 QE_w = object[['QE_w']],
                 QM = object[['QM']],
                 dimensions = dimensions,
                 digits = digits,
                 estimator = object[['estimator']],
                 #variance = variances,
                 #covariance = covariances,
                 I2 = list('I2_within' = object[['I2_within']],
                           'I2_jackson' = object[['I2_jackson']],
                           'I2_between' = object[['I2_between']]
                           ),
                 I2_fixed = i2_fe_part2
  ), class = "summary.mars")


}

#' @importFrom utils packageVersion
#' @export
print.summary.mars <- function(x,
                               signif.stars = getOption("show.signif.stars"),
                               tidy = FALSE, ...) {

  if(tidy) {

  } else {
    cat("Results generated with MARS:")
    cat(paste('v', packageVersion("mars"), '\n'))
    cat(format(Sys.time(), "%A, %B %e, %Y"))
    cat('\n\nModel Type: \n')
    if(x[['structure']] %in% c('univariate', 'multilevel')) {
      cat(paste0(x[['structure']]))
    } else {
      cat('multivariate')
    }
    cat("\n\nEstimation Method: \n")
    if(x[['estimator']] == 'REML') {
      cat("Restricted Maximum Likelihood")
    } else if(x[['estimator']] == 'MLE') {
      cat("Maximum Likelihood")
    } else {
      cat("Fixed/Common Effect")
    }
    cat('\n\nModel Formula: \n')
    print(x[['formula']], digits = x[['digits']])
    cat('\nData Summary: \n')
    cat('Number of Effect Sizes: ')
          cat(x[['dimensions']][['number_effectsize']])
          cat('\nNumber of Fixed Effects: ')
          cat(x[['dimensions']][['number_fixed']])
          cat('\nNumber of Random Effects: ')
          cat(x[['dimensions']][['number_random']])
    if(!is.null(x[['random_stats']])) {
      cat("\n\nRandom Data Summary: \n")
      cat(paste(lapply(seq_along(x[['dimensions']][['dim_random']]),
                   function(xx) paste0("Number unique ", names(x[['dimensions']][['dim_random']][xx]), ": ", x[['dimensions']][['dim_random']][xx])),
            collapse = " \n"))
    }
    cat('\n\nRandom Components: \n')
    print(x[['tau2_cor']], row.names = FALSE, digits = x[['digits']])
    # cat('Variance Estimates: \n', x$variance, "\n \n")
    # cat('Covariance Estimates: \n', x$covariance, "\n \n")
    cat('\nFixed Effects Estimates: \n')
    print(x[['beta']], row.names = FALSE, digits = x[['digits']])
    if(!is.null(x[['beta_robust']])) {
      cat('\n\nFixed Effects Estimates with Robust SEs: \n')
      print(x[['beta_robust']], row.names = FALSE, digits = x[['digits']])
    }
    # data.frame(x$fixed_coef[, 1:2], printCoefmat(x$fixed_coef[, 3:ncol(x$fixed_coef)], digits = digits, signif.stars = signif.stars,
    #              has.Pvalue = TRUE, cs.ind = 1:2, tst.ind = 3, P.values = TRUE))
    cat("\nModel Fit Statistics: \n")
    print(x[['fit_measures']], row.names = FALSE, digits = x[['digits']])
    cat("\nQ Error: ")
    cat(paste0(round(x[['QE']][['value']], 3), ' (', x[['QE']][['df']], "), ", x[['QE']][['p_r']]))
    cat("\n\nI2 (General): ")
    if(x[['estimator']] == 'FE') {
      cat("\n")
      cat(round(x[['I2_fixed']], x[['digits']]))
    } else {
      if(x[['structure']] != 'univariate') {
        cat("\n")
        print(x[['I2']][['I2_within']], row.names = FALSE, digits = x[['digits']])
        cat("\nI2 (Jackson): ")
        if(is.data.frame(x[['I2']][['I2_jackson']])) {
          cat("\n")
          print(x[['I2']][['I2_jackson']], row.names = FALSE, digits = x[['digits']])
        } else {
          cat(x[['I2']][['I2_jackson']])
        }
        cat("\nI2 (Between): ")
        cat(x[['I2']][['I2_between']])
      } else {
        cat(x[['I2']][['I2_within']])
      }
    }
  }

}