AI에게 판단을 맡기는 시대, 준비된 청년은 따로 있다

# ── (0) 설정 ────────────────────────────────────────────────────────────────
D103    <- list.files("C:/Claude", pattern = "^103", full.names = TRUE)[1]
PATH_SAV <- list.files(D103, pattern = "raw data.*\\.SAV$", full.names = TRUE,
                       ignore.case = TRUE)[1]
PATH_CB  <- list.files(D103, pattern = "Layout.*\\.xlsx$", full.names = TRUE)[1]
OUTDIR   <- "C:/Claude/108 ai 디지털 소사이어티 논문 공모/R분석"
MISS_THR <- 0.30   # 후보 예측변수 결측률 상한
SEED     <- 2026

# 규정 확인: 1~8차 통합 원자료만 사용
stopifnot(grepl("1~8차", basename(PATH_SAV)), grepl("1~8차", basename(PATH_CB)))

dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE)
suppressPackageStartupMessages({library(haven); library(dplyr); library(readxl)})
set.seed(SEED)

# 사교육 과목별 하위문항: 사교육 경험 자체가 없는 사람(001==2)의 NA → 2(없음)
sk_sub <- intersect(paste0(sprintf("L08H104%03d", 2:30), "_hs"), names(wide))
no_sk  <- !is.na(wide$L08H104001_hs) & wide$L08H104001_hs == 2
n_fill <- sum(vapply(sk_sub, function(v) sum(no_sk & is.na(wide[[v]])), integer(1)))
for (v in sk_sub) wide[[v]][no_sk & is.na(wide[[v]])] <- 2
cat(sprintf("\n▶ 사교육 하위문항 %d개: 비경험자(001==2, %d명) 결측 %d칸 → 2(없음) 채움\n",
            length(sk_sub), sum(no_sk), n_fill))

# ── (5) 파생변수 (105 정본 122~145줄 이식) ──────────────────────────────────
nv <- paste0("L08H102", sprintf("%03d", 65:69), "_hs")   # 고3 1학기 내신 5과목
wide <- wide %>% mutate(
  naesin_avg  = nan2na(rowMeans(across(all_of(nv)), na.rm = TRUE)),
  grade_drop  = naesin_avg - L08H102153_hs,
  dosi_gyumo  = L08F101025_hs,
  is_Seoul    = ifelse(is.na(L08F101024_hs), NA, ifelse(L08F101024_hs == 1, 1, 0)),
  is_Capital  = ifelse(is.na(L08F101024_hs), NA, ifelse(L08F101024_hs %in% c(1, 4, 8), 1, 0)),
  par_edu_max = nan2na(pmax(L08F104021_hs, L08F104051_hs, na.rm = TRUE)),
  fa_income   = L08F104027_hs,
  emp_w7      = as.integer(!is.na(L08J302012_w7)),
  emp_w8      = as.integer(!is.na(L08S109001_w8)))

# ── (6) 종속변수 ────────────────────────────────────────────────────────────
wide <- wide %>% mutate(
  y_rq1  = L08S107011_w7,
  y_rq1b = nan2na(rowMeans(across(c(L08S107011_w7, L08S107012_w7)), na.rm = TRUE)),
  y_rq2  = nan2na(rowMeans(across(c(L08S108007_w8, L08S108008_w8, L08S108009_w8)), na.rm = TRUE)),
  y_use  = ifelse(is.na(L08S108001_w8), NA_integer_, as.integer(L08S108001_w8 == 1)))

# (3-2) 대학 유형 3범주: 4년제(4,5,6 · 기타7 포함=기준) / 전문대(1,2,3) / 미진학
# 미진학 = 대학 유형·계열·소재지가 모두 결측일 때만
mk_univ <- function(type, major, region) {
  none <- is.na(type) & is.na(major) & is.na(region)
  jun  <- ifelse(none, 0, ifelse(is.na(type), NA, as.integer(type %in% c(1, 2, 3))))
  list(junior = jun, none = as.integer(none))
}
U <- if (TARGET_DV == "rq1")
  mk_univ(d$univ_type_w7last, d$univ_major_w7last, d$univ_region_w7last) else
  mk_univ(d$univ_type_last,   d$univ_major_last,   d$univ_region_last)
maj <- if (TARGET_DV == "rq1") d$univ_major_w7last else d$univ_major_last
reg <- if (TARGET_DV == "rq1") d$univ_region_w7last else d$univ_region_last
cat(sprintf("\n  [대학 변수] RQ1=%s | 미진학 %d명 / 전문대 %d명 / 4년제·기타 %d명 / NA %d명\n",
            if (TARGET_DV == "rq1") "univ_*_w7last" else "univ_*_last",
            sum(U$none), sum(U$junior == 1, na.rm = TRUE),
            sum(U$junior == 0 & U$none == 0, na.rm = TRUE), sum(is.na(U$junior))))

dgt <- sprintf("L08S109%03d_w8", 1:6)
d <- d %>% mutate(
  hs_voc   = ifelse(is.na(L08H101002_hs), NA, as.integer(L08H101002_hs %in% c(3, 8))),
  hs_elite = ifelse(is.na(L08H101002_hs), NA, as.integer(L08H101002_hs %in% c(2, 4, 5, 6, 7))),
  univ_junior  = U$junior,
  univ_none    = U$none,
  # 계열 대분류 = 코드 첫자리 (1인문 2사회 3교육 4공학 5자연 6의약 7예체능) → 4·5·6 = STEM
  stem         = ifelse(U$none == 1, 0L,
                   ifelse(is.na(maj), NA, as.integer(as.integer(substr(as.character(maj), 1, 1)) %in% c(4, 5, 6)))),
  univ_capital = ifelse(U$none == 1, 0L,
                   ifelse(is.na(reg), NA, as.integer(reg %in% c(1, 4, 8)))),
  # ※ L08H104150_hs 는 '원' 단위(중위 0 / 75% 450,000 / 최대 30,000,000) → 10만원 단위로 환산
  edu_cost10   = L08H104150_hs / 100000,
  ai_use_idx   = nan2na(rowMeans(across(all_of(sprintf("L08S107%03d_w7", 1:4))), na.rm = TRUE)),
  ai_hope_idx  = nan2na(rowMeans(across(all_of(sprintf("L08S107%03d_w7", 5:8))), na.rm = TRUE)),
  digital_tasks = ifelse(emp_w8 == 0, 0, rowSums(across(all_of(dgt)) == 2, na.rm = TRUE)))

# ── (4) 분석 표본 ────────────────────────────────────────────────────────────
d$yy <- if (TARGET_DV == "rq1") d$y_rq1 else d$y_rq2
full <- d[!is.na(d$yy), ]
full$W1 <- 1.0
schools <- unique(full$Y16AID)
cat(sprintf("\n▶ 분석표본(%s): N=%d | 학교 %d개 | DV 평균 %.3f (SD %.3f, 범위 %.2f~%.2f)\n",
            TARGET_DV, nrow(full), length(schools), mean(full$yy), sd(full$yy),
            min(full$yy), max(full$yy)))

# ── (5) B안 이론 블록 ────────────────────────────────────────────────────────
vars_m1 <- c("GENDER", "par_edu_max", "fa_income", "is_Capital", "dosi_gyumo", "L08F103001_hs")
vars_m2 <- c("naesin_avg", "L08H102153_hs", "L08H102042_hs", "L08H102048_hs",
             "hs_voc", "hs_elite", "edu_cost10", "L08H104019_hs", "L08H104010_hs",
             "L08H102030_hs", "L08H102035_hs", "L08H102087_hs", "L08H102125_hs",
             "L08H102126_hs", "L08L102005_hs", "L08G101004_hs")
vars_m3 <- c("univ_junior", "univ_none", "stem", "univ_capital", "emp_w7")
if (TARGET_DV == "rq2") vars_m3 <- c(vars_m3, "emp_w8", "digital_tasks")
vars_m4 <- c(sprintf("L08S107%03d_w7", c(1:10, 13, 14)))
if (TARGET_DV == "rq2") vars_m4 <- c(vars_m4, sprintf("L08S107%03d_w7", c(11, 12, 15)))
all_theory <- unique(c(vars_m1, vars_m2, vars_m3, vars_m4))
miss_theory <- setdiff(all_theory, names(full))
cat(sprintf("\n▶ B안 이론변수: M1 %d / M2 %d / M3 %d / M4 %d (합 %d)\n",
            length(vars_m1), length(vars_m2), length(vars_m3), length(vars_m4), length(all_theory)))
if (length(miss_theory)) cat("  *** 데이터에 없어 제외:", paste(miss_theory, collapse = ", "), "\n") else
  cat("  데이터에 없는 이론변수: 없음\n")

# ── (6) A안 후보 ─────────────────────────────────────────────────────────────
DERIV_NEW <- c("hs_voc", "hs_elite", "univ_junior", "univ_none", "stem", "univ_capital", "edu_cost10")
if (TARGET_DV == "rq2") DERIV_NEW <- c(DERIV_NEW, "emp_w8", "digital_tasks")
# 재코딩으로 대체된 원변수는 후보에서 뺀다(명목코드 오해석 방지)
RECODED <- c("L08H101002_hs", "L08H104150_hs",
             "univ_type_last", "univ_major_last", "univ_region_last", "univ_adm_last",
             "univ_type_w7last", "univ_major_w7last", "univ_region_w7last")

cand0 <- DICT$var[DICT$block_M %in% c("M1", "M2", "M3", "M4", "CTRL_RISK")]
cand0 <- setdiff(union(cand0, DERIV_NEW), RECODED)
if (TARGET_DV == "rq1") {
  n_w8 <- sum(grepl("_w8$", cand0))
  cand0 <- cand0[!grepl("_w8$", cand0)]
  cat(sprintf("\n▶ [누수검사] RQ1: '_w8' 변수 %d개 제거 → 잔여 '_w8' 후보 %d개 (0이어야 함)\n",
              n_w8, sum(grepl("_w8$", cand0))))
} else {
  cand0 <- cand0[!grepl("^L08S108", cand0)]          # DV 계열
  cat(sprintf("\n▶ [누수검사] RQ2: 'L08S108*' 후보 %d개 (0이어야 함)\n", sum(grepl("^L08S108", cand0))))
}
cand0 <- intersect(cand0, names(full))
cdf     <- full[, cand0, drop = FALSE]
ok_miss <- vapply(cdf, function(x) mean(is.na(x)) < MISS_THR, logical(1))
ok_var  <- vapply(cdf, function(x) isTRUE(var(x, na.rm = TRUE) > 0), logical(1))
candA   <- cand0[ok_miss & ok_var]
cat(sprintf("▶ A안 후보: %d개 (탈락 — 결측률 %d / 무분산 %d)\n",
            length(candA), sum(!ok_miss), sum(ok_miss & !ok_var)))

admin <- c("Y16SID", "Y16AID", "W1", "yy")
if (TARGET_DV == "rq2") admin <- c(admin, "y_use")
varsB <- drop_const(intersect(all_theory, names(full)), full)
nonA  <- full[, candA, drop = FALSE]
nonB  <- full[, varsB, drop = FALSE]
adm   <- full[, admin, drop = FALSE]
if (TARGET_DV == "rq2") adm$y_use <- full$y_use
adfA  <- bind_cols(adm, nonA); adfB <- bind_cols(adm, nonB)
set.seed(SEED); test_sch <- sample(schools, floor(0.5 * length(schools)))
adfA$is_test <- adfA$Y16AID %in% test_sch
cat(sprintf("▶ 학교 50:50 분할: 훈련 %d개교(N=%d) + 검증 %d개교(N=%d) = %d개교\n",
            length(schools) - length(test_sch), sum(!adfA$is_test),
            length(test_sch), sum(adfA$is_test), length(schools)))

# ── (7) 결측 대체 ────────────────────────────────────────────────────────────
cat(sprintf("\n▶ 결측대체: %s\n", ifelse(USE_MICE, "MICE(CART, m=1, maxit=5)", "중앙값")))
impute_df <- function(df, keep) {
  X <- df %>% select(-any_of(keep))
  if (USE_MICE) {
    pm <- mice::quickpred(X, mincor = 0.1)
    im <- mice::mice(X, m = 1, method = "cart", predictorMatrix = pm, seed = SEED, printFlag = FALSE)
    X  <- mice::complete(im, 1)
  }
  X <- X %>% mutate(across(everything(), ~ ifelse(is.na(.), median(., na.rm = TRUE), .)))
  bind_cols(df %>% select(any_of(keep)), X)
}
keepA <- c(admin, "is_test")
t0 <- proc.time()["elapsed"]
adfA <- impute_df(adfA, keepA); adfB <- impute_df(adfB, admin)
cat(sprintf("  대체 완료 %s | 잔여 결측 A=%d / B=%d\n", fmt_sec(proc.time()["elapsed"] - t0),
            sum(is.na(adfA)), sum(is.na(adfB))))

remove_high_cor <- function(M, priority, thr = 0.7) {
  v <- colnames(M); keep <- c(priority[priority %in% v], setdiff(v, priority))
  cm <- cor(M[, keep, drop = FALSE], use = "pairwise.complete.obs"); rm_ <- character(0)
  for (i in seq_along(keep)) if (!keep[i] %in% rm_)
    for (j in seq_along(keep)) if (i < j && !keep[j] %in% rm_ &&
        !is.na(cm[keep[i], keep[j]]) && abs(cm[keep[i], keep[j]]) > thr) rm_ <- c(rm_, keep[j])
  setdiff(keep, rm_)
}
remove_high_vif_lm <- function(fit0, design, outcome, thr = 5, maxit = 30) {
  fit <- fit0; it <- 0
  repeat { it <- it + 1; vf <- tryCatch(car::vif(fit), error = function(e) NULL)
    if (is.null(vf)) break; if (is.matrix(vf)) vf <- vf[, 1]
    if (max(vf) <= thr || it > maxit) break
    worst <- names(which.max(vf)); tm <- setdiff(attr(terms(fit), "term.labels"), worst)
    if (!length(tm)) break
    fit <- svyglm(as.formula(paste(outcome, "~", paste(tm, collapse = "+"))),
                  design = design, family = gaussian())
  }; fit
}

# ── (8) A안 Post-LASSO ───────────────────────────────────────────────────────
cat("\n===== A안: Post-LASSO (gaussian) =====\n")
X <- as.matrix(adfA %>% select(-any_of(keepA)))
X <- X[, apply(X, 2, function(c) var(c, na.rm = TRUE) > 0), drop = FALSE]
y <- adfA$yy; tf <- adfA$is_test
cv  <- cv.glmnet(X[!tf, ], y[!tf], family = "gaussian", alpha = 1, nfolds = 10, standardize = TRUE)
cf  <- coef(cv, s = cv$lambda.min)
sel <- setdiff(rownames(cf)[as.matrix(cf)[, 1] != 0], "(Intercept)")
predte <- as.numeric(predict(cv, newx = X[tf, ], s = "lambda.min"))
cvr2 <- 1 - min(cv$cvm) / var(y[!tf]); ter2 <- wr2(y[tf], predte); termse <- rmse(y[tf], predte)
cat(sprintf("▶ LASSO 선택(λ.min) %d개 | 훈련 CV R²≈%.4f | 검증표본 R²=%.4f, RMSE=%.4f\n",
            length(sel), cvr2, ter2, termse))

sigA <- data.frame()
if (length(sel) >= 2) {
  prio  <- names(sort(abs(as.matrix(cf)[sel, 1]), decreasing = TRUE))
  surv  <- remove_high_cor(as.data.frame(X[tf, sel, drop = FALSE]), prio)
  surv2 <- drop_const(surv, adfA)
  desA  <- svydesign(id = ~Y16AID, weights = ~W1, data = adfA, nest = TRUE)
  f0    <- svyglm(as.formula(paste("yy ~", paste(surv2, collapse = "+"))),
                  design = subset(desA, is_test), family = gaussian())
  fA    <- remove_high_vif_lm(f0, subset(desA, is_test), "yy")
  co    <- summary(fA)$coefficients
  sigA  <- data.frame(변수 = rownames(co), 라벨 = sapply(rownames(co), labof),
    계수 = round(co[, 1], 4), SE = round(co[, 2], 4), t = round(co[, 3], 3),
    p = round(co[, 4], 4), row.names = NULL) %>% filter(변수 != "(Intercept)") %>% arrange(p)
  sigA$q_BH <- round(p.adjust(sigA$p, method = "BH"), 4)   # 후선택 추론(BH-FDR), 105:88 참조
  cat(sprintf("▶ r>.7 제거 후 %d개 → VIF>5 제거 후 %d개 | 검증표본 유의(p<.05) %d개\n",
              length(surv2), length(attr(terms(fA), "term.labels")), sum(sigA$p < 0.05)))
  print(head(sigA[sigA$p < 0.05, ], 25), row.names = FALSE)
}

# ── (9) B안 위계 선형회귀 ────────────────────────────────────────────────────
cat("\n===== B안: 위계 선형회귀 (gaussian, 집락 Y16AID·비가중) =====\n")
desB <- svydesign(id = ~Y16AID, weights = ~W1, data = adfB, nest = TRUE)
av <- function(v) drop_const(intersect(v, names(adfB)), adfB)
a1 <- av(vars_m1); a2 <- av(vars_m2); a3 <- av(vars_m3); a4 <- av(vars_m4)
blocks <- list(M1 = a1, M2 = c(a1, a2), M3 = c(a1, a2, a3), M4 = c(a1, a2, a3, a4))
hier <- data.frame(); prevR <- 0; fitM4 <- NULL
for (nm in names(blocks)) {
  vv <- unique(blocks[[nm]]); if (!length(vv)) next
  fit <- svyglm(as.formula(paste("yy ~", paste(vv, collapse = "+"))), design = desB, family = gaussian())
  if (nm == "M4") fitM4 <- fit
  yh <- as.numeric(predict(fit, type = "response"))
  R2 <- wr2(adfB$yy, yh)
  hier <- rbind(hier, data.frame(모형 = nm, 투입변수 = length(vv), R2 = round(R2, 4),
    수정R2 = round(adjr2(R2, nrow(adfB), length(vv)), 4), ΔR2 = round(R2 - prevR, 4),
    RMSE = round(rmse(adfB$yy, yh), 4)))
  prevR <- R2
}
print(hier, row.names = FALSE)

tagof <- function(v) ifelse(v %in% vars_m1, "M1.가정배경",
                     ifelse(v %in% vars_m2, "M2.고교시기",
                     ifelse(v %in% vars_m3, "M3.대학·노동",
                     ifelse(v %in% vars_m4, "M4.AI이용·기대", "-"))))
sigB <- data.frame(); fullB <- data.frame()
if (!is.null(fitM4)) {
  co <- summary(fitM4)$coefficients; ci <- confint(fitM4)
  fullB <- data.frame(변수 = rownames(co), 라벨 = sapply(rownames(co), labof),
    블록 = tagof(rownames(co)), 계수 = round(co[, 1], 4), SE = round(co[, 2], 4),
    t = round(co[, 3], 3), p = round(co[, 4], 4),
    CI_low = round(ci[, 1], 4), CI_high = round(ci[, 2], 4), row.names = NULL) %>%
    filter(변수 != "(Intercept)")
  sigB <- fullB %>% arrange(p)
  cat(sprintf("\n▶ B안 M4 유의(p<.05) %d개 / 전체 %d개\n", sum(sigB$p < 0.05), nrow(sigB)))
  print(head(sigB[sigB$p < 0.05, ], 30), row.names = FALSE)
}

# ── (10) Multisplit 안정성 ───────────────────────────────────────────────────
cat(sprintf("\n===== Multisplit %d회 안정성 (회차별 훈련표본 중앙값 대체) =====\n", N_REPEAT))
acc <- character(0); t0 <- proc.time()["elapsed"]
Xr0 <- as.data.frame(nonA); yr <- full$yy; aid <- full$Y16AID
for (r in seq_len(N_REPEAT)) {
  set.seed(2000 + r)
  tsr <- aid %in% sample(schools, floor(0.5 * length(schools)))
  med_tr <- vapply(Xr0, function(col) median(col[!tsr], na.rm = TRUE), numeric(1))
  Xc <- Xr0
  for (j in seq_along(Xc)) { m <- med_tr[j]; if (is.na(m)) m <- 0
    v <- Xc[[j]]; v[is.na(v)] <- m; Xc[[j]] <- v }
  M <- as.matrix(Xc); M <- M[, apply(M, 2, function(c) var(c[!tsr]) > 0), drop = FALSE]
  cvr <- cv.glmnet(M[!tsr, ], yr[!tsr], family = "gaussian", alpha = 1, nfolds = 10, standardize = TRUE)
  pr  <- setdiff(rownames(coef(cvr, s = cvr$lambda.min))[as.matrix(coef(cvr, s = cvr$lambda.min))[, 1] != 0],
                 "(Intercept)")
  if (length(pr)) acc <- c(acc, pr)
  if (r %% 10 == 0) cat(sprintf("  [%02d/%02d] 누적 %s\n", r, N_REPEAT, fmt_sec(proc.time()["elapsed"] - t0)))
}
stab <- as.data.frame(table(acc), stringsAsFactors = FALSE) %>%
  rename(변수 = acc, 선택횟수 = Freq) %>%
  mutate(생존율 = round(선택횟수 / N_REPEAT * 100, 1), 라벨 = sapply(변수, labof)) %>%
  filter(생존율 >= SURV_PCT) %>% arrange(desc(선택횟수))
cat(sprintf("▶ 생존율 %d%% 이상: %d개\n", SURV_PCT, nrow(stab)))
print(head(stab, 30), row.names = FALSE)

# ── (11) 교차검증 핵심변수 (A안 검증유의 ∩ B안 M4 유의 ∩ 생존율≥60%) ────────
sA <- if (nrow(sigA)) gsub("\\.", "", sigA$변수[sigA$p < 0.05]) else character(0)
sB <- if (nrow(sigB)) sigB$변수[sigB$p < 0.05] else character(0)
crossv <- Reduce(intersect, list(sA, sB, stab$변수))
cat(sprintf("\n===== 교차검증 핵심 결정요인 %d개 (A∩B∩생존%d%%) =====\n", length(crossv), SURV_PCT))
if (length(crossv)) for (v in crossv) cat(sprintf("  - %-22s %s\n", v, labof(v))) else
  cat("  (없음 — 기준 미조정, 있는 그대로 보고)\n")

# ── (6) 표 G2: 표준화 SES 로지스틱/선형 (105 방식) ───────────────────────────
cat("\n============ 표 G2: 표준화 SES 효과 (per +1SD, 집락 Y16AID) ============\n")
OUTC <- list(
  list(nm = "011 동의(≥4)",   y = d$ag011,   fam = "bin"),
  list(nm = "011 평균",        y = d$y_rq1,   fam = "gau"),
  list(nm = "rq2 상위3분위",   y = d$top_rq2, fam = "bin"),
  list(nm = "rq2 평균",        y = d$y_rq2,   fam = "gau"),
  list(nm = "AI앱 사용(use)",  y = d$y_use,   fam = "bin"))
AXES2 <- list(`부모 최고학력` = d$par_edu_max, `부 월소득` = d$fa_income)
MODS  <- list(단순 = "z", `통제(성별+수도권거주)` = "z + GENDER + is_Capital")

G2 <- list()
for (o in OUTC) for (an in names(AXES2)) for (mn in names(MODS)) {
  dd <- data.frame(Y16AID = d$Y16AID, W1 = 1, y = o$y, z = zstd(AXES2[[an]]),
                   GENDER = d$GENDER, is_Capital = d$is_Capital)
  dd <- dd[complete.cases(dd), ]
  est <- NA; lo <- NA; hi <- NA; pp <- NA
  if (nrow(dd) >= 50) try({
    des <- svydesign(id = ~Y16AID, weights = ~W1, data = dd, nest = TRUE)
    fit <- svyglm(as.formula(paste("y ~", MODS[[mn]])), design = des,
                  family = if (o$fam == "bin") quasibinomial() else gaussian())
    co  <- summary(fit)$coefficients["z", ]
    if (o$fam == "bin") { est <- exp(co[1]); lo <- exp(co[1] - 1.96 * co[2]); hi <- exp(co[1] + 1.96 * co[2]) }
    else                { est <- co[1];      lo <- co[1] - 1.96 * co[2];      hi <- co[1] + 1.96 * co[2] }
    pp <- co[4]
  }, silent = TRUE)
  G2[[length(G2) + 1]] <- data.frame(결과 = o$nm, SES축 = an, 모형 = mn,
    지표 = ifelse(o$fam == "bin", "OR per +1SD", "beta per +1SD"), N = nrow(dd),
    추정치 = round(est, 3), CI_low = round(lo, 3), CI_high = round(hi, 3),
    p = round(pp, 4), Bonferroni통과_005 = ifelse(is.na(pp), NA, pp < 0.05 / BONF_M),
    row.names = NULL)
}
G2 <- bind_rows(G2)

# =============================================================================
# S5. 두 기울기 차이 검정 (권고)
# =============================================================================
cat(sprintf("\n============ S5: |β₂|−|β₁| 집락 부트스트랩 %d회 ============\n", N_BOOT))
ov <- d0[!is.na(d0$y_rq1) & !is.na(d0$y_rq2) & !is.na(d0$par_edu_max) & !is.na(d0$Y16AID), ]
cat(sprintf("  7·8차 중첩 + 부모학력 비결측 N=%d | 학교 %d\n", nrow(ov), length(unique(ov$Y16AID))))
ov$z1 <- zstd(ov$y_rq1); ov$z2 <- zstd(ov$y_rq2); ov$zp <- zstd(ov$par_edu_max); ov$W1 <- 1
desO <- svydesign(id = ~Y16AID, weights = ~W1, data = ov, nest = TRUE)
c1 <- summary(svyglm(z1 ~ zp, design = desO, family = gaussian()))$coefficients["zp", ]
c2 <- summary(svyglm(z2 ~ zp, design = desO, family = gaussian()))$coefficients["zp", ]
dstat <- abs(c2[1]) - abs(c1[1])
cat(sprintf("  β₁(판단권 전망) = %.4f (SE %.4f, p=%.4g) | β₂(리터러시) = %.4f (SE %.4f, p=%.4g)\n",
            c1[1], c1[2], c1[4], c2[1], c2[2], c2[4]))

set.seed(SEED)
sch <- unique(ov$Y16AID); idxby <- split(seq_len(nrow(ov)), ov$Y16AID)
bs <- numeric(N_BOOT); b1s <- numeric(N_BOOT); b2s <- numeric(N_BOOT)
t0 <- proc.time()["elapsed"]
for (b in seq_len(N_BOOT)) {
  s  <- sample(sch, length(sch), replace = TRUE)
  ii <- unlist(idxby[as.character(s)], use.names = FALSE)
  o  <- ov[ii, ]
  p  <- zstd(o$par_edu_max); u1 <- zstd(o$y_rq1); u2 <- zstd(o$y_rq2)
  b1 <- coef(lm(u1 ~ p))[2]; b2 <- coef(lm(u2 ~ p))[2]
  b1s[b] <- b1; b2s[b] <- b2; bs[b] <- abs(b2) - abs(b1)
}
ci <- quantile(bs, c(.025, .975), na.rm = TRUE)
cat(sprintf("  |β₂|−|β₁| = %.4f | 부트스트랩 95%% CI [%.4f, %.4f] | 0 포함 %s | %s\n",
            dstat, ci[1], ci[2], ifelse(ci[1] <= 0 & ci[2] >= 0, "예", "아니오"),
            fmt_sec(proc.time()["elapsed"] - t0)))
S5 <- data.frame(
  항목 = c("β₁: 판단권 이양 전망(z) ~ 부모학력(z)", "β₂: AI 리터러시(z) ~ 부모학력(z)",
           "|β₂| − |β₁|"),
  추정치 = round(c(c1[1], c2[1], dstat), 4),
  SE_집락 = c(round(c1[2], 4), round(c2[2], 4), NA),
  p_집락 = c(round(c1[4], 4), round(c2[4], 4), NA),
  부트_CI_low = c(round(quantile(b1s, .025), 4), round(quantile(b2s, .025), 4), round(ci[1], 4)),
  부트_CI_high = c(round(quantile(b1s, .975), 4), round(quantile(b2s, .975), 4), round(ci[2], 4)),
  N = nrow(ov), 학교수 = length(sch), 부트반복 = N_BOOT, row.names = NULL)
print(S5, row.names = FALSE)
sc(S5, "S5_slope_diff.csv")

impute_multi <- function(B, m, keep_dv, seed = SEED) {
  drop_cols <- if (keep_dv) setdiff(B$admin, "yy") else B$admin
  X  <- B$adfB_raw %>% select(-any_of(drop_cols))
  pm <- mice::quickpred(X, mincor = 0.1)
  me <- mice::make.method(X); me[me != ""] <- "cart"
  if (keep_dv) stopifnot(sum(is.na(X$yy)) == 0, me[["yy"]] == "")
  im <- mice::mice(X, m = m, method = me, predictorMatrix = pm, maxit = 5,
                   seed = seed, printFlag = FALSE)
  lapply(seq_len(m), function(i) {
    cm <- med_fill(mice::complete(im, i))
    out <- bind_cols(B$adfB_raw[, B$admin, drop = FALSE], cm %>% select(-any_of("yy")))
    out$par_edu_z <- zstd(out$par_edu_max)
    out
  })
}
rubin_manual <- function(fits) {
  Q <- sapply(fits, coef); U <- sapply(fits, function(f) diag(vcov(f)))
  m <- length(fits)
  qbar <- rowMeans(Q); ubar <- rowMeans(U); b <- apply(Q, 1, var)
  tt <- ubar + (1 + 1 / m) * b
  lam <- (1 + 1 / m) * b / tt; lam[!is.finite(lam)] <- 0; lam <- pmax(lam, 1e-8)
  dfcom <- df.residual(fits[[1]])
  df_old <- (m - 1) / lam^2
  df_obs <- (dfcom + 1) / (dfcom + 3) * dfcom * (1 - lam)
  dfa <- df_old * df_obs / (df_old + df_obs)
  se <- sqrt(tt); tv <- qbar / se
  data.frame(term = names(qbar), estimate = qbar, std.error = se,
             statistic = tv, p.value = 2 * pt(-abs(tv), dfa), fmi = lam, row.names = NULL)
}
pool_fits <- function(fits) {
  out <- NULL
  try({ ps <- summary(mice::pool(mice::as.mira(fits)))
    out <- data.frame(term = as.character(ps$term), estimate = ps$estimate, std.error = ps$std.error,
                      statistic = ps$statistic, p.value = ps$p.value, fmi = NA_real_, row.names = NULL)
    attr(out, "how") <- "mice::pool" }, silent = TRUE)
  if (is.null(out)) { out <- rubin_manual(fits); attr(out, "how") <- "Rubin 수식 직접구현" }
  rownames(out) <- out$term
  out
}

cat("\n▶ 분석자료 구성\n")
B1 <- build_dv("rq1"); B2 <- build_dv("rq2")
cat(sprintf("  RQ1 N=%d | RQ2 N=%d\n", B1$N, B2$N))
cat(sprintf("  M4a(345번 이용빈도, 4점) = %s\n", paste(B1$m4a, collapse = ", ")))
cat(sprintf("  M4b RQ1(346~348번 태도, 6점) = %s\n", paste(B1$m4b, collapse = ", ")))
cat(sprintf("  M4b RQ2 = 위 + %s\n", paste(setdiff(B2$m4b, B1$m4b), collapse = ", ")))
cat("\n  ※ [중대 1 관련 설문 원문 기록] 7차 설문 348번 도입문은 '다음에 대해 얼마나 동의하십니까?'이고,\n")
cat("     1) L08S107013 = '새로운 과학기술의 발견은 미래 인류의 삶의 질을 향상시킬 것이다'\n")
cat("     2) L08S107014 = '새로운 과학기술이 가져다 주는 혜택은 소수에게 독점될 것이다'\n")
cat("     3) L08S107015 = '인간은 인공지능이나 로봇이 대체할 수 없는 능력을 갖고 있다'\n")
cat("     즉 013·014 는 AI 를 지칭하지 않는 '새로운 과학기술' 일반 문항이다(015 만 AI 를 지칭).\n")

# ── (c) 011 의 잔차 ─────────────────────────────────────────────────────────
cat("\n---- (c) 011 을 009·010·014 에 회귀한 잔차의 계층 기울기 ----\n")
VR  <- sprintf("L08S107%03d_w7", c(9, 10, 14))
okr <- complete.cases(d0[, c("y_rq1", VR)])
fr  <- lm(as.formula(paste("y_rq1 ~", paste(VR, collapse = "+"))), data = d0[okr, ])
d0$resid011 <- NA_real_; d0$resid011[okr] <- as.numeric(residuals(fr))
cat(sprintf("  잔차 산출 lm: N=%d | R²=%.4f | 계수 %s\n", sum(okr), summary(fr)$r.squared,
            paste(sprintf("%s=%.3f", VR, coef(fr)[VR]), collapse = ", ")))
rc1 <- zslope(d0, "resid011")
rc2 <- zslope(d0, "resid011", extra = c("GENDER", "is_Capital"))
S23 <- rbind(S23,
  row_slope("④ 011 잔차(009·010·014 통제)", "resid011", "잔차, 통제 없음", rc1, -.063),
  row_slope("④ 011 잔차(009·010·014 통제)", "resid011_gc", "잔차, 성별·수도권 거주 추가 통제", rc2, -.052))
print(tail(S23, 2)[, c("변수", "라벨", "N", "표준화기울기", "SE", "t", "p", "위원장값", "차이", "일치_pm002")],
      row.names = FALSE)

# ── r(emp_w8, digital_tasks) — rq2 분석표본 원자료 ──────────────────────────
cat("\n---- r(emp_w8, digital_tasks): RQ2 분석표본(원자료, 10_main.R 산식) ----\n")
F2 <- A2$full
r_raw  <- cor(F2$emp_w8, F2$digital_tasks, use = "complete.obs")
r_imp  <- cor(A2$adfA$emp_w8, A2$adfA$digital_tasks)
r_test <- cor(A2$adfA$emp_w8[A2$tf], A2$adfA$digital_tasks[A2$tf])
cat(sprintf("  원자료 전체 분석표본 N=%d → r=%.4f\n", sum(complete.cases(F2[, c("emp_w8", "digital_tasks")])), r_raw))
cat(sprintf("  대체본 전체 N=%d → r=%.4f | 대체본 검증표본 N=%d → r=%.4f (필터가 실제로 보는 값)\n",
            nrow(A2$adfA), r_imp, sum(A2$tf), r_test))
cat(sprintf("  위원장 값 .7645 와의 차이: 원자료 %+.4f / 검증표본 %+.4f\n", r_raw - .7645, r_test - .7645))
S24 <- rbind(S24, data.frame(구분 = "③ r(emp_w8, digital_tasks)", DV = "RQ2",
  항목 = c("원자료 분석표본", "대체본 전체", "대체본 검증표본(필터 적용 대상)", "위원장 재계산값"),
  라벨 = "8차 재직 여부 ↔ 업무 디지털 활용 과업 수",
  값 = round(c(r_raw, r_imp, r_test, .7645), 4), SE = NA_real_, p = NA_real_,
  N = c(sum(complete.cases(F2[, c("emp_w8", "digital_tasks")])), nrow(A2$adfA), sum(A2$tf), NA),
  비고 = "", row.names = NULL))

mk_ipw <- function(yv) {
  f  <- svyglm(as.formula(paste(yv, "~ par_edu3 + fa_inc3 + hs3 + gr4 + GENDER + is_Capital")),
               design = desL, family = quasibinomial())
  ph <- as.numeric(predict(f, type = "response"))
  sw <- mean(L[[yv]]) / ph                    # 안정화 가중치
  sw[L[[yv]] == 0] <- NA                      # 탈락자는 분석표본에 없다
  q  <- quantile(sw, TRIM, na.rm = TRUE)
  swt <- pmin(pmax(sw, q[1]), q[2])           # 1~99분위 절단
  cat(sprintf("  [%s] 잔존확률 범위 %.3f~%.3f | 가중치 평균 %.3f SD %.3f 최소 %.3f 최대 %.3f (절단 전 최대 %.3f)\n",
              yv, min(ph), max(ph), mean(swt, na.rm = TRUE), sd(swt, na.rm = TRUE),
              min(swt, na.rm = TRUE), max(swt, na.rm = TRUE), max(sw, na.rm = TRUE)))
  setNames(swt, L$Y16SID)
}
w7 <- mk_ipw("ret7"); w8 <- mk_ipw("ret8")
d$ipw7 <- unname(w7[as.character(d$Y16SID)])
d$ipw8 <- unname(w8[as.character(d$Y16SID)])