## ==============================================================================
## JTHFT — Dimensionality and Internal Consistency
## Affected vs. Non-Affected Hand
## Full analytic pipeline: PMM imputation -> EFA (Principal Axis Factoring) ->
## Internal consistency -> Sensitivity analyses
## ==============================================================================

library(psych)

set.seed(42)

## ------------------------------------------------------------------------
## 1. DATA IMPORT
## ------------------------------------------------------------------------
# Expects "JTHFT_final_dataset_complete.xlsx" / equivalent CSV with columns:
# ID, Lateralidade, Mao_Afetada_D1_E2, JT_I1_A...JT_I7_A, JT_I1_NA...JT_I7_NA
# Raw values of 0 indicate task abandonment (not missing data collection).

df <- read.csv("JTHFT_raw_data.csv")          # raw data (zeros = abandoned trials)
lat <- read.csv("JTHFT_laterality.csv")        # Lateralidade, Mao_Afetada_D1_E2

cols_A  <- c("JT_I1_A","JT_I2_A","JT_I3_A","JT_I4_A","JT_I5_A","JT_I6_A","JT_I7_A")
cols_NA <- c("JT_I1_NA","JT_I2_NA","JT_I3_NA","JT_I4_NA","JT_I5_NA","JT_I6_NA","JT_I7_NA")

## ------------------------------------------------------------------------
## 2. IMPUTATION — Predictive Mean Matching (Little, 1988, Section 4.1, Eq. 2)
## ------------------------------------------------------------------------
# For each abandoned task: regress the task on the remaining six tasks using
# completing participants (donors); predict values for donors and the
# recipient; match the recipient to the single nearest donor (K = 1) by
# predicted value; impute the donor's observed value.

little_pmm_impute <- function(data, cols) {
  out <- data
  for (target in cols) {
    if (any(out[[target]] == 0)) {
      other_cols <- setdiff(cols, target)
      donors     <- out[out[[target]] != 0, ]
      recipients <- out[out[[target]] == 0, ]

      form <- as.formula(paste(target, "~", paste(other_cols, collapse = " + ")))
      fit  <- lm(form, data = donors)
      donor_pred <- predict(fit, newdata = donors)

      for (i in seq_len(nrow(recipients))) {
        recip_pred <- predict(fit, newdata = recipients[i, , drop = FALSE])
        nearest_donor <- donors[which.min(abs(donor_pred - recip_pred)), ]
        out[rownames(recipients)[i], target] <- nearest_donor[[target]]
      }
    }
  }
  out
}

full_A  <- little_pmm_impute(df, cols_A)
full_NA <- little_pmm_impute(df, cols_NA)

X_A  <- full_A[, cols_A]
X_NA <- full_NA[, cols_NA]

## ------------------------------------------------------------------------
## 3. DESCRIPTIVES (Mean / SD)
## ------------------------------------------------------------------------
describe_hand <- function(X, label) {
  cat("\n===", label, "===\n")
  print(round(rbind(mean = colMeans(X), sd = apply(X, 2, sd)), 2))
  tot <- rowSums(X)
  cat(sprintf("Total: mean = %.2f, sd = %.2f\n", mean(tot), sd(tot)))
}
describe_hand(X_A,  "Affected hand")
describe_hand(X_NA, "Non-affected hand")

## ------------------------------------------------------------------------
## 4. NORMALITY (Shapiro-Wilk) -> justifies Spearman correlation matrix
## ------------------------------------------------------------------------
cat("\n=== Shapiro-Wilk ===\n")
for (X in list(Affected = X_A, NonAffected = X_NA)) {
  for (col in colnames(X)) {
    sw <- shapiro.test(X[[col]])
    cat(sprintf("%s: W = %.3f, p = %.4f\n", col, sw$statistic, sw$p.value))
  }
}

## ------------------------------------------------------------------------
## 5. PARALLEL ANALYSIS (Horn, 1965) — Spearman-based, decides # of factors
## ------------------------------------------------------------------------
spearman_R <- function(X) cor(X, method = "spearman")

parallel_analysis_spearman <- function(n, k, n_iter = 5000) {
  sims <- matrix(0, n_iter, k)
  for (i in 1:n_iter) {
    rand <- matrix(rnorm(n * k), n, k)
    sims[i, ] <- sort(eigen(cor(rand, method = "spearman"))$values, decreasing = TRUE)
  }
  colMeans(sims)
}

pa_ref <- parallel_analysis_spearman(n = 60, k = 7)
cat("\nParallel analysis (Spearman) simulated mean eigenvalues:\n")
print(round(pa_ref, 3))

## ------------------------------------------------------------------------
## 6. SAMPLING ADEQUACY — KMO, Bartlett, item-level MSA, determinant
## ------------------------------------------------------------------------
run_adequacy <- function(X, label) {
  cat("\n===", label, "— Sampling adequacy ===\n")
  R <- spearman_R(X)
  print(KMO(R))                      # overall KMO + item-level MSA
  print(cortest.bartlett(R, n = nrow(X)))
  cat("Determinant:", format(det(R), scientific = TRUE), "\n")

  eig <- eigen(R)$values
  cat("Eigenvalues (full matrix):\n"); print(round(eig, 3))
  n_retain <- sum(eig > pa_ref)
  cat(sprintf("Factors retained (parallel analysis): %d\n", n_retain))
}
run_adequacy(X_A,  "Affected hand")
run_adequacy(X_NA, "Non-affected hand")

## ------------------------------------------------------------------------
## 7. EXPLORATORY FACTOR ANALYSIS — Principal Axis Factoring (single factor)
## ------------------------------------------------------------------------
# PAF (not ML) is used because most items violate normality (Fabrigar et al.,
# 1999). Single factor retained per parallel analysis; no rotation applied.

run_paf <- function(X, label) {
  cat("\n===", label, "— EFA (Principal Axis Factoring) ===\n")
  R <- spearman_R(X)
  fit <- fa(R, nfactors = 1, fm = "pa", rotate = "none", n.obs = nrow(X))
  print(fit$loadings, cutoff = 0.40)
  cat("Communalities:\n"); print(round(fit$communality, 3))
  cat(sprintf("Mean communality: %.3f\n", mean(fit$communality)))
  cat(sprintf("Variance explained: %.1f%%\n", fit$Vaccounted[2, 1] * 100))
  fit
}
fit_A  <- run_paf(X_A,  "Affected hand")
fit_NA <- run_paf(X_NA, "Non-affected hand")

## ------------------------------------------------------------------------
## 8. INTERNAL CONSISTENCY — Cronbach's alpha (+ CI), inter-item correlation
## ------------------------------------------------------------------------
alpha_ci <- function(a, n, k, conf = 0.95) {
  df1 <- n - 1; df2 <- (n - 1) * (k - 1)
  c(lo = 1 - (1 - a) * qf(1 - (1 - conf) / 2, df1, df2),
    hi = 1 - (1 - a) * qf((1 - conf) / 2, df1, df2))
}

run_alpha <- function(X, label) {
  cat("\n===", label, "— Cronbach's alpha ===\n")
  a <- psych::alpha(X)
  print(a$total); print(a$alpha.drop)
  ci <- alpha_ci(a$total$raw_alpha, nrow(X), ncol(X))
  cat(sprintf("95%% CI: [%.3f, %.3f]\n", ci["lo"], ci["hi"]))

  for (i in seq_along(X)) {
    a_i  <- psych::alpha(X[, -i])$total$raw_alpha
    ci_i <- alpha_ci(a_i, nrow(X), ncol(X) - 1)
    cat(sprintf("without %s: alpha = %.3f [%.3f, %.3f]\n",
                colnames(X)[i], a_i, ci_i["lo"], ci_i["hi"]))
  }

  R <- spearman_R(X)
  off_diag <- R[upper.tri(R)]
  cat(sprintf("Average inter-item correlation (Spearman): %.3f\n", mean(off_diag)))
}
run_alpha(X_A,  "Affected hand")
run_alpha(X_NA, "Non-affected hand")

## ------------------------------------------------------------------------
## 9. McDONALD'S OMEGA (from PAF loadings) + bootstrap CI
## ------------------------------------------------------------------------
omega_paf <- function(X) {
  R <- spearman_R(X)
  fit <- fa(R, nfactors = 1, fm = "pa", rotate = "none", n.obs = nrow(X))
  load1 <- as.numeric(fit$loadings)
  uniq  <- 1 - fit$communality
  (sum(load1)^2) / (sum(load1)^2 + sum(uniq))
}

bootstrap_omega <- function(X, n_boot = 10000) {
  n <- nrow(X)
  boots <- numeric(n_boot)
  b <- 1
  while (b <= n_boot) {
    idx <- sample(1:n, n, replace = TRUE)
    Xb <- X[idx, ]
    if (all(apply(Xb, 2, sd) > 0)) {
      boots[b] <- tryCatch(omega_paf(Xb), error = function(e) NA)
      if (!is.na(boots[b])) b <- b + 1
    }
  }
  quantile(boots, c(0.025, 0.975))
}

for (nm in c("A", "NA_")) {
  X <- if (nm == "A") X_A else X_NA
  label <- if (nm == "A") "Affected hand" else "Non-affected hand"
  obs <- omega_paf(X)
  ci  <- bootstrap_omega(X)
  cat(sprintf("\n%s: omega (PAF) = %.3f, 95%% CI [%.3f, %.3f]\n",
              label, obs, ci[1], ci[2]))
}

## ------------------------------------------------------------------------
## 10. PAIRED BOOTSTRAP — Between-hand comparison (alpha and omega)
## ------------------------------------------------------------------------
# Resamples participants (rows), preserving within-subject pairing, to test
# whether alpha and omega differ significantly between the affected and
# non-affected hands.

paired_bootstrap_diff <- function(X_A, X_NA, stat_fun, n_boot = 10000) {
  n <- nrow(X_A)
  diffs <- numeric(n_boot)
  b <- 1
  while (b <= n_boot) {
    idx <- sample(1:n, n, replace = TRUE)
    Xb_A  <- X_A[idx, ]
    Xb_NA <- X_NA[idx, ]
    if (all(apply(Xb_A, 2, sd) > 0) && all(apply(Xb_NA, 2, sd) > 0)) {
      diffs[b] <- tryCatch(stat_fun(Xb_A) - stat_fun(Xb_NA), error = function(e) NA)
      if (!is.na(diffs[b])) b <- b + 1
    }
  }
  ci <- quantile(diffs, c(0.025, 0.975))
  p  <- 2 * min(mean(diffs <= 0), mean(diffs >= 0))
  list(ci = ci, p = p)
}

alpha_fun <- function(X) psych::alpha(X)$total$raw_alpha

diff_alpha_obs <- alpha_fun(X_A) - alpha_fun(X_NA)
diff_omega_obs <- omega_paf(X_A) - omega_paf(X_NA)

boot_alpha <- paired_bootstrap_diff(X_A, X_NA, alpha_fun)
boot_omega <- paired_bootstrap_diff(X_A, X_NA, omega_paf)

cat("\n=== Paired bootstrap: between-hand difference ===\n")
cat(sprintf("Alpha difference: %.3f, 95%% CI [%.3f, %.3f], p = %.4f\n",
            diff_alpha_obs, boot_alpha$ci[1], boot_alpha$ci[2], boot_alpha$p))
cat(sprintf("Omega difference: %.3f, 95%% CI [%.3f, %.3f], p = %.4f\n",
            diff_omega_obs, boot_omega$ci[1], boot_omega$ci[2], boot_omega$p))

## ------------------------------------------------------------------------
## 11. SENSITIVITY ANALYSIS 1 — Imputed vs. complete-case (Task 1 / Task 7)
## ------------------------------------------------------------------------
cronbach <- function(X) {
  k <- ncol(X)
  (k / (k - 1)) * (1 - sum(apply(X, 2, var)) / var(rowSums(X)))
}

sens_task <- function(full_data, cols, target, label) {
  X_imp <- full_data[, cols]
  complete <- full_data[df[[target]] != 0, cols]   # exclude abandoners for this task
  cat(sprintf("\n--- %s sensitivity (imputed n=%d vs complete-case n=%d) ---\n",
              label, nrow(X_imp), nrow(complete)))
  cat("Alpha (item removed), imputed:      ", round(cronbach(X_imp[, setdiff(cols, target)]), 3), "\n")
  cat("Alpha (item removed), complete-case:", round(cronbach(complete[, setdiff(cols, target)]), 3), "\n")
}
sens_task(full_A,  cols_A,  "JT_I1_A",  "Affected hand — Task 1")
sens_task(full_A,  cols_A,  "JT_I7_A",  "Affected hand — Task 7")
sens_task(full_NA, cols_NA, "JT_I1_NA", "Non-affected hand — Task 1")

## ------------------------------------------------------------------------
## 12. SENSITIVITY ANALYSIS 2 — Testing order (dominance confound)
## ------------------------------------------------------------------------
merged <- merge(full_A, full_NA, by = "ID")
merged <- merge(merged, lat, by = "ID")
merged$total_A  <- rowSums(merged[, cols_A])
merged$total_NA <- rowSums(merged[, cols_NA])

g1 <- merged[merged$Mao_Afetada_D1_E2 == 1, ]  # affected = dominant, tested SECOND
g2 <- merged[merged$Mao_Afetada_D1_E2 == 2, ]  # affected = non-dominant, tested FIRST

cat("\n=== Sensitivity analysis: testing order ===\n")
cat(sprintf("Affected total  — tested 2nd (n=%d): Mdn=%.2f | tested 1st (n=%d): Mdn=%.2f\n",
            nrow(g1), median(g1$total_A), nrow(g2), median(g2$total_A)))
cat(sprintf("Non-affected total — tested 2nd (n=%d): Mdn=%.2f | tested 1st (n=%d): Mdn=%.2f\n",
            nrow(g1), median(g1$total_NA), nrow(g2), median(g2$total_NA)))
print(wilcox.test(g1$total_A, g2$total_A))
print(wilcox.test(g1$total_NA, g2$total_NA))

## ------------------------------------------------------------------------
## 13. FIGURES — Scree plots (with parallel analysis) + PAF loadings
## ------------------------------------------------------------------------
task_labels <- paste0("Task ", 1:7)
col_navy <- "#2c5282"; col_red <- "#c53030"; col_grey <- "#718096"

png("JTHFT_Figure1.png", width = 2400, height = 2000, res = 200)
par(mfrow = c(2, 2), mar = c(5, 5, 4, 2))

for (X in list(X_A, X_NA)) {
  label <- if (identical(X, X_A)) "Affected Hand" else "Non-Affected Hand"
  eig <- eigen(spearman_R(X))$values
  plot(1:7, eig, type = "b", pch = 19, col = col_navy, lwd = 2, cex = 1.3,
       ylim = c(0, max(eig) * 1.3), xlab = "Component", ylab = "Eigenvalue",
       main = paste("Scree Plot -", label))
  lines(1:7, pa_ref, type = "b", pch = 15, col = col_grey, lty = 2)
  abline(h = 1, col = col_red, lty = 3)
  legend("topright", legend = c("Observed", "Parallel analysis", "Kaiser criterion"),
         col = c(col_navy, col_grey, col_red), lty = c(1, 2, 3),
         pch = c(19, 15, NA), bty = "n", cex = 0.8)
}

for (X in list(X_A, X_NA)) {
  label <- if (identical(X, X_A)) "Affected Hand" else "Non-Affected Hand"
  fit <- fa(spearman_R(X), nfactors = 1, fm = "pa", rotate = "none", n.obs = nrow(X))
  load1 <- as.numeric(fit$loadings)
  barplot(rev(load1), horiz = TRUE, names.arg = rev(task_labels), col = col_navy,
          xlim = c(0, 1), las = 1, xlab = "Loading (PAF)",
          main = paste(label, "- Single Factor (PAF)"))
  abline(v = 0.40, col = col_red, lty = 2)
}

dev.off()
cat("\nFigure saved as JTHFT_Figure1.png\n")

## ==============================================================================
## END OF SCRIPT
## ==============================================================================