# ============================================================
# 유아교육연구의 양적접근 | 2026-1학기
# 6주차: 경로모형 실습 (R lavaan)
# 강사: 서울대학교 의학연구원 김주현
# 모형: 사회적지지 3종류 → 우울(MDEP) → 양육행동(PAR_SOCIAL_M)
# + 3종류 → 양육행동 직접효과 포함
# ============================================================
##############################################################
# ────────────────────────────────────────────────────────────
# STEP 0. 패키지 설치
# ────────────────────────────────────────────────────────────
##############################################################
library(lavaan)
## This is lavaan 0.6-21
## lavaan is FREE software! Please report any bugs.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(magrittr)
##############################################################
# ────────────────────────────────────────────────────────────
# STEP 1. 데이터 불러오기
# ────────────────────────────────────────────────────────────
##############################################################
data <- read.delim("C:/Users/Lenovo/Desktop/path_analysis/w2.dat", row.names=1)
data[data == 99999999] <- NA
df_sel <- data %>%
select(SS_EMO, # 사회적지지: 정서적 평균 (X1)
SS_INS, # 사회적지지: 도구적 평균 (X2)
SS_SOC, # 사회적지지: 사교적 평균 (X3)
MDEP, # 어머니 우울 평균 (M)
PAR_SOCIAL_M) # 양육행동 사회적 평균 (Y)
# ────────────────────────────────────────────────────────────
# STEP 3-1. 모형 수정 및 추정
# ────────────────────────────────────────────────────────────
model_rev <- '
# ── 구조 방정식 (회귀식) ──────────────────────────────────
MDEP ~ a1 * SS_EMO # X1 → M
MDEP ~ a2 * SS_INS # X2 → M
MDEP ~ a3 * SS_SOC # X3 → M
PAR_SOCIAL_M ~ b * MDEP # M → Y
PAR_SOCIAL_M ~ SS_SOC #사교적 지지 ->양육행동 추가
# ── 간접효과 정의 ─────────────────────────────────────────
ind_EMO := a1 * b # SS_EMO → MDEP → PAR_SOCIAL_M
ind_INS := a2 * b # SS_INS → MDEP → PAR_SOCIAL_M
ind_SOC := a3 * b # SS_SOC → MDEP → PAR_SOCIAL_M
'
fit_rev <- sem(model_rev,
data = df_sel,
estimator = "ML",
missing = "fiml",
fixed.x = FALSE)
## Warning: lavaan->lav_data_full():
## some cases are empty and will be ignored: 16 17 18 21 28 47 56 64 67 81 99
## 105 108 110 113 125 127 129 136 137 161 166 173 175 182 198 204 224 225
## 227 235 252 261 262 266 282 294 298 303 311 319 320 324 330 337 345 348
## 349 353 361 366 370 379 388 389 392 405 419 423 439 445 467 468 480 486
## 489 495 503 521 535 547 551 559 583 587 590 610 612 630 632 634 653 671
## 674 676 693 699 700 710 721 722 731 747 756 763 765 771 772 780 781 782
## 787 790 793 829 836 840 851 853 864 869 882 886 892 893 897 898 902 914
## 916 919 922 924 928 930 931 937 939 941 948 952 973 976 978 979 1054 1081
## 1123 1128 1146 1151 1162 1173 1182 1192 1241 1282 1299 1301 1306 1307 1310
## 1316 1329 1338 1350 1358 1360 1363 1375 1378 1390 1394 1402 1407 1409 1414
## 1418 1422 1433 1435 1440 1444 1451 1454 1467 1489 1532 1538 1539 1542 1552
## 1555 1562 1564 1575 1580 1582 1589 1592 1604 1606 1607 1611 1612 1619 1623
## 1649 1656 1667 1673 1684 1688 1690 1700 1705 1708 1728 1737 1738 1748 1750
## 1765 1776 1791 1794 1820 1831 1837 1842 1846 1860 1870 1877 1882 1888 1895
## 1904 1905 1945 1948 1951 1954 1958 1960 1963 1976 1978 1982 1989 1991 2006
## 2015 2017 2050 2056 2059 2088 2095.
summary(fit_rev,
fit.measures = TRUE,
standardized = TRUE,
rsquare = TRUE)
## lavaan 0.6-21 ended normally after 40 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 18
##
## Used Total
## Number of observations 1901 2150
## Number of missing patterns 9
##
## Model Test User Model:
##
## Test statistic 26.612
## Degrees of freedom 2
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 593.915
## Degrees of freedom 7
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.958
## Tucker-Lewis Index (TLI) 0.853
##
## Robust Comparative Fit Index (CFI) 0.958
## Robust Tucker-Lewis Index (TLI) 0.852
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7858.705
## Loglikelihood unrestricted model (H1) -7845.399
##
## Akaike (AIC) 15753.410
## Bayesian (BIC) 15853.312
## Sample-size adjusted Bayesian (SABIC) 15796.126
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.080
## 90 Percent confidence interval - lower 0.055
## 90 Percent confidence interval - upper 0.109
## P-value H_0: RMSEA <= 0.050 0.026
## P-value H_0: RMSEA >= 0.080 0.550
##
## Robust RMSEA 0.081
## 90 Percent confidence interval - lower 0.055
## 90 Percent confidence interval - upper 0.110
## P-value H_0: Robust RMSEA <= 0.050 0.026
## P-value H_0: Robust RMSEA >= 0.080 0.561
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.019
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## MDEP ~
## SS_EMO (a1) -0.181 0.029 -6.195 0.000 -0.181 -0.200
## SS_INS (a2) -0.075 0.035 -2.171 0.030 -0.075 -0.076
## SS_SOC (a3) -0.082 0.035 -2.351 0.019 -0.082 -0.082
## PAR_SOCIAL_M ~
## MDEP (b) -0.188 0.015 -12.319 0.000 -0.188 -0.268
## SS_SOC 0.176 0.015 11.511 0.000 0.176 0.251
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## SS_EMO ~~
## SS_INS 0.385 0.016 24.806 0.000 0.385 0.694
## SS_SOC 0.375 0.015 24.538 0.000 0.375 0.683
## SS_INS ~~
## SS_SOC 0.371 0.014 25.897 0.000 0.371 0.741
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .MDEP 3.225 0.092 34.978 0.000 3.225 4.591
## .PAR_SOCIAL_M 3.577 0.074 48.479 0.000 3.577 7.252
## SS_EMO 3.836 0.018 214.186 0.000 3.836 4.921
## SS_INS 3.907 0.016 238.702 0.000 3.907 5.486
## SS_SOC 3.875 0.016 239.507 0.000 3.875 5.503
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .MDEP 0.441 0.014 30.730 0.000 0.441 0.894
## .PAR_SOCIAL_M 0.201 0.007 30.757 0.000 0.201 0.828
## SS_EMO 0.608 0.020 30.767 0.000 0.608 1.000
## SS_INS 0.507 0.016 30.761 0.000 0.507 1.000
## SS_SOC 0.496 0.016 30.755 0.000 0.496 1.000
##
## R-Square:
## Estimate
## MDEP 0.106
## PAR_SOCIAL_M 0.172
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind_EMO 0.034 0.006 5.531 0.000 0.034 0.054
## ind_INS 0.014 0.007 2.138 0.033 0.014 0.020
## ind_SOC 0.015 0.007 2.310 0.021 0.015 0.022
# ────────────────────────────────────────────────────────────
# STEP 4-1. 수정모형 평가
# ────────────────────────────────────────────────────────────
fit_rev_idx <- fitMeasures(fit_rev, c("chisq", "df", "pvalue",
"cfi", "tli",
"rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
"srmr"))
cat("\n=== 모형 적합도 지수 ===\n")
##
## === 모형 적합도 지수 ===
print(round(fit_rev_idx, 3))
## chisq df pvalue cfi tli
## 26.612 2.000 0.000 0.958 0.853
## rmsea rmsea.ci.lower rmsea.ci.upper srmr
## 0.080 0.055 0.109 0.019
# ────────────────────────────────────────────────────────────
# STEP 5. 경로계수
# ────────────────────────────────────────────────────────────
params_rev <- parameterEstimates(fit_rev, standardized = TRUE)
direct_paths_rev <- params_rev[params_rev$op == "~", ]
cat("\n=== 직접 경로계수 ===\n")
##
## === 직접 경로계수 ===
print(direct_paths_rev[, c("lhs","op","rhs","est","se","z",
"pvalue","ci.lower","ci.upper","std.all")])
## lhs op rhs est se z pvalue ci.lower ci.upper std.all
## 1 MDEP ~ SS_EMO -0.181 0.029 -6.195 0.000 -0.238 -0.123 -0.200
## 2 MDEP ~ SS_INS -0.075 0.035 -2.171 0.030 -0.143 -0.007 -0.076
## 3 MDEP ~ SS_SOC -0.082 0.035 -2.351 0.019 -0.150 -0.014 -0.082
## 4 PAR_SOCIAL_M ~ MDEP -0.188 0.015 -12.319 0.000 -0.218 -0.158 -0.268
## 5 PAR_SOCIAL_M ~ SS_SOC 0.176 0.015 11.511 0.000 0.146 0.206 0.251
# ────────────────────────────────────────────────────────────
# STEP 6. Bootstrap 간접효과
# ────────────────────────────────────────────────────────────
set.seed(1234)
fit_rev_boot <- sem(model_rev,
data = df_sel,
estimator = "ML",
missing = "FIML",
se = "bootstrap",
bootstrap = 5000)
## Warning: lavaan->lav_data_full():
## 265 cases were deleted due to missing values in exogenous variable(s),
## while fixed.x = TRUE.
boot_result_rev <- parameterEstimates(fit_rev_boot,
boot.ci.type = "bca.simple",
level = 0.95)
boot_indir_rev <- boot_result_rev[boot_result_rev$op == ":=", ]
cat("\n=== 간접효과 Bootstrap 95% CI ===\n")
##
## === 간접효과 Bootstrap 95% CI ===
print(boot_indir_rev[, c("lhs","est","se","pvalue","ci.lower","ci.upper")])
## lhs est se pvalue ci.lower ci.upper
## 19 ind_EMO 0.034 0.006 0.000 0.022 0.047
## 20 ind_INS 0.014 0.007 0.059 0.000 0.029
## 21 ind_SOC 0.015 0.007 0.052 0.000 0.030