Rstudioで Postgresqlからデータを読み込んで主成分分析する

# 概要
* [RPostgreSQLライブラリ](https://cran.r-project.org/web/packages/RPostgreSQL/index.html)でPostgresqlからRstudioへデータをインポート
* インポートしたデータを[tidyrライブラリ](https://cran.r-project.org/web/packages/tidyr/README.html)を使って主成分分析できるようにピボットする
* [prcompライブラリ](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/prcomp.html)で主成分分析する
    ![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/a86c9b63-bebf-60db-4f1e-fe5f8cb2e6b1.png)

# 環境
* RStudio
* R


# 手順
## Rstudioのインストール
* RstudioとRは下記でインストールする
 * https://www.rstudio.com/products/rstudio/download/#download
 * https://cran.r-project.org/

## デモデータ
* 全国のDPC病院の疾患別・手術別集計を使います
 * https://www.mhlw.go.jp/stf/shingi2/0000150723.html
 * 今回は、MDC01 (神経系疾患) のみを対象としました。 https://www.mhlw.go.jp/file/05-Shingikai-12404000-Hokenkyoku-Iryouka/0000023561.xls
* 今回の解析の目的は： 神経系疾患病院を分類する際の評価軸の検討としました。

## データ前処理
前記集計には、入院期間と患者件数データがあります。今回は、患者件数データを対象とします。
患者件数は様々な治療種類毎のデータが入っています。今回は事前に治療種類毎の単純合計値、主成分解析の解析対象としました。
また、さらに、有効なデータが少ない、病気および病院は除去しました。（詳細は省略）


![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/9ca73b55-450b-9afd-adc5-3a7e964aa4fb.png)

## データをPostgresqlへ投入
* dpcというユーザー名でdpcというデータベースとdpc_brainというテーブルを用意します
* dpc_brainテーブルには、hospital（病院名）とdisease（病気名）とpatient（患者件数合計）というカラムを用意します。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/ff6ee8ab-9af9-8572-7524-ac3d25a0412d.png)
* データを投入します。
  * 元のデータは、病院名, 病気1, 病気2,..というテーブル構造です。それをdpc_brainの構造に変換してimportします。（詳細割愛）
  * import後は下記の様になります
   ![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/6891dcda-bc3a-dd47-52e0-89da9ac28ab0.png)

## Postgresqlのデータを読み込む
* [PostgreSQLのテーブルをRのデータフレームに一瞬で流しこむ方法](https://qiita.com/junkonakajima/items/07ad8a227c6c586e1f99)の通りです。
* dbnameとuserとpasswordは今回のデータベースにあわせます。

```r
install.packages("RPostgreSQL")
require(RPostgreSQL)
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv,dbname="dpc",host="localhost",port=5432,user="dpc",password="")
dpc_brain<- dbReadTable(con,"dpc_brain")
```

## 分析しやすいようにデータを変換する
### 使用するデータだけを抽出
1列目はidなので、2列名以降（hospital, disease, patient）を選択して selected_dpc_brainに入れる

```r
selected_dpc_brain <- dpc_brain[,2:4]
```

### データをピボットする
主成分分析等、統計解析する際は、各列は変数ごとになっていることを想定してなされることがほとんどです。今回の場合もdisease列の各病気名を列名とするように、selected_dpc_brainを[ピボット](https://ja.wikipedia.org/wiki/%E3%83%94%E3%83%9C%E3%83%83%E3%83%88%E3%83%86%E3%83%BC%E3%83%96%E3%83%AB)する必要があります。
そのために[tidyr](https://cran.r-project.org/web/packages/tidyr/README.html)をつかいます。

```r
spreaded_dpc_brain <- tidyr::spread(selected_dpc_brain, key = "disease", value = "patient")
```
その他、tidyrの細かな使い方は、こちらの[tidyr::gather( )とtidyr::spread( )でデータフレームを自在に変形する](https://qiita.com/wakuteka/items/a16783571794f12af2ac)を参照してください。

### ラベルをつける
普通に読むこむと各行のラベルは連番（1, 2, 3...）です。これでも解析上は問題ないのですが、結果をみるとき不便なので、病院名をラベル（行名、rownames）にします。
1列目が病院名なので、それを行名として投入します。

```r
rownames(spreaded_dpc_brain) <- spreaded_dpc_brain[,1]
```

## 解析する
### グラフ内で日本語が表示できるようにする

こちらの[ページ](https://www.karada-good.net/analyticsr/r-58)からコードを抜粋し、実行します。これをやらないと文字が豆腐になってしまいます。

```r
par(family = "HiraKakuProN-W3")
```

### 基礎データ
summaryを見ます。

```r
summary(spreaded_dpc_brain[,2:11])
 くも膜下出血、破裂脳動脈瘤    てんかん       パーキンソン病      脳血管障害         脳梗塞          脳腫瘍       脳脊髄の感染を伴う炎症
 Min.   :  0.000            Min.   :   0.00   Min.   :  0.000   Min.   :  0.00   Min.   :  0.0   Min.   :  0.00   Min.   :  0.000       
 1st Qu.:  0.000            1st Qu.:   0.00   1st Qu.:  0.000   1st Qu.:  0.00   1st Qu.: 19.0   1st Qu.:  0.00   1st Qu.:  0.000       
 Median :  0.000            Median :  16.00   Median :  0.000   Median :  0.00   Median : 54.0   Median :  0.00   Median :  0.000       
 Mean   :  4.446            Mean   :  30.97   Mean   :  6.996   Mean   : 16.46   Mean   : 99.6   Mean   : 19.68   Mean   :  6.393       
 3rd Qu.:  0.000            3rd Qu.:  43.00   3rd Qu.:  0.000   3rd Qu.: 13.00   3rd Qu.:151.2   3rd Qu.: 16.00   3rd Qu.: 11.000       
 Max.   :157.000            Max.   :2428.00   Max.   :363.000   Max.   :649.00   Max.   :882.0   Max.   :561.00   Max.   :123.000       
 非外傷性硬膜下血腫 非外傷性頭蓋内血腫（非外傷性硬膜下血腫以外） 未破裂脳動脈瘤  
 Min.   :  0.000    Min.   :  0.00                               Min.   :  0.00  
 1st Qu.:  0.000    1st Qu.:  0.00                               1st Qu.:  0.00  
 Median :  0.000    Median : 10.00                               Median :  0.00  
 Mean   :  5.188    Mean   : 25.65                               Mean   : 14.72  
 3rd Qu.:  0.000    3rd Qu.: 40.00                               3rd Qu.: 10.00  
 Max.   :109.000    Max.   :308.00                               Max.   :721.00  
```

### 相関係数行列
各病気間の相関係数をみます。

```r
pairs(spreaded_dpc_brain[,2:10])
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/8b59b202-c94c-c64d-a5d4-2f3eb8e75e20.png)

```r
library(corrr)
dpc_brain_correlate <- correlate(spreaded_dpc_brain[,2:11])
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/16fdb2eb-200c-8813-b694-1c1d8c92ce57.png)

### 主成分分析
prcompを使うだけ。

```r
pca <- prcomp(spreaded_dpc_brain[,c(2:11)],scale=T)
```

summaryを見る。

```r
> summary(pca)
Importance of components:
                          PC1    PC2     PC3    PC4     PC5     PC6     PC7     PC8     PC9    PC10
Standard deviation     2.1683 1.0566 0.99355 0.8526 0.82848 0.79096 0.69449 0.60600 0.48495 0.26669
Proportion of Variance 0.4702 0.1116 0.09871 0.0727 0.06864 0.06256 0.04823 0.03672 0.02352 0.00711
Cumulative Proportion  0.4702 0.5818 0.68051 0.7532 0.82185 0.88441 0.93265 0.96937 0.99289 1.00000
```

Cumulative Propotionをみると主成分4ぐらいまで軸がありそうです。


結果をグラフにプロットしてみます。

```r
biplot(pca)
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/124445/a86c9b63-bebf-60db-4f1e-fe5f8cb2e6b1.png)

赤い矢印の軸をみると、脳の病気に関して、病院を分類する際は、4種類ぐらい評価軸がありそうなことがわかりました。


# Extra 主成分分析はRubyでもできるよ!

主成分分析をrubyでやろうとしても、デファクトスタンダードなライブラリがありません。gem pcaで調べると2つぐらいでますが、スター数が少ないので、不安・・・

しかし、主成分分析の実装は簡単です。計算はたった2つなので。

1. 相関係数行列（標準化されたデータの分散共分散行列）を計算する！
2. 上記行列の固有値ベクトルを計算する！

（理論は、　主成分分析 固有値ベクトル とかでググってください）

1,2いずれも普通の行列計算ソフトであれば実装されているものです。
今回は、演算が速いとされている日本製のライブラリ[Numo](https://github.com/ruby-numo)を使った例を示します。（標準の[matrix](https://docs.ruby-lang.org/ja/latest/class/Matrix.html)でもできます。）


まず、正解例としてRを使って計算します。（なお、計算を簡単にするため今回は相関係数行列ではなく、分散共分散行列を使った例です。）


```r
> a <- data.frame(x1=c(1,4,2), x2=c(4,8,3))
> var(a)
         x1  x2
x1 2.333333 3.5
x2 3.500000 7.0
> pca <- prcomp(a,scale=F)
> pca
Standard deviations (1, .., p=2):
[1] 2.9787822 0.6783732

Rotation (n x k) = (2 x 2):
         PC1        PC2
x1 0.4718579  0.8816746
x2 0.8816746 -0.4718579
```

次にruby の例を示します。


```irb
irb(main):002:0> require 'numo/narray'
=> true
irb(main):003:0> require 'numo/linalg'
=> true
irb(main):004:0> a = Numo::DFloat.cast([[1,4,2], [4,8,3]])
=> Numo::DFloat#shape=[2,3]
[[1, 4, 2], 
 [4, 8, 3]]
irb(main):005:0> a.cov
=> Numo::DFloat#shape=[2,2]
[[2.33333, 3.5], 
 [3.5, 7]]
irb(main):007:0> Numo::Linalg.eig(a.cov)
=> [Numo::DComplex#shape=[2]
[0.46019+0i, 8.87314+0i], nil, Numo::DComplex#shape=[2,2]
[[-0.881675+0i, -0.471858+0i], 
 [0.471858+0i, -0.881675+0i]]]

```

値が一致するので、主成分分析ができることがわかりました。

# 参考
* 主成分分析
    * https://logics-of-blue.com/principal-components-analysis/
    * https://www.cis.doshisha.ac.jp/mjin/R/24/24.html
* DPC
    * https://www.mhlw.go.jp/stf/shingi2/0000150723.html