
---
title: "PF DSEM"
output: html_notebook
author: André Kerber
---

Moodpath was a depression screening and psychoeducation app that, in addition to depression symptom tracking, included items on PF (OPD-SQS). In the following, data from moodpath users is preprocessed and analyzed for the interaction of PF, depression symptoms, mood and daily assessed situations.




# Data Preprocessing
## Libraries
```{r}
library(data.table)
library(psych)
library(tidyverse)
library(dplyr)
library(gridExtra)
library(parallel)
library(MplusAutomation)
library(lavaan)
library(kableExtra)
library(openxlsx)
library(ggpubr)

```

## Load data and keep users with sufficient EMA and PF assessments
```{r}
data =  fread("DATA/sven-dump.csv")

#### TIME ####
data[,assessment_time := as.POSIXct(assessment_time)]
# keep only data after addition of PF items
data = data[assessment_time > "2018-11-12"]


# delete rows were users just did smiley scale
#data[which(sapply(data[,mget(questvec)], function(x) all(is.na(x)))) := NULL ,]

questvec = names(data)[grep("question", names(data))]

#availableD =rowSums(is.na(data[,mget(questvec)]))

questlengt = length(questvec)
data[rowSums(is.na(data[,mget(questvec)])) != questlengt ]


#### ASSESSMENTS
# variable with days, day 0 = first assessment
data[,totalday := as.numeric(assessment_time - min(assessment_time))/86400, by = .(identity_id) ]
# variable for the amount of days used (maximum of totalday) 
data[, days.used := max(totalday),  by = .(identity_id) ]
# create variable amount of assessments 
data[, N_ass := length(totalday),  by = .(identity_id) ]
 
# remove users with less than 3 assessments
data = data[N_ass > 2]

# keep only participants that have situation assessments on at least 80% of all days within one month (30 days)  
data[totalday < 31, sit.adj.30 := (sum(situations != "") / length(totalday)),  by = .(identity_id)  ]
ids.30 = unique(data[sit.adj.30 > .8, identity_id]) 
data = data[identity_id %in% ids.30]


gc()


# Available PF data (12 OPD-SFK items in 3 parcels. Highest loading items on the three subfactors are presented first, and if user continues to use the app, additional items from that factor are also presented )
days = 31
data[,func1 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_138","question_193", "question_194", "question_195")]
data[,func2 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_139", "question_199", "question_200", "question_201")]
data[,func3 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_140","question_196","question_197","question_198")]

ids.PF.30 = unique(data[sit.adj.30 > .8,rowMeans(.SD), by = .(identity_id),.SDcols = c("func1","func2","func3")][!is.na(V1),identity_id])

data= data[identity_id %in% ids.PF.30]
save(data, file="DATA/dataids30.rda")

gc()


```

## Load more longitudinal data
```{r}

## Load data for all remaining participants after change of situations and emotions assessment
data.FU <- bind_rows(
  fread("DATA/FU DATA/Data_FU_IDs1.csv"),
  fread("DATA/FU DATA/Data_FU_IDs2.csv"),
  fread("DATA/FU DATA/Data_FU_IDs3.csv"),
  fread("DATA/FU DATA/Data_FU_IDs4.csv"),
  fread("DATA/FU DATA/Data_FU_IDs5.csv"),
  fread("DATA/FU DATA/Data_FU_IDs6.csv"),
  fread("DATA/FU DATA/Data_FU_IDs7.csv"),
  fread("DATA/FU DATA/Data_FU_IDs8.csv"),
  fread("DATA/FU DATA/Data_FU_IDs9.csv")
)



names(data.FU)[3] = "assessment_time"
data.FU = data.FU[,-c(5,6,7)]

# load more FU data from persons not included in previous (but without happiness score because of change in metabase)
data.FU2 <- bind_rows(
  fread("DATA/FU DATA/Data_FU_noch1.csv"),
  fread("DATA/FU DATA/Data_FU_noch2.csv"),
  fread("DATA/FU DATA/Data_FU_noch3.csv"),
  fread("DATA/FU DATA/Data_FU_noch4.csv")
)
names(data.FU2)[3] = "assessment_time"

#combine
data.FU = bind_rows(data.FU, data.FU2)
rm(data.FU2)

# only IDs with func
data.FU= data.FU[identity_id %in% ids.PF.30]

#delete other question answers
data.FU = data.FU[data.FU$question_id %in% c(1:14,16,46,47,138:140,193:201),]

data.FU[,assessment_time := as.POSIXct(assessment_time)]
# variable with days, day 0 = first assessment
data.FU[,totalday := as.numeric(round((assessment_time - min(assessment_time)) / 86400,0)), by = .(identity_id) ]
# variable for the amount of days used (maximum of totalday) 
data.FU[, days.used := max(totalday),  by = .(identity_id) ]
#days_used <- data.FU[, .N, by = days.used]
# create variable for fraction of assessments at which situations were reported
data.FU[, N_ass := length(totalday),  by = .(identity_id) ]

#long to wide
data.FU = spread(data.FU, question_id, intensity)
names(data.FU) = c("identity_id","assessment_time","totalday","days.used","N_ass","question_001","question_002", "question_003", "question_004", "question_005", "question_006",
                   "question_007","question_008","question_009","question_010", "question_011", "question_012","question_013", "question_014", "question_016","question_046","question_047",
                   "question_138","question_139","question_140","question_193", "question_194", "question_195","question_196","question_197","question_198",
                   "question_199", "question_200", "question_201")


data.FU = data.FU[order(data.FU$identity_id,data.FU$assessment_time),]


data.FU[,c("situations","emotions")] =  "text"
data.FU[,c("situations","emotions")] =  lapply(data.FU[,c("situations","emotions")] ,function(x) as.character(x))
data.FU[,c("situations","emotions")] =  NA
data.FU$happiness_score = NA
data.FU$func1 = NA
data.FU$func2 = NA
data.FU$func3 = NA



data.FU = data.FU[,c("identity_id","assessment_time","totalday","days.used","N_ass","situations","emotions","happiness_score","question_001","question_002", "question_003", "question_004", "question_005", "question_006",
                     "question_007","question_008","question_009","question_010", "question_011", "question_012","question_013", "question_014", "question_016","question_046","question_047",
                     "question_138","question_139","question_140","question_193", "question_194", "question_195","question_196","question_197","question_198",
                     "question_199", "question_200", "question_201","func1","func2","func3")]

```


## Merge 
```{r}
# keep only relevant cols
cols = c("identity_id","assessment_time","totalday", "days.used", "N_ass","situations","emotions","happiness_score",
         names(data)[grep("question",names(data))],"func1","func2","func3")
data = data[,mget(cols)]


##merge data
data= rbindlist(list(data, data.FU), fill = TRUE)  
# test = merge(data, data.FU, by = c("identity_id", "assessment_time"), all.x = T)
# test = data %>% right_join(data.FU, by=c("identity_id", "assessment_time"))
#data = bind_rows(data, data.FU)

data= data[order(data$identity_id,data$assessment_time),]

## days used again
data[,totalday := as.numeric(assessment_time - min(assessment_time))/86400, by = .(identity_id) ]
# variable for the amount of days used (maximum of totalday) 
data[, days.used := max(totalday),  by = .(identity_id) ]

 data[, N_ass := length(totalday),  by = .(identity_id) ]

#  
#  # keep only participants that have situation assessments on at least 80% of all days within one month (30 days)
# data[totalday < 31, sit.adj.30 := (sum(situations != "") / length(totalday)),  by = .(identity_id)  ]
# ids.30 = unique(data[sit.adj.30 > .8, identity_id])
# #data = data[identity_id %in% ids.30]

rm(data.FU)
gc()


```



## Situation and Mood assessments
```{r}

data[, sits := situations ]

# Create situation variables where 0 = not chosen at assessment, 1 = chosen at assessment
data[, sit_pwo := 0 ]
data[situations %like% "Schönes mit Anderen" | situations %like% "Good time with someone", sit_pwo := 1 ]
#data[situations == "", sit_pwo := NA ]
data[sit_pwo == 1 , sits := "pos_with_others" ]

data[ , sit_psuc := 0 ]
data[situations %like% "Erfolgserlebnis" | situations %like% "chievement", sit_psuc :=1]
#data[situations == "", sit_psuc := NA]
data[sit_psuc == 1 , sits := "pos_success" ]

data[ , sit_prel := 0 ]
data[situations %like% "Entspannung" | situations %like% "Relaxation", sit_prel :=1]
#data[situations == "", sit_prel := NA]
data[sit_prel == 1 , sits := "relaxation" ]

data[ , sit_ncnf := 0 ]
data[situations %like% "Konflikt mit Anderen" | situations %like% "onflict", sit_ncnf :=1]
#data[situations == "", sit_ncnf := NA]
data[sit_ncnf == 1 , sits := "conflict"]

data[ , sit_novl := 0 ]
data[situations %like% "Überlastung" | situations %like% "Overwhelming task" | situations %like% "workload", sit_novl :=1]
#data[situations == "", sit_novl := NA]
data[sit_novl == 1 , sits := "overload"]

data[ , sit_nemp :=0]
data[situations %like% "Leere oder Langeweile" | situations %like% "Emptiness or boredom", sit_nemp :=1]
#data[situations == "", sit_nemp := NA]
data[sit_nemp == 1 , sits := "emptiness_or_boredom"]

# create variable for positive/negative situation reported at assessment
data[, sit_pos := sit_pwo + sit_psuc + sit_prel]
#data$sit.pos_binary <- recode(data$sit.pos_all, "2=1; 3=1")

data[, sit_neg := sit_ncnf + sit_novl + sit_nemp]
#data$sit.neg_binary <- recode(data$sit.neg_all, "2=1; 3=1")


# create mood and situation scores that are centered around 0
#data$mood = data$happiness_score -2
# data$sits = data$sit_pos - data$sit_neg
#data$sit_pos = data$sit_pos - 1.5
#data$sit_neg = data$sit_neg - 1.5

```


## Average PF and Depression Items

```{r}

# create single PF items for FIML 
days = 31
data[,func1 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_138")]
data[,func2 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_193")]
data[,func3 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_194")]
data[,func4 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_195")]
data[,func5 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_139")]
data[,func6 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_199")]
data[,func7 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_200")]
data[,func8 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_201")]
data[,func9 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_140")]
data[,func10 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_196")]
data[,func11 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_197")]
data[,func12 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_198")]



```



## Depression core symptoms over 14 days for as long as data is available
```{r}
## function that checks if icd depression core symptoms are available and then calculates average score from all dep questions
depfunction = function(data, days){
  if((!is.na(mean(data[totalday >=days & totalday < days+15,question_001],na.rm = T)) | 
      !is.na(mean(data[totalday >=days & totalday < days+15,question_002],na.rm = T))) &
     (!is.na(mean(data[totalday >=days & totalday < days+15,question_003],na.rm = T)) | 
      !is.na(mean(data[totalday >=days & totalday < days+15,question_004],na.rm = T)) |
      !is.na(mean(data[totalday >=days & totalday < days+15,question_005],na.rm = T)) | 
      !is.na(mean(data[totalday >=days & totalday < days+15,question_006],na.rm = T)))){
    dep = mean(rowMeans(data[totalday >days & totalday < days+14, c("question_001","question_002", "question_003", "question_004", "question_005", "question_006",
                                                                    "question_007","question_008","question_009","question_010", "question_011", "question_012","question_013", "question_014", "question_016","question_046","question_047")],na.rm = T),na.rm = T)
    
  }else{dep =NA_real_}
  return(dep)
}

# create baseline depression score
data = data[,Dep_D1 := depfunction(.SD,1), by = .(identity_id)]

# weeks
Ddays = seq(30, max(data$days.used), by = 14)


# initializew numeric columns
data[,paste0("Dep_D",Ddays)] = NA_integer_

for(col in paste0("Dep_D",Ddays)){
    set(data, j = col, value = as.numeric(data[[col]]))
}

# 
# # Function to parallelize the calculation over different time periods
# calculate_depression_scores <- function(data, time_periods) {
#   library(future.apply)
# options(future.globals.maxSize= 22000*1024^2)
#   # Set up parallel processing
#   plan(multisession, workers = availableCores() - 1) # Use all but one core
#   
#   # Apply the depfunction across time periods in parallel
#   results <- future_lapply(time_periods, function(days) {
#     data[, paste0("Dep_D", days) := depfunction(.SD, days), by = .(identity_id)]
#     return(data)
#   })
#   
#   # Combine results back into a single data.table
#   result_data <- rbindlist(results, use.names = TRUE, fill = TRUE)
#   return(result_data)
# }
# 
# 
#   results <- future_lapply(time_periods, function(days) {
#     test[, paste0("Dep_D", days) := depfunction(.SD, days), by = .(identity_id)]
#     return(test)
#   })
# 
# 
# # Parallelize the computation for the rest of the time periods
# test = data[1:50000]
# 
# test <- calculate_depression_scores(test, seq(30, max(data$days.used), by = 14))


#then every two weeks beginning after first month (inefficient loop)
for (days in seq(30,max(data$days.used), by = 14)){
  data = data[,paste0("Dep_D",days) := depfunction(.SD,days), by = .(identity_id)]
}




# how much available data 

availableD= NA
for (i in 1:length(Ddays)){
  availableD[i]=sum(data[!duplicated(identity_id),!is.na(.SD), .SDcols = paste("Dep_D",as.numeric(Ddays[i]), sep="")])
}
availableD = rbind(paste("Dep_D",as.numeric(Ddays), sep=""), availableD)

deps = names(data)[grep("Dep_D", names(data))]

data$dep1y = F
data$dep1y[rowSums(!is.na(data[,mget(c("Dep_D310",
"Dep_D324",
"Dep_D338",
"Dep_D352",
"Dep_D366",
"Dep_D380"))]))>0] = T
data$dep2m = F
data$dep2m[rowSums(!is.na(data[,mget(c("Dep_D44",
"Dep_D58",
"Dep_D72"))]))>0] = T

# create variable with average depression over one year
data$DepAV1Y[data$dep1y == T] = scale(rowMeans(as.data.frame(data)[data$dep1y == T,paste0("Dep_D", seq(30,380, by = 14))],na.rm = T), scale = FALSE)

# create variable with average depression over 2 months
data$DepAV2M[data$dep2m == T] = scale(rowMeans(as.data.frame(data)[data$dep2m == T,paste0("Dep_D", seq(30,72, by = 14))],na.rm = T), scale = FALSE)


#make variable with available Depression assessments per id
data$availableD = rowSums(!is.na(data[,mget(paste0("Dep_D",Ddays))]))
availableD.t = data[!duplicated(identity_id),availableD]


### daten saven
save(data, file="DATA/data.dyn.P.rda")

```


# Results
## prepare descriptives
```{r}
questions = c("question_001" , "question_002" , "question_046" , "question_047","question_003" , "question_004" , "question_005" , "question_006" , "question_007" , "question_008","question_009" , "question_010" , "question_016","question_011" , "question_012" , "question_013" , "question_014", "question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198" )

DEPnames = c("DEP_depr_mood_1","DEP_hopelessness_1","DEP_hopelessness_2","DEP_hopelessness_3","DEP_dimin_interest_1","DEP_dimin_interest_2","DEP_reduced_energy_1", "reduced_energy_2","DEP_dimin_selfconf_1",
"dimin_selfconf_2","DEP_feelings_guilt", "DEP_feelings_worthlessn", "DEP_sucidal_ideation","DEP_concentration_1","DEP_concentration_2","DEP_insomnia","DEP_loss_appetite")
PFnames = c("PF_relationshp_mod_1","PF_relationshp_mod_2","PF_relationshp_mod_3","PF_relationshp_mod_4","PF_self_1","PF_self_2","PF_self_3","PF_self_4", "PF_contact_1","PF_contact_2","PF_contact_3","PF_contact_4" )


data.descr = data[, c("identity_id","totalday", c("question_001" , "question_002" , "question_046" , "question_047","question_003" , "question_004" , "question_005" , "question_006" , "question_007" , "question_008","question_009" , "question_010" , "question_016","question_011" , "question_012" , "question_013" , "question_014", "question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198"))]

names(data.descr) = c("identity_id","totalday","DEP_depr_mood_1","DEP_hopelessness_1","DEP_hopelessness_2","DEP_hopelessness_3","DEP_dimin_interest_1","DEP_dimin_interest_2","DEP_reduced_energy_1", "reduced_energy_2","DEP_dimin_selfconf_1",
"dimin_selfconf_2","DEP_feelings_guilt", "DEP_feelings_worthlessn", "DEP_sucidal_ideation","DEP_concentration_1","DEP_concentration_2","DEP_insomnia","DEP_loss_appetite","PF_relationshp_mod_1","PF_relationshp_mod_2","PF_relationshp_mod_3","PF_relationshp_mod_4","PF_self_1","PF_self_2","PF_self_3","PF_self_4", "PF_contact_1","PF_contact_2","PF_contact_3","PF_contact_4" )




## create table with frequency of assessment per ID for every question
freq.table = data.frame(id = unique(data.descr$identity_id))
freq.table[,c(DEPnames,PFnames)] = 0

#test = data.frame(table(data[id ==freq.table$id[110],answers_question_id]))


pb = txtProgressBar(min = 0, max = dim(freq.table)[1], initial = 0)
for(i in 1:dim(freq.table)[1]){
  temp = data.descr[identity_id ==freq.table$id[i],]

setTxtProgressBar(pb,i)
   if(length(temp) > 0){
      if(dim(temp)[1] >0){
      for (j in 1:length(DEPnames)){
    freq.table[i,DEPnames[j]] = sum(!is.na(temp[totalday <14,get(DEPnames[j])]))
      }
        for (j in 1:length(PFnames)){
    freq.table[i,PFnames[j]] = sum(!is.na(temp[totalday <31,get(PFnames[j])]))
  }
  }
  }
if(i == dim(freq.table)[1]) close(pb)
}


## create table with mean and SD of assessment per ID for every question
descr.table = data.frame(id = unique(data.descr$identity_id))
descr.table[,c(DEPnames,PFnames)] = 0

#test = data.frame(table(data[id ==descr.table$id[110],answers_question_id]))


pb = txtProgressBar(min = 0, max = dim(descr.table)[1], initial = 0)
for(i in 1:dim(descr.table)[1]){
  temp = data.descr[identity_id ==descr.table$id[i],]

setTxtProgressBar(pb,i)
   if(length(temp) > 0){
      if(dim(temp)[1] >0){
      for (j in 1:length(DEPnames)){
    descr.table[i,DEPnames[j]] = mean(temp[totalday <14,get(DEPnames[j])], na.rm = T)
      }
        for (j in 1:length(PFnames)){
    descr.table[i,PFnames[j]] = mean(temp[totalday <31,get(PFnames[j])], na.rm = T)
      
  
  }
  }
  }
if(i == dim(descr.table)[1]) close(pb)
}


save(descr.table, freq.table, file= "itemdescrfreqs.rda")

```



# Depression indicators: two core and 1 additional symptoms indicator
```{r}

# depression parcels Core and additional
# Dep1: depressed mood
# Dep2: reduced interest or pleasure in activities previously enjoyed (anhedonia)
# Dep3: Additional symptoms difficulty concentrating, feelings of worthlessness or guilt, recurrent thoughts of death or suicide, changes in appetite or sleep, and reduced energy

# Dep
days = 14
data[,DepP1 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_001" , "question_002" )]

days = 14
data[,DepP2 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_003" , "question_004" )]

days = 14
data[,DepP3 := mean(rowMeans(.SD[totalday <days], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_012","question_010","question_009","question_016","question_014","question_013","question_005")]



# Dep outcome

data[,ODepP1 := mean(rowMeans(.SD[totalday >31], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_001" , "question_002")]

data[,ODepP2 := mean(rowMeans(.SD[totalday >31], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_003" , "question_004")]

data[,ODepP3 := mean(rowMeans(.SD[totalday >31], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_012","question_010","question_009","question_016","question_014","question_013","question_005")]

```



## PF outcome data
```{r}


## funcoutcome find ids with available func after one year
days = 299
data[,FUNC1Y := mean(rowMeans(.SD[totalday > days], na.rm = T),na.rm = T),
         by = .(identity_id),
         .SDcols = c("question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]


# average over one year PF outcome (for users with available func over one year)
days = 31
data[,Ofunc1 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_138")]
data[,Ofunc2 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_193")]
data[,Ofunc3 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_194")]
data[,Ofunc4 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_195")]
data[,Ofunc5 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_139")]
data[,Ofunc6 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_199")]
data[,Ofunc7 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_200")]
data[,Ofunc8 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_201")]
data[,Ofunc9 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_140")]
data[,Ofunc10 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_196")]
data[,Ofunc11 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_197")]
data[,Ofunc12 := mean(rowMeans(.SD[!is.na(FUNC1Y) & totalday > 31], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_198")]
```



## create mplus frame
```{r}
# create a data frame that contains the time for which we want to model the situation-emotion-dynamics, at least one month usage and at least 30 assessments
data.mplus <- data[totalday < 31]

# exclude MP test users (hundreds of assessments within two weeks)
data.mplus = data.mplus[!(data.mplus$identity_id %in% freq.table[freq.table[,-1]>13,"id"][!is.na(freq.table[freq.table[,-1]>13,"id"])]),]

# more than 30 assessments
data.mplus[, N_ass := length(totalday),  by = .(identity_id) ]
# remove rows without mood or situations assessments
data.mplus <- data.mplus[totalday < 31 & N_ass > 30 & !is.na(happiness_score)]
# adjust N_ass so that it only counts mood assessments
data.mplus[, N_ass := length(totalday),  by = .(identity_id) ]
data.mplus <- data.mplus[totalday < 31 & N_ass > 30 & !is.na(happiness_score)]



# create discrete time variable that advances +2 if night in between (input time in days as decimals)
discrtime = function(asstimes){
  discrtimes = 1
  for (i in 2:length(asstimes)){
    # if time between assessments is < 1 day advance 1 if smaller than .4 or 3 discrete time points when >.4 (when night in between)
    if((asstimes[i]-asstimes[i-1]) <1){
      if((asstimes[i]-asstimes[i-1]) <.4){
        discrtimes[i] = discrtimes[i-1] +1
      }else{
        discrtimes[i] = discrtimes[i-1] +2
      }
      # if more than one day, advance days/.25 times
    } else{
      discrtimes[i] = discrtimes[i-1] + round(((asstimes[i]-asstimes[i-1])/.25), digits = 0)
    }
  }
  return(discrtimes)
}

# temp = data.mplus[,.(identity_id,assessment_time,totalday)]
# 
#  temp[,asstime_discr := discrtime(totalday), by = .(identity_id) ]
# view(temp[,.(assessment_time,totalday,asstime_discr)])


#data.mplus[,asstime_2 := as.numeric(round((assessment_time - min(assessment_time)) / 14400,0)), by = .(identity_id) ]

data.mplus[,asstime_discr := discrtime(totalday), by = .(identity_id)]


# 
# ## check if works
# timecheck = function(asstimes){
#   error = FALSE
#   for (i in 2:length(asstimes)){
#     # if time between assessments is < 1 day advance 1 or 2 discrete time points
#     if((asstimes[i]-asstimes[i-1]) <= 0){error = T
#     cat(i)} 
#     
#   }
#   return(error)
# }
# 
# data.mplus[,checktime := timecheck(asstime_discr), by = .(identity_id)]
# 
# tests = data.mplus[!duplicated(identity_id),.(identity_id,checktime)]
# tests = tests[checktime == T,]
# 
# temp = data.mplus[identity_id == tests$identity_id[5], .(assessment_time,asstime_discr,totalday)]
# 
# #check
# asstimes = data.mplus[identity_id == tests$identity_id[5], asstime_discr]


# 
# # minutes for RDSEM
# data.mplus[,asstime_minutes := as.numeric(assessment_time - min(assessment_time))/60, by = .(identity_id)]
# # 
# # check if still duplicated mood assessments
# data.mplus[,m_ass := length(!duplicated(asstime_minutes)), by = .(identity_id)]
# # data.mplus[, N_ass := length(totalday),  by = .(identity_id) ]
# # 
# # temp = data.mplus[1:1000, .(identity_id,m_ass,N_ass)]
# # 
# sum(data.mplus$N_ass != data.mplus$m_ass)

```

## descriptives for analysis sample
```{r}
descr.dep = data.mplus[totalday < 14, c("question_001" , "question_002" ,"question_003" , "question_004" , "question_005","question_009" , "question_010" , "question_016" , "question_012" , "question_013" , "question_014")]

c("question_001","question_002", "question_003", "question_004", "question_005", "question_012", "question_013","question_014","question_009","question_010","question_016")

names(descr.dep) = c("DEP_depr_mood_1","DEP_hopelessness_1","DEP_dimin_interest_1","DEP_dimin_interest_2","DEP_reduced_energy_1","DEP_feelings_guilt", "DEP_feelings_worthlessn", "DEP_sucidal_ideation","DEP_concentration_2","DEP_insomnia","DEP_loss_appetite" )

descr.PF = data.mplus[totalday < 31, c("question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]

names(descr.PF) = c("PF_relationshp_mod_1","PF_relationshp_mod_2","PF_relationshp_mod_3","PF_relationshp_mod_4","PF_self_1","PF_self_2","PF_self_3","PF_self_4", "PF_contact_1","PF_contact_2","PF_contact_3","PF_contact_4" )


described = rbind(describe(descr.dep, ranges = FALSE, interp = FALSE), 
describe(descr.PF, ranges = FALSE, interp = FALSE))

freq.table = freq.table[freq.table$id %in% unique(data.mplus$id),]
descr.table = descr.table[descr.table$id %in% unique(data.mplus$id),]


items = read.xlsx("Items.xlsx")
itemsquestions = data.frame(questions = questions, itemnames = c(DEPnames,PFnames))


## create frequencieS and descriptives for questions
descr.quests = data.frame(question = names(descr.table)[c(-1)])
for(i in 1:length(names(descr.table)[c(-1)])){
  #frequencies
  descr.quests$freqmin[i] = min(freq.table[,i+1])
  descr.quests$freqmax[i] = max(freq.table[,i+1])
  descr.quests$freqmeansd[i] = paste0(round(mean(freq.table[,i+1], na.rm = T),digits=1), " [",
                                        round(sd(freq.table[,i+1], na.rm = T), digits = 1),"]")
  #descriptives
  descr.quests$descrmeansd[i] = paste0(round(mean(descr.table[,i+1], na.rm = T),digits=1), " [",
                                        round(sd(descr.table[,i+1], na.rm = T), digits = 1),"]")
  descr.quests$skew[i] = round(described$skew[i],digits= 1)
    descr.quests$kurtosis[i] = round(described$kurtosis[i], digits = 1)

        
        
 descr.quests$item[i] = items$q_text_en[items$item == itemsquestions$questions[itemsquestions$itemnames == descr.quests$question[i]]]
}

kable(descr.quests) %>% kable_paper()

write.xlsx(descr.quests, file = "table1.xlsx")
```
## Distributions
```{r}


#plot histograms depression
ggplot(as.data.frame(descr.dep %>% gather(scale, intensity)), aes(intensity, group = scale)) +
  geom_histogram(aes(y = after_stat(density) * 1), binwidth = 1) +
  facet_wrap(~ scale, ncol = 5) +
  ylab("relative frequency") +
  xlab("intensity")



#plot histograms PF
ggplot(as.data.frame(descr.PF %>% gather(scale, intensity)), aes(intensity, group = scale)) +
  geom_histogram(aes(y = after_stat(density) * 1), binwidth = 1) +
  facet_wrap(~ scale, ncol = 5) +
  ylab("relative frequency") +
  xlab("intensity")

```


## How much assessments for longitudinal depression
```{r}
# dep2m
# nassfunc = function(data){
#   return(sum(!is.na(colSums(data, na.rm = T))))
# }
# count = 0
# test = data.frame(data.mplus[identity_id==15651,mget(c("question_001","question_002", "question_003", "question_004", "question_005", "question_006","question_011", "question_012", "question_013","question_014","question_007","question_008","question_009","question_010","question_016","question_046","question_047"))])
# for (i in 1:dim(test)[2]){
#   count = count + sum(!is.na(test[,i]))
# }

# data[dep2m == T,dep2m_Nass := sum(!is.na(.SD)), 
#      by = .(identity_id),
#      .SDcols = c("question_001","question_002", "question_003", "question_004", "question_005", "question_012", "question_013","question_014","question_009","question_010","question_016")]

data[dep1y == T,dep1y_Nass := sum(!is.na(.SD)), 
     by = .(identity_id),
     .SDcols = c("question_001","question_002", "question_003", "question_004", "question_005", "question_012", "question_013","question_014","question_009","question_010","question_016")]

data[dep1y == T,PF1y_Nass := sum(!is.na(.SD)), 
     by = .(identity_id),
     .SDcols = c("question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]


describe(data[identity_id %in% data.mplus$id &!duplicated(identity_id),dep1y_Nass])

describe(data[identity_id %in% data.mplus$id &!duplicated(identity_id),PF1y_Nass])

```


## Graph with amount of total dep, PF, mood and event asssessments
```{r}

data.flow = data[,c("totalday","happiness_score","situations","question_001" , "question_002","question_003" , "question_004","question_012","question_010","question_009","question_016","question_014","question_013","question_005","question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]

data.flow$day = round(data.flow$totalday, digits = 0)
data.flow = data.flow[day > 0 & day < 366]

data.flow = data.flow[order(data.flow$day),]


data.flow[,mood := sum(!is.na(happiness_score)),
     by = .(day)]

data.flow[,events := sum(situations != "", na.rm = T),
     by = .(day)]

data.flow[,PF := sum(!is.na(.SD)),
     by = .(day),
     .SDcols = c("question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]

data.flow[,DEP := sum(!is.na(.SD)),
     by = .(day),
     .SDcols = c("question_001" , "question_002","question_003" , "question_004","question_012","question_010","question_009","question_016","question_014","question_013","question_005")]

data.flow = data.flow[!duplicated(day),.(day,mood,events,PF,DEP)]

data.flow = gather(data.flow, key = "construct", value = "assessments", -day)

library(scales)

# Create plot 
ggplot(data.flow, aes(x = day, y = assessments, colour = construct) ) +
  geom_line( )+
  labs(x = "Days", y = "Number of assessments")  +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 5))
```






# Two months outcome
```{r}

# depression parcels Core and additional
# Dep1: depressed mood
# Dep2: reduced interest or pleasure in activities previously enjoyed (anhedonia)
# Dep3: Additional symptoms difficulty concentrating, feelings of worthlessness or guilt, recurrent thoughts of death or suicide, changes in appetite or sleep, and reduced energy


# Dep outcome (average over 2 months)

data[,ODep2MP1 := mean(rowMeans(.SD[dep2m == T & totalday >31 & totalday < 73], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_001" , "question_002")]

data[,ODep2MP2 := mean(rowMeans(.SD[dep2m == T & totalday >31 & totalday < 73], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_003" , "question_004")]

data[,ODep2MP3 := mean(rowMeans(.SD[dep2m == T & totalday >31 & totalday < 73], na.rm = T),na.rm = T),
     by = .(identity_id),
     .SDcols = c("question_012","question_010","question_009","question_016","question_014","question_013","question_005")]



# average over one year PF outcome (for users with available func over one year)

## funcoutcome find ids with available func after one year
data[,FUNC2M := mean(rowMeans(.SD[totalday > 31 & totalday < 73], na.rm = T),na.rm = T),
         by = .(identity_id),
         .SDcols = c("question_138","question_193", "question_194", "question_195","question_139", "question_199", "question_200", "question_201","question_140","question_196","question_197","question_198")]

days = 31
data[,OPF2M1 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_138")]
data[,OPF2M2 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_193")]
data[,OPF2M3 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_194")]
data[,OPF2M4 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_195")]
data[,OPF2M5 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_139")]
data[,OPF2M6 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_199")]
data[,OPF2M7 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_200")]
data[,OPF2M8 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_201")]
data[,OPF2M9 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_140")]
data[,OPF2M10 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c( "question_196")]
data[,OPF2M11 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_197")]
data[,OPF2M12 := mean(rowMeans(.SD[!is.na(FUNC2M) & totalday > 31 & totalday < 73], na.rm = T),na.rm = T), 
         by = .(identity_id),
         .SDcols = c("question_198")]


```



### Prepare data for Mplus ONE MONTH AVErage PRED ITEMS ONE YEAR AVERAGE OUTCOME ITEMS PARCEL
```{r}

# prepare data, keep only relevant variables:
# id, situations, mood, personality functioning, depression at T1, average depression, assessment time variable
data.mplus <- as.data.frame(data.mplus[,list(identity_id, sits,sit_pos,sit_neg, mood, func1, func2, func3,func4, func5, func6, func7, func8, func9, func10, func11, func12, Ofunc1, Ofunc2, Ofunc3,Ofunc4, Ofunc5, Ofunc6, Ofunc7, Ofunc8, Ofunc9, Ofunc10, Ofunc11, Ofunc12,DepP1,DepP2,DepP3,ODepP1,ODepP2,ODepP3,ODep2MP1,ODep2MP2,ODep2MP3, OPF2M1, OPF2M2, OPF2M3,OPF2M4, OPF2M5, OPF2M6, OPF2M7, OPF2M8, OPF2M9, OPF2M10, OPF2M11, OPF2M12, DepAV2M,DepAV1Y, asstime_discr , dep1y, dep2m,sit_pwo, sit_psuc, sit_prel, sit_ncnf, sit_novl, sit_nemp)])
names(data.mplus) = c("id", "sits","sit_pos","sit_neg", "mood", "func1", "func2", "func3","func4", "func5", "func6", "func7", "func8", "func9", "func10", "func11", "func12", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6", "Ofunc7", "Ofunc8", "Ofunc9", "Ofunc10", "Ofunc11", "Ofunc12" ,"DepP1","DepP2","DepP3","ODepP1","ODepP2","ODepP3","ODep2MP1","ODep2MP2","ODep2MP3","OPF2M1", "OPF2M2", "OPF2M3","OPF2M4", "OPF2M5", "OPF2M6", "OPF2M7", "OPF2M8", "OPF2M9", "OPF2M10", "OPF2M11", "OPF2M12" ,"DepAV2M", "DepAV1Y", "asstimeD" , "dep1y", "dep2m","sit_pwo", "sit_psuc", "sit_prel", "sit_ncnf", "sit_novl", "sit_nemp")
# id = identity variable for each person
# sits = situation variable that is centered around zero
# Mood = mood-variable (based on smiley-scale) that is centered around zero
#  = average on all depression items in the first 14 days
# Dep_AV = average on all depression items from day 56 to day 236
# asstimeD = variable with 4h time intervalls


```




# estimate in Mplus

## Prediction with pos neg sits ESS = 400

## Single situation items not outcome ESS = 400
```{r}


# Bifactor S-1, PF as Ref, no Outcome, sit_pwo 
ModNoOutS1PFsit_pwo = 
  mplusObject(TITLE = "ModNoOutS1PFsit_pwo;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_pwo DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_pwo(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_pwo ON sit_pwo&1;
              m_s1 | mood ON sit_pwo&1;
              s_m1 | sit_pwo ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_pwo; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_pwo@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_pwo on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_pwo","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_pwo, modelout = "ModNoOutS1PFsit_pwo.inp",
             run = F)

runModels("ModNoOutS1PFsit_pwo.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")




# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_psuc
ModNoOutS1PFsit_psuc = 
  mplusObject(TITLE = "ModNoOutS1PFsit_psuc;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_psuc DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_psuc(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_psuc ON sit_psuc&1;
              m_s1 | mood ON sit_psuc&1;
              s_m1 | sit_psuc ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_psuc; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_psuc@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_psuc on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_psuc","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_psuc, modelout = "ModNoOutS1PFsit_psuc.inp",
             run = F)

runModels("ModNoOutS1PFsit_psuc.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_prel
ModNoOutS1PFsit_prel = 
  mplusObject(TITLE = "ModNoOutS1PFsit_prel;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_prel DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_prel(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_prel ON sit_prel&1;
              m_s1 | mood ON sit_prel&1;
              s_m1 | sit_prel ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_prel; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_prel@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_prel on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_prel","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_prel, modelout = "ModNoOutS1PFsit_prel.inp",
             run = F)

runModels("ModNoOutS1PFsit_prel.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")





# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_ncnf
ModNoOutS1PFsit_ncnf = 
  mplusObject(TITLE = "ModNoOutS1PFsit_ncnf;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_ncnf, modelout = "ModNoOutS1PFsit_ncnf.inp",
             run = F)

runModels("ModNoOutS1PFsit_ncnf.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")





# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_novl
ModNoOutS1PFsit_novl = 
  mplusObject(TITLE = "ModNoOutS1PFsit_novl;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_novl DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_novl(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_novl ON sit_novl&1;
              m_s1 | mood ON sit_novl&1;
              s_m1 | sit_novl ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_novl; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_novl@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_novl on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_novl","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_novl, modelout = "ModNoOutS1PFsit_novl.inp",
             run = F)

runModels("ModNoOutS1PFsit_novl.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")




# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_nemp
ModNoOutS1PFsit_nemp = 
  mplusObject(TITLE = "ModNoOutS1PFsit_nemp;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_nemp DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_nemp(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_nemp ON sit_nemp&1;
              m_s1 | mood ON sit_nemp&1;
              s_m1 | sit_nemp ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_nemp; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_nemp@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_nemp on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_nemp","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_nemp, modelout = "ModNoOutS1PFsit_nemp.inp",
             run = F)

runModels("ModNoOutS1PFsit_nemp.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")










# Bifactor S-1, DEP as Ref, no Outcome, sit_pwo 
ModNoOutS1DEPsit_pwo = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_pwo;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_pwo DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_pwo(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_pwo ON sit_pwo&1;
              m_s1 | mood ON sit_pwo&1;
              s_m1 | sit_pwo ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_pwo; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_pwo@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_pwo on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_pwo","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_pwo, modelout = "ModNoOutS1DEPsit_pwo.inp",
             run = F)

runModels("ModNoOutS1DEPsit_pwo.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")




# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_psuc
ModNoOutS1DEPsit_psuc = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_psuc;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_psuc DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_psuc(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_psuc ON sit_psuc&1;
              m_s1 | mood ON sit_psuc&1;
              s_m1 | sit_psuc ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_psuc; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_psuc@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_psuc on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_psuc","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_psuc, modelout = "ModNoOutS1DEPsit_psuc.inp",
             run = F)

runModels("ModNoOutS1DEPsit_psuc.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_prel
ModNoOutS1DEPsit_prel = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_prel;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_prel DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_prel(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_prel ON sit_prel&1;
              m_s1 | mood ON sit_prel&1;
              s_m1 | sit_prel ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_prel; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_prel@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_prel on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_prel","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_prel, modelout = "ModNoOutS1DEPsit_prel.inp",
             run = F)

runModels("ModNoOutS1DEPsit_prel.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")





# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_ncnf
ModNoOutS1DEPsit_ncnf = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_ncnf;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_ncnf, modelout = "ModNoOutS1DEPsit_ncnf.inp",
             run = F)

runModels("ModNoOutS1DEPsit_ncnf.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")





# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_novl
ModNoOutS1DEPsit_novl = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_novl;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_novl DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_novl(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_novl ON sit_novl&1;
              m_s1 | mood ON sit_novl&1;
              s_m1 | sit_novl ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_novl; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_novl@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_novl on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_novl","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_novl, modelout = "ModNoOutS1DEPsit_novl.inp",
             run = F)

runModels("ModNoOutS1DEPsit_novl.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")




# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_nemp
ModNoOutS1DEPsit_nemp = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_nemp;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_nemp DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_nemp(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_nemp ON sit_nemp&1;
              m_s1 | mood ON sit_nemp&1;
              s_m1 | sit_nemp ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_nemp; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_nemp@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_nemp on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_nemp","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_nemp, modelout = "ModNoOutS1DEPsit_nemp.inp",
             run = F)

runModels("ModNoOutS1DEPsit_nemp.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")











###################################################### RERUN with MODELPRIORS for SM1 because of convergence problems

# Bifactor S-1 wo change,  PF as Ref, no Outcome, sit_ncnf WITH MODEL PRIORS
ModNoOutS1PFsit_ncnf = 
  mplusObject(TITLE = "ModNoOutS1PFsit_ncnf;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.01;
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1PFsit_ncnf, modelout = "ModNoOutS1PFsit_ncnf.inp",
             run = F)

runModels("ModNoOutS1PFsit_ncnf.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")






# Bifactor S-1 wo change,  DEP as Ref, no Outcome, sit_ncnf
ModNoOutS1DEPsit_ncnf = 
  mplusObject(TITLE = "ModNoOutS1DEPsit_ncnf;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01847585;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.002;
               !Personality functioning factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 func4 func5
              func6 func7 func8 func9 func10 func11 func12;
              FUNC*; DEP*;
              !covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12", "asstimeD")])

mplusModeler(ModNoOutS1DEPsit_ncnf, modelout = "ModNoOutS1DEPsit_ncnf.inp",
             run = F)

runModels("ModNoOutS1DEPsit_ncnf.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



```



# 1Y Prediction ESS 400
```{r}


# 1 year neg sits Bifactor S-1 pred cor out, PF as ref
Model1YnempS1PCO = 
  mplusObject(TITLE = "Model1YnempS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_nemp DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_nemp(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_nemp ON sit_nemp&1;
              m_s1 | mood ON sit_nemp&1;
              s_m1 | sit_nemp ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_nemp; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_nemp@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              !same for outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_nemp on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_nemp FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_nemp","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YnempS1PCO, modelout = "Model1YnempS1PCO.inp",
             run = F)

runModels("Model1YnempS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")


# 1 year novl sits Bifactor S-1 pred cor out, PF as ref
Model1YnovlS1PCO = 
  mplusObject(TITLE = "Model1YnovlS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_novl DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_novl(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_novl ON sit_novl&1;
              m_s1 | mood ON sit_novl&1;
              s_m1 | sit_novl ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_novl; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_novl@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              !same for outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_novl on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_novl FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_novl","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YnovlS1PCO, modelout = "Model1YnovlS1PCO.inp",
             run = F)

runModels("Model1YnovlS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year ncnf sits Bifactor S-1 pred cor out, PF as ref
Model1YncnfS1PCO = 
  mplusObject(TITLE = "Model1YncnfS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.001;
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              !same for outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_ncnf FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YncnfS1PCO, modelout = "Model1YncnfS1PCO.inp",
             run = F)

runModels("Model1YncnfS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")






# 1 year pos Bifactor S-1 pred cor out CHANGE (baseline DEp and PF predicts PF and DEP -> change), PF as reference
Model1YprelS1PCO = 
  mplusObject(TITLE = "Model1YprelS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_prel DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_prel(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_prel ON sit_prel&1;
              m_s1 | mood ON sit_prel&1;
              s_m1 | sit_prel ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_prel; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_prel@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.001;
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              !same for outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_prel on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_prel FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_prel","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YprelS1PCO, modelout = "Model1YprelS1PCO.inp",
             run = F)

runModels("Model1YprelS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")
# 

# 1 year psuc Bifactor S-1 pred cor out CHANGE (baseline DEp and PF predicts PF and DEP -> change), PF as reference
Model1YpsucS1PCO =
  mplusObject(TITLE = "Model1YpsucS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_psuc DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8
              Ofunc9 Ofunc10 Ofunc11 Ofunc12
              asstimeD;
              LAGGED = mood(1) sit_psuc(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES;
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_psuc ON sit_psuc&1;
              m_s1 | mood ON sit_psuc&1;
              s_m1 | sit_psuc ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_psuc; !Random unique innovation variance for sits
              m_s_lv by mood@1 sit_psuc@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by
              func1* func2 func3 func4 func5@1 func6 func7 func8
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by
              DepP1@1 DepP2 DepP3;
              !same for outcome
              OFUNC by
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 Ofunc8
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression
              ODEP by
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_psuc on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_psuc FUNC DEP;
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_psuc","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YpsucS1PCO, modelout = "Model1YpsucS1PCO.inp",
             run = F)

runModels("Model1YpsucS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year pwo Bifactor S-1 pred cor out CHANGE (baseline DEp and PF predicts PF and DEP -> change), PF as reference
Model1YpwoS1PCO = 
  mplusObject(TITLE = "Model1YpwoS1PCO;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_pwo DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_pwo(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_pwo ON sit_pwo&1;
              m_s1 | mood ON sit_pwo&1;
              s_m1 | sit_pwo ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_pwo; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_pwo@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning reference factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12 DepP1 DepP2 DepP3;
              !Depression specific residual factor
              DEP by 
              DepP1@1 DepP2 DepP3;
              !outcome correlated factors
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_pwo on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_pwo FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_pwo","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YpwoS1PCO, modelout = "Model1YpwoS1PCO.inp",
             run = F)

runModels("Model1YpwoS1PCO.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")






# 1 year novl sits Bifactor S-1 pred cor out, Dep as ref
Model1YnovlS1PCODep = 
  mplusObject(TITLE = "Model1YnovlS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_novl DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_novl(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_novl ON sit_novl&1;
              m_s1 | mood ON sit_novl&1;
              s_m1 | sit_novl ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_novl; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_novl@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_novl on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_novl FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_novl","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YnovlS1PCODep, modelout = "Model1YnovlS1PCODep.inp",
             run = F)

runModels("Model1YnovlS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year nemp sits Bifactor S-1 pred cor out, Dep as ref
Model1YnempS1PCODep = 
  mplusObject(TITLE = "Model1YnempS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_nemp DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_nemp(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_nemp ON sit_nemp&1;
              m_s1 | mood ON sit_nemp&1;
              s_m1 | sit_nemp ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_nemp; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_nemp@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_nemp on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_nemp FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_nemp","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YnempS1PCODep, modelout = "Model1YnempS1PCODep.inp",
             run = F)

runModels("Model1YnempS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year ncnf sits Bifactor S-1 pred cor out, Dep as ref
Model1YncnfS1PCODep = 
  mplusObject(TITLE = "Model1YncnfS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_ncnf DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_ncnf(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062; 
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_ncnf ON sit_ncnf&1;
              m_s1 | mood ON sit_ncnf&1;
              s_m1 | sit_ncnf ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_ncnf; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_ncnf@-1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.002;
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_ncnf on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_ncnf FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_ncnf","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YncnfS1PCODep, modelout = "Model1YncnfS1PCODep.inp",
             run = F)

runModels("Model1YncnfS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year prel sits Bifactor S-1 pred cor out, Dep as ref
Model1YprelS1PCODep = 
  mplusObject(TITLE = "Model1YprelS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_prel DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_prel(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062; 
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_prel ON sit_prel&1;
              m_s1 | mood ON sit_prel&1;
              s_m1 | sit_prel ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_prel; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_prel@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
              s_m1@0.001;
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_prel on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_prel FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_prel","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YprelS1PCODep, modelout = "Model1YprelS1PCODep.inp",
             run = F)

runModels("Model1YprelS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")



# 1 year psuc sits Bifactor S-1 pred cor out, Dep as ref
Model1YpsucS1PCODep = 
  mplusObject(TITLE = "Model1YpsucS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_psuc DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_psuc(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 24; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_psuc ON sit_psuc&1;
              m_s1 | mood ON sit_psuc&1;
              s_m1 | sit_psuc ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_psuc; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_psuc@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_psuc on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_psuc FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_psuc","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YpsucS1PCODep, modelout = "Model1YpsucS1PCODep.inp",
             run = F)

runModels("Model1YpsucS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")




# 1 year pwo sits Bifactor S-1 pred cor out, Dep as ref
Model1YpwoS1PCODep = 
  mplusObject(TITLE = "Model1YpwoS1PCODep;",
              VARIABLE = "
              CLUSTER = id;
               USEVARIABLES = id mood sit_pwo DepP1 DepP2 DepP3 
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3 
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12 
              asstimeD;
              LAGGED = mood(1) sit_pwo(1);
               BETWEEN = DepP1 DepP2 DepP3
              func1 func2 func3 func4 func5 func6 func7 func8 
              func9 func10 func11 func12 ODepP1 ODepP2 ODepP3
              Ofunc1 Ofunc2 Ofunc3 Ofunc4 Ofunc5 Ofunc6 Ofunc7 Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              tinterval = asstimeD(1);",
              ANALYSIS ="type = TWOLEVEL RANDOM; estimator = BAYES; 
               thin = 10; PROCESS = 28; BCONVERGENCE= 0.01571062;   
              BITERATIONS = 100000 20000;",
              MODEL = "%WITHIN%
              m_m | mood ON mood&1;
              s_s | sit_pwo ON sit_pwo&1;
              m_s1 | mood ON sit_pwo&1;
              s_m1 | sit_pwo ON mood&1;
              miv | mood; !Random unique innovation variance for mood
              siv | sit_pwo; !Random unique innovation variance for sits	
              m_s_lv by mood@1 sit_pwo@1; !Latent covariance between innovations
              m_s | m_s_lv; !Random residual covariance
              %BETWEEN%
               !Personality functioning residual factor
              FUNC by 
              func1* func2 func3 func4 func5@1 func6 func7 func8 
              func9 func10 func11 func12;
              !Depression specific reference factor
              DEP by 
              DepP1@1 DepP2 DepP3 func1 func2 func3 
              func4 func5 func6 func7 func8 
              func9 func10 func11 func12;
              !corrfac as outcome
              OFUNC by 
              Ofunc1* Ofunc2 Ofunc3 Ofunc4 Ofunc5@1 Ofunc6 Ofunc7 
              Ofunc8 
              Ofunc9 Ofunc10 Ofunc11 Ofunc12;
              !Depression 
              ODEP by 
              OdepP1@1 ODepP2 ODepP3;
              !factor variances
              OFUNC*; ODEP*;
              FUNC*; DEP*;
              ! pred covariance/correlation fixed to 0
              FUNC with DEP@0;
              m_s miv siv m_m s_s m_s1 s_m1 mood sit_pwo on
              FUNC DEP;
              ODEP OFUNC on m_m s_s m_s1 s_m1 m_s miv siv
              mood sit_pwo FUNC DEP; 
              m_s;",
              OUTPUT = "TECH1 STANDARDIZED TECH4 TECH8 STDYX;",
              #PLOT = "type = PLOT3; factors = ALL;",
              rdata = data.mplus[data.mplus$dep1y == T,c("id","mood","sit_pwo","DepP1", "DepP2", "DepP3", "func1", "func2", "func3","func4", "func5", "func6","func7", "func8", "func9","func10", "func11", "func12","ODepP1", "ODepP2", "ODepP3", "Ofunc1", "Ofunc2", "Ofunc3","Ofunc4", "Ofunc5", "Ofunc6","Ofunc7", "Ofunc8", "Ofunc9","Ofunc10", "Ofunc11", "Ofunc12", "asstimeD")])

mplusModeler(Model1YpwoS1PCODep, modelout = "Model1YpwoS1PCODep.inp",
             run = F)

runModels("Model1YpwoS1PCODep.inp", showOutput = TRUE,Mplus_command = "/opt/mplus/8.11/mplus")







```












# Extract model params
```{r}
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,
                      onlywithin = F) {
  
rmzero <- function(x, digits = NULL) {
  s <- if (is.null(digits)) {
    as.character(x)
  } else {
    format(x, nsmall = digits)
  }
  sub("^(-?)0\\.", "\\1.", s)
}
  if (length(adds1models) == 1) {
    # read mplus model
    mplodel = readModels(target = mplodelnames[1])
    #extract std params
    stdout = mplodel$parameters$stdyx.standardized
    #delete unnecessary rows
    if(!onlywithin){
    stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
    stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
    }else if(onlywithin){
      stdout = stdout[!stdout$BetweenWithin == "Between",]
    }
    
    rownames(stdout) = 1:dim(stdout)[1]
    stdout$est = as.numeric(stdout$est)
    stdout$lower_2.5ci = as.numeric(stdout$lower_2.5ci)
    stdout$upper_2.5ci = as.numeric(stdout$upper_2.5ci)
    stdout$pval = as.numeric(stdout$pval)
    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], "<br> N=", N))
    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 =   as.numeric(temp[temp$paramHeader == "M_S_LV.BY", "est"][1])
      loads =  as.numeric(temp[temp$paramHeader == "M_S_LV.BY", "est"][2])
      # covariance between residual variances
      covms = as.numeric(loads) * as.numeric(loadm)
      # residual variances
      varresidm = as.numeric(temp[temp$paramHeader == "MIV.|", "est"]) + loadm ^ 2
      varresids = as.numeric(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(as.numeric(stdout[grep("M_S_LV", stdout$param), "est"])) ^ .3) * 2.25152 * (N ^ -0.50089)
      stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = as.numeric(stdout[grep("M_S_LV", stdout$param), "est"]) - CI
      stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = as.numeric(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 = ";"),
          "\" >",
          rmzero(round(stdout$est[!noconf], digits = 2)),
          " <br> [",
          rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
          " ",
          rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
          "]</span>"
        )}
       if (html == F) {
        table2[!noconf, 2] = ifelse(
          stdout$pval < .005,
          yes = paste0(
            rmzero(round(stdout$est[!noconf], digits = 2)),
            " <br> [",
            rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
            " ",
            rmzero(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 = ";"),
        "\" >",
        rmzero(round(stdout$est[noconf], digits = 2)),
        "</span>"
      )
    if (html == F) {
      table2[noconf, 2] = ifelse(stdout$pval[noconf] < .005, yes = paste0(rmzero(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
     if(!onlywithin){
      stdout = stdout[-grep("esidual", stdout$paramHeader) , ]
      stdout = stdout[-grep("ntercept", stdout$paramHeader) , ]
       }else if(onlywithin){
      stdout = stdout[!stdout$BetweenWithin == "Between",]}
      
      rownames(stdout) = 1:dim(stdout)[1]
          stdout$est = as.numeric(stdout$est)
    stdout$lower_2.5ci = as.numeric(stdout$lower_2.5ci)
    stdout$upper_2.5ci = as.numeric(stdout$upper_2.5ci)
    stdout$pval = as.numeric(stdout$pval)
      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 =   as.numeric(temp[temp$paramHeader == "M_S_LV.BY", "est"][1])
        loads =  as.numeric(temp[temp$paramHeader == "M_S_LV.BY", "est"][2])
        # covariance between residual variances
        covms = loads * loadm
        # residual variances
        varresidm = as.numeric(temp[temp$paramHeader == "MIV.|", "est"]) + loadm ^ 2
        varresids = as.numeric(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(as.numeric(stdout[grep("M_S_LV", stdout$param), "est"])) ^ .3) * 2.25152 * (N ^ -0.50089)
        stdout[grep("M_S_LV", stdout$param), "lower_2.5ci"] = as.numeric(stdout[grep("M_S_LV", stdout$param), "est"]) - CI
        stdout[grep("M_S_LV", stdout$param), "upper_2.5ci"] = as.numeric(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 = ";"),
          "\" >",
          rmzero(round(stdout$est[!noconf], digits = 2)),
          " <br> [",
          rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
          " ",
          rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
          "]</span>"
        )
       if (html == F) {
        table2[!noconf, i + 1] = ifelse(
          stdout$pval[!noconf] < .005,
          yes = paste0(
            rmzero(round(stdout$est[!noconf], digits = 2)),
            " <br> [",
            rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
            " ",
            rmzero(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 = ";"),
        "\" >",
        rmzero(round(stdout$est[noconf], digits = 2)),
        "</span>"
      )
    if (html == F) {
      table2[noconf, i + 1] = ifelse(stdout$pval[noconf] < .005, yes = paste0(rmzero(round(stdout$est[noconf], digits = 2))), no = "")
    }
      names(table2)[i + 1] = paste0(tablenames[i], "<br> N=", N)
      # 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], "<br> N=", N))
    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 = ";"),
    "\" >",
    rmzero(round(rowMeans(cbind(stdout$est[!noconf],stdouts1$est[!noconf])), digits = 2)),
    " <br> [",
    rmzero(round(rowMeans(cbind(stdout$lower_2.5ci[!noconf],stdouts1$lower_2.5ci[!noconf])), digits = 2)),
    " ",
    rmzero(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(
            rmzero(round(rowMeans(cbind(stdout$est[!noconf],stdouts1$est[!noconf])), digits = 2)),
            " [",
            rmzero(round(rowMeans(cbind(stdout$lower_2.5ci[!noconf],stdouts1$lower_2.5ci[!noconf])), digits = 2)),
            " ",
            rmzero(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) {
      # bivariate unresidualized first in cell
      table2[noconf & stdout$param == "FUNC", 2] = paste0(
        "<span>",
        ifelse(stdout$pval[noconf & stdout$param == "FUNC"] < .005, yes = "<b>", no = ""),
        rmzero(round(stdout$est[noconf & stdout$param == "FUNC"], digits = 2)),
        "</b>",
        ifelse(stdouts1$pval[noconf & stdouts1$param == "FUNC"] < .005, yes = "</b><b>", no = "</b>"),"[",
        rmzero(round(stdouts1$est[noconf & stdouts1$param == "FUNC"], digits = 2)),
        "] ",
        ifelse((!(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"])) | !(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]))) & stdout$pval[noconf & stdout$param == "FUNC"] < .005 & stdouts1$pval[noconf & stdout$param == "FUNC"] < .005, yes = "**", no =  ifelse((!(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"])) | !(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]))) & !(stdout$pval[noconf & stdout$param == "FUNC"] < .005 & stdouts1$pval[noconf & stdout$param == "FUNC"] < .005), yes = "*", no = "")),
        "</span>")
        
        table2[noconf & stdout$param == "DEP", 2] = paste0(
        "<span>",
        ifelse(stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "</b><b>", no = "</b>"),
        rmzero(round(stdouts1$est[noconf & stdout$param == "DEP"], digits = 2)),
        "</b>",
        ifelse(stdout$pval[noconf & stdout$param == "DEP"] < .005, yes = "<b>", no = ""),"[",
        rmzero(round(stdout$est[noconf & stdout$param == "DEP"], digits = 2)),
        "] ",
        ifelse((!(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"])) | !(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]))) & stdout$pval[noconf & stdout$param == "DEP"] < .005 & stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "**", no = ifelse((!(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"])) | !(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]))) & !(stdout$pval[noconf & stdout$param == "DEP"] < .005 & stdouts1$pval[noconf & stdout$param == "DEP"] < .005), yes = "*", no = "")),
        "</span>"
      )}
    if (html == F) {
      table2[noconf, 2] = ifelse(stdout$pval[noconf] < .005, yes = paste0(rmzero(round(stdout$est[noconf], digits = 2)), rmzero(round(stdouts1$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] = 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 = ";"),
    "\" >",
    rmzero(round(rowMeans(cbind(stdout$est,stdouts1$est))[!noconf], digits = 2)),
    " <br> [",
    rmzero(round(rowMeans(cbind(stdout$lower_2.5ci,stdouts1$lower_2.5ci))[!noconf], digits = 2)),
    " ",
    rmzero(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(
            rmzero(round(rowMeans(cbind(stdout$est,stdouts1$est))[!noconf], digits = 2)),
            " <br> [",
            rmzero(round(rowMeans(cbind(stdout$lower_2.5ci,stdouts1$lower_2.5ci))[!noconf], digits = 2)),
            " ",
            rmzero(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 & stdout$param == "FUNC", i+1] = paste0(
        "<span>",
        ifelse(stdout$pval[noconf & stdout$param == "FUNC"] < .005, yes = "<b>", no = ""),
        rmzero(round(stdout$est[noconf & stdout$param == "FUNC"], digits = 2)),
        "</b>",
        ifelse(stdouts1$pval[noconf & stdouts1$param == "FUNC"] < .005, yes = "</b><b>", no = "</b>"),"[",
        rmzero(round(stdouts1$est[noconf & stdouts1$param == "FUNC"], digits = 2)),
        "] ",
        ifelse((!(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"])) | !(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]))) & stdout$pval[noconf & stdout$param == "FUNC"] < .005 & stdouts1$pval[noconf & stdout$param == "FUNC"] < .005, yes = "**", no =  ifelse((!(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"])) | !(pmin(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "FUNC"], stdouts1$est[noconf & stdout$param == "FUNC"]))) & !(stdout$pval[noconf & stdout$param == "FUNC"] < .005 & stdouts1$pval[noconf & stdout$param == "FUNC"] < .005), yes = "*", no = "")),
        "</span>")
        
        table2[noconf & stdout$param == "DEP", i +1] = paste0(
        "<span>",
        ifelse(stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "</b><b>", no = "</b>"),
        rmzero(round(stdouts1$est[noconf & stdout$param == "DEP"], digits = 2)),
        "</b>",
        ifelse(stdout$pval[noconf & stdout$param == "DEP"] < .005, yes = "<b>", no = ""),"[",
        rmzero(round(stdout$est[noconf & stdout$param == "DEP"], digits = 2)),
        "] ",
        ifelse((!(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"])) | !(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]))) & stdout$pval[noconf & stdout$param == "DEP"] < .005 & stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "**", no = ifelse((!(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) < 0) & (0 < pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"])) | !(pmin(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]) > 0) & (0 > pmax(stdout$est[noconf & stdout$param == "DEP"], stdouts1$est[noconf & stdout$param == "DEP"]))) & !(stdout$pval[noconf & stdout$param == "DEP"] < .005 & stdouts1$pval[noconf & stdout$param == "DEP"] < .005), yes = "*", no = "")),
        "</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(rmzero(round(stdout$est[noconf], digits = 2)), rmzero(round(stdouts1$est[noconf], digits = 2))), no = "")
    }}
      names(table2)[i + 1] = paste0(tablenames[i], "<br> N=", N)
      # 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))
}



```




# Create tables
```{r}

## within table dscrete situations without mood
table = table2func(mplodelnames = c("ModWithin_sit_ncnf.out","ModWithin_sit_nemp.out","ModWithin_sit_novl.out","ModWithin_sit_prel.out","ModWithin_sit_psuc.out","ModWithin_sit_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,1,4,6,8,9), html = T,logtransparams = T,onlywithin = 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))


# 
# ## within table dscrete situations without mood
# table = table2func(mplodelnames = c("ModNoOutS1PFsit_ncnf.out","ModNoOutS1DEPsit_nemp.out","ModNoOutS1PFsit_novl.out","ModNoOutS1PFsit_prel.out","ModNoOutS1PFsit_psuc.out","ModNoOutS1PFsit_pwo.out"), tablenames = c("conflict","emptiness<br> or boredom","overwhelming <br> task","relaxation","achievement","good time<br> 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<br> or boredom","overwhelming <br> task","relaxation","achievement","good time<br> 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 old
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))



```



# 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",
  "mxe~>PF", "mxe~>DEP",
  "me1~>PF", "me1~>DEP",
  "em1~>PF", "em1~>DEP",
  "miv~>PF", "miv~>DEP",
  "eiv~>PF", "eiv~>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("PF~>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,3])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfcovrs[[y]]$cormat[1,3], digits = 2)),
          "</b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,3])$p < .005, yes = "<b>", no = "</b>"),
          "[",
          rmzero(round(depcovrs[[y]]$cormat[1,3], digits = 2)),
          "] </b>| ",
          # also put squared semipartial correlations
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,3])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfsqspcors[[y]]$OFUNC[1,2], digits = 3)),
          "</b>",
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,3])$p < .005, yes = "<b>", no = "</b>"),
          "[" ,
          # also put squared semipartial correlations
          rmzero(round(depsqspcors[[y]]$OFUNC[1,2], digits = 3)),
          "]</b>",
          "</span>"
        )
    table4[grep("DEP~>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],depcovrs[[y]]$cormat[2,3])$p < .005, yes = "<b>", no = "</b>"),
          rmzero(round(depcovrs[[y]]$cormat[2,3], digits = 2)),
          "</b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,3])$p < .005, yes = "<b>", no = ""),
          "[",
          rmzero(round(pfcovrs[[y]]$cormat[2,3], digits = 2)),
          "] </b>| ",
          # also put squared semipartial correlations
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,3])$p < .005, yes = "<b>", no = "</b>"),
          rmzero(round(depsqspcors[[y]]$OFUNC[1,3], digits = 3)),
          "</b>",
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,3])$p < .005, yes = "<b>", no = ""),
          "[" ,
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$OFUNC[1,3], digits = 3)),
          "]</b>",
          "</span>"
        )
      
    # for coloring average
    table4r[grep("~>OPF", rownamess),y] = rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,3],depcovrs[[y]]$cormat[1:2,3]))
    # unique variance explanation of unresidualized PF and DEP
    table4sqsp[grep("DEP~>OPF", rownamess),y] = depsqspcors[[y]]$OFUNC[1,3]
    table4sqsp[grep("PF~>OPF", rownamess),y] = pfsqspcors[[y]]$OFUNC[1,2]


    
     # PF OPF DEP ODEP
    table4[grep("PF~>ODEP", 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,4])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfcovrs[[y]]$cormat[1,4], digits = 2)),
          "</b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,4])$p < .005, yes = "<b>", no = "</b>"),
          "[",
          rmzero(round(depcovrs[[y]]$cormat[1,4], digits = 2)),
          "] </b>| ",
          # also put squared semipartial correlations
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,4])$p < .005, yes = "<b>", no = ""),
          rmzero(round(pfsqspcors[[y]]$ODEP[1,2], digits = 3)),
          "</b>",
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,4])$p < .005, yes = "<b>", no = "</b>"),
          "[" ,
          # also put squared semipartial correlations
          rmzero(round(depsqspcors[[y]]$ODEP[1,2], digits = 3)),
          "]</b>",
          "</span>"
        )
    
    table4[grep("DEP~>ODEP", 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],depcovrs[[y]]$cormat[2,4])$p < .005, yes = "<b>", no = "</b>"),
          rmzero(round(depcovrs[[y]]$cormat[2,4], digits = 2)),
          "</b>",
          #bold if significant correlation
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,4])$p < .005, yes = "<b>", no = ""),
          "[",
          rmzero(round(pfcovrs[[y]]$cormat[2,4], digits = 2)),
          "] </b>| ",
          # also put squared semipartial correlations
          ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,4])$p < .005, yes = "<b>", no = "</b>"),
          rmzero(round(depsqspcors[[y]]$ODEP[1,3], digits = 3)),
          "</b>",
          ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,4])$p < .005, yes = "<b>", no = ""),
          "[" ,
          # also put squared semipartial correlations
          rmzero(round(pfsqspcors[[y]]$ODEP[1,3], digits = 3)),
          "]</b>",
          "</span>"
        )
      
    table4r[grep("~>ODEP", rownamess),y] = rowMeans(cbind(pfcovrs[[y]]$cormat[1:2,4],depcovrs[[y]]$cormat[1:2,4]))
    # unique variance explanation of unresidualized PF and DEP
    table4sqsp[grep("DEP~>ODEP", rownamess),y] = depsqspcors[[y]]$ODEP[1,3]
    table4sqsp[grep("PF~>ODEP", rownamess),y] = pfsqspcors[[y]]$ODEP[1,2]

    
  }



```


##  Create prediction 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))


```


# Average R2 for PF and DEP
```{r}
# PF with PF as reference
mean(colMeans(table4sqsp[c("PF~>OPF","PF~>ODEP"),]))

# DEP with DEP as reference
mean(colMeans(table4sqsp[c("DEP~>OPF","DEP~>ODEP"),]))

# ODEP
rowMeans(table4sqsp["PF~>ODEP",])
rowMeans(table4sqsp["DEP~>ODEP",])

# OPF
rowMeans(table4sqsp["PF~>OPF",])
rowMeans(table4sqsp["DEP~>OPF",])




```


