[Multilevel Regression] 다층 회귀 분석 in R (1) 개념, 데이터 전처리, lme4 패키지, Null Model 돌리기

작년 말 나를 괴롭게 했던 다층회귀분석.

석사까지는 어떻게 통계 개념 허술해도 대충 아는 척하고 넘어갈 수 있었는데, 박사 과정에 오고나니 허술하면 들통날 수 밖에 없다. 교수님들과 계속 내 분석을 공유해야하고 그 과정에서 대차게 까일 수 있으니 말이다.. ㅠㅠ 물박사가 될 바엔 학생일 때 좀 괴로운 게 낫다고 생각해서 꾸역꾸역 공부했다. 하고싶은 것만 할 수는 없지!

 

그리고 통계를 제대로 모르면 사실 딥러닝 모델을 이해하고 개선시키는 데에도 큰 한계가 있다는 걸 느꼈다. 요즘 모델 돌리기는 너무 쉽지만, '좋은' 모델을 만드는 것은 여전히 어렵다. 결국 머신러닝 모델을 잘 만들기 위해서는 어떤 요소가 모델링에 중요한가를 알아야하는데, 이걸 알기 위해 정통 통계만큼 정확한 것이 없다. 

 

Jumping right into Machine Learning severely loses explainability. Doable, but at what cost?  - me

 

🤔 Multilevel Regression을 왜 해야할까?

일반적으로 쓰는 회귀분석은 각 데이터 포인트가 independent 하다는 커다란 조건 (i.i.d.) 을 가지고 시작한다. 도시데이터로 분석을 하는 나의 분야에서 실제로 데이터가 모두 independent하지 않은 경우가 종종 있다. 특히, 설문 데이터를 많이 이용하는 경우 설문이 누구를 대상으로 어떻게 뿌려졌는지를 설명하는 Sampling Design 을 자세히 보아야한다. 국가에서 진행하는 (National) 설문조사의 경우, 특정 인구에게 무작위에게 뿌리는 경우는 생각보다 별로 없고, 보통 주(State) 단위로 먼저 나누고 인구 수에 맞춰 응답자 수를 맞추는 경우가 많다.

 

내가 이번에 쓰게 된 RECS (Residential Energy Consumption Survey) 설문 역시 미국을 52개주로 먼저 나눈 후, 각 주에서 목표한 인구를 무작위로 추출했다. 주 내에서는 무작위이지만, 주 외에서는 주마다 sampling weight가 달랐다. 

이와 같이, 데이터에 내재된 다층적인 특성에 따라 일반적인 회귀분석을 쓰지 못할 때, 다층회귀를 써야한다. 

다중회귀는 Heterogeneity Issue (설명변수와 종속변수와의 관계가 클러스터마다 일정하지 않을 때 생기는 문제),  그리고 Robinson effect (aggregation bias  -하나의 클러스터에만 해당하는 것을 다른 클러스터에 일반화하는 것) 를 방지할 수 있다.

 

🤔 Multilevel Regression의 세가지 기본 개념

다층회귀에서는 아래 1,2번 관계로 나누어서 회귀식을 작성하게 되며, 3번 ICC를 통해서 다층회귀가 적절한 것인지 판단할 수 있다.

1) Fixed Effects - 모든 그룹에게 동일한 설명-종속변수 관계. 여기서 쓰는 RECS 데이터의 경우, state와 상관없이 독립변수가 종속변수에 동일하게 주는 영향

2) Random Effects - 그룹마다 다른 설명-종속변수 관계. state에 따라서 다른, 독립변수가 종속변수에게 주는 영향.

3) Intraclass Correlation Coefficient (ICC) - (그룹 간 variance)/(전체 variance), 즉, 전체 분산 중에 그룹별로 다른 관계가 나타나는 것으로 인해 발생하는 분산은 어느정도인지 측정한 값.

> Rule of Thumb! 주로 ICC가 10% 이상일때, 다층회귀를 쓴다. ICC를 검토한 후에, 4단계를 거쳐 최종 모델 식을 도출한다. 

 

 

👩‍💻 lme4 패키지로 다층회귀 시작하기 

이 예시에서는 앞서 말한 RECS 서베이 데이터를 이용해서, 사람들이 집에서 사용하는 에너지 총량 (TOTAL BTU) 값이 어떤 설명변수들에 의해 달라지는지 다층회귀 분석을 해보려고 한다. 이 설문데이터는 미국의 52개주 단위로 나누어 구득이 되었기 때문에 주에 따라 다른 Random Effect, 주와 상관없이 다른 Fixed Effect를 측정할 수 있다. 

 

Step 1. 라이브러리 설치

우선 필요한 라이브러리를 설치, import한다.

# Install packages
install.packages("lme4")
install.packages("lmerTest")

# Multilevel regression
library(lme4) 
library(lmerTest) # to get p-value from t-values in lmer summary

# Data Wrangling
library(dplyr)
library(ggplot2)
library(tidyverse)

 

 

 

Step 2. 데이터 전처리 📊 

사실 가장 많이 시간을 쓰고, 전 과정을 여러번 엎어야했던 이유가 데이터 전처리 과정에서 놓친 점들 때문이었다.

아래 세가지 스텝은 꼭 거쳐야 한다. 

2.1 Outlier 제거하기 -- 가장 먼저, 데이터셋에서 말도 안되는 값을 가진 row들을 제거한다. 

# (1) Delete Outlier
#   For each State (lowest strata level), delete TOTALBTU outliers
#   Outliers can be non-residential (e.g., home businesses), which is not REC
recs_mr2 <- recs_mr %>% 
  group_by(STATE_FIPS) %>%
  mutate(
    Q1 = quantile(TOTALBTU, 0.25),
    Q3 = quantile(TOTALBTU, 0.75),
    IQR = Q3 - Q1,
    lower_bound = Q1 - 1.5 * IQR,
    upper_bound = Q3 + 1.5 * IQR
  ) %>%
  filter(TOTALBTU >= lower_bound & TOTALBTU <= upper_bound) %>%
  dplyr::select(-Q1, -Q3, -IQR, -lower_bound, -upper_bound) # Remove additional columns

dim(recs_mr) # 18496  
dim(recs_mr2) # 18003 -- 493 responses were deleted as they are outliers

 

2.2 Data Type 지정하기 -- Categorical, Dummy (0,1), Ordinal 데이터 처리

가진 데이터의 metadata를 잘 참고해서 NA값이 있거나 missing data 값 (0, -1, -2, 99 와같은) 이 있다면 지우거나, 따로 처리해준다. 아래는 R 예시이다. 

# Categorical Variables 
recs_mr2$TYPEHUQ <- factor(recs_mr2$TYPEHUQ,  levels = c(1,2,3,4,5))
recs_mr2$UATYP10 <- factor(recs_mr2$UATYP10,  levels = c("C", "R", "U"))

# Dummy variables
recs_mr2$CRAWL[recs_mr2$CRAWL == -2] <- 0
recs_mr2$HEATAPT[recs_mr2$HEATAPT == -2] <- 0
recs_mr2$COOLAPT[recs_mr2$COOLAPT == -2] <- 0
recs_mr2$CRAWL <- factor(recs_mr2$CRAWL,  levels = c(1,0))
recs_mr2$HEATHOME <- factor(recs_mr2$HEATHOME,  levels = c(1,0))
recs_mr2$HEATAPT <- factor(recs_mr2$HEATAPT,  levels = c(1,0))
recs_mr2$AIRCOND <- factor(recs_mr2$AIRCOND,  levels = c(1,0))
recs_mr2$COOLAPT <- factor(recs_mr2$COOLAPT,  levels = c(1,0))

# Ordinal Variables
recs_mr2$DRAFTY <- factor(recs_mr2$DRAFTY, ordered = TRUE, levels = c(1:4))
recs_mr2$MONEYPY[recs_mr2$MONEYPY < 9] <- 1
recs_mr2$MONEYPY[recs_mr2$MONEYPY %in% c(9:11)] <- 2
recs_mr2$MONEYPY[recs_mr2$MONEYPY %in% c(12:13)] <- 3
recs_mr2$MONEYPY[recs_mr2$MONEYPY > 13] <- 4
recs_mr2$MONEYPY <- factor(recs_mr2$MONEYPY, ordered = TRUE, levels = c(1:4))
recs_mr2$EDUCATION <- factor(recs_mr2$EDUCATION, ordered = F, levels = c(1:5))

recs_mr2$EQUIPAGE[recs_mr2$EQUIPAGE == -2] <- 0
recs_mr2$ACEQUIPAGE[recs_mr2$ACEQUIPAGE == -2] <- 0
recs_mr2$ACEQUIPAGE[recs_mr2$ACEQUIPM_PUB == -2] <- 0
recs_mr2$EQUIPAGE <- factor(recs_mr2$EQUIPAGE, ordered = TRUE, levels = c(0:6))
recs_mr2$ACEQUIPAGE <- factor(recs_mr2$ACEQUIPAGE, ordered = TRUE, levels = c(0:6))

 

2.3 Continuous data는 Group Mean으로 Standardize 

모든 연속 변수는 클러스터 설명변수 레벨에서 평균값을 중심으로 standardize 한다. 이 예시에서는 주 레벨로 데이터가 구득되었으므로, 주 별 값 평균을 구한 뒤 standardize해준다. 

# (3) Centering Continuous Variables
#   For each state, center continuous variables
#   Group-mean centering

recs_mr3 <- recs_mr2 %>%
  group_by(STATE_FIPS) %>%
  mutate(HDD65_mean = mean(HDD65),
         CDD65_mean = mean(CDD65),
         #MONEYPY_mean = mean(MONEYPY),
         YEARMADERANGE_mean = mean(YEARMADERANGE),
         SQFTEST_mean = mean(SQFTEST),
         NHSLDMEM_mean = mean(NHSLDMEM),
         NUMCHILD_mean = mean(NUMCHILD),
         NUMADULT1_mean = mean(NUMADULT1), 
         NUMADULT2_mean = mean(NUMADULT2)
         ) %>% 
  mutate(CENHDD65 = HDD65 - HDD65_mean,
         CENCDD65 = CDD65 - CDD65_mean,
         # CENMONEYPY = MONEYPY-MONEYPY_mean,
         CENYEARMADERANGE = YEARMADERANGE - YEARMADERANGE_mean,
         CENSQFTEST = SQFTEST - SQFTEST_mean,
         CENNHSLDMEM = NHSLDMEM - NHSLDMEM_mean,
         CENNUMCHILD = NUMCHILD - NUMCHILD_mean,
         CENNUMADULT1 = NUMADULT1 - NUMADULT1_mean,
         CENNUMADULT2 = NUMADULT2 - NUMADULT2_mean
         ) %>%
  ungroup()

 

+ ))) 2.4 Sampling Weight 을 다루는 방식은 사실 다층회귀의 개념보다는 설문데이터의 속성에 의해 필요한 것이기 때문에 여기서 다루지 않겠다. 

 

Step 3. R로 나의 데이터가 다층회귀에 적합한지 확인하기.

Null Model - (클러스터를 만드는 설명변수 - 종속변수)만 넣고 회귀모델 돌리기

첫 스텝은 ICC를 이용해 다층회귀가 적절한 모델인지 확인하는 단계이다. 적절한지는  모델을 돌린 뒤, 모델의 ICC = (Between Group Variance)/(Within Group Variance) 가 0.1 이상이면 적절! 

 

이 예시에서는 주 (State) 변수가 내 데이터에서 클러스터 (그룹)을 결정하므로  이 경우 State Code 변수가 '클러스터를 만드는 설명변수 Group Variable' 이다.  그리고 종속변수는 총 에너지 소비량인 TOTALBTU 변수이다. 따라서 lmer 함수를 이용해 아래와 같이 코드를 돌린다. 

Y: TOTALBTU

X: STATE_FIPS

 

일반적인 회귀식과 같이 Y~X 형식으로 formula를 적으면 되는데, 클러스터 설명변수의 경우, (1| X) 괄호 안에 넣어준다. 

클러스터 설명변수가 2개 이상인 경우에도 (1 | X1 + X2.. ) 이런 식으로, 괄호 안에 적어준다.  모델이 (1 앞의 것은 일반적 설명변수, 뒤의 것은 클러스터 설명변수로 인식한다. 

nullm_a <- lmer(TOTALBTU ~ 1 + (1|STATE_FIPS), data = recs_rescaled, weights = pweights_a, REML=FALSE)
nullm_b <- lmer(TOTALBTU ~ 1 + (1|STATE_FIPS), data = recs_rescaled, weights = pweights_b, REML=FALSE)

 

모델의 적절함을 보기 위해서, ICC계산 및 여러가지 모델 goodness of fit을 측정하는 함수를 만든다. 

measureF <- function(model){
  
  # Check residuals
  #plot(resid(model))
  #qqnorm(resid(model))
  
  # Extracting variance components
  variances <- VarCorr(model) 
  variance_df <- as.data.frame(variances)
  
  # Between-group variance
  between_group_var <- as.numeric(variance_df$vcov[1])
  cat("Between group variance:", between_group_var,"\n")
  
  # Residual (within-group) variance
  residual_var <- attr(variances, "sc")^2
  cat("Withint group variance:", residual_var,"\n")
  
  # Calculating ICC
  icc <- between_group_var / (between_group_var + residual_var)
  cat("ICC:", icc,"\n")
  
  # Extracting AIC and BIC
  aic_value <- AIC(model)
  bic_value <- BIC(model)
  
  # Displaying AIC and BIC
  cat("AIC: ", aic_value, "\n")
  cat("BIC: ", bic_value, "\n")
  
  # Extracting log-likelihood and number of parameters
  logLik_value <- logLik(model)
  k <- attr(logLik_value, "df")
  n <- length(resid(model))
  
  # Calculating AICC
  aicc_value <- -2 * logLik_value + (2 * k * (k + 1))/(n - k - 1)
  
  # Displaying AICC
  cat("AICC: ", aicc_value, "\n")
  
  # Calculating CAIC
  caic_value <- -2 * logLik_value + k * (log(n) + 1)
  
  # Displaying CAIC
  cat("CAIC: ", caic_value, "\n")
  
  # Calculating HQIC
  hqic_value <- -2 * logLik_value + 2 * k * log(log(n))
  
  # Displaying HQIC
  cat("HQIC: ", hqic_value, "\n")
}

 

그 후, measureF함수에 모델을 인풋으로 넣고 돌리면 variance를 비롯해서 ICC, AIC, BIC등의 측정값을 볼 수 있다. 아래처럼 ICC가 0.15, 즉 15% 이면 전체 데이터 분산 중 15%가 클러스터 간의 variance로 인해 설명 가능하다고 해석한다. 

데이터의 속성에 따라서 weight를 두가지를 테스트해봐야할 수 있는데 이 페이퍼를 참고해서 계산했다. 두가지 weight로 모두 돌려본 후 결과가 비슷한 것을 확인해야한다. 

 

 

이번 포스팅에서는 데이터 전처리와 나의 데이터가 다층회귀에 적합한지 알아보는 단계를 거쳤다. 다음 포스팅에서는 Level 2, Level 1 variable을 추가하는 방법, 그리고 최종모델을 도출한 후 모델의 퍼포먼스를 측정하는 단계까지 설명해보겠다. 

 

< 다음 포스팅 개념 >

Step 4. Level 2 Variables - 클러스터 간 다른 설명변수 추가하기

Level 2 variable은 클러스터 내에서는 같지만, 클러스터 간에는 다른 변수를 말한다. 혹시 클러스터에 따라 다른 영향을 주는 변수가 있는지 확인하고, 따로 처리해서 일반회귀의 iid 조건을 맞추기 위해 level 2 variable을 걸러내야 한다. 예시에서 같은 주 내에서는 값이 같지만, 주 간에는 다른 변수들이 있다. 예를 들어, 주 평균 인당 가처분 소득 (state-avg PDPI), 주 평균 에너지 소비자 가격이 그렇다. 주마다 다르지만, 주내에서는 같은 평균값을 가진다. 이 변수들을 하나씩 추가해보고, AIC, BIC가 10 이상 떨어진다면 그 변수를 level 2 variable로 채택한다. (*AIC, BIC는 작을수록 좋은 값이다.) 

 

Step 5. Level 1 Variables - 나머지 변수 추가하기

클러스터와 상관없이 모두 값이 다른 변수들이 있을 것이다. 그 변수를 level 1 variable  candidate으로 생각하고 마찬가지로 하나씩 추가해보면서 AIC, BIC가 10 이상 떨어지는지 확인한다. 

 

Step 6. 최종 모델 도출

모델의 최종 퍼포먼스 Goodness of Fit는 어떻게 확인할까?

 

쉬어가기 >>

2024.01.18 - [대학원 유학일기 US GRAD SCHOOL/유학 일기] - 박사 2학년 2학기 첫주. 널뛰는 나으 마음상태!!

박사 2학년 2학기 첫주. 널뛰는 나으 마음상태!!

 

더 공부하기 >>

2023.06.29 - [컴퓨터 COMPUTER/Data Visualization 데이터시각화] - [지도시각화] Kepler.gl 사용법, 저장 방법 및 접속 링크

[지도시각화] Kepler.gl 사용법, 저장 방법 및 접속 링크