library(MASS)
library(psych)
library(lavaan)
library(semTools)

set.seed(20260907)
N <- 800

# Trzy powiązane, ale rozróżnialne konstrukty.
sigma_latent <- matrix(
  c(1.00, 0.45, -0.35,
    0.45, 1.00, -0.20,
   -0.35, -0.20, 1.00),
  nrow = 3,
  byrow = TRUE
)
eta <- MASS::mvrnorm(N, mu = c(0, 0, 0), Sigma = sigma_latent)
colnames(eta) <- c("samoregulacja", "motywacja", "przeciazenie")

utworz_pozycje <- function(cecha, prefiks, ladunki) {
  wynik <- sapply(seq_along(ladunki), function(i) {
    lambda <- ladunki[i]
    ciagla <- lambda * cecha + sqrt(1 - lambda^2) * rnorm(N)
    as.integer(cut(
      ciagla,
      breaks = c(-Inf, -1.05, -0.35, 0.35, 1.10, Inf),
      labels = FALSE
    ))
  })
  colnames(wynik) <- paste0(prefiks, seq_along(ladunki))
  wynik
}

dane <- data.frame(
  utworz_pozycje(eta[, "samoregulacja"], "s", c(.82, .78, .76, .74, .72, .70)),
  utworz_pozycje(eta[, "motywacja"], "m", c(.80, .77, .74, .71)),
  utworz_pozycje(eta[, "przeciazenie"], "p", c(.81, .78, .75, .72))
)

# Zmienne zewnętrzne do sprawdzenia hipotez trafności teoretycznej.
dane$wytrwalosc <- 0.60 * eta[, "samoregulacja"] +
                    0.15 * eta[, "motywacja"] + rnorm(N, sd = 0.72)
dane$prokrastynacja <- -0.55 * eta[, "samoregulacja"] +
                       0.10 * eta[, "przeciazenie"] + rnorm(N, sd = 0.75)
dane$wynik_nauki <- 0.32 * eta[, "samoregulacja"] +
                     0.22 * eta[, "motywacja"] + rnorm(N, sd = 0.88)

pozycje <- c(paste0("s", 1:6), paste0("m", 1:4), paste0("p", 1:4))
indeks_efa <- sample(seq_len(N), N / 2)
dane_efa <- dane[indeks_efa, ]
dane_cfa <- dane[-indeks_efa, ]

# EFA na pierwszej połowie próby. Liczby czynników nie wybieramy
# na podstawie samego kryterium wartości własnej > 1.
rho_efa <- psych::polychoric(dane_efa[pozycje])$rho
rownolegla <- psych::fa.parallel(
  rho_efa,
  n.obs = nrow(dane_efa),
  fa = "fa",
  fm = "minres",
  n.iter = 100,
  plot = FALSE
)
print(rownolegla$nfact)

efa <- psych::fa(
  r = rho_efa,
  nfactors = 3,
  n.obs = nrow(dane_efa),
  fm = "minres",
  rotate = "oblimin"
)
print(efa$loadings, cutoff = .30, sort = TRUE)

# CFA na niezależnej połowie próby.
model_pomiarowy <- '
  samoregulacja =~ s1 + s2 + s3 + s4 + s5 + s6
  motywacja =~ m1 + m2 + m3 + m4
  przeciazenie =~ p1 + p2 + p3 + p4
'

fit_cfa <- cfa(
  model_pomiarowy,
  data = dane_cfa,
  ordered = pozycje,
  estimator = "WLSMV",
  std.lv = TRUE
)
stopifnot(lavInspect(fit_cfa, "converged"))

print(round(fitMeasures(
  fit_cfa,
  c(
    "chisq.scaled", "df.scaled", "pvalue.scaled",
    "cfi.scaled", "tli.scaled", "rmsea.scaled",
    "rmsea.ci.lower.scaled", "rmsea.ci.upper.scaled", "srmr"
  )
), 3))

print(subset(
  standardizedSolution(fit_cfa),
  op == "=~",
  select = c(lhs, rhs, est.std, se, pvalue)
))

# Rzetelność kompozytowa rho_C i AVE wyprowadzone z całkowicie
# standaryzowanego rozwiązania. Jawny zapis ułatwia sprawdzenie,
# z jakich ładunków i wariancji błędu powstały wyniki.
metryki_pomiaru <- function(fit) {
  ladunki <- subset(
    standardizedSolution(fit),
    op == "=~",
    select = c(lhs, rhs, est.std)
  )

  wynik <- lapply(split(ladunki, ladunki$lhs), function(blok) {
    lambda <- blok$est.std
    theta <- 1 - lambda^2
    data.frame(
      rho_C = sum(lambda)^2 / (sum(lambda)^2 + sum(theta)),
      AVE = sum(lambda^2) / (sum(lambda^2) + sum(theta))
    )
  })
  round(do.call(rbind, wynik), 3)
}

print(metryki_pomiaru(fit_cfa))

# HTMT dotyczy rozróżnialności par konstruktów. Żadna z tych
# liczb nie zastępuje dowodów treściowych ani relacji z kryteriami.
print(semTools::htmt(
  model_pomiarowy,
  data = dane_cfa,
  ordered = pozycje,
  htmt2 = TRUE
))

# Trafność teoretyczna: wcześniej zapisane hipotezy o relacjach
# konstruktu z wytrwałością, prokrastynacją i wynikiem nauki.
model_trafnosci <- paste0(
  model_pomiarowy,
  '\nwytrwalosc ~ samoregulacja + motywacja',
  '\nprokrastynacja ~ samoregulacja + przeciazenie',
  '\nwynik_nauki ~ samoregulacja + motywacja'
)
fit_trafnosci <- sem(
  model_trafnosci,
  data = dane_cfa,
  ordered = pozycje,
  estimator = "WLSMV",
  std.lv = TRUE
)
print(subset(
  parameterEstimates(fit_trafnosci, standardized = TRUE, ci = TRUE),
  op == "~",
  select = c(lhs, rhs, est, se, pvalue, ci.lower, ci.upper, std.all)
))

write.csv(dane, "dane_walidacja_skali.csv", row.names = FALSE)
sessionInfo()