R, Python 분석과 프로그래밍의 친구 (by R Friend)

[R data.table] data.table 의 구문 DT[i, j, by] 에서 행 subset 과 열 select 하고 계산하기

R 분석과 프로그래밍/R 데이터 전처리 2020. 10. 3. 23:49

지난번 포스팅에서는 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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] 여러개의 변수를 처리하는 data.table의 특수 부호 .SD, .SDcols (3)	2020.10.10
[R data.table] data.table의 by 구문으로 그룹별 집계하기 (Aggregations by Group) (4)	2020.10.04
[R data.table] data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문의 비교 (0)	2020.10.02
[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기 (0)	2020.09.27
[R data.table] data.table 은 무엇이고, 왜 data.table 인가? (0)	2020.08.24

Posted by Rfriend

,

[R data.table] data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문의 비교

R 분석과 프로그래밍/R 데이터 전처리 2020. 10. 2. 23:13

지난번 포스팅에서는 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)`

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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table의 by 구문으로 그룹별 집계하기 (Aggregations by Group) (4)	2020.10.04
[R data.table] data.table 의 구문 DT[i, j, by] 에서 행 subset 과 열 select 하고 계산하기 (4)	2020.10.03
[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기 (0)	2020.09.27
[R data.table] data.table 은 무엇이고, 왜 data.table 인가? (0)	2020.08.24
[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수 (6)	2020.02.23

Posted by Rfriend

,

[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기

R 분석과 프로그래밍/R 데이터 전처리 2020. 9. 27. 22:04

지난번 포스팅에서는 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 % 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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table 의 구문 DT[i, j, by] 에서 행 subset 과 열 select 하고 계산하기 (4)	2020.10.03
[R data.table] data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문의 비교 (0)	2020.10.02
[R data.table] data.table 은 무엇이고, 왜 data.table 인가? (0)	2020.08.24
[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수 (6)	2020.02.23
[R] 여러개의 변수를 가진 DataFrame을 무작위 층화 샘플링으로 Train, Test set 분할하고 표준화하기 (0)	2020.01.19

Posted by Rfriend

,

[R data.table] data.table 은 무엇이고, 왜 data.table 인가?

R 분석과 프로그래밍/R 데이터 전처리 2020. 8. 24. 23:00

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 구문 소개 포스팅 링크

데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기
: https://rfriend.tistory.com/563
data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문 비교

: https://rfriend.tistory.com/565
data.table의 구문 DT[i, j, by]에서 행 subset과 열 select하고 계산하기

: https://rfriend.tistory.com/566
data.table의 by 구문으로 그룹별 집계하기 (Aggregation by Group)

: https://rfriend.tistory.com/567
여러개의 변수를 처리하는 data.table의 특수 부호 .SD, .SDcols

: https://rfriend.tistory.com/568
data.table의 키와 빠른 이진 탐색 기반의 부분집합 선택

: https://rfriend.tistory.com/569
data.table에서 Key를 칼럼 j와 그룹 by와 함께 사용하기

: https://rfriend.tistory.com/570
data.table에서 Key 값의 매칭되는 모든 행(mult="all"), 첫번째 행(mult="first"), 마지막 행(mult="last"), 값이 존재하는 행만(nomatch=NULL) 가져오기
: https://rfriend.tistory.com/571
data.table을 참조하여 := 연산자로 data.table의 칼럼 추가, 갱신, 삭제하기
(add/updata/delete column in R data.table)
: https://rfriend.tistory.com/573
data.table을 참조하여 := 연산자로 얕은 복사 시의 부작용, copy() 함수로 깊은 복사하는 방법
: https://rfriend.tistory.com/574
data.table의 효율적인 데이터 재구조화 : melt(), dcast()
: https://rfriend.tistory.com/575

패턴이 있는 여러개의 칼럼을 동시에 길게 녹이고(melt), 옆으로 넓게 주조(cast) 하기
: rfriend.tistory.com/576
그룹별 관측치 개수 별로 Data.Table을 구분해서 생성하기
: rfriend.tistory.com/607

.SDcols 로 일부 칼럼 가져오기 (Column subsetting using .SDcols)
: rfriend.tistory.com/608

회귀 모델의 오른쪽 부분(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

R data.table의 자동 인덱싱 (Auto indexing)
: rfriend.tistory.com/616

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

'를 꾹 눌러주세요.

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

728x90

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문의 비교 (0)	2020.10.02
[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기 (0)	2020.09.27
[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수 (6)	2020.02.23
[R] 여러개의 변수를 가진 DataFrame을 무작위 층화 샘플링으로 Train, Test set 분할하고 표준화하기 (0)	2020.01.19
[R] 장비 On, Off 상태로 가동 시간을 구하고, 장비 별 평균 가동 시간 구하기 (16)	2019.10.24

Posted by Rfriend

,

[R] 자기상관계수 (Autocorrelation Coefficients), 자기상관그림(Autocorrelation Plot)

R 분석과 프로그래밍/R 통계분석 2020. 3. 22. 23:58

자기상관계수(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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 통계분석' 카테고리의 다른 글

[R] 빈도 데이터 분석을 위한 포아송 회귀모형과 대안 회귀모형에 대한 소개 (9)	2019.10.23
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (4/4): 로지스틱 회귀모형 적합 및 모델평가, 해석 (19)	2018.12.18
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (3/4): 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석 (8)	2018.12.17
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (2/4): 탐색적 데이터 분석 및 전처리 (32)	2018.12.16
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (1/4): WDBC Data 소개 (4)	2018.12.16

Posted by Rfriend

,

[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수

R 분석과 프로그래밍/R 데이터 전처리 2020. 2. 23. 20:01

이번 포스팅에서는 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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기 (0)	2020.09.27
[R data.table] data.table 은 무엇이고, 왜 data.table 인가? (0)	2020.08.24
[R] 여러개의 변수를 가진 DataFrame을 무작위 층화 샘플링으로 Train, Test set 분할하고 표준화하기 (0)	2020.01.19
[R] 장비 On, Off 상태로 가동 시간을 구하고, 장비 별 평균 가동 시간 구하기 (16)	2019.10.24
[R] R 패키지 함수의 소스 코드를 볼 수 있는 방법 (how to see function's source codes in R package) (0)	2019.10.22

Posted by Rfriend

,

[R] 여러개의 변수를 가진 DataFrame을 무작위 층화 샘플링으로 Train, Test set 분할하고 표준화하기

R 분석과 프로그래밍/R 데이터 전처리 2020. 1. 19. 13:08

이번 포스팅에서는 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
Price	Horsepower	RPM	Length	Type	Origin
<dbl>	<int>	<int>	<int>	<fct>	<fct>
15.9	140	6300	177	Small	non-USA
33.9	200	5500	195	Midsize	non-USA
29.1	172	5500	180	Compact	non-USA
37.7	172	5500	193	Midsize	non-USA
30.0	208	5700	186	Midsize	non-USA
15.7	110	5200	189	Midsize	USA

table(X$Origin)

USA non-USA 48 45

y <- Cars93$MPG.highway

y

31
25
26
26
30
31
28
25
27
25
25
36
34
28
29
23
20
26
25
28
28
26
33
29
27
21
27
24
33
28
33
30
27
29
30
20
30
26
50
36
31
46
31
33
29
34
27
22
24
23
26
26
37
36
34
24
25
29
25
26
26
33
24
33
30
23
26
31
31
23
28
30
41
31
28
27
28
26
38
37
30
30
43
37
32
29
22
33
21
30
25
28
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
	Price	Horsepower	RPM	Length	Type	Origin
	<dbl>	<int>	<int>	<int>	<fct>	<fct>
44	8.0	81	5500	168	Small	non-USA
2	33.9	200	5500	195	Midsize	non-USA
39	8.4	55	5700	151	Small	non-USA
40	12.5	90	5400	164	Sporty	non-USA
3	29.1	172	5500	180	Compact	non-USA
53	8.3	82	5000	164	Small	non-USA
45	10.0	124	6000	172	Small	non-USA
90	20.0	134	5800	180	Compact	non-USA
42	12.1	102	5900	173	Small	non-USA
16	16.3	170	4800	178	Van	USA
7	20.8	170	4800	200	Large	USA
11	40.1	295	6000	204	Midsize	USA
73	9.0	74	5600	177	Small	USA
12	13.4	110	5200	182	Compact	USA
8	23.7	180	4000	216	Large	USA
23	9.2	92	6000	174	Small	USA
17	16.6	165	4000	194	Van	USA
74	11.1	110	5200	181	Compact	USA
14	15.1	160	4600	193	Sporty	USA

table(X$Origin)

[Out]:
    USA non-USA 
     48      45

table(X_test$Origin)

[Out]:
    USA non-USA 
     10       9

y_test

[Out]: 33
25
50
36
26
37
29
30
46
23
28
25
41
36
25
33
20
31
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
	Price	Horsepower	RPM	Length
1	-0.44621989	-0.1188381	1.7177896	-0.46498935
4	1.80188107	0.5065143	0.3726422	0.64947906
5	1.00782706	1.2100357	0.7089291	0.16189913
41	-0.04403669	0.2720072	0.8770725	-0.60429791
43	-0.28122166	-0.1188381	0.5407856	0.09224485
46	-1.05465089	-1.0568667	0.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
	Price	Horsepower	RPM	Length
44	-1.2608987	-1.2718315	0.3726422	-1.0918778
2	1.4100103	1.0536976	0.3726422	0.7887876
39	-1.2196491	-1.7799303	0.7089291	-2.2760005
40	-0.7968411	-1.0959512	0.2044988	-1.3704949
3	0.9150156	0.5065143	0.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
	Price	Horsepower	RPM	Length
1	0.15596330	0.3248945	0.9259259	0.4615385
4	0.55596330	0.4599156	0.6296296	0.6666667
5	0.41467890	0.6118143	0.7037037	0.5769231
41	0.22752294	0.4092827	0.7407407	0.4358974
43	0.18532110	0.3248945	0.6666667	0.5641026
46	0.04770642	0.1223629	0.6481481	0.3205128

head(X_test_num_scaled)

A matrix: 6 × 4 of type dbl
	Price	Horsepower	RPM	Length
44	0.01100917	0.07594937	0.6296296	0.3461538
2	0.48623853	0.57805907	0.6296296	0.6923077
39	0.01834862	-0.03375527	0.7037037	0.1282051
40	0.09357798	0.11392405	0.5925926	0.2948718
3	0.39816514	0.45991561	0.6296296	0.5000000
53	0.01651376	0.08016878	0.4444444	0.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.Midsize	Type.Small	Origin.non-USA
0	1	1
1	0	1
0	0	1
1	0	1
1	0	1
1	0	0

head(X_test_cat_dummy)

A matrix: 6 × 6 of type dbl
	Type.Large	Type.Midsize	Type.Small	Type.Sporty	Origin.non-USA
75	0	0	0	1	0
76	0	1	0	0	0
77	1	0	0	0	0
78	0	0	0	0	1
79	0	0	1	0	0
80	0	0	1	0	1

(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)

[Out]: Small
Midsize
Compact
Large
Sporty
Van

unique(X_train_cat$Origin)

[Out]: non-USA
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_small	type_midsize	type_compact	type_large	type_sporty	origin_nonusa
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	0	0	0	0	1
0	1	0	0	0	1
0	0	1	0	0	1
0	1	0	0	0	1
0	1	0	0	0	1
0	1	0	0	0	0

# 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_small	type_midsize	type_compact	type_large	type_sporty	origin_nonusa
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
0	0	0	0	1	0
0	1	0	0	0	0
0	0	0	1	0	0
0	0	1	0	0	1
1	0	0	0	0	0
1	0	0	0	0	1

(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
Price	Horsepower	RPM	Length	type_small	type_midsize	type_compact	type_large	type_sporty	origin_nonusa
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
0.1559633	0.3469388	0.9259259	0.4615385	1	0	0	0	0	1
0.4862385	0.5918367	0.6296296	0.6923077	0	1	0	0	0	1
0.3981651	0.4775510	0.6296296	0.5000000	0	0	1	0	0	1
0.5559633	0.4775510	0.6296296	0.6666667	0	1	0	0	0	1
0.4146789	0.6244898	0.7037037	0.5769231	0	1	0	0	0	1
0.1522936	0.2244898	0.5185185	0.6153846	0	1	0	0	0	0

# 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
	Price	Horsepower	RPM	Length	type_small	type_midsize	type_compact	type_large	type_sporty	origin_nonusa
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
75	0.18899083	0.42857143	0.2962963	0.70512821	0	0	0	0	1	0
76	0.20366972	0.59183673	0.4444444	0.69230769	0	1	0	0	0	0
77	0.31192661	0.46938776	0.3703704	0.46153846	0	0	0	1	0	0
78	0.39082569	0.34693878	0.8148148	0.55128205	0	0	1	0	0	1
79	0.06788991	0.12244898	0.4444444	0.44871795	1	0	0	0	0	0
80	0.01834862	0.07346939	0.6666667	0.06410256	1	0	0	0	0	1

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

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

728x90

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table 은 무엇이고, 왜 data.table 인가? (0)	2020.08.24
[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수 (6)	2020.02.23
[R] 장비 On, Off 상태로 가동 시간을 구하고, 장비 별 평균 가동 시간 구하기 (16)	2019.10.24
[R] R 패키지 함수의 소스 코드를 볼 수 있는 방법 (how to see function's source codes in R package) (0)	2019.10.22
[R] 시계열 정수의 순차 3개 묶음 패턴 별 개수를 구하고 내림차순 정렬하기 (16)	2019.10.15

Posted by Rfriend

,

[Python 시계열 자료 분석] 시계열 분해 (Time series Decomposition)

Python 분석과 프로그래밍/Python 통계분석 2020. 1. 2. 23:46

지난번 포스팅에서는 시계열 자료의 구성 요인 (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

저작자표시 비영리 변경금지

'Python 분석과 프로그래밍 > Python 통계분석' 카테고리의 다른 글

[Python] DataFrame에서 여러개의 변수에 대해 일원분산분석 검정하기 (ANOVA test for multiple numeric variables in pandas DataFrame) (0)	2021.05.08
[Python] 샘플 크기가 다른 2개 이상 그룹간 일원분산분석 (one-way ANOVA with different sized samples) (0)	2021.05.07
[Python] 선형회귀모형, 로지스틱 회귀모형에 대한 각 관측치 별 변수별 기여도(민감도) 분석 (Sensitivity analysis of linear regression & Logistic regression per each variables and each observations) (1)	2020.01.12
[Python 시계열 자료 분석] 시계열 패턴별 지수 평활법 (exponential smoothing by time series patterns) (5)	2020.01.05
[Python 시계열 자료 분석] 시계열 구성 요인 (Time series component factors): 추세(trend), 순환(cycle), 계절(seasonal), 불규칙(irregular) 요인 (2)	2020.01.01

Posted by Rfriend

,

[Python, R, PostgreSQL] 그룹별로 행을 내리기, 올리기 (Lag, Lead a row by Group)

Python 분석과 프로그래밍/Python 데이터 전처리 2019. 12. 9. 00:55

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

	dt	id	val
0	2019-12-01	a	5
1	2019-12-02	a	6
2	2019-12-03	a	7
3	2019-12-04	a	8
4	2019-12-01	b	13
5	2019-12-02	b	14
6	2019-12-03	b	15
7	2019-12-04	b	16

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

(a) lag: sort_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'])

	dt	id	val	val_lag_1	val_lead_1
0	2019-12-01	a	5	NaN	6.0
1	2019-12-02	a	6	5.0	7.0
2	2019-12-03	a	7	6.0	8.0
3	2019-12-04	a	8	7.0	NaN
4	2019-12-01	b	13	NaN	14.0
5	2019-12-02	b	14	13.0	15.0
6	2019-12-03	b	15	14.0	16.0
7	2019-12-04	b	16	15.0	NaN

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
dt	id	val
<fct>	<fct>	<dbl>
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

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
dt	id	val	val_lag_1	val_lead_1
<fct>	<fct>	<dbl>	<dbl>	<dbl>
2019-12-01	a	5	NA	6
2019-12-02	a	6	5	7
2019-12-03	a	7	6	8
2019-12-04	a	8	7	NA
2019-12-01	b	13	NA	14
2019-12-02	b	14	13	15
2019-12-03	b	15	14	16
2019-12-04	b	16	15	NA

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

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

728x90

저작자표시 비영리 변경금지

'Python 분석과 프로그래밍 > Python 데이터 전처리' 카테고리의 다른 글

[Python datetime] datetime 모듈의 데이터 유형 (data types in pandas datetime module) (0)	2019.12.21
[Python pandas] resample() 메소드로 시계열 데이터를 10분 단위 구간별로 집계/요약 하기 (8)	2019.12.15
[Python] pandas DataFrame: ValueError: If using all scalar values, you must pass an index 에러 해결 방법 (0)	2019.09.15
[Python] 웹에서 XML 포맷 데이터를 Python으로 읽어와서 DataFrame으로 만들기 (18)	2019.09.04
[Python] 웹으로 부터 JSON 포맷 데이터 읽어와서 pandas DataFrame으로 만들기 (0)	2019.08.31

Posted by Rfriend

,

[R] 빈도 데이터 분석을 위한 포아송 회귀모형과 대안 회귀모형에 대한 소개

R 분석과 프로그래밍/R 통계분석 2019. 10. 23. 23:37

이번 포스팅에서는 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)

표현 예:

해석 :

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

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

만큼 승법적으로 변화(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 model과 complete 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]을 참고하세요.

* Reference: http://cran.at.r-project.org/web/packages/pscl/vignettes/countreg.pdf

[표1] Overview of count regression models

유형	분포	추정 방법	Description	R package & function
GLM	Poisson	ML	Poisson Regression: classical GLM, estimated by maximum likelihood(ML)	{stats} 패키지의 glm()
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

저작자표시 비영리 변경금지

'R 분석과 프로그래밍 > R 통계분석' 카테고리의 다른 글

[R] 자기상관계수 (Autocorrelation Coefficients), 자기상관그림(Autocorrelation Plot) (10)	2020.03.22
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (4/4): 로지스틱 회귀모형 적합 및 모델평가, 해석 (19)	2018.12.18
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (3/4): 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석 (8)	2018.12.17
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (2/4): 탐색적 데이터 분석 및 전처리 (32)	2018.12.16
[R] 로지스틱 회귀분석을 통한 유방암 예측(분류) (1/4): WDBC Data 소개 (4)	2018.12.16

Posted by Rfriend

,

R, Python 분석과 프로그래밍의 친구 (by R Friend)

'r'에 해당되는 글 267건

[R data.table] data.table 의 구문 DT[i, j, by] 에서 행 subset 과 열 select 하고 계산하기

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table의 기본 구문 DT[i, j, by][order] 과 Base R, dplyr, SQL 구문의 비교

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] 데이터를 읽어와서 data.table 만들기, data.frame을 data.table로 변환하기

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R data.table] data.table 은 무엇이고, 왜 data.table 인가?

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R] 자기상관계수 (Autocorrelation Coefficients), 자기상관그림(Autocorrelation Plot)

'R 분석과 프로그래밍 > R 통계분석' 카테고리의 다른 글

[R dplyr] 여러개의 if, else if 조건절을 벡터화해서 처리해주는 case_when() 함수

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[R] 여러개의 변수를 가진 DataFrame을 무작위 층화 샘플링으로 Train, Test set 분할하고 표준화하기

'R 분석과 프로그래밍 > R 데이터 전처리' 카테고리의 다른 글

[Python 시계열 자료 분석] 시계열 분해 (Time series Decomposition)

'Python 분석과 프로그래밍 > Python 통계분석' 카테고리의 다른 글

[Python, R, PostgreSQL] 그룹별로 행을 내리기, 올리기 (Lag, Lead a row by Group)

'Python 분석과 프로그래밍 > Python 데이터 전처리' 카테고리의 다른 글

[R] 빈도 데이터 분석을 위한 포아송 회귀모형과 대안 회귀모형에 대한 소개

'R 분석과 프로그래밍 > R 통계분석' 카테고리의 다른 글

카테고리

태그목록

티스토리툴바