---
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 .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 & 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 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], "
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("π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(
"",
rmzero(round(stdout$est[!noconf], digits = 2)),
"
[",
rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
" ",
rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
"]"
)}
if (html == F) {
table2[!noconf, 2] = ifelse(
stdout$pval < .005,
yes = paste0(
rmzero(round(stdout$est[!noconf], digits = 2)),
"
[",
rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
" ",
rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
"]"
),
no = ""
)
}
table2[noconf, 2] = paste0(
"",
rmzero(round(stdout$est[noconf], digits = 2)),
""
)
if (html == F) {
table2[noconf, 2] = ifelse(stdout$pval[noconf] < .005, yes = paste0(rmzero(round(stdout$est[noconf], digits = 2))), no = "")
}
#
# table2[,2] = paste0("",round(stdout$est,digits = 2),"
[",round(stdout$lower_2.5ci,digits = 2)," ",round(stdout$upper_2.5ci,digits = 2),"]")
#
# table2[,2]= paste(round(stdout$est,digits = 2),paste0("
[",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(
"",
rmzero(round(stdout$est[!noconf], digits = 2)),
"
[",
rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
" ",
rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
"]"
)
if (html == F) {
table2[!noconf, i + 1] = ifelse(
stdout$pval[!noconf] < .005,
yes = paste0(
rmzero(round(stdout$est[!noconf], digits = 2)),
"
[",
rmzero(round(stdout$lower_2.5ci[!noconf], digits = 2)),
" ",
rmzero(round(stdout$upper_2.5ci[!noconf], digits = 2)),
"]"
),
no = ""
)
}
table2[noconf, i + 1] = paste0(
"",
rmzero(round(stdout$est[noconf], digits = 2)),
""
)
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], "
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], "
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("π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(
"",
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)),
"]"
)}
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(
# "",
# round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2),
# ""
# )
# 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(
"",
ifelse(stdout$pval[noconf & stdout$param == "FUNC"] < .005, yes = "", no = ""),
rmzero(round(stdout$est[noconf & stdout$param == "FUNC"], digits = 2)),
"",
ifelse(stdouts1$pval[noconf & stdouts1$param == "FUNC"] < .005, yes = "", no = ""),"[",
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 = "")),
"")
table2[noconf & stdout$param == "DEP", 2] = paste0(
"",
ifelse(stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "", no = ""),
rmzero(round(stdouts1$est[noconf & stdout$param == "DEP"], digits = 2)),
"",
ifelse(stdout$pval[noconf & stdout$param == "DEP"] < .005, yes = "", 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 = "")),
""
)}
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("",round(stdout$est,digits = 2),"
[",round(stdout$lower_2.5ci,digits = 2)," ",round(stdout$upper_2.5ci,digits = 2),"]")
#
# table2[,2]= paste(round(stdout$est,digits = 2),paste0("
[",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(
"",
rmzero(round(rowMeans(cbind(stdout$est,stdouts1$est))[!noconf], digits = 2)),
"
[",
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)),
"]"
)}
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)),
"
[",
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(
# "",
# round(rowMeans(cbind(stdout$est,stdouts1$est)), digits = 2),
# ""
# )
# 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(
"",
ifelse(stdout$pval[noconf & stdout$param == "FUNC"] < .005, yes = "", no = ""),
rmzero(round(stdout$est[noconf & stdout$param == "FUNC"], digits = 2)),
"",
ifelse(stdouts1$pval[noconf & stdouts1$param == "FUNC"] < .005, yes = "", no = ""),"[",
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 = "")),
"")
table2[noconf & stdout$param == "DEP", i +1] = paste0(
"",
ifelse(stdouts1$pval[noconf & stdout$param == "DEP"] < .005, yes = "", no = ""),
rmzero(round(stdouts1$est[noconf & stdout$param == "DEP"], digits = 2)),
"",
ifelse(stdout$pval[noconf & stdout$param == "DEP"] < .005, yes = "", 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 = "")),
""
)}
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], "
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
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
or boredom","overwhelming
task","relaxation","achievement","good time
with someone"),params = c("mm","ee","mxe","me1","em1","miv","eiv"),orders = c(3,5,7,4,6,8,9), html = T,logtransparams = T)
#
# # log back transformation m_s, miv siv
# #exp(m_s) / (sqrt(exp(miv)+exp(m_s))*sqrt(exp(siv)+exp(m_s)))
#
# knitr::kable(table[[1]],escape = FALSE, align = "c",format = "html") %>%
# kable_paper() %>%
# column_spec(2,
# background = spec_color(table[[3]][,2], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
# column_spec(3,
# background = spec_color(table[[3]][,3], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
# column_spec(4,
# background = spec_color(table[[3]][,4], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
#
# column_spec(5,
# background = spec_color(table[[3]][,5], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
#
# column_spec(6,
# background = spec_color(table[[3]][,6], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
#
# column_spec(7,
# background = spec_color(table[[3]][,7], scale_from = c(-.8,.8), palette =colorRampPalette(c("#fde725","white","#21918c"))(256)))%>%
#
# add_header_above(c(" ", "Within-Person-Dynamics" = 6))
#
#
## Table with baseline PF, DEP and events with two averaged S1 mods
table = table2func(mplodelnames = c("ModNoOutS1PFsit_ncnf.out","ModNoOutS1PFsit_nemp.out","ModNoOutS1PFsit_novl.out","ModNoOutS1PFsit_prel.out","ModNoOutS1PFsit_psuc.out","ModNoOutS1PFsit_pwo.out"), adds1models = c("ModNoOutS1DEPsit_ncnf.out","ModNoOutS1DEPsit_nemp.out","ModNoOutS1DEPsit_novl.out","ModNoOutS1DEPsit_prel.out","ModNoOutS1DEPsit_psuc.out","ModNoOutS1DEPsit_pwo.out"), tablenames = c("conflict","emptiness
or boredom","overwhelming
task","relaxation","achievement","good time
with someone"),
params = c("PF~>mm", "DEP~>mm","PF~>ss", "DEP~>ss","PF~>mxe", "DEP~>mxe" , "PF~>me1","DEP~>me1","PF~>em1", "DEP~>em1", "PF~>miv" , "DEP~>miv","PF~>eiv" , "DEP~>eiv","PF~>∅mood" , "DEP~>∅mood" ,"PF~>∅event" , "DEP~>∅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
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","∅mood~>OPF" , "∅mood~>ODEP" ,"∅event~>OPF" , "∅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",
"∅mood~>PF", "∅mood~>DEP",
"∅event~>PF", "∅event~>DEP",
"PF~>OPF¹", "DEP~>OPF¹",
"PF~>ODEP¹", "DEP~>ODEP¹"
)
colnames = c("conflict
N=930","emptiness or boredom
N = 1103","overwhelming task
N = 1155","relaxation
N = 992","achievement
N = 963 ","good time with someone
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(
"",
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)),
""
)
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(
"",
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)),
""
)
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(
'',
#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 = "", 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)),
# "]
",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,3])$p < .005, yes = "", no = ""),
rmzero(round(pfcovrs[[y]]$cormat[1,3], digits = 2)),
"",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,3])$p < .005, yes = "", no = ""),
"[",
rmzero(round(depcovrs[[y]]$cormat[1,3], digits = 2)),
"] | ",
# also put squared semipartial correlations
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,3])$p < .005, yes = "", no = ""),
rmzero(round(pfsqspcors[[y]]$OFUNC[1,2], digits = 3)),
"",
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,3])$p < .005, yes = "", no = ""),
"[" ,
# also put squared semipartial correlations
rmzero(round(depsqspcors[[y]]$OFUNC[1,2], digits = 3)),
"]",
""
)
table4[grep("DEP~>OPF", rownamess),y] =
paste0(
'',
#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 = "", 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)),
# "]
",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,3])$p < .005, yes = "", no = ""),
rmzero(round(depcovrs[[y]]$cormat[2,3], digits = 2)),
"",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,3])$p < .005, yes = "", no = ""),
"[",
rmzero(round(pfcovrs[[y]]$cormat[2,3], digits = 2)),
"] | ",
# also put squared semipartial correlations
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,3])$p < .005, yes = "", no = ""),
rmzero(round(depsqspcors[[y]]$OFUNC[1,3], digits = 3)),
"",
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,3])$p < .005, yes = "", no = ""),
"[" ,
# also put squared semipartial correlations
rmzero(round(pfsqspcors[[y]]$OFUNC[1,3], digits = 3)),
"]",
""
)
# 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(
'',
#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 = "", 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)),
# "]
",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,4])$p < .005, yes = "", no = ""),
rmzero(round(pfcovrs[[y]]$cormat[1,4], digits = 2)),
"",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,4])$p < .005, yes = "", no = ""),
"[",
rmzero(round(depcovrs[[y]]$cormat[1,4], digits = 2)),
"] | ",
# also put squared semipartial correlations
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[1,4])$p < .005, yes = "", no = ""),
rmzero(round(pfsqspcors[[y]]$ODEP[1,2], digits = 3)),
"",
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[1,4])$p < .005, yes = "", no = ""),
"[" ,
# also put squared semipartial correlations
rmzero(round(depsqspcors[[y]]$ODEP[1,2], digits = 3)),
"]",
""
)
table4[grep("DEP~>ODEP", rownamess),y] =
paste0(
'',
#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 = "", 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)),
# "]
",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,4])$p < .005, yes = "", no = ""),
rmzero(round(depcovrs[[y]]$cormat[2,4], digits = 2)),
"",
#bold if significant correlation
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,4])$p < .005, yes = "", no = ""),
"[",
rmzero(round(pfcovrs[[y]]$cormat[2,4], digits = 2)),
"] | ",
# also put squared semipartial correlations
ifelse(psych::r.test(Ns[y],depcovrs[[y]]$cormat[2,4])$p < .005, yes = "", no = ""),
rmzero(round(depsqspcors[[y]]$ODEP[1,3], digits = 3)),
"",
ifelse(psych::r.test(Ns[y],pfcovrs[[y]]$cormat[2,4])$p < .005, yes = "", no = ""),
"[" ,
# also put squared semipartial correlations
rmzero(round(pfsqspcors[[y]]$ODEP[1,3], digits = 3)),
"]",
""
)
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",])
```