R/fit_model.r

In mars: Meta Analysis and Research Synthesis

#' Path Model Function
#'
#' This function fits the path model and returns adjusted standard errors.
#'
#' @param mars_object The mars fitted object.
#' @param model This is model syntax specified in the format by lavaan.
#' @param num_obs Number of observations
#' @param adjust_se Adjust the standard errors to reflect the ...
#' @param method_null Unsure
#' @param ... Currently not used.
#'
#' @details
#' The input is the coefficients and the variance covariance matrix returned
#' from the \link{mars} function.
#'
#'
#' @returns List output with class path; The output is the parameter estimates
#'  from the fitted path model.
#'
#' @export
path_model <- function(mars_object, model, num_obs = NULL, adjust_se = TRUE,
                       method_null = 'sem', ...) {

  if(is.null(num_obs)) {
    N <- mars_object[['sample_size']]
  } else{
    N <- num_obs
  }

  model_stats <- extract_model(mars_object,
                               variable_names = mars_object[['variable_names']])

  fitted_model <- model_fit(model_input = model, R = model_stats[['beta_matrix']],
                            method_null = method_null,
                            N = N)

  coefficients <- format_coefficients(fitted_model[['path_coefficients']])


  if(adjust_se) {

    num_outcomes <- length(coefficients)
    variables <- lapply(coefficients, unique_names)

    data_list <- lapply(seq_along(variables),
                        function(xx)
                          matrix_subset(model_stats[['beta_matrix']], variables[[xx]])
    )

    corr_pairs <- lapply(seq_along(data_list), function(xx)
      data.frame(t(utils::combn(colnames(data_list[[xx]]), 2))))

    corr_vector <- lapply(seq_along(data_list), function(xx)
      corpcor::sm2vec(data_list[[xx]]))

    corr_elements <- lapply(seq_along(data_list), function(xx)
      cbind(corr_pairs[[xx]], corr = corr_vector[[xx]]))

    full_corr_matrix <- cbind(data.frame(t(utils::combn(colnames(model_stats[['beta_matrix']]), 2))),
                              corr = corpcor::sm2vec(model_stats[['beta_matrix']]))
    full_corr_matrix[['num']] <- 1:nrow(full_corr_matrix)

    matches <- lapply(seq_along(corr_elements), function(xx)
      merge(corr_elements[[xx]], full_corr_matrix,
            by = c('X1', 'X2', 'corr')))

    matches <- lapply(seq_along(matches), function(xx)
      matches[[xx]][order(matches[[xx]][['num']]), ]
      )

    var_list <- lapply(seq_along(variables),
                       function(xx)
                         model_stats[['var_matrix']][matches[[xx]][['num']],
                                              matches[[xx]][['num']]]
    )

    computed_se <- lapply(seq_along(data_list), function(xx)
      var_path(data_list[[xx]],
               var_list[[xx]],
               type = 'stdslopes')
    )
  } else {
    computed_se <- NULL
  }

  model_output <- list(beta_matrix = model_stats[['beta_matrix']],
                       parameter_estimates = coefficients,
                       fit_measures = fitted_model,
                       computed_var = computed_se,
                       model = model,
                       est_tau = mars_object[['est_values']][['Tau']])

  class(model_output) <- 'path'

  model_output
}

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

  fixed_coef <- do.call('rbind', object[['parameter_estimates']])

  standard_errors <- lapply(seq_along(object[['computed_var']]), function(xx)
                            sqrt(diag(object[['computed_var']][[xx]])))

  fixed_coef[['standard_errors']] <- do.call('c', standard_errors)
  fixed_coef[['test_statistic']] <- fixed_coef[['estimate']] / fixed_coef[['standard_errors']]
  fixed_coef[['p_value']] <- stats::pnorm(abs(fixed_coef[['test_statistic']]), lower.tail = FALSE) * 2

  if(fit_measures) {
    fit_meas <- object[['fit_measures']]
    fit_meas[['path_coefficients']] <- NULL
  } else {
    fit_meas <- NULL
  }

  structure(list(beta_matrix = data.frame(object[['beta_matrix']]),
                 model = object[['model']],
                 fit_measures = fit_meas,
                 covariance = object[['est_tau']],
                 fixed_coef = fixed_coef
                 ), class = "summary.path")

}

#' @export
print.summary.path <- function(x, digits = max(3, getOption("digits") - 3),
                                   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')
      cat('multivariate')

    cat('\n \nAverage Correlation Matrix: \n')
    print(x[['beta_matrix']], digits = x[['digits']])
    cat('\n \n Model Fitted: \n', x[['model']], "\n \n")
        cat('Fixed Effects: \n')
        print(x[['fixed_coef']], digits = x[['digits']])
        cat('\n \n Fit Statistics: \n')
        print(format_fit_stats(x[['fit_measures']]))
  }
}

format_fit_stats <- function(object, digits = 3) {
  data.frame(
    Type = c('Model Chi-Square', 'Null Model Chi-Square',
             'CFI', "TLI", "RMSEA", "SRMR"),
    Value = c(paste0(round(object[['Model']][['Chi2']], digits = digits), " (", object[['Model']][['df']], "), ", object[['Model']][['pvalue']]),
              paste0(round(object[['NullModel']][['Chi2Null']], digits = digits), " (", object[['NullModel']][['dfNull']], ")"),
              round(object[['CFI']], digits = digits),
              round(object[['TLI']], digits = digits),
              paste0(round(object[['RMSEA']], digits = digits), " [", paste(round(object[['RMSEA_CI']], digits = digits), collapse = ", "), "]"),
              round(object[['SRMR']], digits = digits))
  )
}