어떤 고등학교 특성이 청년기 임금 격차로 이어지는가

# ── (4) select → 라벨 제거·numeric → 임금블록 음수 결측코드 NA → 리네임 → WAVE=9 ─
block <- raw9 %>% select(all_of(c(wage_src, id_src)))
block[wage_src] <- lapply(block[wage_src], naneg)   # 임금블록: 음수→NA
block[id_src]   <- lapply(block[id_src],   num)     # ID·가중치: 라벨 제거 numeric(음수 처리 안 함)
names(block)[match(wage_src, names(block))] <- unname(map[wage_src])  # Y25SH02XXX → L08J302XXX
block$WAVE <- 9L

saveRDS(block, OUTRDS)
cat(sprintf("▶ 저장: %s (%d행 × %d열)\n", OUTRDS, nrow(block), ncol(block)))

# ── (3) 종속변수: 월임금(통합구성, 로그) ─────────────────────────────────────
pay<-naneg(raw$L08J302033); ann<-naneg(raw$L08J302034); mon<-naneg(raw$L08J302035)
wk <-naneg(raw$L08J302036); dy <-naneg(raw$L08J302037); hr <-naneg(raw$L08J302038)
hrs<-naneg(raw$L08J302031); dys<-naneg(raw$L08J302030); inc<-naneg(raw$L08J302039)
mw_form <- dplyr::case_when(
  pay==1 ~ ann/12, pay==2 ~ mon, pay==3 ~ wk*4.345,
  pay==4 ~ dy*ifelse(is.na(dys),5,dys)*4.345,
  pay==5 ~ hr*ifelse(is.na(hrs),40,hrs)*4.345, TRUE ~ NA_real_)
mw <- dplyr::coalesce(mw_form, inc, mon, ann/12, hr*ifelse(is.na(hrs),40,hrs)*4.345)

wdf <- data.frame(Y16SID=raw$Y16SID, WAVE=W, mw=mw, hrs=hrs) %>%
  filter(WAVE==TARGET_WAVE, mw>=WAGE_FLOOR, mw<=WAGE_CEIL)
if (SAMPLE=="fulltime") wdf <- wdf %>% filter(!is.na(hrs), hrs>=35)
q <- quantile(wdf$mw, WINSOR, na.rm=TRUE); wdf$mw <- pmin(pmax(wdf$mw, q[1]), q[2])
wage_df <- wdf %>% mutate(lnwage=log(mw)) %>% select(Y16SID, mw, lnwage)

# ── (5) 결측 대체 (MICE 또는 중앙값) ─────────────────────────────────────────
cat(sprintf("\n▶ 결측대체: %s\n", ifelse(USE_MICE,"MICE(CART, 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)
    Xc <- mice::complete(im, 1)
  } else Xc <- X
  Xc <- Xc %>% mutate(across(everything(), ~ ifelse(is.na(.), median(.,na.rm=TRUE), .)))
  bind_cols(df %>% select(any_of(keep)), Xc)
}
keepA <- c(admin,"is_test")
t0<-proc.time()["elapsed"]; adfA <- impute_df(adfA, keepA); adfB <- impute_df(adfB, admin)

# ── (6) A안 Post-LASSO (gaussian) ────────────────────────────────────────────
cat("\n===== A안: Post-LASSO (gaussian) =====\n")
X <- adfA %>% select(-any_of(keepA)) %>% as.matrix()
X <- X[, apply(X,2,function(c) var(c,na.rm=TRUE)>0), drop=FALSE]
y <- adfA$lnwage; w <- adfA$Y16_WEIGHT; tf <- adfA$is_test
wn <- (w[!tf])/sum(w[!tf])*sum(!tf)
cv <- cv.glmnet(X[!tf,], y[!tf], family="gaussian", alpha=1, weights=wn, 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))

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
}
drop_const <- function(vars, data, k=2)

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)
  dfit <- adfA; names(dfit) <- make.names(names(dfit), unique=TRUE)
  surv2 <- drop_const(make.names(surv,unique=TRUE), dfit)
  desA <- svydesign(id=~Y16AID, weights=~Y16_WEIGHT, data=dfit, nest=TRUE)
  f0 <- svyglm(as.formula(paste("lnwage~",paste(surv2,collapse="+"))), design=desA, family=gaussian())
  fA <- remove_high_vif_lm(f0, desA, "lnwage")
  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)
  cat(sprintf("▶ A안 정제후 %d개 | 유의(p<.05) %d개\n",
    length(attr(terms(fA),"term.labels")), sum(sigA$p<0.05)))
  print(head(sigA[sigA$p<0.05,], 20), row.names=FALSE)

# [중요] 매 회차 '훈련표본 기반 중앙값'으로 결측대체.
#  - pmm/MICE는 314개 공선 변수에서 역행렬이 특이(singular)가 되어 실패하므로 사용하지 않음.
#  - 훈련표본 중앙값만 사용 → 검증표본 정보 누수 차단 + 계산 안정 + 50회 수 초 완료.
acc <- character(0); t0<-proc.time()["elapsed"]
baseA <- bind_cols(adm, nonA)
Xr0 <- baseA %>% select(-any_of(admin)) %>% as.data.frame()
for (r in seq_len(N_REPEAT)){
  set.seed(2000+r)
  tsr <- baseA$Y16AID %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]
  yr <- baseA$lnwage; wr <- baseA$Y16_WEIGHT
  wnr <- (wr[!tsr])/sum(wr[!tsr])*sum(!tsr)
  cvr <- cv.glmnet(M[!tsr,], yr[!tsr], family="gaussian", alpha=1, weights=wnr, nfolds=10, standardize=TRUE)
  cfr <- coef(cvr, s=cvr$lambda.min)
  pr <- setdiff(rownames(cfr)[as.matrix(cfr)[,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)))
}

cat("\n===== B안: 위계 선형회귀 (gaussian) =====\n")
desB <- svydesign(id=~Y16AID, weights=~Y16_WEIGHT, 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("lnwage~",paste(vv,collapse="+"))), design=desB, family=gaussian())
  if (nm=="M4") fitM4 <- fit
  yh <- as.numeric(predict(fit, type="response"))
  R2 <- wr2(adfB$lnwage, yh, adfB$Y16_WEIGHT); aR2 <- adjr2(R2, nrow(adfB), length(vv))
  hier <- rbind(hier, data.frame(모형=nm, 투입변수=length(vv),
    R2=round(R2,4), 수정R2=round(aR2,4), ΔR2=round(R2-prevR,4))); prevR<-R2
}
print(hier, row.names=FALSE)

preds<-intersect(preds,names(ghs))
heck<-ghs%>%filter(Y16SID%in%ids_w9)%>%select(Y16SID,any_of(preds))%>%mutate(across(-Y16SID,~ifelse(num(.)<0,NA,num(.))))
heck<-left_join(heck,wage_df,by="Y16SID");heck$sel<-as.integer(!is.na(heck$lnwage))
heck[preds]<-medfill(heck[preds])
cat(sprintf("Heckman(12변수) universe N=%d, 선택(전일제 임금)=%d (%.1f%%)\n",nrow(heck),sum(heck$sel),mean(heck$sel)*100))
rhs<-paste(preds,collapse="+")
p1<-glm(as.formula(paste("sel ~",rhs)),family=binomial("probit"),data=heck)
xb<-predict(p1,type="link"); heck$imr<-dnorm(xb)/pnorm(xb)
sub<-heck[heck$sel==1,]
o0<-lm(as.formula(paste("lnwage ~",rhs)),data=sub)
o1<-lm(as.formula(paste("lnwage ~",rhs,"+ imr")),data=sub)
co<-summary(o1)$coefficients
cat("\n[Heckman 2단계·12변수] 결과식(역Mills 포함) 계수\n"); print(round(co,4))
imr_p<-co["imr","Pr(>|t|)"]; imr_b<-co["imr","Estimate"]
cat(sprintf("\n역Mills(λ) 계수=%.4f, p=%.4f → %s\n",imr_b,imr_p,
  ifelse(imr_p<0.05,"선택편향 유의","비유의=검정력부족·판단보류(식별 취약)")))

# β°,R°,β̃,R̃ 산출 후 Rmax 3수준 δ·β*
osterbase<-function(target){ ctrl<-intersect(setdiff(m4v,target),names(D))
  fs<-lm(as.formula(paste("lnwage ~",target)),data=D); b0<-coef(fs)[[target]]; R0<-summary(fs)$r.squared
  fl<-lm(as.formula(paste("lnwage ~",paste(c(target,ctrl),collapse="+"))),data=D); bt<-coef(fl)[[target]]; Rt<-summary(fl)$r.squared
  list(b0=b0,R0=R0,bt=bt,Rt=Rt) }
RES<-list()
for(t in c("L08H104018","L08H104019","L08H101009")){ o<-osterbase(t)
  cat(sprintf("\n[%s] β°=%.4f R°=%.4f | β̃=%.4f R̃=%.4f\n", t,o$b0,o$R0,o$bt,o$Rt))
  for(mult in c(1.3,2,3)){ Rmax<-min(1, mult*o$Rt)
    delta<-o$bt*(o$Rt-o$R0)/((o$b0-o$bt)*(Rmax-o$Rt))
    bstar<-o$bt-(o$b0-o$bt)*(Rmax-o$Rt)/(o$Rt-o$R0)
    cat(sprintf("   Rmax=%.2fR̃=%.4f → δ=%.3f, β*(δ=1)=%.4f  [%s]\n",
      mult,Rmax,delta,bstar, ifelse(delta<1,"δ<1 취약","δ>1 강건")))
    RES[[length(RES)+1]]<-data.frame(변수=t, 라벨=lab(t),
      β_무통제=round(o$b0,4), R2_무통제=round(o$R0,4), β_전체통제=round(o$bt,4), R2_전체통제=round(o$Rt,4),
      Rmax수준=sprintf("%.1fR̃",mult), Rmax값=round(Rmax,4),
      delta_δ=round(delta,3), bstar_δ1=round(bstar,4),
      판정=ifelse(delta<1,"δ<1(OVB 취약)","δ≥1(강건)")) }}

# ── (a) BH-FDR on A_LASSO_sig.csv ────────────────────────────────────────────
sig<-tryCatch(read.csv(paste0(OUT,"A_LASSO_sig.csv"),fileEncoding="UTF-8",check.names=FALSE),error=function(e)NULL)
bhq<-setNames(rep(NA_real_,length(report_vars)),report_vars)
sel_p<-setNames(rep(NA_real_,length(report_vars)),report_vars)
if(!is.null(sig) && nrow(sig)){
  vcol<-sig[[1]]; pcol<-as.numeric(sig[[which(tolower(names(sig))=="p")[1]]]); if(all(is.na(pcol))) pcol<-as.numeric(sig[[ncol(sig)]])
  q_all<-p.adjust(pcol, method="BH")
  cat(sprintf("(a) A_LASSO_sig %d변수 BH-FDR 보정\n", length(pcol)))
  for(v in report_vars){ i<-which(vcol==v); if(length(i)){ sel_p[v]<-pcol[i[1]]; bhq[v]<-q_all[i[1]] } }
  print(data.frame(변수=vcol, 선택후p=round(pcol,4), BHq=round(q_all,4)))

# ── (b) 다중분할 Meinshausen 집계 (무누수: 훈련중앙값) ───────────────────────
B<-50; allv<-colnames(nonA); Xr0<-as.data.frame(nonA)
padj<-matrix(1, nrow=length(allv), ncol=B, dimnames=list(allv,NULL))
selsize<-integer(B); okb<-0
for(b in 1:B){
  set.seed(2000+b)
  tf<-aid %in% sample(schools, floor(0.5*length(schools)))   # 검증(추론) 절반
  med<-vapply(Xr0, function(c) median(c[!tf],na.rm=TRUE), numeric(1))
  Xc<-Xr0; for(j in seq_along(Xc)){m<-med[j];if(is.na(m))m<-0;v<-Xc[[j]];v[is.na(v)]<-m;Xc[[j]]<-v}
  M<-as.matrix(Xc); kc<-apply(M,2,function(c) var(c[!tf])>0); M<-M[,kc,drop=FALSE]
  cv<-cv.glmnet(M[!tf,], y[!tf], family="gaussian", alpha=1, nfolds=10, standardize=TRUE)
  cf<-coef(cv,s="lambda.min"); sel<-setdiff(rownames(cf)[as.matrix(cf)[,1]!=0],"(Intercept)")
  sel<-intersect(sel, colnames(M)); selsize[b]<-length(sel); if(!length(sel)) next
  vdf<-as.data.frame(M[tf, sel, drop=FALSE]); names(vdf)<-sel
  vdf$lnwage<-y[tf]; vdf$Y16AID<-aid[tf]; vdf$ones<-1
  fit<-tryCatch(svyglm(as.formula(paste("lnwage ~",paste(sel,collapse="+"))),
    design=svydesign(id=~Y16AID,weights=~ones,data=vdf,nest=TRUE), family=gaussian()),
    error=function(e) NULL)
  if(is.null(fit)) next
  co<-summary(fit)$coefficients; okb<-okb+1
  ps<-setNames(rep(1,length(sel)),sel); for(v in sel) if(v%in%rownames(co)) ps[v]<-co[v,4]
  padj[sel,b]<-pmin(1, ps*length(sel))
}

cat(sprintf("(b) 다중분할 %d회 중 svyglm 성공 %d회, 평균 선택집합 %.1f개\n", B, okb, mean(selsize[selsize>0])))
gmin<-0.05; gg<-seq(gmin,1,by=0.025)
aggp<-sapply(allv, function(v){ pv<-padj[v,]
  Q<-sapply(gg, function(g) min(1, as.numeric(quantile(pv/g, probs=g, names=FALSE, type=1))))
  min(1, (1-log(gmin))*min(Q)) })
sel_freq<-rowMeans(padj<1)   # 선택 빈도(참고)

# ── CSV에서 9차 ΔR² 읽기 (ΔR2 = 5번째 열, 인코딩 무관 인덱스 접근) ─────────────
d9  <- read.csv(paste0(OUT,"B_hierarchy_deltaR2.csv"), fileEncoding="UTF-8")
dr9 <- as.numeric(d9[[5]])
lab <- c("M1\n합리적선택","M2\n+제도","M3\n+학생참여","M4\n+가족SES·지역")

# ── 자가검증: 기대값 부합 (M1≈.064/M2≈.019/M3≈.035/M4≈.028), 다르면 중단 ──────
exp <- c(0.064, 0.019, 0.035, 0.028)
stopifnot(length(dr9) == 4, all(abs(dr9 - exp) < 0.002))
cat("[검증] 9차 ΔR²:",
    paste(sprintf("%s=%.4f", c("M1","M2","M3","M4"), dr9), collapse="  "), "\n")

# ── 그림2 단독: 블록별 ΔR² 9차 단독 막대 (74 그림2부 스타일) ───────────────────
png(paste0(OUT,"Figure2_hierarchy_dR2_w9only.png"), 2200, 1400, res=300)
par(mar=c(4,4.5,3,1))
bp <- barplot(dr9, names.arg=lab, col="#ffc24d", ylim=c(0,max(dr9)*1.35),
   ylab="증분 ΔR² (로그 월임금)", main="위계 영역별 증분 ΔR² (9차·약26세)")
text(bp, dr9, labels=formatC(dr9, format="f", digits=3), pos=3, cex=0.75)
dev.off()
cat("저장: Figure2_hierarchy_dR2_w9only.png\n")