多変量解析

またデータはirisデータセットに戻る。ここからはアドバンスのため、余力がある人やこれらの統計が必要になった人だけ行う形で問題ない。

## Rに収載されているデータ"iris"を呼び出して"iris_data"に格納
## なお、以降このデータを使うため、一度Rstudioを閉じた場合は一番最初から順々に動かして、各データを取ると良い。
iris_data = iris
## 最初の6行のみを表示
head(iris_data)
##   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
## setosa, versicolor, virginicaの各グループごとに分けたデータを作成
iris_sub_set = subset(x = iris_data, subset = iris_data$Species == "setosa")
iris_sub_ver = subset(x = iris_data, subset = iris_data$Species == "versicolor")
iris_sub_vir = subset(x = iris_data, subset = iris_data$Species == "virginica")

ロジスティック回帰

本書を見た場合、なぜ線形回帰分析で重回帰分析を紹介しないのかと思うかもしれない。しかし中でも書いたようにロジスティック回帰の本質は線形回帰であり、そのパラメータを変えるだけで実施可能である。また単純な線形回帰よりも拡張された一般化線形モデルのほうが様々な分布を扱えるため、本稿ではロジスティック回帰を用いて学習していく。

分布名 family(デフォルトのリンク関数)
二項分布 binomial(link = “logit”)
正規分布（ガウス分布） gaussian(link = “identity”)
ガンマ分布 Gamma(link = “inverse”)
逆ガウス分布（ワルド分布） inverse.gaussian(link = “1/mu^2”)
ポアソン分布 poisson(link = “log”)
疑似尤度 quasi(link = “identity”, variance = “constant”)
疑似二項分布 quasibinomial(link = “logit”)
疑似ポアソン分布 quasipoisson(link = “log”)

ロジスティック単回帰

まずは単回帰分析ということで「がくの長さ（Sepal.Length）でアヤメの2品種（versicolorとvirginica）を予測」をしてみる。versicolor = 0, virginica = 1とすると二項分布になる。

## 2品種のデータを結合
ver_vir_data = rbind(iris_sub_ver, iris_sub_vir)
# versicolorを0、virginicaを1というダミー数字を割り当てた列Sを作成
ver_vir_data$S = c(rep(0, 50),rep(1,50))

## この分布を描画
library(ggplot2)
## Warning: パッケージ 'ggplot2' はバージョン 4.2.3 の R の下で造られました
logit_ver_vir_gg = ggplot(data = ver_vir_data, 
                       aes(x = Sepal.Length, y = S)) +
  geom_count(aes(color = Species)) +
  theme_classic()
logit_ver_vir_gg

続いてロジスティック回帰を行ってみる。ロジスティック回帰は一般化線形モデルを実行できるglm関数を用いて行う。

## ロジスティック単回帰を実施
## versicolorかvirginica（つまり0か1か）をがくの長さ（Sepal.Length）から求める確率を算出する
logit_ver_vir = glm(S ~ Petal.Width, family = binomial(link = "logit"),
                    data = ver_vir_data)
summary(logit_ver_vir)
## 
## Call:
## glm(formula = S ~ Petal.Width, family = binomial(link = "logit"), 
##     data = ver_vir_data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.13868  -0.16468  -0.00929   0.13000   2.46891  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -21.126      4.597  -4.596 4.31e-06 ***
## Petal.Width   12.947      2.843   4.554 5.25e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 138.629  on 99  degrees of freedom
## Residual deviance:  33.421  on 98  degrees of freedom
## AIC: 37.421
## 
## Number of Fisher Scoring iterations: 7

Coefficients（係数）を見るとPetal.Width の係数 は 12.947 で、この値は Petal.Width が1単位増加したときに、目的変数 S のログオッズがどれだけ変わるかを示している。この場合、Petal.Width が1増えるごとに、ログオッズが約12.95増加するということを示している。 P値（Pr(>|z|)） は 5.25e-06 で、非常に小さいことから、Petal.Width は統計的に有意であることが示されている。
confint()という関数で、パラメータの信頼区間を算出できる。これを用いてデータを描画してみる。

confint(logit_ver_vir)
## Waiting for profiling to be done...
##                  2.5 %    97.5 %
## (Intercept) -32.144431 -13.68017
## Petal.Width   8.355608  19.79404
## 上で書いた図にgeom_smooth関数を使ってglmを実施した結果で線を引いている
logit_ver_vir_gg = ggplot(data = ver_vir_data, 
                       aes(x = Sepal.Length, y = S)) +
  geom_count(aes(color = Species)) +
  geom_smooth(method = "glm",　
              method.args = list(family = binomial(link = "logit"))) +
  theme_classic()
logit_ver_vir_gg
## `geom_smooth()` using formula = 'y ~ x'

ロジスティック重回帰

単回帰の結果では有意なものの、お互いをきれいに分類することは難しい。では「がくの長さと幅からアヤメの三品種（setosaとversicolorとvirginica）を分類」してみる。まずがくの長さと幅の関係を表してみる。

iris_PL_PW_gg = ggplot(data = iris_data, 
                       aes(x = Sepal.Length, y = Sepal.Width)) +
  geom_point(aes(color = Species)) +
  theme_classic()
iris_PL_PW_gg

ではこれを識別する多項ロジスティック回帰を行っていく。アヤメの三品種（setosaとversicolorとvirginica）を３つの組み合わせで0と1に分類してglmを行っても良いが、書籍のスペースに限りがあり、コードも冗長になるため、nnetパッケージのmultinom関数を用いて解析していく。

## nnetパッケージの呼び出し
library(nnet)
## Warning: パッケージ 'nnet' はバージョン 4.2.3 の R の下で造られました
## 多項ロジスティック回帰を行っている
## 多項分布専用のため、family の指定はない
iris_PL_PW_multi = multinom(Species ~ Sepal.Length + Sepal.Width, data = iris_data)
## # weights:  12 (6 variable)
## initial  value 164.791843 
## iter  10 value 62.715967
## iter  20 value 59.808291
## iter  30 value 55.445984
## iter  40 value 55.375704
## iter  50 value 55.346472
## iter  60 value 55.301707
## iter  70 value 55.253532
## iter  80 value 55.243230
## iter  90 value 55.230241
## iter 100 value 55.212479
## final  value 55.212479 
## stopped after 100 iterations
## 結果の確認
summary(iris_PL_PW_multi)
## Call:
## multinom(formula = Species ~ Sepal.Length + Sepal.Width, data = iris_data)
## 
## Coefficients:
##            (Intercept) Sepal.Length Sepal.Width
## versicolor   -92.09924     40.40326   -40.58755
## virginica   -105.10096     42.30094   -40.18799
## 
## Std. Errors:
##            (Intercept) Sepal.Length Sepal.Width
## versicolor    26.27831     9.142717    27.77772
## virginica     26.37025     9.131119    27.78875
## 
## Residual Deviance: 110.425 
## AIC: 122.425
# モデルの係数に対する有意性のテスト (z値の計算)
z_values <- summary(iris_PL_PW_multi)$coefficients / summary(iris_PL_PW_multi)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z_values)))

# z値とp値の表示
print("Z-values:")
## [1] "Z-values:"
print(z_values)
##            (Intercept) Sepal.Length Sepal.Width
## versicolor   -3.504763     4.419175   -1.461155
## virginica    -3.985588     4.632613   -1.446197
print("P-values:")
## [1] "P-values:"
print(p_values)
##             (Intercept) Sepal.Length Sepal.Width
## versicolor 4.570135e-04  9.90786e-06    0.143973
## virginica  6.731329e-05  3.61079e-06    0.148122

最終的な結果からがくの長さのP値は低いものの、がくの幅のP値は低くない。summaryの結果から、以下のように考える。
1. Coefficients（回帰係数）
Sepal.Length（がくの長さ）の係数：

  • versicolor の係数は 40.40326 で、Sepal.Length が1単位増加すると、versicolor に分類される対数オッズが約40.4増加することを意味する。

  • virginica の係数は 42.30094 で、同様に Sepal.Length の増加が virginica に分類されるオッズを高めることを示している

Sepal.Width（がくの幅）の係数：

  • versicolor の係数は -40.58755 で、Sepal.Width が1単位増加すると versicolor に分類される対数オッズが約40.59減少することを意味する。つまり、Sepal.Width が大きくなるほど versicolor に分類される可能性は低くなる。

  • virginica の係数も同様に -40.18799 で、Sepal.Width が大きくなると virginica に分類されるオッズも低くなる。

2. Std. Errors（標準誤差）

  • Sepal.Length の標準誤差は、versicolor で約 9.14、virginica で約 9.13

  • Sepal.Width の標準誤差は、versicolor で約 27.78、virginica で約 27.79

標準誤差の値が大きい（特に Sepal.Width）ため、この説明変数に対する回帰係数の信頼性がやや低いことがわかる。つまり、がくの幅が分類に与える影響は不確実性が高いことが示唆される。

階層的クラスタリング

ロジスティック回帰ではあらかじめ正解の情報を与えて実施していた。しかし世の中には事前に正解情報が分からない場合も多い。そこで階層的クラスタリングを用いてirisデータセットの各情報の類似度から分類を行ってみる。まず距離行列を作成し、その後各データ間の距離を求めている。

## 品種情報を与えてしまうとそれで分類されてしまうため、がくと花びらの情報のみを使用する
## methodでユークリッド距離を与えている
iris_dist = dist(x = iris_data[, 1:4], method = "euclidean")

## 階層的クラスタリングを実施
## 今回はウォード法(ward.D2)を使用している
iris_hc = hclust(iris_dist, method = "ward.D2")

## 今回は視覚化しやすくするためにggplotのデンドログラム用パッケージを用いている
library(ggdendro)
## Warning: パッケージ 'ggdendro' はバージョン 4.2.3 の R の下で造られました
iris_dd = dendro_data(as.dendrogram(iris_hc), type = "rectangle")

# クラスタリング結果のラベルを抽出し、元の品種と対応付ける
iris_dd$labels$Species = iris_data[as.numeric(iris_dd$labels$label), "Species"]

# 品種ごとに色を割り当ててデンドログラムを描画
iris_hc_gg = ggplot(data = segment(iris_dd)) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_text(data = iris_dd$labels, aes(x = x, y = y, label = label, color = Species), hjust = 1, angle = 90, size = 3) +
  scale_color_manual(values = c("setosa" = "red", "versicolor" = "blue", "virginica" = "green")) +
  labs(title = "Hierarchical Clustering of Iris Data", x = "", y = "") +
  theme_classic()
iris_hc_gg

この結果を見ると教師なし学習でもsetosaは良く分類されており、一つで独立したクラスターを形成している。一方でversicolorとvirginicaもそこそこ分類されているが、やはり一部混ざってしまっており、がくと花びらの情報だけで分類することは難しいことがよくわかる。

主成分分析（Principal Component Analysis; PCA）

では続いてirisデータセットを主成分分析を用いて分類してみる。このままデータを使うとSepal.Lengthが最も大きく、それに引っ張られてしまうため標準化を行う。主成分分析を行うprcomp関数には標準化用のオプションがあるため、それを利用する。まずはデータの第一主成分（PC1）と第二主成分（PC2）を使って視覚化してみる。

## 主成分分析を実行
## center = TRUE, scale. = TRUEとすると標準化されて実行される
iris_pca = prcomp(x = iris_data[, 1:4], center = TRUE, scale. = TRUE)

# 主成分得点をデータフレームに変換し、元の品種情報を追加
iris_pca_data = data.frame(iris_pca$x, Species = iris_data$Species)

## ggplot2を使って主成分分析結果を可視化
iris_pca12_gg = ggplot(data = iris_pca_data, 
                     aes(x = PC1, y = PC2, color = Species)) +
  geom_point(size = 3) +
  labs(title = "PCA of Iris Dataset", x = "PC1", y = "PC2") +
  scale_color_manual(values = c("setosa" = "red", "versicolor" = "blue", "virginica" = "green")) +
  theme_classic()

iris_pca12_gg

まだ分類が上手く行かないため、各主成分毎の寄与率を求めてみる。

## 主成分分析を実行
## center = TRUE, scale. = TRUEとすると標準化されて実行される
iris_pca = prcomp(x = iris_data[, 1:4], center = TRUE, scale. = TRUE)


## 各主成分の寄与率を計算
pca_variance = iris_pca$sdev^2  # 主成分ごとの分散
pca_variance_ratio = pca_variance / sum(pca_variance)  # 寄与率

## 寄与率をデータフレームに変換
pca_df = data.frame(Principal_Component = paste0("PC", 1:length(pca_variance_ratio)), Variance_Explained = pca_variance_ratio)

## 寄与率の棒グラフを描画
pca_df_gg = ggplot(data = pca_df, 
                   aes(x = Principal_Component, y = Variance_Explained)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = round(Variance_Explained, 2)), vjust = -0.3) +
  labs(title = "Variance Explained by Principal Components", x = "Principal Components", y = "Proportion of Variance Explained") +
  theme_classic()

pca_df_gg

この結果を見るとPC1で73%、PC2で23%説明を付けられるが、残り5%がPC3とPC4で説明されることがわかる。そこで試しにPC3との関係性を書いてみる。

## 主成分分析を実行
## center = TRUE, scale. = TRUEとすると標準化されて実行される
iris_pca = prcomp(x = iris_data[, 1:4], center = TRUE, scale. = TRUE)

# 主成分得点をデータフレームに変換し、元の品種情報を追加
iris_pca_data = data.frame(iris_pca$x, Species = iris_data$Species)

## ggplot2を使ってPC1とPC3を可視化
iris_pca13_gg = ggplot(data = iris_pca_data, 
                     aes(x = PC1, y = PC3, color = Species)) +
  geom_point(size = 3) +
  labs(title = "PCA of Iris Dataset", x = "PC1", y = "PC3") +
  scale_color_manual(values = c("setosa" = "red", "versicolor" = "blue", "virginica" = "green")) +
  theme_classic()

iris_pca13_gg

## ggplot2を使ってPC2とPC3を可視化
iris_pca23_gg = ggplot(data = iris_pca_data, 
                     aes(x = PC2, y = PC3, color = Species)) +
  geom_point(size = 3) +
  labs(title = "PCA of Iris Dataset", x = "PC2", y = "PC3") +
  scale_color_manual(values = c("setosa" = "red", "versicolor" = "blue", "virginica" = "green")) +
  theme_classic()

iris_pca23_gg

この結果を見ると意外なことにPC1とPC3の方がよく分かれているように見える。このためPC1-3までのデータを用いれば少なくとも重心は3群に分けることができるが、versicolorとvirginicaを完全に分けることは難しいということがわかる。

興味がある方は3次元空間で見た場合を試してみてほしい。

##plotlyを使って三次元プロットを作成
library(plotly)
## Warning: パッケージ 'plotly' はバージョン 4.2.3 の R の下で造られました
## 
##  次のパッケージを付け加えます: 'plotly'
##  以下のオブジェクトは 'package:ggplot2' からマスクされています:
## 
##     last_plot
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     filter
##  以下のオブジェクトは 'package:graphics' からマスクされています:
## 
##     layout
plot_ly(data = iris_pca_data, x = ~PC1, y = ~PC2, z = ~PC3, color = ~Species, colors = c("red", "blue", "green")) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')),
         title = "3D PCA of Iris Dataset")

Heatmap

引き続きirisデータセットで行ってみる。色々なパッケージがあるが拡張性が高いComplexheatmapパッケージを紹介する。

## irisデータセットの数値データのみを抽出
## 数値データ: Sepal.Length, Sepal.Width, Petal.Length, Petal.Width
## Zスコア化（標準化）

iris_zscore <- as.data.frame(scale(iris_data[, 1:4]))
iris_zscore[iris_zscore < -2] = -2
iris_zscore[iris_zscore > 2] = 2

library(ComplexHeatmap)
##  要求されたパッケージ grid をロード中です
## ========================================
## ComplexHeatmap version 2.14.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## 
##  次のパッケージを付け加えます: 'ComplexHeatmap'
##  以下のオブジェクトは 'package:plotly' からマスクされています:
## 
##     add_heatmap
library(circlize)
## Warning: パッケージ 'circlize' はバージョン 4.2.3 の R の下で造られました
## ========================================
## circlize version 0.4.16
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================
## 品種ごとの注釈データを作成
ann_colors = list(
  Species = c(setosa = "red", versicolor = "blue", virginica = "green"))
ha = HeatmapAnnotation(Species = iris_data$Species, col = ann_colors)

## ヒートマップを作成
ComplexHeatmap::Heatmap(as.matrix(t(iris_zscore)), name = "Z-Score", 
                        top_annotation = ha, 
                        show_row_names = TRUE, 
                        clustering_distance_rows = "euclidean",
                        clustering_method_rows = "ward.D2",
                        show_column_names = TRUE,
                        clustering_distance_columns = "euclidean",
                        clustering_method_columns = "ward.D2",
                        column_names_gp = gpar(fontsize = 10),
                        cluster_rows = TRUE,  # 行クラスタリングを有効にする
                        cluster_columns = TRUE)  # 列クラスタリングを有効にする