practice problem 2

# 사전작업
str(olympic)
str(olympic$score)

str(olympic$tab)
names(olympic$tab)

# 정렬하기
olympic$tab = olympic$tab[,c(1,5,6,10,4,2,8,3,7,9)]
names(olympic$tab)

#### 1번
# na값 파악하기
sum(is.na(olympic$tab))
sum(is.na(olympic$score))

# 산점도 파악
windows()
par(mfrow=c(2,5))
for(i in 1:10){
    with(olympic, plot(score~tab[,i],
         xlab="score", ylab=names(tab)[i]))
}
# 결과: 계주 종목 4개(100, 500, 110, 1500)은 점수와 음의 상관관계가 있음. 더 빨리 뛰어 시간이 줄 수록 점수가 높아짐.
# 또 점프 이벤트 3개는 점수와 양의 상관관계가 있음. 더 높이 점프할수록 점수가 높아짐.
# 던지기 이벤트 3개는 점수와 양의 상관관계가 있음. 더 높게 던질수록 점수가 높아짐.

# 상관관계 파악 1 : 
with(olympic, round(cor(tab, score), 2))
# 결과: 앞서 본 바와 같이 계주는 점수와 음의 상관관계. 점프와 던지기는 점수와 양의 상관관계를 가지고 있음.

# 상관관계 파악 2:
with(olympic, round(cor(tab[, 1:4]), 2))
# 단거리(100, 110, 400)은 서로 상관계수가 높은 반면 장거리(1500)은 단거리와 상관계수가 낮다. 다만 1500과 400은 연관성이 높은 편이다.

# 상관관계 파악 3:
with(olympic, round(cor(tab[, 5:7]), 2))
# 점프 종목끼리는 서로 상관계수가 낮다.

# 상관관계 파악 4:
with(olympic, round(cor(tab[, 8:10]), 2))
# 던지기 종목끼리는 상관계수가 0.44~0.81까지 서로 높은 편이다.   

# boxplot 그리기
# 먼저 연속형인 score을 이산형으로 바꾼다.
summary(olympic$score)
score_cut = cut(olympic$score, breaks=c(6907, 7579, 7869, 8180, 8488), labels=c("G1", "G2", "G3", "G4"), include.lowest = TRUE)

windows()
par(mfrow=c(2,5))
for(i in 1:10){
    with(olympic, boxplot(tab[,i]~score_cut,
                    ylab=names(tab)[i]))
}
# 역시 비슷한 결과를 보여준다.

### 2번
# score_cut 기반 산점도 보기
windows()
par(mfrow=c(2,5))
for(i in 1:10){
    with(olympic, plot(x=score, y=tab[,i],
         xlab="score", ylab=names(tab)[i],
         col=score_cut, pch=19))
}

# 그룹별 상관관계 찾아보기
cor_score_cut = matrix(nrow=10, ncol=4)
for(i in 1:4){
    cor_score_cut[,i] = with(olympic, cor(tab[score_cut == levels(score_cut)[i], ], score[score_cut == levels(score_cut)[i]]))
}
rownames(cor_score_cut) = names(olympic$tab)
colnames(cor_score_cut) = levels(score_cut)

cor_score_cut


# 평균과 분산 찾아보기
mean_score_cut = matrix(nrow=10, ncol=4)
sd_score_cut = matrix(nrow=10, ncol=4)
for(i in 1:10){
    mean_score_cut[,i] = tapply(olympic$tab[,i], score_cut, mean)
    sd_score_cut[,i] = tapply(olympic$tab[,i], score_cut, sd)
}

rownames(mean_score_cut) = levels(score_cut)
rownames(sd_score_cut) = levels(score_cut)

colnames(mean_score_cut) = names(olympic$tab)
colnames(sd_score_cut) = names(olympic$tab)

### 3-a
# 공분산 pca
# 상관계수 pca
pca1 = prcomp(olympic$tab) 
pca2 = prcomp(olympic$tab, scale=T) 

# biplot으로 좋은 기준 발견하기
windows()
par(mfrow=c(1,1))
biplot(pca2)
abline(v=0, col="blue", lty=2)
abline(h=0, col='blue', lty=2)

# biplot 대신 할 수 있는 것
pcs = scale(olympic$tab) %*% pca2$rotation
plot(pcs[,1], pcs[,2], type="n")
text(pcs[,1], pcs[,2], labels=rownames(pcs), cex=0.8)
abline(v=0, col="blue", lty=2)
abline(h=0, col='blue', lty=2)

summary(pca2)

# pca
# biplot 이벤트를 볼 때, pc1의 관점에서 달리기 종목과 다른 종목는 서로 반대 관계임.
# 속도가 높아야 하는 종목들일 수록 pc1의 절댓값이 높아짐. 반대로 pc2는 힘이 좋아야 하는 종목일수록 절댓값이 높아짐.

# 따라서 pc1은 속도이며, pc2는 근육의 힘이다

### 3-b
par(mfrow=c(1,2))
screeplot(pca2)
screeplot(pca2, type='lines')

# 분산과 비율을 같이 계산하고 그리는 방법
par(mfrow=c(1,2))
plot(pca2)

tmp = summary(pca2)
prop = tmp$importance[2,]
plot(prop, main='scree plot', type='b',
     ylab='proportion')

# 모으기
cumsum(pca2$sdev^2 / sum(pca2$scale^2))
tmp$importance[3,]

# 60% 정도의 총 분산이 2개의 주성분으로 설명될 수 있음.

## 3-C. 두 주성분과 최종 score간의 관계는?
par(mfrow=c(1,2))
plot(pca2$x[,1], olympic$score, xlab="pc1", ylab='score')
plot(pca2$x[,2], olympic$score, xlab="pc2", ylab='score')

# pc1이 더 선형적임. 따라서 pc1이 더 중요.

### 4
fa1 = factanal(olympic$tab, factors=2, rotation='none')
fa2 = factanal(olympic$tab, factors=2, rotation='varimax', scores="reg")
fa3 = factanal(olympic$tab, factors=2, rotation='promax')

windows()
par(mfrow=c(1,3))
# 회전 없는 경우
with(fa1, plot(loadings[,1], loadings[,2],
    xlim=c(-1,1), ylim=c(-1,1), 
    main="회전 없음", 
    xlab="요인1", ylab="요인2"))
with(fa1, text(loadings[,1], loadings[,2],
        labels=rownames(loadings), pos=4))
abline(v=0, col='red')
abline(v=0.5, col='blue', lty=2)
abline(v=-0.5, col='blue', lty=2)

abline(h=0, col='red')
abline(h=0.5, col='blue', lty=2)
abline(h=-0.5, col='blue', lty=2)

# varimax 들어간 경우
with(fa2, plot(loadings[,1], loadings[,2],
        xlim=c(-1,1), ylim=c(-1,1), 
        main="varimax", 
        xlab="요인1", ylab="요인2"))
with(fa2, text(loadings[,1], loadings[,2],
        labels=rownames(loadings), pos=4))
abline(v=0, col='red')
abline(v=0.5, col='blue', lty=2)
abline(v=-0.5, col='blue', lty=2)

abline(h=0, col='red')
abline(h=0.5, col='blue', lty=2)
abline(h=-0.5, col='blue', lty=2)

# promax 들어간 경우
with(fa3, plot(loadings[,1], loadings[,2],
        xlim=c(-1,1), ylim=c(-1,1), 
        main="varimax", 
        xlab="요인1", ylab="요인2"))
with(fa3, text(loadings[,1], loadings[,2],
        labels=rownames(loadings), pos=4))
abline(v=0, col='red')
abline(v=0.5, col='blue', lty=2)
abline(v=-0.5, col='blue', lty=2)

abline(h=0, col='red')
abline(h=0.5, col='blue', lty=2)
abline(h=-0.5, col='blue', lty=2)

# 요인1 - 달리기 & 던지기 이벤트와 상관있음.
# 요인2 - 던지기 이벤트랑 상관있음.

fa1$loadings
fa2$loadings
fa3$loadings
# fa1 회전없을 때 : 총분산의 52.3% 설명
# fa2 varimax : 총분산의 52.3% 설명
# 뭐가 더 나은지 찾기 좀 힘들어도, fa1에서 가장 큰 분산 설명력을 가진 varimax가 제일 적합하다. 


## pca와 fa의 비슷한 점 찾기
windows()
par(mfrow=c(1,2))
plot(fa2$scores, main="FA results with 2 factors")
abline(h=0)
abline(v=0)
text(fa2$scores, labels=1:33, pos=4)

biplot(pca2, main="pca")
abline(h=0)
abline(v=0)
# 굉장히 비슷하게 생김.

# final score과 fa 성적
par(mfrow=c(1,2))
plot(fa2$scores[, 1], olympic$score, xlab='요인1', ylab='성적')
plot(fa2$scores[, 2], olympic$score, xlab='요인2', ylab='성적')
# fa1은 final score과 음의 상관.
# fa2는 final score과 양의 상관.
# fa1이 더 강한 관계 지님.


cor(cbind(fa2$scores[, 1:2], olympic$score))
# 성적과 fa1의 상관관계는 양수
# 성적과 fa2의 상관관계는 음수
# fa1이 더 강한 관계(= 더 큰 절댓값)

### 5번
tab_scale = scale(olympic$tab)
tab_scale

a1 = apply(tab_scale[, 1:4], 1, mean)
a2 = apply(tab_scale[, 5:7], 1, mean)
a3 = apply(tab_scale[, 8:10], 1, mean)
mean_tab = data.frame(a1, a2, a3)
names(mean_tab) = c('score.run', 'score.jump', 'score.throw')


# 마지막 점수와 비교
par(mfrow=c(1, 3))
plot(mean_tab[,1], olympic$score)
plot(mean_tab[,2], olympic$score)
plot(mean_tab[,3], olympic$score)

cor(mean_tab, olympic$score)

# k-means 클러스터링: k = 2~4
km2 = kmeans(mean_tab, 2, nstart=5)
km3 = kmeans(mean_tab, 3, nstart=5)
km4 = kmeans(mean_tab, 4, nstart=5)

# 실루엣 그리기
library(cluster)
mean_tab_d = dist(mean_tab)
km2_sil = silhouette(km2$cluster, mean_tab_d)
km3_sil = silhouette(km3$cluster, mean_tab_d)
km4_sil = silhouette(km4$cluster, mean_tab_d)

par(mfrow=c(3,1))
plot(km2_sil, main='실루엣 k = 2')
plot(km3_sil, main='실루엣 k = 3')
plot(km4_sil, main='실루엣 k = 4')
# 실루엣 스코어 모두 0.4 이하
# 3개의 평균 중 2개만 사용할 것.

mean_tab2 = mean_tab[,-2]
library(cluster)
km22 = kmeans(mean_tab2, 2, nstart=5)
km32 = kmeans(mean_tab2, 3, nstart=5)
km42 = kmeans(mean_tab2, 4, nstart=5)

mean_tab_d2 = dist(mean_tab2)
km2_sil2 = silhouette(km22$cluster, mean_tab_d2)
km3_sil2 = silhouette(km32$cluster, mean_tab_d2)
km4_sil2 = silhouette(km42$cluster, mean_tab_d2)

windows()
par(mfrow=c(1, 1))
plot(mean_tab2[, 1],mean_tab2[, 2],type="n")
points(km42$centers, pch=19, col='red')
cols = c('black', 'blue', 'purple', 'red')
text(mean_tab2[, 1], mean_tab2[, 2],
     rownames(mean_tab2),
     col=cols[km42$cluster])
abline(h=0, col='blue', lty=3)
abline(v=0, col='blue', lty=3)

###### 6
# 종목별 기록으로 정리하기기
trans_tab = t(olympic$tab)
id = rownames(olympic$tab)
game_name = colnames(olympic$tab)

data_set2 = data.frame(
    id=rep(id, each = length(game_name)), 
    game_name = rep(game_name, times=length(id)), 
    record=as.vector(trans_tab)
)

head(data_set2)