R/tssem_variables.R

In JConigrave/msemtools: Routines, tables, and figures for metaSEM analyses

#' cormat_list
#'
#' Generate a list of correlation matrices for tssem
#' @param yi correlation effect
#' @param vi correlation variance
#' @param var1 the name of the column containing the first variable
#' @param var1 the name of the column containing the second variable
#' @param cluster the name of the clustering variable
#' @param data data.frame
#' @param transform a function to transform individual correlations
#' @param ni the name of the column containing the sample size
#' @import data.table

cormat_list = function(yi, vi, ni, var1, var2, cluster, data, transform = NULL){

  data = data.table::data.table(data)
  inspect_vars = data[,c(yi,vi,var1,var2,cluster), with = FALSE]
  missing = unlist(lapply(1:nrow(inspect_vars), function(x) any(is.na(inspect_vars[x,]))))
  data = data[!missing, ]
  vars <- unique(c(data[,var1], data[,var2]))

  make_matrix = function(vars) {
    m = matrix(nrow = length(vars), ncol = length(vars))
    rownames(m) = vars
    colnames(m) = vars

    diag(m) = 1
    m
  }

  get_matrix = function(id, vars, data, transform = NULL){
    m = make_matrix(vars)
    dat = data[which(data[,cluster, with = FALSE] == id),]

    for(col in seq_along(colnames(m))){
      for(row in seq_along(rownames(m))){
        # message("row = ",row, "col = ", col)
        if(col == row){
          m[row, col] = 1
          next
        }

        vars = c(colnames(m)[col], rownames(m)[row])
        valid_rows = unlist(dat[,var1, with = FALSE]) %in% vars &
          unlist(dat[,var2, with = FALSE]) %in% vars &
          !is.na(unlist(dat[,yi, with = FALSE])) &
          !is.na(unlist(dat[,vi, with = FALSE]))

        if(sum(valid_rows) == 0){
          m[row,col] = NA
          next
        }

        y <- dat[valid_rows, yi, with = FALSE]
        v <- dat[valid_rows, vi, with = FALSE]

        res <- stats::weighted.mean(y, 1/v)

        if(!is.null(transform)){
          res <- transform(res)
        }
        m[row,col] = res
      }

    }
    rownames(m) = rownames(m)
    colnames(m) = colnames(m)
    m

  }

  unique_ids = unique(unlist(data[,cluster, with = FALSE]))

  dat_list = list()
  dat_list$data = lapply(unique_ids, function(x) get_matrix(x, vars, data, transform))
  names(dat_list$data) <- unique_ids

  get_n = function(id, ni){
    ceiling(mean(unlist(data[unlist(data[,cluster, with = FALSE]) == id, ni, with = FALSE]),na.rm = TRUE))
  }


  dat_list$n = unlist(lapply(unique_ids, function(x) get_n(x, ni)))
  dat_list

}


star_matrix <- function(m, stars) {
  get_stars = function(p, stars) {
    if (is.na(p))
      p <- 1
    n_stars = sum(p < stars)
    paste(rep("*", n_stars), collapse = "")
  }

  s_matrix = m
  s_matrix[] =  sapply(m, function(p)
    get_stars(p, stars = stars))
  return(s_matrix)
}

#' tssem1_table
#'
#' Creates a table from the pooled correlation matrix from tssem1
#' @param model the tssem stage 1 model

tssem1_table = function(model){

  dim_names = model$original.names
  r_mat = coef(model, select = "fixed") %>%
    metaSEM::vec2symMat(diag = FALSE)
  dimnames(r_mat) = list(dim_names,dim_names)

  coefs = summary(model)$coefficients # get all coefs
  p_mat = coefs[grepl("Intercept", rownames(coefs)), "Pr(>|z|)"] %>% metaSEM::vec2symMat(diag = FALSE) #keep intercept coefs (not taus)
  p_mat[upper.tri(p_mat)] = 1 # I don't want stars for the diagonal.
  s_mat = star_matrix(p_mat , stars = c(0.05,0.01,0.001)) #star matrix

  i2 = summary(model)$I2.values[,"Estimate"]*100
  i2_mat = metaSEM::vec2symMat(i2, diag = FALSE) # create i2 matrix
  #r_mat = round(r_mat,2) #round r matrix
  r_mat[upper.tri(r_mat)] = i2_mat[upper.tri(i2_mat)] # merge in i2 matrix
  r_mat = papyr::digits(r_mat, 2)
  r_mat[] = paste0(r_mat,s_mat) #add in significance stars
  diag(r_mat) = "-" # add in diagonal
  r_mat[lower.tri(r_mat)] = gsub("0\\.",".",r_mat[lower.tri(r_mat)]) #get rid of leading zeros
  r_mat = data.frame(r_mat)
  rownames(r_mat) = paste0(seq_along(rownames(r_mat)),". ", rownames(r_mat)) #change rownames
  colnames(r_mat) = seq_along(rownames(r_mat)) #change colnames
  r_mat
}

#' tssem2_table
#'
#' Tabulate tssem2 regressions
#' @param wls a wls model
#' @param ... recode variable names, new_name = old_name
#' @param estimate string, name for estimate variable
#' @param transf transform function to apply to results
#' @param t.name name for transformed results
#' @export

tssem2_table = function(wls, ..., transf = NULL, t.name = NULL, estimate = "Estimate", round = 2){

  epi <- list(...)

  reg <- summary(wls)$coefficients
  vars <- colnames(wls$Cov)
  temp_vars <- paste0("v", 1:length(vars))


  reg <- tibble::rownames_to_column(reg)

  for(i in seq_along(vars)){
    reg$rowname = gsub(vars[i],temp_vars[i],reg$rowname)
  }

  # allow recodes -------------------------

  if(length(epi) > 0){
    for(i in seq_along(epi)){
      vars[grepl(epi[i],vars)] = names(epi[i])
    }
  }

  # ---------------------------------------


  colnames(reg) = c("Predictor", "Estimate","SE","lbound","ubound","z", "p")

  if(!is.null(transf)){
    reg$Estimate <- papyr::digits(transf(reg$Estimate),2)
    reg$lbound <- papyr::digits(transf(reg$lbound),2)
    reg$ubound <- papyr::digits(transf(reg$ubound),2)
    reg$transf <- glue::glue("{reg$Estimate} [{reg$lbound}, {reg$ubound}]")
    reg = reg[,c(1, 8, 2:7)]
  }else{
    est <- papyr::digits(reg$Estimate,2)
    lower <- papyr::digits(reg$lbound,2)
    upper <- papyr::digits(reg$ubound,2)
    reg$Estimate <- glue::glue("{est} [{lower}, {upper}]")
  }

  reg$lbound = NULL
  reg$ubound = NULL

  reg$outcome = gsub("on.*","",reg$Predictor)
  reg$outcome[grepl("with", reg$Predictor)] = "Covariances"
  reg$Predictor = gsub(".*on","", reg$Predictor)

  repLace <- function(x) {
    n <- which(sapply(temp_vars, function(i)
      grepl(i, x)))

    if(sum(n) == 0){
      return("Covariances")
    }

    if(length(n) > 1) {
      one <- vars[n[1]]
      two <- vars[n[2]]
      out <- paste(one, "with", two)
    }else{
      out <- vars[n]
    }
    out
  }

  reg$SE = papyr::digits(reg$SE, 2)
  reg$z = papyr::digits(reg$z, 2)
  reg$p = papyr::round_p(reg$p)

  reg$Predictor <- sapply(reg$Predictor, repLace)
  reg$outcome <- sapply(reg$outcome, repLace)
  papyr:::to_rowhead(reg, outcome)


}

#' report_tssem2
#'
#' Report fit of tssem2
#' @param x wls model
#' @export

report_tssem2 <- function(x, pattern = "The model had {assessment} fit $\\chi^2$({df}) = {chi}; $p$ = {p}; RMSEA = {RMSEA}; SRMR = {SRMR}; TLI = {TLI}") {
  if(!is(x, "wls")) {
    stop("Can only be used with wls objects")
  }
  x <- summary(x)
  l = as.list(x$stat)
  names(l) = rownames(x$stat)

  if (l$`DF of target model` == 0) {
    mess = "The model had zero degrees of freedom and so had perfect fit."
    return(mess)
  }

  RMSEA <- l$RMSEA
  RMSEA_cat <- dplyr::case_when(RMSEA <= 0.01  ~ 1, # excellent (MacCallum, Browne, Sugawara)
                                RMSEA <= 0.05  ~ 2, # good
                                RMSEA <= 0.08  ~ 3, # mediocre
                                RMSEA > 0.08   ~ 4) # poor
  TLI <- l$TLI
  TLI_cat <- dplyr::case_when(TLI > 0.95  ~ 1, # excellent
                              TLI >= .95  ~ 2, # good
                              TLI >= 0.9  ~ 3, # mediocre
                              TLI <  0.9  ~ 4) # poor
  CFI <- l$CFI
  CFI_cat <- dplyr::case_when(CFI > 0.95  ~ 1, # excellent
                              CFI >= .95  ~ 2, # good
                              CFI >= 0.9  ~ 3, # mediocre
                              CFI <  0.9  ~ 4) # poor
  SRMR <- l$SRMR
  SRMR_cat <- dplyr::case_when(SRMR <= 0.04  ~ 1, # excellent (MacCallum, Browne, Sugawara)
                               SRMR <= 0.08  ~ 2, # good
                               SRMR <= 0.1   ~ 3, # mediocre
                               SRMR > .1     ~ 4) # poor



  assessment <- round(mean(c(TLI_cat, RMSEA_cat, SRMR_cat)),0)
  assessment <- c("excellent", "good","mediocre", "poor")[assessment]

  chi = papyr::digits(l$`Chi-square of target model`,2)
  df = papyr::digits(l$`DF of target model`,0)
  p = papyr::round_p(l$`p value of target model`)

  RMSEA = papyr::digits(RMSEA, 3)
  TLI = papyr::digits(TLI, 3)
  CFI = papyr::digits(CFI, 3)
  SRMR = papyr::digits(SRMR,3)

  mess <- glue::glue(pattern)
  mess
}