반응형

 

[SAS] 피어슨 상관 계수 (Pearson's correlation coefficient) - PROC CORR

 

 세상의 많은 지표들은 서로 상관되어 있는 경우가 많다. 비흡연자에 비해 흡연자가 폐암에 잘 걸리고, 부모의 키가 클수록 아이의 키도 크다. 소득 수준이 높은 나라일수록 평균 수명도 길고, 수학 점수가 높은 학생은 물리 점수가 높으며, 식당에서 음식을 많이 시킬수록 계산해야 하는 금액도 커진다.

 그런데, 세상만사 칼같이 딱딱 떨어지는 게 아니다. 폐암에 걸리지 않는 꼴초도 있는가 하면, 부모의 키는 크지만 아이의 키는 작을 수 있다. 수학만 잘하고 물리는 못 할 수도 있다. 이렇게 예외가 있는 관계가 있는가 하면, 예외가 있을 수 없는 관계도 있다. $10,000$원짜리 삼겹살을 $x$인분 시키면 계산해야 하는 금액은 반드시 $10,000x$가 되는 게 일례가 될 수 있다. 

불확실성이 존재하는 관계

 

예외가 없는 관계

 

 예외가 많은 부모와 자녀의 키 문제로 돌아가 보자. 부모의 키로 자녀의 키를 유추하면 틀리는 경우가 많이 생기게 된다. 다시 말해, 불확실성이 높아진다. 반대로 삼겹살 주문량으로 예측한 계산 금액은 틀릴 수가 없다. 우리는 얼마나 정확하게 예측할 수 있는지를 계량하고자 한다. 그 방법 중 하나가 피어슨 상관 계수 (Pearson's correlation coefficient)다. 

 

피어슨 상관 계수 (Pearson's correlation coefficient)는 두 변수 사이의 상관관계의 정도를 나타내는 지표다. 이는 -1과 1 사이의 값을 갖는데 -1 혹은 1에 가까워질수록, 달리 표현하면 상관 계수의 절댓값이 커질수록, 즉 0에서 멀어질수록 "예측할 수 있는 정도"가 높아진다. 1에 가까워지는 것은 양의 상관관계인데 "예측도"가 증가하는 것을 의미하고, -1에 가까워지는 것은 음의 상관관계인데 "예측도"가 증가하는 것을 의미한다. 여기에서 예측하는 방법은 선형 방법이다. 즉, 한 변수의 값 $x$를 알 때, 적절한 실수 $a, b$에 대해 $y$값이 $y=ax+b$ 일 것이라고 예측할 수 있다는 것이다. 여기에서 피어슨 상관 계수에 대해 대표적으로 오해하는 점 두 가지를 들 수 있다. 

 

 1) 피어슨 상관 계수가 크면 기울기($a$값)이 크다. 기울기는 위 삼겹살 예시에서 1인분당 가격을 의미하는 것인데, 1인분 추가 주문 시 증가되는 전체 주문 금액으로 이해할 수 있다. 그런데 피어슨 상관 계수에는 이에 대한 정보가 전혀 없다. 피어슨 상관 계수는 (1) 양의 상관관계인지 음의 상관관계인지, (2) 얼마나 "예측도"가 큰지, 즉 데이터들이 얼마나 직선상에 모여있는지 에 대한 정보만 제공한다. 

 

 2) 피어슨 상관 계수가 0이면, 두 변수 사이에는 어떠한 상관 관계도 없다. 위에서 설명했듯이, 피어슨 상관 계수는 선형 방법으로 예측할 수 있는 정도를 표현한다. 따라서 선형 관계가 아닌 2차 함수 관계, 사인 함수 관계, 원형 함수 관계 등이 존재한다면, 두 변수 사이의 어떤 관련성은 분명히 존재하는 것이지만, 피어슨 상관 계수는 0일 수 있다. 

 

 본격적으로 피어슨 상관 계수를 구해보자. 

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

 

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

 

목표:  혈중 ALT와 AST에는 선형 상관관계가 존재하는가?

피어슨 상관 계수 (Pearson's correlation coefficient)를 구하는 코드는 다음과 같다.

 

코드

PROC CORR DATA=hong.df;
VAR AST ALT;
RUN;

PROC CORR DATA=hong.df; : 상관 분석 (CORRelation)을 시행하되 데이터는 hong 라이브러리의 df를 이용하라.
VAR AST ALT; : 변수는 AST와 ALT다. 둘 사이의 상관계수를 구하라.

 

결과

- 상관계수는 0.78869이다. 1에 가까울수록 양의 상관관계가 존재함을 나타내는데, 0.79 정도면 꽤 높은 "선형" 상관관계가 있다고 할 수 있다.

- 유의확률은 "<0.0001"으로 상당히 작으므로 매우 유의함을 알 수 있다. 따라서 귀무가설을 기각하고 대립 가설을 택해야 하는데, 이 검정의 귀무가설은 "상관 계수는 0"이라는 것이다. 따라서 상관 계수가 0은 아니라고 결론 내릴 수 있다.

 

얼마나 높아야 높다고 할 수 있는가?

위에서 0.79 정도면 꽤 높은 상관관계가 있다고 할 수 있다고 했는데, "얼마나 높은 것이 높은 것인가?"라는 의문이 들 수 있다. 즉, "how high is high?"라는 질문에 대한 답은 사실 정해져 있진 않다. 연구자마다 기준은 다른데, 그중 한 가지는 다음과 같다.

피어슨 상관 계수의 절댓값 $ ( \vert \rho \vert )$ 해석
$0.0<\vert \rho \vert\le0.2$ Very weak
$0.2<\vert \rho \vert\le0.4$ Weak
$0.4<\vert \rho \vert\le0.6$ Moderate
$0.6<\vert \rho \vert\le0.8$ Strong
$0.8<\vert \rho \vert\le1.0$ Very strong

 

피어슨 상관 계수 산출의 가정

 피어슨 상관 계수를 산출할 때, 기본적으로 가정은 존재하지 않는다. 수학적으로 표현하면 공분산을 각각의 표준 편차로 나눠준 것이니, 어떠한 가정도 필요하지 않다. 다만 해석에 유의해야 한다. 피어슨 상관 계수가 1 또는 -1에 가까운 값이 나왔을 때의 해석은 어디까지나 "선형 상관관계가 있다."는 것이다. 심지어 1에 가까운 값이더라도, 선형 상관관계보다 두 변수 사이의 관계를 더 잘 표현하는 관계가 있을 수 있다. 예를 들어, 어떤 변수 x와 y의 피어슨 상관 계수는 0.89로 매우 높은데, 둘 사이의 그래프는  아래와 같이 그려진다고 하자.

둘 사이의 상관 계수가 1(완벽한 선형 예측이 가능)은 아니지만, x와 y사이에는 이차 함수의 관계가 있는 것으로 보인다. "예측도"측면에서 상관 계수가 1은 아니므로 선형 예측은 완벽하지 않지만, 이차 함수의 관계로 보면 완벽하게 예측할 수 있는 관계다. 

 

또한, 피어슨 상관 계수가 0에 가까운 값이 나왔을 때의 해석은 어디까지나 "선형 상관관계가 없다."는 것이다.  예를 들어 상관계수가 0인 다음 두 변수 x와 y의 그래프는 다음과 같이 그려진다고 한다.

"예측도"측면에서 상관 계수가 0이므로 선형 예측이 불가능해야 하지만, 이차 함수의 관계로 보면 완벽하게 예측할 수 있는 관계다. 따라서 피어슨 상관 계수를 "예측도", "관련성"이라는 의미로 사용하려 한다면, 두 변수 사이의 산점도를 꼭 그려보는 것이 좋다.

 

두 변수 사이의 산점도 (플롯) 그리기.

산점도를 그리기위해서는 상관 계수를 구하는 코드에 한 마디만 추가하면 된다.

PROC CORR DATA=hong.df PLOTS=SCATTER;
VAR AST ALT;
RUN;

"PLOTS=SCATTER"라는 문구 하나를 추가하면 다음과 같은 산점도를 그려준다.

 

 

 

둘 사이에는 어떤 선형 관계가 있어 보인다. 그렇다면 우리가 산출한 피어슨 상관 계수는 AST와 ALT의 관련성을 잘 나타낸다고 할 수 있다. 즉, 이차 함수 관계 등 다른 비선형적 관계가 있는데도 불구하고 우리가 선형 상관 계수를 구해서 1보다 낮은 값이 나온 것이 아니라, 원래 선형 관계인데 어느 정도의 불확실성 때문에 1보다는 낮은 값이 나온 것이라고 결론지을 수 있다. 따라서 산점도를 꼭 그려보고 어떠한 관계 속에 있는지 확인하고 상관 계수를 해석하는 것은 매우 중요하다. 

 

가정이 정말 없는가?

 가정이 필요하다는 글, 책들이 정말 많을 것이다. 그럼에도 불구하고 나는 가정이 필요 없다는 말을 계속해왔다. 정말 가정이 없는 것일까? 답은 "그렇다"이다. 정말 없다. 하지만, Pearson's correlation coefficient가 적절할 때, 적절하지 않을 때는 분명히 존재한다. 피어슨 상관 계수는 두 변수가 연속형 변수일 때 좋다. 그리고 이상치(outlier)가 없을 때 정말 좋다. 이상치(outlier)는 피어슨 상관 계수의 값을 크게 변화시킬 수 있으며, 특히 표본의 수가 많지 않을 때 이런 일은 더 빈번하게 나타날 수 있다. 이럴 때에는 스피어만 상관 계수 (Spearman correlation coefficient)나 켄달의 타우 (Kendall's tau)가 더 적절할 수 있다. 이 두 가지 방식은 특히 순서형 변수일 때 더 큰 효용성을 갖는다. 

 아마도 위와 같은 이유로 피어슨 상관 계수 (Pearson's correlation coefficient)의 전제조건으로 정규성을 많이 이야기하는지도 모른다. 정규성을 따른다면 적어도 이상치로부터 어느 정도 자유로울 수 있고, 왜곡된 상관 계수가 산출되지 않을 수 있다.

 

[SAS] 피어슨 상관 계수 (Pearson's correlation coefficient) 정복 완료!

 

작성일: 2023.04.10.

최종 수정일: 2023.10.17.

이용 프로그램: SAS v9.4

운영체제: Windows 11

반응형
반응형

[SAS] 크루스칼 왈리스 검정 (비모수 일원 배치 분산 분석: Kruskal-Wallis Test) - NPAR1WAY()

 

셋 이상의 분포에서 평균이 다른지 확인하는 방법은 ANOVA (ANalysis Of VAriance)였다. 

2022.10.12 - [모평균 검정/SAS] - [SAS] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - PROC GLM, PROC ANOVA

2022.11.26 - [모평균 검정/SAS] - [SAS] 등분산 가정이 성립하지 않는 일원 배치 분산 분석 (Welch's ANOVA) - PROC GLM, PROC MIXED

 

그런데, 여기에는 각각의 분포가 정규성을 따른다는 중요한 가정이 필요하다. 하지만 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 크루스칼-왈리스 검정(Kruskal-Wallis Test)이다.  이번 포스팅에서는 크루스칼 왈리스 검정과 그 사후 분석 (post hoc analysis)까지 알아보고자 한다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

목표:  흡연 상태(SMOK)에 따라 니코틴 중독 점수 (NICOT_ADDICT) 중앙값이 모집단 수준에서 서로 다르다고 말할 수 있는가?

 

 만약 NICOT_ADDICT의 분포가 정규성을 따른다면, ANOVA로 이를 확인할 수 있을 것이다. 따라서 정규성 여부를 먼저 확인해 보자.

정규성 - 코드

 - 정규성 검정에 관한 분석 내용은 다음 글에서 살펴볼 수 있다. 2022.08.12 - [기술 통계/SAS] - [SAS] 정규성 검정 - PROC UNIVARIATE

 - PROC UNIVARIATE명령문에서 "CLASS"옵션을 사용하는 방법은 다음 글에서 살펴볼 수 있다. 2022.09.23 - [기술 통계/SAS] - [SAS] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등) - PROC UNIVARIATE, PROC MEANS

PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS SMOK;
VAR NICOT_ADDICT;
HISTOGRAM NICOT_ADDICT/ NORMAL (MU=EST SIGMA=EST);
RUN;

정규성 - 결과

결과 - 1) 비흡연자 (SMOK=0)

표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 right-skewed 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.

 

결과 - 2) 과거 흡연자 (SMOK=1)

표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 right-skewed 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.

 

결과 - 3) 현재 흡연자 (SMOK=2)

표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 uniform distribution 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.

 

따라서 일원 배치 분산 분석 (ANOVA)을 시행할 수 없고, 크루스칼 왈리스 분석 (Kruskal-Wallis Test)을 시행해야 한다.

 

지금까지 많은 분석을 할 때 분산이 같은지 여부 (등분산성)을 중요하게 여겼었다. 하지만 크루스칼 왈리스 분석은 분산이 달라도 어느 정도 괜찮다. 물론 분산에 아주 큰 차이가 나면 위험해지지만, 어느 정도는 robust 하므로 웬만하면 그냥 써도 괜찮다.

 

크루스칼 왈리스 (Kruskal Wallis) 검정 코드

PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;

코드는 두 군에서의 비모수 중앙값 검정인 맨 휘트니 U 검정과 같다.(2022.12.20 - [모평균 검정/SAS] - [SAS] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) - NPAR1WAY())

 

PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO; : 비모수 검정을 시행한다. 데이터는 hong라이브러리의 df를 사용하고, 크루스칼 왈리스 검정을 시행한다. 연속성 수정은 하지 않는다.(연속성 수정에 관한 내용은 다음 링크를 확인하길 바란다.2022.08.30 - [통계 이론] - [이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity))
CLASS SMOK; : SMOK값에 따라 검정을 시행한다. (흡연 여부에 따라 시행한다.)
VAR NICOT_ADDICT; : NICOT_ADDICT값을 검정한다.

 

결과

P-value가 <0.0001이므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

 귀무가설:  $H_0=$ 세 집단의 모평균은 "모두" 동일하다.

 대립 가설:  $H_1=$ 세 집단의 모평균이 모두 동일한 것은 아니다.

 

우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.

이 말은, 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수는 $_3 C_2 =3$이다.

 (1) 비흡연자 vs 과거 흡연자

 (2) 비흡연자 vs 현재 흡연자

 (3) 과거 흡연자 vs 현재 흡연자

 

비모수 검정인 Kruskal Wallis 검정의 사후 분석으로 쓰이는 방법은 다음 네 가지가 있다.

1) Pairwise Wilcoxon Test

 - 가장 적절하지 않다. 왜냐하면 Kruskal wallis 검정은 기본적으로 순위를 이용한 검정 방법인데, 전체 집단에서 특정 피험자의 순위는 부분 집단에서의 순위와 다르기 때문이다. 예를 들어 전체 집단에서 흡연자인 어떤 피험자의 니코틴 중독 점수의 순위는 (비흡연자 + 현재 흡연자) 집단에서의 순위와 크게 달라진다. 동시에 분산도 크게 달라지기 때문에 적절하지 않다.

2) Conover Test

 - 상기 문제점을 해결하였고, t 분포에 근거해 판단한다.

3) Dunn Test

 - 상기 문제점을 해결하였고, z 분포에 근거해 판단한다.

4) Nemenyi Test

 

하지만 SAS는 위 네 개 중 하나의 방법이 아닌 다른 방법을 사용한다. 

Dwass, Steel, Critchlow-Fligner (DSCF) multiple comparison procedure를 사용하며 이는 pairwise Wilcoxon test를 기반으로 하나, 특수한 방법으로 그 familywise error rate를 보정하고 있다. 이 방법은 Conover Test보다 더 우수한 것으로 받아들여지고 있다. 

 

DSCF 사후 검정 - 코드

PROC NPAR1WAY DATA=hong.df WILCOXON DSCF CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;

일반적인 Kruskal Wallis 검정을 시행하되, 첫 번째 줄에 "DSCF"를 추가해 준 것이 끝이다.

 

결과

여러 값이 제시되지만, 마지막 열인 "Pr>DSCF"인 p-value만 보면 된다. 

 

 

 (1) 과거 흡연자(SMOK=1) vs 현재 흡연자 (SMOK=2) : p-value < 0.0001

 (2) 과거 흡연자 (SMOK=1) vs 비흡연자 (SMOK=0) : p-value < 0.0001

 (3) 현재 흡연자 (SMOK=2) vs 비흡연자 (SMOK=0) : p-value < 0.0001

 

모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.

비흡연자 (0)와 과거 흡연자 (1)  사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

비흡연자 (0)와 현재 흡연자 (2) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

과거 흡연자 (1)와 현재 흡연자 (2)) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

 

어떤 사후 분석을 쓸 것인가

 이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, Kruskal-Wallis test에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.

 

왜 굳이 Kruskal Wallis test를 쓰는 것인가?

이 시점에서 이런 의문이 들 수 있다.

"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 Kruskal Wallis test라는 방법까지 사용하는 것인가?"

 아주 논리적인 의문점이다. 하지만, 반드시 Kruskal-Wallis test를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 $\alpha=0.05$로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $\left( 1-0.05 \right)^3 \approx 0.86$이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)

 그런데, Kruskal-Wallis test의 귀무가설은 "모든 집단의 중앙값이 같다."이다. 따라서 모든 집단의 중앙값이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 $1-0.86=0.14$가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.

 따라서, 각각을 비교해 보는 것이 아니라 한꺼번에 비교하는 Kruskal-Wallis test를 시행해야 함이 마땅하다. 

 

여기에서 한 번 더 의문이 들 수 있다.

 "사후 분석을 할 때에는 1종 오류가 발생하지 않는가?"

그렇다. 1종 오류가 발생할 확률이 있으므로, p-value의 기준을 더 엄격하게 (0.05보다 더 작게) 잡아야 한다. P-value를 보정하는 방법은 일반적으로는 Bonferroni, holm, hochberg, hommel, BH, BY, FDR 등이 있고, Kruskal-Wallis test의 사후 분석 방법으로는 특별히  Pairwise Wilcoxon Test, Conover Test, Dunn Test, Nemenyi Test, DSCF을 사용하고 있다. 

 

코드 정리

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

*정규성 검정하기;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS SMOK;
VAR NICOT_ADDICT;
HISTOGRAM NICOT_ADDICT/ NORMAL (MU=EST SIGMA=EST);
RUN;

*크루스칼 왈리스 검정;
PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;

*크루스칼 왈리스 검정 & 사후검정;
PROC NPAR1WAY DATA=hong.df WILCOXON DSCF CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;

[SAS] 크루스칼 왈리스 검정 (비모수 일원 배치 분산 분석: Kruskal-Wallis Test) 정복 완료!

 

작성일: 2023.04.05.

최종 수정일: 2023.04.05.

이용 프로그램: SAS v9.4

운영체제: Windows 11

반응형
반응형

[R] 선형 회귀 분석의 모델 최적화 : 전진 선택법, 후진 제거법 - lm(), step()

 

 앞선 글에서 여러 개의 변수를 넣고 선형 회귀 분석을 하는 방법까지 소개하였다. 하지만 여기에는 큰 문제가 있는데, 바로 "어떤 교란 변수까지 모델에 넣을 것인가?"라는 질문이다. 이에 대한 방법은 여럿이 소개되어 있다.

1) 단변량 분석에서 유의하게 나온 변수를 모델에 넣는다. 이때 유의 수준은 기존의 0.05보다 덜 보수적인 0.1 혹은 0.2로 책정될 수도 있다.

2) 본인이 갖고 있는 데이터로 시행한 단변량 분석에서 유의하지 않더라도 과거 문헌상 고려되어야만 하는 교란 변수라면 공변량으로 넣는다.

3) 하나의 공변량을 추가로 넣을 때 Estimate를 크게 변화시키지 않으면서 신뢰구간만 넓히게 된다면 굳이 공변량을 넣지 마라.

4) Interaction term을 쓴다면 Hierachically well formulated model을 쓰기 위해 가능한 모든 변수들의 조합을 공변량으로 넣어라.

5) 전진 선택법 혹은 후진 제거법을 이용하여 변수를 선택하라.

 

당장 생각나는 것은 이 정도인데 아마도 더 많은 방법론들이 존재할 것이다. 이번 포스팅에서는 위 중 5번째에 속하는 "전진 선택법(Forward selection)"과 "후진 제거법(Backward elimination)"을 소개할 것이다.

 

전진 선택법과 후진 제거법

 전진 선택법은 어떠한 독립 변수도 없는 모델 (null model)에서 시작하여 가장 유의한 변수부터 하나씩 추가하며 AIC (Akaike's Information Criterion)가 최소화되는 조합을 찾아내는 방법이다. AIC는 회귀 분석에서 독립 변수들이 종속 변수를 얼마나 잘 설명하는지 평가하는 동시에 모델이 얼마나 간단한지도 평가하는 도구인데 자세한 내용은 다음 링크를 확인하길 바란다. (링크 추가 예정) 

 후진 제거법은 분석가가 지정한 모든 독립 변수를 포함한 모델에서 시작하여 제일 쓸모없는 변수를 하나씩 제거해보며 AIC를 최소화하는 조합을 찾아내는 방법이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  심혈관 질환 위험 점수 (CVD_RISK)를 예측하는 변수로 수축기 혈압 (SBP), BMI, 흡연 상태 (SMOK)를 고려하되, 전진 선택법과 후진 제거법을 이용하여 최적의 모델을 구성하는 공변량의 조합을 찾아라.

이전에 SBP, BMI, SMOK를 고려한 선형 회귀 분석에서 SMOK는 유의하지 않은 결과를 냈다(2022.12.23 - [선형 회귀 분석/R] - [R] 다중 선형 회귀 분석 (Multiple linear regression) - lm(), factor()). 전진 선택법과 후진 선택법이 제대로 작동한다면 SMOK는 최종 모델에 들어가지 않는 것이 적절할 것이다. 과연 그런지 확인해보도록 하겠다.

1) 전진 선택법

코드

전진 선택법을 위해서는 어떠한 독립 변수도 없는 모델(null model)을 먼저 지정해주어야 한다.

Null<-lm(CVD_RISK~1, data=df)

Null<-lm(CVD_RISK~1, data=df) : 종속 변수에 CVD_RISK를 넣고 독립 변수에는 어떠한 변수도 없는 선형 회귀 모델을 구축하여 Null이라는 객체에 저장한다. 종속 변수에 1이라고 적으면 어떠한 독립 변수도 존재하지 않는 모델을 의미하는 것이다.

 

#SMOK는 범주형 변수임을 미리 선언한다.
df$SMOK<-factor(df$SMOK)

Forward<-step(Null, direction="forward", scope=(CVD_RISK~SBP+BMI+SMOK))
summary(Forward)

Forward<-step(Null, direction="forward", scope=(CVD_RISK~SBP+BMI+SMOK))

  : Null이라는 모델을 가져와서, 전진 선택법(direction="forward")으로 변수를 선택할 것인데, 아무리 변수가 많이 선택되어도 SBP, BMI, SMOK까지만 들어간다. 최종 모델을 Forward에 저장한다.
summary(Forward) : 최종 모델을 보여달라.

 

결과

Call:
lm(formula = CVD_RISK ~ SBP + BMI, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3967  -3.7096   0.0512   3.4044  15.9551 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.36106    1.83287   7.835 1.19e-14 ***
SBP          0.99389    0.01931  51.478  < 2e-16 ***
BMI          1.08243    0.17174   6.303 4.38e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.079 on 997 degrees of freedom
Multiple R-squared:  0.9393,	Adjusted R-squared:  0.9392 
F-statistic:  7712 on 2 and 997 DF,  p-value: < 2.2e-16

맨 윗줄 "Call:"을 보면, SBP와 BMI만 모델에 들어간 것을 알 수 있다.

 

 

2) 후진 제거법

코드

Full<-lm(CVD_RISK~SBP+BMI+SMOK, data=df)
Back<-step(Full, direction="backward")
summary(Back)

Full<-lm(CVD_RISK~SBP+BMI+SMOK, data=df) : 먼저, 다중 선형 회귀 분석을 이전 포스팅과 같이 시행하겠다. (2022.12.23 - [선형 회귀 분석/R] - [R] 다중 선형 회귀 분석 (Multiple linear regression) - lm(), factor())
Back<-step(Full, direction="backward") :  Full이라는 모델을 가져와서, 후진 제거법(direction="backward")으로 변수를 제거하며 최적 모델을 선정하고 Back에 저장한다.
summary(Back) : 최종 모델을 보여달라.

 

결과

Call:
lm(formula = CVD_RISK ~ SBP + BMI, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3967  -3.7096   0.0512   3.4044  15.9551 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.36106    1.83287   7.835 1.19e-14 ***
SBP          0.99389    0.01931  51.478  < 2e-16 ***
BMI          1.08243    0.17174   6.303 4.38e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.079 on 997 degrees of freedom
Multiple R-squared:  0.9393,	Adjusted R-squared:  0.9392 
F-statistic:  7712 on 2 and 997 DF,  p-value: < 2.2e-16

전진 선택법과 같은 결과를 내주는 것을 볼 수 있다. 

 

하지만 두 방법이 꼭 같은 결과를 내는 것은 아니다. 서로 다른 결과를 낼 수 있음에 유의하길 바란다. 

 

[R] 선형 회귀 분석의 모델 최적화 : 전진 선택법, 후진 제거법 정복 완료!

작성일: 2022.12.23.

최종 수정일: 2022.12.23.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 다중 선형 회귀 분석 (Multiple linear regression) - lm(), factor()

 

 지난 포스팅에서 연속형 변수의 단순 선형 회귀 분석 (Simple linear regression: 2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()), 범주형 변수의 단순 선형 회귀 분석 (Simple linear regression with categorical variables:2022.12.22 - [선형 회귀 분석/R] - [R] 범주형 변수의 선형 회귀 분석 (Simple linear regression with categorical variable) - lm(), factor())에 대해 알아보았다.  위 분석들은 모두 독립 변수가 1개인 단순 선형 회귀 분석이다. 하지만 세상 일은 그렇게 단순하지 않다. 

 

교란 변수 (Confounder, confounding variable)

 우리는 앞서  연속형 변수의 단순 선형 회귀 분석 (Simple linear regression: 2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()) 포스팅에서 수축기 혈압과 심혈관 질환의 관련성을 찾아보았다. 그런데, 비만인 사람은 혈압이 대체적으로 높은 것으로 알려져 있다. 더구나 비만인 사람은 심혈관 질환의 위험도도 높아지는 것으로 알려져 있다. 그렇다면, "수축기 혈압이 높았던 사람에게서 높은 심혈관 질환 위험 점수가 발견되었던 것은 그냥 그 사람이 비만도가 높아 당연히 수축기 혈압도 높고, 심혈관 질환 위험 점수도 높았던 것이 아닐까?"라는 의문이 든다. 즉, "우리가 이전 단순 선형 회귀 분석에서 발견한 과학적 현상(수축기 혈압이 높으면 심혈관 질환 위험 점수가 높아진다.)은 사실은 비만 때문이 아니었을까?"라는 의문이 들게 된다. 이렇게 독립변수와 종속변수 모두의 원인이 되는 변수를 교란 변수라고 하며, 이런 현상(교란)은 통계 분석을 할 때 반드시 교정을 해야 한다. (사실 독립변수는 인과적 관계가 아니라 연관 관계이기만 해도 충분하다.)

교란 현상의 해결 방법: 다중 회귀 분석

 교란 현상을 해결하는 방법은 보정(adjustment)과 층화(stratification)가 있다. 층화도 좋은 방법이고, 때로는 보정을 할 수 없어 층화를 해야만 할 때도 있지만, 보통 표본의 수가 부족해지는 문제가 있기 때문에, 층화가 아니라 보정을 시행하곤 한다. 보정을 시행할 수 있는 방법으로는 이전에 다룬 편상관 분석(2022.12.18 - [상관분석/R] - [R] 편상관 계수, 부분 상관 계수 (Partial correlation coefficient) - pcor.test())이 있다. 하지만 이는 상관 계수만 구할 수 있어 "예측할 수 있는 정도"만 표현할 수 있다. 그 대신 비만도의 영향을 보정했을 때, 수축기 혈압이 1단위 증가할 때 심혈관 질환 위험 점수가 얼마나 증가하는지(기울기)를 알 수 있는 방법으로는 본 포스팅에서 다룰 다중 선형 회귀 분석이 있다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  BMI의 영향을 배제 혹은 교정했을 때 수축기 혈압(SBP)이 1단위 증가할 때 심혈관 질환 위험 점수(CVD_RISK)는 얼마나 증가하는가?

 

단순 선형 회귀 분석과 같이 lm()함수를 사용하며, 코드는 다음과 같다.

코드

LR1<-lm(CVD_RISK~SBP+BMI, data=df)
summary(LR1)

LR1<-lm(CVD_RISK~SBP+BMI, data=df) : df 데이터를 사용하여 선형 회귀 분석을 시행한다. 종속 변수는 CVD_RISK로, 독립 변수에는 SBP와 BMI를 사용하라. 그리고 그 결과를 LR1에 저장하라. SBP와 BMI의 순서가 바뀌는 것은 아무 상관이 없다. 즉, BMI를 보정했을 때의 SBP의 영향과 SBP를 보정했을 때의 BMI의 영향을 동시에 파악하는 것이다.
summary(LR1) : LR1을 요약하여 보여달라.

 

결과

Call:
lm(formula = CVD_RISK ~ SBP + BMI, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3967  -3.7096   0.0512   3.4044  15.9551 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.36106    1.83287   7.835 1.19e-14 ***
SBP          0.99389    0.01931  51.478  < 2e-16 ***
BMI          1.08243    0.17174   6.303 4.38e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.079 on 997 degrees of freedom
Multiple R-squared:  0.9393,	Adjusted R-squared:  0.9392 
F-statistic:  7712 on 2 and 997 DF,  p-value: < 2.2e-16

해석의 방법은 단순 선형 회귀 분석에서와 같다. (2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm())

 

Coefficients:
                    Estimate   Std. Error t value     Pr(>|t|)    
(Intercept)  14.36106    1.83287   7.835  1.19e-14 ***
SBP              0.99389    0.01931  51.478   < 2e-16 ***
BMI               1.08243     0.17174   6.303 4.38e-10 ***

 

 BMI의 영향을 보정하더라도, SBP가 1단위 증가할 때 CVD_RISK는 0.99389만큼 증가하며 이에 대한 p-value는 $2 \times 10^{-16}$보다 작아 매우 유의함을 알 수 있다. 

 또한, SBP의 영향을 보정하더라도, BMI가 1단위 증가할 때 CVD_RISK는 1.08243만큼 증가하며 이에 대한 p-value는 $4.38 \times 10^{-10}$보다 작아 매우 유의함을 알 수 있다. 

 

즉, 단순 회귀 분석에서 나타난 수축기 혈압과 심혈관 질환 위험 점수 사이의 유의한 관련성은 BMI의 영향을 보정하더라도 유의함을 알 수 있다. 또한, 단순 회귀 분석에서의 회귀 계수  1.101935보다 작으므로, 일부분의 설명력은 BMI에게 빼았겼음을 알 수 있다.

 

 

보정 변수가 두 개 이상일 때 코드

 BMI외에 흡연 상태로도 보정하고 싶을 수 있다. 이럴 땐 어떻게 해야 할까? 답은, 단순히 BMI 뒤에 추가적인 보정변수를 '+'연산자를 이용하여 나열하면 된다. 그런데, 선형 회귀 분석 전에 연속형 변수와 범주형 변수의 전처리가 다르다는 것을 일전에 소개한 적이 있다. 범주형 변수는 분석 전에 범주형 자료임을 선언해주어야 한다. 내용은 다음 링크를 확인하면 된다.

연속형 변수: 2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()

범주형 변수: 2022.12.22 - [선형 회귀 분석/R] - [R] 범주형 변수의 선형 회귀 분석 (Simple linear regression with categorical variable) - lm(), factor()

 그런데, 흡연 상태 (SMOK)는 범주형 자료다. 따라서 선형 회귀 분석을 시행하기 전에 범주형임을 선언해주어야 한다.

df$SMOK<-factor(df$SMOK)
LR2<-lm(CVD_RISK~SBP+BMI+SMOK, data=df)
summary(LR2)

 

결과

Call:
lm(formula = CVD_RISK ~ SBP + BMI + SMOK, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.280  -3.663   0.133   3.454  16.141 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  15.4269     1.9747   7.812 1.42e-14 ***
SBP           0.9842     0.0205  48.009  < 2e-16 ***
BMI           1.0832     0.1717   6.310 4.21e-10 ***
SMOK1         0.1361     0.4189   0.325    0.745    
SMOK2         0.7824     0.5036   1.554    0.121    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.077 on 995 degrees of freedom
Multiple R-squared:  0.9394,	Adjusted R-squared:  0.9392 
F-statistic:  3859 on 4 and 995 DF,  p-value: < 2.2e-16

SBP, BMI는 유의한 결과를 보이고 있지만, SMOK는 1(과거 흡연자), 2(현재 흡연자) 모두 유의한 결과를 보이고 있지 않다. 

"BMI와 흡연 상태로 보정을 하였을 때, 수축기 혈압 1 단위 증가할 때, CVD_RISK는 0.9842점 증가한다"라고 할 수 있지만

"SBP와 BMI로 보정을 하였을 때, 비흡연자와 과거 흡연자 혹은 비흡연자와 현재 흡연자 사이의 CVD_RISK의 차이는 존재하지만 유의하다고 할 수 없다."라고 결론 내리게 된다.

 

 코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

##다중 선형 회귀 분석 1
LR1<-lm(CVD_RISK~SBP+BMI, data=df)
summary(LR1)

##다중 선형 회귀 분석 2
df$SMOK<-factor(df$SMOK)
LR2<-lm(CVD_RISK~SBP+BMI+SMOK, data=df)
summary(LR2)

 

[R] 다중 선형 회귀 분석 (Multiple linear regression) 정복 완료!

작성일: 2022.12.23.

최종 수정일: 2022.12.23.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 범주형 변수의 선형 회귀 분석 (Simple linear regression with categorical variable) - lm(), factor()

 

지난 포스팅에서 우리는 수축기 혈압(SBP)과 심혈관 질병 위험 점수 (CVD_RISK)의 선형 회귀 분석에 대해 알아보았다.(2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()) 이 분석에 쓰인 독립 변수인 수축기 혈압 (SBP)는 연속형 변수다. 만약 범주형 변수로 선형 회귀 분석을 시행하려 한다면 어떻게 해야 할까? 예를 들어 흡연 여부(SMOK: 비흡연자 or 과거 흡연자 or  현재 흡연자)로 심혈관 위험 점수를 예측하려 할 땐 어떻게 해야 할까? 이번 포스팅에서는 이에 대해 알아볼 것이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  흡연 여부(SMOK)로 심혈관 위험 점수(CVD_RISK)를 예측할 수 있는가?

이전 단순 선형 회귀 분석에서 배운 대로 선형 회귀 분석을 시행하면 어떻게 될까? 자료의 유형 (연속형 vs 범주형)을 고려하지 않은 책 우선 분석을 시행해보자. 분석 방법은 이전 포스팅(2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm())을 보면 알 수 있다.

 

코드

LR_SMOK_CVD_ft<-lm(CVD_RISK~SMOK, data=df)
summary(LR_SMOK_CVD_ft)

 

결과

Call:
lm(formula = CVD_RISK ~ SMOK, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.581 -11.818  -0.802  11.130  54.669 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 159.9106     0.7216  221.61   <2e-16 ***
SMOK         15.2974     0.6421   23.82   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.45 on 998 degrees of freedom
Multiple R-squared:  0.3626,	Adjusted R-squared:  0.3619 
F-statistic: 567.6 on 1 and 998 DF,  p-value: < 2.2e-16

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept) 159.9106     0.7216  221.61   <2e-16 ***
SMOK           15.2974     0.6421   23.82   <2e-16 ***

 

주목해야할 건 "SMOK"이다. 이전 연속형 독립 변수로 시행한 선형 회귀 분석의 결과 해석은 "독립 변수가 1 단위 증가할 때 종속 변수의 변화량"="기울기"이다. 그러면 이 결과를 해석하면 어떻게 될까? 참고로 SMOK는 다음과 같은 변수다.

SMOK 해석
0 비흡연자
1 과거 흡연자
2 현재 흡연자

 SMOK가 1단위 증가할 때 CVD_RISK의 변화량은 "비흡연자에 비해 과거 흡연자가", "과거 흡연자에 비해 현재 흡연자가" 얼마나 더 높은 CVD_RISK의 값을 갖는지 표현하는 것이다. 근데 이 값은 15.2974로 같다. 즉, 비흡연자에서 과거 흡연자로 변하는 것과, 과거 흡연자에서 현재 흡연자로 변화하는 것이 같은 수준이라고 가정하고 있는 것이다. 이 말은 SMOK를 등간척도로 보고 있다는 말이다. 등간척도는 연속 변수의 일종이니, R은 현재 SMOK를 연속변수로 확인하고 있다는 것이다. 과연 R이 SMOK를 연속변수로 확인하고 있는지 확인해 보자. (변수의 유형을 확인하는 방법은 다음 링크에서 확인할 수 있다.2022.11.21 - [통계 프로그램 사용 방법/R] - [R] 변수의 유형 (타입, type) 확인 및 변경 - as.factor(), as.numeric(), as.character(), str())

코드

str(df$SMOK)

결과

num [1:1000] 1 2 0 0 2 1 1 1 1 0 ...

역시나 numeric (연속형 변수)로 인지하고 있다. 그도 그럴 것이 0,1,2로 구성된 변수를 R이 무슨 수로 범주형 변수인지, 연속형 변수인지 알겠는가. 따라서 우리가 직접 지정해주어야만 한다. 범주형 변수로 변환시켜 주겠다. 변수의 유형을 변환하는 방법에 대해서는 다음 링크(2022.11.21 - [통계 프로그램 사용 방법/R] - [R] 변수의 유형 (타입, type) 확인 및 변경 - as.factor(), as.numeric(), as.character(), str())에서 확인할 수 있다. 이 링크의 말미에 간단히 언급한 factor() 함수를 사용할 것이다.

 

코드

df$SMOK<-factor(df$SMOK)

 

코드

이제 다시 단순 선형 회귀 분석을 시행해보자.

LR_SMOK_CVD<-lm(CVD_RISK~SMOK, data=df)
summary(LR_SMOK_CVD)

결과

Call:
lm(formula = CVD_RISK ~ SMOK, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-47.33 -11.61  -0.49  10.87  55.23 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 159.3539     0.7629  208.88   <2e-16 ***
SMOK1        17.6012     1.2229   14.39   <2e-16 ***
SMOK2        30.0866     1.3021   23.11   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.42 on 997 degrees of freedom
Multiple R-squared:  0.3657,	Adjusted R-squared:  0.3644 
F-statistic: 287.4 on 2 and 997 DF,  p-value: < 2.2e-16

Coefficients:
                      Estimate  Std. Error t value Pr(>|t|)    
(Intercept)  159.3539     0.7629  208.88   <2e-16 ***
SMOK1           17.6012     1.2229     14.39   <2e-16 ***
SMOK2         30.0866     1.3021      23.11   <2e-16 ***

 

결과를 보면 범주형 변수로 바꾸기 전과 다르게, SMOK에 대한 결괏값이 2개가 나온다. SMOK1과 SMOK2다. SMOK는 0,1,2 세 개의 값을 갖는데, SMOK=1인 그룹과 SMOK=2인 그룹에 대한 통계 검정 결과를 나타내는 것이다. 그렇다면 왜 0은 없을까? 0은 참조값(reference)으로 쓰인 것이다. 즉 비흡연자(SMOK=0)에 비해 과거 흡연자 (SMOK=1)는 심혈관 질환 위험 점수 (CVD_RISK)가 17.6012점 높으며, p-value는 매우 작으므로 이는 유의미한 차이라고 할 수 있다. 또한, 비흡연자(SMOK=0)에 비해 현재 흡연자 (SMOK=2)는 심혈관 질환 위험 점수 (CVD_RISK)가 30.0866점 높으며, p-value는 매우 작으므로 이는 유의미한 차이라고 할 수 있다. 

 

 R은 가장 작은 값을 자동으로 참조값으로 잡는다. 그래서 SMOK=0인 비흡연자가 reference로 잡혔다. 만약, 참조값(reference)으로 비흡연자가 아닌 다른 군을 지정하고 싶다면 어떻게 할까? 바로, 범주형 변수로 지정할 때, levels라는 argument의 맨 앞자리에 원하는 숫자를 적으면 된다. 만약 과거 흡연자 (SMOK=1)를 reference로 잡고자 한다면, levels=c(1,0,2)와 같이 적으면 된다. 마지막 0과 2의 순서는 바뀌어도 reference와는 아무 상관이 없고, 다만 선형 회귀 분석 결과가 나타나는 순서를 지정하는 역할만 한다. 다음 예시를 보자.

 

코드

#범주형 변수 지정
df$SMOK<-factor(df$SMOK, levels=c(1,0,2))

#선형회귀분석 시행
LR_SMOK_CVD<-lm(CVD_RISK~SMOK, data=df)
summary(LR_SMOK_CVD)

 

결과

Call:
lm(formula = CVD_RISK ~ SMOK, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-47.33 -11.61  -0.49  10.87  55.23 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 176.9552     0.9557  185.15   <2e-16 ***
SMOK0       -17.6012     1.2229  -14.39   <2e-16 ***
SMOK2        12.4853     1.4237    8.77   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.42 on 997 degrees of freedom
Multiple R-squared:  0.3657,	Adjusted R-squared:  0.3644 
F-statistic: 287.4 on 2 and 997 DF,  p-value: < 2.2e-16

SMOK0과 SMOK2에 대한 내용만 있고, SMOK1은 없어졌다. 만약 levels=c(1,2,0)이라고 적었다면 두 행의 순서가 뒤바뀌어 나왔을 것이다.

 

 

이분형 변수는 범주형 지정이 필요 없다.

이분형 변수는 범주형 변수로의 변환이 굳이 필요는 없다. 왜냐하면 어차피 값이 두 개밖에 존재하지 않기 때문에, 연속형이든 범주형이든 같은 결과를 내기 때문이다. 하지만 이 경우 반드시 reference가 0이 되므로 만약 1을 reference로 잡고 싶은 경우, 범주형 변수로 변경하며 reference를 지정해주면 된다.

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

##범주형 지정
df$SMOK<-factor(df$SMOK)

##선형 회귀 분석 시행
LR_SMOK_CVD<-lm(CVD_RISK~SMOK, data=df)
summary(LR_SMOK_CVD)

 

[R] 범주형 변수의 선형 회귀 분석 (Simple linear regression with categorical variable) 정복 완료!

작성일: 2022.12.22.

최종 수정일: 2022.12.22.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 선형 회귀 분석의 전제 조건 - lm(), plot(), gvlma()

 

 지난 포스팅에서 단순 선형 회귀 분석을 소개하며, 전제 조건(가정)에 대한 논의가 필요하다고 했다. (단순 선형 회귀 분석:2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()) 보통 네 가지 전제조건을 의미하며, 그중 통계적으로 검정하는 조건은 세 가지다.

 

1) 정규성(Normality) : 잔차의 정규성

2) 선형성(Linearity) : 독립변수와 종속변수의 선형성

3) 잔차의 등분산성 (Homoscedasticity) : 잔차의 분산은 종속변수의 값에 따라 달라지지 않는다.

4) 독립성 (Independence): 관측 간의 독립성

 

이 중 4번 "독립성"은 통계적으로 검정하진 않는다. 독립성이란 개별 관측 간에는 연관되어있지 않다는 뜻이다. 대표적인 사례가 반복 측정이다. 같은 사람에서 약물 복용 전 간기능 수치와 약물 복용 후 간기능 수치를 비교한다면 두 개의 값은 "한 사람"에서 나왔기 때문에 독립이라고 볼 수 없다. 이런 경우 독립성을 만족한다고 볼 수 없다. 이런 경우 독립 표본 t 검정이 아닌 대응 표본 t 검정을 사용하는 것과 일맥상통한다. 유전적 정보를 공유하는 가족이나, 같은 경험을 하는 같은 반 학생들은 (항상 그런 것은 아니지만) 가끔 일반적인 선형 회귀 분석이 적절하지 않을 수 있다.

 

나머지 세 가지 전제조건에 대해 검증을 해보도록 하겠다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  수축기 혈압(SBP)과 심혈관 질환 위험 점수 (CVD_RISK) 사이의 선형 회귀 분석의 전제조건이 성립하는가?

 

먼저 선형 회귀 분석을 실시한다. (2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm())

LR_SBP_CVD<-lm(CVD_RISK~SBP, data=df)

 

R은 참 고맙게도, plot() 함수 안에 선형 회귀 분석을 실시한 객체를 넣어주면 전제 조건 확인을 위한 그래프를 바로 보여준다.

 

코드

plot(LR_SBP_CVD)

실시 후 엔터(Window) 혹은 Return(Mac)을 누르면 4개의 그래프를 보여주며, 우리는 마지막 그림을 제외하고 앞 세 개의 그래프만 사용할 것이다.

 

결과

1) 잔차의 정규성 (Normality)

2번째 나오는 그림을 이용한다.

일반 Q-Q plot처럼 해석하면 된다. (Q-Q plot내용은 다음 포스팅을 참고하길 바란다.2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()) 대부분의 점이 45도의 점선 상에 있는 것을 알 수 있다. (비록 $x$축과 $y$축의 scale이 달라 45도처럼 보이지 않으나, $y=x$직선이다.) 따라서 잔차는 정규성을 따른다고 할 수 있다.

 

2) 선형성 (Linearity)

첫 번째 나오는 그림을 이용한다.

해석 방법: Residual이 0 위아래에 고루 분포하며 특정한 형태나 패털을 보이지 않는다. 즉, 점선 위나 아래에 몰려있는 공간이 있으면 안 된다. 예를 들어 Fitted Values가 140인 근방에서는 Residuals이 0보다 큰 게 많고, 200인 근방에서는 0보다 작은 게 많은 상황이 발생했다면 선형성이 망가진 것이다. 

 

 여기서 유의할 것은 절대로 일반적인 산포도로 선형성을 확인하면 안 된다는 것이다. 이유로는 다음 두 가지를 들 수 있다.

 1) 독립변수가 여러 개인 다중 선형 회귀 분석에서는 차원의 문제로 산포도로 선형성을 확인하는 것이 쉽지 않다.

 2) 단순한 변수의 산포도를 확인하면 안 된다. 예를 들어 A와 B의 선형 회귀 분석을 시행할 것인데 A는 로그 변환하여 사용하기로 했다면 산포도는 A와 B의 산포도가 아니라, log(A)와 B의 산포도를 확인해야 한다. 

 

3) 등분산성 (Homoscedasticity)

세 번째 나오는 그림을 사용한다.

해석 방법: 점들이 어떤 직사각형 안에 고루 퍼져있다면 등분산성을 만족한다. 그러므로 본 사례는 등분산성을 만족한다.

 설명의 용이성을 위해 단순 선형 회귀 분석에서 "잔차의 등분산성"의 의미는 다음과 같다. "각 SBP값에 대해 잔차의 분산은 같다." 그림으로 설명하면 다음 그림에서, 주황 박스 안의 잔차들의 분산이 파란 박스 안의 잔차들의 분산과 같다는 것이다. 물론 박스는 저 두 곳에만 위치하는 것이 안이라 전 구간에 걸쳐 존재하게 된다. 

 따라서, 점들이 어떤 직사각형 안에 고루 퍼져있지 않고, 특정 Fitted value에서만 넓게 퍼져있거나 좁게 분포해 있다면 등분산성 가정에 위배되게 될 것이다.

 

 

만약 단순 선형 회귀 분석이 아니라 다중 선형 회귀 분석인 경우 "각 SBP(독립 변수) 값에 대해 잔차의 분산은 같다."가 아니라 "예측된 종속 변수(CVD_RISK) 값에 대해 잔차의 분산은 같다."가 될 것이다.

 

주관적이지 않나?

 위 분석들은 모두 주관적으로 결론 내리고 있는 것이 사실이다. 통계적으로 어떤 p-value를 내는 검정 방법이 충분히 존재하겠지만, 현재로서는 대부분은 주관적으로만 판단을 내리고 있는 것이 현실이다. 세상이 바뀌고 객관적인 결정을 내리기를 요구하는 시대가 도래할지도 모르지만, 만약 그런 세대가 온다면, 아마도 위 전제 조건들을 엄격하게 모두 만족하는 데이터는 그리 많지 않을 것이다. 아마 현재로서도 이러한 현실적인 이유 때문에 그렇게까지 가정을 검정하라고 하지는 않는 것 같다.

 하지만 제안된 통계적 방법은 존재한다. 바로 Edsel의 Global Validation of Linear Models Assumptions(GVLMA)다. 이는 gvlma패키지의 gvlma() 함수로 시행할 수 있다.

코드

#gvlma 패키지 설치
install.packages("gvlma")
library("gvlma")

#gvlma 시행
gv<-gvlma(LR_SBP_CVD)
summary(gv)

gv<-gvlma(LR_SBP_CVD) : 위에서 시행한 선형 회귀 분석 결과를 gvlma함수에 넣고 그 내용을 gv객체에 저장한다.
summary(gv) : gv 내용을 보여달라

 

결과

                    Value p-value                Decision
Global Stat        1.8578  0.7619 Assumptions acceptable.
Skewness           0.0394  0.8427 Assumptions acceptable.
Kurtosis           0.1602  0.6890 Assumptions acceptable.
Link Function      1.1687  0.2797 Assumptions acceptable.
Heteroscedasticity 0.4895  0.4841 Assumptions acceptable.

                                             Value p-value                                Decision

Global Stat                    1.8578  0.7619 Assumptions acceptable.

  : 아래 네 개를 종합적으로 판단한 결과다. p-value>0.05이므로 종합적으로 선형 회귀 분석의 가정이 성립한다고 본다.

Skewness                      0.0394  0.8427 Assumptions acceptable.

  : 정규성 판정 중 Skewness기준으로 본 것이다. 정규분포의 Skewness는 0이므로 0에 가까울수록 정규성을 따른다고 보며, p-value>0.05이므로 이 항목은 선형 회귀 분석의 가정이 성립한다고 본다.

Kurtosis                         0.1602  0.6890 Assumptions acceptable.

  : 정규성 판정 중 Kurtosis기준으로 본 것이다. 정규분포의 Kurtosis는 0이므로 0에 가까울수록 정규성을 따른다고 보며, p-value>0.05이므로 이 항목은 선형 회귀 분석의 가정이 성립한다고 본다.

Link Function              1.1687  0.2797 Assumptions acceptable.

  : link function을 제대로 설정했는가에 대한 항목으로 p-value>0.05이므로 이 항목은 선형 회귀 분석의 가정이 성립한다고 본다.

Heteroscedasticity   0.4895  0.4841 Assumptions acceptable.

  : 등분산성에 대한 항목이고 p-value>0.05이므로 이 항목은 선형 회귀 분석의 가정이 성립한다고 본다.

 

 

그런데, 왜 하필이면 이런 가정이 필요한걸까? 이는 선형회귀분석으로 구한 회귀계수들이 참값이기 위한 가정이며, 이 가정들을 가우스-마르코프 정리 (Gauss-Markov theorem)라고 한다. 관련 내용은 다음 링크에서 확인할 수 있다. 2023.06.21 - [통계 이론] - [이론] 가우스-마르코프 정리 (Gauss-Markov Theorem)

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")


##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

#회귀 분석 시행
LR_SBP_CVD<-lm(CVD_RISK~SBP, data=df)

#가정 확인 방법 1
plot(LR_SBP_CVD)

#가정 확인 방법 2
#gvlma 패키지 설치
install.packages("gvlma")
library("gvlma")

#gvlma 시행
gv<-gvlma(LR_SBP_CVD)
summary(gv)

[R] 선형 회귀 분석의 전제 조건 정복 완료!

작성일: 2022.12.22.

최종 수정일: 2023.10.17.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] [plot() 함수] 추세선 그리기 - plot(), abline()

 

앞선 포스팅에서 산점도를 그리는 방법을 소개하였다.(2022.12.07 - [[R] 그래프 작성] - [R] [plot() 함수] 산점도 그리기 (1) - plot(), colors())산점도를 그리고 난 뒤에는 추세선을 같이 그리기도 하는데, 추세선을 그리는 방법은 여러 가지가 있다. 여기에서는 선형 추세선을 그리는 법에 대해 소개하겠다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표: SBP와 CVD_RISK의 산점도를 그리고, 가장 잘 표현하는 추세선을 그려라

산점도 그리기- 코드

plot(df$SBP, df$CVD_RISK)

결과

 

이 둘은 어떤 선형 관계가 있는 것 같고 그 관계식을 구하는 법은 선형 회귀분석이다. 이는 다음 링크를 확인하길 바란다. 2022.12.22 - [선형 회귀 분석/R] - [R] 단순 선형 회귀 분석 (Simple linear regression) - lm()

우리는 다음 선형 회귀 분석을 통해 y절편이 23.036415, 기울기가 1.101935임을 알 수 있다.

선형 회귀 분석 코드

LR_SBP_CVD<-lm(CVD_RISK~SBP, data=df)
LR_SBP_CVD$coefficients

결과

(Intercept)         SBP 
  23.036415    1.101935

 

추세선 그리기

선형 추세선은 abline()이라는 함수를 이용한다. abline의 첫 번째 자리에는 y절편을, 두 번째 자리에는 x절편을 넣으면 된다. 

abline(23.036415, 1.101935)

결과

그런데, 색깔이 표본의 색과 검은색으로 같아 눈에 잘 띄지 않는다. 색깔은 파란색으로 바꿔주도록 한다. 색깔에 관한 내용은 다음 포스팅을 참고하길 바란다.2022.12.07 - [[R] 그래프 작성] - [R] [plot() 함수] 산점도 그리기 (1) - plot(), colors()

 

그리고, 너무 얇아 티가 나지 않으므로 "lwd"라는 옵션을 사용하여 선의 굵기 (Line WiDth)를 조정하겠다. 숫자를 바꿔가며 본인에게 맞는 굵기를 정하면 된다. 필자는 3으로 하겠다.

 

코드

abline(23.036415,1.101935, col="blue", lwd=3)

결과

 

 

추세선을 그리는 다른 방법

"23.036415","1.101935"라는 숫자를 적기가 너무 귀찮을 수도 있다. 이 값들은 각각 LR_SBP_CVD$coefficients의 첫 번째, 두 번째 값이므로 다음과 같이 표현할 수도 있다.

abline(LR_SBP_CVD$coefficients[1],LR_SBP_CVD$coefficients[2], col="blue", lwd=3)

혹은 그냥 단순 선형 회귀 분석을 저장한 객체를 넣는 것도 방법이 된다.

abline(LR_SBP_CVD, col="blue", lwd=3)

위 두개 코드는 같은 결과를 내준다.

 

[R] [plot() 함수] 추세선 그리기 정복 완료

작성일: 2022.12.22.

최종 수정일: 2022.12.23.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

 

[R] 단순 선형 회귀 분석 (Simple linear regression) - lm()

 

 지난 상관관계 분석에서 우리는 서로 다른 두 변수가 얼마나 관련되어 있는지를 알아보았다. 즉, 한 변수의 값으로 다른 변수의 값을 예측할 수 있는 정도인 "예측도"의 관점에서 상관 계수를 이해하였다.(2022.12.15 - [상관분석/R] - [R] 피어슨 상관 계수 (Pearson's correlation coefficient) - cor.test()) 그런데, 상관 계수의 단점은 "예측도"만 확인할 수 있을 뿐, 실제로 예측을 할 수는 없다. 가령, 수축기 혈압(SBP)이 120인 사람의 심혈관 질환 위험 점수 (CVD_RISK)가 얼마인지는 알 수 없는 것이다. 이걸 가능하게 하는 것이 선형 회귀 분석 (linear regression)이다. 본 포스팅에서는 독립 변수가 1개인 단순 선형 회귀 분석 (simple linear regression)을 시행해보도록 하겠다. 이는 독립 변수가 1개라는 점에서 univariate linear regression이라고도 부른다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.18)

 

분석용 데이터 (update 22.12.18)

2022년 12월 18일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 상관분석 - 반복 측정 자료

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  수축기 혈압(SBP)으로 심혈관 질환 위험 점수 (CVD_RISK)를 예측할 수 있는가?

 

선형 회귀 분석은 lm()이라는 함수를 사용한다.

코드

LR_SBP_CVD<-lm(CVD_RISK~SBP, data=df)
summary(LR_SBP_CVD)

LR_SBP_CVD<-lm(CVD_RISK~SBP, data=df) : 데이터는 df를 사용한다. 독립변수를 SBP로, 종속변수를 CVD_RISK로 하는 선형 회귀 분석을 시행하고 이 결과를 LR_SBP_CVD에 저장하라.
summary(LR_SBP_CVD) : LR_SBP_CVD의 내용을 요약해서 보여달라.

 

결과

Call:
lm(formula = CVD_RISK ~ SBP, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.6754  -3.6753   0.0733   3.4071  15.9155 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 23.036415   1.233590   18.67   <2e-16 ***
SBP          1.101935   0.009055  121.69   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.176 on 998 degrees of freedom
Multiple R-squared:  0.9369,	Adjusted R-squared:  0.9368 
F-statistic: 1.481e+04 on 1 and 998 DF,  p-value: < 2.2e-16

하나씩 살펴보도록 하자.

 

Call: lm(formula = CVD_RISK ~ SBP, data = df)

: 당신이 요청한 분석은, 데이터는 df를, 종속변수는 CVD_RISK를, 독립변수는 SBP를 사용하는 선형 회귀 분석이다.

 

Residuals:      Min            1Q   Median         3Q        Max 

               -16.6754  -3.6753   0.0733   3.4071  15.9155 

: 잔차의 최솟값은 -16.6754, 중앙값은 0.0733, 최댓값은 15.9155다.

잔차는 실제 관측값에서 예측값을 뺀 것이고, 예측값을 계산하는 방법은 아래에 나올 것이다.

 

Coefficients:

                         Estimate   Std. Error  t value    Pr(>|t|)    

(Intercept) 23.036415   1.233590    18.67   <2e-16 ***

SBP                 1.101935  0.009055  121.69   <2e-16 ***

--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

SBP에 관한 것부터 보자.

 추정치(Estimate)는 1.101935다. 이는 SBP가 1 증가할 때, CVD_RISK가 1.101935 증가한다는 의미다. 즉, 이는 $xy$평면에서 기울기를 의미하고, 양의 상관관계가 있다는 말이다. 이 값은 논문에서 종종 $\beta$값으로 불린다.

 추정치의 표준 오차는 0.009055다. 이는 추정치(Estimate)의 변동가능성을 의미하고, 이보다 오른쪽에 있는 Pr(>|t|)를 구할 때 사용된다.

Pr(>|t|)는 2e-16보다 작다. 이는 p-value를 의미하는데, $2\times 10^{-16}$보다 작으므로 매우 유의한 결과임을 알 수 있다. 그러므로 귀무가설을 기각하면 되는데, 선형 회귀분석에서는 귀무가설과 대립가설이 무엇일까?

 귀무가설: $H_0=$ 추정치(Estimate)는 0이다. 즉, 이 표본이 기원한 모집단에서 SBP와 CVD_RISK는 아무런 관련성이 없다.

 대립가설: $H_1=$ 추정치(Estimate)는 0이 아니다. 즉, 이 표본이 기원한 모집단에서 SBP와 CVD_RISK는 모종의 관련성이 있다.

따라서, 우리는  SBP와 CVD_RISK 사이의 양의 상관관계가 유의미하며, 추정치 1.101935는 유의미하다고 결론내린다.

 

하지만, 우리가 선형 회귀 분석을 시작한 목적은 SBP의 값으로 CVD_RISK의 값을 예측하기 위함이었는데, 기울기만 알아서는 $xy$평면에서 직선을 그릴 수 없다. 우리는 $y$절편 또한 알아야 한다. 따라서 Intercept값을 보아야 한다.

Intercept는 23.036415이다. 사실 이 값에 대한 p-value는 큰 의미가 없다. 왜냐하면 우리는 이 값이 0이든 아니든 관심이 없기 때문이다. 

 

어찌 되었든 우리는 SBP와 CVD_RISK사이의 관계식이 다음과 같음을 알 수 있게 되었다.

$CVD\_RISK=23.036415+1.101935\times SBP$

 

(Residual standard error부터의 내용은 조금 뒤에 다시 설명하겠다.)

 

그런데 SBP와 CVD_RISK의 산점도를 그리고 위 식을 같이 그려 넣으면 다음과 같은 plot을 얻을 수 있다.

산점도 그리는 법:2022.12.07 - [[R] 그래프 작성] - [R] [plot() 함수] 산점도 그리기 (1) - plot(), colors()

추세선 그리는 법:2022.12.22 - [[R] 그래프 작성] - [R] [plot() 함수] 추세선 그리기 - plot(), abline()

#산점도 그리기
plot(df$SBP, df$CVD_RISK)

#추세선 그리기 방법 1
abline(23.036415,1.101935, col="blue", lwd=3)

#추세선 그리기 방법 2
abline(LR_SBP_CVD$coefficients[1],LR_SBP_CVD$coefficients[2], col="blue", lwd=3)

#추세선 그리기 방법 3
abline(LR_SBP_CVD, col="blue", lwd=3)

추세선 그리기 방법 2는 다음과 같이 이해할 수 있다.

LR_SBP_CVD$coefficients는 위의 summary(LR_SBP_CVD)중 다음 부분만을 불러온 것이다.

Coefficients:

                         Estimate

(Intercept)   23.036415

SBP                 1.101935 

따라서 LR_SBP_CVD$coefficients의 첫 번째 항목을 의미하는 LR_SBP_CVD$coefficients[1]은 Intercept을 의미한다.

LR_SBP_CVD$coefficients의 첫 번째 항목을 의미하는 LR_SBP_CVD$coefficients[2]는 기울기를 의미한다. 

 

결과

추세선처럼 데이터에서 SBP가 증가할수록 CVD_RISK가 증가하고 있다. 하지만 모든 데이터가 완벽히 선 위에 있는 것은 아니다. 선 위의 값들은 예측된 값이므로 실제값과는 차이가 있다. 이 차이를 Residual라고 한다. 

 

그럼 원래 summary(LR_SBP_CVD)의 내용으로 돌아가보자.

 

Residual standard error: 5.176 on 998 degrees of freedom :

이는 위에서 설명한 residual의 표준오차가 5.176이라는 뜻이다. 자유도(degrees of freedom)가 언급된 이유는, 표준오차를 계산할 때 자유도를 이용하기 때문이다.

Multiple R-squared:  0.9369 : R squared는 SBP가 CVD_RISK를 어느 정도 설명하는지를 의미한다. 여기에서는 SBP가 CVD_RISK를 93.69% 설명한다고 말할 수 있다.

Adjusted R-squared:  0.9368 : 본 분석은 독립 변수가 1개이므로 별로 의미는 없는 지표이다. 그럼에도 불구하고 Adjusted R-squared를 설명하면, 다음과 같다. 독립 변수가 많아지면 많아질수록 종속 변수를 설명하는 능력은 증가되기 마련이다. 하지만 마냥 변수가 많아진다고 실제로 설명력이 높아지는 것은 아니다. 따라서 어느 정도의 페널티를 주어야 하는데, 이 패널티를 고려한 것이 Adjusted R-squared이다.

F-statistic: 1.481e+04 on 1 and 998 DF,  p-value: < 2.2e-16 : 이 내용도 독립 변수가 1개인 본 분석에는 별로 의미가 없다. 그럼에도 불구하고 설명하면 이 p-value의 귀무가설은 "모든 독립 변수의 Estimate가 0이다"이다. p-value가 0.05보다 작으므로 적어도 한 개 이상의 Estimate가 유의미한 값을 갖는다는 결론을 내리게 된다. 하지만 여기에서는 독립 변수가 원래 1개이므로 그 독립변수, SBP가 유의미하다는 것을 알 수 있다.

 

 

가정

 단순 선형 회귀 분석 자체는 간단하지만, 단순 선형 회귀 분석 시행의 전제조건을 따지는 것이 이보다 까다롭고 더 중요하다. 이는 내용이 꽤 많으므로 다음 포스팅을 참고하길 바란다.2022.12.22 - [선형 회귀 분석/R] - [R] 선형 회귀 분석의 전제 조건 - lm(), plot()

  

 

독립 변수의 종류: 연속형

 필자는 연속형 변수인 SBP를 독립 변수로 사용하였다. 만약 범주형 변수를 사용하고자 한다면 분석 방법이 약간은 달라지게 된다. 이에 관한 내용은 다음 링크에서 확인할 수 있다.2022.12.22 - [선형 회귀 분석/R] - [R] 범주형 변수의 선형 회귀 분석 (Simple linear regression with categorical variables) - lm(), factor()

 

[R] 단순 선형 회귀 분석 (Simple linear regression) 정복 완료!

 

작성일: 2022.12.22.

최종 수정일: 2023.05.14.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형

+ Recent posts