---
title: "R Notebook"
output: html_notebook
---

# extract correlation and covariance matrices from mplus output
```{r}
covrmatfunc = function(mplodel){
  
  moreout = MplusAutomation::readModels(target = mplodel)
  moreout = moreout$output[grep("ESTIMATES DERIVED FROM THE MODEL FOR BETWEEN", moreout$output) : length(moreout$output)]
  cormat = moreout[grep("ESTIMATED CORRELATION MATRIX FOR THE LATENT VARIABLES", moreout) : length(moreout)]
  
  covmat = moreout[grep("ESTIMATED COVARIANCE MATRIX FOR THE LATENT VARIABLES", moreout) : grep("ESTIMATED CORRELATION MATRIX FOR THE LATENT VARIABLES", moreout) -1]
  
  
  extract_covrmat <- function(txt) {
    # txt: character vector, one element per line of the printed output
    
    # helper to trim whitespace
    trim <- function(x) gsub("^\\s+|\\s+$", "", x)
    
    block_cols <- list()
    rows_list  <- list()
    
    i <- 1
    while(i <= length(txt)) {
      # look for a header line: non‐underscore names separated by 2+ spaces,
      # followed immediately by an "____" line
      if (grepl("^\\s*\\S+", txt[i]) && grepl("\\s{2,}", txt[i]) &&
          i + 1 <= length(txt) && grepl("^\\s*_{2,}", txt[i+1])) {
        
        # parse column names from the header
        this_cols <- strsplit(trim(txt[i]), "\\s{2,}")[[1]]
        block_cols[[length(block_cols) + 1]] <- this_cols
        
        # now read the numeric rows until a blank line or end
        j <- i + 2
        while(j <= length(txt) && nzchar(txt[j])) {
          line <- txt[j]
          parts <- strsplit(trim(line), "\\s+")[[1]]
          var   <- parts[1]
          vals  <- as.numeric(parts[-1])
          rows_list[[length(rows_list) + 1]] <-
            list(var = var, vals = vals, cols = this_cols)
          j <- j + 1
        }
        i <- j
      } else {
        i <- i + 1
      }
    }
    
    # flatten column order (preserve block order), drop cols and rows with "&"
    all_cols <- unique(unlist(block_cols))
    all_cols <- all_cols[-grep("&", all_cols)]
    
    # initialize empty matrix
    M <- matrix(NA_real_,
                nrow = length(all_cols),
                ncol = length(all_cols),
                dimnames = list(all_cols, all_cols))
    diag(M) <- 1
    
    # fill in off‐diagonals
    for (rec in rows_list) {
      var  <- rec$var
      if (grepl("&", var)) next
      vals <- rec$vals
      cols <- rec$cols
      # skip all‐zero rows
      if (all(vals == 0)) next
      
      for (k in seq_along(cols)) {
        colk <- cols[k]
        if (grepl("&", colk)) next
        M[var, colk] <- vals[k]
        M[colk, var] <- vals[k]
      }
    }
    
    # return as data.frame
    as.data.frame(M, row.names = all_cols)[-1,-1]
  }

  
 list(cormat = extract_covrmat(cormat), covmat = extract_covrmat(covmat))
  
}
```


# calculate (squared semipartial) correlations and put in table
```{r}
sitnames = c("conflict","boredeom","overwhelm","relax","achievem","positive_social")

pfmods = c("Model1YncnfS1PCO.out","Model1YnempS1PCO.out","Model1YnovlS1PCO.out","Model1YprelS1PCO.out","Model1YpsucS1PCO.out","Model1YpwoS1PCO.out")

#create list with correlation and covariance matrices from PF S1 models
pfcovrs = list()
for (i in 1 : length(pfmods)){
  pfcovrs[[i]] = covrmatfunc(pfmods[i])
}

names(pfcovrs) = sitnames

depmods = c("Model1YncnfS1PCODep.out","Model1YnempS1PCODep.out","Model1YnovlS1PCODep.out","Model1YprelS1PCODep.out","Model1YpsucS1PCODep.out","Model1YpwoS1PCODep.out")

# depmods = c("Model1YncnfS1PCODep.out","Model1YnempS1PCODep.out","Model1YnovlS1PCODep.out","Model1YprelS1PCODep.out","Model1YpwoS1PCODep.out","Model1YpwoS1PCODep.out")


#create list with correlation and covariance matrices from DEP S1 models
depcovrs = list()
for (i in 1 : length(depmods)){
  depcovrs[[i]] = covrmatfunc(depmods[i])
}

names(depcovrs) = sitnames



# create list with squared semipartial correlations


## FUNC as ref
pfsqspcors = list()
for (i in 1 : length(pfcovrs)){
  # vars for prediction
predvarnames = names(pfcovrs[[i]][[1]])[c(-3,-4)]
  
##correlations OFUNC
cor_mat.OFUNC = pfcovrs[[i]]$cormat[c("OFUNC",predvarnames),
                         c("OFUNC",predvarnames)]
##correlations ODEP
cor_mat.ODEP = pfcovrs[[i]]$cormat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]


### calculate squared semipartial correlations in predicting ODEP and OFUNC
# using covariances
## OFUNC
cov_mat.OFUNC = pfcovrs[[i]]$covmat[c("OFUNC",predvarnames),
                              c("OFUNC",predvarnames)]



sqsemipcor.OFUNC <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.OFUNC),cov = as.matrix(cov_mat.OFUNC)))^2,digits = 4)


## ODEP
cov_mat.ODEP = pfcovrs[[i]]$covmat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]

sqsemipcor.ODEP <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.ODEP),cov = as.matrix(cov_mat.ODEP)))^2,digits = 4)

pfsqspcors[[i]]  = list(OFUNC = sqsemipcor.OFUNC, ODEP = sqsemipcor.ODEP)
}

names(pfsqspcors) = sitnames




## DEP as ref
depsqspcors = list()
for (i in 1 : length(depcovrs)){
  # vars for prediction
predvarnames = names(depcovrs[[i]][[1]])[c(-3,-4)]
  
##correlations OFUNC
cor_mat.OFUNC = depcovrs[[i]]$cormat[c("OFUNC",predvarnames),
                         c("OFUNC",predvarnames)]
##correlations ODEP
cor_mat.ODEP = depcovrs[[i]]$cormat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]


### calculate squared semipartial correlations in predicting ODEP and OFUNC
# using covariances
## OFUNC
cov_mat.OFUNC = depcovrs[[i]]$covmat[c("OFUNC",predvarnames),
                              c("OFUNC",predvarnames)]



sqsemipcor.OFUNC <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.OFUNC),cov = as.matrix(cov_mat.OFUNC)))^2,digits = 4)


## ODEP
cov_mat.ODEP = depcovrs[[i]]$covmat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]

sqsemipcor.ODEP <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.ODEP),cov = as.matrix(cov_mat.ODEP)))^2,digits = 4)

depsqspcors[[i]]  = list(OFUNC = sqsemipcor.OFUNC, ODEP = sqsemipcor.ODEP)
}

names(depsqspcors) = sitnames



```


## create table with cors and squared semipartial cors, for PF and Dep with both s1 mods
```{r}
# table data frame
rownamess = c(
  "mm~>PF", "mm~>DEP",
  "ss~>PF", "ss~>DEP",
  "me1~>PF", "me1~>DEP",
  "em1~>PF", "em1~>DEP",
  "miv~>PF", "miv~>DEP",
  "eiv~>PF", "eiv~>DEP",
  "mxe~>PF", "mxe~>DEP",
  "&empty;mood~>PF", "&empty;mood~>DEP",
  "&empty;event~>PF", "&empty;event~>DEP",
  "PF~>OPF", "DEP~>OPF",
  "PF~>ODEP", "DEP~>ODEP"
)

colnames = c("conflict <br> N=930","emptiness or boredom <br> N = 1103","overwhelming task <br> N = 1155","relaxation <br> N = 992","achievement <br> N = 963 ","good time with someone <br> N = 1197")

Ns = c(930,1103,1155,992,963,1197)

table4 = data.frame(matrix(,nrow = length(rownamess), ncol = length(colnames)),row.names = rownamess)
names(table4) = colnames

rmzero <- function(x, digits = NULL) {
  s <- if (is.null(digits)) {
    as.character(x)
  } else {
    format(x, nsmall = digits)
  }
  sub("^(-?)0\\.", "\\1.", s)
}
table4r = table4
table4sqsp = table4

# for(x in 1: nrow(table4)){
  for(y in 1: ncol(table4)){
    # dynamics and OFUNC 
    table4[grep("~>PF", rownamess),y] = 
      paste0(
          "<span style=\"",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[5:13,3])$p < .005, yes = "font-weight: bold;", no = ";"),
          "\" >",
          rmzero(round(pfcovrs[[y]]$cormat[5:13,3], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$OFUNC[1,4:12], digits = 3)),
          "]</span>"
        )
    table4r[grep("~>PF", rownamess),y] = pfcovrs[[y]]$cormat[5:13,3]
    table4sqsp[grep("~>PF", rownamess),y] = t(pfsqspcors[[y]]$OFUNC[1,4:12])
    
    #dynamics and ODEP
    table4[grep("~>DEP", rownamess),y] = 
      paste0(
          "<span style=\"",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[5:13,4])$p < .005, yes = "font-weight: bold;", no = ";"),
          "\" >",
          rmzero(round(pfcovrs[[y]]$cormat[5:13,4], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$ODEP[1,4:12], digits = 3)),
          "]</span>"
        )
    table4r[grep("~>DEP", rownamess),y] = pfcovrs[[y]]$cormat[5:13,4]
    table4sqsp[grep("~>DEP", rownamess),y] = t(pfsqspcors[[y]]$ODEP[1,4:12])

    
    # PF OPF DEP ODEP
    table4[grep("~>OPF", rownamess),y] = 
      paste0(
          '<span style=";" >',
          #average from two S1 mods
          # ifelse(psych::r.test(Ns[y],rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,3],depcovrs[[y]]$cormat[1:2,3])))$p < .005, yes = "<b>", no = ""),
          # rmzero(round(rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,3],depcovrs[[y]]$cormat[1:2,3])), digits = 2)),
          # "[",
          # # also put squared semipartial correlations
          # rmzero(round(colMeans(rbind(pfsqspcors[[y]]$OFUNC[1,2:3],depsqspcors[[y]]$OFUNC[1,2:3])), digits = 3)),
          # "]</b> <br>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1:2,3])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfcovrs[[y]]$cormat[1:2,3], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$OFUNC[1,2:3], digits = 3)),
          "] </b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1:2,3])$p < .005, yes = "<b>", no = ""),
          rmzero(round(depcovrs[[y]]$cormat[1:2,3], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(depsqspcors[[y]]$OFUNC[1,2:3], digits = 3)),
          "]</b>",
          "</span>"
        )
    table4r[grep("~>OPF", rownamess),y] = rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,3],depcovrs[[y]]$cormat[1:2,3]))
    table4sqsp[grep("~>OPF", rownamess),y] = colMeans(rbind(pfsqspcors[[y]]$OFUNC[1,2:3],depsqspcors[[y]]$OFUNC[1,2:3]))

    
     # PF OPF DEP ODEP
    table4[grep("~>ODEP", rownamess),y] = 
      paste0(
          '<span style=";" >',
          # #average from two S1 mods
          # ifelse(psych::r.test(Ns[y],rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,4],depcovrs[[y]]$cormat[1:2,4])))$p < .005, yes = "<b>", no = ""),
          # rmzero(round(rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,4],depcovrs[[y]]$cormat[1:2,4])), digits = 2)),
          # "[",
          # # also put squared semipartial correlations
          # rmzero(round(colMeans(rbind(pfsqspcors[[y]]$ODEP[1,2:3],depsqspcors[[y]]$ODEP[1,2:3])), digits = 3)),
          # "]</b> <br>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1:2,4])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfcovrs[[y]]$cormat[1:2,4], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$ODEP[1,2:3], digits = 3)),
          "] </b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1:2,4])$p < .005, yes = "<b>", no = ""),
          rmzero(round(depcovrs[[y]]$cormat[1:2,4], digits = 2)),
          "[",
          # also put squared semipartial correlations
          rmzero(round(depsqspcors[[y]]$ODEP[1,2:3], digits = 3)),
          "]</b>",
          "</span>"
        )
    table4r[grep("~>ODEP", rownamess),y] = rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,4],depcovrs[[y]]$cormat[1:2,4]))
    table4sqsp[grep("~>ODEP", rownamess),y] = colMeans(rbind(pfsqspcors[[y]]$ODEP[1,2:3],depsqspcors[[y]]$ODEP[1,2:3]))

    
  }



```


##  Create table with cors and squared semipartial cors, for PF and Dep with both s1 mods
```{r}
library(tidyverse)
library(kableExtra)

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


```




# Rest
```{r}






covr1y = covrmatfunc("Model1YncnfS1PCO.out")

vars = names(covr1y[[1]])

# without outcome
predvarnames = vars[c(-3,-4)]

##correlations OFUNC
cor_mat.OFUNC = covr1y$cormat[c("OFUNC",predvarnames),
                         c("OFUNC",predvarnames)]

##correlations ODEP
cor_mat.ODEP = covr1y$cormat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]



#### Regression model based on correlatoin matrix
library(metafor)
# # with cov
# regpsfs1 = matreg("psf_av_fu", x=predvarnames, cov_mat.psf.s1, cov = T, n=27173)

#with cor
regOPF = matreg("OFUNC", x=predvarnames, cor_mat.OFUNC, n=930, digits = 2)
regODEP = matreg("ODEP", x=predvarnames, cor_mat.ODEP, n=930, digits = 2)



### calculate squared semipartial correlations in predicting ODEP and OFUNC
# using covariances
## OFUNC
cov_mat.OFUNC = covr1y$covmat[c("OFUNC",predvarnames),
                              c("OFUNC",predvarnames)]



sqsemipcor.OFUNC <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.OFUNC),cov = as.matrix(cov_mat.OFUNC)))^2,digits = 4)


## ODEP
cov_mat.ODEP = covr1y$covmat[c("ODEP",predvarnames),
                             c("ODEP",predvarnames)]

sqsemipcor.ODEP <- round(data.frame(correlation::cor_to_spcor(cor = as.matrix(cor_mat.ODEP),cov = as.matrix(cov_mat.ODEP)))^2,digits = 4)

```


