R Notebook

Extract model params

library(MplusAutomation)
# function for table with all params from DSEM models, 2months and one year
# mplodelnames: names of models for single events (situations) to fill the table with; adds1models: names of additional models for single events (situations) to fill the table with, if second model no confidence intervals; tablenames: names of columns, orders = rownames of  parameters of interest that should be extracted from mplus output; params: rownames/names of parameters for table; html: Output table in html; logtransparams: calculate betas for moodxevent based on hamaker and olivoto procedure; predtable: is it a table with prediction using between person independent variables or is it just within; noconf: no confidence intervals
table2func = function(mplodelnames,
                      adds1models = NA,
                      tablenames,
                      orders,
                      orders1,
                      params = NA ,
                      html = T,
                      logtransparams = F,
                      predtable = F,
                      noconf = NA) {
  if (length(adds1models) == 1) {
    # read mplus model
    mplodel = readModels(target = mplodelnames[1])
    #extract std params
    stdout = mplodel$parameters$stdyx.standardized
    #delete unnecessary rows
    stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
    stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
    rownames(stdout) = 1:dim(stdout)[1]
    temp = stdout
    # reorder within, between, pred and extract only variables we are interested in
    stdout = stdout[orders, ]
    # check available N after removal of cases wo variance and or
    nocases = unlist(regmatches(unlist(mplodel$warnings), gregexpr('[0-9]+', unlist(
      mplodel$warnings
    ))))
    nocases = nocases[nchar(nocases) > 1]
    N = mplodel$data_summary$overall$NClusters - length(nocases)
    # create table
    table2 = stdout[, c("paramHeader", "est")]
    names(table2) = c("param", paste0(tablenames[1], " N=", N, "<br> stdyx."))
    if (length(is.na(params)) == 1) {
      #put params
      table2$param = paste(stdout$paramHeader, stdout$param)
    } else{
      table2$param = params
    }
    # create table with boldness when nonoverlapping ci
    if (predtable == T) {
      tablebolds = table2
      tablebolds[, 2] = FALSE
      # bold for non-overlapping PF and Dep conf int
      i = 1
      while (i < dim(stdout)[1]) {
        a1 = stdout$lower_2.5ci[i]
        a2 = stdout$upper_2.5ci[i]
        b1 = stdout$lower_2.5ci[i + 1]
        b2 = stdout$upper_2.5ci[i + 1]
        if (!((a1 <= b2) && (b1 <= a2))) {
          tablebolds[i, 2] = T
          tablebolds[i + 1, 2] = T
        }
        i = i + 2
      }
    } else{
      tablebolds = table2
      tablebolds[, 2:ncol(tablebolds)] = F
    }
    if (logtransparams == T) {
      #old procedure
      #  rawout =  mplodel$parameters$unstandardized
      # rawout = rawout[grep("ntercept",rawout$paramHeader) ,]
      # # exchange m_s, miv or siv params with log-backtransformed values from raw output
      # rawout = rawout[rawout$param %in% c("MIV","SIV","M_S"),]
      # stdout[grep("M_S_LV", stdout$param),"est"] = exp(rawout$est[rawout$param == "M_S"]) / (sqrt(exp(rawout$est[rawout$param == "MIV"])+exp(rawout$est[rawout$param == "M_S"]))*sqrt(exp(rawout$est[rawout$param == "SIV"])+exp(rawout$est[rawout$param == "M_S"])))
      # # if negative loading
      # if(temp < 0){stdout[grep("M_S_LV", stdout$param),"est"] = stdout[grep("M_S_LV", stdout$param),"est"] *-1}
      ###### procedure advised by Ellen Hamaker  #######
      # unique part of residual variances
      loadm =   temp[temp$paramHeader == "M_S_LV.BY", "est"][1]
      loads =  temp[temp$paramHeader == "M_S_LV.BY", "est"][2]
      # covariance between residual variances
      covms = loads * loadm
      # residual variances
      varresidm = temp[temp$paramHeader == "MIV.|", "est"] + loadm ^ 2
      varresids = temp[temp$paramHeader == "SIV.|", "est"] + loads ^ 2
      # correlation = cov / sqrt(varresidm * varresids)
      stdout[grep("M_S_LV", stdout$param), "est"] = covms / sqrt(varresidm * varresids)
      #confidence interval according to Olivoto et al. (2018)
      CI = (abs(stdout[grep("M_S_LV", stdout$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
      stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] - CI
      stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] + CI
      # stdout[grep("SIV.", stdout$paramHeader),"est"] =exp(rawout$est[rawout$param == "SIV"])
      # stdout[grep("SIV.", stdout$paramHeader),"lower_2.5ci"] =exp(rawout$lower_2.5ci[rawout$param == "SIV"])
      # stdout[grep("SIV.", stdout$paramHeader),"upper_2.5ci"] =exp(rawout$upper_2.5ci[rawout$param == "SIV"])
      # if("&pi;m" %in% table2$param){
      #  stdout[grep("MIV.", stdout$paramHeader),"est"] =exp(rawout$est[rawout$param == "MIV"])
      #  stdout[grep("MIV.", stdout$paramHeader),"lower_2.5ci"] =exp(rawout$lower_2.5ci[rawout$param == "MIV"])
      # stdout[grep("MIV.", stdout$paramHeader),"upper_2.5ci"] =exp(rawout$upper_2.5ci[rawout$param == "MIV"])
      # }
    }
    
    # create noconf boolean
    if(length(noconf) == 1){ if(sum(is.na(noconf)) > 0) {noconf = rep(F,length(orders))} else if(noconf){noconf = rep(T,length(orders))}} 
  if (html == T){
        table2[!noconf, 2] = paste0(
          "<span style=\"",
          ifelse(stdout$pval[!noconf] < .005, yes = "font-weight: bold;", no = ";"),
          "\" >",
          round(stdout$est[!noconf], digits = 2),
          " <br> [",
          round(stdout$lower_2.5ci[!noconf], digits = 2),
          " ",
          round(stdout$upper_2.5ci[!noconf], digits = 2),
          "]</span>"
        )}
       if (html == F) {
        table2[!noconf, 2] = ifelse(
          stdout$pval < .005,
          yes = paste0(
            round(stdout$est[!noconf], digits = 2),
            " <br> [",
            round(stdout$lower_2.5ci[!noconf], digits = 2),
            " ",
            round(stdout$upper_2.5ci[!noconf], digits = 2),
            "]"
          ),
          no = ""
        )
      }
   
      table2[noconf, 2] = paste0(
        "<span style=\"",
        ifelse(stdout$pval[noconf] < .005, yes = "font-weight: bold;", no = ";"),
        "\" >",
        round(stdout$est[noconf], digits = 2),
        "</span>"
      )
    if (html == F) {
      table2[noconf, 2] = ifelse(stdout$pval[noconf] < .005, yes = paste0(round(stdout$est[noconf], digits = 2)), no = "")
    }
    #
    # table2[,2] = paste0("<span style=\"", ifelse(stdout$pval <.005,yes = "font-weight: bold; display:block;width:100%;background-color:yellow;", no = "display:block;width:100%;background-color:yellow;"),   "\" >",round(stdout$est,digits = 2)," <br> [",round(stdout$lower_2.5ci,digits = 2)," ",round(stdout$upper_2.5ci,digits = 2),"]</span>")
    #
    # table2[,2]= paste(round(stdout$est,digits = 2),paste0("<br>[",round(stdout$lower_2.5ci,digits = 2), " ",round(stdout$upper_2.5ci,digits = 2),"]"))
    # make table with significance
    table2p = table2[, 1:2]
    names(table2p)[2] = tablenames[1]
    table2p[, 2] = stdout$pval
    # make table with regressionweight
    table2r = table2[, 1:2]
    names(table2r)[2] = tablenames[1]
    table2r[, 2] = stdout$est
    # make table with all models
    for (i in 2:length(mplodelnames)) {
      # read mplus model
      mplodel = MplusAutomation::readModels(mplodelnames[i])
      #extract unstd params for miv, siv, m_e
      #extract std params
      stdout = mplodel$parameters$stdyx.standardized
      #delete unnecessary rows
      stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
      stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
      rownames(stdout) = 1:dim(stdout)[1]
      temp = stdout
      # reorder within, between, pred
      stdout = stdout[orders, ]
      # check available N after removal of cases wo variance and or
      nocases = unlist(regmatches(unlist(mplodel$warnings), gregexpr(
        '[0-9]+', unlist(mplodel$warnings)
      )))
      nocases = nocases[nchar(nocases) > 1]
      N = mplodel$data_summary$overall$NClusters - length(nocases)
      if (logtransparams == T) {
        ###### procedure advised by Ellen Hamaker  #######
        # unique part of residual variances
        loadm =   temp[temp$paramHeader == "M_S_LV.BY", "est"][1]
        loads =  temp[temp$paramHeader == "M_S_LV.BY", "est"][2]
        # covariance between residual variances
        covms = loads * loadm
        # residual variances
        varresidm = temp[temp$paramHeader == "MIV.|", "est"] + loadm ^ 2
        varresids = temp[temp$paramHeader == "SIV.|", "est"] + loads ^ 2
        # correlation = cov / sqrt(varresidm * varresids)
        stdout[grep("M_S_LV", stdout$param), "est"] = covms / sqrt(varresidm * varresids)
        #confidence interval according to Olivoto et al. (2018)
        CI = (abs(stdout[grep("M_S_LV", stdout$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
        stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] - CI
        stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] + CI
      }
      #put params
      #table2$param = paste(stdout$paramHeader, stdout$param)
        table2[!noconf, i + 1] = paste0(
          "<span style=\"",
          ifelse(stdout$pval[!noconf] < .005, yes = "font-weight: bold;", no = ";"),
          "\" >",
          round(stdout$est[!noconf], digits = 2),
          " <br> [",
          round(stdout$lower_2.5ci[!noconf], digits = 2),
          " ",
          round(stdout$upper_2.5ci[!noconf], digits = 2),
          "]</span>"
        )
       if (html == F) {
        table2[!noconf, i + 1] = ifelse(
          stdout$pval[!noconf] < .005,
          yes = paste0(
            round(stdout$est[!noconf], digits = 2),
            " <br> [",
            round(stdout$lower_2.5ci[!noconf], digits = 2),
            " ",
            round(stdout$upper_2.5ci[!noconf], digits = 2),
            "]"
          ),
          no = ""
        )
      }
   
      table2[noconf, i + 1] = paste0(
        "<span style=\"",
        ifelse(stdout$pval[noconf] < .005, yes = "font-weight: bold;", no = ";"),
        "\" >",
        round(stdout$est[noconf], digits = 2),
        "</span>"
      )
    if (html == F) {
      table2[noconf, i + 1] = ifelse(stdout$pval[noconf] < .005, yes = paste0(round(stdout$est[noconf], digits = 2)), no = "")
    }
      names(table2)[i + 1] = paste0(tablenames[i], " N=", N, "<br> stdyx.")
      # create table with boldness when nonoverlapping ci
      if (predtable == T) {
        # bold for non-overlapping PF and Dep conf int
        tablebolds[, i + 1] = FALSE
        d = 1
        while (d < dim(stdout)[1]) {
          a1 = stdout$lower_2.5ci[d]
          a2 = stdout$upper_2.5ci[d]
          b1 = stdout$lower_2.5ci[d + 1]
          b2 = stdout$upper_2.5ci[d + 1]
          if (!((a1 <= b2) && (b1 <= a2))) {
            tablebolds[d, i + 1] = T
            tablebolds[d + 1, i + 1] = T
          }
          d = d + 2
        }
      } else{
        tablebolds[, i + 1] = FALSE
      }
      # put in table with significance
      table2p[, i + 1] = stdout$pval
      names(table2p)[i + 1] = tablenames[i]
      # put in table with regression
      table2r[, i + 1] = stdout$est
      names(table2r)[i + 1] = tablenames[i]
    }
    rownames(table2) = NULL
    rownames(table2p) = NULL
    rownames(table2r) = NULL
  } else {
    
    
    #####################################################################################################
    #put results of two models in one cell
    # read mplus model
    mplodel = readModels(target = mplodelnames[1])
    #extract std params
    stdout = mplodel$parameters$stdyx.standardized
    # take unstandardized if std is not available
    if (is.null(stdout)) {
      stdout = mplodel$parameters$unstandardized
      # large p value for unstandardized so that it gets marked in final table
      stdout$pval = 1
    }
    #delete unnecessary rows
    stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
    stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
    rownames(stdout) = 1:dim(stdout)[1]
    # for mxe calc
    temp = stdout
    # reorder within, between, pred and extract only variables we are interested in
    stdout = stdout[orders, ]
    if (length(adds1models) > 1) {
      # read additional S1 model
      # read mplus model
      adds1model = readModels(target = adds1models[1])
      #extract std params
      stdouts1 = adds1model$parameters$stdyx.standardized
      # take unstandardized if std is not available
      if (is.null(stdouts1)) {
        stdouts1 = adds1model$parameters$unstandardized
        # large p value for unstandardized so that it gets marked in final table
        stdouts1$pval = 1
        # stdouts1 = stdout
        # stdouts1[,3:ncol(stdouts1)] = NA
      }
      #delete unnecessary rows
      stdouts1 = stdouts1[-grep("esidual", stdouts1$paramHeader) , ]
      stdouts1 = stdouts1[-grep("ntercept", stdouts1$paramHeader) , ]
      rownames(stdouts1) = 1:dim(stdouts1)[1]
      # for mxe calc
      temps1 = stdouts1
      # reorder within, between, pred and extract only variables we are interested in
      stdouts1 = stdouts1[orders1, ]
    }
    # check available N after removal of cases wo variance and or
    nocases = unlist(regmatches(unlist(mplodel$warnings), gregexpr('[0-9]+', unlist(
      mplodel$warnings
    ))))
    nocases = nocases[nchar(nocases) > 1]
    N = mplodel$data_summary$overall$NClusters - length(nocases)
    # create table
    table2 = stdout[, c("paramHeader", "est")]
    names(table2) = c("param", paste0(tablenames[1], " N=", N, "<br> stdyx."))
    if (length(is.na(params)) == 1) {
      #put params
      table2$param = paste(stdout$paramHeader, stdout$param)
    } else{
      table2$param = params
    }
    # create table with boldness when nonoverlapping ci
    if (predtable == T) {
      tablebolds = table2
      tablebolds[, 2] = FALSE
      # bold for non-overlapping PF and Dep conf int
      i = 1
      while (i < dim(stdout)[1]) {
        a1 = stdout$lower_2.5ci[i]
        a2 = stdout$upper_2.5ci[i]
        b1 = stdout$lower_2.5ci[i + 1]
        b2 = stdout$upper_2.5ci[i + 1]
        if (!((a1 <= b2) && (b1 <= a2))) {
          tablebolds[i, 2] = T
          tablebolds[i + 1, 2] = T
        }
        i = i + 2
      }
    } else{
      tablebolds = table2
      tablebolds[, 2:ncol(tablebolds)] = F
    }
    if (logtransparams == T) {
      #old procedure
      #  rawout =  mplodel$parameters$unstandardized
      # rawout = rawout[grep("ntercept",rawout$paramHeader) ,]
      # # exchange m_s, miv or siv params with log-backtransformed values from raw output
      # rawout = rawout[rawout$param %in% c("MIV","SIV","M_S"),]
      # stdout[grep("M_S_LV", stdout$param),"est"] = exp(rawout$est[rawout$param == "M_S"]) / (sqrt(exp(rawout$est[rawout$param == "MIV"])+exp(rawout$est[rawout$param == "M_S"]))*sqrt(exp(rawout$est[rawout$param == "SIV"])+exp(rawout$est[rawout$param == "M_S"])))
      # # if negative loading
      # if(temp < 0){stdout[grep("M_S_LV", stdout$param),"est"] = stdout[grep("M_S_LV", stdout$param),"est"] *-1}
      ###### procedure advised by Ellen Hamaker  #######
      # unique part of residual variances
      loadm =   temp[temp$paramHeader == "M_S_LV.BY", "est"][1]
      loads =  temp[temp$paramHeader == "M_S_LV.BY", "est"][2]
      # covariance between residual variances
      covms = loads * loadm
      # residual variances
      varresidm = temp[temp$paramHeader == "MIV.|", "est"] + loadm ^ 2
      varresids = temp[temp$paramHeader == "SIV.|", "est"] + loads ^ 2
      # correlation = cov / sqrt(varresidm * varresids)
      stdout[grep("M_S_LV", stdout$param), "est"] = covms / sqrt(varresidm * varresids)
      #confidence interval according to Olivoto et al. (2018)
      CI = (abs(stdout[grep("M_S_LV", stdout$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
      stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] - CI
      stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] + CI
      # stdout[grep("SIV.", stdout$paramHeader),"est"] =exp(rawout$est[rawout$param == "SIV"])
      # stdout[grep("SIV.", stdout$paramHeader),"lower_2.5ci"] =exp(rawout$lower_2.5ci[rawout$param == "SIV"])
      # stdout[grep("SIV.", stdout$paramHeader),"upper_2.5ci"] =exp(rawout$upper_2.5ci[rawout$param == "SIV"])
      # if("&pi;m" %in% table2$param){
      #  stdout[grep("MIV.", stdout$paramHeader),"est"] =exp(rawout$est[rawout$param == "MIV"])
      #  stdout[grep("MIV.", stdout$paramHeader),"lower_2.5ci"] =exp(rawout$lower_2.5ci[rawout$param == "MIV"])
      # stdout[grep("MIV.", stdout$paramHeader),"upper_2.5ci"] =exp(rawout$upper_2.5ci[rawout$param == "MIV"])
      # }
      # also for additional s1
      if (length(adds1models) > 1) {
        ###### procedure advised by Ellen Hamaker  #######
        # unique part of residual variances
        loadm =   temps1[temps1$paramHeader == "M_S_LV.BY", "est"][1]
        loads =  temps1[temps1$paramHeader == "M_S_LV.BY", "est"][2]
        # covariance between residual variances
        covms = loads * loadm
        # residual variances
        varresidm = temps1[temps1$paramHeader == "MIV.|", "est"] + loadm ^ 2
        varresids = temps1[temps1$paramHeader == "SIV.|", "est"] + loads ^ 2
        # correlation = cov / sqrt(varresidm * varresids)
        stdouts1[grep("M_S_LV", stdouts1$param), "est"] = covms / sqrt(varresidm * varresids)
        #confidence interval according to Olivoto et al. (2018)
        CI = (abs(stdouts1[grep("M_S_LV", stdouts1$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
        stdouts1[grep("M_S_LV", stdouts1$param), "lower_2.5ci"] = stdouts1[grep("M_S_LV", stdouts1$param), "est"] - CI
        stdouts1[grep("M_S_LV", stdouts1$param), "upper_2.5ci"] = stdouts1[grep("M_S_LV", stdouts1$param), "est"] + CI
      }
    }
   # create noconf boolean
 if(length(noconf) == 1){ if(sum(is.na(noconf)) > 0) {noconf = rep(F,length(orders))} else if(noconf){noconf = rep(T,length(orders))}} 
if (html == T) {
 table2[!noconf, 2] = paste0(
    "<span style=\"",
    ifelse(stdout$pval[!noconf] < .005 & stdouts1$pval[!noconf] < .005, yes = "font-weight: bold;", no = ";"),
    "\" >",
    round(rowMeans(cbind(stdout$est[!noconf],stdouts1$est[!noconf])), digits = 2),
    " <br> [",
    round(rowMeans(cbind(stdout$lower_2.5ci[!noconf],stdouts1$lower_2.5ci[!noconf])), digits = 2),
    " ",
    round(rowMeans(cbind(stdout$upper_2.5ci[!noconf],stdouts1$upper_2.5ci[!noconf])), digits = 2),
    "]</span>"
)}
if (html == F) {
    table2[!noconf, 2] = ifelse(
        stdout$pval[!noconf] < .005 & stdouts1$pval[!noconf] < .005,
        yes = paste0(
            round(rowMeans(cbind(stdout$est[!noconf],stdouts1$est[!noconf])), digits = 2),
            " [",
            round(rowMeans(cbind(stdout$lower_2.5ci[!noconf],stdouts1$lower_2.5ci[!noconf])), digits = 2),
            " ",
            round(rowMeans(cbind(stdout$upper_2.5ci[!noconf],stdouts1$upper_2.5ci[!noconf])), digits = 2),
            "]"
        ),
        no = ""
    )
}
# 
# table2[noconf, 2] = paste0(
#     "<span style=\"",
#     ifelse(stdout$pval < .005, yes = "font-weight: bold;", no = ";"),
#     "\" >",
#     round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2),
#     "</span>"
# )
# if (html == F) {
#     table2[noconf, 2] = ifelse(stdout$pval < .005, yes = paste0(round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2)), no = "")
# }
    
    # if cells without confidence interval then use estimates of the two S-1 models in these cells
    if(sum(noconf) > 0){
    if (html == T) {
      table2[noconf, 2] = paste0(
        "<span style=\"",
        ifelse(stdout$pval[noconf] < .005 & stdouts1$pval[noconf] < .005 & !((pmin(stdout$est[noconf], stdouts1$est[noconf]) <= 0) & (0 <= pmax(stdout$est[noconf], stdouts1$est[noconf]))) , yes = "font-weight: bold;", no = ";"),
        "\" >",
        ifelse(!(pmin(stdout$est[noconf], stdouts1$est[noconf]) <= 0) & (0 <= pmax(stdout$est[noconf], stdouts1$est[noconf])), yes = "<b>", no = ""),
        round(rowMeans(cbind(stdout$est[noconf], stdouts1$est[noconf]), na.rm = T), digits = 2),
        "</b><br> [",
        ifelse(stdout$pval[noconf] < .005, yes = "<b>", no = ""),
        round(stdout$est[noconf], digits = 2),
        " ",
        ifelse(stdouts1$pval[noconf] < .005, yes = "</b><b>", no = "</b>"),
        round(stdouts1$est[noconf], digits = 2),
        "</b>]",
        "</span>"
      )}
    if (html == F) {
      table2[noconf, 2] = ifelse(stdout$pval[noconf] < .005, yes = paste0(round(rowMeans(cbind(stdout$est[noconf], stdouts1$est[noconf]), na.rm = T), digits = 2)), no = "")
    }
    }
    #
    # table2[,2] = paste0("<span style=\"", ifelse(stdout$pval <.005,yes = "font-weight: bold; display:block;width:100%;background-color:yellow;", no = "display:block;width:100%;background-color:yellow;"),   "\" >",round(stdout$est,digits = 2)," <br> [",round(stdout$lower_2.5ci,digits = 2)," ",round(stdout$upper_2.5ci,digits = 2),"]</span>")
    #
    # table2[,2]= paste(round(stdout$est,digits = 2),paste0("<br>[",round(stdout$lower_2.5ci,digits = 2), " ",round(stdout$upper_2.5ci,digits = 2),"]"))
    # make table with significance
    table2p = table2[, 1:2]
    names(table2p)[2] = tablenames[1]
    table2p[, 2] = apply(cbind(stdout$pval, stdouts1$pval), 1, function (x)
      max(x, na.rm = T))
    # make table with regressionweight
    table2r = table2[, 1:2]
    names(table2r)[2] = tablenames[1]
    table2r[, 2] = rowMeans(cbind(stdout$est, stdouts1$est), na.rm = T)
    # make table with all models
    for (i in 2:length(mplodelnames)) {
      # read mplus model
      mplodel = readModels(target = mplodelnames[i])
      #extract std params
      stdout = mplodel$parameters$stdyx.standardized
      # take unstandardized if std is not available
      if (is.null(stdout)) {
        stdout = mplodel$parameters$unstandardized
        # large p value for unstandardized so that it gets marked in final table
        stdout$pval = 1
        stdout$lower_2.5ci = -1
        stdout$upper_2.5ci = 1
      }
      #delete unnecessary rows
      stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
      stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
      rownames(stdout) = 1:dim(stdout)[1]
      # for mxe calc
      temp = stdout
      # reorder within, between, pred and extract only variables we are interested in
      stdout = stdout[orders, ]
      if (length(adds1models) > 1) {
        # read additional S1 model
        # read mplus model
        adds1model = readModels(target = adds1models[i])
        #extract std params
        stdouts1 = adds1model$parameters$stdyx.standardized
        # take unstandardized if std is not available
        if (is.null(stdouts1)) {
          stdouts1 = adds1model$parameters$unstandardized
          # large p value for unstandardized so that it gets marked in final table
          stdouts1$pval = 1
          stdouts1$lower_2.5ci = -1
          stdouts1$upper_2.5ci = 1
          # stdouts1 = stdout
          # stdouts1[,3:ncol(stdouts1)] = NA
        }
        #delete unnecessary rows
        stdouts1 = stdouts1[-grep("esidual", stdouts1$paramHeader) , ]
        stdouts1 = stdouts1[-grep("ntercept", stdouts1$paramHeader) , ]
        rownames(stdouts1) = 1:dim(stdouts1)[1]
        # for mxe calc
        temps1 = stdouts1
        # reorder within, between, pred and extract only variables we are interested in
        stdouts1 = stdouts1[orders1, ]
      }
      # check available N after removal of cases wo variance and or
      nocases = unlist(regmatches(unlist(mplodel$warnings), gregexpr(
        '[0-9]+', unlist(mplodel$warnings)
      )))
      nocases = nocases[nchar(nocases) > 1]
      N = mplodel$data_summary$overall$NClusters - length(nocases)
      if (logtransparams == T) {
        ###### procedure advised by Ellen Hamaker  #######
        # unique part of residual variances
        loadm =   temp[temp$paramHeader == "M_S_LV.BY", "est"][1]
        loads =  temp[temp$paramHeader == "M_S_LV.BY", "est"][2]
        # covariance between residual variances
        covms = loads * loadm
        # residual variances
        varresidm = temp[temp$paramHeader == "MIV.|", "est"] + loadm ^ 2
        varresids = temp[temp$paramHeader == "SIV.|", "est"] + loads ^ 2
        # correlation = cov / sqrt(varresidm * varresids)
        stdout[grep("M_S_LV", stdout$param), "est"] = covms / sqrt(varresidm * varresids)
        #confidence interval according to Olivoto et al. (2018)
        CI = (abs(stdout[grep("M_S_LV", stdout$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
        stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] - CI
        stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = stdout[grep("M_S_LV", stdout$param), "est"] + CI
        # also for additional s1
        if (length(adds1models) > 1) {
          ###### procedure advised by Ellen Hamaker  #######
          # unique part of residual variances
          loadm =   temps1[temps1$paramHeader == "M_S_LV.BY", "est"][1]
          loads =  temps1[temps1$paramHeader == "M_S_LV.BY", "est"][2]
          # covariance between residual variances
          covms = loads * loadm
          # residual variances
          varresidm = temps1[temps1$paramHeader == "MIV.|", "est"] + loadm ^ 2
          varresids = temps1[temps1$paramHeader == "SIV.|", "est"] + loads ^ 2
          # correlation = cov / sqrt(varresidm * varresids)
          stdouts1[grep("M_S_LV", stdouts1$param), "est"] = covms / sqrt(varresidm * varresids)
          #confidence interval according to Olivoto et al. (2018)
          CI = (abs(stdouts1[grep("M_S_LV", stdouts1$param), "est"]) ^ .3) * 2.25152 * (N ^ -0.50089)
          stdouts1[grep("M_S_LV", stdouts1$param), "lower_2.5ci"] = stdouts1[grep("M_S_LV", stdouts1$param), "est"] - CI
          stdouts1[grep("M_S_LV", stdouts1$param), "upper_2.5ci"] = stdouts1[grep("M_S_LV", stdouts1$param), "est"] + CI
        }
      }
      
#if there are cells with confints
  if(sum(!noconf)>0){
      
      if (html == T) {
      #put params
      table2[!noconf, i+1] = paste0(
    "<span style=\"",
    ifelse(stdout$pval[!noconf] < .005 & stdouts1$pval[!noconf] < .005, yes = "font-weight: bold;", no = ";"),
    "\" >",
    round(rowMeans(cbind(stdout$est,stdouts1$est))[!noconf], digits = 2),
    " <br> [",
    round(rowMeans(cbind(stdout$lower_2.5ci,stdouts1$lower_2.5ci))[!noconf], digits = 2),
    " ",
    round(rowMeans(cbind(stdout$upper_2.5ci,stdouts1$upper_2.5ci))[!noconf], digits = 2),
    "]</span>"
)}
if (html == F) {
    table2[!noconf, i +1] = ifelse(
        stdout$pval[!noconf] < .005 & stdouts1$pval[!noconf] < .005,
        yes = paste0(
            round(rowMeans(cbind(stdout$est,stdouts1$est))[!noconf], digits = 2),
            " <br> [",
            round(rowMeans(cbind(stdout$lower_2.5ci,stdouts1$lower_2.5ci))[!noconf], digits = 2),
            " ",
            round(rowMeans(cbind(stdout$upper_2.5ci,stdouts1$upper_2.5ci))[!noconf], digits = 2),
            "]"
        ),
        no = ""
    )
}}
# 
# table2[noconf, i +1] = paste0(
#     "<span style=\"",
#     ifelse(stdout$pval < .005, yes = "font-weight: bold;", no = ";"),
#     "\" >",
#     round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2),
#     "</span>"
# )
# if (html == F) {
#     table2[noconf, 2] = ifelse(stdout$pval < .005, yes = paste0(round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2)), no = "")
# }
      
  # if cells without confint use estimates from the two S-1 models    
   if(sum(noconf) > 0){
     if (html == T) {
      table2[noconf, i +1] = paste0(
        "<span style=\"",
        ifelse(stdout$pval[noconf] < .005 & stdouts1$pval[noconf] < .005 & !((pmin(stdout$est[noconf], stdouts1$est[noconf]) <= 0) & (0 <= pmax(stdout$est[noconf], stdouts1$est[noconf]))) , yes = "font-weight: bold;", no = ";"),
        "\" >",
        ifelse(!(pmin(stdout$est[noconf], stdouts1$est[noconf]) <= 0) & (0 <= pmax(stdout$est[noconf], stdouts1$est[noconf])), yes = "<b>", no = ""),
        round(rowMeans(cbind(stdout$est, stdouts1$est), na.rm = T)[noconf], digits = 2),
        "</b><br> [",
        ifelse(stdout$pval[noconf] < .005, yes = "<b>", no = ""),
        round(stdout$est[noconf], digits = 2),
        " ",
        ifelse(stdouts1$pval[noconf] < .005, yes = "</b><b>", no = "</b>"),
        round(stdouts1$est[noconf], digits = 2),
        "</b>]",
        "</span>"
      )}
    if (html == F) {
      table2[noconf, i +1] = ifelse(stdout$pval[noconf] < .005 & stdouts1$pval[noconf] < .005 & !((pmin(stdout$est[noconf], stdouts1$est[noconf]) <= 0) & (0 <= pmax(stdout$est[noconf], stdouts1$est[noconf]))),
                                    yes = paste0(round(rowMeans(cbind(stdout$est, stdouts1$est), na.rm = T)[noconf], digits = 2)), no = "")
    }}
      names(table2)[i + 1] = paste0(tablenames[i], " N=", N, "<br> stdyx.")
      # create table with boldness when nonoverlapping ci
      if (predtable == T) {
        # bold for non-overlapping PF and Dep conf int
        tablebolds[, i + 1] = FALSE
        d = 1
        while (d < dim(stdout)[1]) {
          a1 = stdout$lower_2.5ci[d]
          a2 = stdout$upper_2.5ci[d]
          b1 = stdout$lower_2.5ci[d + 1]
          b2 = stdout$upper_2.5ci[d + 1]
          if (!((a1 <= b2) && (b1 <= a2))) {
            tablebolds[d, i + 1] = T
            tablebolds[d + 1, i + 1] = T
          }
          d = d + 2
        }
      } else{
        tablebolds[, i + 1] = FALSE
      }
      # put in table with significance
      table2p[, i + 1] = apply(cbind(stdout$pval, stdouts1$pval), 1, function (x)
        max(x, na.rm = T))
      names(table2p)[i + 1] = tablenames[i]
      # put in table with regression
      table2r[, i + 1] = rowMeans(cbind(stdout$est, stdouts1$est), na.rm = T)
      names(table2r)[i + 1] = tablenames[i]
    }
    rownames(table2) = NULL
    rownames(table2p) = NULL
    rownames(table2r) = NULL
  }
  return(list(table2, table2p, table2r, tablebolds))
}

Tables

## within table dscrete situations without mood
table = table2func(mplodelnames = c("ModNoOutS1PFsit_ncnf.out","ModNoOutS1PFsit_nemp.out","ModNoOutS1PFsit_novl.out","ModNoOutS1PFsit_prel.out","ModNoOutS1PFsit_psuc.out","ModNoOutS1PFsit_pwo.out"), tablenames = c("conflict","emptiness or boredom","overwhelming <br> task","relaxation","achievement","good time with someone"),params = c("mm","ee","mxe","me1","em1","miv","eiv"),orders = c(3,5,7,4,6,8,9), html = T,logtransparams = T)

# log back transformation m_s, miv siv
#exp(m_s) / (sqrt(exp(miv)+exp(m_s))*sqrt(exp(siv)+exp(m_s)))

knitr::kable(table[[1]],escape = FALSE, align = "c",format = "html") %>%
      kable_paper() %>% 
  column_spec(2,
              background = spec_color(table[[3]][,2], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(3,
              background = spec_color(table[[3]][,3], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(4,
              background = spec_color(table[[3]][,4], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 

  column_spec(5,
              background = spec_color(table[[3]][,5], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
 
  column_spec(6,
              background = spec_color(table[[3]][,6], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  
   column_spec(7,
              background = spec_color(table[[3]][,7], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 

   add_header_above(c(" ", "Within-Person-Dynamics" = 6))
  



## Table with baseline PF, DEP and events with two averaged S1 mods
table = table2func(mplodelnames = c("ModNoOutS1PFsit_ncnf.out","ModNoOutS1PFsit_nemp.out","ModNoOutS1PFsit_novl.out","ModNoOutS1PFsit_prel.out","ModNoOutS1PFsit_psuc.out","ModNoOutS1PFsit_pwo.out"), adds1models = c("ModNoOutS1DEPsit_ncnf.out","ModNoOutS1DEPsit_nemp.out","ModNoOutS1DEPsit_novl.out","ModNoOutS1DEPsit_prel.out","ModNoOutS1DEPsit_psuc.out","ModNoOutS1DEPsit_pwo.out"), tablenames = c("conflict","emptiness or boredom","overwhelming <br> task","relaxation","achievement","good time with someone"),
params = c("PF~>mm",   "DEP~>mm","PF~>ss",   "DEP~>ss","PF~>mxe",   "DEP~>mxe" , "PF~>me1","DEP~>me1","PF~>em1",  "DEP~>em1", "PF~>miv" , "DEP~>miv","PF~>eiv" , "DEP~>eiv","PF~>&empty;mood" , "DEP~>&empty;mood" ,"PF~>&empty;event" , "DEP~>&empty;event"),orders = c(34,35,36,37,28,29,38,39,40,41,30,31,32,33,42,43,44,45), orders1 = c(43, 44, 45, 46, 37, 38, 47, 48, 49, 50, 39, 40, 41, 42, 51, 52, 53, 54), html = T, predtable = F, noconf = T)

table[[1]]%>% kable(escape = FALSE, align = "c",format = "html") %>%
kable_paper() %>%
column_spec(2,bold = table[[4]][,2],
background = spec_color(table[[3]][,2], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(3,bold = table[[4]][,3],
background = spec_color(table[[3]][,3], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(4,bold = table[[4]][,4],
background = spec_color(table[[3]][,4], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(5,bold = table[[4]][,5],
background = spec_color(table[[3]][,5], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(6,bold = table[[4]][,6],
background = spec_color(table[[3]][,6], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(7,bold = table[[4]][,7],
background = spec_color(table[[3]][,7], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
add_header_above(c(" ", "Prediction of averaged between-person differences in within-person affect-sitation dynamics" = 6))


## Prediction table
table2 = table2func(mplodelnames = c("Model1YncnfS1PCO.out","Model1YnempS1PCO.out","Model1YnovlS1PCO.out","Model1YprelS1PCO.out","Model1YpsucS1PCO.out","Model1YpwoS1PCO.out"), adds1models = c("Model1YncnfS1PCODep.out","Model1YnempS1PCODep.out","Model1YnovlS1PCODep.out","Model1YprelS1PCODep.out","Model1YpsucS1PCODep.out","Model1YpwoS1PCODep.out"), tablenames = c("conflict","emptiness or boredom","overwhelming <br> task","relaxation","achievement","good time with someone"),
params = c("mm~>OPF",   "mm~>ODEP","ss~>OPF",   "ss~>ODEP","mxe~>OPF",   "mxe~>ODEP" , "me1~>OPF","me1~>ODEP","em1~>OPF",  "em1~>ODEP", "miv~>OPF" , "miv~>ODEP","eiv~>OPF" , "eiv~>ODEP","&empty;mood~>OPF" , "&empty;mood~>ODEP" ,"&empty;event~>OPF" , "&empty;event~>ODEP","PF~>OPF *","DEP~>OPF *","PF~>ODEP *","DEP~>ODEP *"),orders = c(66,57,67,58,70,61,68,59,69,60,71,62,72,63,77,75,78,76,73,74,64,65), orders1 = c(75, 66, 76, 67, 79, 70, 77, 68, 78, 69, 80, 71, 81, 72, 86, 84, 87, 85, 82, 83, 73, 74), html = T, predtable = F, noconf = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, T, T, T))


table2[[1]]%>% kable(escape = FALSE, align = "c",format = "html") %>%
kable_paper() %>%
column_spec(2,#bold = table2[[4]][,2],
background = spec_color(table2[[3]][,2], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(3,#bold = table2[[4]][,3],
background = spec_color(table2[[3]][,3], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(4,#bold = table2[[4]][,4],
background = spec_color(table2[[3]][,4], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(5,#bold = table2[[4]][,5],
background = spec_color(table2[[3]][,5], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(6,#bold = table2[[4]][,6],
background = spec_color(table2[[3]][,6], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
column_spec(7,#bold = table2[[4]][,7],
background = spec_color(table2[[3]][,7], scale_from = c(-.95,.95), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
add_header_above(c(" ", "Prediction of between-person differences in PF and depression (averaged over one year) using within-person affect-sitation dynamics" = 6))





## prediction within IN EACH CELL BOTH S1 PARAMS, CIs and correlated factors in SUPPLEMENT

table = table2func(mplodelnames = c("ModNoOutS1sit_ncnf.out","ModNoOutS1DEPsit_ncnf.out","ModNoOutS1sit_nemp.out","ModNoOutS1DEPsit_nemp.out","ModNoOutS1sit_novl.out","ModNoOutS1DEPsit_novl.out","ModNoOutS1sit_prel.out","ModNoOutS1DEPsit_prel.out","ModNoOutS1sit_psuc.out","ModNoOutS1DEPsit_psuc.out","ModNoOutS1sit_pwo.out","ModNoOutS1DEPsit_pwo.out"), tablenames = c("conflict","emptiness or boredom","overwhelming <br> task","relaxation","achievement","good time with someone"),
params = c("PF~>&Psi;me",   "DEP~>&Psi;me" , "PF~>&Phi;ee","DEP~>&Phi;ee","PF~>&Phi;me1",  "DEP~>&Phi;me1", "PF~>&Phi;em1" , "DEP~>&Phi;em1","PF~>&pi;e" , "DEP~>&pi;e" ,"PF~>&empty;SITS" , "DEP~>&empty;SITS"),orders = c(10,11,14,15,16,17,18,19,22,23,26,27), html = F, predtable = T)

table[[1]]%>% kable(escape = FALSE, align = "c",format = "html") %>%
      kable_paper() %>% 
  column_spec(2,bold = table[[4]][,2],
              background = spec_color(table[[3]][,2], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(3,bold = table[[4]][,3],
              background = spec_color(table[[3]][,3], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(4,bold = table[[4]][,4],
              background = spec_color(table[[3]][,4], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(5,bold = table[[4]][,5],
              background = spec_color(table[[3]][,5], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(6,bold = table[[4]][,6],
              background = spec_color(table[[3]][,6], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(7,bold = table[[4]][,7],
              background = spec_color(table[[3]][,7], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
   add_header_above(c(" ", "Prediction of averaged between-person differences in within-person affect-sitation dynamics" = 6))  %>% save_kable("table3.pdf")




## prediction within only correlated factors

table = table2func(mplodelnames = c("ModNoOutCORsit_ncnf.out","ModNoOutCORsit_nemp.out","ModNoOutCORsit_novl.out","ModNoOutCORsit_prel.out","ModNoOutCORsit_psuc.out","ModNoOutCORsit_pwo.out"), tablenames = c("conflict","emptiness or boredom","overwhelming <br> task","relaxation","achievement","good time with someone"),
params = c("PF~>&Psi;me",   "DEP~>&Psi;me" , "PF~>&Phi;ee","DEP~>&Phi;ee","PF~>&Phi;me1",  "DEP~>&Phi;me1", "PF~>&Phi;em1" , "DEP~>&Phi;em1","PF~>&pi;e" , "DEP~>&pi;e" ,"PF~>&empty;SITS" , "DEP~>&empty;SITS"), orders = c(10,11,14,15,16,17,18,19,22,23,26,27), html = F, predtable = T)

rownames(table[[1]]) = c("PF~>mxe",   "DEP~>mxe" , "PF~>ee","DEP~>ee","PF~>me1",  "DEP~>me1", "PF~>em1" , "DEP~>em1","PF~>&eiv" , "DEP~>eiv" ,"PF~>&empty;e" , "DEP~>&empty;e")

table[[1]]%>% kable(escape = FALSE, align = "c",format = "html") %>%
      kable_paper() %>% 
  column_spec(2,bold = table[[4]][,2],
              background = spec_color(table[[3]][,2], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(3,bold = table[[4]][,3],
              background = spec_color(table[[3]][,3], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(4,bold = table[[4]][,4],
              background = spec_color(table[[3]][,4], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(5,bold = table[[4]][,5],
              background = spec_color(table[[3]][,5], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(6,bold = table[[4]][,6],
              background = spec_color(table[[3]][,6], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
  column_spec(7,bold = table[[4]][,7],
              background = spec_color(table[[3]][,7], scale_from = c(-.6,.6), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>% 
   add_header_above(c(" ", "Prediction of averaged between-person differences in within-person affect-sitation dynamics" = 6))  %>% save_kable("table3.pdf")




# %>% 
#     DT::datatable(
#     extensions = 'Buttons',
#     options = list(dom = 'Bfrtip',
#                    buttons = list('excel', "csv")))
# 




  

# Function to apply bold formatting based on logical table
bold_if_true <- function(value, condition) {
  if (condition) {
    return(cell_spec(value, bold = TRUE))
  } else {
    return(value)
  }
}

# Apply the function to the data table
formatted_table <- table2[[1]]
for (col in 2: length(names(table2[[1]]))) {
  formatted_table[[col]] <- mapply(bold_if_true, formatted_table[[col]], table2[[2]][[col]] < .005)
}

# Generate the kable object with kableExtra
kable(formatted_table,escape = FALSE) %>%
      kable_classic_2()