library(lavaan)

set.seed(20260907)
N <- 600

# Dane syntetyczne: dwie cechy stałe i odchylenia wewnątrzosobowe.
ri <- MASS::mvrnorm(N, mu = c(0, 0),
                    Sigma = matrix(c(1, -0.25, -0.25, 1), 2, 2))
w_stres <- w_sen <- matrix(NA_real_, N, 3)
w_stres[, 1] <- rnorm(N, 0, 0.75)
w_sen[, 1] <- rnorm(N, 0, 0.75)

for (t in 2:3) {
  w_stres[, t] <- 0.38 * w_stres[, t - 1] +
                   0.00 * w_sen[, t - 1] + rnorm(N, 0, 0.55)
  w_sen[, t] <- 0.44 * w_sen[, t - 1] -
                 0.21 * w_stres[, t - 1] + rnorm(N, 0, 0.55)
}

dane <- data.frame(
  stres1 = ri[, 1] + w_stres[, 1], sen1 = ri[, 2] + w_sen[, 1],
  stres2 = ri[, 1] + w_stres[, 2], sen2 = ri[, 2] + w_sen[, 2],
  stres3 = ri[, 1] + w_stres[, 3], sen3 = ri[, 2] + w_sen[, 3]
)

ri_clpm <- '
  # Random intercepts: stable between-person differences.
  RI_stres =~ 1*stres1 + 1*stres2 + 1*stres3
  RI_sen   =~ 1*sen1   + 1*sen2   + 1*sen3

  # Latent within-person deviations.
  ws1 =~ 1*stres1; ws2 =~ 1*stres2; ws3 =~ 1*stres3
  wn1 =~ 1*sen1;   wn2 =~ 1*sen2;   wn3 =~ 1*sen3
  stres1 ~~ 0*stres1; stres2 ~~ 0*stres2; stres3 ~~ 0*stres3
  sen1 ~~ 0*sen1; sen2 ~~ 0*sen2; sen3 ~~ 0*sen3

  # Equal dynamic paths across intervals: an assumption to test.
  ws2 ~ a*ws1 + b*wn1
  ws3 ~ a*ws2 + b*wn2
  wn2 ~ c*wn1 + d*ws1
  wn3 ~ c*wn2 + d*ws2

  RI_stres ~~ RI_sen
  RI_stres ~~ 0*ws1 + 0*wn1
  RI_sen   ~~ 0*ws1 + 0*wn1
  ws1 ~~ wn1; ws2 ~~ wn2; ws3 ~~ wn3
'

fit_ri <- lavaan(ri_clpm, data = dane, estimator = 'MLR',
                 meanstructure = TRUE, int.ov.free = TRUE,
                 auto.var = TRUE, auto.fix.first = FALSE,
                 auto.cov.lv.x = FALSE)

lavInspect(fit_ri, 'converged')
fitMeasures(fit_ri, c('chisq.scaled', 'df.scaled', 'pvalue.scaled',
                       'cfi.scaled', 'tli.scaled', 'rmsea.scaled',
                       'rmsea.ci.lower.scaled', 'rmsea.ci.upper.scaled', 'srmr'))
subset(parameterEstimates(fit_ri, ci = TRUE),
       op == '~' & lhs %in% c('ws2', 'wn2'))

# Classic CLPM fitted to the same data for a sensitivity comparison.
clpm <- '
  stres2 ~ a*stres1 + b*sen1
  stres3 ~ a*stres2 + b*sen2
  sen2 ~ c*sen1 + d*stres1
  sen3 ~ c*sen2 + d*stres2
  stres1 ~~ sen1
  stres2 ~~ sen2
  stres3 ~~ sen3
'
fit_clpm <- sem(clpm, data = dane, estimator = 'MLR', meanstructure = TRUE)
fitMeasures(fit_clpm, c('chisq.scaled', 'df.scaled', 'pvalue.scaled',
                         'cfi.scaled', 'tli.scaled', 'rmsea.scaled', 'srmr'))
subset(parameterEstimates(fit_clpm, ci = TRUE),
       op == '~' & lhs %in% c('stres2', 'sen2'))
sessionInfo()