library(dplyr)
library(haven)
library(psych)
library(PerformanceAnalytics)
library(corrplot)
library(Hmisc)
library(lavaan)
library(semTools)

#In the first step, a CFA analysis was conducted for the Polish sample model, 
#and the data were prepared for analysis in RStudio.


CRS_dec <- read_sav("CRS_dec.sav") %>%
  filter(COUNTRY == 3) %>%
  select(starts_with("CRS_"))

#Originally, 4 models were selected for testing:
# 1) Short version of CRS in standard form (CRS_Brief)
# 2) Short version of CRS with the inclusion of the Bi-factor, which is RW (CRS_Brief_RW)
# 3) Full version of the CRS method in standard form (CRS_full)
# 4) Full version of the CRS method with the inclusion of the Bi-factor, which is RW (CRS_full_RW)

CRS_Brief <- 'Brief_CRS =~ CRS_1 + CRS_2 + CRS_4 + CRS_5 + CRS_6 +
                            CRS_9 + CRS_16 + CRS_20 + CRS_22 + CRS_24 + 
                            CRS_25 + CRS_27 + CRS_33 + CRS_34'

CRS_Brief_RW <- 'Brief_CRS =~ CRS_1 + CRS_2 + CRS_4 + CRS_5 + CRS_6 +
                            CRS_9 + CRS_16 + CRS_20 + CRS_22 + CRS_24 + 
                            CRS_25 + CRS_27 + CRS_33 + CRS_34
                RW =~ CRS_5 + CRS_9 + CRS_16 + CRS_20 + CRS_22 + CRS_33 + CRS_34
RW ~~ 0*Brief_CRS'

CRS_full <- 'agree =~ CRS_6 + CRS_9r + CRS_11r + CRS_15r
          close =~ CRS_2 + CRS_17 + CRS_24 + CRS_28r + CRS_30
          exp =~ CRS_31 + CRS_32 + CRS_33 + CRS_34 + CRS_35
          supp =~ CRS_3 + CRS_10 + CRS_19 + CRS_25 + CRS_26 + CRS_27
          und =~ CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
          par =~ CRS_1 + CRS_4 + CRS_7r + CRS_14 + CRS_18 + CRS_23 + CRS_29r
          div =~ CRS_5r + CRS_20r'

CRS_full_RW <- 'agree =~ CRS_6 + CRS_9r + CRS_11r + CRS_15r
          close =~ CRS_2 + CRS_17 + CRS_24 + CRS_28r + CRS_30
          exp =~ CRS_31 + CRS_32 + CRS_33 + CRS_34 + CRS_35
          supp =~ CRS_3 + CRS_10 + CRS_19 + CRS_25 + CRS_26 + CRS_27
          und =~ CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
          par =~ CRS_1 + CRS_4 + CRS_7r + CRS_14 + CRS_18 + CRS_23 + CRS_29r
          div =~ CRS_5r + CRS_20r 
          RW =~ CRS_7r + CRS_9r +
          CRS_11r + CRS_15r + CRS_28r + CRS_29r + CRS_31 + CRS_32 + 
          CRS_33 + CRS_34 + CRS_35 + CRS_20r + CRS_5r +
          CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
RW ~~ 0*agree
RW ~~ 0*close
  RW ~~ 0*exp
  RW ~~ 0*supp
  RW ~~ 0*und
  RW ~~ 0*par
  RW ~~ 0*div
'
#Here, a theoretical validity/goodness-of-fit analysis of the model to the data was conducted for the four proposed models.
res.CRS_Brief <- lavaan::cfa(CRS_Brief, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_Brief, fit = TRUE, standardized = TRUE)

res.CRS_Brief_RW <- lavaan::cfa(CRS_Brief_RW, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_Brief_RW, fit = TRUE, standardized = TRUE)

res.CRS_full <- lavaan::cfa(CRS_full, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_full, fit = TRUE, standardized = TRUE)

res.CRS_full_RW <- lavaan::cfa(CRS_full_RW, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_full_RW, fit = TRUE, standardized = TRUE)

# The goodness-of-fit indices analysis indicated that in the case of the short version of the 
# method, satisfactory parameters (CFI, TLI, RMSEA, SRMR) were not achieved to accept this 
# solution. The analysis of the full method in its original version (without the Bi-factor) also
# does not suggest an optimal level of fit.
#
# An attempt to analyze using the method factor (Bi-factor; RW) allowed us to achieve the
# desired goodness-of-fit properties. Due to the report of obtaining a negative covariance 
# matrix for latent variables, we decided to conduct further model inspection to verify the 
# above issue.

lavInspect(res.CRS_full_RW, "cov.lv")
lavInspect(res.CRS_full_RW, "cor.lv")
phi <- lavInspect(res.CRS_full_RW, "cov.lv")
eigen(phi)

lavInspect(res.CRS_full, "cov.lv")
lavInspect(res.CRS_full, "cor.lv")
phi <- lavInspect(res.CRS_full, "cov.lv")
eigen(phi)

# The obtained indices in the form of a correlation matrix, covariance, as well as eigenvalues, 
# allowed us to observe that the division of labor factor exhibited an unusual level of 
# correlation with the other latent variables. We decided to conduct a detailed inspection of 
# the variance generated by the items comprising this factor. It turned out that item20 
# (CRS_20) had a high variance and standard error. This prompted us to consider three 
# models that could address this issue:
#
# 1) Removing the division of labor factor while keeping the test items CRS_5 and CRS_20 in the method factor (RW).
#    Model name: CRS_full_RW_withoutDIV
#
# 2) Removing the division of labor factor and also removing test items CRS_5 and CRS_20 in the method factor (RW).
#    Model name: CRS_full_RW_withoutDIV2
#
# 3) Stabilizing the variance of item CRS_20 by normalizing (stabilizing) its variance to the
#    lowest possible level among the other variables in the model (for this purpose, CRS_34 was selected).
#    Model name: CRS_full_RW_fixed_CRS20_variance
#

#Model 1)

CRS_full_RW_withoutDIV <- 'agree =~ CRS_6 + CRS_9r + CRS_11r + CRS_15r
          close =~ CRS_2 + CRS_17 + CRS_24 + CRS_28r + CRS_30
          exp =~ CRS_31 + CRS_32 + CRS_33 + CRS_34 + CRS_35
          supp =~ CRS_3 + CRS_10 + CRS_19 + CRS_25 + CRS_26 + CRS_27
          und =~ CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
          par =~ CRS_1 + CRS_4 + CRS_7r + CRS_14 + CRS_18 + CRS_23 + CRS_29r
          RW =~ CRS_7r + CRS_9r + CRS_5r + CRS_20r +
          CRS_11r + CRS_15r + CRS_28r + CRS_29r + CRS_31 + CRS_32 + 
          CRS_33 + CRS_34 + CRS_35 +
          CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
RW ~~ 0*agree
RW ~~ 0*close
  RW ~~ 0*exp
  RW ~~ 0*supp
  RW ~~ 0*und
  RW ~~ 0*par'

#Model 2)

CRS_full_RW_withoutDIV2 <- 'agree =~ CRS_6 + CRS_9r + CRS_11r + CRS_15r
          close =~ CRS_2 + CRS_17 + CRS_24 + CRS_28r + CRS_30
          exp =~ CRS_31 + CRS_32 + CRS_33 + CRS_34 + CRS_35
          supp =~ CRS_3 + CRS_10 + CRS_19 + CRS_25 + CRS_26 + CRS_27
          und =~ CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
          par =~ CRS_1 + CRS_4 + CRS_7r + CRS_14 + CRS_18 + CRS_23 + CRS_29r
          RW =~ CRS_7r + CRS_9r +
          CRS_11r + CRS_15r + CRS_28r + CRS_29r + CRS_31 + CRS_32 + 
          CRS_33 + CRS_34 + CRS_35 +
          CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
RW ~~ 0*agree
RW ~~ 0*close
  RW ~~ 0*exp
  RW ~~ 0*supp
  RW ~~ 0*und
  RW ~~ 0*par'

#Model 3)

CRS_full_RW_fixed_CRS20_variance <- 'agree =~ CRS_6 + CRS_9r + CRS_11r + CRS_15r
          close =~ CRS_2 + CRS_17 + CRS_24 + CRS_28r + CRS_30
          exp =~ CRS_31 + CRS_32 + CRS_33 + CRS_34 + CRS_35
          supp =~ CRS_3 + CRS_10 + CRS_19 + CRS_25 + CRS_26 + CRS_27
          und =~ CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
          par =~ CRS_1 + CRS_4 + CRS_7r + CRS_14 + CRS_18 + CRS_23 + CRS_29r
          div =~ CRS_5r + CRS_20r 
          RW =~ CRS_7r + CRS_9r +
          CRS_11r + CRS_15r + CRS_28r + CRS_29r + CRS_31 + CRS_32 + 
          CRS_33 + CRS_34 + CRS_35 + CRS_20r + CRS_5r +
          CRS_8 + CRS_12 + CRS_13 + CRS_16 + CRS_21 + CRS_22
RW ~~ 0*agree
RW ~~ 0*close
  RW ~~ 0*exp
  RW ~~ 0*supp
  RW ~~ 0*und
  RW ~~ 0*par
  RW ~~ 0*div
CRS_34~~a*CRS_34
CRS_20r~~a*CRS_20r'

###LEGEND###
# agree = Coparenting agreement
# close = Coparenting closeness
# exp = Exposure to conflict
# supp = Coparenting support
# und = 	Coparenting undermining	
# par = Endorsement partner’s parenting
# RW = Reverse worded items

#Conducting analyses for the above three models:

#Model 1)
res.CRS_full_RW_withoutDIV <- lavaan::cfa(CRS_full_RW_withoutDIV, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_full_RW_withoutDIV, fit = TRUE, standardized = TRUE)

lavInspect(res.CRS_full_RW_withoutDIV, "cov.lv")
lavInspect(res.CRS_full_RW_withoutDIV, "cor.lv")
phi <- lavInspect(res.CRS_full_RW_withoutDIV, "cov.lv")
eigen(phi)

#Model 2)

res.CRS_full_RW_withoutDIV2 <- lavaan::cfa(CRS_full_RW_withoutDIV2, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_full_RW_withoutDIV2, fit = TRUE, standardized = TRUE)

lavInspect(res.CRS_full_RW_withoutDIV2, "cov.lv")
lavInspect(res.CRS_full_RW_withoutDIV2, "cor.lv")
phi <- lavInspect(res.CRS_full_RW_withoutDIV2, "cov.lv")
eigen(phi)

#Model 3)

res.CRS_full_RW_fixed_CRS20_variance <- lavaan::cfa(CRS_full_RW_fixed_CRS20_variance, CRS_dec, estimator = "MLR", mimic="MPLUS")
summary(res.CRS_full_RW_fixed_CRS20_variance, fit = TRUE, standardized = TRUE)

resid(res.CRS_full_RW_fixed_CRS20_variance, type='normalized')

lavInspect(res.CRS_full_RW_fixed_CRS20_variance, "cov.lv")
lavInspect(res.CRS_full_RW_fixed_CRS20_variance, "cor.lv")
phi <- lavInspect(res.CRS_full_RW_fixed_CRS20_variance, "cov.lv")
eigen(phi)

# In the above cases, convergence was observed, and a positive covariance matrix for the latent variables was defined.
# Further analysis of model properties revealed that model (1) had the lowest fit indices compared to the others (2 and 3).

# Model (3), on the other hand, despite setting the variance level of item CRS_20, had a high level of standard error,
# making it difficult to sensibly interpret the division of labor factor. Consequently, model (2) was chosen,
# which exhibited high and acceptable model fit indices and lacked artifacts related to item standard errors.

# Model (2) was used for the measurement invariance analysis with respect to gender (script: CRS_2023_MI_sex).