1-1.결측치 제거

• We have 6 variables, and there seems to be 2 missing values(NA).
• Amelia 패키지의 missingness map을 이용하여 결측치를 시각화해보니, 2개의 관측값이 누락된 것 같다.

• 누락된 결측값은 DH0722와 JL0766이다.
• 샘플수가 충분히 많으므로, 두 개의 결측값은 제외하더라도 분석에 큰 영향을 미치지 않을 것으로 보인다. 따라서 두 결측값을 제외하고 분석을 진행하기로 한다.

    RowNames p4075 p4076 p4077 p4078 p4079 p4080
286   DH0722    NA    NA    NA    NA    NA    NA
357   JL0776    NA    NA    NA    NA    NA    NA

> cnp6.narm <- na.omit(cnp6)

1.2 데이터의 모양(분포파악)

• scatterplot matrix을 이용하여 각 변수의 선형적 상관관계를 살펴보도록 하자.
• 변수 p4079와 p4075가 0.59라는 제일 큰 상관관계를 보인다.
• 변수 p4078 ~ p2080의 상관계수도 모두 0.5보다 큰 것으로 보아 상관관계가 있어 보인다.
• 반면 p4076 ~ p4078은 상관관계가 거의 없고 무작위로 흩어져있는 느낌이 강하다.
• 2차원 상의 산포도를 보는 것만으로는 군집의 개형이 보이지 않는다.

2-1. 다양한 연결법 비교

• 어떤 방법을 사용해야 가장 좋은 결과를 얻을 수 있을까?
  
  

2-2.최적의 연결법 선정하기

• 위 8개의 덴드로그램을 종합하면 3개 군집으로 나누는 것이 가장 좋아보인다. 그렇다면 어떤 방법을 사용 했을 때 3개 군집이 가장 적절하게 만들어질 수 있을까?

• canberra거리행렬을 이용한 연결법은 각 군집에 균일한 할당이 안 된다.

• 반면, maximum거리행렬을 이용한 연결법은 전체적으로 군집을 더 균일하게 나누는 양상을 보인다.

• 그 중에서도 maximum거리행렬을 이용한 평균연결법(제일 왼쪽)이 군집화가 가장 적절하게 되는 것으로 보인다. (관측값이 각 군집별로 균일하게 할당됨)

2-3 결과 시각화하기

아래 방법을 따라 군집화한 결과가 가장 좋았다. - number of clusters : 3 - distance method : maximum - clustering method : average

maximum거리행렬과 average linkage method를 이용하여 3개 군집으로 나누었을 때,
군집들 사이의 거리도 충분히 멀고, 각 군집 별로 관측수도 1/3 내외로 균일하게 잘 분배된다. 따라서 제일 좋은 방법이라 생한다.

3-1. 군집개수 정하기 (with Elbow Method)

• 가장 간단한 방법을 사용하자.
• 숫자 3 이후로 기울기가 급격하게 완만해지므로, 3개의 군집으로 나누는 것이 좋겠다.

3-2. 결과 시각화하기(with scatterplot matrices)

• 군집의 개수는 3개로 놓고 군집분석을 수행해보자.
  
  

• 군집이 순서대로 오른쪽에서 왼쪽으로 이어지는 양상을 볼 수 있다.

• 검은색 군집 대체적으로 낮은 값을 가지고(왼쪽 아래), 빨간색 군집이 중간 값을 가지고, 초록색 군집이 대체적으로 큰 값(오른쪽 위)을 가지는 양상을 볼 수 있다.

3-3. 결과 시각화하기(with PCA score plot)

• scatterplot matrix 상에서는 군집이 구분되는게 잘 보이지 않아 주성분분석을 하여 보기로 한다.
• 첫번째와 두번째 주성분을 축으로 두어 그래프 위에 표시하였다.
• 그랬더니 scatterplot matrix서는 볼 수 없었던 깔끔한 결과를 볼 수 있게 되었다.
• 각 군집들이 겹치지 않게 배타적인 집단으로 나뉜 것을 볼 수 있다.
• 3번째 주성분까지 포함하면 80%의 데이터를 파악할 있는데 그러기 위해서는 3차원을 이용해야하므로 복잡하다. 따라서 생략하도록 한다.
• 아래는 pca분석 결과이다.
• 2번째 주성분까지 포함하면 전체 데이터의 66.4%를 설명할 수 있으며, 3번째 주성분까지 합하면 80%의 데이터를 파악할 수 있다.

Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6
Standard deviation     1.6687 1.0986 0.8956 0.69001 0.63377 0.57334
Proportion of Variance 0.4641 0.2011 0.1337 0.07935 0.06694 0.05479
Cumulative Proportion  0.4641 0.6652 0.7989 0.87827 0.94521 1.00000

4-1. 결과의 시각화

• mixture model을 이용해 군집 분석한 결과는 아래와 같다.

4-2. 가정 확인하기(정규성)

• mixture model에 의한 군집들이 다변량정규분포를 따르는지 확인해보자.
• 각 집단이 다변량정규분포를 따라야만 mixture model의 결과를 신뢰할 수 있기 때문이다.
• 만약 각 집단(군집)이 다변량정규분포를 따르지 않는다면 mixture model을 적용하고 해석하는데 어려움이 있다.

sROC 0.1-2 loaded

   Mardia's Multivariate Normality Test 
--------------------------------------- 
   data : cluster1[, 2:7] 

   g1p            : 4.104177 
   chi.skew       : 58.14251 
   p.value.skew   : 0.396363 

   g2p            : 46.44312 
   z.kurtosis     : -0.7324851 
   p.value.kurt   : 0.4638725 

   chi.small.skew : 60.80059 
   p.value.small  : 0.3071751 

   Result         : Data are multivariate normal. 
--------------------------------------- 

• 첫번째 군집은 정규성을 따르나(위 결과), 두번째 군집(아래 결과)은 정규성가정을 만족하지 못한다. 따라서 Mixture model의 결과를 받아들이기 어렵다.

   Mardia's Multivariate Normality Test 
--------------------------------------- 
   data : cluster2[, 2:7] 

   g1p            : 1.229773 
   chi.skew       : 65.17797 
   p.value.skew   : 0.1877657 

   g2p            : 45.80559 
   z.kurtosis     : -1.996947 
   p.value.kurt   : 0.0458309 

   chi.small.skew : 65.97012 
   p.value.small  : 0.1701646 

   Result          : Data are not multivariate normal. 
--------------------------------------- 

분포에 대한 가정이 없는 Kmeans나 Hierarchical clustering을 이용하는 것이 더 정확한 군집분석 결과를 낼 것이다.

5.1 군집분석 방법들의 비교

5.2 최종 결과

6-1. 데이터 준비

#데이터 살펴보기
cnp6 <- read.csv("./cnp6.csv")
head(cnp6)
str(cnp6)
summary(cnp6)

# 결측값 탐색하기
library(Amelia)
library(RColorBrewer)
missmap(cnp6[2:7], legend=TRUE, col=brewer.pal(2, "Pastel1" ))

# 결측값 제거하기
cnp6[!complete.cases(cnp6),]
cnp6.narm <- na.omit(cnp6)

# 산점도 행렬 그리기
library(ggplot2)
panel.lm <- function(x, y, col=par("col"), bg=NA, pch=par("pch"),
                     cex=1, col.smooth="black", ...) {
  points(x, y, pch = 20, col = alpha("midnightblue", 0.2), bg = bg, cex = cex)
  abline(stats::lm(y~x), col = col.smooth, ...)
}

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
  usr <- par("usr")
  on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if (missing(cex.cor)) cex.cor <- 0.8 / strwidth(txt)
  text(.5, 0.5, txt, cex = cex.cor * r)
}


panel.hist <- function(x, ...) {
  usr <- par("usr")
  on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5))
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks
  nB <- length(breaks)
  y <- h$counts
  y <- y / max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "darkblue", ...)
}


pairs(
  cnp6.narm[,2:7],
  lower.panel = panel.lm,
  upper.panel = panel.cor,
  diag.panel = panel.hist,
  pch = 20,
  main = "Scatterplot/Histogram/Correlation coeff. of cnp6 data"
)

6-2. Hierarchical clustering

# 표준화하기
cnp.scaled <- cbind(cnp6.narm$RowNames, scale(cnp6.narm[, 2:7]))

library(NbClust)
library(factoextra)

# 거리행렬 구하기
cnp.dist.canberra <- dist(cnp.scaled, method = "canberra")
cnp.dist.maximum <- dist(cnp.scaled, method = "maximum")
cnp.dist.mink <- dist(cnp.scaled, method = "minkowski")

# 다양한 연결법 시도해보기
cnp.hclust.man <- hclust(cnp.dist.man, method = "average")
cnp.hclust.euc <- hclust(cnp.dist.euc, method = "average")
cnp.hclust.can <- hclust(cnp.dist.can, method = "average")
cnp.hclust.max1 <- hclust(cnp.dist.max, method = "average") #less good
cnp.hclust.can.max <- hclust(cnp.dist.canberra, method = "complete")
cnp.hclust.can.war <- hclust(cnp.dist.canberra, method = "ward")
cnp.hclust.can.sin <- hclust(cnp.dist.canberra, method = "single")
cnp.hclust.can.ave <- hclust(cnp.dist.canberra, method = "average")
cnp.hclust.max.max <- hclust(cnp.dist.maximum, method = "complete")
cnp.hclust.max.war <- hclust(cnp.dist.maximum, method = "ward")
cnp.hclust.max.sin <- hclust(cnp.dist.maximum, method = "single")
cnp.hclust.max.ave <- hclust(cnp.dist.maximum, method = "average") # BEST RESULT

# 덴드로그램 그려서 살펴보기
plot(cnp.hclust.euc)
plot(cnp.hclust.man)
plot(cnp.hclust.can)
plot(cnp.hclust.euc)
plot(cnp.hclust.max1)
plot(cnp.hclust.can.max, cex=0.1, hang=0.1)
plot(cnp.hclust.can.war, cex=0.1, hang=0.1)
plot(cnp.hclust.can.sin, hang=.1, cex=0.01)
plot(cnp.hclust.can.ave, hang=.01 ,cex=0.01)
plot(cnp.hclust.max.max)
plot(cnp.hclust.max.war)
plot(cnp.hclust.max.sin, cex=.1, hang=0.1)
plot(cnp.hclust.max.ave, cex=.1, hang=0.1)

# 덴드로그램 그리기
library(RColorBrewer)
fviz_dend(
  cnp.hclust.max.ave, k = 3, cex = 0.3,
  k_colors = brewer.pal(3, "Set1"),
  color_labels_by_k = TRUE, rect = TRUE
)

6-3. Kmeans clustering

# 군집 개수 정하기
fviz_nbclust(cnp.scaled, kmeans, method = "wss") +
  geom_vline(xintercept = 3, linetype = 2) +
  labs(subtitle = "Elbow method")

# 군집분석 수행하기
m <- 20
set.seed(100)

cnp.km.3 <- kmeans(
  cnp.scaled[, 2:7],
  centers = 3, nstart = 40
)

#산접도 행렬을 이용해 결과 확인해보기
panel.lm.km <- function(x, y, col=par("col"), bg=NA, pch=par("pch"),
                        cex=0.2, col.smooth="black", ...) {
  points(x, y, pch = cnp.km.3$cluster, col = cnp.km.3$cluster, bg = bg, cex = cex)
  abline(stats::lm(y~x), col = col.smooth, ...)
}

pairs(
  cnp.scaled[, 2:7],
  lower.panel = panel.lm.km,
  upper.panel = panel.cor,
  diag.panel = panel.hist,
  pch = 20,
  main = "scatter plot with Kmeans clustering"
)

# score plot 상에서 결과 확인해보기
fviz_cluster(cnp.km.3, data = cnp.scaled[, 2:7])

# pca 결과 
pr.out <- prcomp(cnp6.narm[2:7], scale=TRUE)
summary(pr.out)

6-4. Mixture Model

# mixture model 돌리기
library(mclust)
fit <- Mclust(cnp6.narm[, 2:7])
summary(fit)기

# 시각화
plot(fit, what="classification")

# 정규성 검정 
cluster1 <- cnp.scaled[fit$classification==1,]
cluster2 <- cnp.scaled[fit$classification==2,]

library(MVN)
mardiaTest(cluster1[,2:7])
mardiaTest(cluster2[,2:7])#정규성 위반