② 초미세먼지 지역간 이동 PSM 분석

10m 역궤적 기반 당진·서산 경유 여부에 따른 화성 PM2.5 차이 검증

1 개요

활용한 데이터는 모두 공개 자료로부터 추출한 것이다. 그 과정은 너무 복잡해서 여기서 설명하기 어렵지만, 화성시 8개 측정소의 시간자료에, 도착 지점까지의 10m 고도 역궤적이 특정 지자체(당진·서산)를 경유했는지 여부를 결합한 데이터로, “경유한 날의 PM2.5가 경유하지 않은 날보다 유의미하게 높은가” 를 검증한다.

단순 평균 비교(t-test)로는 이를 제대로 볼 수 없다. 경유하지 않은 날의 수가 경유한 날보다 월등히 많아, 실질적으로 미미한 차이도 통계적으로 유의하게 나올 수 있기 때문이다. 그래서 날씨 조건이 유사한 “쌍둥이 날”끼리 짝지어 비교하는 성향점수 짝짓기 (PSM) 분석방법을 사용해야 한다. 실험설계 조건은 아래와 같다.

노트네 가지 실험설계

케이스	처치(경유)	대조군
Case1	당진 only	무경유(neither)
Case2	당진 any (only+both)	무경유(neither)
Case3	서산 any (only+both)	무경유(neither)
Case4	통합 treated any	무경유(neither)

2 데이터 로드와 황사일 제외

hs <- read.csv("./data/psm_dataset_full.csv")
dplyr::glimpse(hs)   # 전체 250,281건 · 32개 변수

Rows: 250,281
Columns: 32
$ station_code      <int> 131551, 131551, 131551, 131551, 131551, 131551, 1315…
$ pm25_obs          <dbl> 11, 9, 10, 9, 12, 7, 6, 9, 11, 19, 14, 14, 10, 11, 1…
$ u_aws             <dbl> 1.3972, -0.1282, -0.3488, -0.0861, 0.0000, -0.2072, …
$ v_aws             <dbl> -0.0879, -0.3789, -0.8297, -0.1805, 0.0000, -0.3421,…
$ u_aws_grid_avg    <dbl> 0.4448, 0.1238, -0.0702, 0.0693, 0.2659, 0.0688, -0.…
$ v_aws_grid_avg    <dbl> -0.0454, -0.1747, -0.2609, -0.0265, -0.1390, -0.1839…
$ mixdepth          <dbl> 572.8, 523.1, 473.4, 413.1, 352.7, 292.2, 276.6, 261…
$ hm                <dbl> 48, 47, 49, 59, 62, 70, 70, 68, 53, 44, 36, 32, 32, …
$ si_w_mean_24h     <dbl> 128.2407, 128.2407, 128.2407, 128.2407, 128.2407, 12…
$ arrival_only_aod  <dbl> 0.0778670, 0.0778670, 0.0852834, 0.0926997, 0.092699…
$ falt10_pm25_avg   <dbl> 10.7593, 9.3588, 9.4780, 9.9565, 11.1737, 10.7528, 1…
$ falt1000_aod_avg  <dbl> 0.0688176, 0.0717999, 0.0743557, 0.0789165, 0.083262…
$ str_kst           <int> 2022010101, 2022010102, 2022010103, 2022010104, 2022…
$ str_utc           <int> 2021123116, 2021123117, 2021123118, 2021123119, 2021…
$ group             <chr> "neither", "neither", "neither", "neither", "neither…
$ seosan_frac       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ dangjin_frac      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ intersect_seosan  <chr> "False", "False", "False", "False", "False", "False"…
$ intersect_dangjin <chr> "False", "False", "False", "False", "False", "False"…
$ n_points          <int> 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, …
$ pm25_traj         <dbl> 11, 9, 10, 9, 12, 7, 6, 9, 11, 19, 14, 14, 10, 11, 1…
$ has_trajectory    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ is_seosan         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ is_dangjin        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ is_both           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ neither           <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ seosan_any        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ dangjin_any       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ treated_any       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ is_dust_day       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ year              <int> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022…
$ year_month        <int> 202201, 202201, 202201, 202201, 202201, 202201, 2022…

황사일은 지자체가 통제할 수 없는 불가항력적 외부 현상이므로 분석에서 제외한다.

hs |>
  group_by(is_dust_day) |>
  summarize(count = n(), .groups = "drop")   # 황사일 10,697건

# A tibble: 2 × 2
  is_dust_day  count
        <int>  <int>
1           0 239584
2           1  10697

hs_no_dust <- hs |> filter(is_dust_day == 0)

dim(hs_no_dust)                     # 239,584건

[1] 239584     32

sum(complete.cases(hs_no_dust))     # 결측 없는 완전한 데이터

[1] 239584

실제 데이터의 구성을 다시 한번 확인해 본다.

hs_no_dust |>
  group_by(group) |>
  summarize(count = n(), .groups = "drop")

# A tibble: 4 × 2
  group         count
  <chr>         <int>
1 both          36263
2 dangjin_only  18704
3 neither      182412
4 seosan_only    2205

3 예비 분석: 다중회귀로 주요 영향변수 파악

본격적인 PSM에 앞서, PM2.5에 영향을 주는 요인들을 한 번에 저울에 올려 상대적 영향력과 다중공선성을 확인한다.

노트주요 활용 변수 개요

변수명	의미	비고
falt10_pm25_avg	10m 역궤적이 지난 시공간의 초미세먼지 농도 평균	데이터 누출(leakage)에 주의 필요
falt1000_aod_avg	1000m 역궤적이 지나간 시공간의 에어로졸 광학두께 평균	상동
mixdepth	혼합고	초미세먼지 농도와 반비례하는 관계임
si_w_mean_24h	일사량의 24시간 이동 평균	야간시간대 보간을 위해 24시간 이동평균 활용
hm	습도	지역을 대표하는 위치에 있는 종관기상 관측소(ASOS) 데이터 활용
u_aws	풍향/풍속의 동서(東西) 성분	관측소에서 가장 가까운 방재기상 관측소(AWS) 데이터 활용
v_aws	풍향/풍속의 남북(南北) 성분	상동

hs_pm25_reg <- lm(
  pm25_obs ~ u_aws + v_aws + mixdepth + hm + si_w_mean_24h +
             arrival_only_aod + falt10_pm25_avg + falt1000_aod_avg,
  data = hs_no_dust
)
summary(hs_pm25_reg)   # 결정계수(R²) ≈ 0.72

Call:
lm(formula = pm25_obs ~ u_aws + v_aws + mixdepth + hm + si_w_mean_24h + 
    arrival_only_aod + falt10_pm25_avg + falt1000_aod_avg, data = hs_no_dust)

Residuals:
    Min      1Q  Median      3Q     Max 
-70.143  -4.625  -0.688   3.773 129.545 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -1.585e+00  1.202e-01  -13.18   <2e-16 ***
u_aws             1.665e-01  1.209e-02   13.77   <2e-16 ***
v_aws             1.530e-01  1.453e-02   10.53   <2e-16 ***
mixdepth         -4.208e-03  5.544e-05  -75.90   <2e-16 ***
hm                4.124e-02  1.126e-03   36.61   <2e-16 ***
si_w_mean_24h     1.096e-02  2.287e-04   47.92   <2e-16 ***
arrival_only_aod  1.081e+00  9.003e-02   12.00   <2e-16 ***
falt10_pm25_avg   1.020e+00  1.447e-03  704.86   <2e-16 ***
falt1000_aod_avg  2.133e+00  1.106e-01   19.28   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.402 on 239575 degrees of freedom
Multiple R-squared:  0.7188,    Adjusted R-squared:  0.7188 
F-statistic: 7.655e+04 on 8 and 239575 DF,  p-value: < 2.2e-16

단위가 서로 다른 요인들의 영향력을 비교하기 위해 표준화 계수를 구한다.

standardize_parameters(hs_pm25_reg)

# Standardization method: refit

Parameter        | Std. Coef. |         95% CI
----------------------------------------------
(Intercept)      |  -5.04e-15 | [ 0.00,  0.00]
u aws            |       0.02 | [ 0.01,  0.02]
v aws            |       0.01 | [ 0.01,  0.01]
mixdepth         |      -0.11 | [-0.11, -0.10]
hm               |       0.05 | [ 0.05,  0.05]
si w mean 24h    |       0.06 | [ 0.05,  0.06]
arrival only aod |       0.02 | [ 0.02,  0.02]
falt10 pm25 avg  |       0.81 | [ 0.81,  0.81]
falt1000 aod avg |       0.03 | [ 0.03,  0.03]

vif(hs_pm25_reg)   # 통상 10 이상이면 다중공선성 의심

           u_aws            v_aws         mixdepth               hm 
        1.169167         1.087349         1.644898         1.689297 
   si_w_mean_24h arrival_only_aod  falt10_pm25_avg falt1000_aod_avg 
        1.193321         2.105500         1.126975         2.030046 

4 성향점수매칭 (Case1: 당진 only)

먼저 Case1(당진 only vs 무경유)로 매칭 절차와 품질 진단을 자세히 본다. 당진·서산 both와 서산 only를 제외해 당진 only를 표상하는 변수(is_dangjin)에 해당하는 데이터만 남도록 전처리한 후 분석에 사용한다.

중요PSM 공변량 선택 기준

성향점수 공변량으로는 순수 기상변수(u_aws, v_aws, mixdepth, hm, si_w_mean_24h)만 사용한다. arrival_only_aod, falt10_pm25_avg, falt1000_aod_avg는 처치(경유)의 결과에 가까운 사후(post-treatment) 성격이라 공변량에서 제외한다. 결과 성격의 변수를 넣으면 편향이 생기기 때문이다.

hs_case1 <- hs_no_dust |>
  filter(is_both != 1 & is_seosan != 1)
dim(hs_case1)

[1] 201116     32

m1.out <- matchit(
  is_dangjin ~ u_aws + v_aws + mixdepth + hm + si_w_mean_24h,
  data = hs_case1, method = "nearest", distance = "glm"
)
summary(m1.out)   # 매칭 전후 balance

Call:
matchit(formula = is_dangjin ~ u_aws + v_aws + mixdepth + hm + 
    si_w_mean_24h, data = hs_case1, method = "nearest", distance = "glm")

Summary of Balance for All Data:
              Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance             0.1615        0.0860          0.5815     3.4281    0.2325
u_aws                0.3004        0.2357          0.0477     0.7727    0.0148
v_aws                0.6603       -0.1526          0.6370     1.2286    0.1298
mixdepth           409.0972      491.7705         -0.2396     0.6970    0.0446
hm                  74.0501       67.8241          0.3310     0.9067    0.0681
si_w_mean_24h      180.7542      157.7917          0.2665     1.1911    0.0695
              eCDF Max
distance        0.3377
u_aws           0.0366
v_aws           0.2506
mixdepth        0.0751
hm              0.1459
si_w_mean_24h   0.1551

Summary of Balance for Matched Data:
              Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance             0.1615        0.1576          0.0298     1.2665    0.0001
u_aws                0.3004        0.2893          0.0082     0.7900    0.0160
v_aws                0.6603        0.6109          0.0387     1.2783    0.0139
mixdepth           409.0972      401.7473          0.0213     0.8502    0.0155
hm                  74.0501       74.0893         -0.0021     1.0497    0.0052
si_w_mean_24h      180.7542      184.2909         -0.0410     1.0580    0.0112
              eCDF Max Std. Pair Dist.
distance        0.0124          0.0298
u_aws           0.0405          1.0940
v_aws           0.0277          0.5743
mixdepth        0.0319          1.0508
hm              0.0173          1.0258
si_w_mean_24h   0.0470          1.0753

Sample Sizes:
          Control Treated
All        182412   18704
Matched     18704   18704
Unmatched  163708       0
Discarded       0       0

4.1 매칭 품질 진단 (balance)

매칭의 성패는 처치효과가 아니라 공변량이 두 집단에서 얼마나 고르게 맞았는가로 판단한다. love plot에서 매칭 후(파란 점) 모든 변수가 ±0.1 기준선 안으로 들어오면 성공이다.

love.plot(m1.out, thresholds = 0.1)

그림 1: 공변량 balance(love plot). 파란 점이 모두 ±0.1 안 → 날씨 조건이 균형있게 맞춰짐

개별 변수의 분포가 매칭 전후로 얼마나 겹치는지, 성향점수 공통지지(common support)와 함께 확인한다.

bal.plot(m1.out, var.name = "mixdepth", which = "both") /
  bal.plot(m1.out, var.name = "hm", which = "both") /
  bal.plot(m1.out, var.name = "distance", which = "both")

그림 2: 매칭 전후 분포 겹침. distance = 성향점수 공통지지 확인

4.2 처치효과(ATT) 추정

매칭된 데이터에서 처치효과를 추정한다. 매칭쌍(subclass) 구조를 반영한 클러스터-강건 표준오차로 유의성을 본다.

m1_data <- match.data(m1.out)

fit <- lm(pm25_obs ~ is_dangjin, data = m1_data, weights = weights)

# 매칭쌍 클러스터 보정 표준오차
coeftest(fit, vcov. = vcovCL, cluster = ~subclass)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 20.00481    0.11534 173.443 < 2.2e-16 ***
is_dangjin   2.53042    0.16373  15.454 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

노트

matchit의 기본 추정 대상은 ATT(처치집단에 대한 평균 처치효과)이다. 즉 “실제로 당진을 거쳐온 날들에서” 당진 경유가 PM2.5를 얼마나 높였는가를 뜻한다.

5 네 가지 실험설계 자동화

Case2~4는 “처치변수 이름”과 “제외 필터”만 다르므로 함수로 묶어 일괄 처리한다.

run_psm <- function(data, treat, drop_expr) {
  d  <- data |> filter({{ drop_expr }})
  f  <- as.formula(paste(treat, "~ u_aws + v_aws + mixdepth + hm + si_w_mean_24h"))
  m  <- matchit(f, data = d, method = "nearest", distance = "glm")  # balance 우수 → caliper 불필요
  md <- match.data(m)
  fit <- lm(pm25_obs ~ md[[treat]], data = md, weights = weights)
  ct <- coeftest(fit, vcov. = vcovCL, cluster = ~subclass)
  list(match = m, fit = fit,
       att   = unname(coef(fit)[2]),      # 처치효과(㎍/㎥)
       se    = ct[2, "Std. Error"],
       base  = unname(coef(fit)[1]))      # 대조군 평균
}

# Case1: 당진 only vs neither
r1 <- run_psm(hs_no_dust, "is_dangjin",  is_both != 1 & is_seosan != 1)
# Case2: 당진 any vs neither (both 포함, 서산 only 제외)
r2 <- run_psm(hs_no_dust, "dangjin_any", is_seosan != 1)
# Case3: 서산 any vs neither (both 포함, 당진 only 제외)
r3 <- run_psm(hs_no_dust, "seosan_any",  is_dangjin != 1)
# Case4: treated any vs neither
r4 <- run_psm(hs_no_dust, "treated_any", TRUE)

6 결과 종합

네 케이스의 처치효과를 한 표로 모은다. 95% 신뢰구간은 ATT ± 1.96 × SE로 계산한다.

results <- data.frame(
  경로       = c("당진 only", "당진 any", "서산 any", "통합(treated any)"),
  대조군평균 = round(c(r1$base, r2$base, r3$base, r4$base), 2),
  ATT        = round(c(r1$att, r2$att, r3$att, r4$att), 2),
  SE         = round(c(r1$se, r2$se, r3$se, r4$se), 3)
)
results$CI_하한 <- round(results$ATT - 1.96 * results$SE, 2)
results$CI_상한 <- round(results$ATT + 1.96 * results$SE, 2)
results$상승률  <- round(results$ATT / results$대조군평균 * 100, 1)
saveRDS(results, "./data/results.rds")
results

               경로 대조군평균  ATT    SE CI_하한 CI_상한 상승률
1         당진 only      20.00 2.53 0.164    2.21    2.85   12.7
2          당진 any      20.30 1.37 0.094    1.19    1.55    6.7
3          서산 any      20.29 1.08 0.110    0.86    1.30    5.3
4 통합(treated any)      20.36 1.39 0.092    1.21    1.57    6.8

ggplot(results, aes(x = reorder(경로, ATT), y = ATT)) +
  geom_col(width = 0.6, fill = "#4C78A8") +
  geom_errorbar(aes(ymin = CI_하한, ymax = CI_상한), width = 0.2) +
  geom_text(aes(label = sprintf("+%.2f", ATT)), hjust = -0.4, vjust=-1.2) +
  coord_flip() +
  labs(x = NULL, y = "PM2.5 상승분 (㎍/㎥)",
       title = "경로별 화성 초미세먼지 기여도(ATT)") +
  theme_minimal(base_size = 14)

그림 3: 경로별 화성 PM2.5 상승분(ATT). 오차막대 = 95% 신뢰구간

7 소결

• 네 실험설계 모두에서 경유일의 PM2.5가 무경유일보다 통계적으로 확실하게 높다.
• 특히 당진 only(+2.53 ㎍/㎥) 의 상승분이 가장 크다.
• 매칭으로 화성 지점의 바람(u/v)을 통제했으므로, 이 차이는 “그날 바람 탓”이 아니라 경로 자체의 효과로 해석할 수 있다.
• 결과에 대한 해석과 편익 분석내용은 ③ 결과 해석 및 편익 분석을 참조하기 바란다.