#======================================================
#	9 クラスター分析、主成分分析及び因子分析
#======================================================
#======================================================
#	9.1 クラスター分析
#======================================================
#パッケージのロード
library(cluster)

#架空のデータ例

set.seed(456)

v1 <- rnorm(10)                               #連続変数
v2 <- runif(10)                               #連続変数
v3 <- sample(c("A", "B", "C"), 10, replace=T) #名義尺度（因子）
v4 <- ordered(sample(c("不満", "普通", "満足"), 10, replace=T), 
	levels=c("不満", "普通", "満足"))      #順序因子
v5 <- sample(0:1, 10, replace=T)              #二値変数（対称）
v6 <- sample(0:1, 10, replace=T)              #二値変数（非対称）
X <- data.frame(v1, v2, v3, v4, v5, v6)

#データの表示

X

#非類似度行列

D <- daisy(X, type=list(ordratio=4, symm=5, asymm=6))
print(D, digits=3)

#---------------<階層型手法>---------------
#---------------<都道府県別データ集計>---------------

obj.DF <- data.frame(
    log人口密度 = log10(都道府県$人口密度),
    高齢者人口率 = 都道府県$人口65歳以上/都道府県$人口総数,
    単独世帯率 = 都道府県$単独世帯/都道府県$一般世帯,
    高齢単身者率 = 都道府県$高齢単身者/都道府県$一般世帯,
    ３世代世帯率 = 都道府県$３世代世帯/都道府県$一般世帯,
    就業率 = 都道府県$就業者数/(都道府県$"人口15〜64歳" + 都道府県$人口65歳以上),
    第３次産業率 = 都道府県$第３次産業/都道府県$就業者数,
    農業と林業率 = 都道府県$農業と林業/都道府県$就業者数,
    製造業率 = 都道府県$製造業/都道府県$就業者数,
    医療と福祉率 = 都道府県$医療と福祉/都道府県$就業者数
)
row.names(obj.DF) <- コード$都道府県

#---------------<都道府県別データ：散布図>---------------
#図9.1
pairs(obj.DF)

#---------------<凝集型階層：agnes>---------------

ken.agnes <- agnes(obj.DF, stand=T)                        #当てはめ
#図9.2
windows(10, 7)
plot(ken.agnes, which=2, main="", xlab="", ylab="非類似度")#プロット

#---------------<区分型階層：diana>---------------

ken.diana <- diana(obj.DF, stand=T)         #当てはめ
#図9.3
windows(10, 7)
plot(ken.diana, which=2, main="", xlab="", ylab="非類似度")#プロット

#---------------<分割型手法>---------------
#---------------<メドイド近傍分割：pam>---------------

# k=2
ken.pam2 <- pam(obj.DF, k=2, stand=T)
plot(ken.pam2, which=1, main="k=2", labels=2, col.txt="blue", cex.txt=0.7)
plot(ken.pam2, which=2, main="k=2")
#ave(s(i))=0.18

# k=3：おそらく、これが最良
ken.pam3 <- pam(obj.DF, k=3, stand=T)
#図9.4
plot(ken.pam3, which=1, main="k=3", labels=2, col.txt="blue", cex.txt=0.7)
#図9.5
plot(ken.pam3, which=2, main="k=3")
#ave(s(i))=0.33
summary(ken.pam3)  #要約

# k=4
ken.pam4 <- pam(obj.DF, k=4, stand=T)
plot(ken.pam4, which=1, main="k=4", labels=2, col.txt="blue", cex.txt=0.7)
plot(ken.pam4, which=2, main="k=4")
#ave(s(i))=0.27

# k=5
ken.pam5 <- pam(obj.DF, k=5, stand=T)
plot(ken.pam5, which=1, main="k=5", labels=2, col.txt="blue", cex.txt=0.7)
plot(ken.pam5, which=2, main="k=5")
#ave(s(i))=0.25

# k=6
ken.pam6 <- pam(obj.DF, k=6, stand=T)
plot(ken.pam6, which=1, main="k=6", labels=2, col.txt="blue", cex.txt=0.7)
plot(ken.pam6, which=2, main="k=6")
#ave(s(i))=0.25

#======================================================
#	9.2 主成分分析
#======================================================

#パッケージのロード
library(rrcov)  #主成分分析のため  
library(MASS)   #関数mvrnormのため

#---------------<分散共分散行列と相関行列>---------------

#分散共分散行列
cov.ken <- cov(obj.DF[ ,-c(1, 4, 7)])
round(cov.ken, 2)

#相関行列
cor.ken <- cor(obj.DF[ ,-c(1, 4, 7)])
round(cor.ken, 2)
 
#---------------<固有値分解による主成分分析>---------------
#---------------<主成分分析の概念図>---------------

set.seed(456)

#多変量正規分布に従う乱数

X <- mvrnorm(50, mu=c(2.5, 2), Sigma=matrix(c(1.1, 0.85, 0.85, 1.5), nrow=2))
dimnames(X)[[2]] <- c("x1", "x2")              #変数名を付ける。
n <- nrow(X)                                   #観測値の数

x.eig <- eigen(cov(X))　　　　                 #固有値分解
pc <- X %*% x.eig$vectors                      #固有ベクトルによる回転
pc.x1 <- lm((pc[ ,1]-mean(pc[ ,1])) ~ pc[ ,1]) 
pred.x1 <- predict(pc.x1) + mean(pc[ ,1])      #縮約された座標軸上の値

#図9.6
plot(X, xlim=c(-1, 8), ylim=c(-3, 6), pch=4, cex=0.8) #散布図
abline(h=0, v=0, col="green4", lty=2, lwd=2)          #元々の座標軸  
points(pc, col="tomato", pch=16)                      #回転後のデータ
legend("topright", legend=c("原データ", "回転後のデータ"), 
	pch=c(4, 16), col=c("black", "tomato"))

#図9.7
plot(X, xlim=c(-1, 8), ylim=c(-3, 6), pch=4, cex=0.8,
        xlab="PC1", ylab="") #散布図
abline(h=0, v=0, col="green4", lty=2, lwd=2) #元々の座標軸
points(pc, col="tomato", pch=16)             #回転後のデータ
abline(h=mean(pc[ ,2]),lwd=3)                #縮約された座標軸
points(pred.x1, rep(mean(pc[ ,2]), n), 
        pch=16, col="steelblue")             #縮約された座標軸上の点
segments(x0=pc[ ,1], y0=pc[ ,2], x1=pred.x1, y1=mean(pc[ ,2]))
legend("topright", legend=c("縮約された座標軸", "座標軸上の点"), 
	lty=c(1, NA), pch=c(NA, 16), col=c("black", "steelblue"))

#---------------<適用例>---------------

#主成分分析

ken.pca <- PcaClassic(obj.DF, scale=T)
summary(ken.pca)　#要約

#固有値の確認

ken.eig <- eigen(cor(obj.DF)) #固有値分解
Lambda <- ken.eig$value       #固有値
sqrt(Lambda)                  #標準偏差
Lambda/sum(Lambda)            #分散
#図9.8
screeplot(ken.pca, main="")   #スクリープロット
cumsum(Lambda)/10             #累積分散

#主成分負荷量の確認
ken.pca@loadings            #主成分負荷量
ken.eig$vectors             #固有ベクトル

#主成分負荷量の図示
plot.loadings <- function(pca, k=2) {
    load <- getLoadings(pca)
    op <- par(mfrow=c(k, 1), mar=c(3, 4, 2, 2))
    for (i in 1:k) {
        PC <- load[ ,i]
        ord <- rev(order(abs(PC)))
        barplot(PC[ord], ylab=paste("PC", i, sep=""))
    }
    par(op)
}

#図9.9
windows(10, 7)
plot.loadings(ken.pca, k=4) #４つの主成分まで表示

#主成分得点の確認
head(ken.pca@scores)                   #主成分得点
X <- scale(obj.DF)                     #原データの標準化
head(as.matrix(X) %*% ken.eig$vectors) #確認

#バイプロット
#図9.10
biplot(ken.pca, cex=0.8)

#---------------<特異値分解による主成分分析>---------------

ken.pca2 <- prcomp(obj.DF, scale=T)
summary(ken.pca2)         #要約

#確認
ken.svd <- svd(X)         #特異値分解
ken.svd$d/sqrt(nrow(X)-1) #標準偏差
ken.svd$v                 #主成分負荷量

#---------------<ロバストな主成分分析>---------------

ken.pca3 <- PcaHubert(X, k=4) #ロバストな主成分分析
#図9.11
biplot(ken.pca3, cex=0.8)     #バイプロット

#======================================================
#	9.3 因子分析
#======================================================

ken.fa <- factanal(X, factors=4, scores="r")      #因子分析
ken.fa
ken.fa2 <- factanal(X, factors=3, scores="r")      #因子分析
ken.fa2

plot.loadings2 <- function(fa, k=2) {
    load <- loadings(fa)
    op <- par(mfrow=c(k, 1), mar=c(3, 4, 2, 2))
    for (i in 1:k) {
        Factor <- load[ ,i]
        ord <- rev(order(abs(Factor)))
        barplot(Factor[ord], ylab=paste("Factor", i, sep=""))
    }
    par(op)
}

#図9.12
windows(10, 7)
plot.loadings2(ken.fa2, k=3) #因子負荷量

#バイプロット

#第1因子と第2因子
#図9.13
biplot(ken.fa2$scores[ ,1:2], ken.fa2$loadings[ ,1:2], cex=0.8)

#第2因子と第3因子
#図9.14
biplot(ken.fa2$scores[ ,2:3], ken.fa2$loadings[ ,2:3], cex=0.8)

