1 모형과 가설검정 1 2
모형으로 다양한 통계검정을 이해하게 되면 비모수 통계검정을 비롯하여 표본집단에 대한 구분부터 각각을 개별적으로 이해하고 암기할 필요가 없으며 단 하나의 수식 \(y = \beta_0 + \beta_1 \times x\) 이라는 단순한 공식으로 상당수 통계검정을 처리할 수 있다. 예를 들어 상관분석, 일원 분산분석, t-검정, \(\chi^2\)-검정 등이 여기에 포함된다.
평균에 대한 가설검정은 다음 작업흐름도에 따라 어떻게 보면 기계적인 과정을 거쳐 판정에 이를 수 있다. 특히, Common statistical tests are linear models (or: how to teach stats) 블로그에 관련 내용이 잘 정리되어 있다.
1.1 검정 작업흐름도
webshot2::webshot("http://www.rebeccabarter.com/blog/2018-12-04_hypothesis_testing/",
selector = "#content > p:nth-child(3) > img",
"fig/hypothesis-testing.png")
1.2 선형모형과 통계검정
1.3 수식으로 일별하기
비모수 통계검정도 유사하게 확장이 가능하지만 우선 평균에 대한 통계검정을 수식으로 표현하면 다음과 같다.
- 일표본
t-검정: \(y = \beta_0 \qquad \mathcal{H}_0: y = 0\) - 이표본
t-검정: \(y = \beta_0 + \beta_1 x_1 \qquad \mathcal{H}_0: \beta_1 = 0\) - 다표본 분산분석: \(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +... \qquad \mathcal{H}_0: y = \beta_0\)
2 NBA 선수선발 데이터
실제 NBA 선수선발 데이터를 가지고 앞서 언급한 사례를 살펴보자. NBA 선수선발 데이터는 LeBron James 부터 선수선발 순위를 포함한 선수 역량과 관련된 정보가 정리되어 있다.
library(tidyverse)
download.file(url="http://khannay.com/data/NBA_Draft_Data.csv", destfile = "data/NBA_Draft_Data.csv")
nba_df <- read_csv("data/NBA_Draft_Data.csv") %>%
janitor::clean_names()
nba_df %>%
arrange(-pts) %>%
select(tm, player, pts, pick_number, everything()) %>%
# sample_n(100) %>%
DT::datatable()
3 피어슨 상관관계
피어슨 상관관계를 모형으로 표현하면 다음과 같다.
\(y = \beta_0 + \beta_1 \times x\) 선형모형에서 \(\beta_1 = 0\) 인지 검정하는 것과 유사하다.
# Fixed correlation
D_correlation <- data.frame(MASS::mvrnorm(30, mu = c(0.9, 0.9),
Sigma = matrix(c(1, 0.8, 1, 0.8), ncol = 2),
empirical = TRUE)) %>%
as_tibble() # Correlated data
# Add labels (for next plot)
D_correlation <- D_correlation %>%
mutate(label_num = sprintf('(%.1f,%.1f)', X1, X2),
label_rank = sprintf('(%i,%i)', rank(X1), rank(X2)))
# Plot it
fit <- lm(I(X2 * 0.5 + 0.4) ~ I(X1 * 0.5 + 0.2), D_correlation)
intercept_pearson <- coefficients(fit)[1]
P_pearson <- ggplot(D_correlation, aes(x=X1*0.5+0.2, y=X2*0.5+0.4)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, lwd=2, aes(colour='beta_1')) +
geom_segment(x=-100, xend=100,
y=intercept_pearson, yend=intercept_pearson,
lwd=2, aes(color="beta_0")) +
scale_color_manual(name=NULL, values=c("blue", "red"), labels=c(bquote(beta[0]*" (intercept)"), bquote(beta[1]*" (slope)")))
theme_axis(P_pearson, legend.position = c(0.4, 0.9))
앞선 이론적인 배경을 바탕으로 MP: Minutes Played와 PTS: Points를 산점도로 파악하고 이에 대한 상관관계를 통계검정과 선형모형으로 비교해보자. 즉, 인과관계는 모르겠지만 MP와 PTS간의 상관관계는 존재하는 파악하고자 한다.
\[ \hat {\beta} = {\rm cor}(Y_i, X_i) \cdot \frac{ {\rm SD}(Y_i) }{ {\rm SD}(X_i) } \]
\({\rm SD}(Y_i) = {\rm SD}(X_i)\) 동일한 경우에만 \(\hat{\beta}\)이 상관계수와 동일한 관계를 갖는 점을 명심하면 cor.test() 내장함수로 구현된 상관계수 통계검정과 glm() 선형모형을 비교해보면 약간의 차이는 있지만 t 통계량은 물론 p-값도 유사함을 발견할 수 있다.
Pearson's product-moment correlation
data: nba_df$mp and nba_df$pts
t = 59.585, df = 735, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8969551 0.9218351
sample estimates:
cor
0.9102128 # A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -3.43 0.225 -15.3 5.83e- 46
2 mp 0.569 0.00955 59.6 1.30e-283# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -1.12e-17 0.0153 -7.35e-16 1.00e+ 0
2 scale(mp) 9.10e- 1 0.0153 5.96e+ 1 1.30e-2834 일표본 t-검정
GLM 모형과 t-검정을 비교해보자. 즉, 수식은 다음과 같고 의미하는 바는 평균이 0과 다름이 있느냐를 따지게 된다. \(y = \beta_0 + \beta_1*x\)에서 \(x=0\)이라 절편만 검정하게 된다.
\[y = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\]
t-검정
# A tibble: 1 x 8
estimate statistic p.value parameter conf.low conf.high method alternative
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 9.19 50.5 1.90e-241 736 8.83 9.54 One Sam~ two.sided 선형 모형
# A tibble: 1 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 9.19 0.182 50.5 1.90e-2414.1 대응표본 t-검정
반복측정, 예를 들어 동일한 선수가 동일한 조건에서 약물복용 이후 약물복용 효과를 측정하고자 하는 경우 대응표본 t-검정(Paired t-test)이 되고 수식으로 다음과 같이 표현된다.
\(y_2-y_1 = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)
먼저 앞선 경우와 비교하기 위해서 시각화를 한다. sleep 데이터는 최면약(soporific drug)으로 통제집단과 비교하여 효과가 있는지를 10명의 환자에게 투여하여 비교해서 자체 내장된 데이터셋이다.
library(ggthemes)
sleep %>%
ggplot(aes(x = group, y = extra, color=ID)) +
geom_point(show.legend = FALSE) +
geom_line(aes(group = ID), show.legend = FALSE) +
theme_tufte()
대응표본 t-검정
# A tibble: 1 x 8
estimate statistic p.value parameter conf.low conf.high method alternative
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 -1.58 -4.06 0.00283 9 -2.46 -0.700 Paired t-~ two.sided 통계모형
# A tibble: 1 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -1.58 0.389 -4.06 0.002835 독립 두 표본 t-검정
독립된 두 표본(two-sample) 평균검정은 선형모형으로 다음과 같이 나타낼 수 있다. 여기서 \(x_i\)는 0 혹은 1 을 갖는 범주형 변수로 두 표본을 표식하는데 사용된다. 이를 확장하면 다집단에 대한 평균검정도 할 수 있는데 이 경우 ANOVA로 확장된다.
\(y_i = \beta_0 + \beta_1 x_i \qquad \mathcal{H}_0: \beta_1 = 0\)
이를 위해서 선수선발 16 이상은 하위, 1 - 15 순위는 상위권으로 두어 두개의 집단으로 나눠 평균득점에 차이가 있는 살펴보자.
nba_df <- nba_df %>%
mutate(pick_binary = if_else(pick_number > 15, "하위", "상위")) %>%
mutate(pick_binary = factor(pick_binary, levels = c("하위", "상위")))독립표본 t-검정
# A tibble: 1 x 9
estimate1 estimate2 statistic p.value parameter conf.low conf.high method
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 7.04 11.0 -12.0 2.00e-30 735 -4.66 -3.35 " Two~
# ... with 1 more variable: alternative <chr>범주를 갖는 회귀모형
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 7.04 0.244 28.8 7.70e-123
2 pick_binary상위 4.00 0.334 12.0 2.00e- 306 분산분석(ANOVA)
\(x_i\)는 (\(x=0\) or \(x=1\))을 갖는 지시변수로, 다른 모든 변수가 \(x_i=0\)이 될 때 \(x_i=1\)는 1을 갖는다. 수식으로는 다음과 같이 표현되고 그룹간에 평균적인 차이가 있는지 검정하는 것으로 집단이 2개이상될 때 사용한다.
\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +... \qquad \mathcal{H}_0: y = \beta_0\)
2표본 이상 다집단에 대한 평균에 대한 차이를 분석하기 위해서 다음과 같이 팀을 3개 뽑아내서 새로운 데이터프레임으로 준비를 한다. type-II sum of square 혹은 type-III sum of square 에 따라 차이가 있기는 하지만 glm 모형으로 비교를 할 수 있다.
ANOVA 분산분석
# A tibble: 2 x 5
term sumsq df statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 tm 77.3 2 1.59 0.210
2 Residuals 2142. 88 NA NA 가변수 회귀모형
glm(pts ~ tm, data=anova_df) %>%
broom::glance(.) %>%
mutate(sumsq = null.deviance - deviance) %>%
select(null.deviance, sumsq, deviance, everything())# A tibble: 1 x 8
null.deviance sumsq deviance df.null logLik AIC BIC df.residual
<dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int>
1 2219. 77.3 2142. 90 -273. 554. 564. 887 교호작용(interaction)
앞선 일원분산분석을 확장하여 이원분산분석도 가능한데 이런 경우 교호작용을 파악하는 것이 중요하다. 이를 위해서 두가지 모형에 대해서 비교를 해보자. 수식으로는 다음과 같고 두가지 요인 하나는 팀(tm), 다른 하나는 상위 혹은 하위 선수선발 가변수(pick_binary)로 두 변수간에 상승 혹은 역효과가 평균득점에 효과를 나타내는지 살펴보자.
\(y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1 X_2 \qquad \mathcal{H}_0: \beta_3 = 0\)
교호작용 검정
# A tibble: 4 x 5
term sumsq df statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 pick_binary 614. 1 34.2 0.0000000902
2 tm 20.4 2 0.568 0.569
3 pick_binary:tm 1.76 2 0.0490 0.952
4 Residuals 1526. 85 NA NA interaction 항을 갖는 회귀모형
anova_null = lm(pts ~ 1 + tm + pick_binary, anova_df)
anova_full = lm(pts ~ 1 + tm + pick_binary + tm:pick_binary, anova_df)
anova(anova_null, anova_full) %>%
broom::tidy(.)# A tibble: 2 x 6
res.df rss df sumsq statistic p.value
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 87 1528. NA NA NA NA
2 85 1526. 2 1.76 0.0490 0.9528 ANCOVA
평균득점같은 연속형 변수에 대한 검정을 할 때 다른 요인을 고려한 후에 특정 요인에 대한 효과를 살펴보고자 하는 경우가 상식적이다. 이런 경우 나이(age)를 고려한 후에 다양한 효과를 파악하는 것도 가능하다. \(x_i\)는 앞서와 같은 지시변수로 0 혹은 1을 갖는다는 점에서 동일하다.
\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_3 age\)
Ancova 통계검정
# A tibble: 4 x 5
term sumsq df statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 273. 1 17.3 0.0000735
2 pick_binary 512. 1 32.5 0.000000156
3 age 162. 1 10.2 0.00191
4 Residuals 1387. 88 NA NA Ancova GLM
ancova_full_lm <- lm(pts ~ 1 + pick_binary + age, data=anova_df)
ancova_reduced_lm <- lm(pts ~ 1 + pick_binary, data=anova_df)
anova(ancova_full_lm, ancova_reduced_lm, test="LRT") %>%
broom::tidy(.) # A tibble: 2 x 5
res.df rss df sumsq p.value
<dbl> <dbl> <dbl> <dbl> <dbl>
1 88 1387. NA NA NA
2 89 1549. -1 -162. 0.00137Kevin Hannay (Jan 13 2020), “Everything is a Regression:In search of unifying paradigms in statistics”↩
Rebecca Barter (December 4, 2018), “Which hypothesis test should I use? A flowchart: A flowchart to decide what hypothesis test to use.”↩
Phillip M. Alday (May 2018), “Explicit GLM(M) Equivalents for Standard Tests”↩