지난번 포스팅에서는 Base R, dplyr, SQL 구문과 R data.table 의 기본구문을 비교해봄으로써 data.table 의 기본 구문이 상대적으로 간결하여 코딩을 읽고 쓰기에 좋다는 점을 소개하였습니다. 


이번 포스팅에서는 R data.table 의 기본 구문 DT[i, j, by] 에서 한걸음 더 나아가서 i 행을 선별하고 정렬하기, j 열을 선택하고 계산하기, by 그룹 별로 집계하기 등을 하는데 있어 소소한 팁을 추가로 설명하겠습니다. 


R data.table 의 기본 구문 DT[i, j, by][order]



MASS 패키지에 있는 Cars93 data.frame을 data.table로 변환하고 앞에서부터 8개 변수만 가져와서 예제로 사용하겠습니다. 



library(data.table)


# getting Cars93 dataset

library(MASS)

DT <- as.data.table(Cars93) # convert data.frame to data.table

DT <- DT[,1:8]  # use the first 8 variables


> str(DT)

Classes 'data.table' and 'data.frame': 93 obs. of  8 variables:

 $ Manufacturer: Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...

 $ Model       : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...

 $ Type        : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...

 $ Min.Price   : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...

 $ Price       : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...

 $ Max.Price   : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...

 $ MPG.city    : int  25 18 20 19 22 22 19 16 19 16 ...

 $ MPG.highway : int  31 25 26 26 30 31 28 25 27 25 ...

 - attr(*, ".internal.selfref")=<externalptr> 

>

> head(DT)

   Manufacturer   Model    Type Min.Price Price Max.Price MPG.city MPG.highway

1:        Acura Integra   Small      12.9  15.9      18.8       25          31

2:        Acura  Legend Midsize      29.2  33.9      38.7       18          25

3:         Audi      90 Compact      25.9  29.1      32.3       20          26

4:         Audi     100 Midsize      30.8  37.7      44.6       19          26

5:          BMW    535i Midsize      23.7  30.0      36.2       22          30

6:        Buick Century Midsize      14.2  15.7      17.3       22          31

 




  (1) DT[i, j, by] : i 행 선별하기 (Subset rows in i)


DT[i, j, by] 기본 구문에서 i 행 조건문을 이용하여 [ 차종(Type)이 소형("Small")이고 & 제조사(Manufacturer)가 현대("Hyundai") 인 관측치 ] 를 선별해보겠습니다. 



> # Subset rows in i

> # Get all the cars with "Small" as the Type of "Hyundai" manufacturer

> # No need of prefix 'DT$' each time.

> small_hyundai <- DT[Type == "Small" & Manufacturer == "Hyundai"]

>

> print(small_hyundai)

   Manufacturer   Model  Type Min.Price Price Max.Price MPG.city MPG.highway

1:      Hyundai   Excel Small       6.8     8       9.2       29          33

2:      Hyundai Elantra Small       9.0    10      11.0       22          29

 



이때 Base R 대비 소소한 차이점이 있는데요, Base R에서는 dataframe_nm$variable_nm 처럼 접두사로서 원천이 되는 "dataframe_nm$" 을 명시해주는데요, data.table은 위의 예에서 처럼 DT[...] 의 [...] 안에서는 'DT$variable_nm' 처럼 'DT$' 를 명시 안해줘도 되고, 아니면 아래 예에서처럼 DT[...] 대괄호 안에서 'DT$'를 명시해줘도 됩니다. ('DT$'를 안해줘도 잘 작동하는데 굳이 'DT$' 접두사를 꼬박꼬박 붙일 이유는 없겠지요).  



> # Prefix of 'DT$' works fine. 

> DT[DT$Type == "Small" & DT$Manufacturer == "Hyundai"]

   Manufacturer   Model  Type Min.Price Price Max.Price MPG.city MPG.highway

1:      Hyundai   Excel Small       6.8     8       9.2       29          33

2:      Hyundai Elantra Small       9.0    10      11.0       22          29




Base R 과의 또 하나 차이점은 data.table의 경우 위의 예에서 처럼 i 조건절만 써주고 뒤에 '콤마(,)'를 안찍어 줘도 모든 열(all columns)을 자동으로 선택해온다는 점입니다. 물론 아래의 예와 같이 Base R 구문처럼 i 행 조건절 다음에 '콤마(,)'를 명시적으로 찍어줘서 '모든 열(all columns)을 선택'하게끔 해도 정상 작동하기는 합니다. 



> # a comma after the condition in i works fine. 

> DT[Type == "Small" & Manufacturer == "Hyundai", ]

   Manufacturer   Model  Type Min.Price Price Max.Price MPG.city MPG.highway

1:      Hyundai   Excel Small       6.8     8       9.2       29          33

2:      Hyundai Elantra Small       9.0    10      11.0       22          29

 



슬라이싱을 이용해서 [ 1번째부터 ~ 5번째 행까지의 데이터 가져오기 ] 를 할 수 있습니다. (이때 행 조건 다음에 콤마를 사용하지 않은 DT[1:5] 는 콤마(,)를 추가한 DT[1:5, ] 와 결과는 동일합니다. )



# Get the first 5 rows.

DT_5 <- DT[1:5]  # = DT[1:5, ]


> print(DT_5)

   Manufacturer   Model    Type Min.Price Price Max.Price MPG.city MPG.highway

1:        Acura Integra   Small      12.9  15.9      18.8       25          31

2:        Acura  Legend Midsize      29.2  33.9      38.7       18          25

3:         Audi      90 Compact      25.9  29.1      32.3       20          26

4:         Audi     100 Midsize      30.8  37.7      44.6       19          26

5:          BMW    535i Midsize      23.7  30.0      36.2       22          30

 




이번에는 정렬을 해볼 텐데요, [ 차종(Type) 기준 오름차순으로 정렬(sorting by Type in ascending order)한 후 가격(Price) 기준으로 내림차순으로 정렬(sorting by Price in descending order)하기 ] 를 해보겠습니다. 


data.table 의 정렬은 DT[order()] 구문을 사용하는데요, 단 data.table 은 정렬을 할 때  Base R의 'order' 함수가 아니라 data.table 의 'forder()' 함수를 이용함으로써 훨씬 빠른 속도로 정렬을 수행합니다. 


오름차순(ascending order)이 기본 설정이며, 내림차순(descending order)을 하려면 변수 앞에 '마이너스(-)'를 붙여주면 됩니다. 



# Sort DT by Type in ascending order, and then by Price in descending order

# using data.table's internal fast radix order 'forder()', not 'base::order'.

DT_ordered <- DT[order(Type, -Price)] # '-' for descending order

> head(DT_ordered)

    Manufacturer  Model    Type Min.Price Price Max.Price MPG.city MPG.highway

1: Mercedes-Benz   190E Compact      29.0  31.9      34.9       20          29

2:          Audi     90 Compact      25.9  29.1      32.3       20          26

3:          Saab    900 Compact      20.3  28.7      37.1       20          26

4:         Volvo    240 Compact      21.8  22.7      23.5       21          28

5:    Volkswagen Passat Compact      17.6  20.0      22.4       21          30

6:        Subaru Legacy Compact      16.3  19.5      22.7       23          30





  (2) DT[i, j, by] : j 열 선택하기 (Select column(s) in j)


열을 선택할 때는 DT[i, j, by] 에서 j 부분에 선택하려는 변수(칼럼) 이름을 넣어주면 됩니다. 이때 변수 이름만 넣어주면 결과로 벡터가 반환되며, 아니면 list(column_nm) 처럼 list() 로 싸주면 data.table 이 반환됩니다. 


  • Price 변수를 선택하여 벡터로 반환하기 (Select 'Price' but return it as a Vector)


# Select columns(s) in j

# Select 'Price' column with all rows but return it as a vector

price_vec <- DT[, Price]


> head(price_vec)

[1] 15.9 33.9 29.1 37.7 30.0 15.7

 



  • Price 변수를 선택하여 data.table로 반환하기 (Select 'Price' but return it as a Data.Table)
list(Price) 처럼 변수 이름을 list() 로 감싸주면 data.table 자료형태로 반환합니다. 


# Select 'Price' column with all rows, but return it as a 'data.table'

# wrap the variables within list()

# '.()' is an alias to 'list()'

price_dt <- DT[, list(Price)]


> head(price_dt)

   Price

1:  15.9

2:  33.9

3:  29.1

4:  37.7

5:  30.0

6:  15.7

 



data.table 에서 .() 는 list() 와 동일하게 계산이나 수행 결과를 data.table 로 반환합니다. 많은 경우 코딩을 더 간결하게 하기 위해 list() 보다는 .() 를 더 많이 사용하는 편입니다. (반면, .() 를 처음보는 분은 '이게 도대체 무엇일까?' 하고 궁금할 것 같긴합니다.)



# '.()' is an alias to 'list()'

price_dt2 <- DT[, .(Price)] # .() = list()


> head(price_dt2)

   Price

1:  15.9

2:  33.9

3:  29.1

4:  37.7

5:  30.0

6:  15.7

 




변수를 선택해 올 때 DT[, .(새로운 이름 = 원래 이름)] 형식의 구문을 사용하면 변수 이름을 다른 이름으로 변경 (rename) 할 수도 있습니다.



# Select 'Model' and 'Price' and rename them to 'model_2', 'price_2'

model_price_dt2 <- DT[, .(model_2 = Model, price_2 = Price)]


> head(model_price_dt2)

   model_2 price_2

1: Integra    15.9

2:  Legend    33.9

3:      90    29.1

4:     100    37.7

5:    535i    30.0

6: Century    15.7

 




data.table의 1번째~4번째 칼럼까지를 선택해올 때는 Base R과 동일하게 DT[, 1:4] 또는 직접 칼럼 이름을 입력해서 DT[, c("Manufacturer", "Model", "Type", "Min.Price")] 의 두 가지 방법 모두 가능합니다. 



> head(DT[, 1:4])

   Manufacturer   Model    Type Min.Price

1:        Acura Integra   Small      12.9

2:        Acura  Legend Midsize      29.2

3:         Audi      90 Compact      25.9

4:         Audi     100 Midsize      30.8

5:          BMW    535i Midsize      23.7

6:        Buick Century Midsize      14.2

>

> head(DT[, c("Manufacturer", "Model", "Type", "Min.Price")])

   Manufacturer   Model    Type Min.Price

1:        Acura Integra   Small      12.9

2:        Acura  Legend Midsize      29.2

3:         Audi      90 Compact      25.9

4:         Audi     100 Midsize      30.8

5:          BMW    535i Midsize      23.7

6:        Buick Century Midsize      14.2





  (3) DT[i, j, by] : by 그룹별로 j 열을 계산하거나 실행하기 

                       (Compute or do in j by the group)


DT[i, j, by] 에서 j 열에 대해 by 그룹별로 함수 연산, 계산, 처리를 수행할 수 있습니다. 아래 예에서는
[ 차종(Type) 그룹 별로 가격(Price)의 평균을 구해서 data.table로 반환하기 ] 를 해보겠습니다. 



# Compute or do in j

# calculate mean Price by Type


> DT[, .(mean_price = mean(Price)), Type]

      Type mean_price

1:   Small   10.16667

2: Midsize   27.21818

3: Compact   18.21250

4:   Large   24.30000

5:  Sporty   19.39286

6:     Van   19.10000





DT[i, j, by] 에서 i 행을 선별하고 & j 열에 대해 연산/집계를 by 그룹별로 수행하는 것을 해보겠습니다. 

이때 data.table 은 DT[i, j, by] 처럼 i, j, by 가 DT[...] 안에 모두 들어 있기 때문에 i 행 선별이나 j 열 선택이나 by 그룹별 연산을 각각 분리해서 평가하는 것이 아니라 동시에 고려하여 내부적으로 수행 최적화를 하기 때문에 속도도 빠르고 메모리 효율도 좋습니다


아래 예에서는 i 행 선별할 때 

  (a) 차종이 소형인 행 선별: Type == "Small"

  (b) 차종이 소형이 아닌 행 선별: Type != "Small"

  (c) 차종이 소형, 중형, 컴팩트형인 행 선별: Type %in% c("Small", "Midsize", "Compact")

처럼 조건 유형에 맞게 '==', '!=', '%in%' 등의 연산자(operator)를 사용하였습니다. 



> # Subset in i and do in j

> # Because the three main components of the query (i, j and by) are together inside [...], 

> # data.table can see all three and optimize the query altogether before evaluation, not each separately. 

> # We are able to therefore avoid the entire subset

> DT[Type == "Small", .(mean_price = mean(Price)), Type]

    Type mean_price

1: Small   10.16667

> DT[Type != "Small", .(mean_price = mean(Price)), Type]

      Type mean_price

1: Midsize   27.21818

2: Compact   18.21250

3:   Large   24.30000

4:  Sporty   19.39286

5:     Van   19.10000

> DT[Type %in% c("Small", "Midsize", "Compact"), .(mean_price = mean(Price)), Type]

      Type mean_price

1:   Small   10.16667

2: Midsize   27.21818

3: Compact   18.21250

 




아래 예에서는 조건을 만족하는 관측치의 개수를 세는 두 가지 방법을 소개하겠습니다. 첫번째 방법은 length(변수이름) 을 사용해서 변수 벡터의 길이를 반환하도록 하는 방법이구요, 또다른 방법은 .N 이라는 data.table의 특수 내장변수를 사용해서 관측치 개수를 세는 방법입니다. 아래의 Q1, Q2에 대해 length()와 .N 을 사용한 방법을 차례대로 소개합니다. (역시 .N 이 더 간결하므로 많이 쓰입니다.)


[ Q1: 차종별 관측치 개수 ]

[ Q2: 차종이 소형("Small")이고 가격이 11.0 미만(Price < 11.0) 인 관측치 개수 ]




  • length() 를 사용해서 관측치 개수 세기


> # How many Small cars with price less than 11.0

> DT[, .(N = length(Price)), Type]

      Type  N

1:   Small 21

2: Midsize 22

3: Compact 16

4:   Large 11

5:  Sporty 14

6:     Van  9

> DT[Type == "Small" & Price < 11.0, length(Price)]

[1] 14




  • .N 을 이용하여 관측치 개수 세기


> # .N is a special built-in variable that holds the number of observations in the current group

> DT[, .N, Type]

      Type  N

1:   Small 21

2: Midsize 22

3: Compact 16

4:   Large 11

5:  Sporty 14

6:     Van  9

> DT[Type == "Small" & Price < 11.0, .N]

[1] 14

 



[ Reference ]

* Introduction to R data.table: https://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.html


다음번 포스팅에서는 data.table 의 DT[i, j, by] 에서 by 구문을 사용하여 그룹별로 집계하는 방법을 소개하겠습니다. 


많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. 



728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 R data.table이 무엇이고 또 왜 좋은지, 그리고 data.table을 생성하는 4가지 방법을 소개하였습니다. 


이번 포스팅에서는 R data.table 의 일반 구문(general syntax of data.table)을 Base R, dplyr, SQL 구문과 비교해가면서 설명해보겠습니다. 그러면 상대적으로 data.table 이 매우 간결하고 일관성있으며, 또 SQL 구문과도 비슷한 면이 있어서 익히기 쉽겠구나 하는 인상을 받으실 겁니다. 




R data.table 은 DT[i, j, by][order] 의 기본 구문 형태를 따릅니다. 여기서 DT 는 data.table을 말하여, i, j, by, order 는 

  • i : 데이터 조작/집계의 대상이 되는 행(rows)이 무엇인지를 지정합니다. Base R의 subset() 또는 which() 함수dplyr의 filter() 함수, SQL의 where 구문과 같은 역할을 합니다. 

  • j : 일반적으로 조작/집계의 실행이 일어나는 열(columns)과 실행(action, function)을 말합니다. 기존 열을 수정할 수도 있고, 새로운 조작/집계를 통해 새로운 열을 만들 수도 있습니다. 어떤 R 패키지의 어떤 함수 표현도 모두 사용할 수 있습니다. Base R, dplyr 의 function(column_name), 또는 SQL의 select 구문과 같은 역할을 합니다. 

  • by : j에서 지정한 열(columns)과 실행(action, function)을 그룹 별로(by group) 수행하라고 지정해줍니다. Base R의 tapply()나 by() 함수, dplyr의 group_by() 함수, SQL 의 group by 구문과 같은 역할을 합니다. 

  • order : 결과 테이블을 기준 변수에 따라 오름차순(ascending) 또는 내림차순(descending)으로 정렬할 때 사용합니다. Base R의 sort() 또는 order() 함수, dplyr의 arrange() 함수, SQL의 order by 구문과 같은 역할을 합니다. 


아래에는 소속 그룹을 나타내는 'g' 변수와 정수 값을 가지는 'x' 변수로 구성된 'data'라는 이름의 data.frame, data.table에서, 'x' 가 8 이하인 관측치를 대상으로, 그룹 'g' 별로 'x'의 합을 구한 후에 그룹 'g'를 기준으로 내림차순으로 그 결과를 제시하라는 예제입니다. 


Base R, dplyr, data.table, 그리고 SQL 을 각각 비교해보겠습니다. 


먼저, 예제로 사용할 간단한 data.frame 과 data.table 을 만들어보겠습니다. 변수 'g'는 (그룹 "a" 5개, 그룹 "b" 5개)의 관측치로 이루어져있으며, 변수 'x'는 1~10 까지 순차적인 정수로 이루어져 있습니다.  



# loading packages

library(data.table)

library(dplyr)


# creating data.frame with 2 variables, 10 observations

> data <- data.frame(g = c(rep('a', 5), rep('b', 5)), x = 1:10)

> data

   g  x

1  a  1

2  a  2

3  a  3

4  a  4

5  a  5

6  b  6

7  b  7

8  b  8

9  b  9

10 b 10

>

>

>

# converting data.frame to data.table

> data_dt <- as.data.table(data)

> str(data_dt)

Classes 'data.table' and 'data.frame': 10 obs. of  2 variables:

 $ g: Factor w/ 2 levels "a","b": 1 1 1 1 1 2 2 2 2 2

 $ x: int  1 2 3 4 5 6 7 8 9 10

 - attr(*, ".internal.selfref")=<externalptr> 




이제 아래 과업(task)을 Base R, dplyr, data.table, SQL 로 각각 수행해보겠습니다. 


# ===== 과업 (Task) =====

'data' 라는 데이터셋에서 

'x'가 8 이하인 관측치를 대상으로 

'x'의 합을 구하되

그룹 변수 g 별로 (그룹 'a', 그룹 'b') 구분하여 구하여 

결과를 그룹 변수 g의 내림차순으로 정렬(그룹 'b' 먼저, 그룹 'a' 나중에) 하시오. 


 Base R

dplyr 

data.table 

data2 <- subset(data, 

                subset = (x <= 8))


sort(tapply(data2$x, 

            data2$g, 

            sum), 

     decreasing=TRUE)

data %>% 

  filter(x <= 8) %>% 

  group_by(g) %>% 

  summarise(sum(x)) %>% 

  arrange(desc(g))

data_dt[x<=8, sum(x), by=g][

  order(-g)]

# Out

 b  a 

21 15

# Out

# A tibble: 2 x 2

  g     `sum(x)`

  <fct>    <int>

1 b           21

2 a           15 

# Out

   g V1

1: b 21

2: a 15 



Base R 은 data2$x, data2$g 처럼 각 칼럼 앞에 접두사로 'data2$' 가 필요했는데요, data.table은 각 칼럼별로 접두사 'DT$' 가 필요없습니다. DT[ ] 안에 각 변수 이름을 그냥 써주면 DT의 변수라고 인식을 합니다. 


Base R의 경우 x 가 8 이하인 관측치만 먼저 subset()을 하고 나서, 이후 tapply()로 그룹별로 x에 대해서 sum을 하고, sort()로 내림차순 정렬하는 식으로 하다보니 코드가 길고 복잡합니다. 반면에 data.table의 경우 DT[i, j, by][order] 구문 형식에 따라 i 부분에서 x<=8 로 조건을 만족하는 관측치만 가져온 후, j 부분에서 sum(x) 로 합계를 구하되, by 부분에서 by=g 로 그룹별로 수행하라고 설정하고, [order(-g)] 를 뒤에 설정해주어서 내림차순으로 제시하라고 하였습니다. Base R 대비 무척 간결하고 직관적입니다! 


dplyr 도 Base R 보다는 간결하기는 합니다만, data.table이 조금 더 간결합니다. 


만약 SQL 구문에 이미 익숙한 분이라면 data.table 의 구문이 좀더 수월하게 다가올 수도 있겠습니다. 아래 표의 왼쪽 칸에 SQL 구문으로 똑같은 과업을 수행했을 때의 query 인데요, 오른쪽의 data.table 구문과 비교해보면 어떠신가요? 좀 비슷하지요?!


 SQL

data.table 

SELECT 

g, 

SUM(x) AS x_sum

FROM data_dt

WHERE x <= 8

GROUP BY g 

ORDER BY g DESC;

data_dt[x<=8, 

        sum(x), 

        by=g][

          order(-g)]



이번 포스팅에서는 data.table을 Base R, dplyr, SQL과 비교해보면서 data.table의 구문이 얼마나 간결하고 깜끔한지에 대한 이해를 도모했다면, 다음 포스팅에서는 data.table에 대해서만 DT[i, j, by][order] 기본 구문의 i, j, by, order 별로 좀더 다양한 사용방법에 대해서 한번 더 소개하겠습니다. 


많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. 



728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 R data.table 이 무엇이고, 왜 data.table 이 좋은지, 어떻게 패키지를 설치하는지에 대해서 알아보았습니다. 


그러면 이번 포스팅에서는 샘플 데이터를 가지고 data.table 을 만들어보는 몇 가지 방법을 소개하겠습니다. 


-- 로컬에서 csv 파일, txt 파일을 읽어와서 data.table 만들기

(1) fread() 함수로 csv 파일을 빠르게 읽어와서 R data.table 로 만들기

(2) fwrite() 함수로 R data.table을 csv 파일로 빠르게 쓰기


-- 메모리 상에서 data.table 만들기

(3) data.table() 함수로 R data.table 만들기

(4) data.table() 함수로 기존의 data.frame을 data.table로 변환하기

(5) rbindlist() 함수로 여러개로 쪼개져 있는 파일들을 하나의 data.table로 합치



먼저, data.table의 fread() 함수를 사용해서 로컬 머신에 있는 csv file 을 빠르게(Fast!!) 읽어와서 R data.table로 만들거나, 혹은 반대로 fwrite() 함수를 사용해서 data.table을 빠르게(Fast!!!) csv file로 쓰는 방법을 소개하겠습니다. 




  (1) fread() 함수로 csv 파일을 빠르게 읽어와서 R data.table 로 만들기


data.table의 fread() 함수는 Base R의 read.csv() 함수 또는 read.table() 함수와 유사한 역할을 합니다만, 대신에 fread() 라는 이름에서 짐작할 수 있듯이 매우 빠릅니다(Fast Read)!  벤치마킹 테스트를 해보면 fread() 가 read.table() 보다 40배 이상 빠릅니다!!


파일을 빠르게 읽어와서 data.table 자료로 만들 때 로컬 file path (directory/file_name)를 입력해줘도 되고, https:// 로 시작하는 url 을 입력해줘도 됩니다. 아래 예에서는 UCI Machine Learning Repository에 오픈되어 있는 abalone.data 데이터셋을 url을 통해 읽어와서 data.table로 만들어보았습니다. 


Base R 의 data.frame에서는 문자형 칼럼의 경우 디폴트가 요인형(factor)으로 인식하는데요, data.table은 그렇지 않습니다. 만약 data.table에서 문자형 칼럼을 요인형(factor)로 인식하게 하려면 stringsAsFactors = TRUE 라고 명시적으로 설정을 해줘야 합니다. 



#install.packages("data.table")

library(data.table)


# (1) Fast reading csv file using fread()

url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data'


abalone_dt <- fread(url, 

                    sep=",", 

                    stringsAsFactors=TRUE, 

                    header = FALSE, 

                    col.names=c("sex", "length", "diameter", "height", "whole_weight", 

                                "shucked_weight", "viscera_weight", "shell_weight", "rings"))


> head(abalone_dt)

   sex length diameter height whole_weight shucked_weight viscera_weight shell_weight rings

1:   M  0.455    0.365  0.095       0.5140         0.2245         0.1010        0.150    15

2:   M  0.350    0.265  0.090       0.2255         0.0995         0.0485        0.070     7

3:   F  0.530    0.420  0.135       0.6770         0.2565         0.1415        0.210     9

4:   M  0.440    0.365  0.125       0.5160         0.2155         0.1140        0.155    10

5:   I  0.330    0.255  0.080       0.2050         0.0895         0.0395        0.055     7

6:   I  0.425    0.300  0.095       0.3515         0.1410         0.0775        0.120     8

 


* 참고: https://www.rdocumentation.org/packages/data.table/versions/1.13.0/topics/fread




  (2) fwrite() 함수로 R data.table을 csv 파일로 빠르게 쓰기


data.table의 fwrite() 함수는 Base R의 write.csv() 와 유사한 역할을 합니다만, 역시 fwrite()의 이름에서 짐작할 수 있듯이 매우 매우 빠릅니다!  Base R의 write.csv()보다 data.table의 fwrite()가 약 30~40배 더 빠릅니다!!!  여러개의 CPU가 있을 경우 fwrite()는 이들 CPU를 모두 이용해서 쓰기를 하기 때문입니다. 



# (2) fwrite(): Fast CSV Writer

fwrite(x = abalone_dt, 

       file = "/Users/ihongdon/Downloads/abalone.csv", 

       append = FALSE, 

       sep = ",",

       na = "",

       row.names = FALSE, 

       col.names = TRUE

       )

 


-- [Terminal] command line

(base) ihongdon@lhongdon-a01 ~ % cd Downloads

(base) ihongdon@lhongdon-a01 Downloads % 

(base) ihongdon@lhongdon-a01 Downloads % ls

abalone.csv

(base) ihongdon@lhongdon-a01 Downloads % 

(base) ihongdon@lhongdon-a01 Downloads % 

(base) ihongdon@lhongdon-a01 Downloads % cat abalone.csv 

sex,length,diameter,height,whole_weight,shucked_weight,viscera_weight,shell_weight,rings

M,0.455,0.365,0.095,0.514,0.2245,0.101,0.15,15

M,0.35,0.265,0.09,0.2255,0.0995,0.0485,0.07,7

F,0.53,0.42,0.135,0.677,0.2565,0.1415,0.21,9

...



* 참고: https://www.rdocumentation.org/packages/data.table/versions/1.13.0/topics/fwrite





  (3) data.table() 함수로 R data.table 만들기


data.table 자료구조는 data.frame 의 특성을 그대로 이어받아서 확장한 데이터 구조입니다. 따라서 Base R 에서 사용하는 함수를 모두 사용할 수 있습니다. 


아래 예에서는 2개의 변수와 10개의 관측치를 가진 데이터셋에 대해 data.table() 함수를 사용해서 data.table 클래스 객체를 만들어보았습니다. 



# (3) constructing data.table

data_dt <- data.table(g = c(rep('a', 5), rep('b', 5)), x = 1:10)


>

> str(data_dt)

Classes 'data.table' and 'data.frame': 10 obs. of  2 variables:

 $ g: chr  "a" "a" "a" "a" ...

 $ x: int  1 2 3 4 5 6 7 8 9 10

 - attr(*, ".internal.selfref")=<externalptr> 

>

> head(data_dt)

   g x

1: a 1

2: a 2

3: a 3

4: a 4

5: a 5

6: b 6



* 참고: https://www.rdocumentation.org/packages/data.table/versions/1.13.0/topics/data.table-package




  (4) data.table() 함수로 기존의 data.frame을 data.table로 변환하기


기존 Base R의 data.frame을 data.table() 함수를 사용하면 간단하게 data.table 자료구조로 변환할 수 있습니다. 



> # DataFrame of base R

> data_df <- data.frame(g = c(rep('a', 5), rep('b', 5)), x = 1:10)

> str(data_df)

'data.frame': 10 obs. of  2 variables:

 $ g: Factor w/ 2 levels "a","b": 1 1 1 1 1 2 2 2 2 2

 $ x: int  1 2 3 4 5 6 7 8 9 10

> # (4) converting data.frame to data.table

> data_dt2 <- data.table(data_df)

> str(data_dt2)

Classes 'data.table' and 'data.frame': 10 obs. of  2 variables:

 $ g: Factor w/ 2 levels "a","b": 1 1 1 1 1 2 2 2 2 2

 $ x: int  1 2 3 4 5 6 7 8 9 10

 - attr(*, ".internal.selfref")=<externalptr> 

> head(data_dt2)

   g x

1: a 1

2: a 2

3: a 3

4: a 4

5: a 5

6: b 6





  (5) rbindlist() 함수로 여러개로 쪼개져 있는 파일들을 하나의 data.table로 합치


여러개의 파일을 폴더에서 차례대로 읽어들여서 하나의 data.table로 합치는 방법은 rbindlist() 함수를 사용하면 무척 간단하게 할 수 있습니다. 


예를 들어, 아래처럼 'file_1.txt', 'file_2.txt', 'file_3.txt' 의 3개의 텍스트 파일이 있다고 해보겠습니다. 이럴 경우 먼저 list.files() 함수로 특정 폴더에 들어있는 이들 3개의 파일 이름을 리스트로 만들어 놓습니다. 다음으로 읽어온 파일을 하나씩 리스트로 쌓아놓을 빈 templist 를 만들어놓고, for loop 반복문을 사용하여 텍스트 파일 이름이 들어있는 리스트 f_list 로부터 하나씩 차례대로 파일이름을 읽어와서 paste0(base_dir, f_list[i]) 로 전체의 파일 경로(디렉토리/파일이름.txt)를 만들고, fread() 함수로 하나씩 읽어와서 templist 에 차곡차곡 쌓아놓습니다. 마지막에 rbindlist(templist) 로 모든 리스트를 합쳐서 하나의 data.table 을 만듭니다. 





> # (5) combining multiple files into a data.table using rbindlist()

> base_dir <- '/Users/ihongdon/Downloads/files/'

> f_list <- list.files(base_dir)

> print(f_list)

[1] "file_1.txt" "file_2.txt" "file_3.txt"

> paste0(base_dir, f_list[1])

[1] "/Users/ihongdon/Downloads/files/file_1.txt"

> templist <- list() # black list to store multiple data.tables

> for(i in 1:length(f_list)) {

+   f_path <- paste0(base_dir, f_list[i])

+   templist[[i]] <- fread(f_path)

+ }

> all_dt <- rbindlist(templist)

> str(all_dt)

Classes 'data.table' and 'data.frame': 9 obs. of  3 variables:

 $ x1: int  1 4 7 10 13 16 19 22 25

 $ x2: int  2 5 8 11 14 17 20 23 26

 $ x3: int  3 6 9 12 15 18 21 24 27

 - attr(*, ".internal.selfref")=<externalptr> 

> print(all_dt)

   x1 x2 x3

1:  1  2  3

2:  4  5  6

3:  7  8  9

4: 10 11 12

5: 13 14 15

6: 16 17 18

7: 19 20 21

8: 22 23 24

9: 25 26 27




https://rfriend.tistory.com/225 포스팅에 보면 (a) do.call(rbind, data), (b) ldply(data, rbind), (c) rbind.fill(data), (d) rbindlist(data) 의 4가지 패키지/함수별로 비교를 했을 때 data.table의 rbindlist() 가 월등히 빠르다고 소개를 하였습니다.(비교가 안되게 빠름!!!) 


다음번 포스팅에서는 data.table의 기본 구문에 대해서 소개하겠습니다. 



많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요.



728x90
반응형
Posted by Rfriend
,

R 의 가장 대표적인 데이터 자료 구조 중의 하나를 꼽으라고 하면 행(row)과 열(column)로 구성된 2차원 배열 형태의 테이블인 Data Frame 일 것입니다. 

 

그리고 Data Frame 자료 구조를 조작, 관리, 전처리 할 수 있는 패키지로는 Base R과 dplyr 패키지가 많이 사용이 되는데요, 이번 포스팅에서는 Data Frame의 확장 데이터 구조인 data.table 데이터 구조와 이를 조작, 관리, 처리하는데 사용하는 data.table 패키지에 대해서 소개하려고 합니다. 

 

(1) R data.table은 무엇인가? (What is R data.table?)

(2) 왜 R data.table 인가? (Why R data.table?)

(3) R data.table 설치하기 (How to install R data.table)

(4) R data.table 참고할 수 있는 사이트 (Reference sites of R data.table)

(5) R data.table 구문 소개 포스팅 링크

 

 

  (1) R data.table 은 무엇인가? (What is R data.table?)

 

R data.table 은 (a) 행과 열로 구성된 2차원의 테이블 형태를 가지는 'data.frame'의 확장 데이터 구조(data.table: Extension of 'data.frame')이면서, (b) data.table을 생성, 조작, 처리, 집계할 수 있는 패키지(package) 입니다. 

 

아래에 data.table 패키지를 사용하여 문자열을 가진 'g' 변수와 숫자형 데이터를 가진 'x' 변수에 대해 10개의 관측치를 가지고 있는 간단한 data.table 자료구조를 하나 만들어 봤는데요, str() 로 자료 구조에 대해 확인해보면 'Classes 'data.table' and 'data.frame': 10 obs. of 2 variables:' 라는 설명이 나옵니다. 'data.table'이면서 'data.frame'이라고 나옵니다. 

 

View(data_dt) 로 자료를 확인해봐도 기존에 우리가 알고 있던 data.frame 과 data.table 이 구조가 같다는 것을 눈으로 확인할 수 있습니다. 

 

 

library(data.table)

 

> data <- data.table(g = c(rep('a', 5), rep('b', 5)), x = 1:10)

> str(data)

Classes 'data.table' and 'data.frame': 10 obs. of  2 variables:

 $ g: Factor w/ 2 levels "a","b": 1 1 1 1 1 2 2 2 2 2

 $ x: int  1 2 3 4 5 6 7 8 9 10

 - attr(*, ".internal.selfref")=<externalptr>

 

 

 

> View(data)

 

 

 

R CRAN에는 data.table 패키지에 대해서, 

  • 대용량 데이터(예: 100GB)에 대한 빠른 집계
    (Fast aggregation of large data (e.g. 100GB in RAM))
  • 빠른 정렬된 조인
    (Fast ordered joins)
  • 복사를 사용하지 않고 그룹별로 칼럼에 대한 빠른 합치기/수정/삭제
    (Fast add/modify/delete of columns by group using no copies at all)
  • 사용자 친화적이고 빠른 문자 구분 값의 읽기와 쓰기
    (Friendly and fast character-separated-value read/write)
  • 빠른 개발을 위한 자연스럽고 유연한 구문을 제공
    (Offers a natural and flexible syntax, for faster development.)

한다고 소개하고 있습니다. 

 

 

 

  (2) 왜 R data.table 인가? (Why R data.table?)

 

data.table 의 GitHub 페이지에서 '왜 data.table 인가 (Why data.table?)' 에 대해서 6가지 이유를 말하고 있습니다. (위의 data.table 소개 내용과 유사합니다만, 이유 항목별로 순서가 조금 다르고, 항목도 조금 다르긴 합니다.)

  • 빠르게 쓰고 읽을 수 있는 간결한 구문 (concise syntax: fast to type, fast to read)
  • 빠른 속도 (fast speed)
  • 효율적인 메모리 사용 (memory efficient)
  • 세심한 API 라이프사이클 관리 (careful API lifecycle management)
  • 커뮤니티 (community)
  • 풍부한 기능 (feature rich)

 

첫번째 data.frame 이 좋은 이유로 "빠르게 쓰고 읽을 수 있는 간결한 구문 (concise systax: fast to type, fast to read)"을 들었는데요, 이에 대해서는 Base R, dplyr, data.table 의 세 개 패키지별로 위의 data 에 대해서 그룹 'g' 별로 변수 'x'의 평균(mean) 을 구해보는 구문을 비교해보면 이해하기 쉬울 것 같습니다. 

 

자, 아래 구문을 보면 어떠신가요? 제 눈에는 data.table 이 제일 간결해보이네요!

 

 Base R

dplyr 

data.table 

 tapply(

    data$x, 

    data$y, 

    mean)

 data %>%

    group_by(g) %>%

    summarise(mean(x))

 data[, mean(x), by = 'g']

 

 

두번째 data.frame 의 장점으로 "빠른 속도 (fast speed)" 를 들었는데요, H2Oai 에서 R data.table과 R dplyr 뿐만 아니라 Python pandas, Spark, Julia 등 다른 Database-like 언어의 패키지까지 벤치마킹 테스트를 해놓은 자료가 있어서 아래에 소개합니다. 

 

질문은 5GB 데이터셋에 대해 'id1' 그룹별로 'v1' 칼럼에 대해 합(sum)을 구하라는 집계를 수행하는 것인데요, R data.table 이 12초, R dplyr은 156초, Python pandas는 100초가 걸렸네요. data.table이 겁나게 빠른 것을 알 수 있습니다. 

 

 

* 출처: https://h2oai.github.io/db-benchmark/

 

data.table이 이처럼 빠른데는 여러가지 이유가 있는데요, 그중에서도 핵심은 data.table은 테이블을 수정할 때 레퍼런스(Reference)를 이용할 수 있으므로, 새로운 객체를 다시 만들 필요없이, 원래의 객체를 바로 수정할 수 있다는 점입니다.  

 

 

세번째 data.table의 장점으로 "효율적인 메모리 사용 (memory efficient)"을 꼽았는데요, 역시 H2O ai 에서 위와 동일 질문에 대해 이번에는 '50GB' 의 대용량 데이터에 대해서 수행해보았습니다. 

 

특기할 점은 R data.table은 122초, Spark은 374초가 걸려서 에러없이 수행이 된 반면, R dplyr, Python pandas는 "Out of Memory"가 났다는 점입니다. 

 

* 출처: https://h2oai.github.io/db-benchmark/

 

 

 

네번째로 data.table의 장점으로 "세심한 API 라이프사이클 관리(careful API lifecycle management)"를 들고 있습니다. 

 

data.table 패키지는 Base R에만 의존성(dependency)을 가지고 있습니다. 폐쇄망 환경에서 R 패키지 설치하거나 버전 업그레이드 할 때 의존성 확인하면서 하나 하나 설치하려면 보통 일이 아닌데요, data.table 패키지는 Base R에만 의존성이 있으므로 패키지 설치 및 관리가 무척 수월하겠네요. 

 

 

다섯째 data.table의 장점으로 "커뮤니티 (Community)"를 들고 있습니다. data.table 패키지는 Matt Dowle 외 100여명의 contributor 가 제작 및 관리하고 있는 패키지로서, 소스코드는 https://github.com/Rdatatable/data.table 에서 확인하거나 소스 코드에 기여할 수 있습니다. 

 

 

여섯째 data.table의 장점으로 "풍부한 기능 (feature rich)"을 들고 있는데요, 앞으로 하나씩 같이 알아가보시지요. 

 

 

이렇게 data.table 은 굉장히 매력적인 장점을 가지고 있는데요, 한가지 단점이 있다면 기존에 익숙한 Base R, dplyr 을 놔두고 새로 data.table 패키지의 구문을 배워야 한다는 점입니다. 

 

Base R 쓰다가 dplyr 을 처음 봤을 때의 당혹스러움만큼이나 data.table 구문을 처음 보면 아마도 당황스러울 것입니다. '어, 이게 R 맞나? 내가 지금 무슨 프로그래밍 언어를 보고 있는 거지?' 하고 말이지요. 

 

그래도 다행이라면 구문이 간결해서 배우기 (상대적으로) 쉽다는 점이겠습니다. R data.table은 충분히 시간과 노력을 투자해서 배워볼 만한 값어치가 있다고 생각합니다. 

 

 

 

  (3) R data.table 설치하기 (How to install R data.table)

 

패키지 설치 및 로딩은 다른 패키지와 동일합니다. 

 

 

# install data.table package

install.packages("data.table")

 

# latest development version:

data.table::update.dev.pkg()

 

# loading data.table package

library(data.table)

 

 

 

 

  (4) R data.table 참고할 수 있는 사이트 (Reference sites of R data.table)

 

- data.table 패키지 소소코드: https://github.com/Rdatatable/data.table

- R package 다운로드, 매뉴얼, 튜토리얼: https://cran.r-project.org/web/packages/data.table/ 

 

 

많은 도움이 되었기를 바랍니다. 

다음번 포스팅에서는 R data.table 의 데이터셋을 만드는 3가지 방법을 소개하겠습니다. 

 

 

 

  (5) R data.table 구문 소개 포스팅 링크

 

 

  • 패턴이 있는 여러개의 칼럼을 동시에 길게 녹이고(melt), 옆으로 넓게 주조(cast) 하기
    :
    rfriend.tistory.com/576

  • 그룹별 관측치 개수 별로 Data.Table을 구분해서 생성하기
    :
    rfriend.tistory.com/607

 

 

  • 회귀 모델의 오른쪽 부분(model's right-hand side)의 변수 조합을 일괄 다루기
    : rfriend.tistory.com/609

 

  • R data.table의 조건이 있는 상태에서 Key를 기준으로 데이터셋 합치기(Conditional Joins)
    :
    rfriend.tistory.com/610

 

  • R data.table의 .SD[], by 를 사용해서 그룹별로 부분집합 가져오기 (Group Subsetting)

    : rfriend.tistory.com/611

 

  • R data.table 그룹별 선형회귀모형 적합하고 회귀계수 구하기 (Grouped Regression)
    :
    rfriend.tistory.com/614

 

  • R data.table 이차 인덱스(Secondary indices)를 활용한 빠른 탐색 기반 subsetting
    : rfriend.tistory.com/615

 

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감~

'를 꾹 눌러주세요.

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

 

 

728x90
반응형
Posted by Rfriend
,

자기상관계수(Autocorrelation coefficients), 자기상관그림(Autocorrelation Plot) (a) 데이터가 시간에 의존하는 것 없이 무작위성 (randomness, no time dependence in the data) 을 띠는지 확인하는데 사용됩니다. 그리고 (b) 시계열분석에서 Box-Jenkins ARIMA (Autoregressive Integrated Moving Average) 모형 식별 단계(model identification stage)에서 AR 차수를 구하는데도 사용됩니다. 


이번 포스팅에서는 R을 사용해서 


(1) 자기상관계수(Autocorrelation coefficients)를 계산하고, 

(2) 자기상관계수의 95% 신뢰구간(95% confidence interval)을 계산하고, 

(3) 자기상관그림(Autocorrelation Plot)을 그리는 방법


을 소개하겠습니다. 


본론에 들어가기 전에 먼저, 상관계수(correlation coefficients)와 자기상관계수(autocorrelation coefficients)를 비교해보겠습니다. 


상관계수와 자기상관계수의 유사한 점은 두 연속형 변수 간의 관계를 분석하다는 점입니다. 기본 개념은 공분산을 표준편차로 나누어서 표준화해주며, (자기)상관계수가 -1 ~ 1 사이의 값을 가지게 됩니다. 


상관계수와 자기상관계수가 서로 다른 점은, 상관계수는 특정 동일 시점을 횡단면으로 해서 Y와 다른 X1, X2, ... 변수들 간의 관계를 분석합니다. 반면에 자기상관계수는 동일한 변수(Yt, Yt-1, Yt-2, ...)의 서로 다른 시간 차이 (time lag) 를 두고 관계를 분석하는 것입니다. 


기존에 cross-sectional 관점의 Y와 X 변수들 간의 상관관계 분석에 많이 익숙해져 있는 경우에, 시계열 분석을 공부할 때 보면 자기 자신의 시간 차이에 따른 자기상관관계 관점이 처음엔 헷갈리고 잘 이해가 안가기도 합니다. 관점이 바뀐것일 뿐 어려운 개념은 아니므로 정확하게 이해하고 가면 좋겠습니다. 




  (1) 자기상관계수(Autocorrelation coefficients) 구하고, 자기상관그림 그리기


자기상관계수는 아래의 공식과 같이  로 계산합니다. 



간단한 예로서, 1 ~ 50까지의 시간 t에 대해서 싸인곡선 형태의 주기적인 파동을 띠는 값에 정규확률분포 N(0, 1) 에서 추출한 난수를 더하여 생성한 Y 데이터셋에 대해서 R로 위의 공식을 이용해서 각 시차(time lag)별로 자기상관계수를 구해보겠습니다. 



# sine wave with noise

set.seed(123)

noise <- 0.5*rnorm(50, mean=0, sd=1)

t <- seq(1, 50, 1)

y <- sin(t) + noise


plot(y, type="b")





자기상관계수를 구하고 자기상관그림을 그리는 가장 간단한 방법은 R stats 패키지에 들어있는 acf() 함수를 사용하는 것입니다. 



z <- acf(y, type=c("correlation"), plot=TRUE)


# autocorrelation coefficients per time lag k

round(z$acf, 4)


         [,1]

 [1,]  1.0000

 [2,]  0.3753

 [3,] -0.3244

 [4,] -0.5600

 [5,] -0.4216

 [6,]  0.2267

 [7,]  0.5729

 [8,]  0.3551

 [9,] -0.0952

[10,] -0.5586

[11,] -0.4582

[12,]  0.0586

[13,]  0.3776

[14,]  0.4330

[15,]  0.0961

[16,] -0.4588

[17,] -0.4759


# autocorrelation plot



좀 복잡하기는 합니다만, 위의 자기상관계수를 구하는 공식을 이용해서 R로 직접 사용자 정의 함수를 만들고, for loop 반복문을 이용해서 시차(time lag) k 별로 자기상관계수를 구할 수도 있습니다. 자기공분산을 분산으로 나누어서 표준화해주었습니다. 시차 0(time lag 0) 일 경우에는 자기자신과의 상관관계이므로 자기상관계수는 '1'이 되며, 표준화를 해주었기 때문에 자기상관계수는 -1 <= autocorr(Y, Yt-k) <= 1  사이의 값을 가집니다. 



# User Defined Function of Autocorrelation

acf_func <- function(y, lag_k){

  # y: input vector

  # lag_k : Lag order of k

  

  N = length(y) # total number of observations

  y_bar = mean(y)

  

  # Variance

  var = sum((y - y_bar)^2) / N

  

  # Autocovariance

  auto_cov = sum((y[1:(N-lag_k)] - y_bar) * (y[(1+lag_k):(N)] - y_bar)) / N

  

  # Autocorrelation coefficient = Autocovariance / Variance

  r_k = auto_cov / var

  

  return(r_k)

}



# Compute Autocorrelation per lag (from 0 to 9 in this case)

acf <- data.frame()

for (k in 0:(length(y)-1)){

  acf_k <- round(acf_func(y, lag_k = k), 4)

  acf[k+1, 'lag'] = k

  acf[k+1, 'ACF'] = acf_k

}


> print(acf)

   lag     ACF

1    0  1.0000

2    1  0.3753

3    2 -0.3244

4    3 -0.5600

5    4 -0.4216

6    5  0.2267

7    6  0.5729

8    7  0.3551

9    8 -0.0952

10   9 -0.5586

11  10 -0.4582

12  11  0.0586

13  12  0.3776

14  13  0.4330

15  14  0.0961

16  15 -0.4588

17  16 -0.4759

18  17 -0.0594

19  18  0.2982

20  19  0.3968

21  20  0.1227

22  21 -0.2544

23  22 -0.3902

24  23 -0.1981

25  24  0.1891

26  25  0.3701

27  26  0.1360

28  27 -0.1339

29  28 -0.2167

30  29 -0.1641

31  30  0.0842

32  31  0.2594

33  32  0.1837

34  33  0.0139

35  34 -0.1695

36  35 -0.1782

37  36 -0.0363

38  37  0.1157

39  38  0.1754

40  39  0.0201

41  40 -0.1428

42  41 -0.1014

43  42  0.0000

44  43  0.0750

45  44  0.0641

46  45 -0.0031

47  46 -0.0354

48  47 -0.0405

49  48 -0.0177

50  49 -0.0054

 



이번 예의 경우 시간의 흐름에 따라 싸인 파동 (sine wave) 형태를 보이는 데이터이므로, 당연하게도 시차를 두고 주기적으로 큰 양(+)의 자기상관, 큰 음(-)의 자기 상관을 보이고 있습니다. 




  (2) 자기상관계수의 95% 신뢰구간(95% confidence interval) 계산


자기상관계수의 95% 신뢰구간은 아래의 공식으로 간단하게 구할 수 있습니다. 


자기상관계수의 95% 신뢰구간 = 


이때 N은 샘플 개수, z 는 표준정규분포의 누적분포함수, 는 유의수준 5% 입니다. 


이미 익숙하게 알고 있다시피 유의수준 =0.05 에서의 z = 1.96 입니다. 이번 예에서의 관측치 개수 N = 50 입니다. 따라서 자기상관계수의 95% 신뢰구간은  이 됩니다. 



> # z quantile of 95% confidence level

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)

[1] 1.959964

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)

[1] -1.959964

> # sample size

> N <- length(y)

> print(N)

[1] 50

> # 95% confidence interval of autocorrelation

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] 0.2771808

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] -0.2771808




위의 acf() 함수로 그린 자기상관그림(autocorrelation plot)에서 ACF 축 +0.277 과 -0.277 위치에 수평으로 그어진 점선이 바로 95% 신뢰구간의 상, 하한 선입니다. 


위의 자기상관그림에서 보면 95% 신뢰구간의 상, 하한 점선을 주기적으로 벗어나므로 이 Y 데이터셋은 무작위적이지 않으며 (not random), 시계열적인 자기상관(autocorrelation)이 존재한다고 판단할 수 있습니다. 


만약, 무작위적인 데이터셋에 대해서 자기상관계수를 계산하고 자기상관그림을 그린다면, 시차(time lag)에 따른 자기상관계수가 모두 95% 신뢰구간 상, 하한 점선의 안쪽에 위치할 것입니다. 


많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 R dplyr 패키지의 case_when() 함수를 이용해서 연속형 변수를 여러개의 범주로 구분하여 범주형 변수를 만들어보겠습니다. dplyr 패키지의 case_when() 함수를 사용하면 여러개의 if, else if 조건절을 사용하지 않고도 벡터화해서 쉽고 빠르게 처리를 할 수 있습니다. R dplyr 의 case_when() 함수는 SQL의 case when 절과 유사하다고 보면 되겠습니다. 




간단한 예제로 1~10 까지의 양의 정수를 "2 이하", "3~5", "6~8", "9 이상" 의 4개 범주로 구분을 해보겠습니다. 

(dplyr::case_when()에서 dplyr:: 는 생략해도 되며, dplyr 패키지의 함수를 이용하다는 의미입니다)


case_when(

조건 ~ 할당값, 

조건 ~ 할당값, 

TRUE ~ 할당값)

의 형식으로 작성합니다. 


아래의 예에서는 조건절이 총 4개 사용되었는데요, if, else if, else if, else 등의 조건절문 없이 case_when() 함수의 괄호안에 바로 조건을 나열했고, 마지막에는 앞의 조건절에 모두 해당 안되는 나머지(else)에 대해서 TRUE ~ "9~" 로 지정을 해주었습니다. 



library(dplyr)


x <- 1:10

x


[1]  1  2  3  4  5  6  7  8  9 10



dplyr::case_when(

  x <= 2 ~ "~2",

  x <= 5 ~ "3~5",

  x <= 8 ~ "6~8",

  TRUE ~ "9~"

)


 [1] "~2"  "~2"  "3~5" "3~5" "3~5" "6~8" "6~8" "6~8" "9~"  "9~" 





이때 조건절의 순서가 중요합니다. 복수의 조건절을 나열하면 앞에서 부터 순서대로(in order) 조건에 해당하는 관측치에 대해 값을 할당하게 됩니다. 따라서 만약 TRUE ~ "9~"를 case_when(() 조건절의 제일 앞에 사용하게 되면 1~10까지의 모든 값에 대해 "9~" 를 할당하게 됩니다. 따라서 조건절의 처리 순서를 반드시 고려해서 조건절을 작성해줘야 합니다. 



# order matters!!!

case_when(

  TRUE ~ "9~",

  x <= 2 ~ "~2",

  x <= 5 ~ "3~5",

  x <= 8 ~ "6~8",

)


[1] "9~" "9~" "9~" "9~" "9~" "9~" "9~" "9~" "9~" "9~"

 




case_when() 조건절의 오른쪽(right hand side)의 데이터 유형이 모두 동일해야 합니다. 만약 데이터 유형이 다를 경우 error를 발생합니다. 가령, 아래 예에서는 오른쪽에 character를 반환하게끔 되어있는데 logical 인 NA 가 포함되는 경우 Error가 발생합니다. 이때는 'NA_character_' 를 사용해서 NA가 character로 반환되게끔 해주면 됩니다. 


  • 오른쪽에 문자형(character) 반환하는 경우 NA 값으로는 NA_character_ 사용

 잘못된 사용 예 (오른쪽 데이터 유형 다름)

 올바른 사용 예 (오른쪽 데이터 유형 같음)


# error as NA is logical not character

case_when(

  x <= 2 ~ "~2",

  x <= 5 ~ "3~5",

  x <= 8 ~ "6~8",

  TRUE ~ NA

)


Error: must be a character vector, not a logical vector

Call `rlang::last_error()` to see a backtrace


# use NA_character_

case_when(

  x <= 2 ~ "~2",

  x <= 5 ~ "3~5",

  x <= 8 ~ "6~8",

  TRUE ~ NA_character_

)


[1] "~2"  "~2"  "3~5" "3~5" "3~5" "6~8" "6~8" "6~8" NA    NA 



  • 오른쪽에 숫자형(numeric)을 반환하는 경우 NA 값으로는 NA_real_ 사용

  잘못된 사용 예 (오른쪽 데이터 유형 다름)

 올바른 사용 예 (오른쪽 데이터 유형 같음)


# error as NA is logical not numeric

case_when(

  x <= 2 ~ 2,

  x <= 5 ~ 5,

  x <= 8 ~ 8,

  TRUE ~ NA

)


Error: must be a double vector, not a logical vector

Call `rlang::last_error()` to see a backtrace


# use NA_real_

case_when(

  x <= 2 ~ 2,

  x <= 5 ~ 5,

  x <= 8 ~ 8,

  TRUE ~ NA_real_

)


[1]  2  2  5  5  5  8  8  8 NA NA




dplyr의 case_when() 함수는 mutate() 함수와 함께 사용하면 매우 강력하고 편리하게 여러개의 조건절을 사용해서 새로운 변수를 만들 수 있습니다. 아래는 mtcars 데이터셋의 cyl (실린더 개수)  와 hp (자동차 마력) 의 두 개 변수를 사용해  첫번째 "or" 조건절로 "big" 유형으로 찾고, 두번째 "and" 조건절로 "medium" 유형을 찾으며, 마지막으로 나머지에 대해서는 "small" 유형을 명명해본 예입니다. 



mtcars$name <- row.names(mtcars)


mtcars %>% 

  select(name, mpg, cyl, hp) %>% 

  mutate(

    type = case_when(

      cyl >= 8 | hp >= 180 ~ "big",          # or

      cyl >= 4 & hp >= 120 ~ "medium", # and

      TRUE ~ "small"

    )

  )


                 name  mpg cyl  hp   type

1            Mazda RX4 21.0   6 110  small

2        Mazda RX4 Wag 21.0   6 110  small

3           Datsun 710 22.8   4  93  small

4       Hornet 4 Drive 21.4   6 110  small

5    Hornet Sportabout 18.7   8 175    big

6              Valiant 18.1   6 105  small

7           Duster 360 14.3   8 245    big

8            Merc 240D 24.4   4  62  small

9             Merc 230 22.8   4  95  small

10            Merc 280 19.2   6 123 medium

---- 이하 생략 ----

 




위에서 R dplyr의 case_when() 함수로 진행했던 내용을 PostgreSQL, Greenplum DB에서 하려면 SQL CASE WHEN 문을 아래처럼 사용하면 됩니다. 참고하세요. 



-- PostgreSQL CASE WEHN


SELECT 

   name, 

   mpg, 

   cyl, 

   hp, 

   CASE 

      WHEN (cyl >= 8) OR (hp >= 180) THEN "big"

      WHEN (cyl >= 4) AND (hp >= 120) THEN "median"

      ELSE "small"

   END AS type

FROM mtcars

 



많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. 



728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 R을 사용하여 예측이나 분류 모델링을 할 때 기본적으로 필요한 두가지 작업인 


(1) DataFrame을 Train set, Test set 으로 분할하기 

     (Split a DataFrame into Train and Test set)

   - (1-1) 무작위 샘플링에 의한 Train, Test set 분할 

             (Split of Train, Test set by Random Sampling)

   - (1-2) 순차 샘플링에 의한 Train, Test set 분할 

             (Split of Train, Test set by Sequential Sampling)

   - (1-3) 층화 무작위 샘플링에 의한 Train, Test set 분할 

             (Split of Train, Test set by Stratified Random Sampling)


(2) 여러개의 숫자형 변수를 가진 DataFrame을 표준화하기 

      (Standardization of Numeric Data)

   - (2-1) z-변환 (z-transformation, standardization)

   - (2-2) [0-1] 변환 ([0-1] transformation, normalization)


(3) 여러개의 범주형 변수를 가진 DataFrame에서 가변수 만들기 

      (Getting Dummy Variables)


에 대해서 소개하겠습니다. 



예제로 사용할 Cars93 DataFrame을 MASS 패키지로 부터 불러오겠습니다. 변수가 무척 많으므로 예제를 간단하게 하기 위해 설명변수 X로 'Price', 'Horsepower', 'RPM', 'Length', 'Type', 'Origin' 만을 subset 하여 가져오고, 반응변수 y 로는 'MPG.highway' 변수를 사용하겠습니다. 



# get Cars93 DataFrame from MASS package

library(MASS)

data(Cars93)

str(Cars93)

'data.frame': 93 obs. of 27 variables: $ Manufacturer : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4... $ Model : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1... $ Type : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3... $ Min.Price : num 12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ... $ Price : num 15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ... $ Max.Price : num 18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3... $ MPG.city : int 25 18 20 19 22 22 19 16 19 16 ... $ MPG.highway : int 31 25 26 26 30 31 28 25 27 25 ... $ AirBags : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2... $ DriveTrain : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3... $ Cylinders : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4... $ EngineSize : num 1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ... $ Horsepower : int 140 200 172 172 208 110 170 180 170 200 ... $ RPM : int 6300 5500 5500 5500 5700 5200 4800 4000 4800... $ Rev.per.mile : int 2890 2335 2280 2535 2545 2565 1570 1320 1690... $ Man.trans.avail : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1... $ Fuel.tank.capacity: num 13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ... $ Passengers : int 5 5 5 6 4 6 6 6 5 6 ... $ Length : int 177 195 180 193 186 189 200 216 198 206 ... $ Wheelbase : int 102 115 102 106 109 105 111 116 108 114 ... $ Width : int 68 71 67 70 69 69 74 78 73 73 ... $ Turn.circle : int 37 38 37 37 39 41 42 45 41 43 ... $ Rear.seat.room : num 26.5 30 28 31 27 28 30.5 30.5 26.5 35 ... $ Luggage.room : int 11 15 14 17 13 16 17 21 14 18 ... $ Weight : int 2705 3560 3375 3405 3640 2880 3470 4105 3495... $ Origin : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1... 

$ Make : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...



X <- subset(Cars93, select=c('Price', 'Horsepower', 'RPM', 'Length', 'Type', 'Origin'))

head(X)

A data.frame: 6 × 6
PriceHorsepowerRPMLengthTypeOrigin
<dbl><int><int><int><fct><fct>
15.91406300177Smallnon-USA
33.92005500195Midsizenon-USA
29.11725500180Compactnon-USA
37.71725500193Midsizenon-USA
30.02085700186Midsizenon-USA
15.71105200189MidsizeUSA



table(X$Origin)

USA non-USA 48 45



y <- Cars93$MPG.highway

y

  1. 31
  2.  
  3. 25
  4.  
  5. 26
  6.  
  7. 26
  8.  
  9. 30
  10.  
  11. 31
  12.  
  13. 28
  14.  
  15. 25
  16.  
  17. 27
  18.  
  19. 25
  20.  
  21. 25
  22.  
  23. 36
  24.  
  25. 34
  26.  
  27. 28
  28.  
  29. 29
  30.  
  31. 23
  32.  
  33. 20
  34.  
  35. 26
  36.  
  37. 25
  38.  
  39. 28
  40.  
  41. 28
  42.  
  43. 26
  44.  
  45. 33
  46.  
  47. 29
  48.  
  49. 27
  50.  
  51. 21
  52.  
  53. 27
  54.  
  55. 24
  56.  
  57. 33
  58.  
  59. 28
  60.  
  61. 33
  62.  
  63. 30
  64.  
  65. 27
  66.  
  67. 29
  68.  
  69. 30
  70.  
  71. 20
  72.  
  73. 30
  74.  
  75. 26
  76.  
  77. 50
  78.  
  79. 36
  80.  
  81. 31
  82.  
  83. 46
  84.  
  85. 31
  86.  
  87. 33
  88.  
  89. 29
  90.  
  91. 34
  92.  
  93. 27
  94.  
  95. 22
  96.  
  97. 24
  98.  
  99. 23
  100.  
  101. 26
  102.  
  103. 26
  104.  
  105. 37
  106.  
  107. 36
  108.  
  109. 34
  110.  
  111. 24
  112.  
  113. 25
  114.  
  115. 29
  116.  
  117. 25
  118.  
  119. 26
  120.  
  121. 26
  122.  
  123. 33
  124.  
  125. 24
  126.  
  127. 33
  128.  
  129. 30
  130.  
  131. 23
  132.  
  133. 26
  134.  
  135. 31
  136.  
  137. 31
  138.  
  139. 23
  140.  
  141. 28
  142.  
  143. 30
  144.  
  145. 41
  146.  
  147. 31
  148.  
  149. 28
  150.  
  151. 27
  152.  
  153. 28
  154.  
  155. 26
  156.  
  157. 38
  158.  
  159. 37
  160.  
  161. 30
  162.  
  163. 30
  164.  
  165. 43
  166.  
  167. 37
  168.  
  169. 32
  170.  
  171. 29
  172.  
  173. 22
  174.  
  175. 33
  176.  
  177. 21
  178.  
  179. 30
  180.  
  181. 25
  182.  
  183. 28
  184.  
  185. 28




  (1) DataFrame을 Train set, Test set 으로 분할하기 (Split a DataFrame into Train and Test set)


(1-1) 무작위 샘플링에 의한 Train, Test set 분할 (Split of Train, Test set by Random Sampling)


간단하게 일회성으로 무작위 샘플링 하는 것이면 sample() 함수로 난수를 생성해서 indexing을 해오면 됩니다. 

(* 참고 : https://rfriend.tistory.com/58)



# (1) index for splitting data into Train and Test set

set.seed(1004) # for reprodicibility

train_idx <- sample(1:nrow(X), size=0.8*nrow(X), replace=F) # train-set 0.8, test-set 0.2

test_idx <- (-train_idx)


X_train <- X[train_idx,]

y_train <- y[train_idx]

X_test <- X[test_idx,]

y_test <- y[test_idx]


print(paste0('X_train: ', nrow(X_train)))

print(paste0('y_train: ', length(y_train)))

print(paste0('X_test: ', nrow(X_test)))

print(paste0('y_test: ', length(y_test)))

[Out]:

[1] "X_train: 74" [1] "y_train: 74" [1] "X_test: 19" [1] "y_test: 19"





(1-2) 순차 샘플링에 의한 Train, Test set 분할 (Split of Train, Test set by Sequential Sampling)


시계열 분석을 할 경우 시간 순서(timestamp order)를 유지하는 것이 필요하므로 (1-1)의 무작위 샘플링을 하면 안되며, 시간 순서를 유지한 상태에서 앞서 발생한 시간 구간을 training set, 뒤의(미래의) 시간 구간을 test set 으로 분할합니다. 



# sequential sampling

test_size <- 0.2

test_num <- ceiling(nrow(X) * test_size)

train_num <- nrow(X) - test_num


X_train <- X[1:train_num,]

X_test <- X[(train_num+1):nrow(X),]

y_train <- y[1:train_num]

y_test <- y[(train_num+1):length(y)]





(1-3)  층화 무작위 샘플링에 의한 Train, Test set 분할 (Split of Train, Test set by Stratified Random Sampling)


위의 (1-1)과 (1-2)에서 소개한 무작위 샘플링, 순차 샘플링을 사용한 train, test set split 을 random_split() 이라는 사용자 정의함수(user-defined function)으로 정의하였으며, 층화 무작위 샘플링(stratified random sampling)을 사용한 train_test_split() 사용자 정의 함수도 이어서 정의해 보았습니다. (python sklearn의 train_test_split() 함수의 인자, 반환값이 유사하도록  정의해보았습니다) (* 참고 : https://rfriend.tistory.com/58)



# --- user-defined function of train_test split with random sampling

random_split <- function(X, y

                         , test_size

                         , shuffle

                         , random_state) {

    

    test_num <- ceiling(nrow(X) * test_size)

    train_num <- nrow(X) - test_num

    

    if (shuffle == TRUE) {

        # shuffle == True

        set.seed(random_state) # for reprodicibility

        test_idx <- sample(1:nrow(X), size=test_num, replace=F)

        train_idx <- (-test_idx)

            

        X_train <- X[train_idx,]

        X_test <- X[test_idx,]

        y_train <- y[train_idx]

        y_test <- y[test_idx]

    } else {

        # shuffle == False

        X_train <- X[1:train_num,]

        X_test <- X[(train_num+1):nrow(X),]

        y_train <- y[1:train_num]

        y_test <- y[(train_num+1):length(y)]

    }

    

    return (list(X_train, X_test, y_train, y_test))

}



# --- user defined function of train_test_split() with statified random sampling

train_test_split <- function(X, y

                             , test_size=0.2

                             , shuffle=TRUE

                             , random_state=2004

                             , stratify=FALSE, strat_col=NULL){

                        

    if (stratify == FALSE){ # simple random sampling

        split <- random_split(X, y, test_size, shuffle, random_state)

        X_train <- split[1]

        X_test  <- split[2]

        y_train <- split[3]

        y_test  <- split[4]

    } else { # --- stratified random sampling

        strata <- unique(as.character(X[,strat_col]))

        X_train <- data.frame()

        X_test  <- data.frame()

        y_train <- vector()

        y_test  <- vector()

        for (stratum in strata){

            X_stratum <- X[X[strat_col] == stratum, ]

            y_stratum <- y[X[strat_col] == stratum]

            split_stratum <- random_split(X_stratum, y_stratum, test_size, shuffle, random_state)

            X_train <- rbind(X_train, data.frame(split_stratum[1]))

            X_test  <- rbind(X_test,  data.frame(split_stratum[2]))

            y_train <- c(y_train, unlist(split_stratum[3]))

            y_test  <- c(y_test,  unlist(split_stratum[4]))

        }

    }

    return (list(X_train, X_test, y_train, y_test))

}

 



위에서 정의한 train_test_splie() 사용자 정의 함수를 사용하여 'Origin' ('USA', 'non-USA' 의 두 개 수준을 가진 요인형 변수) 변수를 사용하여 층화 무작위 샘플링을 통한 train, test set 분할 (split of train and test set using stratified random sampling in R) 을 해보겠습니다, 



split_list <- train_test_split(X, y

                               , test_size=0.2

                               , shuffle=TRUE

                               , random_state=2004

                               , stratify=TRUE, strat_col='Origin')


X_train <- data.frame(split_list[1])

X_test  <- data.frame(split_list[2])

y_train <- unlist(split_list[3])

y_test  <- unlist(split_list[4])



print(paste0('Dim of X_train: ', nrow(X_train), ', ', ncol(X_train)))

print(paste0('Dim of X_test:  ', nrow(X_test), ', ', ncol(X_test)))

print(paste0('Length of y_train: ', length(y_train)))

print(paste0('Length of y_test:  ', length(y_test)))

[Out]:
[1] "Dim of X_train: 74, 6"
[1] "Dim of X_test:  19, 6"
[1] "Length of y_train: 74"
[1] "Length of y_test:  19"



X_test

A data.frame: 19 × 6
PriceHorsepowerRPMLengthTypeOrigin
<dbl><int><int><int><fct><fct>
448.0815500168Smallnon-USA
233.92005500195Midsizenon-USA
398.4555700151Smallnon-USA
4012.5905400164Sportynon-USA
329.11725500180Compactnon-USA
538.3825000164Smallnon-USA
4510.01246000172Smallnon-USA
9020.01345800180Compactnon-USA
4212.11025900173Smallnon-USA
1616.31704800178VanUSA
720.81704800200LargeUSA
1140.12956000204MidsizeUSA
739.0745600177SmallUSA
1213.41105200182CompactUSA
823.71804000216LargeUSA
239.2926000174SmallUSA
1716.61654000194VanUSA
7411.11105200181CompactUSA
1415.11604600193SportyUSA


table(X$Origin)

[Out]: USA non-USA 48 45



table(X_test$Origin)

[Out]: USA non-USA 10 9


y_test

  1. [Out]: 33
  2.  
  3. 25
  4.  
  5. 50
  6.  
  7. 36
  8.  
  9. 26
  10.  
  11. 37
  12.  
  13. 29
  14.  
  15. 30
  16.  
  17. 46
  18.  
  19. 23
  20.  
  21. 28
  22.  
  23. 25
  24.  
  25. 41
  26.  
  27. 36
  28.  
  29. 25
  30.  
  31. 33
  32.  
  33. 20
  34.  
  35. 31
  36.  
  37. 28





참고로 (1-1) 무작위 샘플링에 의한 Train, Test set 분할을 위의 (1-3)에서 정의한 train_test_split() 사용자 정의 함수를 사용해서 하면 아래와 같습니다. (shuffle=TRUE)



# split of train, test set by random sampling using train_test_split() function

split_list <- train_test_split(X, y

                               , test_size=0.2

                               , shuffle=TRUE

                               , random_state=2004

                               , stratify=FALSE)


X_train <- data.frame(split_list[1])

X_test  <- data.frame(split_list[2])

y_train <- unlist(split_list[3])

y_test  <- unlist(split_list[4])




참고로 (1-2) 순차 샘플링에 의한 Train, Test set 분할을 위의 (1-3)에서 정의한 train_test_split() 사용자 정의 함수를 사용해서 하면 아래와 같습니다. (shuffle=FALSE)



# split of train, test set by sequential sampling using train_test_split() function

split_list <- train_test_split(X, y

                               , test_size=0.2

                               , shuffle=FALSE

                               , random_state=2004

                               , stratify=FALSE)


X_train <- data.frame(split_list[1])

X_test  <- data.frame(split_list[2])

y_train <- unlist(split_list[3])

y_test  <- unlist(split_list[4])

 




  (2) 여러개의 숫자형 변수를 가진 DataFrame을 표준화하기 (Standardization of Nuemric Data)


(2-1) z-변환 (z-transformation, standardization)


X_train, X_test 데이터셋에서 숫자형 변수(numeric variable)와 범주형 변수(categorical varialble)를 구분한 후에, 숫자형 변수로 이루어진 DataFrame 에 대해서 z-표준화 변환 (z-standardization transformation)을 해보겠습니다. (* 참고 : https://rfriend.tistory.com/52)


여러개의 변수를 가진 DataFrame이므로 X_mean <- apply(X_train_num, 2, mean) 로 Train set의 각 숫자형 변수별 평균을 구하고, X_stddev <- apply(X_train_num, 2, sd) 로 Train set의 각 숫자형 변수별 표준편차를 구했습니다. 


그리고 scale(X_train_num, center=X_mean, scale=X_stddev) 로 Train set의 각 숫자형 변수를 z-표준화 변환을 하였으며, scale(X_test_num, center=X_mean, scale=X_stddev) 로 Test set의 각 숫자형 변수를 z-표준화 변환을 하였습니다. 


이때 조심해야 할 것이 있는데요, z-표준화 변환 시 사용하는 평균(mean)과 표준편차(standard deviation)는 Train set으로 부터 구해서 --> Train set, Test set 에 적용해서 z-표준화를 한다는 점입니다. 왜냐하면 Test set는 미래 데이터(future data), 볼 수 없는 데이터(unseen data) 이므로, 우리가 알 수 있는 집단의 평균과 표준편차는 Train set으로 부터만 얻을 수 있기 때문입니다.  (많은 분석가가 그냥 Train, Test set 구분하기 전에 통채로 scale() 함수 사용해서 표준화를 한 후에 Train, Test set으로 분할을 하는데요, 이는 엄밀하게 말하면 잘못된 순서입니다)



# split numeric, categorical variables

X_train_num <- X_train[, c('Price', 'Horsepower', 'RPM', 'Length')]

X_train_cat <- X_train[, c('Type', 'Origin')]

X_test_num  <- X_test[ , c('Price', 'Horsepower', 'RPM', 'Length')]

X_test_cat  <- X_test[ , c('Type', 'Origin')]


# (1) Z Standardization

# (1-1) using scale() function

X_mean   <- apply(X_train_num, 2, mean)

X_stddev <- apply(X_train_num, 2, sd)


print('---- Mean ----')

print(X_mean)

print('---- Standard Deviation ----')

print(X_stddev)

[Out]:
[1] "---- Mean ----"
     Price Horsepower        RPM     Length 
  20.22703  146.08108 5278.37838  183.67568 
[1] "---- Standard Deviation ----"
     Price Horsepower        RPM     Length 
  9.697073  51.171149 594.730345  14.356620 



X_train_scaled <- scale(X_train_num, center=X_mean, scale = X_stddev)

head(X_train_num_scaled)

A matrix: 6 × 4 of type dbl
PriceHorsepowerRPMLength
1-0.44621989-0.11883811.7177896-0.46498935
41.801881070.50651430.37264220.64947906
51.007827061.21003570.70892910.16189913
41-0.044036690.27200720.8770725-0.60429791
43-0.28122166-0.11883810.54078560.09224485
46-1.05465089-1.05686670.4567139-1.23118639


# note that 'mean' and 'stddev' are calculated using X_train_num dataset (NOT using X_test_num)

X_test_scaled <- scale(X_test_num, center=X_mean, scale = X_stddev)

head(X_test_num_scaled)

A matrix: 6 × 4 of type dbl
PriceHorsepowerRPMLength
44-1.2608987-1.27183150.3726422-1.0918778
21.41001031.05369760.37264220.7887876
39-1.2196491-1.77993030.7089291-2.2760005
40-0.7968411-1.09595120.2044988-1.3704949
30.91501560.50651430.3726422-0.2560265
53-1.2299615-1.2522893-0.4680750-1.3704949



# combine X_train_scaled, X_train_cat

X_train_scaled <- cbind(X_train_num_scaled, X_train_cat)


# combine X_trest_scaled, X_test_cat

X_test_scaled <- cbind(X_test_num_scaled, X_test_cat)





(2-2) [0-1] 변환 ([0-1] transformation, normalization)


각 숫자형 변수별 최소값(min)과 최대값(max)을 구해서 [0-1] 사이의 값으로 변환해보겠습니다. 

(* 참고 : https://rfriend.tistory.com/52)



# (2) [0-1] Normalization

# 0-1 transformation

X_max <- apply(X_train_num, 2, max)

X_min <- apply(X_train_num, 2, min)

X_train_num_scaled <- scale(X_train_num, center = X_min, scale = (X_max - X_min))

X_test_num_scaled <- scale(X_test_num, center = X_min, scale = (X_max - X_min))


head(X_train_num_scaled)

A matrix: 6 × 4 of type dbl
PriceHorsepowerRPMLength
10.155963300.32489450.92592590.4615385
40.555963300.45991560.62962960.6666667
50.414678900.61181430.70370370.5769231
410.227522940.40928270.74074070.4358974
430.185321100.32489450.66666670.5641026
460.047706420.12236290.64814810.3205128



head(X_test_num_scaled)

A matrix: 6 × 4 of type dbl
PriceHorsepowerRPMLength
440.011009170.075949370.62962960.3461538
20.486238530.578059070.62962960.6923077
390.01834862-0.033755270.70370370.1282051
400.093577980.113924050.59259260.2948718
30.398165140.459915610.62962960.5000000
530.016513760.080168780.44444440.2948718


# combine X_train_scaled, X_train_cat

X_train_scaled <- cbind(X_train_num_scaled, X_train_cat)


# combine X_trest_scaled, X_test_cat

X_test_scaled <- cbind(X_test_num_scaled, X_test_cat)





 (3) 여러개의 범주형 변수를 가진 DataFrame에서 가변수 만들기 (Getting Dummy Variables) 


(3-1) caret 패키지의 dummyVars() 함수를 이용하여 DataFrame 내 범주형 변수로부터 가변수 만들기



library(caret)


# fit dummyVars()

dummy <- dummyVars(~ ., data = X_train_cat, fullRank = TRUE)


# predict (transform) dummy variables

X_train_cat_dummy <- predict(dummy, X_train_cat)

X_test_cat_dummy <- predict(dummy, X_test_cat)


head(X_train_cat_dummy)

A matrix: 6 × 6 of type dbl
Type.LargeType.MidsizeType.SmallType.SportyType.VanOrigin.non-USA
001001
010001
000001
010001
010001
010000


head(X_test_cat_dummy)

A matrix: 6 × 6 of type dbl
Type.LargeType.MidsizeType.SmallType.SportyType.VanOrigin.non-USA
75000100
76010000
77100000
78000001
79001000
80001001





(3-2) 조건문 ifelse() 함수를 이용하여 수작업으로 가변수 만들기 

        (creating dummy variables manually using ifelse())


아무래도 (3-1)의 caret 패키지를 이용하는 것 대비 수작업으로 할 경우 범주형 변수의 개수와 범주형 변수 내 class 의 종류 수가 늘어날 수록 코딩을 해야하는 수고가 기하급수적으로 늘어납니다. 그리고 범주형 변수나 class가 가변적인 경우 데이터 전처리 workflow를 자동화하는데 있어서도 수작업의 하드코딩의 경우 에러를 야기하는 문제가 되거나 추가적인 비용이 될 수 있다는 단점이 있습니다. 


범주형 변수 내 범주(category) 혹은 계급(class)이 k 개가 있으면 --> 가변수는 앞에서 부터 k-1 개 까지만 만들었습니다. (회귀모형의 경우 dummy trap 을 피하기 위해)



# check level (class) of categorical variables

unique(X_train_cat$Type)

  1. [Out]: Small
  2.  
  3. Midsize
  4.  
  5. Compact
  6.  
  7. Large
  8.  
  9. Sporty
  10.  
  11. Van

unique(X_train_cat$Origin)

  1. [Out]: non-USA
  2.  
  3. USA


# get dummy variables from train set

X_train_cat_dummy <- data.frame(

    type_small = ifelse(X_train_cat$Type == "Small", 1, 0)

    , type_midsize = ifelse(X_train_cat$Type == "Midsize", 1, 0)

    , type_compact = ifelse(X_train_cat$Type == "Compact", 1, 0)

    , type_large = ifelse(X_train_cat$Type == "Large", 1, 0)

    , type_sporty = ifelse(X_train_cat$Type == "Sporty", 1, 0)

    , origin_nonusa = ifelse(X_train_cat$Origin == "non-USA", 1, 0)

)


head(X_train_cat_dummy)

A data.frame: 6 × 6
type_smalltype_midsizetype_compacttype_largetype_sportyorigin_nonusa
<dbl><dbl><dbl><dbl><dbl><dbl>
100001
010001
001001
010001
010001
010000


# get dummy variables from test set

X_test_cat_dummy <- data.frame(

    type_small = ifelse(X_test_cat$Type == "Small", 1, 0)

    , type_midsize = ifelse(X_test_cat$Type == "Midsize", 1, 0)

    , type_compact = ifelse(X_test_cat$Type == "Compact", 1, 0)

    , type_large = ifelse(X_test_cat$Type == "Large", 1, 0)

    , type_sporty = ifelse(X_test_cat$Type == "Sporty", 1, 0)

    , origin_nonusa = ifelse(X_test_cat$Origin == "non-USA", 1, 0)

)


head(X_test_cat_dummy)

A data.frame: 6 × 6
type_smalltype_midsizetype_compacttype_largetype_sportyorigin_nonusa
<dbl><dbl><dbl><dbl><dbl><dbl>
000010
010000
000100
001001
100000
100001





  (4) 숫자형 변수와 범주형 변수 전처리한 데이터셋을 합쳐서 Train, Test set 완성하기



# combine X_train_scaled, X_train_cat

X_train_preprocessed <- cbind(X_train_num_scaled, X_train_cat_dummy)

head(X_train_preprocessed)

A data.frame: 6 × 10
PriceHorsepowerRPMLengthtype_smalltype_midsizetype_compacttype_largetype_sportyorigin_nonusa
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
0.15596330.34693880.92592590.4615385100001
0.48623850.59183670.62962960.6923077010001
0.39816510.47755100.62962960.5000000001001
0.55596330.47755100.62962960.6666667010001
0.41467890.62448980.70370370.5769231010001
0.15229360.22448980.51851850.6153846010000


 

# combine X_trest_scaled, X_test_cat

X_test_preprocessed <- cbind(X_test_num_scaled, X_test_cat_dummy)

head(X_test_preprocessed)

A data.frame: 6 × 10
PriceHorsepowerRPMLengthtype_smalltype_midsizetype_compacttype_largetype_sportyorigin_nonusa
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
750.188990830.428571430.29629630.70512821000010
760.203669720.591836730.44444440.69230769010000
770.311926610.469387760.37037040.46153846000100
780.390825690.346938780.81481480.55128205001001
790.067889910.122448980.44444440.44871795100000
800.018348620.073469390.66666670.06410256100001




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 시계열 자료의 구성 요인 (time series component factors)으로서 추세 요인, 순환 요인, 계절 요인, 불규칙 요인에 대해서 소개하였습니다. 

 

지난번 포스팅에서는 가법 모형으로 가상의 시계열 자료를 만들었다면(time series composition), ==> 이번 포스팅에서는 반대로 시계열 분해(time series decomposition)를 통해 시계열 자료를 추세(순환)(Trend), 계절성(Seasonality), 잔차(Residual)로 분해를 해보겠습니다. 

 

시계열 분해는 직관적으로 이해하기 쉽고 구현도 쉬워서 시계열 자료 분석의 고전적인 방법론이지만 지금까지도 꾸준히 사용되는 방법론입니다. 시계열 분해의 순서와 방법은 대략 아래와 같습니다. 

 

(1) 시도표 (time series plot)를 보고 시계열의 주기적 반복/계절성이 있는지, 가법 모형(additive model, y = t + s + r)과 승법 모형(multiplicative model, y = t * s * r) 중 무엇이 더 적합할지 판단을 합니다. 

 

(가법 모형을 가정할 시)

(2) 시계열 자료에서 추세(trend)를 뽑아내기 위해서 중심 이동 평균(centered moving average)을 이용합니다. 

 

(3) 원 자료에서 추세 분해값을 빼줍니다(detrend). 그러면 계절 요인과 불규칙 요인만 남게 됩니다. 

 

(4) 다음에 계절 주기 (seasonal period) 로 detrend 이후 남은 값의 합을 나누어주면 계절 평균(average seasonality)을 구할 수 있습니다. (예: 01월 계절 평균 = (2020-01 + 2021-01 + 2022-01 + 2023-01)/4, 02월 계절 평균 = (2020-02 + 2021-02 + 2022-02 + 2023-02)/4). 

 

(5) 원래의 값에서 추세와 계절성 분해값을 빼주면 불규칙 요인(random, irregular factor)이 남게 됩니다. 

 

 

시계열 분해 후에 추세와 계절성을 제외한 잔차(residual, random/irregular factor) 가 특정 패턴 없이 무작위 분포를 띠고 작은 값이면 추세와 계절성으로 모형화가 잘 되는 것이구요, 시계열 자료의 특성을 이해하고 예측하는데 활용할 수 있습니다. 만약 시계열 분해 후의 잔차에 특정 패턴 (가령, 주기적인 파동을 그린다거나, 분산이 점점 커진다거나 등..) 이 존재한다면 잔차에 대해서만 다른 모형을 추가로 적합할 수도 있겠습니다. 

 

 

예제로 사용할 시계열 자료로서 '1차 선형 추세 + 4년 주기 순환 + 1년 단위 계절성 + 불규칙 noise' 의 가법 모형 (additive model)으로 시계열 데이터를 만들어보겠습니다. 

 



import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


dates = pd.date_range('2020-01-01', periods=48, freq='M')
dates
[Out]:
DatetimeIndex(['2020-01-31', '2020-02-29', '2020-03-31', '2020-04-30', '2020-05-31', '2020-06-30', '2020-07-31', '2020-08-31', '2020-09-30', '2020-10-31', '2020-11-30', '2020-12-31', '2021-01-31', '2021-02-28', '2021-03-31', '2021-04-30', '2021-05-31', '2021-06-30', '2021-07-31', '2021-08-31', '2021-09-30', '2021-10-31', '2021-11-30', '2021-12-31', '2022-01-31', '2022-02-28', '2022-03-31', '2022-04-30', '2022-05-31', '2022-06-30', '2022-07-31', '2022-08-31', '2022-09-30', '2022-10-31', '2022-11-30', '2022-12-31', '2023-01-31', '2023-02-28', '2023-03-31', '2023-04-30', '2023-05-31', '2023-06-30', '2023-07-31', '2023-08-31', '2023-09-30', '2023-10-31', '2023-11-30', '2023-12-31'], dtype='datetime64[ns]', freq='M')


# additive model: trend + cycle + seasonality + irregular factor


timestamp = np.arange(len(dates))
trend_factor = timestamp*1.1
cycle_factor = 10*np.sin(np.linspace(0, 3.14*2, 48))
seasonal_factor = 7*np.sin(np.linspace(0, 3.14*8, 48))
np.random.seed(2004)
irregular_factor = 2*np.random.randn(len(dates))


df = pd.DataFrame({'timeseries': trend_factor + cycle_factor + seasonal_factor + irregular_factor, 
                   'trend': trend_factor, 
                   'cycle': cycle_factor, 
                   'trend_cycle': trend_factor + cycle_factor,
                   'seasonal': seasonal_factor, 
                   'irregular': irregular_factor},
                   index=dates)


df
[Out]:
  timeseries trend cycle trend_cycle seasonal irregular
2020-01-31 2.596119 0.0 0.000000 0.000000 0.000000 2.596119
2020-02-29 6.746160 1.1 1.332198 2.432198 3.565684 0.748278
2020-03-31 8.112100 2.2 2.640647 4.840647 6.136825 -2.865371
2020-04-30 8.255941 3.3 3.902021 7.202021 6.996279 -5.942358
2020-05-31 16.889655 4.4 5.093834 9.493834 5.904327 1.491495
2020-06-30 16.182357 5.5 6.194839 11.694839 3.165536 1.321981
2020-07-31 14.128087 6.6 7.185409 13.785409 -0.456187 0.798865
2020-08-31 11.943313 7.7 8.047886 15.747886 -3.950671 0.146099
2020-09-30 9.728095 8.8 8.766892 17.566892 -6.343231 -1.495567
2020-10-31 12.483489 9.9 9.329612 19.229612 -6.966533 0.220411
2020-11-30 12.141808 11.0 9.726013 20.726013 -5.646726 -2.937480
2020-12-31 15.143334 12.1 9.949029 22.049029 -2.751930 -4.153764
2021-01-31 21.774516 13.2 9.994684 23.194684 0.910435 -2.330604
2021-02-28 28.432892 14.3 9.862164 24.162164 4.318862 -0.048134
2021-03-31 32.350583 15.4 9.553832 24.953832 6.522669 0.874082
2021-04-30 30.596556 16.5 9.075184 25.575184 6.907169 -1.885797
2021-05-31 32.510523 17.6 8.434753 26.034753 5.365118 1.110653
2021-06-30 30.425519 18.7 7.643955 26.343955 2.326624 1.754939
2021-07-31 24.300958 19.8 6.716890 26.516890 -1.360813 -0.855119
2021-08-31 20.450917 20.9 5.670082 26.570082 -4.668691 -1.450475
2021-09-30 18.870881 22.0 4.522195 26.522195 -6.674375 -0.976939
2021-10-31 21.326310 23.1 3.293690 26.393690 -6.818438 1.751059
2021-11-30 22.902448 24.2 2.006469 26.206469 -5.060699 1.756678
2021-12-31 26.620578 25.3 0.683478 25.983478 -1.891426 2.528526
2022-01-31 27.626499 26.4 -0.651696 25.748304 1.805404 0.072791
2022-02-28 31.858923 27.5 -1.975253 25.524747 4.998670 1.335506
2022-03-31 35.930469 28.6 -3.263598 25.336402 6.797704 3.796363
2022-04-30 30.177870 29.7 -4.493762 25.206238 6.700718 -1.729087
2022-05-31 30.016165 30.8 -5.643816 25.156184 4.734764 0.125217
2022-06-30 26.591729 31.9 -6.693258 25.206742 1.448187 -0.063200
2022-07-31 21.118481 33.0 -7.623379 25.376621 -2.242320 -2.015820
2022-08-31 16.636031 34.1 -8.417599 25.682401 -5.307397 -3.738973
2022-09-30 17.682613 35.2 -9.061759 26.138241 -6.892132 -1.563496
2022-10-31 21.163298 36.3 -9.544375 26.755625 -6.554509 0.962182
2022-11-30 22.455672 37.4 -9.856844 27.543156 -4.388699 -0.698786
2022-12-31 26.919529 38.5 -9.993595 28.506405 -0.998790 -0.588086
2023-01-31 33.964623 39.6 -9.952191 29.647809 2.669702 1.647112
2023-02-28 37.459776 40.7 -9.733369 30.966631 5.593559 0.899586
2023-03-31 40.793766 41.8 -9.341031 32.458969 6.957257 1.377540
2023-04-30 43.838415 42.9 -8.782171 34.117829 6.380433 3.340153
2023-05-31 41.301780 44.0 -8.066751 35.933249 4.023975 1.344556
2023-06-30 39.217866 45.1 -7.207526 37.892474 0.545147 0.780245
2023-07-31 35.125502 46.2 -6.219813 39.980187 -3.085734 -1.768952
2023-08-31 33.841926 47.3 -5.121219 42.178781 -5.855940 -2.480916
2023-09-30 38.770511 48.4 -3.931329 44.468671 -6.992803 1.294643
2023-10-31 37.371216 49.5 -2.671356 46.828644 -6.179230 -3.278198
2023-11-30 46.587633 50.6 -1.363760 49.236240 -3.642142 0.993536
2023-12-31 46.403326 51.7 -0.031853 51.668147 -0.089186 -5.175634

 

 

 

  (1) Python을 이용한 시계열 분해 (Time series decomposition using Python)

 

Python의 statsmodels 라이브러리를 사용해서 가법 모형(additive model) 가정 하에 시계열 분해를 해보겠습니다. 

 



from statsmodels.tsa.seasonal import seasonal_decompose


ts = df.timeseries
result = seasonal_decompose(ts, model='additive')


plt.rcParams['figure.figsize'] = [12, 8]
result.plot()
plt.show()



 

 

 

원래의 시계열 구성요소(추세+순환, 계절성, 불규칙 요인)와 시계열 분해(time series decomposition)를 통해 분리한 추세(&순환), 계절성, 잔차(불규칙 요인)를 겹쳐서 그려보았습니다. (즉, 원래 데이터의 추세요인과 시계열 분해를 통해 분리한 추세를 겹쳐서 그려보고, 원래 데이터의 계절요인과 시계열 분해를 통해 분리한 계절을 겹쳐서 그려보고, 원래 데이터의 불규칙 요인과 시계열 분해를 통해 분리한 잔차를 겹쳐서 그려봄) 

 

원래의 데이터와 얼추 비슷하게, 그럴싸하게 시계열 분해를 한 것처럼 보이지요?   

 



# ground truth & timeseries decompostion all together
# -- observed data
plt.figure(figsize=(12, 12))
plt.subplot(4,1, 1)
result.observed.plot()
plt.grid(True)
plt.ylabel('Observed', fontsize=14)


# -- trend & cycle factor
plt.subplot(4, 1, 2)
result.trend.plot()        # from timeseries decomposition
df.trend_cycle.plot()     # ground truth
plt.grid(True)
plt.ylabel('Trend', fontsize=14)


# -- seasonal factor
plt.subplot(4, 1, 3)
result.seasonal.plot()  # from timeseries decomposition
df.seasonal.plot()        # ground truth
plt.grid(True)
plt.ylabel('Seasonality', fontsize=14)


# -- irregular factor (noise)
plt.subplot(4, 1, 4)
result.resid.plot()    # from timeseries decomposition
df.irregular.plot()    # ground truth
plt.grid(True)
plt.ylabel('Residual', fontsize=14)


plt.show()

 

 

 

원래의 관측치(observed), 추세(trend), 계절성(seasonal), 잔차(residual) 데이터 아래처럼 시계열 분해한 객체에서 obsered, trend, seasonal, resid 라는 attributes 를 통해서 조회할 수 있습니다. 

 

Observed  Trend ( & Cycle)
 print(result.observed)
[Out]:
2020-01-31 2.596119
2020-02-29 6.746160
2020-03-31 8.112100
2020-04-30 8.255941
2020-05-31 16.889655
2020-06-30 16.182357
2020-07-31 14.128087
2020-08-31 11.943313
2020-09-30 9.728095
2020-10-31 12.483489
2020-11-30 12.141808
2020-12-31 15.143334
2021-01-31 21.774516
2021-02-28 28.432892
2021-03-31 32.350583
2021-04-30 30.596556
2021-05-31 32.510523
2021-06-30 30.425519
2021-07-31 24.300958
2021-08-31 20.450917
2021-09-30 18.870881
2021-10-31 21.326310
2021-11-30 22.902448
2021-12-31 26.620578
2022-01-31 27.626499
2022-02-28 31.858923
2022-03-31 35.930469
2022-04-30 30.177870
2022-05-31 30.016165
2022-06-30 26.591729
2022-07-31 21.118481
2022-08-31 16.636031
2022-09-30 17.682613
2022-10-31 21.163298
2022-11-30 22.455672
2022-12-31 26.919529
2023-01-31 33.964623
2023-02-28 37.459776
2023-03-31 40.793766
2023-04-30 43.838415
2023-05-31 41.301780
2023-06-30 39.217866
2023-07-31 35.125502
2023-08-31 33.841926
2023-09-30 38.770511
2023-10-31 37.371216
2023-11-30 46.587633
2023-12-31 46.403326
Freq: M, Name: timeseries,

dtype: float64
print(result.trend)
[Out]
2020-01-31 NaN
2020-02-29 NaN
2020-03-31 NaN
2020-04-30 NaN
2020-05-31 NaN
2020-06-30 NaN
2020-07-31 11.994971
2020-08-31 13.697685
2020-09-30 15.611236
2020-10-31 17.552031
2020-11-30 19.133760
2020-12-31 20.378094
2021-01-31 21.395429
2021-02-28 22.173782
2021-03-31 22.909215
2021-04-30 23.658616
2021-05-31 24.475426
2021-06-30 25.402005
2021-07-31 26.124056
2021-08-31 26.510640
2021-09-30 26.802553
2021-10-31 26.934270
2021-11-30 26.812893
2021-12-31 26.549220
2022-01-31 26.256876
2022-02-28 25.965319
2022-03-31 25.756854
2022-04-30 25.700551
2022-05-31 25.675143
2022-06-30 25.668984
2022-07-31 25.945528
2022-08-31 26.442986
2022-09-30 26.878992
2022-10-31 27.650819
2022-11-30 28.690242
2022-12-31 29.686565
2023-01-31 30.796280
2023-02-28 32.096818
2023-03-31 33.692393
2023-04-30 35.246385
2023-05-31 36.927214
2023-06-30 38.744537
2023-07-31 NaN
2023-08-31 NaN
2023-09-30 NaN
2023-10-31 NaN
2023-11-30 NaN
2023-12-31 NaN 

Freq: M, Name: timeseries,
dtype: float64 

 

 

 Seasonality   Residual (Noise)
print(result.seasonal)
[Out]:
2020-01-31 1.501630
2020-02-29 5.701170
2020-03-31 8.768065
2020-04-30 6.531709
2020-05-31 5.446174
2020-06-30 2.002476
2020-07-31 -1.643064
2020-08-31 -6.011071
2020-09-30 -7.807785
2020-10-31 -5.858728
2020-11-30 -5.849710
2020-12-31 -2.780867
2021-01-31 1.501630
2021-02-28 5.701170
2021-03-31 8.768065
2021-04-30 6.531709
2021-05-31 5.446174
2021-06-30 2.002476
2021-07-31 -1.643064
2021-08-31 -6.011071
2021-09-30 -7.807785
2021-10-31 -5.858728
2021-11-30 -5.849710
2021-12-31 -2.780867
2022-01-31 1.501630
2022-02-28 5.701170
2022-03-31 8.768065
2022-04-30 6.531709
2022-05-31 5.446174
2022-06-30 2.002476
2022-07-31 -1.643064
2022-08-31 -6.011071
2022-09-30 -7.807785
2022-10-31 -5.858728
2022-11-30 -5.849710
2022-12-31 -2.780867
2023-01-31 1.501630
2023-02-28 5.701170
2023-03-31 8.768065
2023-04-30 6.531709
2023-05-31 5.446174
2023-06-30 2.002476
2023-07-31 -1.643064
2023-08-31 -6.011071
2023-09-30 -7.807785
2023-10-31 -5.858728
2023-11-30 -5.849710
2023-12-31 -2.780867
Freq: M, Name: timeseries,

dtype: float64
print(result.resid)
[Out]:
2020-01-31 NaN
2020-02-29 NaN
2020-03-31 NaN
2020-04-30 NaN
2020-05-31 NaN
2020-06-30 NaN
2020-07-31 3.776179
2020-08-31 4.256699
2020-09-30 1.924644
2020-10-31 0.790186
2020-11-30 -1.142242
2020-12-31 -2.453893
2021-01-31 -1.122544
2021-02-28 0.557940
2021-03-31 0.673303
2021-04-30 0.406231
2021-05-31 2.588922
2021-06-30 3.021039
2021-07-31 -0.180034
2021-08-31 -0.048653
2021-09-30 -0.123887
2021-10-31 0.250769
2021-11-30 1.939265
2021-12-31 2.852225
2022-01-31 -0.132007
2022-02-28 0.192434
2022-03-31 1.405550
2022-04-30 -2.054390
2022-05-31 -1.105152
2022-06-30 -1.079730
2022-07-31 -3.183983
2022-08-31 -3.795884
2022-09-30 -1.388594
2022-10-31 -0.628793
2022-11-30 -0.384861
2022-12-31 0.013830
2023-01-31 1.666713
2023-02-28 -0.338212
2023-03-31 -1.666692
2023-04-30 2.060321
2023-05-31 -1.071608
2023-06-30 -1.529146
2023-07-31 NaN
2023-08-31 NaN
2023-09-30 NaN
2023-10-31 NaN
2023-11-30 NaN
2023-12-31 NaN 

Freq: M, Name: timeseries,
dtype: float64 

 

 


# export to csv file
df.to_csv('ts_components.txt', sep=',', index=False)
 

 

 

 

 

  (2) R을 이용한 시계열 분해 (Time series Decomposition using R)

 

위에서 가법 모형을 적용해서 Python으로 만든 시계열 자료를 text 파일로 내보낸 후, 이를 R에서 읽어서 시계열 분해 (time series decomposition)를 해보겠습니다. 

 

ts_components.txt
다운로드

 



# read time series text file
df <- read.table('ts_components.txt', sep=',', header=T)
head(df)
A data.frame: 6 × 6
timeseries trend cycle trend_cycle seasonal irregular
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
2.596119 0.0 0.000000 0.000000 0.000000 2.5961193
6.746160 1.1 1.332198 2.432198 3.565684 0.7482781
8.112100 2.2 2.640647 4.840647 6.136825 -2.8653712
8.255941 3.3 3.902021 7.202021 6.996279 -5.9423582
16.889655 4.4 5.093834 9.493834 5.904327 1.4914946
16.182357 5.5 6.194839 11.694839 3.165536 1.3219812
 

 

 

이렇게 불러와서 만든 df DataFrame의 칼럼 중에서 시계열 분해를 할 'timeseries' 칼럼만을 가져와서 ts() 함수를 사용하여 1년 12개월 이므로 frequency =  12로 설정해 R의 시계열 자료 형태로 변환합니다. 

 

그 다음에 decompose() 함수를 사용하여 시계열 분해를 하는데요, 이때 가법 모형 (additive model)을 적용할 것이므로 decompose(ts, "additive") 라고 설정해줍니다. 

 

시계열 분해를 한 결과를 모아놓은 리스트 ts_decompose 객체를 프린트해보면 원래의 값 $x, 계절 요인 $seasonal, 추세(&순환) 요인 $trend, 불규칙 요인 $random 분해값이 순서대로 저장되어 있음을 알 수 있습니다. 

 



# transforming data to time series with 12 months frequency
ts <- ts(df$timeseries, frequency = 12) # 12 months


# time series decomposition
ts_decompose <- decompose(ts, "additive") # additive model


# decomposition results
ts_decompose
[Out]:
$x Jan Feb Mar Apr May Jun Jul
1 2.596119 6.746160 8.112100 8.255941 16.889655 16.182357 14.128087
2 21.774516 28.432892 32.350583 30.596556 32.510523 30.425519 24.300958
3 27.626499 31.858923 35.930469 30.177870 30.016165 26.591729 21.118481
4 33.964623 37.459776 40.793766 43.838415 41.301780 39.217866 35.125502

Aug Sep Oct Nov Dec
1 11.943313 9.728095 12.483489 12.141808 15.143334
2 20.450917 18.870881 21.326310 22.902448 26.620578
3 16.636031 17.682613 21.163298 22.455672 26.919529
4 33.841926 38.770511 37.371216 46.587633 46.403326

$seasonal Jan Feb Mar Apr May Jun Jul
1 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064
2 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064
3 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064
4 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064

Aug Sep Oct Nov Dec
1 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867
2 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867
3 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867
4 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867

$trend Jan Feb Mar Apr May Jun Jul Aug
1 NA NA NA NA NA NA 11.99497 13.69769
2 21.39543 22.17378 22.90922 23.65862 24.47543 25.40200 26.12406 26.51064
3 26.25688 25.96532 25.75685 25.70055 25.67514 25.66898 25.94553 26.44299
4 30.79628 32.09682 33.69239 35.24639 36.92721 38.74454 NA NA

Sep Oct Nov Dec
1 15.61124 17.55203 19.13376 20.37809
2 26.80255 26.93427 26.81289 26.54922
3 26.87899 27.65082 28.69024 29.68657
4 NA NA NA NA

$random Jan Feb Mar Apr May Jun
1 NA NA NA NA NA NA
2 -1.12254398 0.55793990 0.67330316 0.40623129 2.58892205 3.02103851
3 -0.13200705 0.19243403 1.40555039 -2.05439035 -1.10515213 -1.07973023
4 1.66671294 -0.33821202 -1.66669164 2.06032097 -1.07160802 -1.52914637

Jul Aug Sep Oct Nov Dec
1 3.77617915 4.25669908 1.92464365 0.79018590 -1.14224232 -2.45389325
2 -0.18003389 -0.04865275 -0.12388741 0.25076875 1.93926482 2.85222467
3 -3.18398336 -3.79588443 -1.38859434 -0.62879275 -0.38486059 0.01383049
4 NA NA NA NA NA NA

$figure
[1] 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064
[8] -6.011071 -7.807785 -5.858728 -5.849710 -2.780867

$type
[1] "additive" attr(,"class")
[1] "decomposed.ts"

 

 

 

 

위의 분해 결과가 숫자만 잔뜩 들어있으니 뭐가 뭔지 잘 눈에 안들어오지요?  그러면 이제 원래의 값 (observed)과 시계열 분해된 결과인 trend, seasonal, random 을 plot() 함수를 사용하여 다같이 시각화해보겠습니다. 

 



# change plot in jupyter
library(repr)


# Change plot size to 12 x 10
options(repr.plot.width=12, repr.plot.height=10)


plot(ts_decompose)
 

 

 

위의 분해 결과를 Trend (ts_decompose$trend), Seasonal (ts_decompose$seasonal), Random (ts_decompose$random) 의 각 요소별로 나누어서 시각화해볼 수도 있습니다. 

 



# change the plot size
options(repr.plot.width=12, repr.plot.height=5)


# Trend
plot(as.ts(ts_decompose$trend))


# Seasonality
plot(as.ts(ts_decompose$seasonal))
# Random (Irregular factor)
plot(as.ts(ts_decompose$random))


 

 

 

많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~

'를 꾹 눌러주세요. :-)

 

 

 

728x90
반응형
Posted by Rfriend
,

Lag, Lead window function은 시계열 데이터를 처리할 때 많이 사용하는 매우 유용한 함수입니다. 


이번 포스팅에서는 PostgreSQL, Python (pandas), R (dplyr) 을 이용해서 그룹별로 행을 하나씩 내리기, 올리기 (lag or lead a row by group using PostgreSQL, Python, R) 하는 방법을 소개하겠습니다. 





  1. PostgreSQL로 그룹별로 특정 칼럼의 행을 하나씩 내리기, 올리기 

    (lag, lead a row by group using PostgreSQL lag(), lead() window function)


연월일(dt), 그룹ID(id), 측정값(val) 의 세 개 칼럼을 가진 시계열 데이터의 테이블을 PostgreSQL DB에 만들어보겠습니다. 



DROP TABLE IF EXISTS ts;

CREATE TABLE ts (

    dt date not null

    , id text not null  

    , val numeric not null

);


INSERT INTO ts VALUES 

  ('2019-12-01', 'a', 5)

, ('2019-12-02', 'a', 6)

, ('2019-12-03', 'a', 7)

, ('2019-12-04', 'a', 8)

, ('2019-12-01', 'b', 13)

, ('2019-12-02', 'b', 14)

, ('2019-12-03', 'b', 15)

, ('2019-12-04', 'b', 16);


SELECT * FROM ts ORDER BY id, dt;




PostgreSQL 의 LAG(value, offset, default), LEAD(value, offset, default) Window function을 이용해서 그룹ID('id') 별로 측정값('val')의 행을 하나씩 내리기(lag), 올리기(lead) 해보겠습니다. 행을 내리거나 올린 후에 빈 셀의 값은 'NULL'로 지정해주었습니다. 


LAG(), LEAD() 함수를 사용할 때 그룹ID('id')별로 년월일('dt') 을 기준으로 내림차순 정렬(OVER(PARTITIO BY id ORDER BY dt)) 을 해줍니다. 



-- lead() windows function

SELECT 

    *

    , LAG(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lag_1

    , LEAD(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lead_2

FROM ts;

 



lag(), lead() 함수를 사용해서 lag_1, lead_2 라는 새로운 칼럼을 추가한 'ts_lag_lead' 라는 이름의 테이블을 만들어보겠습니다. 



DROP TABLE IF EXISTS ts_lag_lead;

CREATE TABLE ts_lag_lead AS (

SELECT 

    *

    , LAG(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lag_1

    , LEAD(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lead_2

FROM ts

);


SELECT * FROM ts_lag_lead ORDER BY id, dt;

 





  2. Python pandas 로 DataFrame 내 그룹별 특정 칼럼의 행을 하나씩 내리기, 올리기 

     (shift a row by group using Python pandas library)


위에서 PostgreSQL의 lag(), lead() window function과 똑같은 작업을 Python pandas 를 가지고 수행해보겠습니다. 


먼저 dt, id, val의 칼럼을 가진 pandas DataFrame 시계열 데이터를 만들어보겠습니다. 



import pandas as pd


ts = pd.DataFrame({'dt': ['2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04', 

                          '2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04'], 

                  'id': ['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b'], 

                  'val': [5, 6, 7, 8, 13, 14, 15, 16]})


ts

dtidval
02019-12-01a5
12019-12-02a6
22019-12-03a7
32019-12-04a8
42019-12-01b13
52019-12-02b14
62019-12-03b15
72019-12-04b16

 



shift() 함수를 쓰기 전에 sort_values() 함수로 정렬을 해주는데요, lag 는 내림차순 정렬, lead는 오름차순 정렬임에 주의해야 합니다. (PostgreSQL, R 대비 Python이 좀 불편하긴 하네요 -,-;)


(a) lagsort_values() 함수를 이용해서 년월일('dt')를 기준으로 내림차순 정렬 (ascending=True) 한 후, 'id' 그룹별로 'val' 값을 하나씩 내려기 groupby('id')['val'].shift(1)


(b) lead: sort_values() 함수를 이용해서 년월일('dt')를 기준으로 오름차순 정렬 (ascending=False) 한 후, 'id' 그룹별로 'val' 값을 하나씩 올리기 groupby('id')['val].shift(1)



# lag a row by group 'id'

ts['val_lag_1'] =  ts.sort_values(by='dt', ascending=True).groupby('id')['val'].shift(1)


# lead a row by group 'id'

ts['val_lead_1'] = ts.sort_values(by='dt', ascending=False).groupby('id')['val'].shift(1)


ts.sort_values(by=['id', 'dt'])

dtidvalval_lag_1val_lead_1
02019-12-01a5NaN6.0
12019-12-02a65.07.0
22019-12-03a76.08.0
32019-12-04a87.0NaN
42019-12-01b13NaN14.0
52019-12-02b1413.015.0
62019-12-03b1514.016.0
72019-12-04b1615.0NaN

 





  3. R dplyr 로 dataframe 내 그룹별 특정 칼럼의 행을 하나씩 내리기, 올리기 

     (lag, lead a row by group using R dplyr library)



위에서 PostgreSQL의 lag(), lead() window function과 똑같은 작업을 R dplyr library를 가지고 수행해보겠습니다. 


먼저 dt, id, val의 칼럼을 가진 R DataFrame 시계열 데이터를 만들어보겠습니다. 



#install.packages("dplyr")

library(dplyr)


dt <- c(rep(c('2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04'), 2))

id <- c(rep('a', 4), rep('b', 4)) 

val <- c(5, 6, 7, 8, 13, 14, 15, 16)


ts <- data.frame(dt, id, val)

ts

A data.frame: 8 × 3
dtidval
<fct><fct><dbl>
2019-12-01a5
2019-12-02a6
2019-12-03a7
2019-12-04a8
2019-12-01b13
2019-12-02b14
2019-12-03b15
2019-12-04b16

 



R은 Postgresql 처럼 lag(), lead() window function을 가지고 있고 dplyr library의 chain operator를 써서 arrange() 함수로 'dt' 기준 내림차순 정렬하고, group_by(id)를 써서 그룹ID('id')별로 lag(), lead()를 무척 편리하게 적용해서 새로운 변수를 생성(mutate)할 수 있습니다. 



ts <- ts %>% 

    arrange(dt) %>%

    group_by(id)  %>% 

    mutate(val_lag_1 = lag(val, 1), 

          val_lead_1 = lead(val, 1))

 


arrange(ts, id, dt)

A grouped_df: 8 × 5
dtidvalval_lag_1val_lead_1
<fct><fct><dbl><dbl><dbl>
2019-12-01a5NA6
2019-12-02a657
2019-12-03a768
2019-12-04a87NA
2019-12-01b13NA14
2019-12-02b141315
2019-12-03b151416
2019-12-04b1615NA




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요.



728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 STEFANY COXE, STEPHEN G. WEST, AND LEONA S. AIKEN, 2009, "The Analysis of Count Data: A Gentle Introduction to Poisson Regression and Its Alternatives" 논문을 요약해서 한글로 번역한 내용을 소개하겠습니다. 

 

빈도(count data) 데이터에 대한 분석을 할 때 어떤 분포를 가정하고, 어떤 회귀모형을 적합하는 것이 좋을지 분석 방향을 잡는데 도움이 될 것입니다. 

 

 

 

 


  왜 포아송 회귀모형이 필요한가(Why Poisson Regression)?

 

특정 시간 동안 발생한 사건의 건수에 대한 도수 자료(count data, 음수가 아닌 정수)를 목표변수, 특히 평균이 10 미만일 경우, 최소제곱법(Ordinary Least Squares) 회귀모형을 적합하면 표준 오차와 유의수준이 편향되는 문제가 발생한다. OLS 회귀모형은 오차의 조건부 정규성(conditional normality), 등분산성(homoscedasiticity), 독립성(independence)을 가정한다. 하지만 평균이 작은 도수 데이터는 0 이상의 정수만 있고, 작은 계급에 많은 관측치가 몰려있으며, 우측으로 길게 치우친 분포를 띠어 회귀모형의 가정을 위배하고 음수값도 결과값으로 반환하는 문제가 있다. 종속변수가 도수 데이터이면서 오차가 정규분포, 등분산이 아닌 경우 포아송 회귀모형을 사용한다.

 

 

  • 포아송 회귀모형은 무엇인가? (Overview of Poisson Regression)

포아송 회귀모형은 일반화선형모형(Generalized Linear Model, GLM)의 한 종류로서,

 

-       Random component: Poisson Distribution, 

 

 

-       Systematic component: Mixed linear relationship

 

-       Link function: Log Link (ln)

 

 

 

표현 예

 

해석

 이며, X1을 한 단위 증가시켰을 때 효과를 보면 

  이므로, 다른 변수들을 통제했을 때 X1이 한단위 증가하면 도수의 추정치

 만큼 승법적으로 변화(multiplicative changes)한다.

 

 

선형회귀모형이 최소제곱법(OLS)로 모수를 추정한다면, 포아송 회귀모형은 최대가능도추정(MLE, Maximum Likelihood Estimation)을 통해 모수를 추정한다.

 

 

선형회귀모형은 모델에 의해 설명되는 잔차 제곱합 SS (Sum of Squares)를 총 잔차제곱합(TSS, Total Sum of Squares)로 나눈 값인  R^2로 모델 적합도(Goodness of Fit of the model) 평가하지만, 포아송 회귀모형은 설명변수의 추가에 따른 이탈도의 감소 비율(proportional reduction in deviance)로 표현되는 pseudo-R^2를 사용하여 완전모델(perfect model)에 얼마나 가깝게 적합이 되었는지를 평가한다

 

 

선형회귀모형은 F-test를 통해 모델의 유의성을 검정하지만, 포아송 회귀모형은 base modelcomplete model 간의 이탈도 차이에 대한 Chi-square test를 통해 모델 유의성을 검정한다 

선형회귀모형은 모형 타당성 평가 (Assessing model adequacy)를 위해 잔차도(residual plot)를 그려서 확인해보는 반면, 포아송 회귀모형은 이탈도 잔차도(deviance residuals vs. predicted values)를 그려보거나 실측치와 예측치를 비교해보는 것이다.

 

 

* [참고] 포아송 분포 : https://rfriend.tistory.com/101

 

 

  과대산포(Over-dispersion) 문제는 무엇이며어떻게 해결할 수 있나?

 

포아송분포는 평균과 분산이 동일(equidispersion)하다고 가정하므로, 분산이 평균보다 큰 과대산포 (Overdispersion) 시에는 포아송 회귀모형을 사용할 수 없다. 과대산포 시에 적용할 수 있는 대안 모델로는 과대산포 포아송 회귀(Overdispersed Poisson Regression), 음이항회귀(Negative Binomial Regression)가 있다.

 

 

-   Overdispersed Poisson Regression: 과대산포를 모델링하기 위해 두번째 모수인 조건부 분산 overdispersion scaling parameter

를 추정하며, 과대산포 포아송 회귀모형의 오차 분포는 

 를 가진다.

 

 

-   Negative Binomial Regression Models: 음이항회귀모형은 평균에는 영향이 없으면서 과대산포를 유발하는, 설명이 되지 않는 추가적인 가변성이 있다고 가정한다. 음이항 회귀모형은 똑같은 설명변수 값을 가지는 관측치가 다른 평균 모수를 가지고 포아송 회귀모형에 적합될 수 있도록 해준다. 오차 함수는 포아송 분포와 감마분포의 혼합 분포이다.

 

 

과대산포 검정은 가능도비 검정(Likelihood ratio test) 또는 score 검정을 사용한다.

 

모델 선택을 위해서는 AIC(Akaike Information Criterion), BIC(Bayesian Information Criterion) 을 사용하며, AIC, BIC의 값이 작은 모델을 선택하는데, 샘플 크기가 작을 때는 간단한 모델을 선택하고, 샘플 크기가 클 때는 좀더 복잡한 모델을 선택한다.

 

 

 

  과대영(Excess Zeros) 문제는 무엇이며어떻게 해결할 수 있나?

 

특정 평균을 모수로 가지는 포아송 분포에서 나타나는 ‘0’보다 많은 ‘0’을 가지는 샘플의 경우 과대영(excess zeros)이 발생했다고 말한다. 과대영은 (a) 특정 행동을 나타내지 않은 그룹의 포함 (: 비흡연자 그룹), (b) 특정 행동을 나타내지만 샘플링 계획에서 제외 (: 병원 입원한 사람만 연구) 하는 경우에 발생할 수 있다. 과대영일 경우 ZIP 회귀모형을 대안으로 사용할 수 있다.

 

-   Zero-Inflated Poisson(ZIP) Regression: 두 부분으로 나뉘는데, (1) 로지스틱 회귀모형을 사용하여 항상 ‘0’인 집단(structural zeros) 또는 항상 ‘0’은 아니지만 조사 시점에 ‘0’이라고 응답한 집단에 속할 확률을 구하고, (2) 포아송 회귀 또는 음이항 회귀모형으로 structural zeros를 포함하지 않은 나머지 관측치를 추정하는 모델을 적합한다.

 

 


  음이항분포(Negative Binomial Distribution) 

 

* Reference : https://en.wikipedia.org/wiki/Negative_binomial_distribution


 
확률 이론과 통계에서, 음 이항 분포는 특정 (무작위) 실패 횟수 (r로 표시)가 발생하기 전에 독립적이고 동일하게 분산 된 베르누이 시행 순서에서 성공 횟수의 이산 확률 분포이다. 예를 들어 1을 실패로, 1이 아닌 모든 것을 성공으로 정의하고 1이 세 번째로 나타날 때까지 반복적으로 던지면 (r = 3 번의 실패), 출현한 1이 아닌 수(non-1)의 확률 분포는 음 이항 분포가 된다.

 

 파스칼 분포 (Blaise Pascal 이후) Polya 분포 (George Pólya의 경우)는 음 이항 분포의 특수한 경우이다. 엔지니어, 기후 학자 및 다른 사람들 사이의 협약은 정수 값 정지 시간 매개 변수 r의 경우 "음 이항" 또는 "파스칼"을 사용하고, 실수 값의 경우에는 "Polya"를 사용한다. 토네이도 발생과 같은 "전염성 있는" 이산 이벤트가 발생하는 경우 Polya 분포를 사용하여 Poisson과 달리 평균 및 분산을 다르게 하여 Poisson 분포보다 더 정확한 모델을 제공 할 수 있다. "전염성 있는" 사건은 양의 공분산으로 인해 사건이 독립적이었던 경우보다 양의 상관관계가 있어 Poisson 분포에 의한 분산보다 더 큰 분산을 일으킨다. Poisson 분포는 평균과 분산이 동일하다고 가정하기 때문에 평균보다 분산이 큰 분포(Overdispersed)의 경우 적절하지 않으며, 2개의 모수를 가져서 과대산포분포를 적합할 수 있는 음 이항 분포를 사용할 수 있다.

 

 일련의 독립적 인 베르누이 재판이 있다고 가정하자. 따라서 각 임상 시험은 "성공" "실패"라는 두 가지 결과를 나타낸다. 각 시험에서 성공 확률은 p이고 실패 확률은 (1 - p)이다. 미리 정의 된 실패 횟수 r이 발생할 때까지의 서열을 관찰한다. 그러면 우리가 본 성공의 무작위 수 X는 음 이항 (또는 파스칼) 분포를 가질 것이다

 

 



 R을 활용한 Poisson model, Negative binomial model, Zero-Inflated model

 

Poisson model, Negative binomial model, Zero-Inflated model 의 유형, 분포, 모수 추정 방법, R패키지와 함수를 표로 정리하면 아래 [1]을 참고하세요. 
 

 

 

[1] Overview of count regression models

 

유형
분포
추정
방법
Description
R package & function
GLM
Poisson
ML
Poisson Regression: classical GLM, estimated by maximum likelihood(ML)
{stats} 패키지의 glm()
NB
ML
NB Regression: exteded GLM, estimated by ML including additional shape parameter
{stats} 패키지의 glm(), {MASS} 패키지의 glm.nb() 함수
Zero-augmented
Poisson
ML
Zero-Inflated Poission(ZIP)
{pscl} 패키지의 zeroinfl() 함수

 

각 회귀모형별로 R 패키지와 함수의 활용 방법은 아래와 같다.

-       (1) R을 활용한 Poisson model : glm() in the {stats} package

 



glm(formula, data, subset, na.action, weights, offset,


    family = poisson, start = NULL, control = glm.control(...),


    model = TRUE, y = TRUE, x = FALSE, ...)


 

 

-       (2) R을 활용한 Negative binomial model
   :
glm() in the {stats} package

 



glm(formula, data, subset, na.action, weights, offset,


    family = negative.binomial, start = NULL, control = glm.control(...),


    model = TRUE, y = TRUE, x = FALSE, ...)


 

 

   : glm.nb() in the {MASS} package

 



glm.nb(formula, data, weights, subset, na.action,


       start = NULL, etastart, mustart,


       control = glm.control(...), method = "glm.fit",


       model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, ...,


       init.theta, link = log)


 

-       (3) R을 활용한 Zero-inflated regression model: zeroinfl() in the {pscl} package

 



zeroinfl(formula, data, subset, na.action, weights, offset,


dist = "poisson", link = "logit", control = zeroinfl.control(...),


model = TRUE, y = TRUE, x = FALSE, ...)


 

 

많은 도움이 되었기를 바랍니다. 

 

 

728x90
반응형
Posted by Rfriend
,