상관관계 분석은 두 변수 간의 관계를 이해하는 가장 기본적인 통계 기법 중 하나입니다. R에서 cor.test()를 사용한 기본적인 분석은 익숙할 것입니다.

result <- cor.test(dataset$var1, dataset$var2, method = "pearson")
result

위 코드는 두 변수 간의 피어슨 상관계수를 계산하고 결과를 출력합니다. 하지만 탐색적 데이터 분석(EDA) 단계에서는 어떤 변수들이 서로 관련 있는지 사전에 알기 어렵습니다. 수백 개의 변수를 가진 데이터셋에서 특정 변수와 관련된 다른 변수들을 빠르게 찾아야 한다면 어떨까요? rstatix 패키지의 cor_test() 함수는 이러한 상황에서 강력한 도구가 됩니다.

데이터 준비

예시로 세계은행(World Bank)의 세계개발지표(WDI) 데이터를 사용하겠습니다. WDI는 전 세계 개발 관련 지표를 제공하는 오픈 데이터베이스입니다. WDI R 패키지를 통해 데이터를 불러올 수 있습니다.

install.packages("WDI")
library(WDI)

# 전체 데이터베이스 다운로드
raw_data <- WDIbulk(timeout = 600)

우리의 관심사는 "경제 규모 대비 무역량이 많은 국가들의 특징"이며, 데이터는 2020년 기준으로 분석합니다. 먼저 데이터를 정리합니다.

# 연간 데이터만 필터링
annual_series <- raw_data$Series %>% filter(Periodicity == "Annual")

# 2020년 데이터에서 결측치가 100개 미만인 변수는 제외
valid_vars <- raw_data$Data %>%
  filter(Indicator.Code %in% annual_series$Series.Code,
         year == 2020,
         !is.na(value)) %>%
  group_by(Indicator.Code) %>%
  summarise(n = n()) %>%
  filter(n > 100)

분석 실행

이제 약 790개의 변수와 무역량(GDP 대비 무역 비율, 코드: NE.TRD.GNFS.ZS) 간의 상관관계를 한 번에 분석합니다. cor_test()는 tidyverse 스타일의 파이프 연산을 지원하며, 결과를 tibble 형태로 반환합니다.

# 국가만 추출 (지역 데이터 제외)
countries <- raw_data$Country %>% filter(!Region == "") %>% as_tibble()

# 상관분석 실행
cor_results <- raw_data$Data %>%
  filter(Indicator.Code %in% valid_vars$Indicator.Code,
         year == 2020,
         Country.Code %in% countries$Country.Code) %>%
  select(-Indicator.Name) %>%
  pivot_wider(names_from = Indicator.Code, values_from = value) %>%
  cor_test(NE.TRD.GNFS.ZS, method = "pearson", use = "complete.obs")

cor_results

결과는 변수 쌍별로 상관계수(r), t-통계량, p-값, 신뢰구간 하한/상한을 포함한 깔끔한 테이블입니다. p < 0.05 수준에서 유의한 변수 중 상관계수가 높은 순서로 정렬하여 살펴봅시다.

# 변수명 매핑 테이블 생성
var_labels <- raw_data$Series %>%
  select(Series.Code, Indicator.Name, Short.definition) %>%
  as_tibble()

# 유의한 결과만 추출하여 상관계수 내림차순 정렬
cor_results %>%
  left_join(var_labels, by = c("var2" = "Series.Code")) %>%
  filter(p < 0.05) %>%
  arrange(desc(cor)) %>%
  View()

예상대로 서비스 무역, 상품 무역과의 상관성이 높았지만, 의외로 공적개발원조(ODA) 수혜액과도 중간 정도의 양의 상관관계(r = 0.43)를 보였습니다.

그룹별 분석

이 관계를 더 깊이 파고들어 봅시다. 예를 들어, 시간이 지남에 따라 이 상관성이 어떻게 변하는지 확인하려면 어떻게 해야 할까요? cor_test()는 group_by()와 함께 사용할 수 있어 매우 편리합니다.

# 무역량과 ODA 지표만 추출하여 연도별 상관분석
trend_analysis <- raw_data$Data %>%
  filter(Indicator.Code %in% c("NE.TRD.GNFS.ZS", "DT.ODA.ODAT.GI.ZS"),
         Country.Code %in% countries$Country.Code,
         year < 2021) %>%
  select(-Indicator.Name) %>%
  pivot_wider(names_from = Indicator.Code, values_from = value) %>%
  group_by(year) %>%
  cor_test(NE.TRD.GNFS.ZS, DT.ODA.ODAT.GI.ZS,
           method = "pearson", use = "complete.obs")

# 결과 시각화
library(ggplot2)
trend_analysis %>%
  mutate(significant = ifelse(p < 0.05, "Yes", "No")) %>%
  ggplot(aes(x = year, y = cor)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line() +
  geom_point(aes(color = significant)) +
  labs(y = "상관계수 (r)") +
  theme_minimal()

시각화 결과, 과거에는 두 변수 간 상관성이 거의 없었지만, 최근 몇 년간 유의한 양의 상관관계가 나타나는 것을 확인할 수 있습니다. 이는 추가 연구를 위한 흥미로운 가설을 제시합니다. 상관관계가 인과관계를 의미하지는 않지만, 데이터 기반 가설 생성에 유용합니다.

rstatix 패키지의 장단점

rstatix의 cor_test()는 탐색적 분석에 매우 효율적이지만, 몇 가지 제한점이 있습니다.

• 장점: tidyverse 스타일의 파이프 연산 지원, 빠른 실행 속도, tibble 형태의 깔끔한 출력
• 단점: 피어슨, 스피어만, 켄달 상관계수만 지원 (다양한 방법론 부족), 신뢰구간은 피어슨 r에만 제공, 표본 크기나 자유도는 출력되지 않음

그럼에도 불구하고, 초기 탐색 단계에서는 이러한 한계가 크게 문제되지 않습니다. rstatix는 cor_test() 외에도 다양한 통계 검정 함수를 제공하므로, EDA 작업 시 유용하게 활용할 수 있습니다.