제 9강: 통계적 토대

데이터과학 입문

Author

Affiliation

원중호

서울대학교 통계학과

Published

April 2024

시작하기 전에

다음의 패키지가 설치되어 있지 않으면 설치한다.

# install.packages("tidyverse")
# install.packages("mdsr")
# install.packages("nycflights13")
library(tidyverse)
library(mdsr)
library(nycflights13)
sessionInfo()

R version 4.3.3 (2024-02-29)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Asia/Seoul
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] nycflights13_1.0.2 mdsr_0.2.7         lubridate_1.9.3    forcats_1.0.0     
 [5] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
 [9] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.0      tidyverse_2.0.0   

loaded via a namespace (and not attached):
 [1] gtable_0.3.4      jsonlite_1.8.8    compiler_4.3.3    tidyselect_1.2.0 
 [5] scales_1.3.0      yaml_2.3.8        fastmap_1.1.1     R6_2.5.1         
 [9] generics_0.1.3    skimr_2.1.5       knitr_1.45        htmlwidgets_1.6.4
[13] munsell_0.5.0     pillar_1.9.0      tzdb_0.4.0        rlang_1.1.3      
[17] utf8_1.2.4        stringi_1.8.3     repr_1.1.6        xfun_0.42        
[21] timechange_0.3.0  cli_3.6.2         withr_3.0.0       magrittr_2.0.3   
[25] digest_0.6.34     grid_4.3.3        rstudioapi_0.15.0 base64enc_0.1-3  
[29] hms_1.1.3         lifecycle_1.0.4   vctrs_0.6.5       evaluate_0.23    
[33] glue_1.7.0        fansi_1.0.6       colorspace_2.1-0  rmarkdown_2.25   
[37] tools_4.3.3       pkgconfig_2.0.3   htmltools_0.5.7

들어가며

데이터과학의 궁극적 목표: 데이터에서 의미 추출
- 랭글링: 데이터를 해석하기 쉽도록 사례와 변수를 재구성
- 시각화: 의미를 탐색할 수 있도록 사람의 마음과 데이터 연결
시각화와 인간의 시각적 인지 능력
- 인간은 임의의 노이즈에 의해 부분적으로 가려져 있어도 패턴을 보는 데 매우 능숙
- 패턴이 없는 경우에도 패턴을 보는 데 매우 능숙 — 노이즈에 나타나는 우연적 패턴에 쉽게 현혹 (별자리?)

통계적 방법
- 패턴이 강하고 견고해서 단순한 우연이 아니라고 확신할 수 있을 때 식별
- 통계적 방법은 패턴과 그 강도를 정량화
- 복잡하거나 다면적이어서 시각적으로 확인하기 어려운 패턴을 찾는 데도 매우 중요

표본과 모집단

이전 장에서는 데이터를 고정된 것으로 간주
통계학에서는 우리가 가지고 있는 데이터는 고정되어 있지만 방법론은 훨씬 더 큰 잠재적 사례 집합(모집단)에서 표본을 추출한다고 가정
모집단에서 추출되었을 수 있는 다른 표본을 상상
모집단에서 측정할 수 있는 추가 변수가 있었을 수도 있다고 상상
- 이 변수와 관련하여 나타나는 모든 패턴은 무작위적이고 우연적인 것이어야 한다
도구: 확률

모집단에서의 표본 추출

뉴욕에 거주하는 사람의 샌프란시스코 출장을 위한 출장 계획 세우기
- 특정 시각에 회의가 있다고 가정

이 시각 보다 얼마나 일찍 도착해야 항공편 지연으로 인해 회의에 늦지 않을 수 있을까?

nycflights13 패키지의 336,776개 항공편을 모집단으로 가정

샌프란시스코 공항(SFO)로 들어오는 항공편의 모집단에서 샘플을 추출하여 수집된 데이터(표본)를 기반으로 출장 계획을 세우는 상황을 시뮬레이션 (표본 크기 \(n\)=25)

SF <- flights %>%
  filter(dest == "SFO", !is.na(arr_delay))
set.seed(101)
sf_25 <- SF %>%
  slice_sample(n = 25)

간단한 방법: 가장 긴 항공편 지연을 찾아 이 지연에 대처할 수 있도록 여행 일정 조정

sf_25 %>% skim(arr_delay)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
arr_delay	25	0	0.12	29.79	-38	-23	-5	14	103

최대 지연은 103분. 이만큼 일찍 도착하면 될까?

SF %>% skim(arr_delay)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
arr_delay	13173	0	2.67	47.67	-86	-23	-8	12	1007

모집단에서의 최대 지연 시간은 1007분.
17시간 정도 일찍 도착?
- 하루를 더 여행하면 숙박, 식사 등의 추가비용 발생
- 1007분 이상 지연이 발생하지 않는다는 보장 없음

합리적인 계획: 비용 절약과 지각할 가능성을 절충 — 100번 같은 출장을 간다면 그 중 2번만 지각해도 괜찮음

  sf_25 %>% summarize(q98 = quantile(arr_delay, p = 0.98))

# A tibble: 1 × 1
    q98
  <dbl>
1  67.5

90분 일찍 도착하는 계획을 세웠다면 회의에 2% 이하의 확률로만 지각하겠다는 의도를 달성하는 데 얼마나 효과가 있었을까?

SF %>%
  group_by(arr_delay < 90) %>%
  count() %>%
  mutate(pct = n / nrow(SF))

# A tibble: 2 × 3
# Groups:   arr_delay < 90 [2]
  `arr_delay < 90`     n    pct
  <lgl>            <int>  <dbl>
1 FALSE              640 0.0486
2 TRUE             12533 0.951

90분 정책은 훨씬 더 나쁜 5%의 빈도로 목표를 달성하지 못함

2%의 확률을 정확하게 맞추려면 정책을 90분에서 어떤 값으로 늘려야 할까?

SF %>% summarize(q98 = quantile(arr_delay, p = 0.98))

# A tibble: 1 × 1
    q98
  <dbl>
1   153

대부분의 실제 환경에서는 모집단 데이터에 접근할 수 없고 표본 데이터만 있을 뿐
샘플에서 얻은 결과가 98% 목표를 충족하기에 충분한지 판단하기 위해 샘플을 어떻게 사용할 수 있을까?
충분하지 않다면 표본의 크기는 어느 정도여야 만족할 만한 결과를 얻을 수 있을까?

표본통계량

표본분포

항공편 지연 시간 표본의 98백분위수는 (별로 좋지 않은) 통계량
모집단에서 새로운 표본을 추출하는 경우, 새로운 표본의 통계량이 원래 표본에서 계산된 동일한 통계량과 얼마나 유사할까?
즉, 모집단에서 크기 \(n\)인 표본을 여러 번 추출하고 각 표본에 대한 통계량을 계산하면 모든 이들 통계량은 얼마나 유사할까?

n <- 25
num_trials <- 500
sf_25_means <- 1:num_trials %>%
  map_dfr(
    ~ SF %>%
      slice_sample(n = n) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)

head(sf_25_means)

# A tibble: 6 × 2
  mean_arr_delay     n
           <dbl> <dbl>
1          -5.92    25
2           9.24    25
3          -6.28    25
4          14.2     25
5          -2.84    25
6           0.6     25

sf_25_means %>% skim(mean_arr_delay)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
mean_arr_delay	500	0	3.09	10.16	-20.12	-3.72	1.54	9	56.68

평균이 중앙값의 두 배 이상이므로 표본분포가 오른쪽으로 치우쳐(right-skewed) 있음
다른 표본통계량

sf_25_means %>%
  summarize(
    x_bar = mean(mean_arr_delay),
    se = sd(mean_arr_delay)
  ) %>%
  mutate(
    ci_lower = x_bar - 2 * se, # approximately 95% of observations 
    ci_upper = x_bar + 2 * se  # are within two standard errors
  )

# A tibble: 1 × 4
  x_bar    se ci_lower ci_upper
  <dbl> <dbl>    <dbl>    <dbl>
1  3.09  10.2    -17.2     23.4

표본 크기

신뢰성 있는 결과를 얻으려면 얼마나 큰 표본이 필요한가?

n <- 100
sf_100_means <- 1:500 %>%
  map_dfr(
    ~ SF %>%
      slice_sample(n = n) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)

sf_25_means %>%
  bind_rows(sf_100_means) %>%
  ggplot(aes(x = mean_arr_delay)) + 
  geom_histogram(bins = 30) + 
  facet_grid( ~ n) + 
  xlab("Sample mean")

두 표본분포는 중심은 거의 일치함
표본 크기가 클수록 표준 오차가 더 작아짐. 즉, 표본 크기가 클수록 작은 표본 크기보다 더 신뢰할 수 있음
표본 크기가 큰 경우, 표본분포의 모양은 종 모양이 되는 경향이 있음

부트스트랩

재표집과 부트스트랩

앞 예제에서는 모집단 데이터에 접근할 수 있었기 때문에 모집단에서 반복적으로 표본을 추출하여 표본분포를 구할 수 있었음
현실에서는 전체 모집단이 아닌 하나의 표본만 있음
부트스트랩은 모집단에 대한 접근 없이도 표본분포를 근사할 수 있는 통계적 방법
표본 자체를 모집단처럼 생각, 복원추출을 통해 “새로운 표본”을 많이 생성, 각 표본마다 통계량 계산

표본

three_flights <- SF %>%
  slice_sample(n = 3, replace = FALSE) %>%
  select(year, month, day, dep_time)
three_flights

# A tibble: 3 × 4
   year month   day dep_time
  <int> <int> <int>    <int>
1  2013     8     5     1327
2  2013     3     3     1526
3  2013     4    11      728

재표집

three_flights %>% slice_sample(n = 3, replace = TRUE)

# A tibble: 3 × 4
   year month   day dep_time
  <int> <int> <int>    <int>
1  2013     4    11      728
2  2013     3     3     1526
3  2013     8     5     1327

n <- 200
orig_sample <- SF %>%        # sample
  slice_sample(n = n, replace = TRUE) 
# resample
sf_200_bs <- 1:num_trials %>%    # num_trials == 500
  map_dfr(
    ~orig_sample %>%
      slice_sample(n = n, replace = TRUE) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)
# bootstrap statistics
sf_200_bs %>%
  skim(mean_arr_delay)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
mean_arr_delay	500	0	8.52	4.04	-2.04	5.79	8.42	11.28	24.11

이 예제에서는 모집단을 알고 있으므로 부트스트랩 추정량을 크기 \(n\)=200인 가상의 표본들과 비교

sf_200_pop <- 1:num_trials %>%
  map_dfr(
    ~SF %>%
      slice_sample(n = n, replace = FALSE) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)

sf_200_pop %>%
  skim(mean_arr_delay)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
mean_arr_delay	500	0	2.67	3.3	-6.24	0.38	2.5	4.9	17.22

부트스트랩(sf_200_bs)에는 모집단이 사용되지 않고 원 표본만 사용되었다는 점에 주의
부트스트랩을 사용하여 계산된 표준오차(3.1분)가 모집단에서 반복적으로 표본을 추출하여 계산된 표본분포의 표준 오차(3.3분)에 대한 괜찮은 근사치
부트스트랩 분포 = 부트스트랩 통계량의 분포
- 표본 분포와 정확히 같지는 않지만, 표본 크기가 적당히 크고 부트스트랩 시행 횟수가 충분한 경우 표준 오차 및 분위수 등 중요한 표본 분포의 특성에 수렴 (Efron and Tibshirani 1993).

출장 예제

크기 \(n\)=200인 표본으로부터 98백분위를 계산

orig_sample %>% summarize(q98 = quantile(arr_delay, p = 0.98))

# A tibble: 1 × 1
    q98
  <dbl>
1  209.

부트스트랩을 이용하여 신뢰도를 계산

n <- nrow(orig_sample)
sf_200_bs <- 1:num_trials %>%
  map_dfr(
    ~orig_sample %>%
      slice_sample(n = n, replace = TRUE) %>%
      summarize(q98 = quantile(arr_delay, p = 0.98))
  )

sf_200_bs %>% skim(q98)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
q98	500	0	206.36	57.18	75	153	208.86	251.02	344.32

부트스트랩 표준 오차는 약 29분 — 95% 신뢰구간은 대략 140 ± 58분
140분 일찍 가기 계획은 그리 좋아 보이지 않음

표본 크기를 늘리면 나을까?

set.seed(1001)
n_large <- 10000
sf_10000_bs <- SF %>% 
  slice_sample(n = n_large, replace = FALSE)

sf_200_bs <- 1:num_trials %>%
  map_dfr(~sf_10000_bs %>%
        slice_sample(n = n_large, replace = TRUE) %>%
        summarize(q98 = quantile(arr_delay, p = 0.98))
  )

sf_200_bs %>% skim(q98)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
q98	500	0	153.88	4.14	139.02	151.02	153.02	156.02	169

154 ± 8분은 훨씬 좋아 보임
데이터가 많을수록 추정량을 모수에 더 가깝게 할 수 있으며, 특히 꼬리 부분의 추정이 더 쉬워짐

이상치

이상치 식별

더 많은 데이터가 도움이 되는 분야 중 하나
항공편이 7시간(420분) 이상 지연되는 경우를 극단적인 상황으로 간주해 보자
- 420분은 자의적인 선택이나 심각한 지연을 나타내는 지표로 유용할 수 있음

SF %>%
  filter(arr_delay >= 420) %>% 
  select(month, day, dep_delay, arr_delay, carrier)

# A tibble: 7 × 5
  month   day dep_delay arr_delay carrier
  <int> <int>     <dbl>     <dbl> <chr>  
1    12     7       374       422 UA     
2     7     6       589       561 DL     
3     7     7       629       676 VX     
4     7     7       653       632 VX     
5     7    10       453       445 B6     
6     7    10       432       433 VX     
7     9    20      1014      1007 AA

극단적 지연의 대부분이 7월에 발생 (7건 중 5건)
가장 빈번한 위반 항공사는 Virgin America (VX).
“7월에 일찍 도착하고 VX를 피하라.”
이상값 자체가 오해의 소지가 있을 수 있음
- 2013년 뉴욕에서 샌프란시스코로 가는 항공편 중 극히 일부
- 정시에 회의에 도착하는 데 실패하는 비율을 2%로 낮춘다는 목표에 비하면 아주 작은 부분
- 2013년 7월 샌프란시스코 공항에서는 아시아나항공 214편의 착륙 사고라는 더욱 드문 사건이 발생

Note

이상값을 처리하는 방법은 그 원인에 따라 다름
데이터 불규칙성이나 오류로 인한 이상값은 수정
다른 이상값은 문제에 대한 중요한 직관을 제공할 수 있음
명확한 근거가 없는 이상 이상값은 절대 삭제해서는 안 됨
이상값이 삭제된 경우에는 이를 명확하게 보고해야 함

이상치 제거한 히스토그램

SF %>% 
  filter(arr_delay < 420) %>%
  ggplot(aes(arr_delay)) + 
  geom_histogram(binwidth = 15) + 
  labs(x = "Arrival delay (in minutes)")

대부분의 항공편이 지연 없이 도착하거나 60분 미만의 지연이 발생
장시간 지연이 발생할 가능성이 높은 시기를 예측할 수 있는 패턴을 파악할 수 있을까?
이상값은 해당 월 또는 항공사가 장시간 지연과 관련이 있을 수 있음을 시사

SF %>% 
  mutate(long_delay = arr_delay > 60) %>%
  group_by(month, long_delay) %>%
  count() %>%
  pivot_wider(names_from = month, values_from = n) %>%
  data.frame()

  long_delay  X1  X2  X3  X4   X5  X6  X7   X8   X9  X10  X11  X12
1      FALSE 856 741 812 993 1128 980 966 1159 1124 1177 1107 1093
2       TRUE  29  21  61 112   65 209 226   96   65   36   51   66

6월과 7월이 문제

SF %>% 
  mutate(long_delay = arr_delay > 60) %>%
  group_by(carrier, long_delay) %>%
  count() %>%
  pivot_wider(names_from = carrier, values_from = n) %>%
  data.frame()

  long_delay   AA  B6   DL   UA   VX
1      FALSE 1250 934 1757 6236 1959
2       TRUE  148  86   91  492  220

두 가지 간단한 분석으로 여행자에게 6월과 7월에 더 일찍 도착할 계획을 세우고 SFO로 비행할 항공사로 델타항공(DL)을 고려하도록 조언할 것을 암시

통계 모형

변동성 설명

지금까지의 노력은 항공편마다의 연착 변동성을 일부 설명하려 한 것
통계 모형은 변수를 서로 연관시키는 방법을 제공
관심 있는 현상의 이해 증진

출발예정 시각이 항공편 지연의 기대치에 어떤 영향을 미칠까?

항공편 지연이 하루 동안 연쇄적으로 발생하는 경향이 있기 때문에 이른 시간에 출발하는 항공편이 지연될 가능성이 적다고 생각됨
이 이론이 데이터로 뒷받침되는가?

시간대를 고려하는 것부터 시작하자.
hour: nycflights13 패키지에서 flights 데이터프레임에는 출발예정 시각을 지정하는 변수

SF %>%
  group_by(hour) %>%
  count() %>%
  pivot_wider(names_from = hour, values_from = n) %>%
  data.frame()

  X5  X6   X7  X8  X9  X10 X11 X12 X13 X14 X15 X16  X17  X18 X19 X20 X21
1 55 663 1696 987 429 1744 413 504 476 528 946 897 1491 1091 731 465  57

이른 오전에서 오전 중반, 늦은 오후부터 이른 저녁까지는 많은 항공편이 예정
오전 5시 이전이나 오후 10시 이후에는 예정된 항공편이 없음

연착과 시간대와의 관계 탐구

SF %>%
  ggplot(aes(x = hour, y = arr_delay)) +
  geom_boxplot(alpha = 0.1, aes(group = hour)) +   # 연착 시간의 분포
  geom_smooth(method = "lm") + 
  xlab("Scheduled hour of departure") + 
  ylab("Arrival delay (minutes)") + 
  coord_cartesian(ylim = c(-30, 120))

평균 연착 시간은 하루 동안 시간이 지날수록 증가
추세선 자체는 선형회귀 모형으로 생성

mod1 <- lm(arr_delay ~ hour, data = SF)
broom::tidy(mod1)

# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)   -22.9     1.23       -18.6 2.88e- 76
2 hour            2.01    0.0915      22.0 1.78e-105

이 모형은 실제 항공편의 도착이 시간대별로 2.0 ± 0.18 분씩 증가한다고 95%의 신뢰도로 이야기함
매우 작은 유의확률: 시간대와 항공편 지연 사이에 연관성이 없다면 이렇게 극단적인 결과가 나올 가능성이 매우 낮다는 것을 의미

항공편 지연을 설명하는 데 도움이 될 수 있는 추가 요인은 무엇이 있을까?
- 출발 공항, 항공사, 연도, 월, 요일 등
- 또한 6월과 7월이 지연이 긴 달이라는 사실에 대해 이미 알고 있음

SF <- SF %>% 
  mutate(
    day = as.Date(time_hour), 
    dow = as.character(wday(day, label = TRUE)),  # day of the week
    season = ifelse(month %in% 6:7, "summer", "other month")
  )
mod2 <- lm(arr_delay ~ hour + origin + carrier + season + dow, data = SF)

broom::tidy(mod2)

# A tibble: 14 × 5
   term         estimate std.error statistic   p.value
   <chr>           <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept)   -24.6      2.17     -11.3   1.27e- 29
 2 hour            2.08     0.0898    23.2   1.44e-116
 3 originJFK       4.12     1.00       4.10  4.17e-  5
 4 carrierB6     -10.3      1.88      -5.49  4.07e-  8
 5 carrierDL     -18.4      1.62     -11.4   5.88e- 30
 6 carrierUA      -4.76     1.48      -3.21  1.31e-  3
 7 carrierVX      -5.06     1.60      -3.17  1.54e-  3
 8 seasonsummer   25.3      1.03      24.5   5.20e-130
 9 dowMon          1.74     1.45       1.20  2.28e-  1
10 dowSat         -5.60     1.55      -3.62  2.98e-  4
11 dowSun          5.12     1.48       3.46  5.32e-  4
12 dowThu          3.16     1.45       2.18  2.90e-  2
13 dowTue         -1.65     1.45      -1.14  2.53e-  1
14 dowWed         -0.884    1.45      -0.610 5.42e-  1

적합된 계수(estimate)는 기준 그룹인 EWR 대신 JFK에서 출발하는 경우 평균 지연에 4.1분을 더해야 한다는 것을 알려줌.
델타항공은 다른 항공사보다 연착이 적음
6월과 7월(25분), 일요일(5분)에 평균적으로 지연 시간이 더 긺
- 아시아나 착륙 사고가 7월
토요일은 평균 연착 시간보다 6분 짧음 (각 요일은 알파벳순으로 먼저 오기 때문에 기준 그룹으로 선택한 금요일과 비교)
표준 오차: 추정량의 정밀도를 알려줌
유의확률: 변수 간에 체계적인 연관성이 없더라도 개별 패턴이 우연히 발생할 것으로 예상되는 것과 일치할 확률
다른 모형: “항공편이 100분 이상 연착할 확률은 얼마나 될까?”

교락

다른 요인에 의한 설명

상관관계가 인과관계를 의미하지는 않는다.

상관관계가 인과 관계를 암시하는 경우가 많음 (임의화검정에서뿐만은 아님)
관측 연구에서의 주요 관심사는 두 요인 간의 겉보기 관계를 실제로 결정할 수 있는 다른 요인에 의해 실제 연관성이 왜곡되고 있는 것은 아닌지의 여부
- 이들 요인은 조사 중인 상관(또는 인과)관계에 혼선을 줄 수 있음 (confounder, 교락)

임의화 검정

증거 기반 연구의 표준
온라인 업계에서는 A/B 테스트라고도 불림
일반적으로 처치의 효과를 비교하기 위해 수행(예: 두 가지 형태의 웹 페이지)
새로운 개입을 받는 사람과 대조군(또는 표준 처치)을 받는 사람을 통제함으로써 다른 모든 요인이 두 그룹 간에 평균적으로 균형을 이루도록 함
시험이 끝날 때 측정된 결과에 차이가 있는 경우, 이는 치료법의 적용에 기인한 것이라고 결론을 내릴 수 있음
주의: 피험자가 처치를 따르지 않거나 추적 관찰에서 사라지는 경우 임의화 검정도 교락을 일으킬 수 있음

관측 연구

실용성의 문제
- 담배가 폐암을 유발하는지 여부를 확인하기 위해 일부 어린이에게는 담배를 피우고 다른 어린이에게는 담배를 피우지 않도록 무작위로 배정
- 커피가 장기적으로 건강에 영향을 미치는지 확인하기 위해 성인을 무작위로 커피를 마시게 하거나 마시지 않게 함
관측 연구만이 중요한 질문에 답할 수 있는 유일한 방법일 수 있음
데이터과학자는 주로 관측자료로 작업

미국 수학능력시험(SAT) 자료

미국 50개 주의 평균 교사 급여(2010년 기준)와 평균 총 SAT 점수에 대한 관측 자료

교사 급여가 높을수록 주 차원에서 시험 결과가 더 좋은 것과 관련이 있을까? 그렇다면 시험 성적을 개선하기 위해 급여를 조정해야 하는가?

SAT_2010 <- SAT_2010 %>%
  mutate(Salary = salary/1000)
SAT_plot <- ggplot(data = SAT_2010, aes(x = Salary, y = total)) + 
  geom_point() + 
  geom_smooth(method = "lm") + 
  ylab("Average total score on the SAT") + 
  xlab("Average teacher salary (thousands of USD)")
SAT_plot

SAT_mod1 <- lm(total ~ Salary, data = SAT_2010)
broom::tidy(SAT_mod1)

# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  1871.      113.       16.5  1.86e-21
2 Salary         -5.02      2.05     -2.45 1.79e- 2

교락 요인

숨어 있는 또 다른 중요한 요소: 각 주에서 SAT를 치르는 학생의 비율은 매우 다양(2010년에는 3%에서 93%까지)

SAT_2010 %>% skim(sat_pct)

Variable type: numeric

var	n	na	mean	sd	p0	p25	p50	p75	p100
sat_pct	50	0	38.52	31.99	3	6	27	68	93

주(state)를 두 그룹으로 나누는 SAT_grp 층화 변수를 만들자.

SAT_2010 <- SAT_2010 %>%
  mutate(SAT_grp = ifelse(sat_pct <= 27, "Low", "High"))
SAT_2010 %>%
  group_by(SAT_grp) %>% count()

# A tibble: 2 × 2
# Groups:   SAT_grp [2]
  SAT_grp     n
  <chr>   <int>
1 High       25
2 Low        25

SAT_plot %+% SAT_2010 + 
  aes(color = SAT_grp) + 
  scale_color_brewer("% taking\nthe SAT", palette = "Set2")

7장의 반복법을 사용하여 두 개의 개별 그룹에 적합한 선형회귀 모형의 계수 도출

SAT_2010 %>%
  group_by(SAT_grp) %>%
  group_modify(~broom::tidy(lm(total ~ Salary, data = .x)))

# A tibble: 4 × 6
# Groups:   SAT_grp [2]
  SAT_grp term        estimate std.error statistic  p.value
  <chr>   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 High    (Intercept)  1428.       62.4     22.9   2.51e-17
2 High    Salary          1.16      1.06     1.09  2.85e- 1
3 Low     (Intercept)  1583.      141.      11.2   8.52e-11
4 Low     Salary          2.22      2.75     0.809 4.27e- 1

심슨 역설

심슨 역설의 정량적 버전
- SAT 응시 비율이 낮은 주에서는 교사 급여와 SAT 점수가 양의 상관관계를 보임
- SAT 응시 비율이 높은 주에서는 교사 급여와 SAT 점수가 양의 상관관계를 보임
- 응시 비율을 무시하면 교사 급여와 SAT 점수가 음의 상관관계를 보임.
교락변수가 측정되는 경우 문제를 해결하는 것은 간단
- 층화
- 다중 회귀

회귀 모형에 sat_pct 변수를 예측 변수로 추가

SAT_mod2 <- lm(total ~ Salary + sat_pct, data = SAT_2010)
broom::tidy(SAT_mod2)

# A tibble: 3 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  1589.      58.5       27.2  2.16e-30
2 Salary          2.64     1.15       2.30 2.62e- 2
3 sat_pct        -3.55     0.278    -12.8  7.11e-17

sat_pct를 통제했을 때 Salary의 기울기가 양수이고 통계적으로 유의
- SAT_grp로 층화분석한 결과와 일치
교사 급여가 SAT 점수 상승의 원인이라고 단정할 수는 없음
- 그러나 교락 요인을 고려한 후 관찰한 연관성은 그렇지 않은 연관성보다 더 신뢰할 수 있음
관측자료에서 의미를 찾으려 할 때는 관측된 연관성을 왜곡할 수 있는 잠재적인 교락 요인을 찾아내는 것이 중요

유의확률의 위험성

통계적 가설검정의 위험성

유의확률(p값)은 데이터의 패턴이 실제로 무작위 우연의 결과라는 귀무가설 하에서 통계량이 관찰된 통계량보다 더 극단적인 값으로 나타날 확률
항공사 데이터 예제
- “예측 변수와 항공편 지연 사이에 연관성이 없다”
SAT 성적과 연봉 예제
- “실제(모집단) 회귀 계수(기울기)가 0”
통상적으로 유의확률이 유의수준(예: \(\alpha\)=0.05) 이하이면 귀무가설을 기각

“전부 아니면 전무” 사고방식의 문제
- 귀무가설이 참으로 설정된 세계를 시뮬레이션하여 유의확률 계산
- 유의확률이 크면 데이터가 시뮬레이션과 일치함을 나타낸다.
- 유의확률이 매우 작으면 시뮬레이션이 관찰된 패턴의 메커니즘을 설명하는 것과 관련이 없음을 의미
  - 이 수치만으로는 어떤 종류의 가설이 패턴과 관련성이 있는지 알기 어려움
- “유의한 결과”는 귀무가설을 거부할 수 있다는 뜻 — 어떤 가설을 받아들여야 하는지는 알려주지 않음

미국 통계학회의 유의확률에 대한 성명서(Wasserstein and Lazar 2016)

유의확률은 데이터가 지정된 통계 모형과 얼마나 호환되지 않는지를 나타낼 수 있다.
유의확률은 연구 가설이 참일 확률이나 데이터가 무작위적인 우연에 의해 생성되었을 확률을 측정하지 않는다.
과학적 결론 또는 정책 결정은 유의확률이 특정 임계값을 통과했는지 여부에만 근거해서는 안 된다.
적절한 추론을 위해서는 완전한 보고와 투명성이 필요하다.
유의확률 또는 통계적 유의성은 효과의 크기나 결과의 중요성을 측정하는 것이 아니다.
유의확률은 그 자체로는 모형이나 가설에 대한 증거의 좋은 척도를 제공하지 않는다.

Note

항상 결정(귀무가설 기각 또는 기각 실패)만 보고하지 말고 실제 유의확률(또는 p <0.0001과 같은 문구)을 보고하라.
신뢰 구간은 종종 더 해석하기 쉬우므로 함께 보고하라.

다중 검정

분석에는 하나의 가설 검정이 아니라 수십 개 이상의 가설이 포함될 수 있음
이러한 상황에서는 유의수준이 낮아도 데이터와 귀무가설 사이의 불일치를 입증하지 못함
임의화검정에서는 주결과 및 보조 결과를 미리 지정 가능
- 본페로니 보정 등
관찰연구에서는 사전 지정된 프로토콜이 없으므로 어떤 수준의 본페로니 보정이 적절한지 알기 어려움

갈림길의 정원 (Andrew Gelman)

http://www.stat.columbia.edu/~gelman/research/unpublished/p_hacking.pdf

대부분의 분석에는 최종 분석이 설정되기 전에 데이터 코딩, 중요한 요인 결정, 모형 상정 및 수정 등에 대한 많은 결정이 수반
데이터를 살펴보고 간결한 표현을 구성하는 작업이 포함
- 연속형 예측 변수를 임의의 그룹으로 잘라 해당 예측 변수와 결과 간의 관계를 평가
- 탐색적 자료분석 과정에서 특정 변수를 회귀 모형에 포함하거나 제외
다른 옵션보다 더 많은 신호(또는 더 작은 유의확률)를 산출하는 결정을 선택 — 귀무가설 기각 편향
임상시험: 분석 계획을 미리 지정하고 게시해야 하므로 갈림길의 정원 문제가 덜 일반적
대부분의 데이터과학 문제: 재현 가능한 결과에 대한 의문으로 이어질 수 있음