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

이번 포스팅에서는 누적(cumulative)과 관련된 함수로서 {base} package에 기본으로 내장되어 있는 함수들인

 - 누적 합 (cumulative sums) : cumsum()

 - 누적 곱 (cumulative products) : cumprod()

 - 누적 최소값 (cumulative minima) : cummin()

 - 누적 최대값 (cumulative maxima) : cummax()

함수들에 대해서 아주 간단한 벡터(vector) 예제를 들어서 살펴보도록 하겠습니다.

이 함수들을 알고 있으면 달랑 한줄이면 끝인데요, 이 함수들을 모르고 있으면 누적 합, 누적 곱, 누적 최소값, 누적 최대값 구하는 프로그래밍하기가 은근히 어렵습니다. 

(1) 누적 합 (cumulative sums) : cumsum()

(1-1) 벡터 누적 합 (cumulative sum of vector)

아래와 같이 1~10까지의 정수 벡터가 있다고 합시다. 그리고 이에 대해 오른쪽으로 하나씩 이동해가면서 누적 합을 구해보면 아래의 표와 같습니다.

이걸 R base package의 cumsum() 함수를 이용해서 계산하면 아래와 같습니다.

> ##----------------------------------------------------------
> ## cumulative Sums, Products, and Minima or Maxima (Extremes)
> ##----------------------------------------------------------
> # {base} package
> 
> # (1) cumulative sums : cumsum()
> # (1-1) cumulative sums (vector) : cumsum()
> c(1:10)
 [1]  1  2  3  4  5  6  7  8  9 10
> 
> cumsum(1:10)
 [1]  1  3  6 10 15 21 28 36 45 55

이를 cumsum() 함수를 쓰지 않고 while 반복문을 사용해서 프로그래밍할 수도 있습니다.

(그런데, 계산 결과 제시되는 포맷이 위의 cumsum()과는 다르군요. ^^; )

> # cf. cumulative sums : while loop programing
> x <- 0
> i <- 1
> while( i <= 10) {
+         x = x + i
+         print(x) 
+         i = i + 1
+       }
[1] 1
[1] 3
[1] 6
[1] 10
[1] 15
[1] 21
[1] 28
[1] 36
[1] 45
[1] 55

누적 합(cumulative sums)은 그래도 종종 사용하는 만큼 벡터와 더불어서 데이터 프레임(dataframe) 예제도 곁들어서 살펴보겠습니다.

(1-2) 칼럼 1개만의 누적합을 구하는 방법과, (1-3) (numeric) 칼럼 모두의 누적합을 구하는 방법으로 나누어서 예를 들어보겠습니다.

(1-2) 데이터 프레임 : 변수 1개의 누적 합 (dataframe - cumulative sum of a column)

1~10까지의 정수인 x1 에 대한 누적합이므로 위의 (1-1) 벡터 예제와 결과가 동일하게 나옵니다.

결과 제시 포맷만 데이터 프레임의 새로운 변수(cum_sum_x1)로 생성된게 다를 뿐입니다.

> # (1-2) cumulative sums (dataframe - only 1 variable)
> #  : mydata$cum_sum <- cumsum(madata$var1)
> x1 <- c(1:10)
> x2 <- c(rep(1:5,2))
> x3 <- c(rep(2, 10))
> 
> x123 <- data.frame(x1, x2, x3) # making dataframe
> 
> x123
   x1 x2 x3
1   1  1  2
2   2  2  2
3   3  3  2
4   4  4  2
5   5  5  2
6   6  1  2
7   7  2  2
8   8  3  2
9   9  4  2
10 10  5  2
> 
> 
> x123$cum_sum_x1 <- cumsum(x123$x1) # cumulative sum of x1
> 
> x123
   x1 x2 x3 cum_sum_x1
1   1  1  2          1
2   2  2  2          3
3   3  3  2          6
4   4  4  2         10
5   5  5  2         15
6   6  1  2         21
7   7  2  2         28
8   8  3  2         36
9   9  4  2         45
10 10  5  2         55

(1-3) 데이터 프레임 : 숫자형 모든 변수의 누적 합 (dataframe - cumulative sum of all variables)

cumsum() 함수에 apply() 함수를 응용하면 소기의 목적을 달성할 수 있습니다. 

아래 예제의 1번째 행(row), 2번째 행, 3번째 행이 어떻게 구해졌는지만 풀어서 설명을 해보면요,

1st across_cum_sum = 1 + 1 + 2 + 1 = 5

2nd across_cum_sum = 5(1st across_cum_sum) + (2 + 2 + 2 + 3)(2nd row) = 14

3rd across_cum_sum = 14(2nd across_cum_sum) + (3 + 3 + 2 + 6)(3rd row) = 28

> # (1-3) cumulative sums (dataframe - all variables)
> #  : mydata$across_cum_sum <- cumsum(apply(mydata, 1, sum))
> x123$across_cum_sum <- cumsum(apply(x123, 1, sum))
> 
> x123
   x1 x2 x3 cum_sum_x1 across_cum_sum
1   1  1  2          1              5
2   2  2  2          3             14
3   3  3  2          6             28
4   4  4  2         10             48
5   5  5  2         15             75
6   6  1  2         21            105
7   7  2  2         28            144
8   8  3  2         36            193
9   9  4  2         45            253
10 10  5  2         55            325

(2) 누적 곱 (cumulative products) : cumprod()

1~10까지의 정수 벡터를 가지고 누적 곱 예를 들어보겠습니다.

숫자 하나씩 오른쪽으로 이동하면서 곱해나가면 됩니다. 

(아직 실전에서 써먹어 본적은 없음 -,-;)

> # (2) cumulative products (vector) : cumprod()
> c(1:10)
 [1]  1  2  3  4  5  6  7  8  9 10
> 
> cumprod(1:10)
 [1]       1       2       6      24     120     720    5040   40320  362880 3628800

(3) 누적 최소값 (cumulative minima) : cummin()

정수 (3, 2, 1, 2, 1, 0, 4, 3, 2) 의 순서로 이루어진 벡터에 대해서 누적 최소값을 구해보겠습니다. 

오른쪽으로 하나씩 숫자를 이동하면서 과거 지나온 숫자를 기억하고 있다가, 매번 평가를 통해 가장 작은 값(cumulative minima)을 반환하게 됩니다.

> # (3) cumulative minima (vector) : cummin()
> c(3:1, 2:0, 4:2)
[1] 3 2 1 2 1 0 4 3 2
> 
> cummin(c(3:1, 2:0, 4:2))
[1] 3 2 1 1 1 0 0 0 0

(4) 누적 최대값 (cumulative maxima) : cummax()

정수 (3, 2, 1, 2, 1, 0, 4, 3, 2) 의 순서로 이루어진 벡터에 대해서 누적 최대값을 구해보겠습니다.

오른쪽으로 하나씩 숫자를 이동하면서 과거 지나온 숫자를 기억하고 있다가, 매번 평가를 통해 가장 큰 값(cumulative maxima)을 반환하게 됩니다.

> # (4) cumulative maxima (vector) : cummax()
> c(3:1, 2:0, 4:2)
[1] 3 2 1 2 1 0 4 3 2
> 
> cummax(c(3:1, 2:0, 4:2))
[1] 3 3 3 3 3 3 4 4 4

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

참고로, R dplyr 패키지의 Window funciton 중에서
 - Cumulative aggregates : cumall() 함수, cumany() 함수, cummean() 함수
에 대해서는 여기 ( ☞ http://rfriend.tistory.com/246 )를 참고하세요.

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

지난번 포스팅에서는 분할적 군집화(Partitional Clustering) 중에서 프로토타입 기반(Prototype-based)의 군집화 기법인 K-중심 군집(K-Centroid Clustering)에 대해서 알아보았습니다.

이번 포스팅에서는 분할적 군집화(Partitional Clustering) 중에서 프로토타입 기반(Prototype-based) 군집화 기법의 두번째로 퍼지 군집 (Fuzzy Clustering)에 대해서 알아보겠습니다.

퍼지 군집 (Fuzzy Clustering)은 혼합 분포 군집(Mixture Distribution Clustering)과 함께 "Soft Clustering"이라고도 하는데요, 각 관측치가 (단 하나의 군집에만 속하는 것이 아니라) 여러 군집에 속할 수 있으며, 이를 각 군집에 속할 가능성(possibility), 확률(probability)로 제시해줍니다.

(가령, 관측치 3번은 군집 1 속할 가능성이 0.7, 군집 2에 속할 가능성이 0.3)

[ Partitional Clustering > Prototype-based > Fuzzy Clustering ]

퍼지 군집에 대해 자세히 살펴보기 전에 먼저 퍼지 논리(Fuzzy Logic)에 대해서 간략히 살펴보겠습니다.

Fuzzy를 영어사전에서 찾아보면 '애매 모호함' 이라는 뜻입니다. 

(뜻이) (be) vague (idea), obscure, fuzzy, inexplicit, puzzling, ambiguous (meaning), hazy (notion), dim (memory), doubtful, evasive, equivocal, elusive

Fuzzy Set Theory, Fuzzy Logic은 '애매 모호한 대상을 다루는 논리'입니다. 

고전 논리 연산(Boolean Logic)에서는 1(참, Truth) 아니면 0(거짓, False), 모 아니면 도, 아군 아니면 적군으로 단순 명쾌하게 구분을 합니다.

반면에 Fuzzy 논리(Fuzzy Logic)에서는 진리값(Truth Value)이 0~1 사이의 실수값(real nummber between 0 and 1)으로 확장하여 소속도(degrees of membership)에 대한 애매모호한 정도, 가능성의 정도(degrees of possibility)를 표현합니다.

[ Boolean Logic vs. Fuzzy Logic ]

위의 '키(Height)'를 예를 들어보면요, 고전 논리(Boolean Logic)에서는 키 175 cm 를 기준으로 해서, 175cm 미만은 키가 작은 집단 ('0', short group), 175cm 이상은 '키가 큰 집단' ('1', tall group)으로 구분하고 있습니다.  따라서 고전 논리에 따르면 키 174.999cm인 사람은 '키가 작은 집단'('0', short group)에 속하게 되며, 키 175.001 cm인 사람 (키 174.999 cm인 사람보다 단지 0.002 cm 큼)은 '키가 큰 집단' ('1', tall group)에 속하게 됩니다. 0.002 cm 차이가 두 사람이 속하는 그룹을 나누었는데요, 뭔가 불합리하다는 생각이 들지요?

반면에, 퍼지 논리 (Fuzzy Logic)에서는 소속도(degrees of membership)의 가능성(Possibility)을 진리 값(Truth Value)이 0~1 사이의 연속된 실수값으로 확장해서 나타내준다고 했지요?! 위의 키 예를 살펴보면요, 174.999 cm 인 사람이 '키가 큰 집단'에 속할 가능성이 0.499 이고 & '키가 작은 집단'에 속할 가능성(possibilty)이 0.501 입니다.  키가 175.001 cm 인 사람이 '키가 큰 집단'에 속할 가능성이 0.501 이고 & '키가 작은 집단'에 속할 가능성은 0.499 입니다. 

키나 몸무게, 나이, 시력, 청력... 등 연속된 개념에 대해서는 아무래도 고전 논리(Boolean Logic)보다는 퍼지 논리(Fuzzy Logic)이 애매 모호함의 정도(degrees)를 나타내기에 더 적합해 보입니다.  사실, 인간이 사용하는 용어, 표현, 개념의 "상당 부분이 애매모호한것 같습니다" (<= 이 표현 자체가 바로 애매 모호함 그 자체이지요. "상당 부분"? "같습니다"?).  퍼지 논리, 퍼지 집합론(Fuzzy set theory)을 이용하면 컴퓨터가 인간이 생각하고 표현하는 애매 모호함을 인식하고 연산할 수 있습니다.  

군집분석으로 다시 돌아와서 생각해보면요, K-중심 군집(K-Centroid Clustering)에서는 각 관측치가 특정 군집에 속하거나('1') 혹은 아니거나('0')의 둘 중 하나였습니다(each data point can only belong to exactly one cluster).  반면에, 퍼지 군집(Fuzzy Clustering)에서는 퍼지 이론(Fuzzy set theory)에 기반해서 각 관측치가 여러 군집에 동시에 속할 수 있으며(data points can potentially belong to multiple clusters), 각 군집별로 속할 가능성(degrees of possibility, probability)을 제시해줍니다.

[ Soft Clustering compared to Hard(Non-fuzzy) Clustering ]

관측치 중에서 각 군집과 군집의 중간 사이에 위치해서 특정 군집에 할당하기에 애매 모호한 관측치가 많이 있는 데이터셋이라면 '집단에 속할 가능성'을 나타내주는 퍼지 군집(Fuzzy Clustering)이 '단순 무식하게 여기 아니면 저기'로 배타적으로(exclusively) 나누어주는 K-중심군집(K-Centroid Clustering)보다 상대적으로 더 적합하다고 볼 수 있습니다.  (연산량이 많아지는 단점이 있긴 합니다만, 요즘엔 컴퓨팅 파워가 높으므로 예전 대비 문제가 덜 한편이죠)

개념에 대해 소개했으니, 이제 알고리즘으로 한단계 내려가 보겠습니다.

퍼지 군집 알고리즘으로 가장 많이 사용되는 것으로 Fuzzy C-means(FCM) Clustering Algorithm 입니다. FCM 알고리즘은 1973년 J.C.Dunn이 개발하였고, 1981년 J.C.Bezdek 이 발전시켰습니다. (* source : https://en.wikipedia.org/wiki/Fuzzy_clustering)  이번 포스팅에서는 Fuzzy C-means Clustering Algorithm을 가지고 설명하겠습니다.

Fuzzy C-means Clustering Algorithm은 K-means Algorithm 과 매우 유사합니다.  군집 내 관측치들 간의 유사성을 최대화하고 (즉, 군집의 중심과 관측치 간 거리의 합을 최소화), 군집 간 비유사성을 최대화 (즉, 군집의 중심 간 거리의 합을 최대화) 하는 최적 해(optimal solution)를 반복적인 연산을 통해 찾는 개념은 똑같습니다.  K 개의 군집을 분석가가 사전적으로 지정해주는 것도 같습니다.  유사성의 측도로서 거리(distance)를 사용하는 것도 같습니다 (정확히는 비유사성(dis-similarity)의 측도).

다만, 퍼지 군집에서는 각 관측치가 특정 군집에 속할 가능성, 가중치 w 를 계산하는 것이 다를 뿐입니다. (아래에는 군집 개수를 K로 통일해서 Fuzzy K-means Clustering 으로 사용하겠음)

[ Fuzzy K-means Clustering Algorithm ]

위의 알고리즘 내용을 수식으로 표현하기 전에 표기법부터 정리해보죠. 군집화를 하려고 하는 데이터셋 집합이 m개의 변수, n개의 관측치, K개의 군집이 있고, 각 군집에 속할 가능성을 가중치 값 w 로 표기한다고 하면 아래와 같습니다.  

[ 표기법 (Notation) ]

퍼지 군집 모형은 아래 두 개의 분할 조건을 만족합니다.

1) 데이터 가 각 군집 에 속할 가능성의 가중값 합은 1이다. (즉, 확률 합이 1)

    ---- (식 1)

2) 각 군집 는 하나 이상의 데이터가 0이 아닌 가중값을 가지며, 그 군집의 모든 가중값이 1이 될 수는 없다. (모두 0이거나 모두 1이면 군집분석 효용성 없음)

   ---- (식 2)

퍼지 군집을 하는 원리가 '군집 내 유사성 최대화 (즉, 관측치와 퍼지 군집 중심과의 거리 합 최소화)', '군집 간 비유사성 최대화 (즉, 군집의 중심 간 거리 합 최대화)' 하는 것이라고 했습니다 (K-means clustering 과 원리 동일). 

'군집 내 유사성 최대화 (즉, 관측치와 퍼지 군집 중심과의 거리[d(xi, ck)] 합 최소화)'를 달리 표현하면 '퍼지 군집 내 오차 제곱 합(SSE : Sum of Squared Error) 최소화'라고 할 수 있습니다.

  ---- (식 3)

위의 식에서 p는 1보다 큰 상수로서, 가중값의 정도를 조절하는 매개변수(parameter) 입니다. p 값이 1이면 K-means Clustering과 결과가 비슷해집니다. (K-means Clustering을 Fuzzy Clustering의 특수한 경우이며,  Fuzzy Clustering이 더 포괄적이라고 생각할 수 있습니다.)

위의 (식 3)에서 SSE를 최소로 하는 각 군집의 평균 를 구하기 위해서 (식 3)을 에 대해 편미분(partial derivative with respect to ) 한 후 '0'으로 놓고 연립방정식을 풀면 아래와 같은 해를 얻습니다.  (K-means clustering 에서 각 군집의 중심 평균을 구할 때는 해당 군집의 관측치 값만 사용하는데 반해서, Fuzzy K-means Clustering 에서는 전체 관측치에다가 각 군집에 속할 가능성인 가중치를 곱한 값을 사용함)

   ---- (식 4)

위의 (식 3)에서 SSE를 최소로 하는 가중값 를 구하기 위해서 (식 3)을 에 대해 편미분(partial derivative with respect to ) 한 후 '0'으로 놓고 연립방정식을 풀면 아래와 같은 해를 얻게 됩니다.

     ---- (식 5)

p 값이 커지면 커질수록 각 군집의 평균이 전체 평균에 가까워져서 한 군집으로 데이터를 분류하는 것이 점점 더 모호(fuzzier)해집니다. 일반적으로 위의 (식 5)의 가중값 의 재계산식을 간편히 하기 위해 p=2를 많이 사용합니다.

p=2 를 (식 5)에 대입하면 (식 5)가 아래와 같이 계산하기 용이하게 정리됩니다.

       ---- (식 6)

위의 (식 6)의 분자를 살펴보면, 관측치 가 군집 에 속할 가능성, 가중치 는 관측치 와 군집의 중심 의 거리(distance) 제곱에 반비례함을 알 수 있습니다. 즉 관측치와 군집의 중심 간 거리가 짧을 수록 (즉, 유사할 수록, 그 군집에 속할 가능성이 높을 수록) 는 커지게 됩니다. 

(* 출처 : 이정진)

이때 분모는 각 데이터와 군집 1 ~ K 까지의 거리 제곱의 역수의 합이며, 이것으로 분자(특정 군집 K와 관측치 간의 거리 제곱의 역수)를 나누어주게 되면 여러 군집에 속하는 가중값의 합이 1이 되도록 해주는 표준화 상수 역할을 하게 됩니다.

간단한 예제를 가지고 Fuzzy K-means Clustering을 반복 수행해보겠습니다. 

n = 4 (관측치 4개),

m = 2 (변수 2개),

K = 2 (군집 2개),

p = 2 (가중치 계산 상수 parameter),

군집 중심과 관측치간의 거리는 유클리드 제곱 거리(Squared euclidean distance) 인

간단한 예제이며, 엑셀로 수식 걸어서 반복 수행하였습니다. (아래 첨부한 엑셀 파일 참조하세용~)

Fuzzy_Clustering_example.xlsx

[ 데이터셋 ]

[ x1, x2 축 기준 관측치 산점도 ]

무작위로 obs 2, obs 4에 군집 1 가중값(wi1)으로 0.8, 군집 2 가중값(wi2)으로 0.2를 할당하고, obs 1, obs 3에는 군집 1 가중값(wi1)으로 0.2, 군집 2 가중값(wi2)으로 0.8을 할당하였습니다. (각 관측치의 모든 군집이 가중값 합의 1이 되어야 하므로 군집 1 가중값이 정해지면 나머지 군집 2의 가중값은 1-wi1 으로 자동으로 정해짐)

관측치가 군집 에 속할 가능성인 가중값  는 (식 6)에 의해서, 군집 의 중심  는 (식 4)에 의해서 구했습니다.

가중값 가 수렴할 때까지 위 계산을 반복을 합니다.

(위의 표 'Fuzzy K-means Clustering Algorithm' 참조)

5번째 반복과 6번째 반복의 가중값의 변화가 거의 없으므로 (즉, 수렴하였으므로) 퍼지 군집 알고리즘을 종료합니다.

결과적으로 obs 1 과 obs 2 는 군집 2 (cluster 2)에 속하고, obs 3 과 obs4 는 군집 1로 분류가 되었습니다.

(obs 1의 w11 = 0.02, w12 = 0.98,

 obs 2의 w21 = 0.02, w22 = 0.98,

 obs 3의 w31 = 0.96, w32 = 0.04,

 obs 4이 w41 = 0.94, w42 = 0.06  이므로)

위에 제시한 예제의 산점도를 보면 obs 1과 obs 2가 서로 인접해 있으며, obs 3과 obs 4가 서로 인접해 있으므로 군집화가 제대로 된 셈이네요. 

비록 처음에 무작위로 가중값을 부여했을 때 obs 1과 obs 3을 군집2로 가중치를 0.8 할당, obs 2와 obs 4를 군집1로 가중치를 0.8 할당하였습니다만, 6차례의 반복을 거치면서 각 관측치별 가중치도 새로 계산하고, 군집 1과 군집 2의 중심(K-평균)도 새로 계산하면서 군집화(clustering)을 반복하다보니 인접한(유사한) 관측치끼리 군집으로 잘 묶였습니다.

처음에 가중값 부여할 때 무작위로 한다고 했는데요, 여기서 무작위로 부여하는 숫자가 바뀌면 (무작위 이므로 뭐가 될지 모름 -_-;) 물론 군집화의 결과가 바뀔 수 있다는 점은 알고 계시구요. (K-means clustering 도 초기 중심값이 무작위로 할당되다보니 분석 돌릴 때마다 군집화 결과가 바뀔 수 있다는 건 동일함)

아래에 R 예제에서는 초기 rational starting point를 행렬로 입력할 수 있는 기능이 있으므로, 관측값별 초기 가중값을 합리적으로 부여할 수 있는 상황이라면 이용할 만 하겠습니다.

R의 'fclust' Package를 가지고 위에서 소개했던 예제에 대해서 Fuzzy clustering을 수행해보겠습니다.

1) 먼저, fclust Package 를 설치하고 로딩해보겠습니다.

> ##---------------------------------------------
> ## Fuzzy K-means Clustering : R fclust Package
> ##---------------------------------------------
> install.packages("fclust")
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.3/fclust_1.1.2.zip'
Content type 'application/zip' length 197520 bytes (192 KB)
downloaded 192 KB

package ‘fclust’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\Administrator\AppData\Local\Temp\RtmpkzHe81\downloaded_packages
> library(fclust)

2) Dataset 준비

  위의 이론 설명에서 사용했던 2개 변수(x1, x2), 4개 관측치 예제 데이터를 똑같이 사용하겠습니다.

> # dataset
> x1 <- c(-1, -2, 1, 3)
> x2 <- c(-2, -1, 3, 2)
> x1_x2_d.f <- data.frame(x1, x2)
> x1_x2_d.f
  x1 x2
1 -1 -2
2 -2 -1
3  1  3
4  3  2

3) rational starting point Matrix U

위의 이론 설명에서 사용했던 값을 그대로 사용해서요, 관측값 1번이 군집1에 속할 가중값은 0.2, 군집2에 속할 가중값은 0.8, ... , 관측값 4번이 군집1에 속할 가중값이 0.8, 군집2에 속할 가중값은 0.2 의 초기값을 행렬(Matrix)로 만들었습니다.  이걸 아래의 FKM() 함수의 startU 옵션의 할당값으로 지정해줄겁니다.

> # rational starting point Maxtix U
> rational_starting_point <- matrix(c(0.2, 0.8, 0.2, 0.8, 0.8, 0.2, 0.8, 0.2), 
+                                   nrow = 4, ncol = 2, byrow = F)
> rational_starting_point
     [,1] [,2]
[1,]  0.2  0.8
[2,]  0.8  0.2
[3,]  0.2  0.8
[4,]  0.8  0.2

4) Fuzzy K-means Clustering using FKM() function of fclust package

fclust Package이 FKM() 함수의 매개변수 및 옵션에 대해서 소개하자면 아래와 같으며, 데이터셋(행렬 또는 데이터 프레임) 할당해주는 X, 군집의 개수를 지정해주는 k, 퍼지 매개변수(parameter) m (위의 '식 5'번의 p) 에 대해서만 옵션을 설정해보겠으며, 나머지는 default 설정 그대로 사용하겠습니다.

(같은 척도의 데이터이므로 표준화 필요 없음. default 가 no standardization)

( * source : Paolo Giordani, Maria Brigida Ferraro)

퍼지군집 분석 결과는 아래와 같습니다.

> # Fuzzy K-means clustering with FKM() fuctnion of fclust package
> x1_x2_FKM <- FKM(X = x1_x2_d.f, # Matrix or data.frame
+                  k = 2, # Number of clusters (default: 2)
+                  m = 2, # Parameter of fuzziness (default: 2)
+                  startU = rational_starting_point) 
> # startU : Rational starting point for the membership degree matrix U
> 
> 



> # Fuzzy K-means clustering results
> x1_x2_FKM

 Fuzzy clustering object of class 'fclust' 
 
 Number of objects: 
 4
 
 Number of clusters: 
 2
 
 Closest hard clustering partition: 
1 2 3 4 
2 2 1 1 

 Membership degree matrix (rounded): 
  Clus 1 Clus 2
1   0.02   0.98
2   0.02   0.98
3   0.95   0.05
4   0.96   0.04

 Available components: 
 [1] "U"     "H"     "clus"  "value" "cput"  "iter"  "k"     "m"     "stand" "Xca"  
[11] "X"     "call" 

R의 fclust Package 의 FKM() 함수로 퍼지 군집화를 한 결과와, 위에서 엑셀을 가지고 6번 반복해서 푼 결과가 일치함을 알 수 있습니다. (obs3, obs4의 degree of membership, ie, weight 가 소숫점 두째자리에서 약간 다르기는 하지만 무시할만 합니다. 엑셀로는 반복을 6번하고 멈추었는데요, R은 14번 반복을 했네요.  이는 R fclust Package의 수렴 기준(Convergence criterion)의 default 값이 '1e-9' 으로서 매우 매우 작기 때문에 반복을 좀더 많이 했습니다.

> x1_x2_FKM$iter  # number of iteration
Start 1 
     14

fclust Package FKM() 함수의 분석결과 객체에 대해 소개하자면 아래와 같으며, 필요한 정보는 indexing 해서 사용하면 유용합니다.

( * source : Paolo Giordani, Maria Brigida Ferraro)

예를 들어, 소속 정도의 가중값을 알고 싶다면 'U' 객체를 indexing 해오면 됩니다.

> # Membership degree matrix
> x1_x2_FKM$U
      Clus 1     Clus 2
1 0.01684418 0.98315582
2 0.01734390 0.98265610
3 0.95397165 0.04602835
4 0.96353078 0.03646922

각 퍼지 군집의 중심 위치도 궁금하지요?  이때 쓰는게 'H' 객체입니다.  Cluster 1 은 x1 중심좌표가 '2.008850', x2 중심좌표가 '2.493750' 으로서 우상단에 위치하고 있으며, Cluster 2는 x1 중심좌표가 '-1.493918', x2 중심좌표는 '-1.492924' 로서 좌하단에 위치하고 있군요.

> # Prototype matrix : H
> x1_x2_FKM$H
              x1        x2
Clus 1  2.008850  2.493750
Clus 2 -1.493918 -1.492924

변수가 2개인 2차원 데이터이므로 산점도 그래프로 그려보면 좀더 명확하게 이해할 수 있겠네요.  아래의 그래프를 보면 2개의 군집별로 색깔이 다르게 나와있구요, '*' 표시는 각 군집의 중심을 의미합니다. 바로 위해서 x1_x2_FKM$H 로 indexing했던 바로 그 좌표값입니다.

> # Fuzzy K-means Clustering plot
> plot(x1_x2_FKM)

이상으로 퍼지 군집(Fuzzy K-means Clustering Algorithm)에 대한 소개를 마치겠습니다.

다음번 포스팅에서는 '혼합분포군집(Mixture Distribution Clustering)' 모형에 대해서 알아보겠습니다.

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

[Reference]

- Michael Negnevitsky, "인공지능 개론(Artificial Intelligence)' (2nd Edition), 한빛아카데미, 2013

- 이정진, "R, SAS, MS-SQL을 활용한 데이터마이닝", 자유아카데미, 2011

- Fuzzy C-means Clustering Algorithm : https://en.wikipedia.org/wiki/Fuzzy_clustering

- R fuzzy clusterin package 'fclust' : Paolo Giordani, Maria Brigida Ferraro,  https://cran.r-project.org/web/packages/fclust/fclust.pdf

지난번 포스팅에서는 (1) 응집형 계층적 군집화(Agglomerative Hierarchical Clustering) 방법 5가지(단일연결법, 완전연결법, 평균연결법, 중심연결법, Ward연결법) 중에서, 오차제곱합의 증분으로 군집 간 (비)유사성을 측정해서 군집화를 하는 Ward 연결법에 대해서 알아보았습니다.

이번 포스팅부터는 (2) 분할적 군집화(Partitional Clustering) 기법 중에서 프로토타입 기반(Prototype-based) 의 K-중심 군집화(K-centroid clustering)에 대해서 알아보도록 하겠습니다.

[ 분할적 군집 (Partitional Clustering, Non-hierarchical Clustering) ]

복습하는 차원에서 한번 더 복기를 하자면,

- 응집형 계층적 군집화(Agglomerative Hierarchical Clustering)은 각각의 객체에서 시작해서 유사성 척도(proximity measure)에 의거해 유사한 (거리가 짧은) 객체들을 Tree 형태의 계층적 군집으로 차근 차근 묶어가는 기법입니다. (a set of nested clusters organized as a hierarchical trees).  일단 한번 군집으로 묶이면 그 객체는 계속 그 군집에 속하게 되고 변동이 없으며, 계층 구조의 Tree에서 어느 수준을 기준으로 하느냐에 따라서 군집이 달라지게 됩니다.

- 분할적 군집화(Partitional Clustering)은 객체가 하나의 군집에 exclusive하게 속하도록 군집을 형성합니다. (A division data objects into non-overlapping subsets such that each data object is in exactly one subset). 분할 방법에는 프로토타입 기반(Prototype-based), 분포 기반(distribution-based), 밀도 기반(Density-based), 그래프 기반(Graph-based) 기법이 있습니다.

참고로,

- Hard clustering은 객체별로 어느 군집에 속할지를 명시적으로 할당하는 기법이며,

K-중심군집은 Hard clustering에 속합니다.

- Soft clustering은 각 객체가 어느 군집에 속할지를 가중치(weight)나 확률(probability)로서 가능성 정도를 나타내주는 기법으로서, Fuzzy Clustering과 혼합분포군집(Mixture Distribution Clustering)이 이에 속합니다.

이번 포스팅에서 소개할 분할적 군집화는 이중에서 프로토타입 기반(Prototype-based) 기법 중에서도 K-중심군집(K-centroid Clustering) 모형이 되겠습니다.

프로토타입 기반 군집화(Prototype-based Clustering)는 미리 정해놓은 각 군집의 프로토타입에 각 객체가 얼마나 유사한가 (혹은 가까운가)를 가지고 군집을 형성하는 기법입니다. 

K-중심군집에서는 연속형 데이터의 경우 평균(Mean)이나 중앙값(Median)을 그 군집의 프로토타입으로 하며, 이산형 데이터인 경우는 최빈값(Mode)이나 메도이드(Medoid)라고 해서 해당 군집을 가장 잘 표현할 수 있는 측도를 정해서 프로토타입으로 정하게 됩니다.

보통 군집분석을 공부한다고 했을 때 가장 많이 회자되고, 가장 처음에 배우는 기법이 아마도 'K-평균 군집화(K-means Clustering)이 아닐까 싶습니다.  그런데 앞서 소개드린 것처럼 군집분석 기법에는 정말 많은 알고리즘이 있습니다. K-평균 군집은 그 중에서도 일부에 해당하는 기법일 뿐이며, 프로토타잎도 데이터 형태에 따라서 '평균(Mean)'을 쓰는 K-means Clustering, '중앙값(Median)'을 쓰는 K-median Clustering, '메도이드(Medoid)'를 쓰는 K-medoid Clustering 등으로 세분화된다는 점은 알아두시면 좋겠습니다.  이들을 모두 묶어서 'K-중심군집(K-centroid Clustering)'이라고 합니다.

여기서 'K'는 '군집의 수(number of clusters)'를 나타내는 모수(parameter)로서, 분석가가 사전에 정해주어야 합니다.  참고로, 군집의 수 K를 미리 지정해주어야 하는 군집분석 기법으로는 이번 포스팅의 주제인 K-중심군집(K-centroid Clustering), 그리고 퍼지군집(Fuzzy Clustering), 혼합분포 군집(Mixture Distribution Clustering) 등이 있습니다.

군집의 수 K를 정하는 문제가 참 중요한데요, 좀 어렵고 주관적이고 애매한 부분이 있습니다.

군집 개수별 전체 분산에서 설명가능한 분산의 비율의 꺽인 부분을 보는 The elbow method 나, 동일 군집 내 객체 간의 거리와 다른 군집 간의 거리를 고려한 실루엣 방법(Silhouette method)을 통해서 적정 군집의 개수를 선택하는 방법은 https://rfriend.tistory.com/585 또는 Wikipedia 를 참고하세요.

저는 보통은 위의 방법과 더불어 업에 대한 이해를 바탕으로 분석 목적을 감안하여 복수의 k를 지정해서 군집분석을 수행한 후에, 군집에 대한 profiling을 해보고, 가장 적합하다고 판단되는 k를 정하곤 했습니다.  다분히 분석가의 업에 대한 경험/이해도와 주관이 많이 들어가고, Biz 활용 목적과 현실적 제약조건도 고려해야 하며, 또 시행착오와 profiling을 통한 오랜 탐색이 필요합니다. 

[ 군집의 수 K 정하기의 애미모호함 ?  어려움?? ]

이번 포스팅에서는 일반적으로 가장 많이 사용하는 'K-평균 군집(K-means Clustering)'에 대해서만 대표로 설명을 하도록 하겠습니다.

K-평균 군집(K-means Clustering)의 원리를 알아보면요,

1) 군집 내 응집도 최대화(maximizing cohesion within cluster) : 군집 內 중심(centroid)과 해당 군집의 각 객체 간 거리의 합 최소화

2) 군집 간 분리도 최대화(maxizing separation between clusters) : 각 군집의 중심(centroid) 間 거리 합 최대화 

의 두 가지 목적함수를 만족하는 해(解, solution)를 찾는 것입니다.

[ K-평균 군집의 개념 및 원리 (Concept and Principle of K-means Clustering) ]

즉, 군집분석은 결국 위의 두 목적함수에 대한 최적화 (optimization of global objective function) 문제임을 알 수 있습니다.  복잡도(complexity)를 살펴보면, 군집의 수가 k, 차원의 수가 d, 객체의 수가 n 일 때 입니다. (* 출처 : https://en.wikipedia.org/wiki/K-means_clustering)  기본적으로 객체의 수(n)가 많을 수록 시간이 오래걸리며, 특히 변수의 수(d)와 군집의 수(k)가 늘어날 수록 지수적으로 계산 시간이 증가함을 알 수 있습니다.  따라서 허접한 변수들 몽땅 때려넣고 군집화하라고 컴퓨터한테 일 시킬 것이 아니라, 똘똘하고 핵심적인 변수를 선별해서 차원을 줄인 후에 군집분석을 실행하는 것이 연산시간을 줄이는 측면에서나, 또 Biz. 목적에 맞게 잘 군집화가 되도록 하는 측면에서나 중요합니다.

다음으로, 'K-중심군집 알고리즘(K-centroid Clustering Algorithm)'에 대해서 알아보겠습니다.  알고리즘은 알고나면 허무할 정도로 정말 쉽습니다. ^^' 컴퓨터야 중심과 객체 간 거리를 반복적으로 계산하느라 죽어나겠지만 말이지요.

[ 'K-중심군집 알고리즘(K-centroid Clustering Algorithm) ]

위 알고리즘을 도식화해서 예를 들어보면 아래와 같습니다. 매번의 반복(iteration) 마다 군집의 중심이 새로 계산되고, 새로 바뀐 중심과 각 객체간 거리가 다시 계산되면서 군집이 계속 동적으로 바뀌다가 더 이상 변동이 없이 수렴될 때까지 반복을 하게 됩니다. 

[ K-평균군집화의 개념적인 예시 (Conceptual exmaple of K-means Clustering) ]

위의 개념적 예시 그림에서도 짐작할 수 있듯이, 초기 중심값(initial centroid)가 바뀌면 군집 묶이는게 바뀔 수 있습니다.  이런 점때문에 B2C에서 고객세분화(customer segmentation)에 K-means clustering을 적용하게 되면 매번 군집화를 시도할 때마다 군집에 묶이는 고객이 달라질 수 있어서 문제가 될 소지가 있습니다.  (참고로, R에서 set.seed = 1234 처럼 set.seed() 함수를 사용하면 무작위수가 매번 실행때마다 동일하게 할 수 있음)

이번 포스팅은 K-평균 군집(K-means Clustering)을 다룬다고 했으므로, 중심(centroid)는 '평균(Mean)'이 되겠지요.

그리고 유사도(proximity, similarity)는 보통은 유클리드 거리(Euclidean distance,

) 혹은 유클리드 제곱 거리(SSE, SUM of Squared Error, R에서 사용하는 거리)를 사용합니다.  데이터 특성에 따라서 유사도 측도로서 적합한 Measure를 선택해서 분석하시면 되겠습니다. (☞ 유사도 측정 Distance Measures 참고 : http://rfriend.tistory.com/199)

거리를 가지고 유사성을 측정한다고 했는데요, 이러다 보니 K-평균군집(K-means Clustering)은 노이즈나 이상치에 민감(Sensitive to Noise and Outlier)한 단점이 있습니다.  평균보다는 중앙값이 이상치에 덜 민감하므로 이상치로 인한 왜곡이 우려되면 K-중앙값군집(K-median Clustering)이 대안이 될 수 있겠네요.  아니면 탐색적 분석 단계에서 이상치를 제거하는 것도 방법이 될 수 있겠고요.

마지막으로, 만약 여러 변수들의 계측 단위(scale)이 서로 다르다면 사전에 표준화(standardization)를 해줄 필요가 있습니다.  안그러면 측정 단위가 큰 특정 한, 두개의 변수에 의해 군집화가 휘둘릴 수 있기 때문입니다.  보통 표준정규분포로 z 표준화를 많이 사용합니다. (R의 scale() 함수)

K-중심 군집(K-Centroid Clustering) 이론에 대해서는 왠만큼 소개를 한듯 하니, 이제 데이터셋을 가지고 R script 를 써가면서 실습을 해보겠습니다.  R script 가 무척 짧아서 당황하실수도 있다는 점 미리 안내드립니다. ㅋㅋ

실습 데이터셋은 iris 입니다. 

Sepal.Length, Sepal.Width, Petal.Length, Petal.Width 의 4개의 변수를 가지고 있고, 150개의 관측치를 가지고 있는 데이터프레임입니다.

결측값은 없이 깨끗하네요.

측정 척도의 단위(scale)가 4개 변수 모두 길이를 재는 동일 단위이므로 iris 데이터셋의 경우 별도로 표준화르 할 필요는 없겠네요. 

> ##-----------------------------------
> ## K-centroid clustering
> ##-----------------------------------
> 
> # dataset : iris
> str(iris)
'data.frame':	150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
> 
> # checking missing value
> colSums(is.na(iris))
Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
           0            0            0            0            0

산포도(Scatter Plot)를 그려보면 아래와 같습니다.

> # scatter plot of iris
> panel.fun <- function(x, y, ...) {
+   horizontal <- (par("usr")[1] + par("usr")[2]) / 2; 
+   vertical <- (par("usr")[3] + par("usr")[4]) / 2; 
+   text(horizontal, vertical, format(abs(cor(x,y)), digits=2)) 
+ }
> 
> pairs(iris[1:4], 
+       pch = 21, bg = c("red","green3","blue")[unclass(iris$Species)], 
+       upper.panel=panel.fun, 
+       main = "Statter Plot of Iris Dataset")

K-means Clustering이 중심과의 거리를 가지고 군집을 묶는 방법이다보니, 위의 산포도를 보면 Petal.Length와 Petal.Width 의 두개의 변수를 가지고 군집화(Clustering)를 하는 것이 제일 좋을 것 같군요.(아마 Sepal.Width와 Petal.Length 의 두개 변수를 사용해서 K-means Clustering을 돌리면 좌측 상단의 두 Species가 반반씩 잘못 섞여서 군집화가 될겁니다. 왜 그럴지는 한번 생각해보시길...) 

Petal.Length와 Petal.Width 로 산점도를 그려보면 아래와 같습니다.

[plot 1 : original scatter plot]

> # Original Scatter Plot of Iris Petal.Length & Petal.Width 
> install.packages("ggplot2")
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.3/ggplot2_2.1.0.zip'
Content type 'application/zip' length 2001758 bytes (1.9 MB)
downloaded 1.9 MB

package ‘ggplot2’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\Administrator\AppData\Local\Temp\Rtmpi20q05\downloaded_packages
> library(ggplot2)
> iris_plot <- ggplot(data=iris, aes(x=Petal.Length, y=Petal.Width, colour=Species)) + 
+   geom_point(shape=19, size=4) + 
+   ggtitle("Original Scatter Plot of Iris Petal.Length & Petal.Width")
>   
> iris_plot

> # text annotation with Species Name
> iris_plot_2 <- iris_plot +
+   annotate("text", x=1.5, y=0.7, label="Setosa", size=5) + # text annotation
+   annotate("text", x=3.5, y=1.5, label="Versicolor", size=5) + 
+   annotate("text", x=6, y=2.7, label="Virginica", size=5)
> 
> iris_plot_2
> # adding shadowed box
> iris_plot_3 <- iris_plot_2 +
+   annotate("rect", xmin=0, xmax=2.6, ymin=0, ymax=0.8, alpha=0.1, fill="red") + 
+   annotate("rect", xmin=2.6, xmax=4.9, ymin=0.8, ymax=1.5, alpha=0.1, fill="green") + 
+   annotate("rect", xmin=4.8, xmax=7.2, ymin=1.5, ymax=2.7, alpha=0.1, fill="blue")
> 
> iris_plot_3

Noise나 이상값(outlier)는 없어보이므로 별도 전처리 없이 그대로 데이터 가져다 쓰면 되겠네요.

자, 이제 드디어 K-평균군집(K-means Clustering)을 R로 돌려보겠습니다.  군집의 개수 K = 3 으로 하겠습니다. (iris 품종으로 Setosa, Versicolor, Virginica 의 3종류가 있다고 우리가 이미 알고 있기 때문에 K 3으로 한 것임.  일반적으로 K-평균군집분석을 할 때 Y Label에 대해서 모르거나 없는 상태에서 비지도학습(Un-supervised Learning)으로서 데이터 속에 숨겨진 패턴을 컴퓨터 보고 한번 찾아보라고 마이닝을 시키게 됨.  다시 한번 말하지만, 이번 iris 데이터셋은 Y Lable을 알고 있는 상태에서 K-평균군집분석 결과를 좀더 이해하기 쉽도록 비교해서 보여주는 것 뿐이며, Y Lable 모른 채, 혹은 없는 상태에서 군집분석을 수행하게 됨)

K-평균군집의 R script는..... 음..... 아래 한줄이 끝입니다. -_- 

kmeans(dataset, k)  # k : number of clusters

> # K-means clustering with Petal.Length and Petal.Width
> iris_k_means <- kmeans(iris[,c("Petal.Length", "Petal.Width")], 3) 

K-평균군집의 객체 iris_k_means 를 호출하면 아래와 같이 K-평균 군집 결과를 일목요연하게 볼 수 있습니다. 

- K-means clustering with 3 clusters of sizes 50, 52, 48

   : 군집의 개수(k)가 3개, 군집 1/2/3별 크기가(개체 개수) 50개, 52개, 48개

- Cluster means

   : 군집 1/2/3 별 두개의 변수 Petal.Length, Petal.Width의 평균 좌표 

     (=> profiling 하기에 good!)

- Clustering vector

   : 각 개체별 군집 벡터
- Within cluster sum of squares by cluster

   : 각 군집의 중심(centroid)와 각 군집에 속한 개체간 거리의 제곱합

- Available components

   : 군집분석 결과의 구성요소
     => 필요한거 있으면 이 객체(object)를 indexing해서 쓰면 요긴함
 

> iris_k_means
K-means clustering with 3 clusters of sizes 50, 52, 48

Cluster means:
  Petal.Length Petal.Width
1     1.462000    0.246000
2     4.269231    1.342308
3     5.595833    2.037500

Clustering vector:
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
 [52] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3
[103] 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3

Within cluster sum of squares by cluster:
[1]  2.02200 13.05769 16.29167
 (between_SS / total_SS =  94.3 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"        
[8] "iter"         "ifault"   

아래처럼 names() 함수를 써서 군집분석 결과의 객체들에 무엇이 있는지 확인해볼 수도 있습니다.

- cluster : 각 개체별 할당된 군집 번호, 1부터 k번까지 군집 숫자 

             (A vector of integers (from 1:k) indicating the cluster to which each point is allocated.)

- centers : 군집의 중심 (centroid) 좌표 행렬 (A matrix of cluster centres.)

- totss : 제곱합의 총합 (total sum of squares)

- withinss : 군집 내 군집과 개체간 거리의 제곱합 벡터.

               (Vector of within-cluster sum of squares, one component per cluster)

- tot.withinss : 군집 내 군집과 개체간 거리의 제곱합의 총합, 즉, sum(withinss)

                   (Total within-cluster sum of squares, i.e. sum(withinss))

- betweenss : 군집과 군집 간 중심의 거리 제곱합

                   (The between-cluster sum of squares, i.e. totss-tot.withinss)

- size : 각 군집의 개체의 개수 (The number of points in each cluster)

- iter : 반복 회수 (The number of (outer) iterations)

- ifault : 전문가용의 발생 가능한 알고리즘 문제의 indocator (?? 이거 뭔지 잘 모르겠음... -_-??)
           (integer: indicator of a possible algorithm problem – for experts)

> # objects of k-means clustering in R
> names(iris_k_means)
[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"        
[8] "iter"         "ifault"

그럼, 군집의 크기(객체의 개수)를 한번 확인해볼까요?  table() 함수를 써도 되고 size 객체를 가져다가 확인해도 됩니다. 두 가지 방법 모두 아래에 소개하였습니다.

> # cluster size
> table(iris_k_means$cluster)

 1  2  3 
50 52 48 

> iris_k_means$size
[1] 50 52 48
> 
> prop.table(iris_k_means$size)
[1] 0.3333333 0.3466667 0.3200000

각 군집 1, 2, 3의 중심(centroid)도 확인해보겠습니다.

> # centroid of clusters
> iris_k_means$centers
  Petal.Length Petal.Width
1     1.462000    0.246000
2     4.269231    1.342308
3     5.595833    2.037500 

마지막으로, 각 개체들을 군집별로 색깔을 달리해서 산점도(scatter plot) 그려보겠습니다.  덤으로 각 군집의 중심(centroid)을 검정색 점으로 해서 덮입혀보았습니다.

위의 [plot 1 : original scatter plot]에서 iris species 로 색깔을 구분해서 그린 산점도와 아래의 군집(cluster)별로 색깔을 구분해서 그린 산점도가 거의 유사하지요?

> # Scatter plot of Iris' K-menas clusters
> cluster <- iris_k_means$cluster
> iris_k_means_x_y <- cbind(cluster, iris[,c("Petal.Length", "Petal.Width")])
> head(iris_k_means_x_y)
  cluster Petal.Length Petal.Width
1       1          1.4         0.2
2       1          1.4         0.2
3       1          1.3         0.2
4       1          1.5         0.2
5       1          1.4         0.2
6       1          1.7         0.4
> 
> sapply(iris_k_means_x_y, class)
     cluster Petal.Length  Petal.Width 
   "integer"    "numeric"    "numeric" 
> iris_k_means_x_y <- transform(iris_k_means_x_y, 
+                               cluster = as.factor(cluster))
> 
> sapply(iris_k_means_x_y, class)
     cluster Petal.Length  Petal.Width 
    "factor"    "numeric"    "numeric" 
> 
> 
> library(ggplot2)
> iris_k_means_x_y_plot <- ggplot(data=iris_k_means_x_y, 
+                                 aes(x=Petal.Length, y=Petal.Width, colour=cluster)) + 
+   geom_point(shape=19, size=4) + 
+   ggtitle("Scatter Plot of Iris' K-means clusters")
> 
> iris_k_means_x_y_plot
> 
> 
> # adding centroid points of cluster 1, 2, 3
>   # centers(centroids) by cluster 1, 2, 3
> iris_k_means_centers <- iris_k_means$centers
> iris_k_means_centers
  Petal.Length Petal.Width
1     1.462000    0.246000
2     4.269231    1.342308
3     5.595833    2.037500
> 
> iris_k_means_x_y_plot_2 <- iris_k_means_x_y_plot +
+   annotate("point", x = 1.462, y = 0.246, size = 6, color = "black") + 
+   annotate("point", x = 5.595, y = 2.037, size = 6, color = "black") + 
+   annotate("point", x = 4.269, y = 1.342, size = 6, color = "black") + 
+   
+   annotate("text", x=1.462, y=0.4, label="Cluster 1", size=5) + 
+   annotate("text", x=5.595, y=2.2, label="Cluster 2", size=5) + 
+   annotate("text", x=4.269, y=1.5, label="Cluster 3", size=5)
>  
> iris_k_means_x_y_plot_2

이상으로 K-중심 군집의 하나인 K-평균군집(K-means Clustering)에 대해서 알아보았습니다.

다음번 포스팅에서는 프로토타입 기반(Prototype-based) 군집분석의 두번째 기법으로 퍼지 군지(Fuzzy Clustering)에 대해서 알아보겠습니다.

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

[Reference]

(1) "Introduction to Data Mining", Pang-Ning Tan(Michigan State University), Michael Steinbach(University of Minnesota), Vipin Kumar(University of Minnesota), Addison-Wesley Companion Book

(2) "Clustering Algorithm", Ana Fred, INSTITUTO SUPERIOR TECNICO, Universidade Techica de Lisboa, 2009

(3) "R, SAS, MS-SQL을 활용한 데이터마이닝", 이정진 지음, 자유아카데미, 2011

(4) "Data Mining Cluster Analysis : Basic Concepts and Algorithms", Tan, Steinbach, Kumar, 2004

(5) Wikipedia :

    https://en.wikipedia.org/wiki/K-means_clustering

    https://en.wikipedia.org/wiki/Determining_the_number_of_clusters_in_a_data_set

응집형 계층적 군집화(agglomerative hierarchical clustering) 방법 중에서

 - 지난번 포스팅에서는 중심 연결법 (Centroid Linkage Method)를 다루었으며,

 - 이번 포스팅에서는 Ward 연결법 (Ward Linkage Method) 에 대해서 알아보겠습니다. 

(참고로, Ward는 이 방법을 제시했던 학자인 Joe H. Ward 의 이름을 딴 것입니다)

그동안의 응집형 계층적 군집화에서 다루었던 연결법들로 단일 연결법(single linkage), 완전 연결법(complete linkage), 평균 연결법(average linkage), 중심 연결법(centroid linkage)은 모두 '유클리드 제곱거리(euclidean squared distance)'에 기반해서 군집을 형성하는 방법들 이었습니다.

반면에, 이번에 소개하는 Ward 연결법(ward linkage)은 두 군집 간의 유사성을 두 군집이 합쳐졌을 때의 오차 제곱합(ESS : error sum of squares)의 증가분에 기반해서 측정합니다. (Similarity of two clusters is based on the increase in squared error when two clusters are merged).  즉, 거리 행렬(distance matrix)를 구할 때 오차제곱합의 증분(increase of ESS)을 두 군집 사이의 거리로 측정하게 됩니다.

단일 연결법이 노이즈나 이상치(noise and outlier)에 민감한 반면에, Ward 연결법은 노이즈나 이상치에 덜 민감한 장점이 있습니다. 그리고 Ward 연결법은 비슷한 크기의 군집끼리 묶어주는 경향이 있습니다.

중심 연결법과 Ward 연결법의 유사성 측정 수식이 비슷한데요, 중심 연결법의 유사성 측도 대비 Ward 연결법에는 가중값이 추가되었다는 점이 다릅니다. (분자는 같고, 분모가 서로 다름)

기본 컨셉은 다음번 포스팅의 프로토타입 기반 분할적 군집화의 하나인 "K 평균 군집화(K-means clustering)"와 유사한 측면이 있습니다.  

예제를 가지고 풀어보면서 설명을 이어가 보겠습니다.

(미리 말씀드리자면, 손으로 푼 결과와 R로 계산한 결과가 서로 다르게 나왔습니다. 어디서 차이가 생기는 건지 파악을 못했습니다.  혹시 이 블로그를 보시고 제가 잘못 계산한 부분을 찾으셨다면 댓글로 알려주시면 감사하겠습니다. ㅜ_ㅜ)

R script도 같이 제시하면서 설명하겠습니다.  먼저, 데이터 입력 및 plotting (↑ 위의 산점도 그래프) 입니다.

> ##--------------------------------------------
> ## (1) Agglomerative Hierarchical Clustering 
> ##   (b) Prototype-based
> ##    (1-5) Ward Linkage
> ##--------------------------------------------
> 
> x <- c(1, 2, 2, 4, 5)
> y <- c(1, 1, 4, 3, 4)
> 
> xy <- data.frame(cbind(x, y))
> 
> xy
  x y
1 1 1
2 2 1
3 2 4
4 4 3
5 5 4
> 
> # scatter plot of xy
> my_par = par(no.readonly = TRUE)
> par(mfrow = c(1, 1))
> plot(xy, pch = 19, xlab = c("x coordinate"), ylab = c("y coordinate"), 
+      xlim = c(0, 6), ylim = c(0, 6), 
+      main = "scatter plot of xy")
> 
> # adding student label
> text(xy[,1], xy[,2], labels = abbreviate(rownames(xy)), 
+      cex = 0.8, pos = 1, col = "blue") # pos=1 : at the bottom
> 
> 
> # adding dotted line
> abline(v=c(3), col = "gray", lty = 2) # vertical line
> abline(h=c(3), col = "gray", lty = 2) # horizontal line

유클리드 제곱거리를 구하는 R script 입니다. dist(xy, method="euclidean") 에다가 뒤에 "^2"를 붙여서 제곱을 했습니다.

> # proximity matrix : squares of euclidean distance matrix for 6 points
> dist(xy, method = "euclidean")^2
   1  2  3  4
2  1         
3 10  9      
4 13  8  5   
5 25 18  9  2

(여기까지는 단일 연결법, 완전 연결법, 평균 연결법, 중심 연결법과 동일합니다)

3) 군집(P1, P2)의 중심 구하기

[centroid coordinate of clusters - no.1]

아래의 부분집합그림에 보면 군집 (P1, P2) 의 중심(centroid) 위치에 노란색 별이 추가되었습니다.

> # Agglomerative Hierarchical Clustering : Centroid Linkage method
> hc_ward <- hclust(dist(xy)^2, method="ward.D")
> hc_ward

Call:
hclust(d = dist(xy)^2, method = "ward.D")

Cluster method   : ward.D 
Distance         : euclidean 
Number of objects: 5

> # dendrogram
> my_par = par(no.readonly = TRUE)
> par(oma = c(0, 0, 1, 0))
> par(mfrow = c(1, 2))
> plot(hc_ward)
> plot(hc_ward, hang = -1) # hang = -1 : line from the bottom

> # custering information
> rev(hc_ward)
$dist.method
[1] "euclidean"

$call
hclust(d = dist(xy)^2, method = "ward.D")

$method
[1] "ward.D"

$labels
NULL

$order
[1] 1 2 3 4 5

$height
[1]  1.000000  2.000000  8.666667 28.333333

$merge
     [,1] [,2]
[1,]   -1   -2
[2,]   -4   -5
[3,]   -3    2
[4,]    1    3

이번 포스팅에서는 여러개의 데이터 프레임을 한꺼번에 하나의 데이터프레임으로 묶는 몇가지 방법을 알아보고, 성능 측면을 비교해보겠습니다.

이번 포스팅은 andrew 님이 r-bloggers.com 에 썼던 글을 그대로 가져다가 번역을 한 내용입니다. 

[ Source ] Concatenating a list of data frames , June 6, 2014, By andrew

* source : http://www.r-bloggers.com/concatenating-a-list-of-data-frames/

결론 먼저 말씀드리면, data.table package의 rbindlist(data) 함수가 속도 면에서 월등히 빠르네요.

[ R로 여러개의 데이터프레임을 한꺼번에 하나의 데이터프레임으로 묶기 ]

0) 문제 (The problem)

아래처럼 3개의 칼럼으로 구성된 100,000 개의 자잘한 데이터 프레임을 한개의 커다란 데이터 프레임으로 합치는 것이 풀어야 할 문제, 미션입니다.

data = list() 로 해서 전체 데이터 프레임들을 data라는 리스트로 만들어서 아래 각 방법별 예제에 사용하였습니다.

> ###########################################
> ## Concatenating a list of data frames
> ## do.call(rbind, data)
> ## ldply(data, rbind)
> ## rbind.fill(data)
> ## rbindlist(data) ** winner! **
> ###########################################
> 
> ## The Problem
> 
> data = list()
> 
> N = 100000
> 
> for (n in 1:N) {
+   data[[n]] = data.frame(index = n, 
+                          char = sample(letters, 1), 
+                          z = runif(1))
+ }
> 
> data[[1]]
  index char         z
1     1    j 0.2300154

1) The navie solution : do.call(rbind, data)

가장 쉽게 생각할 수 있는 방법으로 base package에 포함되어 있는 rbind() 함수를 do.call 함수로 계속 호출해서 여러개의 데이터 프레임을 위/아래로 합치는 방법입니다. 

이거 한번 돌리니 정말 시간 오래 걸리네요.  @@~  낮잠 잠깐 자고 와도 될 정도로요.

> ## (1) The Naive Solution
> head(do.call(rbind, data))
  index char          z
1     1    j 0.23001541
2     2    f 0.63555284
3     3    d 0.65774397
4     4    y 0.46550511
5     5    b 0.02688307
6     6    u 0.19057217

2-1) plyr package : ldply(data, rbind)

두번째 방법은 plyr package의 ldply(data, rbind) 함수를 사용하는 방법입니다.

> ## (2) Alternative Solutions #1 and #2
> ## (2-1) plyr package : ldply(data, rbind)
> install.packages("plyr")

> library(plyr)

> head(ldply(data, rbind))
  index char          z
1     1    j 0.23001541
2     2    f 0.63555284
3     3    d 0.65774397
4     4    y 0.46550511
5     5    b 0.02688307
6     6    u 0.19057217

2-2) plyr package : rbind.fill(data)

세번째 방법은 plyr package의 rbind.fill(data) 함수를 사용하는 방법입니다.  결과는 앞의 두 방법과 동일함을 알 수 있습니다.

> ## (2-2) plyr package : rbind.fill(data)
> library(plyr)
> head(rbind.fill(data))
  index char          z
1     1    j 0.23001541
2     2    f 0.63555284
3     3    d 0.65774397
4     4    y 0.46550511
5     5    b 0.02688307
6     6    u 0.19057217

3) data.table package : rbindlist(data)

마지막 방법은 data.table package의 rbindlist(data) 함수를 사용하는 방법입니다.

> ## (3) Alternative Solution 
> ## data.table package : rbindlist(data)
> install.packages("data.table")

> library(data.table)

> head(rbindlist(data))
   index char          z
1:     1    j 0.23001541
2:     2    f 0.63555284
3:     3    d 0.65774397
4:     4    y 0.46550511
5:     5    b 0.02688307
6:     6    u 0.19057217

4) 벤치마킹 테스트 (bechmarking test)

> ## Benchmarking (performance comparison)
> install.packages("rbenchmark")

> library(rbenchmark)
> benchmark(do.call(rbind, data),
+           ldply(data, rbind), 
+           rbind.fill(data), 
+           rbindlist(data))

패키지/함수별 성능 비교를 해본 결과 data.table 패키지의 rbindlist(data) 함수가 월등히 빠르다는 것을 알 수 있습니다.  위의 벤치마킹 결과를 보면, 속도가 가장 빨랐던 rbindlist(data)를 1로 놨을 때, 상대적인 속도(relative 칼럼)를 보면 rbind.fill(data)가 86.932로서 rbindlist(data)보다 86배 더 오래걸리고, ldply(data, rbind)가 292.644로서 rbindlist(data)보다 292배 더 오래걸린다는 뜻입니다.  do.call(rbind, data)는 rbindlist(data) 보다 상대적으로 668.692배 더 시간이 걸리는 것으로 나오네요.

rbindlist(data)가 훨등히 속도가 빠른 이유는 두가지인데요,

(1) rbind() 함수가 각 데이터 프레임의 칼럼 이름을 확인하고, 칼럼 이름이 다를 경우 재정렬해서 합치는데 반해, data.table 패키지의 rbindlist() 함수는 각 데이터 프레임의 칼럼 이름을 확인하지 않고 단지 위치(position)를 기준으로 그냥 합쳐버리기 때문이며,

(따라서, rbindlist() 함수를 사용하려면 각 데이터 프레임의 칼럼 위치가 서로 동일해야 함)

(2) rbind() 함수는 R code로 작성된 반면에, data.table 패키지의 rbindlist() 는 C 언어로 코딩이 되어있기 때문입니다.

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

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

큰 자리수의 숫자를 불러오기하면 지수 표기법(exponential notation, scientific notation)으로 되어 있어서 읽기가 힘든 경우가 있습니다.

가령 '900000' 이라는 숫자를 R에서 읽어들였더니 '9e+5'라고 쓰여있다면 좀 어색하지요? 

이럴 경우 '9e+5'라고 프린트되어 있는 표기를 '900000'처럼 고정형 숫자 표기법(fixed notation)으로 바꿔주고 싶을 때 R에서 사용하는 방법 2가지를 소개하고자 합니다.

(1) 하나는 전역적으로(global environment) 영향을 미치는 옵션으로 : options(scipen = xxx) 이구요,

(2) 또 하나는 국소적으로(local) 특정 변수에 대해서 설정을 하는 옵션으로 : format(var, scientific = FALSE)

입니다.

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

가령, "2016년 7월 9일 10시 30분 50초"의 시간을 나타내는 숫자가 "20160709103050" 처럼 입력되어 있는 데이터가 있다고 칩시다.  이를 options("scipen" = -100) 또는 options(scipen = -100)인 상태에서 입력을 받고서 프린트를 해보면 "2.016071e+13" 이라고 나옵니다.

(괄호 안의 숫자는 충분히 큰 숫자를 써주면 됩니다. options("scipen" = -30) 이라고 해도 동일하게 지수형태로 표기될겁니다)

> #################################################
> ## changing from exponential to numeric notation
> ## options("scipen" = 100)
> ## format(df$var, scientific = FALSE)
> #################################################
> 
> options("scipen" = -100)
> number <- c(20160709103050)
> number
[1] 2.016071e+13
> str(number)
 num 2.02e+13

이걸 options("scipen" = 100) 또는 options(scipen = 100)이라고 하고 다시 number를 프린트해보면 아래처럼 "20160709103050" 으로 원하는 형태로 잘 나옵니다. (괄호 안의 숫자는 충분히 큰 숫자를 써주면 됩니다.)

이 옵션은 R 사용환경 전체에 영향을 미칩니다.  즉, 특정 변수는 지수형 표기, 특정 변수는 숫자표기로 할 수는 없다는 뜻입니다.

> options("scipen" = 100)
> number
[1] 20160709103050 

참고로, options(scipen = x) 의 default 값은 '0' 입니다.

날짜/시간 얘기가 나온 김에 조금 더 진도를 나가 보자면요,

시계열분석이나 년/월/일/시간/분/초 단위로 구분해서 연산을 해야 하는 경우에 as.POSIXlt() 함수로 날짜/시간 type으로 바꾸어주면 편합니다. 

아래에 as.POSIXlt로 숫자형 number를 날짜/시간 유형으로 바꿔보려고 했더니 "NA" 값이 나오네요. 

> datetime_1 <- as.POSIXlt(number, # numeric
+                          format = '%Y%m%d%H%M%S', 
+                          origin = "1970-01-01", 
+                          tz ="UTC")
> 
> datetime_1 # NA
[1] NA

이럴경우 특정 변수만 format(var, scientific = FALSE) 옵션을 사용해서 포맷을 바꾸어주면 numeric이었던 것이 character로 바뀌었고, 큰 따옴표도 생겼습니다.

> # scientific = FALSE
> number_2 <- format(number, scientific = FALSE)
> number_2
[1] "20160709103050"
> str(number_2)
 chr "20160709103050"

이렇게 문자형으로 바꾸고 나서 다시 as.POSIXlt(char_var, format = '%Y%m%d%H%M%S', origin="1970-01-01", tz = UTC) 로 날짜/시간 유형으로 변환하면 잘 작동하네요.

> datetime_2 <- as.POSIXlt(number_2, # character
+                          format = '%Y%m%d%H%M%S', 
+                          origin = "1970-01-01", 
+                          tz ="UTC") # UTC : universal time
> 
> datetime_2 # "2016-07-09 10:30:50 UTC"
[1] "2016-07-09 10:30:50 UTC"

format() 대신에 as.character(numeric_var) 를 사용해서 숫자형을 문자형으로 바꾸어주고 나서, as.POSIXlt(char_var, format = '%Y%m%d%H%M%S', origin="1970-01-01", tz = UTC) 로 날짜/시간 유형으로 변환해도 똑같이 잘 작동합니다.

> # as.character
> number_3 <- as.character(number)
> number_3
[1] "20160709103050"
> str(number_3)
 chr "20160709103050"
> 
> datetime_3 <- as.POSIXlt(number_3, # character
+                          format = '%Y%m%d%H%M%S', 
+                          origin = "1970-01-01", 
+                          tz ="UTC") # UTC : universal time
> 
> datetime_3 # "2016-07-09 10:30:50 UTC"
[1] "2016-07-09 10:30:50 UTC" 

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

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

이번 포스팅에서는 방향성 있는 가중 네트워크 시각화 (directed and weighted network visualization)에 대해서 소개하겠습니다.

 - R igraph package 를 사용해서

 - 심리학 저널의 상호인용 빈도 데이터

  (* source : Lattin, Carroll and Green, 2003)

를 대상으로 논문 상호인용 빈도 네트워크를 그려보겠습니다.

먼저, 방향성 있는 가중 네트워크에 대해서 한번 더 소개하자면, 아래와 같이 화살표로 방향이 있고, 연결선에 가중치(weight)가 있는 네트워크를 말합니다.

[ 방향성 있는 가중 네트워크 (Directed and Weighted Network) ]

심리학 저널의 상호인용 빈도를 행렬 형태로 나타내보면 아래와 같습니다.

행렬의 요소 (i, j) 빈도는 저널 i의 논문이 저널 j에 인용된 회수입니다.

 - 예) (1, 2) = 32 : AJP 논문 32편이 JASP 논문에 인용됨

        (2, 1) = 8 : JASP 논문 8편이 AJP 논문에 인용됨

[ 심리학 저널 간 상호 인용 빈도 ]

library()함수로 R의 igraph package를 로딩해서 시각화를 해보겠습니다.

igraph package는 데이터 포맷이 ID1, ID2, Weight or Frequency 의 행렬(matrix) 입니다.

아래에 위의 표를 행렬로 입력해서 불러들였는데요, 빈도(3번째 열)에 +1 을 해주었습니다.

(0이면 에러가 나는지 안되네요)

R script는 'R을 활용한 사회네트워크분석 입문'의 예시를 거의 대부분 사용하였으며,

그래프의 설정 parameter를 이리 저리 숫자를 조절해보면서 살짝 바꿔보았습니다.

> #############################################
> ## Network Analysis - igraph package
> #############################################
> library(igraph)
> 
> psych_edgelist <- matrix(
+   c(0, 1, 32, 0, 2, 2, 0, 3, 35, 0, 4,  6, 0, 5, 4, 
+     0, 6, 125, 0, 7,  2, 1, 0, 8, 1, 2, 8, 1, 3, 8, 
+     1, 4, 116, 1, 5,  9, 1, 6,  19, 1, 7, 5, 2, 0, 4, 
+     2, 1, 16, 2, 3, 0, 2, 4,  11, 2, 5, 7, 2, 6, 6,
+     2, 7, 5, 3, 0,  21, 3, 1, 11, 3, 2, 1, 3, 4, 1, 
+     3, 5, 0, 3, 6, 70, 3, 7, 0, 4, 0, 0, 4, 1, 73, 
+     4, 2, 7, 4, 3,  0, 4, 5, 3, 4, 6,  0, 4, 7, 13, 
+     5, 0, 1, 5, 1,  9, 5, 2,  8, 5, 3,  1, 5, 4, 7,
+     5, 6, 0, 5, 7,  2, 6, 0,  85, 6, 1, 119, 6, 2, 16, 
+     6, 3, 126, 6, 4,  12, 6, 5, 27, 6, 7, 13, 7, 0, 2, 
+     7, 1, 4, 7, 2,  10, 7, 3, 1, 7, 4,  7, 7, 5, 5, 7,6, 15), 
+   byrow = T, ncol = 3) + 1
> 
> psych_edgelist
      [,1] [,2] [,3]
 [1,]    1    2   33
 [2,]    1    3    3
 [3,]    1    4   36
 [4,]    1    5    7
 [5,]    1    6    5
 [6,]    1    7  126
 [7,]    1    8    3
 [8,]    2    1    9
 [9,]    2    3    9
[10,]    2    4    9
[11,]    2    5  117
[12,]    2    6   10
[13,]    2    7   20
[14,]    2    8    6
[15,]    3    1    5
[16,]    3    2   17
[17,]    3    4    1
[18,]    3    5   12
[19,]    3    6    8
[20,]    3    7    7
[21,]    3    8    6
[22,]    4    1   22
[23,]    4    2   12
[24,]    4    3    2
[25,]    4    5    2
[26,]    4    6    1
[27,]    4    7   71
[28,]    4    8    1
[29,]    5    1    1
[30,]    5    2   74
[31,]    5    3    8
[32,]    5    4    1
[33,]    5    6    4
[34,]    5    7    1
[35,]    5    8   14
[36,]    6    1    2
[37,]    6    2   10
[38,]    6    3    9
[39,]    6    4    2
[40,]    6    5    8
[41,]    6    7    1
[42,]    6    8    3
[43,]    7    1   86
[44,]    7    2  120
[45,]    7    3   17
[46,]    7    4  127
[47,]    7    5   13
[48,]    7    6   28
[49,]    7    8   14
[50,]    8    1    3
[51,]    8    2    5
[52,]    8    3   11
[53,]    8    4    2
[54,]    8    5    8
[55,]    8    6    6
[56,]    8    7   16
> 
> psych.w <- graph.edgelist(psych_edgelist[, 1:2])
> psych.w
IGRAPH D--- 8 56 -- 
+ edges:
 [1] 1->2 1->3 1->4 1->5 1->6 1->7 1->8 2->1 2->3 2->4 2->5 2->6
[13] 2->7 2->8 3->1 3->2 3->4 3->5 3->6 3->7 3->8 4->1 4->2 4->3
[25] 4->5 4->6 4->7 4->8 5->1 5->2 5->3 5->4 5->6 5->7 5->8 6->1
[37] 6->2 6->3 6->4 6->5 6->7 6->8 7->1 7->2 7->3 7->4 7->5 7->6
[49] 7->8 8->1 8->2 8->3 8->4 8->5 8->6 8->7
> 
> E(psych.w)$weight <- psych_edgelist[,3]
> psych.w
IGRAPH D-W- 8 56 -- 
+ attr: weight (e/n)
+ edges:
 [1] 1->2 1->3 1->4 1->5 1->6 1->7 1->8 2->1 2->3 2->4 2->5 2->6
[13] 2->7 2->8 3->1 3->2 3->4 3->5 3->6 3->7 3->8 4->1 4->2 4->3
[25] 4->5 4->6 4->7 4->8 5->1 5->2 5->3 5->4 5->6 5->7 5->8 6->1
[37] 6->2 6->3 6->4 6->5 6->7 6->8 7->1 7->2 7->3 7->4 7->5 7->6
[49] 7->8 8->1 8->2 8->3 8->4 8->5 8->6 8->7
> 
> psych.diag <- c(119, 510, 84, 533, 225, 52, 586, 58)
> 
> psych.name <- c("AJP", "JASP", "JAP", "JCPP", "JCP", "JEDP", "JEXP", "PKA")
> 
> plot(psych.w, 
+      layout = layout.circle, 
+      vertex.size = 2, 
+      vertex.shape = "none", 
+      vertex.size = psych.diag, 
+      vertex.label = psych.name, 
+      vertex.label.font = 2, 
+      vertex.label.cex = sqrt(psych.diag)/10, 
+      edge.width=2 + E(psych.w)$weight/10, 
+      edge.arrow.width = E(psych.w)$weight/100
+ )

igraph plot 의 parameter 기능은 아래와 같습니다. 

parameter를 바꿔가면서 최적의 이쁜 그래프를 찾아가 보시기 바랍니다.

1) Layout : 점의 좌표를 정하는 알고리즘

 - layout.circle : 원 배치 (위의 예시)

 - layout.random : 무작위 배치

- layout.fruchterman.reingold : Fruchterman Reingold 배치 알고리즘

 - layout.kamada.kawai : kamada Kawai 배치 알고리즘

 - layout.lgl : 대규모 네트워크를 위한 배치 알고리즘

2) Edge : 선 관련 파라미터

- edge.color : 선 색 지정 (default = "darkgrey")

- edge.width : 선 폭

- edge.arrow.size : 화살 크기

- edge.arrow.width : 화살 폭

- edge.arrow.mode : 화살 머리 유형 (0 : 없음,  1 : 역방향,  2 : 순방향,   3 : 양방향)

                            (무방향 네트워크의 경우 default = 0)

- edge.lty : 선 유형 ("solid", "dashed", "dotted", "dotdash", "longdash", "twodash")

- edge.label : 선 레이블

- edge.label.family : 선 레이블 종류 ("serif", "sans", "mono" 등)

- edge.label.font : 선 레이블 글자형 (1 : plain text, 2 : bold, 3 : italic, 4 : bold italic)

- edge.label.cex : 선 레이블 크기 (default = 1)

- edge.label.color : 선 레이블 색 (default = "navy")

3) Vertex : 점 관련

- vertex.size : 점 크기, vector도 가능 (default = 15)

- vertex.color : 점 색 (default = "SkyBlue2")

- vertex.frame.color : 점 윤곡의 색 (default = "black")

- vertex.shape : 점 형태 ("circle", "square", "rectangle", "none", default = "circle")

- vertex.label : 점 레이블 (vector)

- vertex.label.family : 점 레이블 종류 ("serif", "sans", "mono" 등)

- vertex.label.font : 점 레이블 글자형 (1 : plain text, 2 : bold, 3 : italic, 4 : bold italic)

- vertex.label.cex : 점 레이블 크기 (default = 1)

- vertex.label.dist : 점 중심과 레이블 간 거리 (default = 0)

- vertex.label.degree : 점 레이블 방향(radian) (좌 : 0, 우 : pi, 상 : -pi/2, 하 : pi/2)

- vertex.label.color : 점 레이블 색 (default = "navy")

네트워크 그래프 하나 그리는데 무슨 놈의 파라미터가 종류가 이렇게 많은 건지 놀랍기도하고, 이걸 언제 다 설정하나 부담되기도 하지요?

default 설정 값을 이용해서 한번 얼른 그려보시고요, 원하는 모양이 아니다 싶으면 parameter 종류 중에서 살짝 살짝 손을 좀 봐가면서 몇 번 더 그래프를 그려보시기 바랍니다. 

R이 제공하는 이런 다양한 그래프 옵션이면 못할 것이 없겠지요?!

[Reference]

- R을 활용한 사회네트워크분석 입문, 허명회 저, 자유아카데미, 2012

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

이번 포스팅에서는 방향을 가진 네트워크 상에서 유량, 에너지, 물자 등의 '네트워크(network) 기반'의 흐름(flow)을 시각화하는데 유용한 Sankey Diagram 에 대해서 알아보겠습니다. 

Sankey Diagram 은 얼핏보면 평행좌표그림(parallel coordinate plot)과 비슷한 면이 있습니다.  하지만 Sankey Diagram은 Node -> Edge 로의 관계, 경로가 있다는 점, 경로의 두께를 weight 에 따라서 다르게 한다든지, 색깔을 부여할 때도 반투명(alpha)하게 한다든지 해서 평행좌표그림 대비 다릅니다. (평행좌표그림은 그냥 다변량 변수들을 y축 높이를 같게 해서 옆으로 x를 죽 늘어놓은 형태. 경로 개념 없음. 선 두께 동일) 

아래 이미지는 구글에서 Sankey Diagram 이라고 키워드 검색했을 때 나오는 이미지들을 화면캡쳐한 것입니다.  아래 이미지를 보면 '아, 이거~ ' 싶으시죠?

* 이미지 출처 : Google 검색

자, 그럼

 - R의 riverplot package(* author : January Weiner) 를 가지고

 - minard 데이터셋(* source : Charles Joseph Minard)을 사용해서
   (R riverplot package에 내장되어 있음)

나폴레옹 군대가 러시아로 진군했다가 퇴각했던 경로를 시각화해보겠습니다.

분석의 재미를 더하기 위해 "minard" 데이터셋에 대한 역사적인 배경을 간략히 알아보고 넘어가겠습니다.

때는 바야흐로 1812년, 유럽이 되겠습니다.  프랑스의 나폴레옹은 유럽의 상당 국가를 점령했으며, 영국을 침탈하기 위해 백방으로 쌈질을 걸었지만 번번히 실패하고 있던 상황이었습니다.  이에 화가 난 나폴레옹은 '바다를 건너가서 싸우는 게 승산이 낮으니 차라리 유럽 본토와 영국과의 무역을 봉쇄해서 영국의 피를 말리자. 그래, 바로 이거야. 내 사전에 불가능은 없어!' 라는 계획을 세우게 됩니다.

 

그러나 but,

 

나폴레옹의 프랑스가 이미 과도하게 힘이 세졌기 때문에 견제가 필요하다가 생각한 러시아의 알렉산더 짜르는 '나폴레옹, 내가 니 봉이냐?  영국 다음에 러시아 공격할거지?  내가 누구 좋으라고 니 말을 따라?' 라면서 나폴레옹의 영국과의 무역 폐쇄령을 쌩까고 게기게 됩니다. 

 

이에 발끈한 나폴레옹은 프랑스 대군을 모집해서 러시아 짜르의 못되고 괴씸한(?) 버릇을 고쳐주고자, 본떼를 보여주고자 1812년 10월, 겨울이 코앞인 시점에 러시아 모스코바로 진격을 하게 됩니다. 이때만 해도 정말 나폴레옹 군대는 "진격의 거인" 이었습니다.

 

그러나 but,

 

러시아 알렉산더 짜르가 대책없이 나폴레옹에게 대든게 아니었습니다.  러시아 짜르는 러시아 군과 국민에게 아주 간단한(?) 전략의 명령을 내립니다.  "프랑스군이 진격하는 곳의 모든 것을 태워서 프랑스군이 아무것도 탈취하지 못하도록 하고, 싸우지는 말고 퇴각하라.  전투 전략 끝!"

 

아마, 나폴레옹의 프랑스군은 처음 며칠은 러시아로 무혈입성하는 것에 신이 났을지도 모릅니다.

 

그러나 but,

 

10월이 11월이 되고, 그 담에 12월이 되면서 추위와 배고픔과 질병에 프랑스 군인들의 대부분이 죽어나갔습니다.  진격의 거인 나폴레옹은 러시아 짜르의 "불태우고 후퇴" 전략에 속수무책으로 당하면서 거의 전멸을 당하게 됩니다.

 

아래의 그림이 나폴레옹 프랑스 군대가 러시아로 진격했다가 퇴각한 진로를 지도 상에 표기한 것입니다.

* 출처 : https://robots.thoughtbot.com/analyzing-minards-visualization-of-napoleons-1812-march

여기까지의 이야기를 토대로 숫자와 그래프를 가지고 한눈에 실감할 수 있는 시각화를 Charles Joseph Minard 이라는 분이 아래와 같이 했습니다.

[ Napoleon army march : minard ]

* 출처 : https://robots.thoughtbot.com/analyzing-minards-visualization-of-napoleons-1812-march

위의 2차원의 minard visualization에는 다양한 차원의 정보가 알차게(!) 들어있는데요,

 - (1) 프랑스군의 진격과 퇴각 경로 (advance and retreat path and direction)

        : 연한 색깔이 진격(advance), 검정 색깔이 퇴각(retreat)

        : 도시 이름은 좌표에 따라 text로 표기

        : 프랑스 군이 세갈래(하나의 큰 줄기, 두 개의 얇은 줄기)로 나누어 진격한 것도 선이 갈라지게 표현

 - (2) 프랑스군의 인명 손실 규모 (loss of life at a time and location)
        : 선의 두께(line width), 처음에는 몽둥이처럼 두꺼웠던 선이 퇴각 마무리 시점에는 실처럼 가늘게 됨

 - (3) 온도 (temperature)

        : 그림 하단에 퇴각 시점의 온도를 그리 넣음.  최저 -30도씨까지 떨어졌음.  얼어죽기 딱 좋은 날씨. -_-;

 - (4) 강 (river)

        : 하단에 얇은 수직 선으로 나폴레옹군이 퇴각 시점에 맞닥트려 시련을 더해 준 강(river)을 그려 넣음. 

이렇게 많은 알찬 정보를 저 위의 시각화 하나에 오롯히 담아 냈습니다!!!  이해하기 쉽죠?!

이걸 데이터로 나타내 보면 아래와 같습니다.  @@;  

이게 무슨 소린가, 데이터가 뭘 말해주려고 하나.... 눈 돌아가지요? 

위의 그래프로 보면 단박에 이해되는 것을 아래의 숫자로 보면 한숨 나오고 갑갑하지요? ㅋㅋ

> install.packages("riverplot")

> library(riverplot)

> data( minard )
> minard
$edges
    ID1  ID2  Value direction
1    A1   A2 422000         A
2    A2   A3 400000         A
3    A3   A4 320000         A
4    A4   A5 320000         A
5    A5   A6 300000         A
6    A6   A7 280000         A
7    A7   A8 240000         A
8    A8   A9 210000         A
9    A9  A10 180000         A
10  A10  A11 175000         A
11  A11  A12 145000         A
12  A12  A13 140000         A
13  A13  A14 127100         A
14  A14  A15 100000         A
15  A15  A16 100000         A
16  A16   R1 100000         A
17   R1   R2 100000         R
18   R2   R3  98000         R
19   R3   R4  97000         R
20   R4   R5  96000         R
21   R5   R6  87000         R
22   R6   R7  55000         R
23   R7   R8  37000         R
24   R8   R9  24000         R
25   R9  R10  20000         R
26  R10  R11  50000         R
27  R11  R12  50000         R
28  R12  R13  48000         R
29  R13  R14  20000         R
30  R14  R15  12000         R
31  R15  R16  14000         R
32  R16  R17   8000         R
33  R17  R18   4000         R
34  R18  R19  10000         R
35  R19  R19  10000         R
36   A3 A4.2  60000         A
37 A4.2 A5.2  40000         A
38 A5.2 A6.2  33000         A
39 A6.2  R10  30000         R
40   A2 A3.1  22000         A
41 A3.1  R18   6000         A

$nodes
       ID Longitude Latitude
A1     A1      24.0     54.9
A2     A2      24.5     55.0
A3     A3      25.5     54.5
A4     A4      26.0     54.7
A5     A5      27.0     54.8
A6     A6      28.0     54.9
A7     A7      28.5     55.0
A8     A8      29.0     55.1
A9     A9      30.0     55.2
A10   A10      30.3     55.3
A11   A11      32.0     54.8
A12   A12      33.2     54.9
A13   A13      34.4     55.5
A14   A14      35.5     55.4
A15   A15      36.0     55.5
A16   A16      37.6     55.8
R1     R1      37.7     55.7
R2     R2      37.5     55.7
R3     R3      37.0     55.0
R4     R4      36.8     55.0
R5     R5      35.4     55.3
R6     R6      34.3     55.2
R7     R7      33.3     54.8
R8     R8      32.0     54.6
R9     R9      30.4     54.4
R10   R10      29.2     54.3
R11   R11      28.5     54.2
R12   R12      28.3     54.3
R13   R13      27.5     54.5
R14   R14      26.8     54.3
R15   R15      26.4     54.4
R16   R16      25.0     54.4
R17   R17      24.4     54.4
R18   R18      24.2     54.4
R19   R19      24.1     54.4
A4.2 A4.2      26.6     55.7
A5.2 A5.2      27.4     55.6
A6.2 A6.2      28.7     55.5
A3.1 A3.1      24.6     55.8

$cities
   Longitude Latitude           Name
1       24.0     55.0          Kowno
2       25.3     54.7          Wilna
3       26.4     54.4       Smorgoni
4       26.8     54.3      Moiodexno
5       27.7     55.2      Gloubokoe
7       28.5     54.3     Studienska
8       28.7     55.5        Polotzk
9       29.2     54.4           Bobr
10      30.2     55.3        Witebsk
11      30.4     54.5         Orscha
13      32.0     54.8       Smolensk
14      33.2     54.9    Dorogobouge
15      34.3     55.2          Wixma
16      34.4     55.5          Chjat
17      36.0     55.5        Mojaisk
18      37.6     55.8         Moscou
19      36.6     55.3      Tarantino
20      36.5     55.0 Malo-Jarosewii

* source : R riverplot package, author : January Weiner, data source : Charles Joseph Minard

 Sankey Diagram에서 사용하는 node, edge라는 용어를 이해하기 위해서, 두 개체 간 쌍을 이룬 관계 (mathematical structures used to model pairwise relations between objects)를 다루는 Graph Theory에 대해서 간략히 짚고 넘어가겠습니다.

아래 그래프처럼 점(Node or Point or Vertice)과 선(Edge or Link or Line or Arc)으로 개체 간의 관계를 나타내는 그래프로 나타내어 연구하는 수학, 컴퓨터 과학 분야가 Graph theory입니다.  최적화(optimization) 할 때도 네트워크 그래프 많이 쓰곤 합니다.

Sankey Diagram 그리려면

  - Nodes : 개체들의 ID, Longitude, Latitude, Labels

  - Edges : 개체 간 관계를 나타내는 ID 1, ID 2, Value (or weight) 

정보가 필요합니다.

> str(minard)
List of 3
 $ edges :'data.frame':	41 obs. of  4 variables:
  ..$ ID1      : chr [1:41] "A1" "A2" "A3" "A4" ...
  ..$ ID2      : chr [1:41] "A2" "A3" "A4" "A5" ...
  ..$ Value    : num [1:41] 422000 400000 320000 320000 300000 280000 240000 210000 180000 175000 ...
  ..$ direction: Factor w/ 2 levels "A","R": 1 1 1 1 1 1 1 1 1 1 ...
 $ nodes :'data.frame':	39 obs. of  3 variables:
  ..$ ID       : chr [1:39] "A1" "A2" "A3" "A4" ...
  ..$ Longitude: num [1:39] 24 24.5 25.5 26 27 28 28.5 29 30 30.3 ...
  ..$ Latitude : num [1:39] 54.9 55 54.5 54.7 54.8 54.9 55 55.1 55.2 55.3 ...
 $ cities:'data.frame':	18 obs. of  3 variables:
  ..$ Longitude: num [1:18] 24 25.3 26.4 26.8 27.7 28.5 28.7 29.2 30.2 30.4 ...
  ..$ Latitude : num [1:18] 55 54.7 54.4 54.3 55.2 54.3 55.5 54.4 55.3 54.5 ...
  ..$ Name     : Factor w/ 20 levels "Bobr","Chjat",..: 5 18 15 9 4 16 13 1 19 12 

R riverplot package를 사용해서 드디어 Sankey diagram을 그려보겠습니다.  R script는 riverplot package(* author : January Weiner)에 있는 예제 script를 그대로 인용하였습니다. (날로 먹는 듯한 이 기분..^^;) 

> ############################
> # Sankey diagram 
> # using R riverplot package 
> # minard data (list format)
> ############################
> 
> # install.packages("riverplot")
> library(riverplot)
> data( minard )
> 
> nodes <- minard$nodes
> edges <- minard$edges
> colnames( nodes ) <- c( "ID", "x", "y" )
> colnames( edges ) <- c( "N1", "N2", "Value", "direction" )
> 
> 
> # color the edges by troop movement direction
> edges$col <- c( "#e5cbaa", "black" )[ factor( edges$direction ) ]
> 
> # color edges by their color rather than by gradient between the nodes
> edges$edgecol <- "col"
> 
> # generate the riverplot object and a style
> river <- makeRiver( nodes, edges )
> style <- list( edgestyle= "straight", nodestyle= "invisible" )
> 
> # plot the generated object
> plot( river, lty= 1, default_style= style )
> 
> # Add cities
> with( minard$cities, points( Longitude, Latitude, pch= 19 ) )
> with( minard$cities, text( Longitude, Latitude, Name, adj= c( 0, 0 ) ) )
> # Add title
> title("Sankey Diagram - Napoleon army march, minard")

다음번 포스팅에서는 igraph package를 사용해서 방향성 있는 가중 네트워크를 시각화해보겠습니다.

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

[Reference]

 - R riverplot package manual : https://cran.r-project.org/web/packages/riverplot/riverplot.pdf

 - On minards visualization : https://robots.thoughtbot.com/analyzing-minards-visualization-of-napoleons-1812-march

폴더에 자잘하게 쪼개진 여러개의 파일들이 있을 때, 그리고 이 파일들을 일일이 R로 불러오기 해야 할 때, 더그리고 이들 불러온 파일을 한개의 데이터셋을 합쳐야 할 때 (이쪽 동네 전문용어로) 노가다를 하지 않고 좀더 스마트하게 하는 방법을 소개하겠습니다.

순서는 다음과 같습니다.

• (1) 폴더 경로 객체로 만들기
• (2) 폴더 내 파일들 이름을 list-up 하여 객체로 만들기
• (3) 파일 개수 객체로 만들기
• (4) 폴더 내 파일들을 LOOP 돌려서 불러오기 : read.table()
• (5) 파일을 내보내면서 합치기 : write.table(dataset, APPEND = TRUE)
• (6) 데이터프레임으로 불러오기, 칼럼 이름 넣기 : read.table(dataset_all, col.names = c())

자, 예를 들면서 순서대로 R script 설명하겠습니다.

아래의 화면캡쳐 예시처럼 MyDocuments > R > FILES 폴더 아래에 daily로 쪼개진 10개의 text 파일들이 들어있다고 해봅시다.  (10개 정도야 일일이 불어올 수도 있겠지만, 100개, 1,000개 파일이 들어있다면?)

• (1) 폴더 경로 객체로 만들기

## cleaning up environment
rm(list=ls())

## making directory as an object 
src_dir <- c("C:/Users/Owner/Documents/R/FILES") # 경로 구분 : '\'를 '/'로 바꿔야 함 

src_dir 
#[1] "C:/Users/Owner/Documents/R/FILES"

• (2) 폴더 내 파일들 이름을 list-up 하여 객체로 만들기 : list.files()

# listing up name of files in the directory => object 
src_file <- list.files(src_dir) # list 

src_file 
#[1] "day_20160701.txt" "day_20160702.txt" "day_20160703.txt" "day_20160704.txt" 
#[5] "day_20160705.txt" "day_20160706.txt" "day_20160707.txt" "day_20160708.txt" 
#[9] "day_20160709.txt" "day_20160710.txt"

"C:/Users/Owner/Documents/R/FILES" 디렉토리에 들어있는 파일들을 열어보면 아래와 같은 데이터들이 들어있습니다. (가상으로 만들어 본 것임)  daily로 집계한 데이터들이 들어있네요.

• (3) 파일 개수 객체로 만들기 : length(list)

# counting number of files in the directory => object 
src_file_cnt <- length(src_file)

src_file_cnt 
#[1] 10

여기까지 R을 실행하면 아래와 같이 environment 창에 객체들이 생겼음을 확인할 수 있습니다.

• (4) 폴더 내 파일들을 LOOP 돌려서 불러오기
  => (5) 파일을 내보내면서 합치기 : write.table(dataset, APPEND = TRUE)
  
  : for(i in 1:src_file_cnt) {read.table()
                                   write.table(dataset, append = TRUE)}

## write.table one by one automatiically, using loop program 
for(i in 1:src_file_cnt) {
	# write.table one by one automatiically, using loop program 
    dataset <- read.table(
    	paste(src_dir, "/", src_file[i], sep=""), 
        sep=",", 
        header=F, 
        stringsAsFactors = F) 
        
    # dataset exporting with 'APPEND = TREU' option, filename = dataset_all.txt
    write.table(dataset, 
    	paste(src_dir, "/", "dataset_all.txt", sep=""), 
        sep = ",", 
        row.names = FALSE, 
        col.names = FALSE, 
        quote = FALSE, 
        append = TRUE) # appending dataset (stacking)
        
    # delete seperate datasets
    rm(dataset) 
    
    # printing loop sequence at console to check loop status
    print(i)
} 

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

여기까지 실행을 하면 아래처럼 MyDocuments>R>FILES 폴더 아래에 'dataset_all.txt' 라는 새로운 텍스트 파일이 하나 생겼음을 확인할 수 있습니다. 

새로 생긴 'dataset_all.txt' 파일을 클릭해서 열어보면 아래와 같이 'day_20160701.txt' ~ 'day_20160710.txt'까지 10개 파일에 흩어져있던 데이터들이 차곡차곡 쌓여서 합쳐져 있음을 확인할 수 있습니다.

• (6) 데이터 프레임으로 불러오기 : read.table()
       칼럼 이름 붙이기 : col.names = c("var1", "var2", ...)

# reading dataset_all with column names 
dataset_all_df <- read.table(
	paste(src_dir, "/", "dataset_all.txt", sep=""), 
    sep = ",", 
    header = FALSE, # no column name in the dataset 
    col.names = c("ymd", "var1", "var2", "var3", "var4", "var5", + "var6", "var7", "var8", "var9", "var10"), # input column names 
    stringsAsFactor = FALSE, 
    na.strings = "NA") # missing value : "NA"

우측 상단의 environment 창에서 'dataset_all_df' 데이터 프레임이 새로 생겼습니다.

클릭해서 열어보면 아래와 같이 'day_20160701.txt' ~ 'day_20160710.txt'까지 데이터셋이 합쳐져있고, "ymd", "var1" ~ "var10" 까지 칼럼 이름도 생겼습니다.

프로그래밍을 통한 자동화가 중요한 이유, 우리의 시간은 소중하니깐요~! ^^

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

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

(2018.03.14일 내용 추가)

댓글 질문에 '폴더에 있는 개별 파일을 하나씩 읽어와서 하나씩 DataFrame 객체로 메모리상에 생성하는 방법에 대한 질문이 있어서 코드 추가해서 올립니다. 위에서 소개한 방법과 전반부는 동일하구요, 마지막에 루프 돌릴 때 assign() 함수로 파일 이름을 할당하는 부분만 조금 다릅니다.

#========================================================= # read all files in a folder and make a separate dataframe #========================================================= rm(list=ls()) # clear all # (1) directory src_dir <- c("D:/admin/Documents/R/R_Blog/326_read_all_files") # (2) make a file list of all files in the folder src_file <- list.files(src_dir) src_file # (3) count the number of files in the directory => object src_file_cnt <- length(src_file) src_file_cnt # 5 # (4) read files one by one using looping #     => make a dataframe one by one using assign function for (i in 1:src_file_cnt){   assign(paste0("day_", i),           read.table(paste0(src_dir, "/", src_file[i]),                     sep = ",",                     header = FALSE))   print(i) # check progress } rm(src_dir, src_file, src_file_cnt, i) # delete temp objects ls() # list-up all dataframes

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

(2021.08.24일 추가)

댓글에 "여러개의 파일을 하나로 합칠 때 "파일 이름을 데이터 프레임의 새로운 칼럼에 값으로 추가한 후"에 합치는 방법"에 대한 문의가 있었습니다.  댓글란에 코드 블락을 복사해 넣으면 들여쓰기가 무시되어서 보기가 힘들므로 본문에 예제 코드 추가해 놓습니다. 

간단한 샘플 텍스트 파일 3개 만들어서 for loop 순환문으로 각 파일 읽어온 후, 파일 이름을 새로운 칼람 'z'의 값으로 할당 해주고, blank data.frame 인 'day_all' 에 순차적으로 rbind 해주었습니다. 

##--------------------------------------------------------
## add new column with file name and append all dataframes
##--------------------------------------------------------

## blank data.frame to save all files later
day_all <- data.frame()

## file list
src_dir <- c("/Users/lhongdon/Documents/day")
src_file <- list.files(src_dir)
src_file
# [1] "day_20160701" "day_20160702" "day_20160703"


for (i in 1:length(src_file)){
  # read dataset 1 by 1 sequentially
  day_temp <- read.table(
    paste0(src_dir, "/", src_file[i]), 
    sep=",", 
    header=T, 
    stringsAsFactors=F)
  
  # add filename as a new column
  day_temp$z <- src_file[i]
  
  # rbind day_temp to day_all data.frame
  day_all <- rbind(day_all, day_temp)
  
  #print(i) # for progress check
}

print(day_all)
# x y            z
# 1 a e day_20160701
# 2 b f day_20160701
# 3 c g day_20160701
# 4 q w day_20160702
# 5 e r day_20160702
# 6 t y day_20160702
# 7 u i day_20160703
# 8 o p day_20160703
# 9 k l day_20160703

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

(2021.08.25 일 추가)

댓글에 추가 질문이 달려서 요건에 맞게 코드를 더 추가하였습니다.

중첩 for loop 문에 조건절이 여러개 들어가다 보니 코드가 많이 복잡해졌네요. 

[데이터 전처리 요건 ]

1. 로컬 머신 폴더 내 여러개의 csv 파일을 읽어와서 한개의 R data.frame 으로 통합

2. 이때 개별 csv 파일로 부터 읽어들인 데이터를 특정 개수의 [행 * 열] data.frame 으로 표준화

    - 가령, 3 행 (rows) * 3 열 (columns) 의 data.frame 으로 표준화하기 원한다면 

    - 개별 csv 파일로 부터 읽어들인 데이터의 행(row)의 개수가 3보다 크면 1~3행까지만 가져와서 합치고 나머지는 버림. 반대로 3개 행보다 부족하면 'NA' 결측값으로 처리함. 

    - 개별 csv 파일로 부터 읽어들인 데이터의 열(column)이 타켓 칼럼 이름(가령, "x", "y", "z") 중에서 특정 칼럼이 없다면 그 칼럼의 값은 모두 'NA' 결측값으로 처리함.(가령, csv 파일 내에 "x", "y" 만 있고 "z" 칼럼은 없다면 "z" 칼럼을 만들어주고 대신 값은 모두 'NA' 처리해줌)

3. 'day' 라는 칼럼을 새로 만들어서 파일 이름(day 날짜가 들어가 있음)을 값으로 넣어줌

[ 예제 데이터 ]

##--------------------------------------------------------
## (1) 3 rows & 3 cols DataFrame
## (2) add new column with file name and append all dataframes
##--------------------------------------------------------

## blank data.frame to save all files later
day_all <- data.frame()

## file list
src_dir <- c("/Users/lhongdon/Documents/day")
src_file <- list.files(src_dir)
src_file
# [1] "day_20160701" "day_20160702" "day_20160703" "day_20160704"


## setting target rows & cols
row_num <- 3 # set your target number of rows
col_name <- c("x", "y", "z") # set your target name of columns

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

  # read dataset 1 by 1 sequentially
  day_temp <- read.table(
    paste0(src_dir, "/", src_file[i]), 
    sep=",", 
    header=T, 
    stringsAsFactors=F)
  
  ##-- if the number of rows is less than 3 then 'NA', 
  ##-- if the number of rows is greater than 3 than ignore them
  ##-- if the name of columns is not in col_nm then 'NA'
  
  # blank temp dataframe with 3 rows and 3 columns
  tmp_r3_c3 <- data.frame(matrix(rep(NA, row_num*col_num), 
  								nrow=row_num, 
                                byrow=T))
  names(tmp_r3_c3) <- col_name
  
  tmp_row_num <-  nrow(day_temp)
  tmp_col_name <- colnames(day_temp)
  
  r <- ifelse(row_num > tmp_row_num, tmp_row_num, row_num)
  
  for (j in 1:r) {
    for (k in 1:length(tmp_col_name)) {
      tmp_r3_c3[j, tmp_col_name[k]] <- day_temp[j, tmp_col_name[k]]
    }
  }
  
  # add filename as a new column 'day'
  tmp_r3_c3$day <- src_file[i]
  
  # rbind day_temp to day_all data.frame
  day_all <- rbind(day_all, tmp_r3_c3)
  
  rm(tmp_r3_c3)
  print(i) # for progress check
}

print(day_all)
# x    y  z          day
# 1     a    e  1 day_20160701
# 2     b    f  3 day_20160701
# 3     c    g  5 day_20160701
# 4     q    w NA day_20160702
# 5     e    r NA day_20160702
# 6     t    y NA day_20160702
# 7     u    i  3 day_20160703
# 8     o    p  6 day_20160703
# 9  <NA> <NA> NA day_20160703
# 10    e    a  6 day_20160704
# 11    d    z  5 day_20160704
# 12    c    x  3 day_20160704

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

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

군집 간 거리를 측정하는 방법에 따라서 여러가지 알고리즘이 있는데요, 지난번 포스팅에서는 응집형 계층적 군집화(agglomerative hierarchical clustering) 알고리즘 중에서 그래프 기반(Graph-based)의 (1-1) 단일(최단) 연결법 (single linkage method, MIN), (1-2) 완전(최장) 연결법 (complete linkage method, MAX), (1-3) 평균 연결법 (average linkage method) 에 대해 소개하였습니다.

응집형 계층적 군집화 알고리즘 중에서 프로토타입 기반(Prototype-based) 모형은 미리 정해놓은 각 군집의 프로토타입에 터이터가 얼마나 가까운가로 군집의 형태를 결정합니다. 프로토타입 기반 유사성 측도로서 두가지 방법인 (1-4) 중심 연결법 (centroid linkage method)과 (1-5) Ward 연결법 (Ward linkage method) 중에서 이번 포스팅에서는 (1-4) 중심 연결법에 대해서 소개하겠습니다.

중심 연결법(Centroid Linkage method)은 두 군집 간의 거리를 측정할 때 각 군집의 중심(centroid) 간의 거리를 사용합니다.  아래 그림의 왼쪽의 이미지를 참고하시기 바랍니다.

[ 프로토타입 기반 유사성 측정 (Prototype-based maesure of proximity) ]

[표기 (denotation) 설명]

- d(i+j, k) : 군집 i와 j가 합쳐진 새로운 군집 i+j (cluster i U j)와 군집 k 간의 거리(distance)

- μi+j : 군집 i와 r군집 j의 데이터를 가중평균(weighted average)을 이용해 계산한 새로운 중심
- ni : 군집 i의 데이터 개수,  nj : 군집 j의 데이터 개수

- 빨간점 : 각 군집의 중심(centroid)

> ##--------------------------------------------
> ## (1) Agglomerative Hierarchical Clustering 
> ##   (b) Prototype-based
> ##    (1-4) Centroid Linkage
> ##--------------------------------------------
> 
> x <- c(1, 2, 2, 4, 5)
> y <- c(1, 1, 4, 3, 4)
> 
> xy <- data.frame(cbind(x, y))
> 
> xy
  x y
1 1 1
2 2 1
3 2 4
4 4 3
5 5 4
> 
> # scatter plot of xy
> plot(xy, pch = 19, xlab = c("x coordinate"), ylab = c("y coordinate"), 
+      xlim = c(0, 6), ylim = c(0, 6), 
+      main = "scatter plot of xy")
> 
> # adding student label
> text(xy[,1], xy[,2], labels = abbreviate(rownames(xy)), 
+      cex = 0.8, pos = 1, col = "blue") # pos=1 : at the bottom
> 
> 
> # adding dotted line
> abline(v=c(3), col = "gray", lty = 2) # vertical line
> abline(h=c(3), col = "gray", lty = 2) # horizontal line

> # proximity matrix : squares of euclidean distance matrix for 6 points
> dist(xy, method = "euclidean")^2
   1  2  3  4
2  1         
3 10  9      
4 13  8  5   
5 25 18  9  2

> # Agglomerative Hierarchical Clustering : Centroid Linkage method
> hc_cent <- hclust(dist(xy)^2, method="centroid")
> hc_cent

Call:
hclust(d = dist(xy)^2, method = "centroid")

Cluster method   : centroid 
Distance         : euclidean 
Number of objects: 5 

> # dendrogram
> my_par = par(no.readonly = TRUE)
> par(oma = c(0, 0, 1, 0))
> par(mfrow = c(1, 2))
> plot(hc_cent)
> plot(hc_cent, hang = -1) # hang = -1 : line from the bottom

> # custering information
> rev(hc_cent)
$dist.method
[1] "euclidean"

$call
hclust(d = dist(xy)^2, method = "centroid")

$method
[1] "centroid"

$labels
NULL

$order
[1] 1 2 3 4 5

$height
[1]  1.00000  2.00000  6.50000 11.80556

$merge
     [,1] [,2]
[1,]   -1   -2
[2,]   -4   -5
[3,]   -3    2
[4,]    1    3

> # comparison of height among linkage methods
> hc_sl <- hclust(dist(xy)^2, method="single")
> hc_cl <- hclust(dist(xy)^2, method="complete")
> hc_avg <- hclust(dist(xy)^2, method="average")
> hc_cent <- hclust(dist(xy)^2, method="centroid")
> 
> # dendrogram
> my_par = par(no.readonly = TRUE)
> par(oma = c(0, 0, 1, 0))
> par(mfrow = c(1, 4))
> plot(hc_sl, main = "Single Linkage")
> plot(hc_cl, main = "Complete Linkage")
> plot(hc_avg, main = "Average Linkage")
> plot(hc_cent, main = "Centroid Linkage")