Yeong Chan Lee - 10_Regression and correlation methods

10-0. Load data

library(dplyr)

Warning: package 'dplyr' was built under R version 4.2.3


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(data.table)


Attaching package: 'data.table'

The following objects are masked from 'package:dplyr':

    between, first, last

Fundamentals of Biostatistics에 소개되어 있는 데이터를 활용해 회귀분석을 진행해봅시다. 먼저 estriol 데이터를 불러옵니다.

est <- fread("data/estriol.csv",data.table = FALSE)

library(ggplot2)
ggplot(data = est) +
  geom_point(aes(x=ESTRIOL, y=BIRTHWEIGHT)) + 
  theme_light(base_size=20)

주어진 데이터를 직접 눈으로 확인하여 결측치가 있는지, 이상치(outlier)가 있는지 메인 분석 전에 미리 확인하는 것이 중요합니다. 이를 Exploratory Data Analysis (EDA)라고도 합니다. 지금은 Estriol과 Birthweight간 그래프를 그려 확인해보았습니다. 선형성이 보이나요?

다음으로는 회귀분석을 해보겠습니다. Linear regression은 다음과 같은 모형을 가집니다.

\[ y=\alpha+\beta x+ e \]

회귀분석을 하는 기본 함수는 lm입니다. 그리고 (BIRTHWEIGHT ~ ESTRIOL) 이렇게 써져있는 것처럼 formula에 대해서 익숙해져야합니다. 수식을 쓴다고 생각하면 쉬울 수 있습니다. 내가 만들고자하는 모형에서 y에 해당하는 변수가 물결(~)의 왼쪽, x에 해당하는 변수가 오른쪽으로 나열되면 됩니다. 지금은 BIRTHWEIGHT가 y이기 때문에 왼쪽, ESTRIOL이 x이기 때문에 오른쪽에 위치하게 됩니다.

lm_fit <- lm(BIRTHWEIGHT ~ ESTRIOL, data = est)

b <- (cov(est$ESTRIOL, est$BIRTHWEIGHT)*30)/(var(est$ESTRIOL)*30)
a <- (mean(est$BIRTHWEIGHT) - (b * mean(est$ESTRIOL)))

model을 해석하기 위해서는 summary()함수를 씁니다. 다양한 정보가 있는데 우리는 Coefficients, R-squared, F test에 대한 p-value정도를 확인하면 됩니다.

summary(lm_fit)


Call:
lm(formula = BIRTHWEIGHT ~ ESTRIOL, data = est)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.1200 -2.0381 -0.0381  3.3537  6.8800 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  21.5234     2.6204   8.214 4.68e-09 ***
ESTRIOL       0.6082     0.1468   4.143 0.000271 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.821 on 29 degrees of freedom
Multiple R-squared:  0.3718,    Adjusted R-squared:  0.3501 
F-statistic: 17.16 on 1 and 29 DF,  p-value: 0.0002712

anova(lm_fit)

Analysis of Variance Table

Response: BIRTHWEIGHT
          Df Sum Sq Mean Sq F value    Pr(>F)    
ESTRIOL    1 250.57 250.574  17.162 0.0002712 ***
Residuals 29 423.43  14.601                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estriol 변수에 대한 추정량 \(\beta\)가 0.6082라고 계산이 되었고 p-value가 0.05보다 작기 때문에 통계적으로 유의한 변수라고 결론내릴 수 있습니다.

모델 자체가 의미있는지를 확인하기 위해서는 F test의 p-value를 확인합니다. 지금은 univariable linear regression이기 때문에 베타의 p-value랑 동일하고 역시 유의하다고 볼 수 있습니다. 이 결과를 사용하여 다음과 같은 회귀식을 만들 수 있습니다.

\[ Birthweight=21.5234+0.6082*Estriol \]

그럼 새로운 데이터가 들어왔을 때 예상되는 birthweight를 추정해볼 수 있겠습니다. Estriol이 15인 새로운 데이터가 생겼다고 가정합시다. 새로운 데이터를 모델에 넣어 y가 어떻게 나오는지 확인하고자할 땐 predict()함수를 사용합니다.

predict(lm_fit, data.frame(ESTRIOL=15))

       1 
30.64629

data <- fread("data/hyp_ihd_cohort.csv", data.table = FALSE)

t.test(filter(data, HYP == 1)$TOT_CHOL, filter(data, HYP == 0)$TOT_CHOL)


    Welch Two Sample t-test

data:  filter(data, HYP == 1)$TOT_CHOL and filter(data, HYP == 0)$TOT_CHOL
t = 3.7837, df = 495.11, p-value = 0.0001735
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  4.138246 13.078552
sample estimates:
mean of x mean of y 
 201.5006  192.8922

이번엔 중간고사 데이터를 다시 불러와보겠습니다. 그리고 범주형 변수인 HYP를 x변수로 두고 연속형 변수 TOT_CHOL을 y로 두어 분석해보겠습니다. 원래와 같이 two-sample t test를 하면 다음과 같이 나옵니다.

이전에 배운 two-sample t test와 연결지어서 어떤 공통점이 있는지 확인해보세요. 회귀분석은 등분산성을 가정하고 있기 때문에 two-sample t test에서는 등분산성을 가정한 결과와 동일하게 나옵니다.

lm_fit <- lm(TOT_CHOL ~ factor(HYP), data=data)
summary(lm_fit)


Call:
lm(formula = TOT_CHOL ~ factor(HYP), data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-83.892 -24.559  -2.113  19.983 155.608 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   192.892      1.135 169.921  < 2e-16 ***
factor(HYP)1    8.608      2.268   3.795 0.000155 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 33.5 on 1160 degrees of freedom
Multiple R-squared:  0.01226,   Adjusted R-squared:  0.01141 
F-statistic:  14.4 on 1 and 1160 DF,  p-value: 0.0001553

t.test(filter(data, HYP == 1)$TOT_CHOL, filter(data, HYP == 0)$TOT_CHOL, var.equal = TRUE)


    Two Sample t-test

data:  filter(data, HYP == 1)$TOT_CHOL and filter(data, HYP == 0)$TOT_CHOL
t = 3.7949, df = 1160, p-value = 0.0001553
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  4.15773 13.05907
sample estimates:
mean of x mean of y 
 201.5006  192.8922

Multivariable linear regression - regression과 ANOVA

이번엔 BP_SYS 변수를 활용해 고혈압 그룹을 새롭게 정의하여 세 군으로 만들어보겠습니다. (정상, 고혈압 1기, 고혈압 2기)

data$HYP2 <- ifelse(data$BP_SYS < 140, 0, 
                    ifelse(between(data$BP_SYS, 140, 160), 1, 
                           ifelse(data$BP_SYS > 160, 2, -1)))
table(data$HYP2)


   0    1    2 
1043  103   16

data %>% group_by(HYP2) %>% summarise(MEAN_SBP = mean(BP_SYS), 
                                      SD_SBP = sd(BP_SYS),
                                      MEAN_BMI = mean(BMI), 
                                      SD_BMI = sd(BMI))

# A tibble: 3 × 5
   HYP2 MEAN_SBP SD_SBP MEAN_BMI SD_BMI
  <dbl>    <dbl>  <dbl>    <dbl>  <dbl>
1     0     119.  10.9      23.5   3.18
2     1     145.   4.68     25.2   3.46
3     2     172.   9.12     26.1   2.61

회귀분석도 해보겠습니다. y는 BMI, x는 HYP2로 해봅시다.

lm_fit <- lm(BMI ~ factor(HYP2), data=data)
summary(lm_fit)


Call:
lm(formula = BMI ~ factor(HYP2), data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.8620 -2.2558 -0.1334  1.9380 13.8380 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   23.46202    0.09918 236.549  < 2e-16 ***
factor(HYP2)1  1.77140    0.33084   5.354 1.04e-07 ***
factor(HYP2)2  2.64423    0.80692   3.277  0.00108 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.203 on 1159 degrees of freedom
Multiple R-squared:  0.03188,   Adjusted R-squared:  0.03021 
F-statistic: 19.08 on 2 and 1159 DF,  p-value: 7.017e-09

anova(lm_fit) # 같은 f statistic과 pvalue

Analysis of Variance Table

Response: BMI
               Df  Sum Sq Mean Sq F value    Pr(>F)    
factor(HYP2)    2   391.6 195.798  19.082 7.017e-09 ***
Residuals    1159 11892.1  10.261                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

이번엔 ANOVA와 비교해봅시다. 세 집단의 평균의 차이를 볼 땐 ANOVA를 실시합니다.

귀무가설: 그룹간 평균의 차이가 없다.

대립가설: 적어도 한 그룹은 나머지 그룹의 평균과 차이가 있다.

aov_fit <- aov(BMI ~ factor(HYP2), data=data)
summary(aov_fit)

               Df Sum Sq Mean Sq F value   Pr(>F)    
factor(HYP2)    2    392  195.80   19.08 7.02e-09 ***
Residuals    1159  11892   10.26                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov_fit)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = BMI ~ factor(HYP2), data = data)

$`factor(HYP2)`
         diff        lwr      upr     p adj
1-0 1.7713976  0.9950081 2.547787 0.0000003
2-0 2.6442334  0.7506084 4.537858 0.0030995
2-1 0.8728358 -1.1471253 2.892797 0.5681241

이 때, ANOVA가 결과가 회귀분석에서의 F statistic과 동일합니다. 사실 둘의 분석 방법이 같기 때문입니다. Tukey multiple comparisons으로 어느 그룹끼리의 차이가 있는지도 확인할 수 있습니다.

Hypertension 예시

이번엔 새로운 데이터 hyp_ped를 불러와 multivariable linear regression을 수행할 것입니다. 가설은 “infant SBP가 Age와 Birthweight와 연관되어 있을 것이다.” 입니다. 직접 분석을 수행하고 결론을 내려보세요.

hyp_ped <- fread("data/hyp_ped.csv", data.table = FALSE)

lm_fit <- lm(SBP ~ AGE + BIRTHWEIGHT, data = hyp_ped)
summary(lm_fit)


Call:
lm(formula = SBP ~ AGE + BIRTHWEIGHT, data = hyp_ped)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0438 -1.3481 -0.2395  0.9688  6.6964 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 53.45019    4.53189  11.794 2.57e-08 ***
AGE          5.88772    0.68021   8.656 9.34e-07 ***
BIRTHWEIGHT  0.12558    0.03434   3.657   0.0029 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.479 on 13 degrees of freedom
Multiple R-squared:  0.8809,    Adjusted R-squared:  0.8626 
F-statistic: 48.08 on 2 and 13 DF,  p-value: 9.844e-07

결정계수와 상관계수

바로 위 모델에서 결정계수는 몇인가요? 수정된 결정계수를 보니 0.86정도로 1에 가까워 설명력이 높다고 생각할 수 있습니다. 물론 꼭 그런건 아니고 이것도 분야, 모델 등마다 다르긴 합니다.

estriol 데이터를 활용해서 상관계수에 대해 알아봅시다.

lm_fit <- lm(BIRTHWEIGHT ~ ESTRIOL, data = est)
summary(lm_fit)


Call:
lm(formula = BIRTHWEIGHT ~ ESTRIOL, data = est)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.1200 -2.0381 -0.0381  3.3537  6.8800 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  21.5234     2.6204   8.214 4.68e-09 ***
ESTRIOL       0.6082     0.1468   4.143 0.000271 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.821 on 29 degrees of freedom
Multiple R-squared:  0.3718,    Adjusted R-squared:  0.3501 
F-statistic: 17.16 on 1 and 29 DF,  p-value: 0.0002712

cor(est$ESTRIOL, est$BIRTHWEIGHT) # Pearson's correlation

[1] 0.6097313

corrplot()라는 함수를 활용해 다음과 같은 그래프도 그릴 수 있습니다.

#install.packages("corrplot")
library(corrplot)

Warning: package 'corrplot' was built under R version 4.2.3

corrplot 0.92 loaded

est2 <- select(est, -PERSON_ID)
corrplot(cor(est2), method="number", type="upper")

corrplot(cor(hyp_ped), method="number", type="upper")

univariable analysis에서 계산된 결정계수는 상관계수의 제곱과 동일합니다.

lm_fit_sum <- summary(lm_fit)

lm_fit_sum$r.squared == cor(est$ESTRIOL, est$BIRTHWEIGHT)^2

[1] TRUE

변수선택법

이번엔 중간고사 데이터(data)를 활용하여 모델을 만든 후 변수선택법을 통해 최적의 모델을 선정해볼 것입니다. 변수선택법 중에 잘 쓰이는 stepwise selection을 해볼 것입니다. R로 쉽게 적용할 수 있고 함수는 step입니다. 이 때 예제로 쓰이는 모델의 y가 무엇이고 x가 무엇인가요?

lm_fit <- lm(BMI ~ factor(SEX) + AGE + BP_SYS + TOT_CHOL + factor(HYP2), data=data)
summary(lm_fit)


Call:
lm(formula = BMI ~ factor(SEX) + AGE + BP_SYS + TOT_CHOL + factor(HYP2), 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2829 -1.9999 -0.1485  1.7210 14.0853 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   11.090102   1.120582   9.897  < 2e-16 ***
factor(SEX)2  -0.855703   0.181585  -4.712 2.75e-06 ***
AGE            0.002792   0.005978   0.467    0.641    
BP_SYS         0.079033   0.008867   8.913  < 2e-16 ***
TOT_CHOL       0.016804   0.002639   6.368 2.75e-10 ***
factor(HYP2)1 -0.490342   0.375019  -1.308    0.191    
factor(HYP2)2 -1.379186   0.885499  -1.558    0.120    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.971 on 1155 degrees of freedom
Multiple R-squared:  0.1701,    Adjusted R-squared:  0.1657 
F-statistic: 39.44 on 6 and 1155 DF,  p-value: < 2.2e-16

지금 모델에서 stepwise selection을 적용해 더 나은 모델을 찾아볼 것입니다.

final_model <- step(lm_fit) #forward로 하고 싶으면 direction="forward"을 넣습니다.

Start:  AIC=2537.55
BMI ~ factor(SEX) + AGE + BP_SYS + TOT_CHOL + factor(HYP2)

               Df Sum of Sq   RSS    AIC
- AGE           1      1.93 10197 2535.8
- factor(HYP2)  2     27.74 10222 2536.7
<none>                      10195 2537.6
- factor(SEX)   1    196.01 10391 2557.7
- TOT_CHOL      1    357.97 10553 2575.7
- BP_SYS        1    701.23 10896 2612.8

Step:  AIC=2535.77
BMI ~ factor(SEX) + BP_SYS + TOT_CHOL + factor(HYP2)

               Df Sum of Sq   RSS    AIC
- factor(HYP2)  2     27.38 10224 2534.9
<none>                      10197 2535.8
- factor(SEX)   1    194.13 10391 2555.7
- TOT_CHOL      1    368.91 10566 2575.1
- BP_SYS        1    749.08 10946 2616.1

Step:  AIC=2534.89
BMI ~ factor(SEX) + BP_SYS + TOT_CHOL

              Df Sum of Sq   RSS    AIC
<none>                     10224 2534.9
- factor(SEX)  1    229.57 10454 2558.7
- TOT_CHOL     1    381.29 10605 2575.4
- BP_SYS       1   1051.20 11275 2646.6

final model로 sex와 total cholesterol, SBP 변수만 남고 나머지는 제거된 것을 확인할 수 있습니다.

summary(final_model)


Call:
lm(formula = BMI ~ factor(SEX) + BP_SYS + TOT_CHOL, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.0439 -1.9979 -0.1247  1.7071 13.8900 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.236726   0.893062  13.702  < 2e-16 ***
factor(SEX)2 -0.904062   0.177296  -5.099 3.98e-07 ***
BP_SYS        0.069729   0.006390  10.912  < 2e-16 ***
TOT_CHOL      0.017177   0.002614   6.572 7.51e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.971 on 1158 degrees of freedom
Multiple R-squared:  0.1677,    Adjusted R-squared:  0.1655 
F-statistic: 77.76 on 3 and 1158 DF,  p-value: < 2.2e-16

Linear regression의 통계적 가정

선형회귀분석을 수행하면서 체크해야할 통계적 가정은 크게 독립성, 선형성, 정규성, 등분산성이 있습니다. 관련된 여러 방법들을 살펴봅시다.

#모델
lm_fit <- lm(BMI ~ factor(SEX) + AGE + BP_SYS + TOT_CHOL + factor(HYP2), data=data)

독립성 (independence)

Variation inflation factor (VIF)

다중공선성을 확인하는 VIF가 10이 넘는지 확인합니다.

#install.packages("car")
library(car)

Warning: package 'car' was built under R version 4.2.3

Loading required package: carData

Warning: package 'carData' was built under R version 4.2.3


Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

vif(lm_fit)

                 GVIF Df GVIF^(1/(2*Df))
factor(SEX)  1.085004  1        1.041635
AGE          1.135375  1        1.065540
BP_SYS       2.023403  1        1.422464
TOT_CHOL     1.039765  1        1.019689
factor(HYP2) 1.871674  2        1.169654

선형성 (linearity)

plot() 함수를 사용하면 다양한 통계적 가정에 관련된 그래프를 볼 수 있습니다. 잔차에 대한 선형성을 확인합니다.

plot(lm_fit, 1)

정규성 (normality)

Q-Q plot을 그려 직선의 형태를 띄고 있는지 확인합니다.

plot(lm_fit, 2)

또는 잔차에 대한 정규성 검정을 실시합니다.

shapiro.test(lm_fit$residuals)


    Shapiro-Wilk normality test

data:  lm_fit$residuals
W = 0.97711, p-value = 1.352e-12

등분산성 (homoscedasticity)

마찬가지로 등분산성도 residual 값이 균등하게 퍼져있는질 눈으로 체크하여 판단할 수 있습니다.

plot(lm_fit, 3)

또는 검정을 통해서도 확인할 수 있습니다. 여기서 귀무가설이 분산이 동질하다 (등분산성이다)는 것을 유의하세요.

library(car)
ncvTest(lm_fit)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 3.987694, Df = 1, p = 0.045834