이번 포스팅에서는 카이제곱 적합도 검정(Chi-squared goodness of fit test)를 활용하여 관측치가 특정 분포를 따르는지 아닌지를 검정하는 방법을 소개하겠습니다.

 

다음과 같은 문제가 있다고 해봅시다.

 

 

[문제]  시간이 지남에 따라 일어난 어떤 특정 사건의 발생건수에 대한 결과가 다음과 같이 주어졌다고 하자. 그 때 특정 사건의 발생 회수가 포아송(λ) 분포를 따르는지에 대한 검정을 유의수준 5%에서 실시하시오.

 

 

 발생건수

0회

1회

2회 

3회 

4회 

5회 

6회 

7회 

8회 

9회 이상 

 관찰회수

6

20

36 

50 

52 

40 

25 

12 

 

 

(풀이)

특정 관측치들이 어떤 이론적 분포를 따르고 있는지를 검정하는 데에는 적합도 검정(Chi-squared goodness of fit test) 를 이용합니다.

 

검정해야 할 귀무가설(H0)와 대립가설(H1)을 세워보면,

  • 귀무가설 : 특정 사건의 발생 회수가 포아송 분포(Poisson(λ))를 따른다.
                     (즉, 관측값의 도수와 가정한 분포로부터의 기대관측도수가 같다)
  • 대립가설 : 특정 사건의 발생 회수가 포아송 분포를 따르지 않는다. (Not )

 

이며,

카이제곱 검정 통계량은  입니다. (이 문제에서는 i=0, 1, ..., 9)

이때 k 는 범주(class)의 개수를 나타내며, 는 도수분포의 각 구간에 있는 관측도수,

는 도수분포의 각 구간의 기대도수를 나타냅니다.

 

 

주어진 데이터로부터 가정된 포아송 분포의 평균(λ)의 추정치을 구하면, =3.844가 됩니다.

 

 = (0*6+1*20+2*36+3*50*4*52+5*40+6*25+7*12+8*4+9*5)/250 = 3.844

 

R 로 위의 관측치 데이터셋을 입력하고 가정된 포아송 분포의 평균 추정치를 구해보도록 하겠습니다.

 

 

> #===========================================
> # Chi-squared Goodness-of-fit Test
> # Poisson Distribution, 5% significant level
> #===========================================
> rm(list=ls()) # clear all
>
> # Class of Events
> classnum <- c(0:9)
>
> # Observed Frequency
> obsnum <- c(6, 20, 36, 50, 52, 40, 25, 12, 4, 5)
>
> # Sum of Observed Frequency
> obssum <- sum(obsnum); obssum
[1] 250
>
> # Calculate Lambda
> lamb <- sum(classnum*obsnum) / sum(obsnum); lamb
[1] 3.844 

 

 

 

포아송분포의 모수 λ에 대한 추정치로 3.844를 사용하여 각 class 별 가설적인 확률 (i=0, 1, ..., 9)를 구해보겠습니다.

 

에 i=0, 1, ..., 9 를 대입해서 풀어보면 아래의 표와 같이 가설적인 확률 를 구할 수 있습니다.

 

R 의 dpois() 함수를 사용하여 확률 구해봤는데요, sum(hyp_prob)=0.994로서 전체 합이 1이 안되고 0.994 이네요.

 

> # Calculate the hypothesised probabilities
> hyp_prob = round(dpois(classnum, lamb), 3); hyp_prob
 [1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.011
> sum(hyp_prob) # not 1
[1] 0.994

 

확률의 공리적 정의 상 표본공간의 P(S)=1 이어야 하므로 모든 class별 확률을 전부 더했을 때 1이 되도록 보정을 해줄 필요가 있습니다. 만약 보정을 안한 상태에서 바로 카이제곱 검정을 실시하면 아래와 같이 "확률값들의 총계는 반드시 1 이어야 합니다"와 같은 에러 메시지가 뜹니다.

 

> # Chi-squared goodness-of-fit test
> chisq.gof.test <- chisq.test(obsnum, p=hyp_prob) # error
Error in chisq.test(obsnum, p = hyp_prob) : 
  확률값들의 총계는 반드시 1 이어야 합니다 

 

 

R의 indexing을 사용하여 가장 마지막인 '9회 이상' class의 확률 를 조정하여 전체 class별 확률의 합이 1이 되도록 조정 해보겠습니다.

 

 

> # Change the last hypothesised probability to make sum(hyp_prob)=1
> hyp_prob[length(hyp_prob)] <- 1-sum(hyp_prob[1:(length(hyp_prob)-1)]); hyp_prob
[1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017
> sum(hyp_prob) # 1
[1] 1 

 

 

 

발생확률 Pi를 구했으므로 전체 발생회수 합인 250과 각 class별 가정된 포아송분포로 부터의 확률을 곱한 250*로 부터 기대도수(Expected frequency)를 구할 수 있습니다.

 

 발생건수

0회

1회

2회 

3회 

4회 

5회 

6회 

7회 

8회 

9회 이상 

 계

 발생확률

0.021

0.082

0.158

0.203

0.195

0.150

0.096

0.053

0.025

0.017

 1

기대도수

(Ei=Pi*250)

5.25

20.50

39.50

50.75

48.75

37.50

24.00

13.25

6.25

4.25

250

 관측도수

(Oi)

6

20

36

50

52

40

25

12

4

5

250

 

 

이제 카이제곱 적합도 검정(Chi-squared goodness of fit test)을 할 준비가 되었네요. 가정된 포아송분포(λ=3.844)로 부터의 기대도수(Expected frequency)와 250개 샘플 데이터로 부터 실제 관측한 도수(Observed frequency) 간에 차이가 있다고 볼 수 있는 것인지, 관측치가 포아송분포를 따른다고 볼 수 있는 것인지 가설을 5% 유의수준 하에서 검정해보겠습니다.

 

R의 chisq.test() 함수를 사용하여 카이제곱 적합도 검정을 수행한 결과는 아래와 같습니다.

 

 

> # Chi-squared goodness-of-fit test
> chisq.gof.test <- chisq.test(obsnum, p=hyp_prob)
Warning message:
In chisq.test(obsnum, p = hyp_prob) :
카이제곱 approximation은 정확하지 않을수도 있습니다
> chisq.gof.test

Chi-squared test for given probabilities

data:  obsnum
X-squared = 1.9258, df = 9, p-value = 0.9926 

 

 

검정통계량인 X-squared 통계량이 1.9258이고 P-value 가 0.9926 이므로 유의수준 5% 하에서 관측값이 포아송 분포를 따른다는 귀무가설을 채택(Accept Null Hyperthesis)합니다.

 

 

R로 카이제곱 적합도 검정을 한 결과를 chisq.gof.test 라는 이름의 객체로 저장을 했는데요, names() 함수로 확인을 해보면 chisq.gof.test 객체 리스트 안에 여러 분석 결과 통계량들이 들어있음을 알 수 있습니다.

 

> names(chisq.gof.test)
[1] "statistic" "parameter" "p.value"   "method"    "data.name"
[6] "observed"  "expected"  "residuals" "stdres"  

 

 

이들 중에서 index을 해서 관측도수(Observed frequency)와 기대도수(Expected frequency)를 확인해보겠습니다.

 

 

> chisq.gof.test$observed
[1]  6 20 36 50 52 40 25 12  4  5
> chisq.gof.test$expected
[1]  5.25 20.50 39.50 50.75 48.75 37.50 24.00 13.25  6.25  4.25
> round(hyp_prob, 3)
[1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017 

 

 

 

카이제곱 검정을 할 때 각 class의 기대도수가 너무 작을 경우 검정이 유효하지 않을 수가 있습니다. 위에 R분석 결과에서 빨간색으로 "카이제곱 approximation은 정확하지 않을 수도 있습니다"라는 warning message가 나온 이유가 '9회 이상' class의 기대도수가 4.25로서 작기 때문입니다. 최소 기대도수에 대해서는 일반적인 합의는 없으며, 보통 3, 4, 5개가 종종 사용되며, 3~5개 이하일 경우 인접 class와 병합해서 재분석을 하기도 합니다.

 

이번 문제에서는 최소 기대도수가 4.25이고, 또 P-value가 0.9926으로서 매우 높은 수준으로 귀무가설을 채택하였으므로 '9회 이상' class를 '8회'와 합쳐서 재분석하는 것은 필요없어 보입니다. (병합해서 재분석해도 결과는 동일하게 귀무가설 채택함)

 

마지막으로, 각 class별 관측도수와 기대도수를 그래프로 그려서 시각화로 비교를 해보겠습니다.

 

> # Plot of Poisson Distribution: Observed vs. Expected
> # Observed Frequency
> plot(c(0:9), chisq.gof.test$observed, 
+      main="Poisson Distribution: Observed vs. Expected", 
+      type='b', 
+      pch=0,
+      col='blue',
+      xlab="Number of Events",
+      ylab="Frequency",
+      ylim=c(0,55))
> 
> # Dual Y-axes plot (secondary Y-axis)
> par(new=T)
> 
> # Expected frequency
> plot(c(0:9), chisq.gof.test$expected, 
+      type='b', 
+      pch=1,
+      col='red', 
+      xlab="", 
+      ylab="",
+      ylim=c(0,55))
> 
> legend(x=6.5, y=50, 
+        c("Observed", "Expected"), 
+        pch=c(0,1), 
+        col=c('blue', 'red'))

 

 

 

 

 

Poisson(λ=3.844) 로부터의 각 class별 기대도수(Expected frequency)와 250개 샘플로 부터의 관측도수(Observed frequency)가 매우 유사한 분포를 띠고 있음을 눈으로도 확인할 수 있습니다.

 

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

 

Posted by R Friend R_Friend

이번 포스팅에서는 R ggplot2 로 y축이 2개 있는 이중 축 그래프 그리는 방법을 소개하겠습니다.

 

엑셀로는 이중 축 그래프 그리기가 그리 어렵지 않은데요, R의 ggplot 도 sec.axis 을 사용하면 두번째 y축을 추가할 수 가 있습니다.

 

먼저 예제로 사용할 데이터셋을 만들어보겠습니다.

 

 

#=========================================
# Dual y-axes plot using ggplot2 sec.axis
#=========================================

rm(list=ls()) # clear all

# make a DataFrame
data_val <- matrix(
  c(-4.2922960, 0.00000000000389755, 0.00000000268579,
  -3.9102123, 0.00000000000939471, 0.00000003082191,
  -3.5281286, 0.00000000002264455, 0.00000027610486,
  -3.1460449, 0.00000000005458189, 0.00000193070482,
  -2.7639612, 0.00000000013156254, 0.00001053865139,
  -2.3818776, 0.00000000031711433, 0.00004490363779,
  -1.9997939, 0.00000000076436268, 0.00014935003987,
  -1.6177102, 0.00000000184239690, 0.00038775418083,
  -1.2356265, 0.00000000444085801, 0.00078584150863,
  -0.8535428, 0.00000001070411038, 0.00124319933681,
  -0.4714591, 0.00000002580086522, 0.00153523161114,
  -0.0893754, 0.00000006218962756, 0.00147990680737,
  0.2927083, 0.00000014989999897, 0.00111358189390,
  0.6747920, 0.00000036131440584, 0.00065408965954,
  1.0568757, 0.00000087090114242, 0.00029990225777,
  1.4389594, 0.00000209919260108, 0.00010733701573,
  1.8210431, 0.00000505982314747, 0.00002998793957,
  2.2031268, 0.00001219600184343, 0.00000653989939,
  2.5852105, 0.00002939662278856, 0.00000111332724,
  3.5852105, 0.00029392734366929, 0.00000000000000,
  4.5852105, 0.00293538878457633, 0.00000000000000,
  5.5852105, 0.02896916461589227, 0.00000000000000,
  6.5852105, 0.25470155892952129, 0.00000000000000,
  7.5852105, 0.94711869936516890, 0.00000000000000,
  8.5852105, 0.99999999999982903, 0.00000000000000,
  9.5852105, 1.00000000000000000, 0.0000000000000),
  nrow=26,
  byrow=TRUE)

# convert a matrix into a dataframe
df <- data.frame(data_val)

# column name
colnames(df) <- c("x", "y1", "y2")

 

> head(df)
          x           y1           y2
1 -4.292296 3.897550e-12 2.685790e-09
2 -3.910212 9.394710e-12 3.082191e-08
3 -3.528129 2.264455e-11 2.761049e-07
4 -3.146045 5.458189e-11 1.930705e-06
5 -2.763961 1.315625e-10 1.053865e-05
6 -2.381878 3.171143e-10 4.490364e-05

 

 

 

두 개 y축의 요약 정보를 보니 최대값의 차이가 매우 크다는 것을 알 수 있습니다. y1 축과 y2 축의 비율을 보니 약 651배 차이가 나네요. 이런 상태에서 그래프를 그리면 y2 축의 그래프가 바닥에 쫙 붙어서 그려지므로 두 축 간의 패턴, 특징, 구조를 비교하기가 어렵게 됩니다. 따라서 y2 축의 데이터를 y1축과 비교하기 쉽도록 y2에 'y1과 y2의 최대값의 비율(max_ratio)'를 곱해서 그래프를 그려보겠습니다.

 

 

> # max ratio b/w 2 axes
> summary(df)
       x                 y1                  y2          
Min.   :-4.2923   Min.   :0.0000000   Min.   :0.000e+00 
1st Qu.:-1.9043   1st Qu.:0.0000000   1st Qu.:7.000e-10 
Median : 0.4838   Median :0.0000003   Median :8.539e-06 
Mean   : 1.1492   Mean   :0.1243873   Mean   :3.020e-04 
3rd Qu.: 3.3352   3rd Qu.:0.0002278   3rd Qu.:3.658e-04 
Max.   : 9.5852   Max.   :1.0000000   Max.   :1.535e-03 
> max_ratio <- max(df$y1)/max(df$y2); max_ratio
[1] 651.367

 

 

 

y1 데이터로는 선 그래프를 그리고, y2 데이터로는 막대그래프를 그려보겠습니다.

이때 y2 에 해당하는 두번째 축의 막대그래프를 그릴 때는

 - (1) max(y1)/max(y2) 로 계산한 두 축의 최대값의 비율인 651을 y2에 곱했으며

       (geom_bar(aes(y = y2*max_ratio) 

 - (2) scale_y_continuous(sec.axis = sec_axis(~.*maxratio, name="y2") 를 사용해서 오른쪽의 두번째 축을 추가하였습니다.

 

 

> # Dual y-axes plot using ggplot2 sec.axis
> library(ggplot2)
> 
> g <- ggplot(df, aes(x = x))
>   g <- g + geom_line(aes(y = y1), colour = "red", size = 2)
>   
>   # adding the relative result5.data.A, 
> # transformed to match roughly the range of the B
> g <- g + geom_bar(aes(y = y2*max_ratio),
> fill = "blue",
> stat = "identity")
> > # adding secondary axis > g <- g + scale_y_continuous(sec.axis = sec_axis(~.*max_ratio, name="y2*651")) > > g

 

 

 

해석하기에 편리하도록 데이터의 소스가 무엇인지를 나타내는 텍스트와 화살표를 추가해보겠습니다.

 

 

> # adding text > g <- g + annotate("text", x = -2, y = 0.75, colour="black", label="y2*651") > g <- g + annotate("text", x = 5, y = 0.5, colour="black", label="y1") > > # adding arrow > library(grid) > g <- g + annotate("segment", x = -1.8, y = 0.75, + xend = -1, yend = 0.6, colour="blue", arrow=arrow()) > g <- g + annotate("segment", x = 5.2, y = 0.5, + xend = 6.7, yend = 0.4, colour="red", arrow=arrow()) > g

 

 

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

 

Posted by R Friend R_Friend

이번 포스팅에서는 


(1) R Environment 상에 존재하는 여러개의 DataFrame 중에서 이름이 특정 조건을 만족하는 DataFrame을 선별하여 

 : ls(pattern = "xx")


(2) 하나의 List 로 묶는 방법, 

 : mget()


(3) List로 부터 특정 DataFrame을 Indexing 하는 방법

 : list[[1]], or  list[["name"]]


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



먼저, 이번 포스팅의 예제의 특성 상 Environment 와 Console을 깨끗하게 청소하는 것부터 시작하겠습니다.  rm(list=ls()) 로 Environment에 있는 모든 객체를 삭제하고, cat("\014")로 Console을 백지상태로 만들어주었습니다. 



# To make an environment, console clear 

rm(list=ls()) # clearing environment

cat("\014") # clearing Console




객체 이름에 'data_'를 공통으로 포함한 3개의 DataFrame 예제와, 객체 이름에 'file_'을 포함한 3개의 DataFrame 예제를 만들어보겠습니다. 



# To make several sample DataFrames

data_1 <- data.frame(var1 = c(1, 2), var2 = c(3, 4))

data_2 <- data.frame(var1 = c('a', 'b'), var2 = c('c', 'd'))

data_3 <- data.frame(var1 = c(TRUE, TRUE), var2 = c(FALSE, TRUE))

file_1 <- data.frame(var1 = c(1, 2), var2 = c(3, 4))

file_2 <- data.frame(var1 = c('a', 'b'), var2 = c('c', 'd'))

file_3 <- data.frame(var1 = c(TRUE, TRUE), var2 = c(FALSE, TRUE))

  [DataFrame samples]

 



Environment에 생성된 객체를 확인해보려면 ls() 함수를 사용하면 됩니다. 그리고 특정 문자열을 포함한 객체만을 선별해서 찾아보려면 ls(pattern = "xx") 처럼 pattern = "xx" 를 추가하면 됩니다. 



> # To list up all objects in an environment

> ls()

[1] "data_1" "data_2" "data_3" "file_1" "file_2" "file_3"

> # To list up objects which have a matching 'pattern' in a DataFrame name

> ls(pattern = "data_")

[1] "data_1" "data_2" "data_3"

> ls(pattern = "file_")

[1] "file_1" "file_2" "file_3"

 




다음으로, mget() 함수를 사용하여 ls(pattern = "data_"), ls(pattern = "file_"로 각각 선별한 DataFrame 객체들을 2개의 List로 각각 묶어서 생성해보겠습니다.  



> # To combine all DataFrame into a List using mget()

> list_df_data <- mget(ls(pattern = "data_"))

> list_df_data

$data_1

  var1 var2

1    1    3

2    2    4


$data_2

  var1 var2

1    a    c

2    b    d


$data_3

  var1  var2

1 TRUE FALSE

2 TRUE  TRUE


> list_df_file <- mget(ls(pattern = "file_"))

> list_df_file

$file_1

  var1 var2

1    1    3

2    2    4


$file_2

  var1 var2

1    a    c

2    b    d


$file_3

  var1  var2

1 TRUE FALSE

2 TRUE  TRUE






마지막으로, List에 묶인 DataFrame을 [[ ]] 을 사용하여 위치(숫자) 혹은 이름으로 Indexing 하는 방법을 아래에 소개합니다. 



> # To index one of object in a List using [[ ]]

> list_df_data[[1]] # using 'number'

  var1 var2

1    1    3

2    2    4

> list_df_data["data_1"] # using 'name'

$data_1

  var1 var2

1    1    3

2    2    4 



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

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



Posted by R Friend R_Friend

이번 포스팅은 페이스북의 R User Group 에서 아이디어 토론 주제로 올라왔길레 한번 짜봤던 코드입니다.  


'0'과 '1'로 가지는 긴 벡터를 처음부터 끝까지 쭉 훓어가면서 


(1) '0'이 연속으로 3번 나오면 그 구간을 기준으로 나누어서 

=> (2) 나누어진 구간을 새로운 벡터로 생성하기


입니다.  


'0'이 연속으로 나오는 회수는 분석가가 필요로 하는 회수로 지정할 수 있도록 매개변수(argument)로 지정해서 프로그래밍해보겠습니다. 


간단한 예제 벡터로서, '0'과 '1'을 원소로 해서 30개의 무작위수(이항분포 랜덤 샘플링, 0이 80%, 1이 20%)로 구성된 벡터를 생성해보겠습니다. 


 

> rm(list=ls()) # clean all


# Sample vector


> set.seed(123) # for reproducibility

> vec_raw <- sample(c(0, 1), size=30, replace=TRUE, prob=(c(0.8, 0.2)))

> vec_raw

 [1] 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0

> vec <- vec_raw # copy for manipulation




R 프로그래밍을 하고자 하는 일의 아웃풋 이미지는 아래와 같습니다. 

'0'이 연속 세번 나오는 구간을 노란색으로 표시를 했습니다. 이 구간을 기준으로 원래의 벡터를 구분해서 나눈 후에(split), 각각을 새로운 벡터로 생성하여 분리해주는 프로그램을 짜보는 것이 이번 포스팅의 주제입니다. 




아래에 wihle() 반복문, if(), else if() 조건문, assign() 함수를 사용해서 위의 요건을 수행하는 코드를 짜보았습니다. 원래 벡터에서 제일 왼쪽에 있는 숫자가 '0'이 아니면 하나씩 빼서 빈 벡터에 차곡차곡 쌓아 놓구요, '0'이 3개 나란히 있는 '000'이 나타나면 그동안 쌓아두었던 벡터를 잘라내서 다른 이름의 벡터로 저장하고, '000'도 그 다음 이름의 벡터로 저장하면서 원래 벡터에서 '000'을 빼놓도록 하는 반복문입니다. assign() 함수는 저장하는 객체의 이름에 paste()를 사용해서 변화하는 숫자를 붙여주고자 할 때 사용하는 함수인데요, 자세한 사용법 예제는 여기(=> http://rfriend.tistory.com/108)를 참고하세요. 


아래 프로그램에서 'bin_range'에 할당하는 숫자를 변경하면 원하는 회수만큼 '0'이 반복되었을 때 구간을 분할하도록 조정할 수 있습니다.  가령, '0'이 연속 10번 나오면 분할하고 싶으면 bin_range <- 10  이라고 할당해주면 됩니다. 



##---------------------------------------------

# Vector Bin Split by successive '0's criteria

##---------------------------------------------


# Setting

vec_tmp <- c() # null vector

bin_range <- 3 # successive '0's criterion

vec_idx <- 1

vec_num <- 1


# Zero_Bin_Splitter

while (length(vec) > 0){

  if (sum(vec[1:bin_range], na.rm=T) != 0){

    vec_tmp[vec_idx] <- vec[1]

    vec <- vec[-1]

    vec_idx <- vec_idx + 1

  } else if (is.null(vec_tmp)){

    assign(paste0("vec_split_", vec_num), vec[1:bin_range])

    vec <- vec[-(1:bin_range)]

    vec_idx <- 1 # initialization

    vec_num <- vec_num + 1

  } else {

    assign(paste0("vec_split_", vec_num), vec_tmp)

    vec_tmp <- c() # initialization

    vec_num <- vec_num + 1

    assign(paste0("vec_split_", vec_num), vec[1:bin_range])

    vec <- vec[-(1:bin_range)]

    vec_idx <- 1 # initialization

    vec_num <- vec_num + 1

  }

}


# delete temp vector and arguments

rm(vec, vec_tmp, bin_range, vec_idx, vec_num)




아래는 위의 프로그램을 실행시켰을 때의 결과입니다.  원래 의도했던대로 잘 수행되었네요. 




=====================================================

아래의 코드는 페이스북 R User Group의 회원이신 June Young Lee 님께서 data.table 라이브러리를 사용해서 짜신 것입니다.  코드도 간결할 뿐만 아니라, 위에 제가 짠 코드보다 월등히 빠른 실행속도를 자랑합니다. 대용량 데이터를 빠른 속도로 조작, 처리하려면 data.table 만한게 없지요. 아래 코드는 저도 공부하려고 옮겨 놓았습니다. 


똑같은 일을 하는데 있어서도 누가 어떤 로직으로 무슨 패키지, 함수를 사용해서 짜느냐에 따라서 복잡도나 연산속도가 이렇게 크게 차이가 날 수 있구나 하는 좋은 예제가 될 것 같습니다. 이런게 프로그래밍의 묘미이겠지요? ^^  June Young Lee 님께 엄지척! ^^b



# by June Young Lee


library(data.table)


# 0이 세번 연속으로 들어 있는 아무 벡터 생성

vec <- c(0,0,0,3,10,2,3,0,4,0,0,0,1,2,50,4,0,0,32,1,0,0,0,1,1)


# data.table함수인 rleid이용하여 연속인 행들을 구분하는 idx1 컬럼 생성 

dt <- data.table(v1=vec, idx1=rleid(vec))


# idx1 기준으로 갯수를 세어 N 컬럼 생성 

dt[,N:=.N, by=idx1]


# 0이 세번 연속으로 나온 그룹(=N컬럼변수가 3인)을 idx2로 체크하기 

dt[v1==0&N==3,idx2:=1L]


# split 하기 위한 기준 벡터를 마찬가지로 rleid함수를 이용하여 idx3로 생성


dt[,idx3:=rleid(idx2)]


# data.table의 split 함수 적용하고, lapply 적용하여 v1 칼럼만 추출 

res <- lapply(split(dt, by="idx3"), function(x) x[,v1])

res

# $`1`

# [1] 0 0 0

# $`2`

# [1] 3 10 2 3 0 4

# $`3`

# [1] 0 0 0

# $`4`

# [1] 1 2 50 4 0 0 32 1

# $`5`

# [1] 0 0 0

# $`6`

# [1] 1 1


# 0이 3번 연속으로 나온 놈들을 남겨두도록 작성한 이유는, 

# 위치확인 & 확인용입니다. 

# 필요없으면, 아래 정도의 코드를 추가하면 되겠네요.


res2 <- res[!unlist(lapply(res, function(x) length(x)==3&sum(x)==0))]

res2

# $`2`

# [1] 3 10 2 3 0 4

# $`4`

# [1] 1 2 50 4 0 0 32 1

# $`6`

# [1] 1 1

 


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


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



Posted by R Friend R_Friend

페이스북의 KRUG (Korean R User Group) 에서 2017.11.18일 주말 퀴즈로 아래의 시각화를 푸는 문제가 있어서 재미로 풀어보았습니다. 


처음엔 금방 코드 짤 것으로 예상했는데요, 출제자분께서 중간 중간 장애물을 심어놓으셔서 제 예상보다 좀 더 걸렸네요. 


문제 풀면 KRUG 운영자분께서 스타벅스 기프티콘 주신다고 하는데요, 기대되네요. ^___^



# KRUG's plot quiz, as of 18th NOV. 2017


library(ggplot2)

#install.packages("dplyr")

library(dplyr)

str(mtcars)


# making 'model' variable

mtcars$model <- rownames(mtcars)


# ordering by cyl, hp

mtcars_ord <- arrange(mtcars[,c('model', 'cyl', 'hp')], cyl, desc(hp))


# marking an inflection point within cluster : threshold >= 40

mtcars_cyl <- c(4, 6, 8)

mtcars_ord_2 <- data.frame()


for (i in 1:length(mtcars_cyl)) {

  

  mtcars_tmp <- subset(mtcars_ord, subset = (cyl == mtcars_cyl[i]))

  

  for (j in 1:nrow(mtcars_tmp)) {

    if (j != nrow(mtcars_tmp) & mtcars_tmp$hp[j] - mtcars_tmp$hp[j+1] >= 40) {

      mtcars_tmp$hp_outlier[j] = '1_outlier'

    } else {

      mtcars_tmp$hp_outlier[j] = '2_normal'

    }

  }

  

  mtcars_ord_2 <- rbind(mtcars_ord_2, mtcars_tmp)

  rm(mtcars_tmp)

}


# converting cyl variable type from numeric to factor

mtcars_ord_2$cyl_cd <- paste0("cyl :", mtcars_ord_2$cyl)


model_order <- mtcars_ord_2$model[order(mtcars_ord_2$cyl_cd, 

                                        mtcars_ord_2$hp, 

                                        decreasing = FALSE)]


mtcars_ord_2$model <- factor(mtcars_ord_2$model, levels = model_order)



# drawing cleveland dot plot

ggplot(mtcars_ord_2, aes(x = hp, y = model)) +

  geom_point(size = 2, aes(colour = hp_outlier)) +

  scale_colour_manual(values = c("red", "black")) + 

  theme_bw() +

  facet_grid(. ~ cyl_cd, scales = "free_y", space = "free_y") +

  xlim(0, max(mtcars_ord_2$hp)) +

  geom_hline(yintercept = nrow(mtcars_ord_2[mtcars_ord_2$cyl == 4,]) + 0.5, 

             colour = "black", linetype = "dashed", size = 0.5) +

  geom_hline(yintercept = nrow(mtcars_ord_2) - 

               nrow(mtcars_ord_2[mtcars_ord_2$cyl == 8,]) + 0.5, 

             colour = "black", linetype = "dashed", size = 0.5) +

  theme(legend.position = 'none')

 




R Korea - KRUG(Korean R User Group) 에 문제 출제해주신 정우준님께서 나중에 정답지 올려주신 코드도 아래에 같이 공유합니다. 제가 짠 코드보다 한결 간결하네요. 



library(data.table)
library(dplyr)
library(ggplot2)

mtcars %>%
mutate(car.name=rownames(.)) %>%
arrange(cyl, hp) %>%
mutate(order.key=1:n()) -> data

data %>%
ggplot(aes(x=hp, y=reorder(car.name, order.key))) +
geom_point(
colour=case_when(
data$car.name %in% c('Ferrari Dino','Maserati Bora') ~ 'red', 
TRUE ~ 'black')) +
geom_hline(yintercept = 11.5, linetype='dashed') +
geom_hline(yintercept = 18.5, linetype='dashed') +
facet_wrap(~ cyl, labeller = label_both) +
scale_x_continuous(limits=c(0,max(data$hp))) +
theme_bw() +
theme(axis.title.y=element_blank())

 






Posted by R Friend R_Friend

문자열이나 숫자를 특정 형식으로 길이를 지정해주면 데이터를 출력했을 때 깨끗하게 정리가 되어 보이기 때문에 가독성이 좋아집니다. 


혹은 데이터가 특정 형식(format)으로 DB에 이미 지정이 되어 있어서 데이터 간 병합이나 join을 하기 위해 특정 형식으로 데이터를 표준화 해주어야 할 경우가 있습니다. 


이번 포스팅에서는 {base} package의 sprintf() 함수를 사용해서 


 - (1) 문자열을 매개변수 width 길이로 만들고, 빈 자리는 '0'으로 채우기 : sprintf("%05d", var)


 - (2) 소수점 숫자(numeric)의 자리수를 지정해주기 : sprintf(".5f", var)


하는 방법에 대해서 알아보겠습니다. 


이번 포스팅의 함수 sprintf()는 데이터 전처리할 때 종종 사용하는 편이예요. 



 (1) 문자열을 특정 길이로 만들고, 빈 자리수만큼 '0'을 채우기 : sprintf("%05d", var)





1자리, 2자리, 3자리, 4자리를 가진 데이터를 가지고 예제로 사용할 간단한 DataFrame을 만들어보겠습니다. 



> # making a sample DataFrame

> df <- data.frame(var1 = c(1, 11, 111, 1111))

> df

  var1

1    1

2   11

3  111

4 1111

 




위의 예제 데이터셋을 가지고, 칼럼 var1의 데이터를 


- '1 자리수를 가진 문자열'로 만들되, '1자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '2 자리수를 가진 문자열'로 만들되, '2자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '3 자리수를 가진 문자열'로 만들되, '3자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '4 자리수를 가진 문자열'로 만들되, '4자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '5 자리수를 가진 문자열'로 만들되, '5자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'


를 해보겠습니다. 


만약 매개변수 자리수보다 데이터의 길이가 더 크다면 '0'이 채워지지는 않습니다.  아래의 예제의 결과를 View(df)로 해서 보면 원래의 변수 var1과 sprintf() 함수를 사용해서 만든 var1_01d 변수의 데이터 출력 형식이 다른 것을 알 수 있습니다. 


그리고 class 함수로 데이터 형식을 살펴보니 원래 변수 var1은 숫자형(numeric)이지만 sprintf() 함수로 만든 새로운 변수는 요인형(factor)의 문자열로 바뀌어 있음을 알 수 있습니다. 



> #-------------------------

> # (1) sprintf(%03d, var) : Format number as fixed width, with leading zeros

> df <- transform(df, 

+                 var1_01d = sprintf("%01d", var1), 

+                 var1_02d = sprintf("%02d", var1), 

+                 var1_03d = sprintf("%03d", var1), 

+                 var1_04d = sprintf("%04d", var1), 

+                 var1_05d = sprintf("%05d", var1))

> df

  var1 var1_01d var1_02d var1_03d var1_04d var1_05d

1    1        1       01      001     0001    00001

2   11       11       11      011     0011    00011

3  111      111      111      111     0111    00111

4 1111     1111     1111     1111     1111    01111



> View(df)



> sapply(df, class)

     var1  var1_01d  var1_02d  var1_03d  var1_04d  var1_05d

"numeric"  "factor"  "factor"  "factor"  "factor"  "factor"




 (2) 소수점 숫자(numeric)의 자리수를 지정해주기 : sprintf(".5f", var)


무리수인 자연상수 e의 소수점 10째 자리까지의 수를 대상으로 sprintf("%.5f", e) 함수를 사용해서 소수점의 자리수를 설정해보겠습니다. "%.숫자f"의 숫자 만큼 소수점을 표시해주는데요, 반올림을 해서 표시해줍니다. 아래의 예제를 보시면 금방 이해할 수 있을 것입니다. 



> #-------------------------

> # (3) sprintf("%.5f", x) : formatting decimal point, 

> e <- c(2.7182818284) # mathematical constant, the base of the natural logarithm

> sprintf("%.0f", e)

[1] "3"

> sprintf("%.1f", e)

[1] "2.7"

> sprintf("%.2f", e)

[1] "2.72"

> sprintf("%.3f", e)

[1] "2.718"

> sprintf("%.5f", e)

[1] "2.71828"

> sprintf("%.10f", e)

[1] "2.7182818284"

 




아래의 예시는 sprintf("%숫자.f, e)로 '숫자' 부분에 매개변수로 정수 부분의 자리수를 지정해주는 예시입니다. 소수점의 자리도 모두 포함해서 '숫자' 부분 매개변수만큼의 길이로 표시 형식을 맞추어줍니다. 



> e <- c(2.7182818284) # mathematical constant, the base of the natural logarithm

> sprintf("%1.1f", e)

[1] "2.7"

> sprintf("%2.1f", e)

[1] "2.7"

> sprintf("%3.1f", e)

[1] "2.7"

> sprintf("%5.1f", e)

[1] "  2.7"

> sprintf("%10.1f", e)

[1] "       2.7"

>  



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

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





Posted by R Friend R_Friend

이번 포스팅에서는 R ggplot2 패키지의 coord_fixed(ratio = number) 옵션을 사용해서 그래프의 크기 (가로, 세로 비율)를 조정하는 방법을 소개하겠습니다. 


1 ~ 4의 정수 좌표를 가지는 X 축과 

1 ~ 4의 정수 좌표를 가지는 Y 축을 가지는 

가상의 데이터를 사용해서 예를 들어보겠습니다. 


X축이 Y축이 1:1 비율인 그래프가 되겠습니다. 



[ 가상 데이터셋 생성 ]



#==========================

# Rssizing the plot

#==========================


install.packages("ggplot2")

library(ggplot2)


# X, Y coordinate and segments(A, B category) data set

my.df <- data.frame(XCoord = c(1, 1, 1, 2, 2, 2, 3, 3, 4, 4), 

                    YCoord = c(1, 2, 4, 1, 3, 4, 2, 4, 1, 4), 

                    Seg = c("A", "A", "A", "A", "B", "A", "B", "B", "B", "B"))





X와 Y좌표별로 Seg. 변수의 "A"는 검정색, "B"는 흰색으로 색깔을 지정해서 Heatmap을 그려보겠습니다. 

X좌표와 Y좌표가 1~4 범위를 동일하게 가지고 있으므로 크기에 대한 설정없이 디폴트 세팅으로 그래프를 그리면 아래 처럼 정사각형 1:1 비율로 그래프가 그려집니다. 



[ 그림 1 ] 원본 그래프 (Original Plot) : 가로축과 세로축이 1:1



# Original plot

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  ggtitle("Original Heatmap by X4*Y4 size")






[ 그림 2 ] 가로축과 세로축의 비율을 1:2 로 설정하기 : coord_fixed(ratio = 2)


[그림 1]에서 원본 이미지가 가로축과 세로축이 1:1 비율의 정사각형 그래프였는데요, 이를 가로:세로 비율을 1:2로 세로가 가로의 2배 비율인 그래프(가로:세로 = 1:2)로 바꾸어 주려면 coord_fixed(ratio = 2) 를 설정해주면 됩니다. 



# Resized plot using coord_fixed(ratio = number)

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  coord_fixed(ratio = 2) +

  ggtitle("Heatmap : Resized X:Y axis size with 1:2 ratio")







[ 그림 3 ] 가로축과 세로축의 비율을 2:1 로 설정하기 : coord_fixed(ratio = 0.5)


가로와 세로축 비율이 1:1인 원본 이미지 [그림 1] 을 가로가 세로축의 2배인 그래프 (가로:세로 = 2:1)로 바꾸고 싶다면 coord_fixed(ratio = 0.5) 로 설정해주면 됩니다. 



# Resized plot using coord_fixed(ratio = number)

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  coord_fixed(ratio = 0.5) +

  ggtitle("Heatmap : Resized X:Y axis size with 2:1 ratio")









아래는 EBImage 패키지의 resize() 함수를 사용해서 png 이미지 파일로 출력할 때 이미지 크기를 (a) 특정 가로, 세로 크기로 설정해주는 방법과 (b) 비율로 설정해주는 방법입니다. 

R code는 stackoverflow 의 답변 중에서 aoles 님께서 달아놓은 것인데요, 코드 그대로 인용해서 소개합니다. ( * R code 출처 : https://stackoverflow.com/questions/35786744/resizing-image-in-r )


#==============
# Image resize using EBImage package

# installing EBImage package
source("http://bioconductor.org/biocLite.R")
biocLite("EBImage")


# resizing image using EBImage package's resize() function

library("EBImage")

x <- readImage(system.file("images", "sample-color.png", package="EBImage"))


# width and height of the original image
dim(x)[1:2]


# scale to a specific width and height
y <- resize(x, w = 200, h = 100)


# scale by 50%; the height is determined automatically so that
# the aspect ratio is preserved
y <- resize(x, dim(x)[1]/2)


# show the scaled image
display(y)


# extract the pixel array
z <- imageData(y)


# or
z <- as.array(y)

 



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

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



Posted by R Friend R_Friend

보통 외부 데이터 불러오기 할 때 {utils} package의 read.table() 함수를 사용하곤 합니다. 


크기가 작은 데이터라면 별 문제를 못느낄 텐데요, 만약 데이터 사이즈가 수 메가, 기가 단위의 큰 데이터라면 데이터를 불러들이는데 너무 오랜 시간이 걸려서 문제가 될 수 있습니다. 


대용량 데이터 처리에 아주 뛰어난 성능을 발휘하는 data.table 패키지의 fread() 함수를 사용하면 큰 용량의 외부 데이터도 빠르게 불러올 수 있습니다. 


(물론 R은 메모리에 데이터를 올려놓고 처리/분석을 하므로 하둡에서 말하는 수테라급의 대용량에는 필적을 못하구요, 분산병렬처리도 아니긴 합니다. 이 포스팅에서 말하는 대용량은 책보고 공부할 때 사용하는 수십, 수백개 row 를 가진 예제 데이터 대비 실전에서 사용하는 수십만, 수백만, 수천만 row 데이터를 말하는 것입니다. ^^;)


아래에 간단한 샘플 데이터를 만들어서 {utils} 패키지의 read.table() 함수와 {data.table} 패키지의 fread() 함수의 데이터 불러오는데 소요되는 시간을 비교해보았습니다. 



[ 외부 데이터 읽어오기 : {utils} 패키지 read.table() 함수 vs. {data.table} 패키지 fread(0 함수 ]




1. 샘플 데이터 만들기 : 1 백 만개 row, 변수 2개를 가지는 데이터 프레임



# generating large scaled data

my_data <- data.frame(var1 = rnorm(n = 1000000, mean = 0, sd = 1), 

                      var2 = rnorm(n = 1000000, mean = 2, sd = 3))

 



> str(my_data)

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

 $ var1: num  -0.556 1.787 0.498 -1.967 0.701 ...

 $ var2: num  1.669 0.597 4.452 1.405 6.936 ...

 



# exporting to text file

write.table(my_data, 

            "/Users/Desktop/R/my_data.txt",

            sep = "|",

            row.names = FALSE, 

            quote = FALSE)

 




2. {utils} 패키지의 read.table() 함수를 사용해서 my_data.txt 불러오기


system.time() 함수로 데이터를 불어오는데 소요된 시간을 재어보았더니 7.287초가 나왔습니다. 

(매번 할 때마다 소요 시간이 조금씩 차이가 날 수는 있습니다)



# reading text file : (1) read.table() of {utils} package

> system.time(my_data_readtable <- read.table("/Users/Desktop/R/my_data.txt",

+                                             sep = "|", 

+                                             header = TRUE, 

+                                             stringsAsFactors = FALSE))

   user  system elapsed 

  7.161   0.100   7.287

 




3. {data.table} 패키지의 fread() 함수를 사용해서 my_data.txt 불러오기


data.table 패키지는 기본 패키지가 아니므로 먼저 별도 설치(install.packages) 및 호출(library)이 필요합니다. 

(매번 할 때마다 소요 시간이 조금씩 차이가 날 수는 있습니다)



# reading text file : (2) fread() of {data.table} package

install.packages("data.table")

library(data.table)

 



system.time() 함수로 my_data를 불러오는데 걸린 시간을 재어봤더니 0.256초가 걸렸습니다. 



> system.time(my_data_fread <- fread("/Users/Desktop/R/my_data.txt", 

+                                    sep = "|", 

+                                    header = TRUE, 

+                                    stringsAsFactors = FALSE))

   user  system elapsed 

  0.242   0.014   0.256

 




1백만 행을 가진 데이터프레임을 읽어오는데 있어, 앞서 read.table() 함수가 7.287 초 걸렸던데 반해, fread() 함수는 0.256 초밖에 걸리지 않았습니다.  fread() 함수는 read.table() 함수를 사용했을 때 대비 약 96.5% 정도 속도가 더 적게 걸린 것입니다.  놀랍지요?!!! 



> 0.256/7.287

[1] 0.03513106

 



R 은 대용량 데이터에는 맥을 못춰라고 지레 평가절하하기 보다는 {data.table} 패키지의 fread() 함수로 대용량 데이터 불러오기 속도 문제를 공략해보시지요. 


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

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


Posted by R Friend R_Friend

(X축) 시간의 흐름에 따른 (Y축) 값의 추세, 변화 분석 및 탐색을 하는데 시계열 선 그래프(time series plot, line graph)를 많이 이용합니다. 


 이번 포스팅에서는 R ggplot2 패키지로 시계열 선그래프를 그리고, 거기에 세로선을 추가하는 작업을 해보겠습니다. 


ggplot2 에서 세로선을 추가할 때 geom_vline() 함수를 사용하는데요, 이게 시계열 데이터의 경우는 as.numeric() 함수를 사용해서 시계열 데이터를 숫자형 데이터로 변환을 해주어야 에러가 안나고 제대로 세로선이 그려집니다. 


이거 몰라서 한참을 구글링하면서 애 좀 먹었습니다. ^^;


간단한 시계열 데이터를 만들어서 예를 들어보겠습니다. 




> ##======================================================

> ## adding multiple vertical lines at time-series plot using R ggplot2

> ##======================================================

> # making time series data

> dt <- c("20170609100000", "20170609100100", "20170609100200", 

+         "20170609100300", "20170609100400", "20170609100500")

> val <- c(5.2, 3.4, 3.9, 6.3, 4.7, 5.6)

> dt_val <- data.frame(dt, val)

> dt_val <- transform(dt_val, 

+                     dt = as.POSIXct(dt, 

+                                     format = '%Y%m%d%H%M%S', 

+                                     origin = "1970-01-01", 

+                                     tz = "UTC"))

> dt_val

                   dt val

1 2017-06-09 10:00:00 5.2

2 2017-06-09 10:01:00 3.4

3 2017-06-09 10:02:00 3.9

4 2017-06-09 10:03:00 6.3

5 2017-06-09 10:04:00 4.7

6 2017-06-09 10:05:00 5.6

 




R ggplot2 패키지의 geom_line() 함수를 사용해서 시계열 선그래프를 그려보겠습니다. 



> # making time series plot

> library(ggplot2)

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size=1, color = "blue") + 

+   ggtitle("Time-series plot")

 




R ggplot2로 세로선을 추가할 때는 geom_vline(xintercept = x) 함수를 추가해주면 됩니다. 하지만 xintercept 에 들어가는 값이 날짜, 시간 포맷의 데이터일 경우 아래 처럼 에러가 납니다. 



> # To add vertical line at time series plot

> # Error in Ops.POSIXt((x - from[1]), diff(from)) : '/' not defined for "POSIXt" objects

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size = 1, color = "blue") +

+   geom_vline(xintercept = dt_val$dt[3], color =  "red", linetype = 2) +

+   ggtitle("Adding vertical line at time-series plot using geom_vline()")

Error in Ops.POSIXt((x - from[1]), diff(from)) : 

  '/' not defined for "POSIXt" objects

 




R ggplot2 시계열 선그래프에 X축이 날짜, 시간 포맷의 시계열 데이터인 경우 특정 날짜/시간에 세로선을 추가하기 위해서는 as.numeric(x) 함수를 사용해서 숫자형 데이터로 포맷을 바꾸어 주어야 합니다



> # Use as.numeric() function at xintercept

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size = 1, color = "blue") +

+   geom_vline(xintercept = as.numeric(dt_val$dt[3]), color = "red", linetype = 2) + # as.numeric() transformation

+   ggtitle("Adding vertical line at time-series plot using geom_vline() and as.numeric() transformation")

 







만약 복수의 세로선을 추가하고 싶다면 아래의 예제를 참고하세요. 만약 3번째와 5번째 x변수의 날짜/시간에 세로선을 추가하고 싶다면 dataset$variable[c(3, 5)] 처럼 indexing을 해서 xintercept 에 넣어주면 됩니다. 

(세로선 2개가 그려지기는 했는데요, 하단에 빨간색으로 "HOW_BACKTRACK environmental varialbe"이라는 경고메시지가 떴습니다. -_-; )



> # adding "Multiple" vertical lines

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size = 1, color = "blue") +

+   geom_vline(xintercept = as.numeric(dt_val$dt[c(3,5)]), color = "red", linetype = 2) +

+   ggtitle("Adding multiple vertical lines at time-series plot using geom_vline()")

HOW_BACKTRACE environmental variable.


 




이번에는 세로선을 그릴 기준 날짜/시간 데이터를 다른 데이터프레임에서 가져와야 하는 경우를 예로 들어보겠습니다. 먼저 세로선의 기준이 되는 xintercept 에 들어갈 날짜/시간 정보가 들어있는 data frame 을 만들어보죠. 



> # adding multiple vertical lines with another data frame

> dt_2 <- c("20170609100150", "20170609100430")

> val_2 <- c("yes", "yes")

> dt_val_2 <- data.frame(dt_2, val_2)

> dt_val_2 <- transform(dt_val_2, 

+                       dt_2 = as.POSIXct(dt_2, 

+                                         format = '%Y%m%d%H%M%S', 

+                                         origin = "1970-01-01", 

+                                         tz = "UTC"))

> dt_val_2

                 dt_2 val_2

1 2017-06-09 10:01:50   yes

2 2017-06-09 10:04:30   yes

 




R ggplot2 시계열 선그래프를 그린 원본 데이터프레임(아래 예제에서는 dt_val)과는 다른 데이터프레임(아래 예제에서는 dt_val_2)의 날짜/시간 데이터를 사용해서 복수의 세로선을 그려보겠습니다.  두 개의 세로선이 그려지기는 했는데요, "HOW_BACKTRACE environmental variable"이라는 빨간색 경고 메시지가 떴습니다. 그런데, 예전에는 에러 메시지 뜨면서 안그려졌었는데, 블로그 쓰려고 다시 해보니 그려지기는 하는군요. ^^; 



> # time series plot with multiple vertical lines from another data frame

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size = 1, color = "blue") +

+   geom_vline(xintercept = as.numeric(dt_val_2$dt_2), color = "red", linetype = 2) +

+   ggtitle("Time-seires plot with multiple vertical line from another dataframe")

HOW_BACKTRACE environmental variable.




위의 예제처럼 했는데 혹시 Error: Aesthetics must be either length 1 or the same as the data (6): xintercept 와 같은 에러 메시지가 뜨고 그래프가 안그려진다면 아래처럼 geom_vline(data = dataframe, xintercept = ... ) 처럼 데이터를 가져오는 데이터프레임을 명시해주면 문제가 해결됩니다.  이걸 몰라서 또 한참을 고민하고, 구글링하고, 참 애먹었던 적이 있습니다. -_-;



> # time series plot with multiple vertical lines from another data frame(2)

> ggplot(dt_val, aes(x = dt, y = val)) +

+   geom_line(size = 1, color = "blue") +

+   geom_vline(data = dt_val_2, xintercept = as.numeric(dt_val_2$dt_2), color = "red", linetype = 2) +

+   ggtitle("Time-series plot with multiple vertical lines from another dataframe 2")



 



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


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



Posted by R Friend R_Friend

지난번 포스팅에서는 범주형 특성 데이터, 텍스트 문서의 거리를 측정하는 지표 중에서 


 - 자카드 거리 (Jaccard distance)


 - 코사인 거리 (Cosine distance)


에 대해서 알아보았습니다. 


이번 포스팅에서는 두 문자열의 거리(distance between two strings of characters), 비유사도(dissimilarity)를 측정하는데 사용하는 편집 거리(edit distance), 혹은 다른 이름으로 Levenshtein metric 에 대해서 알아보겠습니다. 


편집거리(edit distance)는 데이터 항목이 놓인 순서(order)가 중요한 문자열(strings of characters, 예: 주소, 전화번화, 이름 스펠링)이나 서열(sequence, 예 : 염색체 염기서열)의 (비)유사도를 측정하는데 유용하게 사용할 수 있습니다. 


편집거리 (edit distance, Levenshtein metric) 는 두 문자열에서 하나의 문자열을 다른 문자열과 똑같게 만들기 위해서 최소로 필요로 하는 편집 회수(문자 추가, 제거, 위치 변경)를 계산합니다. 


아래에 간단한 예를 들어서 설명해보겠습니다. 



[ 편집 거리 예시 (example of edit distance, Levenshtein Metric) ]




아래처럼 사람 이름을 입력한 두 개의 문자열이 있다고 가정해보겠습니다. 

  • 문자열 1 (character string 1): Shawn Henry
  • 문자열 2 (character string 2): Shan Hennyy


'문자열 2'를 편집해서 '문자열 1'로 변환할 때 필요한 최소한의 조치를 생각해보면, 


(편집 조치 1) '문자열 2'의 4번째 위치에 'w'를 추가 (insert a 'w')

(편집 조치 2) '문자열 2'의 8번째 위치에 'n'을 'r'로 교체 (replace an 'n' with a 'r')

(편집 조치 3) '문자열 2'의 10번째 위치에 있는 'y'를 삭제 (delete the last 'y')


와 같이 3번의 편집 조치가 필요합니다.  따라서 '문자열 1'과 '문자열 2'의 편집 거리는 3입니다. 



'문자열 1'을 편집해서 '문자열 2'로 변환할 때 필요한 최소한의 조치를 생각해보면, 


(편집 조치 1) '문자열 1'의 4번째 위치의 'w'를 삭제 (delete a 'w')

(편집 조치 2) '문자열 1'의 9번째 위치의 'r'을 'n'으로 교체 (replace a 'r' with an 'n')

(편집 조치 3) '문자열 1'의 11번째 위치에 'y'를 추가 (insert an 'y')


이므로, 이렇게 계산해도 역시 '문자열 1'과 '문자열 2'의 편집 거리는 3입니다. 





이제 R 의 stringdist package를 사용해서 편집 거리 (edit distance, Levenshtein metric)를 계산해보겠습니다. 


(1) stringdist 패키지 설치 및 불러오기



# installing and loading of stringdist package

install.packages("stringdist")

library(stringdist)

 




(2) 문자열 편집 거리(edit distance, Levenshtein metric) 계산: stringdist()


문자열 "shawn henry"와 "shan hennyy", "show hurry" 문자열 간의 편집거리를 각각 계산해보겠습니다. 



> # to compute string edit distances

> # default method is 'osa', which is Optimal string alignment, (restricted Damerau-Levenshtein distance)

> stringdist(c("shawn henry"), c("shan hennyy", "show hurry"))

[1] 3 4

 



"shawn henry"와 "show hurry"의 편집거리를 계산하기 위해, 'show hurry'를 'shawn henry'로 변환하기 위한 최소 편집 조치를 살펴보면


(편집 조치 1) 'o'를 'a'로 변경 =>  shaw hurry

(편집 조치 2) 'n'을 추가   => shawn hurry

(편집 조치 3) 'u'를 'e'로 변경 => shawn herry

(편집 조치 4) 'r'을 'n'으로 변경 => shawn henry  (끝)


이므로, 편집거리는 4가 됩니다. 




(3) 여러개의 문자열 간의 편집 거리 (edit distance) 계산 결과를 행렬로 만들기: stringdistmatrix()


R stringdist 패키지에 들어있는 간단한 예제를 인용해서 예를 들어보겠습니다. 



> # to compute a dist object of class dist

> # => can be used by clustering algorithms such as stats::hclust

> stringdistmatrix(c("foo","bar","boo","baz"))

  1 2 3

2 3    

3 1 2  

4 3 1 2

> str_dist_mat <- as.matrix(stringdistmatrix(c("foo","bar","boo","baz")))

> str_dist_mat

  1 2 3 4

1 0 3 1 3

2 3 0 2 1

3 1 2 0 2

4 3 1 2 0

 




(4) 가장 유사한 문자열의 위치 찾기: amatch()


stringdist 패키지에는 'hello' 문자열과 편집 거리(edit distance, Levenshtein metric)가 가장 짧은 문자열의 위치를 찾아주는 amatch() 함수가 있습니다.  아래 예처럼 maxDist=2 라고 설정하면 편집 거리가 2를 넘어서는 문자열은 무시하게 됩니다.  동일 최소 편집거리 문자열이 여러개 있으면 앞에 위치한 문자열의 위치를 제시해줍니다. 



> # Approximate string matching

> # amatch returns the position of the closest match of x in table

> # by default, the OSA algorithm is used 

> # : Optimal string aligment (restricted Damerau-Levenshtein distance)

> amatch(c("hello"),c("hillu","hala","hallo", "hi"),maxDist=2)

[1] 3

 


- 'hello' 문자열과 'hillu' 문자열 간 편집 거리는 2 ('i'를 'e'로 교체, 'u'를 'o'로 교체), 

- 'hello' 문자열과 'hala' 문자열 간의 편집 거리는 3 ('a'를 'e'로 교체, 'l' 추가, 'a'를 'o'로 교체, 단, maxDist=2 이므로 고려 대상에서 제외됨), 

- 'hello' 문자열과 'hallo' 문자열 간의 편집 거리는 1 ('a'를 'e'로 교체)


이므로 편집 거리가 가장 짧은 'hallo' 의 3 을 반환합니다. 



* reference: stringdist package manual ( https://cran.r-project.org/web/packages/stringdist/stringdist.pdf)



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


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



Posted by R Friend R_Friend