[R] 편상관 계수, 부분 상관 계수 (Partial correlation coefficient) - pcor.test()

[R] 편상관 계수, 부분 상관 계수 (Partial correlation coefficient) - pcor.test()

 

세상의 많은 지표들은 서로 상관되어있는 경우가 많다. 비흡연자에 비해 흡연자가 폐암에 잘 걸리고, 부모의 키가 클수록 아이의 키도 크다. 소득 수준이 높은 나라일수록 평균 수명도 길고, 수학 점수가 높은 학생은 물리 점수가 높으며, 식당에서 음식을 많이 시킬수록 계산해야 하는 금액도 커진다. 변수의 형태에 따라 피어슨 상관 계수 (Pearson's correlation coefficient)를 구하거나 스피어만 상관 계수 (Spearman's correlation coefficient), 혹은 켄달의 타우 (Kendall's tau)를 구하기도 한다. 관련 내용은 다음 링크에서 확인할 수 있다.

2022.12.15 - [상관분석/R] - [R] 피어슨 상관 계수 (Pearson's correlation coefficient) - cor.test()

2022.12.16 - [상관분석/R] - [R] 스피어만 상관 계수, 켄달의 타우 (Spearman's rank correlation coefficient, Kendall's tau) - cor.test()

 

 그런데 만약 다음과 같은 상황이라면 어떨까?

목표: 수축기 혈압이 심혈관 질환 위험에 미치는 영향을 알아보고자 한다. 

상황: 수축기 혈압이 높은 사람은 일반적으로 BMI도 높다. 그런데, 높은  BMI는 심혈관 질환의 위험 인자다.

문제: 높아진 심혈관 질환 위험이 수축기 혈압 때문인지, BMI 때문인지 구분할 수가 없다.

 

 이런 상황에서 BMI를 교란 요인 (confounding factor, confounding variable, confounder)이라고 하며, 교란요인의 영향을 보정하는 것은 의학 통계 전반에 있어서 아주 중요한 문제다. 

 

 상관 계수 (correlation coefficient)를 구할 때 교란 변수를 보정하는 방법은 편상관 계수 혹은 부분 상관 계수(partial correlation coefficient)를 구하는 것이다. 이에 대해서 알아보도록 하겠다.

 

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

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

분석용 데이터 (update 22.12.18)

 

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

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

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의 영향을 배제 혹은 교정하고도 수축기 혈압은 심혈관 질환 위험과 상관 관계를 갖는가?

 

편상관 계수를 구하는 함수는 pcor.test()인데, 이는 ppcor이라는 패키지에 있으므로 먼저 설치를 한다. 패키지 설치에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()

install.packages("ppcor")
library("ppcor")

 

코드 - 일반적인 피어슨 상관 계수 구하기

먼저 수축기 혈압과 심혈관 질환 위험 간의 상관 계수를 구해보자. 이는 이전 포스팅에서 다룬 cor.test() 함수를 사용하면 된다. 2022.12.15 - [상관분석/R] - [R] 피어슨 상관 계수 (Pearson's correlation coefficient) - cor.test()

cor.test(df$SBP, df$CVD_RISK)

결과

	Pearson's product-moment correlation

data:  df$SBP and df$CVD_RISK
t = 121.69, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9637531 0.9716095
sample estimates:
     cor 
0.967917

상관 계수가 0.967로 매우 높다. 

 

그렇다면, BMI의 영향을 배제하는 코드는 어떻게 될까?

코드 - 편상관 계수 / 부분 상관 계수

pcor.test(df$SBP, df$CVD_RISK, df$BMI)

pcor.test(df$SBP, df$CVD_RISK, df$BMI) : df 데이터의 SBP와 df 데이터의 CVD_RISK간의 상관계수를 구하라. 근데, df 데이터의 BMI의 영향은 배제시킨 편상관 계수 / 부분 상관 계수를 구하라.

 

결과

   estimate       p.value statistic    n gp  Method
1 0.8524241 4.948042e-283  51.47844 1000  1 pearson

       estimate                  p.value     statistic       n gp  Method

1 0.8524241  4.948042e-283  51.47844 1000  1 pearson

 

BMI로 보정하면, SBP와 CVD_RISK사이의 상관 계수는 0.8524241로 보정하기 전의 0.968보다는 소폭 감소한 모양이다. BMI가 CVD_RISK를 증가시키는 설명력을 배제하였기 때문이다. 이때의 p-value는 4.948042e-283로 매우 작으므로 0이 아니라는 대립 가설을 택할 수 있고, 구해진 상관 계수는 Pearson 상관 계수다.

 

두 개 이상의 교란 변수를 보정

만약, 두 개 이상의 교란 변수를 보정하고자 한다면 보정할 변수들을 list로 만들어 넣으면 된다. BMI뿐만 아니라 흡연력인 SMOK도 보정변수로 넣고 싶다면 다음과 같은 코드를 쓰면 된다. (SMOK는 순서형 변수이므로 Pearson이 적절하지 않으므로 Spearman이나 Kendall 사이에서 골라야 하는데, 같은 값을 갖는 사람(tie)이 많으므로 Kendall의 tau를 구하도록 한다.)

pcor.test(df$SBP, df$CVD_RISK, list(df$BMI,df$SMOK), method="kendall")

 

결과

   estimate       p.value statistic    n gp  Method
1 0.6301046 3.236403e-195  29.80631 1000  2 kendall

 

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

상관계수가 "얼마나 높아야 높다고 할 수 있는가?"라는 의문이 들 수 있다. 즉, "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

 

[R] 편상관 계수, 부분 상관 계수 (Partial correlation coefficient) 정복 완료!

작성일: 2022.12.18.

최종 수정일: 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