'모델의 오른쪽 부분 조절하기'에 해당되는 글 1건

  1. 2021.01.31 [R data.table] 선형회귀 모델의 오른쪽 부분(model's right-hand side)의 변수 조합을 간단하게 다루기 4

지난 포스팅에서는 R data.table 패키지에서 .SD 가 무엇이고, .SD와 .SDcols 를 사용해서 특정 칼럼만을 가져오는 방법(을 소개하였습니다. 그리고 lapply()도 같이 사용해서 여러개의 칼럼의 데이터 유형을 변환하는 방법도 같이 설명하였습니다. (rfriend.tistory.com/608)

 

이번 포스팅에서는 지난 포스팅에 이어서, R data.table의 .SD, .SDcols와 lapply(), combn() 함수를 같이 사용하여 '선형회귀 모델의 오른쪽 부분에 여러 변수의 모든 가능한 조합을 적용하여 다루는 방법'을 소개하겠습니다.

 

(1) 설명변수 간 모든 가능한 조합 만들기: combn(x, m, simplify = TRUE, ...)

(2) 설명변수 간 모든 가능한 조합으로 선형회귀모델 적합하여 회귀계수 가져오기

(3) 설명변수 조합별 회귀계수로 막대그래프 그리기: barplot()

 

만약 모델에서 사용하는 칼럼의 개수가 몇 개 안된다면 수작업으로 일일이 변수 이름을 나열해서 적어주면 됩니다. 하지만, 모델에서 사용하는 칼럼 개수가 수십개인 경우(** 몇 개의 칼럼 간의 조합을 사용하는 경우 기하급수적으로 칼럼 개수가 늘어남)에는 키보드를 쳐가면서 일일이 변수 이름(조합)을 나열하기가 시간도 오래 걸리고, 오타 에러를 유발할 위험도 높아집니다. 이럴 경우에 이번에 소개하는 내용이 매우 유용하고 간편하게 사용될 수 있습니다.

 

먼저, data.table 패키지를 불러오고, 예제로 사용할 데이터로 Lahman 패키지에 들어있는 야규 투구 통계 데이터인 'Pitching' 데이터셋을 Data.Table로 만들어보겠습니다.

 

library(data.table)

## Lahman database on baseball
#install.packages("Lahman")
library(Lahman)
data("Pitching")
setDT(Pitching)
str(Pitching)
# Classes 'data.table' and 'data.frame':	47628 obs. of  30 variables:
#   $ playerID: chr  "bechtge01" "brainas01" "fergubo01" "fishech01" ...
# $ yearID  : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1871 ...
# $ stint   : int  1 1 1 1 1 1 1 1 1 1 ...
# $ teamID  : Factor w/ 149 levels "ALT","ANA","ARI",..: 97 142 90 111 90 136 111 56 97 136 ...
# $ lgID    : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ W       : int  1 12 0 4 0 0 0 6 18 12 ...
# $ L       : int  2 15 0 16 1 0 1 11 5 15 ...
# $ G       : int  3 30 1 24 1 1 3 19 25 29 ...
# $ GS      : int  3 30 0 24 1 0 1 19 25 29 ...
# $ CG      : int  2 30 0 22 1 0 1 19 25 28 ...
# $ SHO     : int  0 0 0 1 0 0 0 1 0 0 ...
# $ SV      : int  0 0 0 0 0 0 0 0 0 0 ...
# $ IPouts  : int  78 792 3 639 27 3 39 507 666 747 ...
# $ H       : int  43 361 8 295 20 1 20 261 285 430 ...
# $ ER      : int  23 132 3 103 10 0 5 97 113 153 ...
# $ HR      : int  0 4 0 3 0 0 0 5 3 4 ...
# $ BB      : int  11 37 0 31 3 0 3 21 40 75 ...
# $ SO      : int  1 13 0 15 0 0 1 17 15 12 ...
# $ BAOpp   : num  NA NA NA NA NA NA NA NA NA NA ...
# $ ERA     : num  7.96 4.5 27 4.35 10 0 3.46 5.17 4.58 5.53 ...
# $ IBB     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ WP      : int  7 7 2 20 0 0 1 15 3 44 ...
# $ HBP     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BK      : int  0 0 0 0 0 0 0 2 0 0 ...
# $ BFP     : int  146 1291 14 1080 57 3 70 876 1059 1334 ...
# $ GF      : int  0 0 0 1 0 1 1 0 0 0 ...
# $ R       : int  42 292 9 257 21 0 30 243 223 362 ...
# $ SH      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ GIDP    : int  NA NA NA NA NA NA NA NA NA NA ...
# - attr(*, ".internal.selfref")=<externalptr>

 

(1) 설명변수 간 모든 가능한 조합 만들기: combn(x, m, simplify = TRUE, ...)

 

먼저, 코드 이해를 돕기 위해서 utils 패키지의 combn(x, m, simplify = TRUE, ...) 함수에 대해서 설명하겠습니다.

 

combn() 함수는 투입되는 객체 x 내 원소들 간의 m 개로 이루어진 모든 가능한 조합을 만들어줍니다. 이때 simplify = TRUE (default 설정) 이면 바로 아래의 예처럼 행렬(matrix), 배열(array) 형태로 간소화해서 조합의 결과를 나타내주며, simplify = FALSE 로 설정해주면 두번째 예처럼 원소 간 조합의 결과를 리스트(list) 형태로 해서 위에서 아래로 길게 나타내줍니다.

 

## combn {utils}
## : combn(x, m, FUN = NULL, simplify = TRUE, ...)
## : Generate all combinations of the elements of 'x' taken 'm' at a time
combn(x = letters[1:4], m = 2)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] "a"  "a"  "a"  "b"  "b"  "c" 
# [2,] "b"  "c"  "d"  "c"  "d"  "d"


## simplify = FALSE : returns the combinations in a list format
combn(x = letters[0:4], m = 2, simplify = FALSE)
# [[1]]
# [1] "a" "b"
# 
# [[2]]
# [1] "a" "c"
# 
# [[3]]
# [1] "a" "d"
# 
# [[4]]
# [1] "b" "c"
# 
# [[5]]
# [1] "b" "d"
# 
# [[6]]
# [1] "c" "d"

 

 

앞에서 combn() 함수로 인풋 x 의 모든 원소 간 조합을 만들 수 있다는 것을 알았으므로, 이제 c('yearID', 'teamID', 'G', 'L') 의 4개 설명변수 간 모든 조합을 만들어보겠습니다. 이때 m=0~4 (설명변수 개수) 까지 바꾸어가면서 lapply() 로 combn() 함수에 x = c('yearID', 'teamID', 'G', 'L') 의 4개 원소를 가진 인풋을 적용해서 가능한 모든 조합을 만들어보겠습니다. 그리고 조합의 결과를 unlist() 해서 아래의 models 결과처럼 각 조합이 문자형 벡터(character vector)로 해서 차곡 차곡 리스트에 묶여있도록 하겠습니다.

 

## this generates a list of the 2^k possible extra variables
## for models of the form ERA ~ G + (...)
extra_var <- c('yearID', 'teamID', 'G', 'L')
models <- unlist(
  lapply(0L:length(extra_var), combn, x = extra_var, simplify = FALSE), 
  recursive = FALSE # the function will not recurse beyond the first level items in x
)

models
# [[1]]
# character(0)
# 
# [[2]]
# [1] "yearID"
# 
# [[3]]
# [1] "teamID"
# 
# [[4]]
# [1] "G"
# 
# [[5]]
# [1] "L"
# 
# [[6]]
# [1] "yearID" "teamID"
# 
# [[7]]
# [1] "yearID" "G"     
# 
# [[8]]
# [1] "yearID" "L"     
# 
# [[9]]
# [1] "teamID" "G"     
# 
# [[10]]
# [1] "teamID" "L"     
# 
# [[11]]
# [1] "G" "L"
# 
# [[12]]
# [1] "yearID" "teamID" "G"     
# 
# [[13]]
# [1] "yearID" "teamID" "L"     
# 
# [[14]]
# [1] "yearID" "G"      "L"     
# 
# [[15]]
# [1] "teamID" "G"      "L"     
# 
# [[16]]
# [1] "yearID" "teamID" "G"      "L"

 

 

위의 lapply() 와 combn() 함수를 같이 쓴 코드가 잘 이해가 안가신다면, 바로 아래처럼 m=0~4 까지 일일이 반복적으로 길게 쓴 코드를 lapply() 사용해서 좀더 간결하게 쓴 것이 위의 코드라고 생각하면 되겠습니다.

 

## or, equivalently by manual
unlist(
  list(
    combn(x=extra_var, m=0, simplify=F), 
    combn(x=extra_var, m=1, simplify=F), 
    combn(x=extra_var, m=2, simplify=F), 
    combn(x=extra_var, m=3, simplify=F), 
    combn(x=extra_var, m=4, simplify=F)
  ), 
  recursive = FALSE
)

 

 

 

(2) 설명변수 간 모든 가능한 조합으로 선형회귀모델 적합하여 회귀계수 가져오기

 

R에서 선형회귀모형은 lm(y ~ x, data) 형태의 구문을 사용합니다. 그리고 coef(lm(y ~ x, data)) 로 적합된 회귀모형의 회귀계수에 접근할 수 있습니다.

아래의 Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', rhs)] 코드가 무척 어려워보일 수 있는데요, 하나씩 설명해보자면요, data = .SD, .SDcols = c('W', rhs) 는 Pitching 데이터셋을 참조해서 .SDcols 에 있는 칼럼(즉, 'W' 변수와 (1)번에서 만든 models 설명변수 조합에 해당하는 rhs 변수들)만 선별해서 가져오라는 뜻입니다. 이들 설명변수 조합을 사용해 lm(ERA ~ .) 로 회귀모형을 적합하며, coef(lm(ERA ~ ., data=.SD))['W'] 는 변수 'W'의 적합된 추정 회귀계수를 가져오라는 뜻입니다.

## -- coefficient of 'W' from linear regression model
## using ERA ~ . and data = .SD, then varying which
##   columns are included in .SD allows us to perform this
##   iteration over 16 models succinctly.
##   coef(.)['W'] extracts the W coefficient from each model fit
lm_coef = sapply(models, function(rhs) {
  Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', rhs)]
})



## coefficients of 'W' from the linear regression models above
lm_coef
# W           W           W           W           W           W           W           W 
# -0.20643825 -0.20877606 -0.20763604 -0.11362170 -0.18675857 -0.20689003 -0.09819310 -0.18763789 
# W           W           W           W           W           W           W           W 
# -0.10765791 -0.18214290 -0.12478923 -0.09595359 -0.18197728 -0.11646014 -0.11793832 -0.11143365

 

 

위의 sapply() 함수와 .SD, .SDcols 를 사용하지 않고 각 설명변수 16개 조합별로 일일이 수작업으로 회귀모형을 나열해서 적으려면 아래처럼 단순 작업을 여러번 해야 합니다. 시간도 오래걸리고, 또 오타 에러를 낼 위험이 높습니다.

 

## all combinations of linear regression models by manually
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'G', 'L')]

 

 

(3) 설명변수 조합별 회귀계수로 막대그래프 그리기: barplot()

 

위의 (2)번에서 구한 각 설명변수 16개 조합별 변수 'W'의 회귀계수를 가지고 barplot() 함수를 사용해서 막대그래프(bar plot)를 그려보겠습니다. 

여기서도 'W' 회귀계수가 적합되었을 때 각 모델별로 사용되었던 설명변수 조합을 names.arg = sapply(models, paste, collapse = '/') 로 해서 넣어주었습니다.

 

## -- barplot
# here are 16 visually distinct colors, taken from the list of 20 here:
#   https://sashat.me/2017/01/11/list-of-20-simple-distinct-colors/
col16 = c('#e6194b', '#3cb44b', '#ffe119', '#0082c8',
          '#f58231', '#911eb4', '#46f0f0', '#f032e6',
          '#d2f53c', '#fabebe', '#008080', '#e6beff',
          '#aa6e28', '#fffac8', '#800000', '#aaffc3')

par(oma = c(2, 0, 0, 0))

barplot(lm_coef, names.arg = sapply(models, paste, collapse = '/'), 
        main = "Wins coefficient\n With various covariates", 
        col = col16, las = 2L, cex.names = .8)
        

 

 

참고로, paste() 함수로 'models' 리스트에 들어있는 변수 이름들의 조합을 하나의 합칠 때, collpase = '/' 를 설정해주면 아래처럼 각 리스트 원소 내 여러개의 변수 이름을 하나로 합치면서 사이에 '/'를 넣어줍니다. (paste 함수에서 sep와 collapse 의 차이점은 여기를 참고하세요.)

sapply(models, paste, collapse = '/')
# [1] ""                  "yearID"            "teamID"            "G"                 "L"                
# [6] "yearID/teamID"     "yearID/G"          "yearID/L"          "teamID/G"          "teamID/L"         
# [11] "G/L"               "yearID/teamID/G"   "yearID/teamID/L"   "yearID/G/L"        "teamID/G/L"       
# [16] "yearID/teamID/G/L"

 

 

[ Reference ]

* R data.table vignettes 'Using .SD for Data Analysis'

  : cran.r-project.org/web/packages/data.table/vignettes/datatable-sd-usage.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)

 

 

 

728x90
반응형
Posted by Rfriend
,