# 教師なし学習  ```{r, setup, include=FALSE, echo=FALSE} knitr::opts_chunk$set( collapse = TRUE ) pacman::p_load(tidyverse) ``` ## クラスタリング ## 階層ありクラスタリング ### 距離行列 [ Github Gistのページ ](https://gist.github.com/ctsaran/42728dad3c7d8bd91f1d) からダウンロードしています。`dist` 関数というものを用いて、緯度・経度のデータからユークリッド距離の距離行列を求めることができます。```{r, filename="データの準備"} kansai <- read.csv("./data/pref_lat_lon.csv", header = T, fileEncoding = "CP932") # 関西の都道府県データのみ取得する kansai <- kansai[24:30, ] kansai # 各県の緯度と経度 ggplot( kansai, aes(x = lon, y = lat, color = pref_name, label = pref_name)) + geom_text() # 県庁所在地の位置 ``` ```{r, filename="データを距離行列に変換"} # ユークリッド距離（数値は24から順に三重、滋賀、京都、大阪、兵庫、奈良、和歌山） dist(kansai[, 2:3], method="euclidean") ``` ## ユークリッド距離 ```{r} # ユークリッド距離 #（軸の数値は24から順に三重、滋賀、京都、大阪、兵庫、奈良、和歌山） dist (kansai[, 2 : 3 ], method = "euclidean" ) ``` ## 最大距離 ```{r} # 縦と横の距離の最大値を取る dist (kansai[, 2 : 3 ], method = "maximum" ) ``` ## マンハッタン距離 ```{r} # マンハッタン距離 dist (kansai[, 2 : 3 ], method = "manhattan" ) ``` ## キャンベラ距離 ```{r} # キャンベラ距離 #（原点近くの距離を大きく見積もる、2点の距離を原点からの距離の絶対値の和で割る） dist (kansai[, 2 : 3 ], method = "canberra" )``` ## バイナリ距離 ```{r} # バイナリ距離（距離が0なら1、0でなければ0を返す） dist (kansai[, 2 : 3 ], method = "binary" ) ``` ## ミンコフスキー距離 ```{r} # ミンコフスキー距離 #（ユークリッドとマンハッタンの中間的なもの、pという引数を別途取る。pのデフォルトは2） dist (kansai[, 2 : 3 ], method = "minkowski" ) ``` ### Rでの階層ありクラスタリング `hclust` 関数を用いて行います。`hclust` 関数の第一引数は距離行列で、上に述べた通り`dist` 関数を用いて計算できるものです。`method` 引数には、クラスタリングの方法を選択します。クラスタリングの方法は、`"complete"` （デフォルト）、`"ward.D"` 、`"ward.D2"` 、`"single"` 、`"average"` 、`"mcquitty"` 、`"median"` 、`"centroid"` の8つから選択します。`hclust` 関数の返り値を`plot` することで、クラスタリングの結果を確認することができます。`hclust` 関数の返り値自体を確認することはあまりないため、基本的には`plot` して使うものだと考えてもらうとよいと思います。`plot` 関数において、`hang` 引数を`-1` に指定すると、クラスタリングする要素の下端をそろえることができます。`plot` 関数で表示されるのは、より近いもの同士も線でつないだ樹形上の図です。線でつながった距離が近いほど、データがよく似ていることを示しています。上記の都道府県の緯度・経度でクラスタリングを行うと、距離の近さとクラスタリングの図上の近さが対応しているのがわかるかと思います。## 下揃えなし ```{r, filename="階層ありクラスタリング", fig.width=10} JPll <- read.csv("./data/pref_lat_lon.csv", header = T, row.names = 1, fileEncoding = "CP932") hc <- hclust(dist(JPll), method = "complete") # 引数は距離行列 plot(hc) # 樹形図の作図 ``` ## 下揃えあり ```{r, filename="階層ありクラスタリング：下を揃える", fig.width=10} plot(hc, hang = -1) ``` ## 階層なしクラスタリング ## k-means [ こちらのリンク ](http://tech.nitoyon.com/ja/blog/2013/11/07/k-means/) 、もしくは[ このWebアプリケーション ](https://www.shinyapps.io/admin/#/application/11051496) を見ていただくと、k-meansをよりイメージしやすいかと思います。WebアプリケーションはRで、`if(!require(shiny)){install.packages("shiny")};runGitHub("kmeans_animated", "sb8001at")` を実行することでも利用することができます。## k-mean：下グラフ描画のスクリプト ```{r, eval=FALSE} set.seed(0) # データにはirisのがく片データを使用 data <- iris[, 1:2] colnames(data) <- c("x", "y") # クラスターの中心の初期値 center_means <- data.frame( x = rnorm(3, mean(data[ ,1]), sd(data[ ,1])), y = rnorm(3, mean(data[ ,2]), sd(data[ ,2])), clus = 1:3) # アウトプットを保存する変数 data_for_gganimate <- NULL for(i in 1:20){ # おそらくi=16ぐらいで収束している # 各中心からの距離の計算 dist_d <- data.frame( dist_1 = (data$x - center_means[1, 1]) ^ 2 + (data$y - center_means[1, 2]) ^ 2, dist_2 = (data$x - center_means[2, 1]) ^ 2 + (data$y - center_means[2, 2]) ^ 2, dist_3 = (data$x - center_means[3, 1]) ^ 2 + (data$y - center_means[3, 2]) ^ 2 ) # 最も近いクラスターを探す min_dist <- apply(dist_d, 1, min) clus <- case_when( dist_d$dist_1 == min_dist ~ 1, dist_d$dist_2 == min_dist ~ 2, dist_d$dist_3 == min_dist ~ 3, ) # グラフ用データを準備 data_temp <- data.frame(data, clus, xcenter = center_means[clus,]$x, ycenter = center_means[clus,]$y, timestate = 2 * i - 1) data_for_gganimate <- rbind(data_for_gganimate, data_temp) # 重心でクラスターの中心を更新する center_means2 <- data_temp |> group_by(factor(clus)) |> summarise(x = mean(x), y = mean(y)) center_means$x <- center_means2$x center_means$y <- center_means2$y # グラフ用データを準備 data_temp <- data.frame(data, clus, xcenter = center_means[clus,]$x, ycenter = center_means[clus,]$y, timestate = 2 * i) data_for_gganimate <- rbind(data_for_gganimate, data_temp) } pacman::p_load(gganimate) # gganimateでグラフ化 p <- ggplot()+ geom_point(data = data_for_gganimate, aes(x = x, y = y, color = factor(clus), size = 2), alpha = 0.5) + geom_point(data = data_for_gganimate, aes(x = xcenter, y = ycenter, color = factor(clus), size = 3), alpha = 0.7)+ geom_segment(data = data_for_gganimate, aes(x = x, y = y, xend = xcenter, yend = ycenter, color = factor(clus)), alpha = 0.5)+ theme(legend.position = "none")+ transition_states(timestate, transition_length = 1, state_length = 1, wrap = TRUE)+ ease_aes('linear') anim_save("./image/p_kmeans.gif", animate(p, renderer = gifski_renderer())) ```  `kmeans` 関数を用いて行うことができます。`kmeans` 関数の引数はデータフレームとクラスターの数（`centers` ）です。以下の例では、`iris` のデータを3つのクラスターに分けています。k-meansのクラスター中心の初期値は乱数で決められますので、クラスターの結果は計算のたびに少しずつ変化します。また、初期値を特定の値として`centers` 引数に行列で与えることもできます。```{r, filename="kmeans関数でクラスタリング"} ## k-means set.seed(0) # k-meansは乱数計算なので、シードを設定 k.cent3 <- kmeans(iris[ ,1:4], centers = 3) # 3つのクラスタに分ける k.cent3$cluster # 計算結果 ``` ```{r, filename="クラスターの中心を指定"} k.cents <- kmeans(iris[ ,1:4], centers = iris[c(1, 51, 101), 1:4]) # 初期値を各種の始めの値に設定 k.cents$cluster ``` `iris` の`Species` の関係を示しています。k-meansでクラスター分けした結果が`Species` とほぼ重なることを確認できます。```{r, filename="クラスターとirisの種の比較", warning=FALSE} spvec <- iris$Species |> levels() k.cent3$cluster <- k.cent3$cluster - 1 k.cent3$cluster <- if_else(k.cent3$cluster == 0, 3, k.cent3$cluster) cbind(iris, cluster = spvec[k.cent3$cluster]) |> gather(tag, value, 5:6) |> ggplot(aes(x = Sepal.Length, y = Sepal.Width, color = value))+ geom_point(size = 3)+ facet_wrap(~ tag) ``` ## c-means [ `e1071` ](https://cran.r-project.org/web/packages/e1071/index.html) パッケージ [ @e1071_bib ] の`cmeans` 関数で行うことができます。`cmeans` 関数も`kmeans` 関数と同様に、データフレームと中心の数（`centers` ）を引数に取ります。クラスターの中心の初期値は行列で与えることもできます。```{r, filename="c-means法", warnings=FALSE, message=FALSE} pacman::p_load(e1071) # クラスタリングの計算（中心を指定） iris.c <- cmeans(iris[,1:4], centers = iris[c(1, 51, 101), 1:4]) # クラスターとirisの種の比較 cbind(iris, cluster = spvec[iris.c$cluster]) |> gather(tag, value, 5:6) |> ggplot(aes(x = Sepal.Length, y = Sepal.Width, color = value))+ geom_point(size = 3)+ facet_wrap(~ tag) ``` ## スペクトラルクラスタリング [ `kernlab` ](https://cran.r-project.org/web/packages/kernlab/index.html) パッケージ [ @kernlab_bib1; @kernlab_bib2 ] の`specc` 関数を用いて行うことができます。`specc` 関数の使い方は`kmeans` 関数や`cmeans` 関数とほぼ同じで、データフレームと`centers` を引数に取ります。また、返り値を`plot` 関数の引数にすることで、簡単にクラスターをグラフにすることもできます。```{r, filename="スペクトラルクラスタリング"} pacman::p_load(kernlab) # データにはkernlabのspiralsを用いる data(spirals) head(spirals) sc <- specc(spirals, centers = 2) # スペクトラルクラスタリング plot(spirals, col = sc) # 結果の表示（色がクラスターを示す） ``` `iris` のクラスタリングに用いると、3クラスターを指定しても、ほとんど2クラスターしか出てこなくなってしまいます。クラスター同士の距離が近く、接していると正しくクラスタリングできないようです。スペクトラルクラスタリングを用いる時にはデータの構造をよく確認する必要があります。```{r, filename="irisでスペクトラルクラスタリング"} # 3クラスタとなるよう設定 sc.iris <- specc(iris[,1:4] |> as.matrix(), centers = 3) # 結果に1がほとんどなく、ほぼ2クラスターになっている sc.iris@.Data ``` ## 主成分分析 [ 26章 ](./chapter26.html) で説明した主成分回帰では説明変数同士が相関しているときに起こる問題、多重共線性を避けることができます。## 主成分分析における主成分の軸のイメージ（スクリプト） ```{r, eval=FALSE} iris_temp <- iris |> filter(Species != "setosa") |> _[ ,1:3] iris.pc2 <- prcomp(iris_temp) pc1 <- iris.pc2$x[ ,1] pc2 <- iris.pc2$x[ ,2] pc3 <- iris.pc2$x[ ,2] center1 <- iris.pc2$center[1] center2 <- iris.pc2$center[2] center3 <- iris.pc2$center[3] iris_temp$pc1_Sepal.Length <- pc1 * iris.pc2$rotation[1, 1] + center1 iris_temp$pc1_Sepal.Width <- pc1 * iris.pc2$rotation[2, 1] + center2 iris_temp$pc1_Petal.Length <- pc1 * iris.pc2$rotation[3, 1] + center3 iris_temp$pc2_Sepal.Length <- pc2 * iris.pc2$rotation[1, 2] + center1 iris_temp$pc2_Sepal.Width <- pc2 * iris.pc2$rotation[2, 2] + center2 iris_temp$pc2_Petal.Length <- pc2 * iris.pc2$rotation[3, 2] + center3 iris_temp$pc3_Sepal.Length <- pc3 * iris.pc2$rotation[1, 3] + center1 iris_temp$pc3_Sepal.Width <- pc3 * iris.pc2$rotation[2, 3] + center2 iris_temp$pc3_Petal.Length <- pc3 * iris.pc2$rotation[3, 3] + center3 pacman::p_load(plotly) fig <- plot_ly( iris_temp, x =~ Sepal.Length, y =~ Sepal.Width, z =~ Petal.Length, type = "scatter3d", mode = "markers", name = 'data', marker = list(size = 2)) |> add_trace( x =~ pc1_Sepal.Length, y =~ pc1_Sepal.Width, z =~ pc1_Petal.Length, type = "scatter3d", mode = "markers+lines", name = 'PC1', marker = list(color = "red", size = 0.1), line = list(color = "red", width = 5)) |> add_trace( x =~ pc2_Sepal.Length, y =~ pc2_Sepal.Width, z =~ pc2_Petal.Length, type = "scatter3d", mode = "markers+lines", name = 'PC2', marker = list(color = "red", size = 0.1), line = list(color = "blue", width = 5)) |> add_trace( x =~ pc3_Sepal.Length, y =~ pc3_Sepal.Width, z =~ pc3_Petal.Length, type = "scatter3d", mode = "markers+lines", name = 'PC3', marker = list(color="red", size=0.1), line = list(color="green", width=5)) fig ``` ```{r, echo=FALSE} iris_temp <- iris |> filter(Species != "setosa") |> _[,1:3] iris.pc2 <- prcomp(iris_temp) pc1 <- iris.pc2$x[,1] pc2 <- iris.pc2$x[,2] pc3 <- iris.pc2$x[,2] center1 <- iris.pc2$center[1] center2 <- iris.pc2$center[2] center3 <- iris.pc2$center[3] iris_temp$pc1_Sepal.Length <- pc1 * iris.pc2$rotation[1, 1] + center1 iris_temp$pc1_Sepal.Width <- pc1 * iris.pc2$rotation[2, 1] + center2 iris_temp$pc1_Petal.Length <- pc1 * iris.pc2$rotation[3, 1] + center3 iris_temp$pc2_Sepal.Length <- pc2 * iris.pc2$rotation[1, 2] + center1 iris_temp$pc2_Sepal.Width <- pc2 * iris.pc2$rotation[2, 2] + center2 iris_temp$pc2_Petal.Length <- pc2 * iris.pc2$rotation[3, 2] + center3 iris_temp$pc3_Sepal.Length <- pc3 * iris.pc2$rotation[1, 3] + center1 iris_temp$pc3_Sepal.Width <- pc3 * iris.pc2$rotation[2, 3] + center2 iris_temp$pc3_Petal.Length <- pc3 * iris.pc2$rotation[3, 3] + center3 pacman::p_load(plotly) fig <- plot_ly( iris_temp, x=~Sepal.Length, y=~Sepal.Width, z=~Petal.Length, type="scatter3d", mode="markers", name = 'data', marker=list(color="orange", size=3)) |> add_trace( x=~pc1_Sepal.Length, y=~pc1_Sepal.Width, z=~pc1_Petal.Length, type="scatter3d", mode="markers+lines", name = 'PC1', marker=list(color="red", size=0.1), line=list(color="red", width=5)) |> add_trace( x=~pc2_Sepal.Length, y=~pc2_Sepal.Width, z=~pc2_Petal.Length, type="scatter3d", mode="markers+lines", name = 'PC2', marker=list(color="red", size=0.1), line=list(color="blue", width=5)) |> add_trace( x=~pc3_Sepal.Length, y=~pc3_Sepal.Width, z=~pc3_Petal.Length, type="scatter3d", mode="markers+lines", name = 'PC3', marker=list(color="red", size=0.1), line=list(color="green", width=5)) fig ``` ```{r, echo=FALSE} data.frame(pc1, pc2) |> ggplot(aes(x=pc1, y=pc2))+ geom_point(size=2, color="orange")+ geom_hline(yintercept=0, color="red")+ geom_vline(xintercept=0, color="blue")+ labs(title="PC1とPC2の2次元グラフ") ``` [ 論文 ](https://www.nature.com/articles/jhg2012114#citeas) が発表されています[ @1523106605673180672 ] 。この論文を見ていただければわかる通り、主成分分析を用いれば膨大なSNPのデータを2次元に集約し、わかりやすい表現で説明することができます。### Rで主成分分析 `prcomp` 関数で主成分分析の計算を行うことができます。`prcomp` 関数は引数に行列やデータフレームを取ります。`prcomp` 関数の返り値は各主成分方向の標準偏差（standard deviations）と、回転（Rotation）です。回転は、データを主成分に変換するときの係数を示します。```{r, filename="主成分分析"} iris.pc <- prcomp(iris[ ,1:4]) iris.pc ``` `prcomp` 関数の返り値を`summary` 関数の引数に取ると、標準偏差、分散の配分（Proportion of Variance）、積算の分配の配分（Cumulative Proportion）が示されます。このうち、分散の配分はその主成分によってばらつきをどれだけ説明できているかを示すものです。下の例では、PC1で92.5%程度のばらつきを示すことができている、つまり、PC1でデータの差を十分説明できていることがわかります。```{r, filename="summaryの結果"} iris.pc |> summary() ``` `$x` で呼び出すことができます。以下のように、`iris` の種が同じ、つまりよく似たデータであれば、比較的近い位置に変換後のデータが集まるのがわかります。```{r, filename="第一主成分と第二主成分をプロット"} iris.pc <- prcomp(iris[ ,1:4]) irispcd <- iris.pc$x |> as.data.frame() irispcd$Species <- iris$Species ggplot(irispcd, aes(x = PC1, y = PC2, color = Species)) + geom_point(size = 2) ``` ## 因子分析 `factanal` 関数を用いて行うことができます。以下の例では、鹿児島大学が[ 成績サンプルデータ ](https://estat.sci.kagoshima-u.ac.jp/data/cgi-bin/data/whats_data/data/img/932722923_9821.xls) として提供しているデータを用いて因子分析を行っています。`factanal` 関数はデータフレームまたは行列を第一引数に取り、因子の数（`factors` ）、回転の方法（`rotation` ）、算出するスコアのタイプ（`scores` ）を引数として設定して用います。`factors` には、出力として得たい因子の数を指定します。以下の例では、理系・文系科目の得意・不得意を評価する因子を作成するため、因子数を`2` としています。`rotation` には`"varimax"` か`"promax"` のどちらかを選択するのが一般的です。`varimax` 回転はデータをそのまま回転させる方法（直交回転）、`promax` 回転はデータの軸の角度を変換させて回転させる方法（斜交回転）です。`promax` 回転では因子間に相関が生じるのが特徴です。`scores` には計算して得られる数値（因子スコア）の計算方法を指定します。方法には`"none"` 、`"Bartlett"` 、`"regression"` のいずれかを選択します。`"none"` を指定するとスコアが算出されないので、`"Bartlett"` か`"regression"` のいずれかを選択します。`factanal` 関数の返り値のうち、Loadingが各データの寄与率を示す値となります。以下のvarimax変換の例では、Factor1（因子1）では数学、物理、化学、Factor2では国語、英語の値が高いため、Factor1は理系科目の、Factor2は文系科目の評価を反映していることがわかります。```{r, filename="因子分析：バリマックス変換"} m1 <- read.csv("./data/testresult.csv", fileEncoding = "CP932") m1v <- factanal(m1[ ,2:6], factors = 2, scores = "Bartlett") # varimax変換 m1v ``` ```{r, filename="因子分析：プロマックス変換"} m1p <- factanal(m1[,2:6], factors = 2, rotation = "promax", scores = "Bartlett") # promax変換 m1p ``` `$scores` で呼び出すことができます。varimax回転でもpromax回転でも、概ね理系・文系科目の得意さを因子で評価できています。## varimax回転 ```{r} <- m1v$ scores |> as.data.frame () # 計算したスコアをデータフレームに変換 $ student <- 1 : nrow (m1vs) # 生徒のIDを追加 ggplot (m1vs, aes (x = Factor1, y = Factor2, color = student, label = student))+ geom_text (size = 3 )+ theme (legend.position = "none" )+ labs (title = "バリマックス回転" , x = "理系度" , y = "文系度" )``` ## promax回転 ```{r} <- m1p$ scores |> as.data.frame ()$ student <- 1 : nrow (m1ps)ggplot (m1ps, aes (x = Factor1, y = Factor2, color = student, label = student))+ geom_text (size = 3 )+ theme (legend.position = "none" )+ labs (title = "プロマックス回転" , x = "理系度" , y = "文系度" )``` ## 多次元尺度法 `dist` 関数で距離行列に変換し、距離行列を多次元尺度法で位置の情報に変換しています。```{r, filename="緯度・経度を距離行列に変換"} JPD <- read.csv("./data/pref_lat_lon.csv", header = T, fileEncoding = "CP932") prefecture <- JPD[,1] %>% unlist() JPD <- JPD[,-1] %>% dist() # 距離行列を計算 ``` `MASS` パッケージの`sammon` 関数を用いて行うことができます。`sammon` 関数の引数は距離行列です。距離行列を`sammon` 関数で変換すると、`x` と`y` という、位置を示す2つの変数が求まります。```{r} :: p_load (MASS)<- JPD |> sammon () |> _$ pointshead (JPpos)``` `sammon` 関数の返り値をプロットしたものを示します。`x` と`y` の単位はありませんが、都道府県の位置関係を正確に反映していることが見て取れると思います。ただし、距離行列には方角のデータが含まれていないため、東西南北の方向は回転したり反転したりすることになります。```{r, filename="多次元尺度法の結果をプロットする"} JPpos <- as.data.frame(JPpos) colnames(JPpos) <- c("x", "y") JPpos$prefecture <- prefecture ggplot(JPpos, aes(x = -x, y = -y, label = prefecture, color = prefecture))+ geom_text(size = 3)+ theme(legend.position = "none") ``` `cmdscale` 関数を用いても多次元尺度法の計算を行うことができます。```{r, filename="cmdscale関数で多次元尺度法"} JPpos2 <- JPD %>% cmdscale %>% as.data.frame colnames(JPpos2) <- c("x", "y") JPpos2$prefecture <- prefecture ggplot(JPpos2, aes(x = x, y = y, label = prefecture, color = prefecture))+ geom_text(size = 3)+ theme(legend.position = "none") ``` 
