이 문서에는 2017년 6월 9일에 발표된 dplyr 패키지 0.7.0 버젼이 사용되었다. 이 문서의 코드를 실행하는데 문제가 있는 경우 새 버젼으로 업데이트할 것을 권장한다.

연속형변수를 기준으로 같은 크기를 갖는 여러 개의 그룹으로 나누기

moonBook 패키지의 acs 데이터는 857명의 급성관동맥증후군 환자의 데이터이다. age열에는 환자의 만 나이가 기록되어 있는데 이를 기준으로 같은 크기를 갖는 여러 개의 그룹으로 나누어보자. 단 같은 값을 같는 경우 같은 그룹으로 나누기로 한다.

첫번째 방법

먼저 rank()함수로 순위를 정한 후 그 순위에 따라 분류하는 방법이다. moonBook 패키지의 rank2group함수는 이 방법을 사용한다.

require(moonBook)  # acs데이터의 사용을 위해

rank2group=function(y, k = 4) {

    count = length(y)
    z = rank(y, ties.method = "min")
    return(floor((z - 1)/(count/k)) + 1)
}

이 방법을 써서 분류를 해서 group1이라는 열에 저장한다.

acs$group1=rank2group(acs$age,5)

이 방법에 의해 제대로 나누어 졌는지 확인해본다.

require(dplyr)

k=5
acs %>% group_by(group1) %>%
    summarize(
             n=n(),
             min=min(age),
             max=max(age),
             ratio=scales::percent(n/nrow(acs)),
             diff=round(abs(n*100/nrow(acs)-100/k),2)
    ) %>% 
    mutate(sumdiff=sum(diff))
# A tibble: 5 x 7
  group1     n   min   max ratio  diff sumdiff
   <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1      1   190    28    53 22.2%  2.17    7.05
2      2   183    54    61 21.4%  1.35    7.05
3      3   151    62    67 17.6%  2.38    7.05
4      4   168    68    74 19.6%  0.40    7.05
5      5   165    75    91 19.3%  0.75    7.05

다섯개의 그룹으로 나누는 경우 이상적으로는 모두 20%로 나누어지는 것이 좋겠지만 같은 값이 있으므로 조금 달라질 수 있다. 이 경우 20%를 기준으로 약 7.05%에서 차이가 나는 것을 알 수 있다.

cut2를 이용하는 방법

최근 구글링을 해보니 Hmisc::cut2함수가 같은 기능을 한다고 한다. 사용하는 방법은 다음과 같다.

acs$group2=Hmisc::cut2(acs$age,g=5)

정확도를 보기 위해 처음에 사용했던 방법을 함수로 만들어 정확도를 보았다.

validateGroup=function(df,group,var,k=5){
     group=enquo(group)    
     var=enquo(var)    
     df %>% 
         group_by(!!group)%>% 
         summarize(
             n=n(),
             min=min(!!var),
             max=max(!!var),
             ratio=scales::percent(n/nrow(df)),
             diff=round(abs(n*100/nrow(df)-100/k),2)
         ) %>% 
         mutate(sumdiff=sum(diff))
}

위 함수의 소스를 보고 quo()함수나 UQ() 함수(또는 !!)가 생소한 분들은 시간을 내어 dplyr을 이용한 프로그래밍 비니에트(http://dplyr.tidyverse.org/articles/programming.html) 를 읽어보시기 바란다.

validateGroup(acs,group2,age)
# A tibble: 5 x 7
   group2     n   min   max ratio  diff sumdiff
   <fctr> <int> <int> <int> <chr> <dbl>   <dbl>
1 [28,54)   190    28    53 22.2%  2.17    7.05
2 [54,62)   183    54    61 21.4%  1.35    7.05
3 [62,68)   151    62    67 17.6%  2.38    7.05
4 [68,75)   168    68    74 19.6%  0.40    7.05
5 [75,91]   165    75    91 19.3%  0.75    7.05

역시 같은 결과를 보여준다.

Prop.table을 이용하는 방법

정확도를 개선시키기 위해 고민하던 중 prop.table을 이용하여 보았다.

cumsum(prop.table(table(acs$age)))*100
         28          30          35          36          37          38 
  0.1166861   0.2333722   0.4667445   0.7001167   1.1668611   1.5169195 
         39          40          41          42          43          44 
  1.8669778   2.6837806   3.9673279   4.5507585   5.1341890   6.3010502 
         45          46          47          48          49          50 
  6.8844807   8.0513419   9.8016336  11.4352392  13.8856476  14.8191365 
         51          52          53          54          55          56 
 17.5029172  19.7199533  22.1703617  23.6872812  26.7211202  29.5215869 
         57          58          59          60          61          62 
 32.3220537  35.0058343  37.4562427  39.6732789  43.5239207  46.5577596 
         63          64          65          66          67          68 
 48.3080513  51.3418903  54.9591599  57.7596266  61.1435239  64.2940490 
         69          70          71          72          73          74 
 67.0945158  70.8284714  73.3955659  75.8459743  78.4130688  80.7467911 
         75          76          77          78          79          80 
 83.6639440  85.5309218  87.8646441  90.3150525  92.0653442  93.6989498 
         81          82          83          84          85          86 
 95.7992999  96.7327888  97.1995333  97.6662777  98.3663944  98.8331389 
         87          88          89          90          91 
 99.1831972  99.6499417  99.7666278  99.8833139 100.0000000

위의 결과를 보면 처음 20%에 가장 가까운 나이는 52세임을 알 수 있다. 하지만 앞의 두가지 방법에서는 첫 그룹을 53세까지 분류했음을 알 수 있다. 또한 두번쨰 그룹은 40%에 가장 가까운 60세까지 분류하는 것이 바람직하지만 앞의 두 방법은 61세로 분류하였다. 따라서 prop.table의 결과를 보고 이상적인 값에서 가장 가까운 값을 산택하여 분류를 하는 함수를 rank2group2라는 이름으로 제작하였다. 이 함수는 ggiraphExtra 패키지에 포함되어 있다.

rank2group2 = function (x, k = 4) 
{
    temp = cumsum(prop.table(table(x)))
    res = c()
    for (i in 1:(k - 1)) {
        result = which.min(abs(temp - (i/k)))
        res = c(res, result)
    }
    res = as.numeric(names(res))
    res = c(min(x, na.rm = TRUE) - 0.01, res, max(x, na.rm = TRUE))
    temp = cut(x, breaks = res)
    as.numeric(temp)
}

이 함수를 이용하여 분류하고 정확도를 비교해본다.

acs$group3=rank2group2(acs$age,5)    
validateGroup(acs,group3,age)
# A tibble: 5 x 7
  group3     n   min   max ratio  diff sumdiff
   <dbl> <int> <int> <int> <chr> <dbl>   <dbl>
1      1   169    28    52 19.7%  0.28    2.95
2      2   171    53    60   20%  0.05    2.95
3      3   184    61    67 21.5%  1.47    2.95
4      4   168    68    74 19.6%  0.40    2.95
5      5   165    75    91 19.3%  0.75    2.95

앞의 두 방법은 7.05%의 오차를 보였으며 새로운 방법으로 분류하면 오차가 2.95%로 주는 것을 알수 있다.

dplyr 패키지의 ntile 함수를 이용하는 방법

acs$group4=ntile(acs$age,5)
validateGroup(acs,group4,age)
# A tibble: 5 x 7
  group4     n   min   max ratio  diff sumdiff
   <int> <int> <int> <int> <chr> <dbl>   <dbl>
1      1   172    28    53 20.1%  0.07    0.29
2      2   171    53    61   20%  0.05    0.29
3      3   172    61    67 20.1%  0.07    0.29
4      4   171    67    74   20%  0.05    0.29
5      5   171    74    91   20%  0.05    0.29

이 방법을 이용하면 오차가 0.29%로 줄어드나 53세의 경우 group 1로 분류되기도 하고 2로 분류되기도 한다. 마찬가지로 61세, 67세, 74세도 두 그룹으로 나누어 들어간다.

누락된 값이 있는 경우

처음 만들었던 rank2group함수는 누락된 값이 있는 경우 처리를 하지 못한다. acs데이터에 있는 EF열은 심장의 기능을 나타내는 수치로 857명의 환자 중 134명에서 누락되어 있다.

sum(is.na(acs$EF))
[1] 134

이 데이터를 rank2group함수로 분류해 본다.

acs$EFgroup1=rank2group(acs$EF,5)
validateGroup(acs,EFgroup1,EF)
# A tibble: 5 x 7
  EFgroup1     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   175  18.0  50.0 20.4%  0.42    1.22
2        2   168  50.1  57.4 19.6%  0.40    1.22
3        3   173  57.6  61.6 20.2%  0.19    1.22
4        4   170  61.7  67.2 19.8%  0.16    1.22
5        5   171    NA    NA   20%  0.05    1.22

누락이 있는 경우 5번째 그룹으로 분류되는 것을 알 수 있다.

validateGroup(acs[!is.na(acs$EF),],EFgroup1,EF)
# A tibble: 5 x 7
  EFgroup1     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   175  18.0  50.0 24.2%  4.20   29.76
2        2   168  50.1  57.4 23.2%  3.24   29.76
3        3   173  57.6  61.6 23.9%  3.93   29.76
4        4   170  61.7  67.2 23.5%  3.51   29.76
5        5    37  67.4  79.0 5.12% 14.88   29.76

누락이 없는 자료만 대상으로 분류가 정확한지 보면 모두 29.76% 에서 오류가 있는 것을 알 수 있다. 이를 개선하기 위해 rank2group함수가 NA 값이 있어도 제대로 동작하도록 함수를 약간 수정하였다.

rank2group3=function(y, k = 4) {
    x=y[!is.na(y)]
    count = length(x)
    z = rank(x, ties.method = "min")
    y[!is.na(y)]=(floor((z - 1)/(count/k)) + 1)
    y
}

이 함수를 가지고 EF를 기준으로 5그룹으로 분류하여 새로운 변수 EFgroup2에 저장한다.

acs$EFgroup2=rank2group3(acs$EF,5)
validateGroup(acs,EFgroup2,EF)
# A tibble: 6 x 7
  EFgroup2     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   146  18.0  48.2   17%  2.96   19.99
2        2   146  48.3  56.0   17%  2.96   19.99
3        3   145  56.1  59.7 16.9%  3.08   19.99
4        4   145  59.8  63.5 16.9%  3.08   19.99
5        5   141  63.6  79.0 16.5%  3.55   19.99
6       NA   134    NA    NA 15.6%  4.36   19.99
validateGroup(acs[!is.na(acs$EF),],EFgroup2,EF)
# A tibble: 5 x 7
  EFgroup2     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   146  18.0  48.2 20.2%  0.19       1
2        2   146  48.3  56.0 20.2%  0.19       1
3        3   145  56.1  59.7 20.1%  0.06       1
4        4   145  59.8  63.5 20.1%  0.06       1
5        5   141  63.6  79.0 19.5%  0.50       1

누락치는 모두 누락치로 분류하였고 누락치를 제외한 자료의 분류 성능을 보면 오차가 1% 로 줄어들었슴을 알수 있다.

cut2 함수를 사용한 경우

acs$EFgroup3=Hmisc::cut2(acs$EF,g=5)
validateGroup(acs,EFgroup3,EF)
# A tibble: 6 x 7
     EFgroup3     n   min   max ratio  diff sumdiff
       <fctr> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1 [18.0,48.3)   146  18.0  48.2   17%  2.96   19.99
2 [48.3,56.1)   146  48.3  56.0   17%  2.96   19.99
3 [56.1,59.8)   145  56.1  59.7 16.9%  3.08   19.99
4 [59.8,63.6)   145  59.8  63.5 16.9%  3.08   19.99
5 [63.6,79.0]   141  63.6  79.0 16.5%  3.55   19.99
6          NA   134    NA    NA 15.6%  4.36   19.99
validateGroup(acs[!is.na(acs$EF),],EFgroup3,EF)
# A tibble: 5 x 7
     EFgroup3     n   min   max ratio  diff sumdiff
       <fctr> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1 [18.0,48.3)   146  18.0  48.2 20.2%  0.19       1
2 [48.3,56.1)   146  48.3  56.0 20.2%  0.19       1
3 [56.1,59.8)   145  56.1  59.7 20.1%  0.06       1
4 [59.8,63.6)   145  59.8  63.5 20.1%  0.06       1
5 [63.6,79.0]   141  63.6  79.0 19.5%  0.50       1

cut2함수를 사용한 경우의 오차도 1%로 rank2group3함수를 사용한 것과 같다.

rank2group2함수를 사용한 경우

acs$EFgroup4=rank2group2(acs$EF,5)
validateGroup(acs,EFgroup4,EF)
# A tibble: 6 x 7
  EFgroup4     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   146  18.0  48.2   17%  2.96   19.99
2        2   146  48.3  56.0   17%  2.96   19.99
3        3   141  56.1  59.6 16.5%  3.55   19.99
4        4   145  59.7  63.4 16.9%  3.08   19.99
5        5   145  63.5  79.0 16.9%  3.08   19.99
6       NA   134    NA    NA 15.6%  4.36   19.99
validateGroup(acs[!is.na(acs$EF),],EFgroup4,EF)
# A tibble: 5 x 7
  EFgroup4     n   min   max ratio  diff sumdiff
     <dbl> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   146  18.0  48.2 20.2%  0.19       1
2        2   146  48.3  56.0 20.2%  0.19       1
3        3   141  56.1  59.6 19.5%  0.50       1
4        4   145  59.7  63.4 20.1%  0.06       1
5        5   145  63.5  79.0 20.1%  0.06       1

dplyr 패키지의 ntile 함수를 이용하는 방법

acs$EFgroup5=ntile(acs$EF,5)
validateGroup(acs,EFgroup5,EF)
# A tibble: 6 x 7
  EFgroup5     n   min   max ratio  diff sumdiff
     <int> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   145  18.0  48.2 16.9%  3.08      20
2        2   145  48.2  56.0 16.9%  3.08      20
3        3   144  56.0  59.7 16.8%  3.20      20
4        4   145  59.7  63.5 16.9%  3.08      20
5        5   144  63.5  79.0 16.8%  3.20      20
6       NA   134    NA    NA 15.6%  4.36      20
validateGroup(acs[!is.na(acs$EF),],EFgroup5,EF)
# A tibble: 5 x 7
  EFgroup5     n   min   max ratio  diff sumdiff
     <int> <int> <dbl> <dbl> <chr> <dbl>   <dbl>
1        1   145  18.0  48.2 20.1%  0.06    0.34
2        2   145  48.2  56.0 20.1%  0.06    0.34
3        3   144  56.0  59.7 19.9%  0.08    0.34
4        4   145  59.7  63.5 20.1%  0.06    0.34
5        5   144  63.5  79.0 19.9%  0.08    0.34

ntile 함수를 사용한 경우의 오차는 0.34% 이지만 EF의 값이 20.1, 56.0, 59.7, 63.5인 경우 두군에 속하게 된다.

diamonds 데이터

ggplot2패키지에 포함되어 있는 diamonds데이타는 모두 53940개의 다이아몬드에 대한 데이터이다. 이중 price를 기준으로 4개의 군으로 나누어 본다.

먼저 rank2group()함수를 사용한 방법이다.

require(ggplot2)
nrow(diamonds)
[1] 53940
diamonds$priceGroup2=rank2group3(diamonds$price,4)
validateGroup(diamonds,priceGroup2,price,4)
# A tibble: 4 x 7
  priceGroup2     n   min   max ratio  diff sumdiff
        <dbl> <int> <int> <int> <chr> <dbl>   <dbl>
1           1 13490   326   950   25%  0.01    0.06
2           2 13495   951  2401   25%  0.02    0.06
3           3 13470  2402  5324   25%  0.03    0.06
4           4 13485  5325 18823   25%  0.00    0.06

오차의 합이 0.06%를 보이고 있다.

cut2()를 이용한 방법

diamonds$priceGroup1=Hmisc::cut2(diamonds$price,g=4)
validateGroup(diamonds,priceGroup1,price,4)
# A tibble: 4 x 7
   priceGroup1     n   min   max ratio  diff sumdiff
        <fctr> <int> <int> <int> <chr> <dbl>   <dbl>
1 [ 326,  951) 13490   326   950   25%  0.01    0.06
2 [ 951, 2402) 13495   951  2401   25%  0.02    0.06
3 [2402, 5325) 13470  2402  5324   25%  0.03    0.06
4 [5325,18823] 13485  5325 18823   25%  0.00    0.06

마찬가지로 오차의 합이 0.06%를 보이고 있다.

rank2group2()를 이용한 방법

diamonds$priceGroup3=rank2group2(diamonds$price,4)
validateGroup(diamonds,priceGroup3,price,4)
# A tibble: 4 x 7
  priceGroup3     n   min   max ratio  diff sumdiff
        <dbl> <int> <int> <int> <chr> <dbl>   <dbl>
1           1 13483   326   949   25%  0.00    0.04
2           2 13476   950  2400   25%  0.02    0.04
3           3 13496  2401  5324   25%  0.02    0.04
4           4 13485  5325 18823   25%  0.00    0.04

rank2group2() 함수를 이용하면 오차의 합계가 0.04%로 감소하는 것을 알수 있다.

ntile()을 이용하는 방법

diamonds$priceGroup4=ntile(diamonds$price,4)
validateGroup(diamonds,priceGroup3,price,4)
# A tibble: 4 x 7
  priceGroup3     n   min   max ratio  diff sumdiff
        <dbl> <int> <int> <int> <chr> <dbl>   <dbl>
1           1 13483   326   949   25%  0.00    0.04
2           2 13476   950  2400   25%  0.02    0.04
3           3 13496  2401  5324   25%  0.02    0.04
4           4 13485  5325 18823   25%  0.00    0.04

ntile()을 이용하면 rank2group2()와 같은 결과를 얻을 수 있다.

결론

4년전 발표한 moonBook패키지의 rank2group 함수는 결측치가 있는 경우를 생각하지 않고 함수를 만들어 결측치가 있는 경우 제대로 동작하지 않는다. 결측치가 없는 경우 매우 정확한 분류를 하며 그 정확도는 Hmisc패키지의 cut2()함수와 유사하였다. 이번 글을 통해 rank2group() 함수의 버그를 개선한 rank2group3()함수를 소개하였다. 하지만 ggiraphExtra 패키지에 포함되어 있는 prop.table()를 이용한 rank2group2()함수가 가장 오차의 정도가 적었다. 따라서 연속형변수를 기준으로 같은 크기를 갖는 여러 개의 그룹으로 나눌 때에는 rank2group2()함수를 사용하기 바란다. 또한 rank2group2()함수에서 사용한 방법이 가장 좋은 방법은 아닐수 있으니 보다 좋은 결과를 보이는 방법을 아시는 분은 저자에게 연락해주시기 바란다.(cardiomoon@gmail.com)