-------------------------------------------------------------------------------------------------------------------------------------------------
title:  "Exploring Ways of Measuring Childhood Adversity: A Comparison of Presence vs. Absence, Severity and Frequency Indicators "
author: "P. Schlechter, J. Fritz, & P. O. Wilkinson"
data <- read_csv("") #include your path here
-------------------------------------------------------------------------------------------------------------------------------------------------
  
  
 
   
-------------------------------------------------------------------------------------------------------------------------------------------------
####Table of contents:
#0. Collapse variables
#1. descriptive statistics
#2. exploratory factor analyses
#3. Confirmatory factor analyses
#4. IRT-approach
#5. Network models
#6. Validation
#7. supplemenatary analyses
-------------------------------------------------------------------------------------------------------------------------------------------------

  
    

-------------------------------------------------------------------------------------------------------------------------------------------------
####loading required packages:
library (psych)
library (car)
library (lavaan)
library (sirt)
library(readr)
library(dplyr)
library(coin)
library("mirt")
library("ltm")
library("pastecs")
library("qgraph")
library("bootnet")
library("mgm")
-------------------------------------------------------------------------------------------------------------------------------------------------
  
  
  
  
-------------------------------------------------------------------------------------------------------------------------------------------------
#### 0. collapse variables
###Note: to replicate our initial results skip this step 
  #and run the analysis with the original variables
##a) severity
data$yts_1a_full <- ifelse (data$yts_1a_full == 0, 0, 
                    ifelse (data$yts_1a_full == 1 | data$yts_1a_full == 2 | 
                            data$yts_1a_full == 3 | data$yts_1a_full == 4, 1, 2))

data$yts_2a_full <- ifelse (data$yts_2a_full == 0, 0, 
                            ifelse (data$yts_2a_full == 1 | data$yts_2a_full == 2 | 
                                      data$yts_2a_full == 3 | data$yts_2a_full == 4, 1, 2))

data$yts_3a_full <- ifelse (data$yts_3a_full == 0, 0, 
                            ifelse (data$yts_3a_full == 1 | data$yts_3a_full == 2 | 
                                      data$yts_3a_full == 3 | data$yts_3a_full == 4, 1, 2))

data$yts_4a_full <- ifelse (data$yts_4a_full == 0, 0, 
                            ifelse (data$yts_4a_full == 1 | data$yts_4a_full == 2 | 
                                      data$yts_4a_full == 3 | data$yts_4a_full == 4, 1, 2))

data$yts_5a_full <- ifelse (data$yts_5a_full == 0, 0, 
                            ifelse (data$yts_5a_full == 1 | data$yts_5a_full == 2 | 
                                      data$yts_5a_full == 3 | data$yts_5a_full == 4, 1, 2))

data$yts_6a_full <- ifelse (data$yts_6a_full == 0, 0, 
                            ifelse (data$yts_6a_full == 1 | data$yts_6a_full == 2 | 
                                      data$yts_6a_full == 3 | data$yts_6a_full == 4, 1, 2))

data$yts_7a_full <- ifelse (data$yts_7a_full == 0, 0, 
                            ifelse (data$yts_7a_full == 1 | data$yts_7a_full == 2 | 
                                      data$yts_7a_full == 3 | data$yts_7a_full == 4, 1, 2))

data$yts_8a_full <- ifelse (data$yts_8a_full == 0, 0, 
                            ifelse (data$yts_8a_full == 1 | data$yts_8a_full == 2 | 
                                      data$yts_8a_full == 3 | data$yts_8a_full == 4, 1, 2))

data$yts_9a_full <- ifelse (data$yts_9a_full == 0, 0, 
                            ifelse (data$yts_9a_full == 1 | data$yts_9a_full == 2 | 
                                      data$yts_9a_full == 3 | data$yts_9a_full == 4, 1, 2))

data$yts_10a_full <- ifelse (data$yts_10a_full == 0, 0, 
                            ifelse (data$yts_10a_full == 1 | data$yts_10a_full == 2 | 
                                      data$yts_10a_full == 3 | data$yts_10a_full == 4, 1, 2))

data$yts_11a_full <- ifelse (data$yts_11a_full == 0, 0, 
                             ifelse (data$yts_11a_full == 1 | data$yts_11a_full == 2 | 
                                       data$yts_11a_full == 3 | data$yts_11a_full == 4, 1, 2))

data$yts_13a_full <- ifelse (data$yts_13a_full == 0, 0, 
                             ifelse (data$yts_13a_full == 1 | data$yts_13a_full == 2 | 
                                       data$yts_13a_full == 3 | data$yts_13a_full == 4, 1, 2))

##b) frequency 
data$yts_1b_full <- ifelse (data$yts_1b_full == 0, 0, 
                            ifelse (data$yts_1b_full == 1 | data$yts_1b_full == 2 | 
                                      data$yts_1b_full == 3 | data$yts_1b_full == 4, 1, 2))

data$yts_2b_full <- ifelse (data$yts_2b_full == 0, 0, 
                            ifelse (data$yts_2b_full == 1 | data$yts_2b_full == 2 | 
                                      data$yts_2b_full == 3 | data$yts_2b_full == 4, 1, 2))

data$yts_3b_full <- ifelse (data$yts_3b_full == 0, 0, 
                            ifelse (data$yts_3b_full == 1 | data$yts_3b_full == 2 | 
                                      data$yts_3b_full == 3 | data$yts_3b_full == 4, 1, 2))

data$yts_4b_full <- ifelse (data$yts_4b_full == 0, 0, 
                            ifelse (data$yts_4b_full == 1 | data$yts_4b_full == 2 | 
                                      data$yts_4b_full == 3 | data$yts_4b_full == 4, 1, 2))

data$yts_5b_full <- ifelse (data$yts_5b_full == 0, 0, 
                            ifelse (data$yts_5b_full == 1 | data$yts_5b_full == 2 | 
                                      data$yts_5b_full == 3 | data$yts_5b_full == 4, 1, 2))


data$yts_6b_full <- ifelse (data$yts_6b_full == 0, 0, 
                            ifelse (data$yts_6b_full == 1 | data$yts_6b_full == 2 | 
                                      data$yts_6b_full == 3 | data$yts_6b_full == 4, 1, 2))

data$yts_7b_full <- ifelse (data$yts_7b_full == 0, 0, 
                            ifelse (data$yts_7b_full == 1 | data$yts_7b_full == 2 | 
                                      data$yts_7b_full == 3 | data$yts_7b_full == 4, 1, 2))

data$yts_8b_full <- ifelse (data$yts_8b_full == 0, 0, 
                            ifelse (data$yts_8b_full == 1 | data$yts_8b_full == 2 | 
                                      data$yts_8b_full == 3 | data$yts_8b_full == 4, 1, 2))

data$yts_9b_full <- ifelse (data$yts_9b_full == 0, 0, 
                            ifelse (data$yts_9b_full == 1 | data$yts_9b_full == 2 | 
                                      data$yts_9b_full == 3 | data$yts_9b_full == 4, 1, 2))

data$yts_10b_full <- ifelse (data$yts_10b_full == 0, 0, 
                             ifelse (data$yts_10b_full == 1 | data$yts_10b_full == 2 | 
                                       data$yts_10b_full == 3 | data$yts_10b_full == 4, 1, 2))

data$yts_11b_full <- ifelse (data$yts_11b_full == 0, 0, 
                             ifelse (data$yts_11b_full == 1 | data$yts_11b_full == 2 | 
                                       data$yts_11b_full == 3 | data$yts_11b_full == 4, 1, 2))

data$yts_13b_full <- ifelse (data$yts_13b_full == 0, 0, 
                             ifelse (data$yts_13b_full == 1 | data$yts_13b_full == 2 | 
                                       data$yts_13b_full == 3 | data$yts_13b_full == 4, 1, 2))
--------------------------------------------------------------------------------------------------------------
  
  
  
  
--------------------------------------------------------------------------------------------------------------
##### 1. descriptive statistics

### a)yes/no, see table 2 
dataYN <- data.frame(data$yts_1, data$yts_2, data$yts_3, data$yts_4, 
                     data$yts_5, data$yts_6, data$yts_7, data$yts_8, 
                     data$yts_9, data$yts_10,data$yts_11,data$yts_12, 
                     data$yts_13)

describe (dataYN) ####se also supplementary material 2

100 * prop.table(table(dataYN$data.yts_1))
table(dataYN$data.yts_1)

100 * prop.table(table(dataYN$data.yts_2))
table(dataYN$data.yts_2)

100 * prop.table(table(dataYN$data.yts_3))
table(dataYN$data.yts_3)

100 * prop.table(table(dataYN$data.yts_4))
table(dataYN$data.yts_4)

100 * prop.table(table(dataYN$data.yts_5))
table(dataYN$data.yts_5)

100 * prop.table(table(dataYN$data.yts_6))
table(dataYN$data.yts_6)

100 * prop.table(table(dataYN$data.yts_7))
table(dataYN$data.yts_7)

100 * prop.table(table(dataYN$data.yts_8))
table(dataYN$data.yts_8)

100 * prop.table(table(dataYN$data.yts_9))
table(dataYN$data.yts_9)

100 * prop.table(table(dataYN$data.yts_10))
table(dataYN$data.yts_10)

100 * prop.table(table(dataYN$data.yts_11))
table(dataYN$data.yts_11)

100 * prop.table(table(dataYN$data.yts_12))
table(dataYN$data.yts_12)

100 * prop.table(table(dataYN$data.yts_13))
table(dataYN$data.yts_13)


##b) Severity
dataseverity <- data.frame(data$yts_1a_full, data$yts_2a_full, data$yts_3a_full, data$yts_4a_full, 
                           data$yts_5a_full, data$yts_6a_full, data$yts_7a_full, data$yts_8a_full, 
                           data$yts_9a_full, data$yts_10a_full,data$yts_11a_full,data$yts_12a_full, 
                           data$yts_13a_full)

describe (dataseverity)
quantile (data$yts_1a_full, na.rm = T)
quantile (data$yts_2a_full, na.rm = T)
quantile (data$yts_3a_full, na.rm = T)
quantile (data$yts_4a_full, na.rm = T)
quantile (data$yts_5a_full, na.rm = T)
quantile (data$yts_7a_full, na.rm = T)
quantile (data$yts_8a_full, na.rm = T)
quantile (data$yts_9a_full, na.rm = T)
quantile (data$yts_10a_full, na.rm = T)
quantile (data$yts_11a_full, na.rm = T)
quantile (data$yts_12a_full, na.rm = T)
quantile (data$yts_13a_full, na.rm = T)

##c) datafrequency
datafrequency <- data.frame(data$yts_1b_full, data$yts_2b_full, data$yts_3b_full, data$yts_4b_full, 
                            data$yts_5b_full, data$yts_6b_full, data$yts_7b_full, data$yts_8b_full, 
                            data$yts_9b_full, data$yts_10b_full,data$yts_11b_full,data$yts_12b_full,
                            data$yts_13b_full)

describe(datafrequency)

quantile (data$yts_1b_full, na.rm = T)
quantile (data$yts_2b_full, na.rm = T)
quantile (data$yts_3b_full, na.rm = T)
quantile (data$yts_4b_full, na.rm = T)
quantile (data$yts_5b_full, na.rm = T)
quantile (data$yts_7b_full, na.rm = T)
quantile (data$yts_8b_full, na.rm = T)
quantile (data$yts_9b_full, na.rm = T)
quantile (data$yts_10b_full, na.rm = T)
quantile (data$yts_11b_full, na.rm = T)
quantile (data$yts_12b_full, na.rm = T)
quantile (data$yts_13b_full, na.rm = T)


##d) agebin, see table 3
##data agebin: C__1 = birth to 12, c__2 = 12 to 18, c___3 = 18 +

data$yts_1c___2[which(data$yts_1c___2 ==1)] <- 2
data$yts_1c___3[which(data$yts_1c___3 ==1)] <- 3
data$yts_1c <- data$yts_1c___1 + data$yts_1c___2 + data$yts_1c___3
data$yts_1c[which(data$yts_1c ==0)] <- NA


data$yts_2c___2[which(data$yts_2c___2 ==1)] <- 2
data$yts_2c___3[which(data$yts_2c___3 ==1)] <- 3
data$yts_2c <- data$yts_2c___1 + data$yts_2c___2 + data$yts_2c___3
data$yts_2c[which(data$yts_2c ==0)] <- NA

data$yts_3c___2[which(data$yts_3c___2 ==1)] <- 2
data$yts_3c___3[which(data$yts_3c___3 ==1)] <- 3
data$yts_3c <- data$yts_3c___1 + data$yts_3c___2 + data$yts_3c___3
data$yts_3c[which(data$yts_3c ==0)] <- NA

data$yts_4c___2[which(data$yts_4c___2 ==1)] <- 2
data$yts_4c___3[which(data$yts_4c___3 ==1)] <- 3
data$yts_4c <- data$yts_4c___1 + data$yts_4c___2 + data$yts_4c___3
data$yts_4c[which(data$yts_4c ==0)] <- NA

data$yts_5c___2[which(data$yts_5c___2 ==1)] <- 2
data$yts_5c___3[which(data$yts_5c___3 ==1)] <- 3
data$yts_5c <- data$yts_5c___1 + data$yts_5c___2 + data$yts_5c___3
data$yts_5c[which(data$yts_5c ==0)] <- NA

data$yts_6c___2[which(data$yts_6c___2 ==1)] <- 2
data$yts_6c___3[which(data$yts_6c___3 ==1)] <- 3
data$yts_6c <- data$yts_6c___1 + data$yts_6c___2 + data$yts_6c___3
data$yts_6c[which(data$yts_6c ==0)] <- NA

data$yts_7c___2[which(data$yts_7c___2 ==1)] <- 2
data$yts_7c___3[which(data$yts_7c___3 ==1)] <- 3
data$yts_7c <- data$yts_7c___1 + data$yts_7c___2 + data$yts_7c___3
data$yts_7c[which(data$yts_7c ==0)] <- NA

data$yts_8c___2[which(data$yts_8c___2 ==1)] <- 2
data$yts_8c___3[which(data$yts_8c___3 ==1)] <- 3
data$yts_8c <- data$yts_8c___1 + data$yts_8c___2 + data$yts_8c___3
data$yts_8c[which(data$yts_8c ==0)] <- NA
data$yts_8c

data$yts_9c___2[which(data$yts_9c___2 ==1)] <- 2
data$yts_9c___3[which(data$yts_9c___3 ==1)] <- 3
data$yts_9c <- data$yts_9c___1 + data$yts_9c___2 + data$yts_9c___3
data$yts_9c[which(data$yts_9c ==0)] <- NA

data$yts_10c___2[which(data$yts_10c___2 ==1)] <- 2
data$yts_10c___3[which(data$yts_10c___3 ==1)] <- 3
data$yts_10c <- data$yts_10c___1 + data$yts_10c___2 + data$yts_10c___3
data$yts_10c[which(data$yts_10c ==0)] <- NA

data$yts_11c___2[which(data$yts_11c___2 ==1)] <- 2
data$yts_11c___3[which(data$yts_11c___3 ==1)] <- 3
data$yts_11c <- data$yts_11c___1 + data$yts_11c___2 + data$yts_11c___3
data$yts_11c[which(data$yts_11c ==0)] <- NA

data$yts_12c___2[which(data$yts_12c___2 ==1)] <- 2
data$yts_12c___3[which(data$yts_12c___3 ==1)] <- 3
data$yts_12c <- data$yts_12c___1 + data$yts_12c___2 + data$yts_12c___3
data$yts_12c[which(data$yts_12c ==0)] <- NA

data$yts_13c___2[which(data$yts_13c___2 ==1)] <- 2
data$yts_13c___3[which(data$yts_13c___3 ==1)] <- 3
data$yts_13c <- data$yts_13c___1 + data$yts_13c___2 + data$yts_13c___3
data$yts_13c[which(data$yts_13c ==0)] <- NA


agebin <- data.frame (data$yts_1c, data$yts_2c, data$yts_3c, data$yts_4c, data$yts_5c, data$yts_6c,
                      data$yts_7c, data$yts_8c, data$yts_9c, data$yts_10c, data$yts_11c, data$yts_12c,data$yts_13c)

describe (agebin)


agebin$data.yts_1agebin <- ifelse (agebin$data.yts_1c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_1c == 2, "12to18",
                                           ifelse (agebin$data.yts_1c == 3, "18+", 
                                                   ifelse (agebin$data.yts_1c >= 4 , "several", NA))))

agebin$data.yts_2agebin <- ifelse (agebin$data.yts_2c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_2c == 2, "12to18",
                                           ifelse (agebin$data.yts_2c == 3, "18+", 
                                                   ifelse (agebin$data.yts_2c >= 4 , "several", NA))))

agebin$data.yts_3agebin <- ifelse (agebin$data.yts_3c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_3c == 2, "12to18",
                                           ifelse (agebin$data.yts_3c == 3, "18+", 
                                                   ifelse (agebin$data.yts_3c >= 4 , "several", NA))))

agebin$data.yts_4agebin <- ifelse (agebin$data.yts_4c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_4c == 2, "12to18",
                                           ifelse (agebin$data.yts_4c == 3, "18+", 
                                                   ifelse (agebin$data.yts_4c >= 4 , "several", NA))))

agebin$data.yts_5agebin <- ifelse (agebin$data.yts_5c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_5c == 2, "12to18",
                                           ifelse (agebin$data.yts_5c == 3, "18+", 
                                                   ifelse (agebin$data.yts_5c >= 4 , "several", NA))))

agebin$data.yts_6agebin <- ifelse (agebin$data.yts_6c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_6c == 2, "12to18",
                                           ifelse (agebin$data.yts_6c == 3, "18+", 
                                                   ifelse (agebin$data.yts_6c >= 4 , "several", NA))))

agebin$data.yts_7agebin <- ifelse (agebin$data.yts_7c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_7c == 2, "12to18",
                                           ifelse (agebin$data.yts_7c == 3, "18+", 
                                                   ifelse (agebin$data.yts_7c >= 4 , "several", NA))))

agebin$data.yts_8agebin <- ifelse (agebin$data.yts_8c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_8c == 2, "12to18",
                                           ifelse (agebin$data.yts_8c == 3, "18+", 
                                                   ifelse (agebin$data.yts_8c >= 4 , "several", NA))))

agebin$data.yts_9agebin <- ifelse (agebin$data.yts_9c == 1, "birthto12", 
                                   ifelse (agebin$data.yts_9c == 2, "12to18",
                                           ifelse (agebin$data.yts_9c == 3, "18+", 
                                                   ifelse (agebin$data.yts_9c >= 4 , "several", NA))))

agebin$data.yts_10agebin <- ifelse (agebin$data.yts_10c == 1, "birthto12", 
                                    ifelse (agebin$data.yts_10c == 2, "12to18",
                                            ifelse (agebin$data.yts_10c == 3, "18+", 
                                                    ifelse (agebin$data.yts_10c >= 4 , "several", NA))))

agebin$data.yts_11agebin <- ifelse (agebin$data.yts_11c == 1, "birthto12", 
                                    ifelse (agebin$data.yts_11c == 2, "12to18",
                                            ifelse (agebin$data.yts_11c == 3, "18+", 
                                                    ifelse (agebin$data.yts_11c >= 4 , "several", NA))))

agebin$data.yts_12agebin <- ifelse (agebin$data.yts_12c == 1, "birthto12", 
                                    ifelse (agebin$data.yts_12c == 2, "12to18",
                                            ifelse (agebin$data.yts_12c == 3, "18+", 
                                                    ifelse (agebin$data.yts_12c >= 4 , "several", NA))))

agebin$data.yts_13agebin <- ifelse (agebin$data.yts_13c == 1, "birthto12", 
                                    ifelse (agebin$data.yts_13c == 2, "12to18",
                                            ifelse (agebin$data.yts_13c == 3, "18+", 
                                                    ifelse (agebin$data.yts_13c >= 4 , "several", NA))))

table(agebin$data.yts_1agebin)
100 * prop.table(table(agebin$data.yts_1agebin))

table(agebin$data.yts_2agebin)
100 * prop.table(table(agebin$data.yts_2agebin))

table(agebin$data.yts_3agebin)
100 * prop.table(table(agebin$data.yts_3agebin))

table(agebin$data.yts_4agebin)
100 * prop.table(table(agebin$data.yts_4agebin))

table(agebin$data.yts_5agebin)
100 * prop.table(table(agebin$data.yts_5agebin))

table(agebin$data.yts_6agebin)
100 * prop.table(table(agebin$data.yts_6agebin))

table(agebin$data.yts_7agebin)
100 * prop.table(table(agebin$data.yts_7agebin))

table(agebin$data.yts_8agebin)
100 * prop.table(table(agebin$data.yts_8agebin))

table(agebin$data.yts_9agebin)
100 * prop.table(table(agebin$data.yts_9agebin))

table(agebin$data.yts_10agebin)
100 * prop.table(table(agebin$data.yts_10agebin))

table(agebin$data.yts_11agebin)
100 * prop.table(table(agebin$data.yts_11agebin))

table(agebin$data.yts_12agebin)
100 * prop.table(table(agebin$data.yts_12agebin))

table(agebin$data.yts_13agebin)
100 * prop.table(table(agebin$data.yts_13agebin))
--------------------------------------------------------------------------------------------------------------
  
  
  
  
--------------------------------------------------------------------------------------------------------------
##### 2. exploratory factor analyses
  
##buidling correlation matrices
matrix1 <- tetrachoric2(dataYN)$rho#####correlation matrix for Yes/no Items
matrix1
matrixtest <- polychoric2(dataYN)$rho #####double-check whether polychoric yields the same results
matrixtest

##correlation matrix for severity
R <- polychoric2(dataseverity)$rho
R

##correlation matrix for frequency
R2 <- polychoric2(datafrequency)$rho
R2

##Principal compomnents analysis####

nfact <-13 ###at first as many factors as Items
nfact2 <- 13
nfact3 <- 13

factors.pc <- principal(matrixtest, nfact)
factors.pc1 <- principal(R, nfact2)
factors.pc2 <- principal(R2, nfact3)

## Kaiser-Criterion & Scree Plot for yes/no Items
plot(c(1:nfact), factors.pc$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(matrix1,n.obs=nrow(matrixtest),fa="fa",main="Parallel Analysis Scree Plots")

## Kaiser-Criterion & Scree Plot for severity
plot(c(1:nfact2), factors.pc1$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(R,n.obs=nrow(R),fa="fa",main="Parallel Analysis Scree Plots")

## Kaiser-Criterion & Scree Plot for frequency
plot(c(1:nfact3), factors.pc2$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(R2,n.obs=nrow(R2),fa="fa",main="Parallel Analysis Scree Plots")


####factor loadings####
##depending on the suggested solution which was 1 factor
faYN <- fa (matrixtest, 1, fm = "wls", rotate = "promax") ###yes/no Items
print (faYN, digits = 2)

##severity
faseverity <- fa (R, 1, fm = "wls", rotate = "promax")
print (faseverity, digits = 2)
faseverity$e.values

##frequency
fafrequency <- fa (R2, 1, fm = "wls", rotate = "promax")
print (fafrequency, digits = 2)
fafrequency$e.values


##EFA factor analysis without item 6 and 12, see table 4
dataYN <- data.frame(data$yts_1, data$yts_2, data$yts_3, data$yts_4, 
                     data$yts_5, data$yts_7, data$yts_8, 
                     data$yts_9, data$yts_10,data$yts_11, 
                     data$yts_13) ####without item 6 and 12

dataseverity <- data.frame( data$yts_1a_full,data$yts_2a_full, data$yts_3a_full, data$yts_4a_full, 
                           data$yts_5a_full,  data$yts_7a_full, data$yts_8a_full, 
                           data$yts_9a_full, data$yts_10a_full,data$yts_11a_full, 
                           data$yts_13a_full) ####without item 6 and 12

datafrequency <- data.frame( data$yts_1b_full,data$yts_2b_full, data$yts_3b_full, data$yts_4b_full, 
                            data$yts_5b_full,  data$yts_7b_full, data$yts_8b_full, 
                            data$yts_9b_full, data$yts_10b_full,data$yts_11b_full,
                            data$yts_13b_full) ####without item 6 and 12



##buidling correlation matrices
matrix1 <- tetrachoric2(dataYN)$rho#####correlation matrix for Yes/no Items
matrix1
matrixtest <- polychoric2(dataYN)$rho #####double-check whether polychoric yields the same results
matrixtest

##correlation matrix for severity
R <- polychoric2(dataseverity)$rho
R

##correlation matrix for frequency
R2 <- polychoric2(datafrequency)$rho
R2


#### Principal compomnents analysis####

nfact <- 11 ###down to 11 because we dropped item 6 and 12 
nfact2 <- 11
nfact3 <- 11

factors.pc <- principal(matrix1, nfact)
factors.pc1 <- principal(R, nfact2)
factors.pc2 <- principal(R2, nfact3)

## Kaiser-Criterion & Scree Plot for yes/no Items
plot(c(1:nfact), factors.pc$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(matrix1,n.obs=nrow(matrix1),fa="fa",main="Parallel Analysis Scree Plots")

## Kaiser-Criterion & Scree Plot for severity
plot(c(1:nfact2), factors.pc1$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(R,n.obs=nrow(R),fa="fa",main="Parallel Analysis Scree Plots")

## Kaiser-Criterion & Scree Plot for frequency
plot(c(1:nfact3), factors.pc2$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(R2,n.obs=nrow(R2),fa="fa",main="Parallel Analysis Scree Plots")


####factor loadings####
##depending on the suggested solution which was 1 factor
faYN <- fa (matrix1, 1, fm = "wls", rotate = "promax") ###yes/no Items
print (faYN, digits = 2)

##severity
faseverity <- fa (R, 1, fm = "wls", rotate = "promax")
print (faseverity, digits = 2)

##frequency
fafrequency <- fa (R2, 1, fm = "wls", rotate = "promax")
print (fafrequency, digits = 2)
--------------------------------------------------------------------------------------------------------------
  


  
--------------------------------------------------------------------------------------------------------------
#### 3. Confirmatory factor analyses, see table 5
##a) yes/no
fs <- data.frame (data$id, data$yts_1, data$yts_2, data$yts_3, data$yts_4, 
                    data$yts_5, data$yts_7, data$yts_8, 
                    data$yts_9, data$yts_10, data$yts_11,  data$yts_13)

fs <- na.omit (fs)


m.cong <- 'Trauma =~  data.yts_1 + data.yts_2 + data.yts_3 + data.yts_4 + data.yts_5 + data.yts_7 +
                      data.yts_8 + data.yts_9 + data.yts_10 + data.yts_11+  data.yts_13
'


#running a cfa
r.cong <- cfa(m.cong, data=fs, meanstructure=T, std.lv=T, mimic = "Mplus",
              ordered = c( "data.yts_1", "data.yts_2","data.yts_3", "data.yts_4", "data.yts_5", 
                           "data.yts_7", "data.yts_8", "data.yts_9","data.yts_10","data.yts_11" 
                           ,"data.yts_13"))

##results
summary(r.cong, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr"))
fitMeasures(r.cong,c('cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

####Omega H 
parameterEstimates(r.cong)
loading2 <- parameterEstimates(r.cong)$est[1:11] 
Errors2 <- parameterEstimates(r.cong)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))

modindices(r.cong, sort. = T)

####with modification indices
m.cong <- 'Trauma =~  a * data.yts_1 + b* data.yts_2 +  c * data.yts_3 +  d*data.yts_4 +   e*data.yts_5 + f*data.yts_7 +
                      g * data.yts_8 + h*  data.yts_9 + i * data.yts_10 + J* data.yts_11+  k* data.yts_13

data.yts_1 ~~ data.yts_7       #####included after the first CFA
data.yts_2 ~~ data.yts_9

a > 0
b > 0
c > 0
d > 0
e > 0
f > 0
g > 0
h > 0
i > 0
J > 0
k > 0

'

#running a cfa
r.cong <- cfa(m.cong, data=fs, meanstructure=T, mimic = "Mplus", std.lv=T, 
              ordered = c( "data.yts_1", "data.yts_2","data.yts_3", "data.yts_4", "data.yts_5", 
                           "data.yts_7", "data.yts_8", "data.yts_9","data.yts_10","data.yts_11", 
                           "data.yts_13"))


##results
summary(r.cong, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr"))
fitMeasures(r.cong,c('cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

modindices(r.cong, sort. = T)

####get factorscores and include them in origingal data
fs$factorscore1 <- predict(r.cong)
library("dplyr")
names (fs) [1] <- "id"
data <- left_join(data, fs,  by = "id")
data$factorscore1

####Omega H 
parameterEstimates(r.cong)
loading2 <- parameterEstimates(r.cong)$est[1:11] 
Errors2 <- parameterEstimates(r.cong)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))

##comparing with essent. tau model 
m.esstau <- 'Trauma =~ A * data.yts_1 +  A * data.yts_2   + A * data.yts_3  + A * data.yts_4 + 
                       A * data.yts_5 +  A * data.yts_7   + A * data.yts_8  +
                       A * data.yts_9 +  A * data.yts_10  + A * data.yts_11 +  
                       A * data.yts_13

data.yts_1 ~~ data.yts_7      
data.yts_2 ~~ data.yts_9

'

r.esstau <- cfa(m.esstau, data=fs, meanstructure=TRUE,std.lv=TRUE, mimic = "Mplus", 
                ordered = c("data.yts_1", "data.yts_2","data.yts_3", "data.yts_4", "data.yts_5","data.yts_7",
                            "data.yts_8", "data.yts_9","data.yts_10","data.yts_11", "data.yts_13"))

summary(r.esstau, fit.measures=TRUE) 


fitMeasures(r.esstau,c('cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

#### Modellcomp. ####
anova(r.cong, r.esstau)

###################################b) severity###################################
fs2 <- data.frame(data$id, data$yts_1a_full, data$yts_2a_full, data$yts_3a_full, data$yts_4a_full, 
                  data$yts_5a_full, data$yts_7a_full, data$yts_8a_full, 
                  data$yts_9a_full, data$yts_10a_full,data$yts_11a_full, 
                  data$yts_13a_full)

fs2 <- na.omit (fs2)

#defining model structure
m.cong1 <- 'Trauma =~  data.yts_1a_full + data.yts_2a_full + data.yts_3a_full + data.yts_4a_full + 
                       data.yts_5a_full + data.yts_7a_full + data.yts_8a_full + data.yts_9a_full +
                       data.yts_10a_full + data.yts_11a_full  + data.yts_13a_full
'

r.cong1 <- cfa(m.cong1, data=fs2, meanstructure=T, std.lv=T, mimic = "Mplus", 
               ordered =c("data.yts_1a_full","data.yts_2a_full","data.yts_3a_full", "data.yts_4a_full",
                          "data.yts_5a_full", "data.yts_7a_full", "data.yts_8a_full", "data.yts_9a_full",  
                          "data.yts_10a_full","data.yts_11a_full", "data.yts_13a_full"))

##results
summary(r.cong1, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong1,c("chisq","df","pvalue","cfi","tli","rmsea",
                      "wrmr", "cfi.scaled","tli.scaled","rmsea.scaled","rmsea.ci.lower.scaled","rmsea.ci.upper.scaled"))

modindices(r.cong1, sort. = T)


#######Omega H 
parameterEstimates(r.cong1)
loading2 <- parameterEstimates(r.cong1)$est[1:11] 
Errors2 <- parameterEstimates(r.cong1)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))


#####with modifications

#defining model structure
m.cong1 <- 'Trauma =~  a * data.yts_1a_full +  b * data.yts_2a_full + c * data.yts_3a_full + d * data.yts_4a_full + 
                       e * data.yts_5a_full +  f * data.yts_7a_full + g * data.yts_8a_full + h * data.yts_9a_full +
                       i * data.yts_10a_full + J* data.yts_11a_full + k * data.yts_13a_full

data.yts_1a_full ~~ data.yts_7a_full ####Modifications made 
data.yts_2a_full ~~ data.yts_9a_full

a > 0
b > 0
c > 0
d > 0
e > 0
f > 0
g > 0
h > 0
i > 0
J > 0
k > 0

'

r.cong1 <- cfa(m.cong1, data=fs2, meanstructure=T, std.lv=T, mimic = "Mplus", 
               ordered =c("data.yts_1a_full","data.yts_2a_full","data.yts_3a_full", "data.yts_4a_full",
                          "data.yts_5a_full", "data.yts_7a_full", "data.yts_8a_full", "data.yts_9a_full",  
                          "data.yts_10a_full","data.yts_11a_full","data.yts_13a_full"))

##results
summary(r.cong1, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong1,c("chisq","df","pvalue","cfi","tli","rmsea",
                      "wrmr", "cfi.scaled","tli.scaled","rmsea.scaled","rmsea.ci.lower.scaled","rmsea.ci.upper.scaled"))

modindices(r.cong1, sort. = T)


#######get factorscores and include them in origingal data
fs2$factorscore2 <- predict(r.cong1)
names (fs2) [1] <- "id"
data <- left_join(data, fs2,  by = "id")
data$factorscore2

#######Omega H 
parameterEstimates(r.cong1)
loading2 <- parameterEstimates(r.cong1)$est[1:11] 
Errors2 <- parameterEstimates(r.cong1)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))


#####esstau
m.esstau1 <- 'Trauma =~  A * data.yts_1a_full + A * data.yts_2a_full + A * data.yts_3a_full + A * data.yts_4a_full + 
                         A * data.yts_5a_full + A * data.yts_7a_full + A * data.yts_8a_full + A * data.yts_9a_full +
                         A * data.yts_10a_full +A *  data.yts_11a_full + A *  data.yts_13a_full

data.yts_1a_full ~~ data.yts_7a_full
data.yts_2a_full ~~ data.yts_9a_full

'

r.esstau1 <- cfa(m.esstau1, data=fs2, meanstructure=T, std.lv=T, mimic = "Mplus", 
                 ordered =c( "data.yts_1a_full","data.yts_2a_full","data.yts_3a_full", "data.yts_4a_full",
                             "data.yts_5a_full", "data.yts_7a_full", "data.yts_8a_full", "data.yts_9a_full",  
                             "data.yts_10a_full","data.yts_11a_full","data.yts_13a_full"))

##results
summary(r.esstau1, fit.measures=TRUE, standardized = T)

fitMeasures(r.esstau1,c("chisq","df","pvalue","cfi","tli","rmsea",
                        'wrmr", cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))



#### Modellcomp. ####
anova(r.cong1, r.esstau1)


#################c) frequency##################################
########frequency###############
fs3 <- data.frame(data$id, data$yts_1b_full, data$yts_2b_full, data$yts_3b_full, data$yts_4b_full, 
                  data$yts_5b_full, data$yts_7b_full, data$yts_8b_full, data$yts_9b_full, 
                  data$yts_10b_full,data$yts_11b_full, data$yts_13b_full)

fs3 <- na.omit (fs3)


#defining model structure
m.cong2 <- 'Trauma =~ data.yts_1b_full + data.yts_2b_full + data.yts_3b_full + data.yts_4b_full +
data.yts_5b_full + data.yts_7b_full + data.yts_8b_full + data.yts_9b_full + 
data.yts_10b_full+ data.yts_11b_full +data.yts_13b_full

'


#running a cfa
r.cong2 <- cfa(m.cong2, data=fs3, meanstructure=T, std.lv=T, mimic = "Mplus",
               ordered = c("data.yts_1b_full", "data.yts_2b_full","data.yts_3b_full", "data.yts_4b_full",
                           "data.yts_5b_full", "data.yts_7b_full", "data.yts_8b_full","data.yts_9b_full",
                           "data.yts_10b_full","data.yts_11b_full","data.yts_13b_full"))


##results
summary(r.cong2, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong2,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr",
                      'cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

modindices(r.cong2, sort. = T)


######Omega H#######
parameterEstimates(r.cong2)
loading2 <- parameterEstimates(r.cong2)$est[1:11] 
Errors2 <-  parameterEstimates(r.cong2)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))

#defining model structure with modification indices
m.cong2 <- 'Trauma =~ a * data.yts_1b_full + b * data.yts_2b_full + c * data.yts_3b_full + d * data.yts_4b_full +
                      e * data.yts_5b_full + f * data.yts_7b_full + g * data.yts_8b_full + h * data.yts_9b_full + 
                      i * data.yts_10b_full+ J * data.yts_11b_full+ k *  data.yts_13b_full

data.yts_1b_full ~~ data.yts_7b_full
data.yts_2b_full ~~ data.yts_9b_full

a > 0
b > 0
c > 0
d > 0
e > 0
f > 0
g > 0
h > 0
i > 0
J > 0
k > 0

'


#running a cfa
r.cong2 <- cfa(m.cong2, data=fs3, meanstructure=T, std.lv=T, mimic = "Mplus",
               ordered = c("data.yts_1b_full", "data.yts_2b_full","data.yts_3b_full", "data.yts_4b_full",
                           "data.yts_5b_full", "data.yts_7b_full", "data.yts_8b_full","data.yts_9b_full",
                           "data.yts_10b_full","data.yts_11b_full","data.yts_13b_full"))


##results
summary(r.cong2, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong2,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr",
                      'cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

modindices(r.cong2, sort. = T)

#######get factorscores and include them in origingal data

fs3$factorscore3 <- predict(r.cong2)
names (fs3) [1] <- "id"
data <- left_join(data, fs3,  by = "id")
data$factorscore3

######Omega H#######
parameterEstimates(r.cong2)
loading2 <- parameterEstimates(r.cong2)$est[1:11] 
Errors2 <-  parameterEstimates(r.cong2)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))


####tau model###
#defining model structure with modification indices
m.esstau2 <- 'Trauma =~ A * data.yts_1b_full + A * data.yts_2b_full + A * data.yts_3b_full + A * data.yts_4b_full +
                        A * data.yts_5b_full + A * data.yts_7b_full + A * data.yts_8b_full + A * data.yts_9b_full + 
                        A * data.yts_10b_full+ A * data.yts_11b_full+ A * data.yts_13b_full

data.yts_1b_full ~~ data.yts_7b_full
data.yts_2b_full ~~ data.yts_9b_full

'


#running a cfa
r.esstau2 <- cfa(m.esstau2, data=fs3, meanstructure=T, std.lv=T, mimic = "Mplus",
                 ordered = c("data.yts_1b_full", "data.yts_2b_full","data.yts_3b_full", "data.yts_4b_full",
                             "data.yts_5b_full", "data.yts_7b_full", "data.yts_8b_full","data.yts_9b_full",
                             "data.yts_10b_full","data.yts_11b_full","data.yts_13b_full"))

##results
summary(r.esstau2, fit.measures=TRUE, standardized = T)

fitMeasures(r.esstau2,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr",
                        'cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

anova(r.cong2, r.esstau2)
--------------------------------------------------------------------------------------------------------------

  
  
  
--------------------------------------------------------------------------------------------------------------
#### 4. IRT approach
####put yes/no items in a dataframe
fs <- data.frame (data$yts_1, data$yts_2, data$yts_3, data$yts_4, 
                  data$yts_5,  data$yts_7, data$yts_8, 
                  data$yts_9, data$yts_10, data$yts_11,  data$yts_13)

fs <- na.omit (fs)

###define model, one dimension
irtmodel = ltm (fs ~ z1, IRT.param = T)
summary (irtmodel) 
coef (irtmodel)
plot (irtmodel, type = "ICC")
plot (irtmodel, type = "IIC")
plot (irtmodel, type = "IIC", items = 8)
factor.scores (irtmodel)
item.fit(irtmodel)


##b)severity 
fs2 <- data.frame(data$yts_1a_full, data$yts_2a_full, data$yts_3a_full, data$yts_4a_full, 
                  data$yts_5a_full, data$yts_7a_full, data$yts_8a_full, 
                  data$yts_9a_full, data$yts_10a_full,data$yts_11a_full, 
                  data$yts_13a_full)
fs2 <- na.omit (fs2)

polymodel1 = mirt (data=fs2, model = 1, itemtype="graded")
summary (polymodel1)
coef (polymodel1)
itemplot (polymodel1, 1, type = "trace")
itemplot (polymodel1, 4, type = "info")
plot (polymodel1, type = "trace")
plot (polymodel1,  type = "info")
plot (polymodel1)
fscores (polymodel1)
itemfit (polymodel1)


##c)frequency
fs2 <- data.frame(data$yts_1b_full, data$yts_2b_full, data$yts_3b_full, data$yts_4b_full, 
                  data$yts_5b_full, data$yts_7b_full, data$yts_8b_full, 
                  data$yts_9b_full, data$yts_10b_full,data$yts_11b_full, 
                  data$yts_13b_full)
fs2 <- na.omit (fs2)

polymodel2 = mirt (data=fs2, model = 1, itemtype = "graded")
summary (polymodel2)
coef (polymodel2)
itemplot (polymodel2, 1, type = "trace")
itemplot (polymodel2, 1, type = "info")
plot (polymodel2, type = "trace")
plot (polymodel2,  type = "info")
plot (polymodel2)
itemfit (polymodel2)  
--------------------------------------------------------------------------------------------------------------

  
  
  
--------------------------------------------------------------------------------------------------------------
#### 5. Network models
  
## adapt items for the frequency model
data$yts_9b_full_II <- ifelse (data$yts_9b_full == 0, 0, 
                            ifelse (data$yts_9b_full == 1 | data$yts_9b_full == 2 | 
                                      data$yts_9b_full == 3 | data$yts_9b_full == 4, 1, 2)); table (data$yts_9b_full_II)
data$yts_10b_full_II <- ifelse (data$yts_10b_full == 0, 0, 
                             ifelse (data$yts_10b_full == 1 | data$yts_10b_full == 2 | 
                                       data$yts_10b_full == 3 | data$yts_10b_full == 4, 1, 2));table (data$yts_10b_full_II)

data$yts_11b_full_II <- ifelse (data$yts_11b_full == 0, 0, 
                             ifelse (data$yts_11b_full == 1 | data$yts_11b_full == 2 | 
                                       data$yts_11b_full == 3 | data$yts_11b_full == 4, 1, 2));table (data$yts_11b_full_II)
Network_Data <- data

## making the data frame for the presence vs absecense network
pres_vs_abs       <- na.omit(data.frame(Network_Data$yts_1,
                                        Network_Data$yts_2,
                                        Network_Data$yts_3,
                                        Network_Data$yts_4,
                                        Network_Data$yts_5,
                                        Network_Data$yts_7,
                                        Network_Data$yts_8,
                                        Network_Data$yts_9,
                                        Network_Data$yts_10,
                                        Network_Data$yts_11,
                                        Network_Data$yts_13));nrow(pres_vs_abs) # 316

## making the data frame for the severity network
severity          <- na.omit(data.frame(Network_Data$yts_1a_full,
                                        Network_Data$yts_2a_full,
                                        Network_Data$yts_3a_full,
                                        Network_Data$yts_4a_full,
                                        Network_Data$yts_5a_full,
                                        Network_Data$yts_7a_full,
                                        Network_Data$yts_8a_full,
                                        Network_Data$yts_9a_full,
                                        Network_Data$yts_10a_full,
                                        Network_Data$yts_11a_full,
                                        Network_Data$yts_13a_full));nrow(severity) # 311

## making the data frame for the frequency network
frequency         <- na.omit(data.frame(Network_Data$yts_1b_full,
                                        Network_Data$yts_2b_full,
                                        Network_Data$yts_3b_full,
                                        Network_Data$yts_4b_full,
                                        Network_Data$yts_5b_full,
                                        Network_Data$yts_7b_full,
                                        Network_Data$yts_8b_full,
                                        Network_Data$yts_9b_full_II,
                                        Network_Data$yts_10b_full_II,
                                        Network_Data$yts_11b_full_II,
                                        Network_Data$yts_13b_full));nrow(frequency) # 312

## define the layout for thenetwork  model (so that the models are replicable)
Layout <- structure(c(0.0783202661557321, -0.507528453222523, 0.0907182736002627, 
                      -1, -0.448827907351369, 0.540670007136345, -0.564085490971391, 
                      -0.763954236336543, 0.422874400479789, -0.160677666011429, 1, 
                      1, 0.121129059006913, -1, -0.237633154992276, -0.570585801969079, 
                      0.55656607801666, -0.993090905825155, 0.62359416364442, -0.510536133029634, 
                      -0.242636163146609, -0.124438351227588), .Dim = c(11L, 2L))

## presence vs absence network
RFtype<- c(rep("c", 11)); RFtype
RFlev<- c(rep(2,11)); RFlev
set.seed(1)
boot_graph_pva <- estimateNetwork(pres_vs_abs,
                                  default = "mgm",
                                  type = RFtype,
                                  lev = RFlev,
                                  rule = "OR",
                                  order = 2,
                                  binarySign = TRUE,
                                  criterion = "CV",
                                  nFolds = 10)
plot(boot_graph_pva, layout = Layout)

## severity network
RFtype<- c(rep("g", 11)); RFtype
RFlev<- c(rep(1,11)); RFlev
set.seed(1)
boot_graph_sev <- estimateNetwork(severity,
                                  default = "mgm",
                                  type = RFtype,
                                  lev = RFlev,
                                  rule = "OR",
                                  order = 2,
                                  binarySign = TRUE,
                                  criterion = "CV",
                                  nFolds = 10)
plot(boot_graph_sev, layout = Layout)

## frequency network
RFtype<- c(rep("g", 7), rep("c", 3), rep("g", 1)); RFtype
RFlev<- c(rep(1,7),rep(2,3),rep(1,1)); RFlev
set.seed(1)
boot_graph_freq <- estimateNetwork(frequency,
                                   default = "mgm",
                                   type = RFtype,
                                   lev = RFlev,
                                   rule = "OR",
                                   order = 2,
                                   binarySign = TRUE,
                                   criterion = "CV",
                                   nFolds = 10)
plot(boot_graph_freq, layout = Layout)

## accuracy: presence vs absence network
set.seed(1)
pres_vs_abs_accuracy <- bootnet(boot_graph_pva, nBoots = 1000, nCores = 8, type = "nonparametric", statistics = c("edge"))
plot(pres_vs_abs_accuracy, labels = TRUE, order = "id")

## accuracy: severity network
set.seed(1)
severity_accuracy <- bootnet(boot_graph_sev, nBoots = 1000, nCores = 8, type = "nonparametric", statistics = c("edge"))
plot(severity_accuracy, labels = TRUE, order = "id")

## accuracy: frequency network
set.seed(1)
frequency_accuracy <- bootnet(boot_graph_freq, nBoots = 1000, nCores = 8, type = "nonparametric", statistics = c("edge"))
plot(frequency_accuracy, labels = TRUE, order = "id")
--------------------------------------------------------------------------------------------------------------
  
   

  
--------------------------------------------------------------------------------------------------------------
#### 6. Validation
  
#####a) perceived stress
######simple regression models
lm1 <- lm(data$pss_sum ~ data$factorscore1)
summary (lm1)
lm2 <- lm (data$pss_sum ~ data$factorscore2)
summary (lm2)
lm3 <- lm (data$pss_sum ~ data$factorscore3)
summary (lm3)
####stepwise regression models###
auswahl <- na.omit(data.frame(data$pss_sum, data$factorscore1, data$factorscore2, data$factorscore3))
step.input <- lm(data.pss_sum ~ 1, data = auswahl)
step.output <- step(step.input, scope = ~ data.factorscore1 +
                      data.factorscore2 + data.factorscore3, data = auswahl)
summary (step.output)

######b) GHQ
####simple regression models
lm1 <- lm (data$w_sum ~ data$factorscore1)
summary (lm1)
lm2 <- lm (data$w_sum ~ data$factorscore2)
summary (lm2)
lm3 <- lm (data$w_sum ~ data$factorscore3)
summary (lm3)
#####stepwise regression
auswahl <- na.omit(data.frame(data$w_sum, data$factorscore1, data$factorscore2, data$factorscore3))
step.input <- lm(data.w_sum ~ 1, data = auswahl)
step.output <- step(step.input, scope = ~ data.factorscore1 +
                      data.factorscore2 + data.factorscore3, data = auswahl)
summary (step.output)


##### c) CCMS_Scale
###simple regrssion models
lm1 <- lm (data$ccms_sum ~ data$factorscore1)
summary (lm1)
lm2 <- lm (data$ccms_sum ~ data$factorscore2)
summary (lm2)
lm3 <- lm (data$ccms_sum ~ data$factorscore3)
summary (lm3)
####stepwise regresssion
auswahl <- na.omit(data.frame(data$ccms_sum, data$factorscore1, data$factorscore2, data$factorscore3))
step.input <- lm(data.ccms_sum ~ 1, data = auswahl)
step.output <- step(step.input, scope = ~ data.factorscore1 +
                      data.factorscore2 + data.factorscore3, data = auswahl)
summary (step.output)

####d )Medication
###simpe logisitc regression
auswahl <- na.omit(data.frame(data$medication, data$factorscore1,data$factorscore2, data$factorscore3))
glm1 <- glm(data.medication ~ data.factorscore1, family = binomial (logit), data = auswahl)
summary (glm1)
exp(coef(glm1))
exp(confint(glm1))
glm2 <- glm(data.medication ~ data.factorscore2, family = binomial (logit), data = auswahl)
summary (glm2)
exp(coef(glm2))
exp(confint(glm2))
glm3 <- glm(data.medication ~ data.factorscore3, family = binomial (logit), data = auswahl)
summary (glm3)
exp(coef(glm3))
exp(confint(glm3))
##stepwise logisitic regression 
step.input <- glm(data.medication ~ 1, family = binomial (logit), data = auswahl)
step.output <- step(step.input, scope = ~ data.factorscore1 +
                      data.factorscore2 + data.factorscore3, data = auswahl)
summary (step.output)
#####ORs and CI for ORs####
exp(coef(step.output))
exp(confint(step.output))


####e )Therapy
###simple logistic regression
auswahl <- na.omit(data.frame(data$therapy, data$factorscore1, data$factorscore2, data$factorscore3))
glm1 <- glm(data.therapy ~ data.factorscore1, family = binomial (logit), data = auswahl)
summary (glm1)
exp(coef(glm1))
exp(confint(glm1))
glm2 <- glm(data.therapy ~ data.factorscore2, family = binomial (logit), data = auswahl)
summary (glm2)
exp(coef(glm2))
exp(confint(glm2))
glm3 <- glm(data.therapy ~ data.factorscore3, family = binomial (logit), data = auswahl)
summary (glm3)
exp(coef(glm3))
exp(confint(glm3))

###stepwise regression
step.input <- glm(data.therapy ~ 1, data = auswahl, family = binomial (logit))
step.output <- step(step.input, scope = ~ data.factorscore1 +
                      data.factorscore2 + data.factorscore3, data = auswahl)
summary (step.output)
#####ORs and CI for ORs####
exp(coef(step.output))
exp(confint(step.output))


################influence of age and gender on factorscores
gen <- lm (data$factorscore1 ~ data$gender)
summary (gen)
gen1 <- lm (data$factorscore2 ~ data$gender)
summary (gen1)
data$age <- as.factor (data$age)
age <- lm (data$factorscore1 ~ data$age)
summary (age)
age1<- lm (data$factorscore2 ~ data$age)
summary (age1)
--------------------------------------------------------------------------------------------------------------
  


  
--------------------------------------------------------------------------------------------------------------
#### 7. supplemenatary analyses
##Supplementary material 1.1 - factor analysis up to the age of 18
  
new <- data.frame (data$id, data$yts_1, data$yts_2, data$yts_3, data$yts_4, data$yts_5, data$yts_6,
                     data$yts_7, data$yts_8, data$yts_9, data$yts_10, data$yts_11, data$yts_12,data$yts_13,
                     agebin$data.yts_1c, agebin$data.yts_2c, agebin$data.yts_3c, agebin$data.yts_4c, 
                     agebin$data.yts_5c, agebin$data.yts_6c, agebin$data.yts_7c, agebin$data.yts_8c,
                     agebin$data.yts_9c, agebin$data.yts_10c, agebin$data.yts_11c, agebin$data.yts_12c,
                     agebin$data.yts_13c) 

#######additional factor analysis#########
new$data.yts_1 <- ifelse (new$agebin.data.yts_1c == 1 | new$agebin.data.yts_1c == 2
                          | new$data.yts_1 == 0, new$data.yts_1, 0) 

new$data.yts_2 <- ifelse (new$agebin.data.yts_2c == 1 | new$agebin.data.yts_2c == 2
                          | new$data.yts_2 == 0, new$data.yts_2, 0) 

new$data.yts_3 <- ifelse (new$agebin.data.yts_3c == 1 | new$agebin.data.yts_3c == 2
                          | new$data.yts_3 == 0 , new$data.yts_3, 0) 

new$data.yts_4 <- ifelse (new$agebin.data.yts_4c == 1 | new$agebin.data.yts_4c == 2
                          | new$data.yts_4 == 0, new$data.yts_4, 0) 

new$data.yts_5 <- ifelse (new$agebin.data.yts_5c == 1 | new$agebin.data.yts_5c == 2
                          | new$data.yts_5 == 0, new$data.yts_5, 0) 

new$data.yts_6 <- ifelse (new$agebin.data.yts_6c == 1 | new$agebin.data.yts_6c == 2
                          | new$data.yts_6 == 0, new$data.yts_6, 0) 

new$data.yts_7 <- ifelse (new$agebin.data.yts_7c == 1 | new$agebin.data.yts_7c == 2
                          | new$data.yts_7 == 0, new$data.yts_7, 0) 

new$data.yts_8 <- ifelse (new$agebin.data.yts_8c == 1 | new$agebin.data.yts_8c == 2
                          | new$data.yts_8 == 0, new$data.yts_8, 0) 

new$data.yts_9 <- ifelse (new$agebin.data.yts_9c == 1 | new$agebin.data.yts_9c == 2
                          | new$data.yts_9 == 0, new$data.yts_9, 0) 

new$data.yts_10 <- ifelse (new$agebin.data.yts_10c == 1 | new$agebin.data.yts_10c == 2
                           | new$data.yts_10 == 0, new$data.yts_10, 0) 

new$data.yts_11 <- ifelse (new$agebin.data.yts_11c == 1 | new$agebin.data.yts_11c == 2
                           | new$data.yts_11 == 0, new$data.yts_11, 0)

new$data.yts_12 <- ifelse (new$agebin.data.yts_12c == 1 | new$agebin.data.yts_12c == 2
                           | new$data.yts_12 == 0, new$data.yts_12, 0)

new$data.yts_13 <- ifelse (new$agebin.data.yts_13c == 1 | new$agebin.data.yts_13c == 2
                           | new$data.yts_13 == 0, new$data.yts_13, 0) 


new1 <- data.frame(new$data.yts_1,new$data.yts_2, new$data.yts_3, new$data.yts_4, new$data.yts_5, new$data.yts_6, 
                   new$data.yts_7, new$data.yts_8, new$data.yts_9, new$data.yts_10, new$data.yts_11, 
                   new$data.yts_13)

new1 <- na.omit (new1)

matrixnew <- polychoric2(new1)$rho   #####correlation matrix for Yes/no Items
matrixnew

nfact <-12 ###at first as many factors as Items
factors.pc <- principal(matrixnew, nfact)

## Kaiser-Criterion & Scree Plot for yes/no Items
plot(c(1:nfact), factors.pc$values, type = "b", xlab = "Factor", ylab = "Eigenvalue")
fa.parallel(matrixnew,n.obs=nrow(matrixnew),fa="fa",main="Parallel Analysis Scree Plots")

####factor loadings####
######depending on the suggested solution which was 1 factor
faYN <- fa (matrixnew, 1, fm = "wls", rotate = "promax") ###yes/no Items
print (faYN, digits = 2)

m.cong <- 'Trauma =~ new.data.yts_2 + new.data.yts_3 + new.data.yts_4 + new.data.yts_5 + new.data.yts_6 + 
                     new.data.yts_7 + new.data.yts_8 + new.data.yts_9 + new.data.yts_10 + new.data.yts_11 + 
                     new.data.yts_13

new.data.yts_2 ~~ new.data.yts_9
'

#running a cfa
r.cong <- cfa(m.cong, data=new1, meanstructure=T, std.lv=T, mimic = "Mplus", 
              ordered=c("new.data.yts_2","new.data.yts_3", "new.data.yts_4", "new.data.yts_5", 
                        "new.data.yts_6", "new.data.yts_7","new.data.yts_8", "new.data.yts_9", "new.data.yts_10",
                        "new.data.yts_11","new.data.yts_13"))
##results
summary(r.cong, fit.measures=TRUE, standardized = T)

fitMeasures(r.cong,c("chisq","df","pvalue","cfi","tli","rmsea","wrmr"))
fitMeasures(r.cong,c("chisq","df","wrmr","pvalue",'cfi.scaled','tli.scaled','rmsea.scaled','rmsea.ci.lower.scaled','rmsea.ci.upper.scaled'))

#######Omega H 
parameterEstimates(r.cong)
loading2 <- parameterEstimates(r.cong)$est[1:11] 
Errors2 <- parameterEstimates(r.cong)$est[12:22]
(omega_H = (sum(loading2)^2)/(sum(loading2)^2+sum(Errors2)))



############## Supplementary material 3. Descriptive statistics severity items and frequency items for people if event endorsed 'yes'

data$yts_1an <- ifelse (data$yts_1a >= 1, data$yts_1a, NA) 
describe (data$yts_1an)
quantile (data$yts_1an, na.rm = T)

data$yts_2an <- ifelse (data$yts_2a >= 1, data$yts_2a, NA) 
table (data$yts_2an)
describe (data$yts_2an)
quantile (data$yts_2an, na.rm = T)

data$yts_3an <- ifelse (data$yts_3a >= 1, data$yts_3a, NA) 
table (data$yts_3an)
describe (data$yts_3an)
quantile (data$yts_3an, na.rm = T)

data$yts_4an <- ifelse (data$yts_4a >= 1, data$yts_4a, NA) 
table (data$yts_4an)
describe (data$yts_4an)
quantile (data$yts_4an, na.rm = T)

data$yts_5an <- ifelse (data$yts_5a >= 1, data$yts_5a, NA) 
table (data$yts_5an)
describe (data$yts_5an)
quantile (data$yts_5an, na.rm = T)

data$yts_6an <- ifelse (data$yts_6a >= 1, data$yts_6a, NA) 
table (data$yts_6an)
describe (data$yts_6an)
quantile (data$yts_6an, na.rm = T)

data$yts_7an <- ifelse (data$yts_7a >= 1, data$yts_7a, NA) 
table (data$yts_7an)
describe (data$yts_7an)
quantile (data$yts_7an, na.rm = T)

data$yts_8an <- ifelse (data$yts_8a >= 1, data$yts_8a, NA) 
table (data$yts_8an)
describe (data$yts_8an)
quantile (data$yts_8an, na.rm = T)

data$yts_9an <- ifelse (data$yts_9a >= 1, data$yts_9a, NA) 
table (data$yts_9an)
describe (data$yts_9an)
quantile (data$yts_9an, na.rm = T)

data$yts_10an <- ifelse (data$yts_10a >= 1, data$yts_10a, NA) 
table (data$yts_10an)
describe (data$yts_10an)
quantile (data$yts_10an, na.rm = T)

data$yts_11an <- ifelse (data$yts_11a >= 1, data$yts_11a, NA) 
table (data$yts_11an)
describe (data$yts_11an)
quantile (data$yts_11an, na.rm = T)

data$yts_12an <- ifelse (data$yts_12a >= 1, data$yts_12a, NA) 
table (data$yts_12an)
describe (data$yts_12an)
quantile (data$yts_12an, na.rm = T)

data$yts_13an <- ifelse (data$yts_13a >= 1, data$yts_13a, NA) 
table (data$yts_13an)
describe (data$yts_13an)
quantile (data$yts_13an, na.rm = T)

data$yts_1bn <- ifelse (data$yts_1b >= 1, data$yts_1b, NA) 
table (data$yts_1bn)
describe (data$yts_1bn)
quantile (data$yts_1bn, na.rm = T)

data$yts_2bn <- ifelse (data$yts_2b >= 1, data$yts_2b, NA) 
table (data$yts_2bn)
describe (data$yts_2bn)
quantile (data$yts_2bn, na.rm = T)

data$yts_3bn <- ifelse (data$yts_3b >= 1, data$yts_3b, NA) 
table (data$yts_3bn)
describe (data$yts_3bn)
quantile (data$yts_3bn, na.rm = T)

data$yts_4bn <- ifelse (data$yts_4b >= 1, data$yts_4b, NA) 
table (data$yts_4bn)
describe (data$yts_4bn)
quantile (data$yts_4bn, na.rm = T)

data$yts_5bn <- ifelse (data$yts_5b >= 1, data$yts_5b, NA) 
table (data$yts_5bn)
describe (data$yts_5bn)
quantile (data$yts_5bn, na.rm = T)

data$yts_6bn <- ifelse (data$yts_6b >= 1, data$yts_6b, NA) 
table (data$yts_6bn)
describe (data$yts_6bn)
quantile (data$yts_6bn, na.rm = T)

data$yts_7bn <- ifelse (data$yts_7b >= 1, data$yts_7b, NA) 
table (data$yts_7bn)
describe (data$yts_7bn)
quantile (data$yts_7bn, na.rm = T)

data$yts_8bn <- ifelse (data$yts_8b >= 1, data$yts_8b, NA) 
table (data$yts_8bn)
describe (data$yts_8bn)
quantile (data$yts_8bn, na.rm = T)


data$yts_11bn <- ifelse (data$yts_11b >= 1, data$yts_11b, NA) 
table (data$yts_11bn)
describe (data$yts_11bn)
quantile (data$yts_11bn, na.rm = T)

data$yts_12bn <- ifelse (data$yts_12b >= 1, data$yts_12b, NA) 
table (data$yts_12bn)
describe (data$yts_12bn)
quantile (data$yts_12bn, na.rm = T)

data$yts_13bn <- ifelse (data$yts_13b >= 1, data$yts_13b, NA) 
table (data$yts_13bn)
describe (data$yts_13bn)
quantile (data$yts_13bn, na.rm = T)

#####Supplementary material 4 - comp with sumscores
####Note: run it once with collapsiing items beforehand and once without
data$sumscore1w <- data$yts_1 + data$yts_2 + data$yts_3  + data$yts_4  + data$yts_5  +  data$yts_7 +
                   data$yts_8 + data$yts_9 + data$yts_10 + data$yts_11  +  data$yts_13

data$sumscore2w <- data$yts_1a_full  + data$yts_2a_full  + data$yts_3a_full  + data$yts_4a_full +
                   data$yts_5a_full  + data$yts_7a_full  + data$yts_8a_full  + data$yts_9a_full +
                   data$yts_10a_full + data$yts_11a_full +  data$yts_13a_full

data$sumscore3w <- data$yts_1b_full  + data$yts_2b_full  + data$yts_3b_full  + data$yts_4b_full +
                   data$yts_5b_full  + data$yts_7b_full  + data$yts_8b_full  + data$yts_9b_full +
                   data$yts_10b_full + data$yts_11b_full +  data$yts_13b_full

cor.test(data$sumscore1w, data$factorscore1)
cor.test(data$sumscore2w, data$factorscore2)
cor.test(data$sumscore2w, data$factorscore1)
cor.test(data$sumscore1w, data$factorscore2)
cor.test(data$factorscore2, data$factorscore1)
cor.test(data$sumscore1w, data$sumscore2w)
cor.test(data$factorscore1, data$factorscore3)
cor.test(data$factorscore1, data$sumscore3w)
cor.test(data$factorscore2, data$factorscore3)
cor.test(data$factorscore2, data$sumscore3w)
cor.test(data$sumscore1w, data$factorscore3)
cor.test(data$sumscore1w, data$sumscore3w)
cor.test(data$sumscore2w, data$factorscore3)
cor.test(data$sumscore2w, data$sumscore3w)
cor.test(data$sumscore3w, data$sumscore3w)

lm1 <- lm(data$pss_sum ~ data$sumscore1w)
summary (lm1)
lm2 <- lm(data$w_sum ~ data$sumscore1w)
summary (lm2)
lm3 <- lm(data$ccms_sum ~ data$sumscore1w)
summary (lm3)

lm1.1 <- lm(data$pss_sum ~ data$sumscore2w)
summary (lm1.1)
lm2.1 <- lm(data$w_sum ~ data$sumscore2w)
summary (lm2.1)
lm3.1 <- lm(data$ccms_sum ~ data$sumscore2w)
summary (lm3.1)


lm1.2 <- lm(data$pss_sum ~ data$sumscore3w)
summary (lm1.2)
lm2.2 <- lm(data$w_sum ~ data$sumscore3w)
summary (lm2.2)
lm3.2 <- lm(data$ccms_sum ~ data$sumscore3w)
summary (lm3.2)
data$w_sum_o2

glm1 <- glm(data$therapy ~ data$sumscore1w, family = binomial (logit))
summary (glm1)
exp(coef(glm1))
exp(confint(glm1))
glm1.1 <- glm(data$therapy ~ data$sumscore2w, family = binomial (logit))
summary (glm1.1)
exp(coef(glm1.1))
exp(confint(glm1.1))
glm2.1 <- glm(data$therapy ~ data$sumscore3w, family = binomial (logit))
summary (glm2.1)
exp(coef(glm2.1))
exp(confint(glm2.1))

glm2 <- glm(data$medication ~ data$sumscore1w, family = binomial (logit))
summary (glm2)
exp(coef(glm2))
exp(confint(glm2))
glm2.1 <- glm(data$medication ~ data$sumscore2w, family = binomial (logit))
summary (glm2.1)
exp(coef(glm2.1))
exp(confint(glm2.1))
glm3.1 <- glm(data$medication ~ data$sumscore3w, family = binomial (logit))
summary (glm3.1)
exp(coef(glm3.1))
exp(confint(glm3.1))



####### categorize sum scores
data$threecat <- ifelse (data$sumscore1w == 0, 0, 
                         (ifelse (data$sumscore1w == 1 | data$sumscore1w == 2, 1, 2)))

table (data$threecat)

w <- lm (data$pss_sum ~ data$threecat)
summary (w)
w <- lm (data$w_sum ~ data$threecat)
summary (w)
w <- lm (data$ccms_sum ~ data$threecat)
summary (w)
glm3.1 <- glm(data$medication ~ data$threecat, family = binomial (logit))
summary (glm3.1)
exp(coef(glm3.1))
exp(confint(glm3.1))
glm3.1 <- glm(data$therapy ~ data$threecat, family = binomial (logit))
summary (glm3.1)
exp(coef(glm3.1))
exp(confint(glm3.1))
--------------------------------------------------------------------------------------------------------------