Modern Data Science with R_Cpt7

7. Statistical foundations

7.1 Samples and populations

Example: Sampling from the population
당신은 샌프란시스코(공항코드 SFO)에서 비즈니스 미팅이 있는 여행객을 시간 t에 늦지않게 만나도록 비행계획을 짜야 한다.R의 nycflights13 패키지를 사용한다.

library(tidyverse)
library(mosaic)
library(nycflights13)

SF <- nycflights13::flights %>% 
  filter(dest == "SFO", !is.na(arr_delay))

set.seed(101)
Sample25 <- SF %>%
  sample_n(size = 25)

favstats(~arr_delay, data = Sample25)

단순하게는 가장 많이 지연된 경우인 max의 103분을 참고하여 계획 작성하기를 고려할 수 있다. 하지만 이는 25개 샘플의 결과이기 때문에, 다음의 전체 결과를 보고 비교할 수 있다.

favstats(~arr_delay, data = SF)

전체 데이터(모집단, population)의 경우 최대 지연시간은 1007분인 것을 확인할 수 있다. 이 경우 미팅을 놓치지 않기 위해선 적어도 하루 전, 약 17시간 전으로 비행 계획을 작성할 필요가 있어보인다. 더욱 안전하기 위해선 다음 경우도 고려할 필요가 있다.

1. 추가 여행일은 여행자 입장에서 비용이 많이 든다.
   
2. 1,007분 보다 더 비행이 지연될 가능성이 없지 않다.
   
   

비용과 지각 가능성을 고려한 현명한 계획을 짜기 위해 2%의 지각 가능성을 두고 98%는 정시에 도착한다고 가정해보자.

qdata(~arr_delay, p = 0.98, data = Sample25)

  98% 
67.48 

98% 백분위 수에 해당하는 지연 시간이 약 67분으로 나타났다. 사실 98분위가 몇 분인지는 중요한 것이 아니다. 2%의 지각 가능성을 열어둔 결과가 과연 효과적인 계획이었는지가 중요하다. 모집단 데이터가 주어져있다면 이를 확인하는 것이 가능하다.

tally(~arr_delay < 68, data = SF, format = "proportion")

arr_delay < 68
     TRUE     FALSE 
0.9297047 0.0702953 

지연 확률이 2%를 넘은 7%로 나타났다. 처음 계획이 잘못 계산된 것이다. 모집단 데이터로 다시 계산해본다.

qdata(~arr_delay, p = 0.98, data = SF)

98% 
153 

모집단을 통해 계산한 결과, 98%의 도착률을 맞추기 위해서는 68분이 아니라 153분의 여유를 두어서 계획해야 함을 알 수 있다.

하지만 본 예시가 아니라 현실에서는 전체 데이터를 접하기 어렵고 샘플 데이터를 가지고 있을 확률이 높다. 이후로는 샘플 데이터를 가지고 높은 효과를 볼 수 있는 방법에 대해 알아볼 예정이다.

7.2 Sample statistics

표본 집단의 통계가 모집단을 얼마나 잘 반영하는지 알아내는 확립된 기술이 있다.

The sampling distribution
표본 추출을 여러 번 실행하여 각 시도에서 표본 통계를 계산하는 방법.

n <- 25
mean(~arr_delay, data = sample_n(SF, size = n, replace = FALSE))

[1] -5.92

아마 이보다도 더 많은 실행을 해보는 것이 좋을 것이다.

trials <- do(500) *
  mean(~arr_delay, data = sample_n(SF, size = n, replace = FALSE))
head(trials)

500번 표본 통계를 반복 추출하여 평균 도착 지연시간을 계산해보았다. 결과가 얼마나 분산되어 있는지 보자.

favstats(~mean, data = trials)

500번이나 반복하였으니 표본의 크기는 중요하지 않다. 평균이 중앙값의 약 두 배 가량인 것을 보아 분포가 약간 오른쪽으로 치우쳐져있음을 알 수 있다.

mean(~mean, data = trials) + 2*sd(~mean, data = trials) * c(-1, 1)

[1] -17.23359  23.41647

평균과 표준편차의 계산을 통해 95% 신뢰구간을 보는 것 역시 판단에 참고할 수 있다.

기존의 n=25였던 표본보다 더 큰 n=100인 표본을 사용한다면 어떤 결과를 보일까?

trials_100 <- do(500) *
  mean(~arr_delay, data = SF %>% sample_n(size = 100, replace = FALSE))

rbind(trials %>% mutate(n = 25),
      trials_100 %>% mutate(n = 100)) %>% 
  ggplot(aes(x = mean)) +
           geom_histogram(bins = 30) +
           facet_grid(~n) +
           xlab("Sample mean")

(좌) n=25인 표본 분포, (우) n=100인 표본 분포

위 분포도를 통해 알 수 있는 것:
1. 두 표본 분포는 같은 값에서 중앙값을 가진다.
2. 표본의 크기가 큰 분포가 더 작은 표준오차를 가진다. 즉, 더 신뢰도가 높다.
3. 표본의 크기가 큰 경우 종모양 분포를 갖게 되고, 이를 정규분포라고 한다.

7.3 The bootstrap

’bootstrap’이란 모집단에 접근할 수 없는 없는 경우에도 표본 분포를 통해 근사화할 수 있게 하는 통계학적 방법이다. 가지고 있는 표본에서 표본을 추출할 것이며, 이를 ’resampling’이라고 한다. 설명을 위해 flights 데이터를 사용하였다.

small <- sample_n(SF, size = 3, replace = FALSE)
small

위 과정을 중복을 허용하여 resampling한 결과는 다음과 같다.

small %>% sample_n(size = 3, replace = TRUE)

보다시피 부트스트래핑은 데이터를 수집하여 새로운 데이터를 분석하는 것이 아니다. 전에 하던 비행 계획 문제로 돌아가서, 부트스트래핑을 사용하여 크기가 200인 표본에서 계산된 평균 도착 시간의 신뢰도를 추정해보자.

n <- 200
orig_sample <- SF %>% sample_n(size = n, replace = FALSE)

mean(~arr_delay, data = sample_n(orig_sample, size = n, replace = TRUE))

[1] -0.88

위 과정을 반복 resampling하여 평균 도착 지연시간을 계산해보자.

bootstrap_trials <- do(500) * mean(~arr_delay, data = sample_n(orig_sample, size = n, replace = TRUE))

favstats(~mean, data = bootstrap_trials)

부트스트랩을 사용하여 계산한 표준 오차가 모집단에서 반복 샘플을 취하여 계산한 샘플링 분포의 표준 오차에 대한 합리적인 근사치라는 것을 확인할 수 있다.(모집단의 평균이 2.67분이었으므로 1.4분이라는 조금의 차이만 난다.)

회의 시간의 98%를 제 시간에 도착하기로 설정한 항공편 계획 문제로 돌아가보자. 문제에 부트스트래핑을 적용하여 n=200인 표본의 통계 신뢰도를 파악할 수 있다.

qdata(~arr_delay, p = 0.98, data = orig_sample)

   98% 
195.26 

표본은 195분 먼저 도착하도록 비행기를 예약할 것을 제안한다. 본 결과의 신뢰도를 부트스트래핑을 사용하여 확인해보자.

n <- nrow(orig_sample)
num_trials <- 500
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)

결과에 따르면 부트스트랩 표준 오류는 약 182분이고, 95% 신뢰 구간은 182 ± (66*2)분이다. 역시 앞에서 본 195분과 약 13분 정도의 차이가 나는 것을 볼 수 있다.

문제를 해결하는 한 가지 방법은 98번째 백분위수에 대한 보다 신뢰할 수 있는 추정치를 얻기 위해 더 많은 데이터를 수집하는 것이다. n=10,000인 케이스의 표본을 생성하는 작업을 수행해보자.

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)

신뢰구간이 153 ± 9분으로 훨씬 작아졌다. 더 많은 데이터가 있으면 특히 꼬리 부분에서 더 쉽게 추정치를 세분화할 수 있다.

7.4 Outliers

데이터는 비정상적이거나 극단적인 이벤트를 식별할 때에도 도움이 된다. 7시간 이상 지연된 비행을 극단적인 상황으로 간주한다고 가정할 때, 420분의 수치는 심각하게 지연된 항공편의 지표로 유용할 수 있다.

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

비행 지연의 대부분이 7월, VX항공에서 발생했다. 이 사실은 항공 계획을 짤 때 7월에는 더 시간을 앞당겨서 계획하고, 되도록 VX항공편을 피하라고 조언할 수 있다.

하지만 이상치 자체에 오해가 있을 수도 있다. 해당 특이점은 2013년 뉴욕에서 샌프란시스코로 향하는 항공편의 2%라는 극히 일부를 차지한다. 분포의 주요 부분을 더 잘 이해하기 위해 이러한 점을 제거할 수 있다.
대부분의 항공편은 60분 미만의 지연 없이 도착한다. 더 긴 지연이 발생할 가능성이 있을 때 전조가 될 수 있는 패턴을 식별할 수 있을까? 특이값은 month 또는 carrier가 긴 지연과 연결될 수 있음을 시사했다. 대다수의 데이터에서 이것이 어떻게 작동하는지 보자.

SF %>% 
  mutate(long_delay = arr_delay > 60) %>% 
  tally(~long_delay | month, data = .)

          month
long_delay    1    2    3    4    5    6    7    8    9   10   11   12
     TRUE    29   21   61  112   65  209  226   96   65   36   51   66
     FALSE  856  741  812  993 1128  980  966 1159 1124 1177 1107 1093

위 결과에서 60분 이상 지연되는 경우가 6~7월에 가장 많기 때문에, 6월과 7월이 문제의 달임을 알 수 있다.

SF %>% filter(arr_delay < 420) %>% 
  ggplot(aes(arr_delay)) +
  geom_histogram(binwidth = 15)

2013년 NYC공항에서 샌프란시스코행 항공편이 7시간 미만 지연된 항공편의 도착 지연 분포다. 이 분포는 이상값을 제거한 후에도 긴 오른쪽 꼬리를 특징으로 보인다.

SF %>% 
  mutate(long_delay = arr_delay > 60) %>% 
  tally(~long_delay | carrier, data = .)

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

델타 항공(DL)은 합리적인 성능을 제공한다. 이 두 가지 간단한 분석을 통해 여행자에게 6월과 7월에 더 일찍 도착할 계획을 세우고 델타항공을 SFO 여행 항공사로 고려하도록 조언할 수 있다.

7.5 Statistical models: Explaining variation

지금까지 비행 지연의 가능성을 낮추기 위해 month와 carrier를 이용했다. 이처럼 변수를 서로 관련시켜 시스템을 이해하도록 하는 것을 통계 모델링이라고 한다. 다른 예를 들면, 일반적으로 사람들은 비행 지연이 하루 종일 연속되는 경향이 있기 때문에 이른 시간대의 비행이 지연될 가능성이 적다고 생각한다. 모델링을 설명하기 위해 출발 시간이 항공편 지연에 어떤 영향을 미치는지 알아보기로 했다.

tally(~hour, data = SF)

hour
   5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20   21 
  55  663 1696  987  429 1744  413  504  476  528  946  897 1491 1091  731  465   57 

결과를 통해 항공편이 이른 아침부터 정오까지, 그리고 늦은 오후부터 이른 저녁까지 예정되어 있음을 알 수 있다. 또한 새벽 5시 전이나 밤 9시 이후에는 항공편이 없다는 것을 알 수 있다.

도착 지연이 시간에 따라 어떻게 달라지는지 살펴보기 위해 두 가지 방법을 사용한다.

1. 도착 지연의 분포를 보여주기 위해 표준 box-and-whisker 그래프를 사용한다.
   
2. 하루 동안의 평균 도착 지연을 추적할 수 있는 선형 모델이라는 일종의 통계 모델을 사용한다.
   
   

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)
msummary(mod1)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -22.93267    1.23275  -18.60   <2e-16 ***
hour          2.01487    0.09154   22.01   <2e-16 ***

Residual standard error: 46.82 on 13171 degrees of freedom
Multiple R-squared:  0.03548,   Adjusted R-squared:  0.03541 
F-statistic: 484.5 on 1 and 13171 DF,  p-value: < 2.2e-16

Estimate 아래의 숫자는 도착 지연이 시간 당 2분 더 높을 것으로 예상됨을 나타낸다. 가장 오른쪽의 p값은 매우 작은 2e-16이며 이는 비행 출발시간과 비행 지연 사이에 연관성이 없다면 극단적인 결과를 볼 가능성이 매우 낮다는 것을 의미한다.

비행 지연을 설명하는 데 더 도움이 될 수 있는 추가 요소는 무엇일까? 데이터 랭글링을 통해 연도, 월, 일에서 요일(dow)을 추출할 수 있다. 또한 6, 7월처럼 지연이 긴 달을 season 변수로 만들 수 있다. 이들은 응답 변수인 도착 지연을 설명하기 위한 설명 변수로 사용될 것이다.

library(lubridate)
SF <- SF %>% 
  mutate(day = ymd(paste0(year, "-", month, "-", day)),
         dow = as.character(wday(day, label = TRUE)),
         season = ifelse(month %in% 6:7, "summer", "other month"))

mod2 <- lm(arr_delay ~ hour + origin + carrier + season + dow, data = SF)
msummary(mod2)

              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -24.54081    2.17451 -11.286  < 2e-16 ***
hour           2.06423    0.08981  22.985  < 2e-16 ***
originJFK      4.19889    1.00442   4.180 2.93e-05 ***
carrierB6    -10.33223    1.87970  -5.497 3.94e-08 ***
carrierDL    -18.40107    1.61464 -11.396  < 2e-16 ***
carrierUA     -4.78247    1.48083  -3.230 0.001243 ** 
carrierVX     -5.03651    1.59792  -3.152 0.001626 ** 
seasonsummer  25.32723    1.03068  24.573  < 2e-16 ***
dow목          3.33587    1.44610   2.307 0.021082 *  
dow수         -0.50140    1.44913  -0.346 0.729347    
dow월          1.44377    1.44444   1.000 0.317553    
dow일          5.53721    1.47092   3.764 0.000168 ***
dow토         -5.94596    1.56172  -3.807 0.000141 ***
dow화         -1.84872    1.45023  -1.275 0.202412    

Residual standard error: 45.44 on 13159 degrees of freedom
Multiple R-squared:  0.09218,   Adjusted R-squared:  0.09129 
F-statistic: 102.8 on 13 and 13159 DF,  p-value: < 2.2e-16

Estimate 열의 숫자는 JFK공항에서 출발할 경우 평균 지연 시간에 4.1분이 더 소요될 것을 의미한다. 델타항공(DL)은 다른 항공사에 비해 평균 지연이 더 좋다(-18). 비행 지연 시간은 평균적으로 ’summer’인 6월과 7월에 길고, 요일 중에선 일요일에 가장 길다. 이 예제에서는 lm()을 사용하여 선형 모델을 구성했다.

7.6 Confounding and accounting for other factors

통계를 공부하며 상관관계가 인과관계를 암시할 때가 많지만, 두 경우가 서로 다르다는 것을 이해해야 한다. 예를 들어, 평균 교사 급여(2010년)와 미국 50개 주 각각의 평균 총 SAT 점수에 대한 관찰 데이터를 사용하여 상관관계와 인과관계가 혼란스러운 경우의 예시를 보자. SAT는 대학 입학에 사용되는 고난도 시험이다. 더 높은 교사 급여가 더 나은 시험 결과와 관련이 있을까? 그렇다면 시험 성적 향상을 위해 급여를 조정해야 하는가? 증명해보았다.

library(mdsr)
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

Figure 9.4: 2010년 미국 50개 지역의 평균 SAT 점수 대 평균 교사 급여(단위: 천 달러)의 산점도.

install.packages(“remotes”)
remotes::install_github(“mdsr-book/mdsr”)
mdsr패키지가 현재 R버전에서 단종되었기에 위 과정을 통해 패키지를 다운로드받아야 한다.

SAT_mod1 <- lm(total ~ Salary, data = SAT_2010)
msummary(SAT_mod1)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1871.104    113.141  16.538   <2e-16 ***
Salary        -5.019      2.048  -2.451   0.0179 *  

Residual standard error: 111.2 on 48 degrees of freedom
Multiple R-squared:  0.1113,    Adjusted R-squared:  0.09273 
F-statistic: 6.008 on 1 and 48 DF,  p-value: 0.01793

위와 같은 분석이 가능하지만, 각 주에서 SAT를 보는 학생들의 비율이 매우 차이난다는 점을 간과할 수 없다. 따라서 SAT_grp변수를 만들어 분석에 포함시킨다.

SAT_2010 <- SAT_2010 %>% 
  mutate(SAT_grp = ifelse(sat_pct <= 27, "Low", "High"))
tally(~SAT_grp, data = SAT_2010)

SAT_grp
High  Low 
  25   25 

SAT_plot %+% SAT_2010 + aes(color = SAT_grp)

2010년 미국 50개 주에서 평균 SAT 점수 대 평균 교사 급여의 산포도(단위: 천달러)
각 주에서 SAT를 치른 학생의 백분율로 계층화됨

coef(lm(total ~ Salary, data = filter(SAT_2010, SAT_grp == "Low")))

(Intercept)      Salary 
1583.265572    2.223863 

coef(lm(total ~ Salary, data = filter(SAT_2010, SAT_grp == "High")))

(Intercept)      Salary 
1428.375819    1.159223 

각 그룹에서 평균 교사 연봉은 SAT점수와 양의 관계에 있다. 정리하면 다음과 같다.

SAT 응시율이 낮은 주에서는 교사 급여와 SAT 점수가 긍정적으로 연관되어 있다.
SAT 응시율이 높은 주에서는 교사 급여와 SAT 점수가 긍정적으로 연관되어 있다.
모든 주에서 급여와 SAT 점수는 음의 상관 관계가 있다.

앞서서 SAT_grp를 고려하지 않았을 때는 음의 상관관계를 보였기 때문이다. 이처럼 혼란변수를 측정하면 문제를 해결하는 것이 간단해진다. 예측 변수로 sat_pct 변수를 회귀 모델에 추가해보자.

SAT_mod2 <- lm(total ~ Salary + sat_pct, data = SAT_2010)
msummary(SAT_mod2)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1589.0065    58.4707  27.176   <2e-16 ***
Salary         2.6370     1.1488   2.295   0.0262 *  
sat_pct       -3.5526     0.2785 -12.756   <2e-16 ***

Residual standard error: 53.18 on 47 degrees of freedom
Multiple R-squared:  0.8008,    Adjusted R-squared:  0.7924 
F-statistic: 94.49 on 2 and 47 DF,  p-value: < 2.2e-16

이제 Salary의 기울기가 양수이고 sat_pct를 제어할 때 통계적으로 유의함을 알 수 있다. 여전히 교사 급여가 SAT점수 향상에 영향을 준다고 결론을 내릴 수 없지만, 혼란 요인을 고려한 후 관찰한 연관성은 고려하기 전의 결과보다 더 신뢰할 수 있다.