1. 先週のショート課題から: 19名


2. 主成分分析

 いくつか(p個)の変量の値を情報の損失をできるだけ少なくして、 少数変量(m個、m<p)の総合的指標(主成分)で代表させる方法として 主成分分析(Principal Component Analysis, PCA)と 因子分析(Factor Analysis, FA)がある。 いくつかのテストの成績を総合した総合的成績、 いろいろな症状を総合した総合的な重症度、 種々の財務指標に基づく企業の評価 等を求めたいといった場合に用いられる。 p変量(p次元)の観測値をm個(m次元)の主成分に縮約させるという意味で、 次元を減少させる(reduce)方法と言うこともでき、 多変量データを要約する有力な方法である。

 両者は似た目的に使われるが、元になっている考え方は異なるので 利用する場面では注意が必要である。

2.1. アイディア

 ここにあるクラスの学生の2科目の試験成績があったとしよう。以下はその散布図である。 このとき、「総合的な成績の指標」として2つの科目の「単純加算(合計点)」 が用いられることがよく行われている。 例えば、この合計点がある閾値を超えていた場合に進級出来る等の状況が想起される。 大学入試の一般選抜における合否ラインもこの考え方で用いられている。

 この「単純加算」とは散布図で考えるとどういう状況を示しているのであろうか? 容易に想像できるであろうが、「傾き1の直線」に各点からの垂線を下ろした時の 足の位置で学生の成績を代表させていることになる。

 つまり、この「傾き1の直線」を新たな軸\(z\)と定めて、その数直線上の位置で 学生の成績を表現することができた。 つまり、2次元の情報を1次元の情報で代表(次元を縮約)させることができた。

\[ z = x_{1} + x_{2} \]

 では、どのようなデータに対しても、この「傾き1の直線」を新たな軸\(z\)に用いることが 妥当であろうか? 学生の成績をより分解して表現したいと考えた場合、 「傾き1の直線」ではうまく表現できない状況もあるのではないか。 例えば、教科\(x_{1}\)の作題に失敗し、 \(x_{1}\)側の散らばり(分散)が小さいような場合の新たな軸\(z\)は どのように取るのが妥当と考えるか?

 つまり、データに応じてそれぞれの軸\(x_{j}\)に対する重み\(a_{j}\)を適切に取って 新たな軸\(z\)を定義することにより、次元を縮約した時により多くの情報を 表現する軸とすることができることを解るであろう。

\[ z = a_{1}x_{1} + a_{2}x_{2} \]

 では、この新たな軸\(z\)の決め方は統計の概念で考えると、どのように表現することになるのか。 それは、データの散らばりが最大になる軸を見つける作業に相当することが理解できるか? つまり、「分散最大」となる軸を見つけ出すことに相当する。

 あくまでも理解を促進するための模式的な説明になるが、 データ空間に「ノギス」を持ち込んで、その最大幅となる軸を見つけ出す行為を 想像すると理解し易いのではないか。ただし、データには密度があり、 これを考慮した上での分散であるので、「単なる端点同士」を 見るけることではないことに注意する必要がある。

[参考:立体の測定: イメージとして]

2.2. 前準備(Rでの利用に際して)

 回帰分析の時と同様に、この後のことを考えて、欠損値を除外した変数を予め作っておこう。 今回は、身長と体重のデータをStudHWとし、それに胸囲を加えたデータをStudHWCとする。 以下の出力から判る通り、元のデータは792サンプルだったが、 StudHWは691サンプル、StudHWCは242サンプルである。 変量が多くなると欠損値を含むサンプルが多くなり、 分析対象とするデータのサンプル数が減る傾向にある。

# データ切り出し
dim(Student24)
## [1] 792   9
StudHW<-na.omit(Student24[,3:4]) # 欠損値を含むデータを削除
dim(StudHW)
## [1] 691   2
# colnames(StudHW)
StudHW[1:5,]
##   Height Weight
## 1  145.0     38
## 2  145.5     42
## 3  146.7     41
## 4  148.0     42
## 5  148.0     43
StudHWC<-na.omit(Student24[,3:5]) # 欠損値を含むデータを削除
dim(StudHWC)
## [1] 242   3
# colnames(StudHWC)
StudHWC[1:5,]
##   Height Weight Chest
## 2  145.5     42    76
## 3  146.7     41    85
## 5  148.0     43    80
## 8  150.0     43    82
## 9  150.0     46    86

3. 2変量の場合の主成分分析 : 理解を助けるため

  1. 定式化 :

    • 変数 : \(x_{1}, x_{2}\)
    • 重み(係数) : \(a_{1}, a_{2}\)
    • 合成変量(線形結合) : \(z = a_{1}x_{1} + a_{2}x_{2}\)
    • 「よく代表する」ように、\(a_{1}, a_{2}\)を決める。
    • より広がって測定できる軸に沿うと情報量が多い。===> 広がり ===> 分散を最大に
      分散 \(V(z) \;\;\; { \to max} \;\;\;\;\; (制約条件 : a_{1}^2+a_{2}^2=1 の下で)\)


      なお、 \[ V(z) = \frac{1}{n} \sum_{i=1}^{n}{(z_{i}-\bar{z})^2} =a_{1}^2 s_{11} + 2 a_{1} a_{2} s_{12} + a_{2}^2 s_{23} \;。 \]

      ただし、\(x_{1}, x_{2}\) の分散を \(s_{11}, s_{23}\)、共分散を \(s_{12}\) とする。

    • 全測定値の分散を最大化する身長と体重(\(x_{1}, x_{2}\))の 総合指標を示す係数\(a_{1}, a_{2}\)を求めよう。
    • 求まったベクトルを「軸」とか「主成分」と呼ぶ。

  2. 出力結果 :
    • Importance of components:
      各軸の標準偏差(Standard deviation, 固有値の平方根)、 比率(Proportion of Variance, 寄与率, 全体に占める各軸の重み)、 累積寄与率(Cumulative Proportion, 前者の累積和)
    • 各軸ごとに
      • 固有ベクトル(係数\(a_{1}, a_{2}\)) : 変量の特長の解釈に使う
      • 主成分得点 : 各個人の得点(\(z\))
    • 身長と体重の散布図
    • 第1軸と第2軸の主成分得点(合成変数)の散布図(x軸が第1軸)

    • Rでは計算結果(PCAresultHW)に対して、
      • $rotation =各主成分軸の固有ベクトル
      • $center =解析に用いた各項目の値の平均
      • $sdev =標準偏差 (固有値の正の平方根)
      • $x =主成分得点
      等が格納される.

    # 主成分分析(2変量)
    PCAresultHW<-prcomp(StudHW, scale=F)  # 分散共分散行列に基づく
    
    # 分析結果要約
    summary(PCAresultHW)
    ## Importance of components:
    ##                            PC1    PC2
    ## Standard deviation     11.0494 4.8732
    ## Proportion of Variance  0.8372 0.1628
    ## Cumulative Proportion   0.8372 1.0000

    ◎主成分軸の固有ベクトル(各軸の特長を把握する)と回転の中心店の座標。

    PCAresultHW$rotation
    ##              PC1        PC2
    ## Height 0.5744243  0.8185577
    ## Weight 0.8185577 -0.5744243
    PCAresultHW$center
    ##    Height    Weight 
    ## 169.53965  59.94356

    ◎主成分得点(合成変数z)の先頭10ケースを表示させてみる。

    PCAresultHW$x[1:10,]
    ##          PC1       PC2
    ## 1  -32.05824 -7.482208
    ## 2  -28.49680 -9.370626
    ## 3  -28.62605 -7.813932
    ## 4  -27.06074 -7.324232
    ## 5  -26.24218 -7.898656
    ## 7  -24.03064 -8.228947
    ## 8  -25.09333 -6.261541
    ## 9  -22.63766 -7.984813
    ## 10 -21.81910 -8.559238
    ## 11 -25.33747 -4.868559

    ◎仕組みを理解してもらうために元データの散布図と、主成分得点(合成変数)のデータの散布図を描いてみる。

    summary(StudHW)
    ##      Height          Weight      
    ##  Min.   :145.0   Min.   : 35.00  
    ##  1st Qu.:165.0   1st Qu.: 53.75  
    ##  Median :170.0   Median : 60.00  
    ##  Mean   :169.5   Mean   : 59.94  
    ##  3rd Qu.:175.0   3rd Qu.: 65.00  
    ##  Max.   :187.0   Max.   :100.00
    plot(StudHW[,2],StudHW[,1])
    abline(v=seq(40,100,10), lty=3)
    abline(h=seq(150,180,10), lty=3)
    abline(h=169.4, lty=1)
    abline(v=59.8, lty=1)

    plot(PCAresultHW$x[,1], PCAresultHW$x[,2])
    abline(h=seq(-20,15,5), lty=3)
    abline(v=seq(-30,40,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    # 別表記(何も指定をしないと第1軸と第2軸を描画)
    plot(PCAresultHW$x)
    abline(v=seq(-30,400,10), lty=3)
    abline(h=seq(-20,20,5), lty=3)
    abline(v=0, lty=1)
    abline(h=0, lty=1)

    ◎ちょっと検算してみましょうか。求められた \(a_{1}, a_{2}\) を用いてそれぞれの軸の主成分得点が計算できていることを確認せよ。

    StudHW_PCA<-cbind(StudHW, PCAresultHW$x)  # 列側(横方向、column bind)に連結(結合)
    StudHW_PCA[1:10,]
    ##    Height Weight       PC1       PC2
    ## 1   145.0     38 -32.05824 -7.482208
    ## 2   145.5     42 -28.49680 -9.370626
    ## 3   146.7     41 -28.62605 -7.813932
    ## 4   148.0     42 -27.06074 -7.324232
    ## 5   148.0     43 -26.24218 -7.898656
    ## 7   149.0     45 -24.03064 -8.228947
    ## 8   150.0     43 -25.09333 -6.261541
    ## 9   150.0     46 -22.63766 -7.984813
    ## 10  150.0     47 -21.81910 -8.559238
    ## 11  151.0     42 -25.33747 -4.868559
    PCAresultHW$center                      # 回転の中心点
    ##    Height    Weight 
    ## 169.53965  59.94356
    PCAresultHW$rotation                    # 各軸の重み
    ##              PC1        PC2
    ## Height 0.5744243  0.8185577
    ## Weight 0.8185577 -0.5744243
    例えば、1番目のサンプル(145.0, 38.0)の主成分得点は、
    第1軸(PC1): -31.6 -32.1=(145.0-169.5)x0.574+(38.0-59.4 59.9)x0.819。
    第2軸(PC2): -7.8 -7.5=(145.0-169.5)x0.819+(38.0-59.4 59.9)x(-0.574)。

    ◎散布図でも確かめてみよう。理解を促進するために本来の丸マークではなく、データ番号で散布図を描く。 元データは身長順、体重順に並べ替えてあり、主成分得点の小さいものが第1軸の左方向に位置していることが理解できる。 回帰分析時にも話題になった体重が重い7-8サンプルが右下に位置している。

    plot(StudHW_PCA[,3], StudHW_PCA[,4], type="n") # マークを描かない
    text(StudHW_PCA[,3], StudHW_PCA[,4])           # その場所に番号を表示
    abline(h=seq(-20,20,5), lty=3)
    abline(v=seq(-30,40,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

  3. 解釈方法
    • データサイズ: 691名、2変量(身長、体重)
    • 寄与率(Proportion of Variance) : その軸がどの程度説明力を持っているか :
      • 第1軸だけで十分(83.7%)。第2軸に含まれる説明力は小さい(16.3%)。
    • 固有ベクトル($rotation) : その軸の特徴を示している :
      • 身長と体重の重みはほぼ同等(0.57と0.82)だが、体重がやや大きめに効いている(第1軸)
    • 主成分得点と散布図 : 各個人がどこに付置されているか
      • 第1軸 : 全体的な体格の指標。身長と体重を足したような指標。
      • ちなみに第2軸 : 全体的な体格の指標だが、体重が負に効く。
    • 元の散布図が回転していることが判るであろうか?
      • 回転の中心点はそれぞれの変量の平均値。

  4. 【テクニック】
    • 事前に欠損値を除外したデータを生成(抽出)し、 そのデータを格納した変数に対して統計手法(回帰分析、主成分分析等)を適用する。

4. 3変量以上の主成分分析


  1. 定式化 :
    • 2変量の時の考え方を拡張したもの(p変量)
    • 合成変量(線形結合) : \(z=\sum_{j=1}^{p}{a_{j} x_{j}}\)
    • より広がって測定できる軸に沿うと情報量が多い。===> 広がり ===> 分散を最大に

      分散 \(V(z) \;\;\; { \to max} \;\;\;\;\; (制約条件 : \sum_{j=1}^{p}{a_{j}^2}=1 の下で)\)

    • 合成変量の分散を最大化する軸を決定する。

    • よく代表するように、\({a_{1},a_{2},...}\) を決める。

  2. 身長、体重、胸囲の3変量での総合指標を求めてみよう :

  3. 出力結果 :
    • Importance of components:
      各軸の標準偏差(Standard deviation, 固有値の平方根)、 比率(Proportion of Variance, 寄与率, 全体に占める各軸の重み)、 累積寄与率(Cumulative Proportion, 前者の累積和)
    • 各軸ごとに
      • 固有ベクトル(係数\(a_{1}, a_{2},...\)) : 変量の特長の解釈に使う
      • 主成分得点 : 各個人の得点(\(z\))
    • 第1軸~第3軸の散布図
    • 各軸の主成分得点(合成変数)

    • Rでは計算結果に対して、
      • $rotation =各主成分軸の固有ベクトル
      • $center =解析に用いた各項目の値の平均
      • $sdev =標準偏差 (固有値の正の平方根)
      • $x =主成分得点
      等が格納される.

    # 主成分分析(3変量)
    PCAresultHWC<-prcomp(StudHWC, scale=F)  # 分散共分散行列に基づく
    
    # 分析結果要約
    summary(PCAresultHWC)
    ## Importance of components:
    ##                            PC1    PC2    PC3
    ## Standard deviation     13.2074 7.9230 4.7961
    ## Proportion of Variance  0.6704 0.2412 0.0884
    ## Cumulative Proportion   0.6704 0.9116 1.0000
    PCAresultHWC$rotation
    ##              PC1        PC2        PC3
    ## Height 0.4535489 -0.4041618  0.7943215
    ## Weight 0.7199210 -0.3592390 -0.5938528
    ## Chest  0.5253639  0.8411900  0.1280321
    PCAresultHWC$center
    ##    Height    Weight     Chest 
    ## 169.15083  61.06983  86.08347
    PCAresultHWC$x[1:10,]
    ##          PC1        PC2        PC3
    ## 2  -29.75307   7.927275  -8.752693
    ## 3  -25.20046  15.372230  -6.053366
    ## 5  -25.79782   9.922391  -6.848614
    ## 8  -23.84000  10.796448  -5.003907
    ## 9  -19.57878  13.083490  -6.273337
    ## 16 -25.19957   8.965851  -3.018845
    ## 17 -31.31909   8.656086   0.695398
    ## 23 -17.33281  12.532575  -4.059267
    ## 25 -15.94175   1.908334 -10.259304
    ## 26 -35.51647 -32.457743 -16.568294
    plot(PCAresultHWC$x[,1], PCAresultHWC$x[,2])
    abline(h=seq(-40,40,10), lty=3)
    abline(v=seq(-30,40,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC$x[,2], PCAresultHWC$x[,3])
    abline(h=seq(-15,15,5), lty=3)
    abline(v=seq(-40,40,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC$x[,3], PCAresultHWC$x[,1])
    abline(h=seq(-20,40,20), lty=3)
    abline(v=seq(-15,10,5), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    ◎理解を促進するために本来の丸マークではなく、データ番号で散布図を描く。

    # 番号で描画してみる。
    plot(PCAresultHWC$x[,1], PCAresultHWC$x[,2], type="n") # マークを描かない
    text(PCAresultHWC$x[,1], PCAresultHWC$x[,2])           # その場所に番号を表示
    abline(h=seq(-40,40,10), lty=3)
    abline(v=seq(-30,40,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC$x[,2], PCAresultHWC$x[,3], type="n") # マークを描かない
    text(PCAresultHWC$x[,2], PCAresultHWC$x[,3])           # その場所に番号を表示
    abline(h=seq(-15,10,5), lty=3)
    abline(v=seq(-40,10,10), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC$x[,3], PCAresultHWC$x[,1], type="n") # マークを描かない
    text(PCAresultHWC$x[,3], PCAresultHWC$x[,1])           # その場所に番号を表示
    abline(h=seq(-20,40,20), lty=3)
    abline(v=seq(-15,10,5), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

  4. 解釈方法
    • データサイズ: 242名、3変量(身長、体重、胸囲)
    • 寄与率 : 2軸まで取れば十分のようだ(累積の項 91.1%)。
    • 各軸の意味付け
      • 第1軸 : 全体的な体格の因子。やや体重が多めに効いている。
      • 第2軸 : 太さの因子(?)。胸囲がで大きい、身長・体重は

      • ちなみに第3軸 : 華奢さの因子(?)。無視しても良い軸であるが。(8.8%)。

5. 相関行列を使う理由

  1. 相関行列を用いて体格の総合指標を求めてみよう :

  2. 出力結果 :
    • Importance of components:
      各軸の標準偏差(Standard deviation, 固有値の平方根)、 比率(Proportion of Variance, 寄与率, 全体に占める各軸の重み)、 累積寄与率(Cumulative Proportion, 前者の累積和)
    • 各軸ごとに
      • 固有ベクトル(係数\(a_{1}, a_{2},...\)) : 変量の特長の解釈に使う
      • 主成分得点 : 各個人の得点(\(z\))
    • 第1軸~第3軸の散布図
    • 各軸の主成分得点(合成変数)

    • Rでは計算結果に対して、
      • $rotation =各主成分軸の固有ベクトル
      • $center =解析に用いた各項目の値の平均
      • $sdev =標準偏差 (固有値の正の平方根)
      • $x =主成分得点
      等が格納される.

    # 主成分分析(3変量, 相関行列を用いて)
    PCAresultHWC2<-prcomp(StudHWC, scale=T)  # 相関行列に基づく(scale=Tがそのことを示す)
    
    # 分析結果要約
    summary(PCAresultHWC2)
    ## Importance of components:
    ##                           PC1    PC2     PC3
    ## Standard deviation     1.4071 0.8546 0.53837
    ## Proportion of Variance 0.6599 0.2434 0.09661
    ## Cumulative Proportion  0.6599 0.9034 1.00000
    PCAresultHWC2$rotation
    ##              PC1        PC2        PC3
    ## Height 0.5973048  0.4933666 -0.6323104
    ## Weight 0.6417864  0.1787820  0.7457527
    ## Chest  0.4809752 -0.8512499 -0.2098486
    PCAresultHWC2$center
    ##    Height    Weight     Chest 
    ## 169.15083  61.06983  86.08347
    # PCAresultHWC2$x[1:10,]
    
    # 別表記
    # plot(PCAresultHWC2$x)
    # abline(h=seq(-5,1,1), lty=3)
    # abline(v=seq(-4,4,1), lty=3)
    # abline(h=0, lty=1)
    # abline(v=0, lty=1)
    
    plot(PCAresultHWC2$x[,1], PCAresultHWC2$x[,2])
    abline(h=seq(-2,5,1), lty=3)
    abline(v=seq(-4,4,1), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC2$x[,2], PCAresultHWC2$x[,3])
    abline(h=seq(-1,2,0.5), lty=3)
    abline(v=seq(-2,5,1), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

    plot(PCAresultHWC2$x[,3], PCAresultHWC2$x[,1])
    abline(h=seq(-4,4,2), lty=3)
    abline(v=seq(-1,2,0.5), lty=3)
    abline(h=0, lty=1)
    abline(v=0, lty=1)

  3. 解釈方法
    • データサイズ: 242名、3変量(身長、体重、胸囲)
    • 寄与率 : 2軸まで取れば十分のようだ(累積の項 90.3%)。
    • 各軸の意味付け
      • 第1軸 : 全体的な体格の因子。
      • 第2軸 : 細さの因子。
    • 分散共分散行列を使ったときよりも第1軸の固有ベクトル同士が近い値になった。しかし、軸の解釈に違いはない。その理由はこの例では 3変量のスケールや分散にそれほどの違いがないためと想像される。

[演習5-1] 分散共分散行列と相関行列を使ったときの違いを見てみよう。身長のみを mm 単位で測定したと考えて、身長のみ10倍したデータを生成後、両者(分散共分散行列と相関行列)の出力を比較してみよ。

6. 主成分の数の決定基準

これでなければならないというようなルールが定まっているわけではないが、以下のような基準が一般的に 用いられている。また、結果の解釈の都合上、多少増減させることもある。

[ノウハウ] 軸の特性を把握するにはラインマーカーが重宝する。 一つの軸内で大きい(もしくは小さい)固有ベクトルにマークすると、 何がグルーピングされているかが理解し易くなる。 ただ、ここまでの例示(学生体格データ)では3軸までしか出現しないので、ご利益は感じにくいとも思う。 多変量(例えば以下の演習)になればその有効性が理解できるであろう。

7. [例題PCA] 食品の嗜好性を探ってみよう

100種類の食品(ごはん、お茶漬け、…、リンゴ、パイ缶)に対する 性・年齢で分割した10グループの嗜好度調査のデータをMoodleに掲載しておく。 グループ1から5は男性、グループ6から10は女性であり、その中の番号の小さい方から順に 15才以下、16~20才以下、21~30才以下、31~40才以下、41才以上の10群を構成している。 また、測定尺度は「1: おらく食べる気がしない」、「2: もし強制されれば食べる」、…、 「8: いつも食べたい」、「9: もっとも好きな食品に入る」までの9段階であり、 各グループごとに尺度の平均値を取ったものが測定値として格納されている。

なお、データのファイル名は「food.xlsx」(エクセル形式)であり、先頭17行には説明が、それに続く行にヘッダーとしてグループ名が格納されている。 分析に供する際にはcsv形式に変換すると共にスキップする必要がある。 縦方向に100種類の食品が、横方向にグループが配置されている。

# ディレクトリの移動。必須ではない。個々人の設定に応じて。
setwd("D:/home_sub3/R_Dir")  # ホームディレクトリに移動(Set Working Directory)
getwd()           # 現在のディレクトリ位置を表示
## [1] "D:/home_sub3/R_Dir"
list.files()      # ファイル・ディレクトリ一覧を表示
##  [1] "DNC_Tsuiseki"           "DNC21"                  "Dragons"               
##  [4] "Food"                   "grain"                  "JEES"                  
##  [7] "KougiDS20"              "KougiDS21"              "KougiDS22"             
## [10] "KougiDS23"              "KougiDS24"              "LibraryInstall2211a.R" 
## [13] "Misc"                   "NitechNSK21_R"          "NitechNSK22_R"         
## [16] "NitechNSK23_R"          "plot1.png"              "PresentationSample.Rmd"
## [19] "PresenTest.html"        "PresenTest.Rmd"         "R_Dir_NKK20"           
## [22] "StatEdu24"              "StatM20"                "StatM21"               
## [25] "StatM22"                "StatM23"                "Terao.zip"             
## [28] "Terao_Lenovo"           "Terao1"                 "Terao2"                
## [31] "Unemployment"           "確立楕円.R"
setwd("KougiDS24")  # ディレクトリを移動
list.files()      # ファイル・ディレクトリ一覧を表示
##  [1] "DNC21a_suugaku1.pdf"           "DNC21a_suugaku1_04.pdf"       
##  [3] "DNC21a_suugaku1_a.pdf"         "DNC21a_suugaku1_a_02.pdf"     
##  [5] "Dragons24.csv"                 "DS01_Questionnaire1.csv"      
##  [7] "DS01_Questionnaire2.csv"       "DS01_Questionnaire3.csv"      
##  [9] "DS2402_1.pdf"                  "DS2402_1c.html"               
## [11] "DS2402_1c.Rmd"                 "DS2402_2.pdf"                 
## [13] "DS2402_2a.html"                "DS2402_2a.Rmd"                
## [15] "DS2402_3.pdf"                  "DS2402_3b.html"               
## [17] "DS2402_3b.Rmd"                 "DS2403_1.pdf"                 
## [19] "DS2403_1a.html"                "DS2403_1a.Rmd"                
## [21] "DS2404_1b.html"                "DS2404_1b.Rmd"                
## [23] "DS2404_1c3a_Ans.html"          "DS2404_1c3a_Ans.Rmd"          
## [25] "DS2404_3b.html"                "DS2404_3b.Rmd"                
## [27] "DS2405_1a.html"                "DS2405_1a.Rmd"                
## [29] "DS2405_7a.html"                "DS2406_7b.html"               
## [31] "DS2406_7b.Rmd"                 "DS2406_7c.html"               
## [33] "DS2406_7c.Rmd"                 "DS2406_RegTopics24a_Haifu.pdf"
## [35] "DS2407_1c.html"                "DS2407_1c.Rmd"                
## [37] "DS2407_1d.BAK"                 "DS2407_1d.html"               
## [39] "DS2407_1d.Rmd"                 "DS2407_1d_files"              
## [41] "DS2407_8a.html"                "DS2407_8a.Rmd"                
## [43] "ExcelStudent24b.jpg"           "Exp"                          
## [45] "food.csv"                      "food.xlsx"                    
## [47] "K2023_Fig2.jpg"                "lec01"                        
## [49] "NewsPaper"                     "OldFiles"                     
## [51] "ReadMe_food.pdf"               "RStudio_Display.jpg"          
## [53] "RStudio_Icon.jpg"              "Seminar24a.pdf"               
## [55] "Seminar24b.pdf"                "StockFiles"                   
## [57] "StudAll23b.csv"                "StudAll24b.csv"               
## [59] "新聞記事"                      "労働男女集約1_抽出2EUC.csv"
# データの読み込み
Food<-read.csv("food.csv", skip=17, header=TRUE)
dim(Food)
## [1] 100  10
colnames(Food)
##  [1] "M1" "M2" "M3" "M4" "M5" "F1" "F2" "F3" "F4" "F5"
Food[1:5,]
##     M1   M2   M3   M4   M5   F1   F2   F3   F4   F5
## 1 7.69 7.31 7.47 7.76 7.87 7.51 7.24 7.70 7.91 7.95
## 2 6.59 5.56 6.21 6.04 5.81 6.64 6.11 6.53 6.44 6.64
## 3 4.55 4.18 4.36 4.25 4.53 4.60 3.66 4.04 3.68 4.43
## 4 6.78 6.11 6.30 5.98 5.56 6.37 6.29 5.43 5.32 5.28
## 5 6.47 6.24 6.02 5.42 5.88 6.00 5.60 4.60 5.40 5.95
# 主成分分析(10変量, 相関行列を用いて)
PCAresultFood<-prcomp(Food, scale=T)  # 相関行列に基づく(scale=Tがそのことを示す)

# 分析結果要約
summary(PCAresultFood)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5    PC6     PC7
## Standard deviation     2.6130 1.3274 0.86859 0.51223 0.34864 0.3130 0.26851
## Proportion of Variance 0.6828 0.1762 0.07545 0.02624 0.01216 0.0098 0.00721
## Cumulative Proportion  0.6828 0.8590 0.93443 0.96067 0.97282 0.9826 0.98983
##                            PC8     PC9    PC10
## Standard deviation     0.20995 0.18908 0.14797
## Proportion of Variance 0.00441 0.00358 0.00219
## Cumulative Proportion  0.99424 0.99781 1.00000
PCAresultFood$rotation
##           PC1         PC2         PC3         PC4         PC5         PC6
## M1 -0.2860327  0.44633451 -0.19351244 -0.42801884  0.16236540  0.01641257
## M2 -0.3313365  0.23984218 -0.33606342 -0.02248782 -0.55959447  0.21236690
## M3 -0.3233449 -0.16633727 -0.44229066  0.43602904 -0.16859426 -0.47692872
## M4 -0.2993294 -0.35862650 -0.37536562 -0.06344873  0.36791199  0.56249120
## M5 -0.2607266 -0.50720850 -0.12741854 -0.37542550  0.14687921 -0.38512320
## F1 -0.3086352  0.40788210  0.08369528 -0.26737495  0.28686604 -0.20987842
## F2 -0.3442715  0.25269723  0.17139988  0.29565475 -0.02504968 -0.13746884
## F3 -0.3478768  0.03231395  0.28908713  0.50750798  0.45237713  0.12839022
## F4 -0.3456363 -0.16436803  0.32223577 -0.04001194 -0.38894366  0.38718897
## F5 -0.3033345 -0.26727324  0.52255945 -0.25126957 -0.19050693 -0.18195525
##            PC7         PC8         PC9        PC10
## M1  0.06213786  0.03849263  0.14161684  0.66805207
## M2 -0.47946501 -0.28332452  0.01373924 -0.22506354
## M3  0.41635370 -0.13615016 -0.08592164  0.16395971
## M4  0.06624504  0.11430118 -0.40371304 -0.06834407
## M5 -0.32530957  0.16753406  0.43783314 -0.14864786
## F1  0.33505785 -0.17613659 -0.09053845 -0.61810688
## F2 -0.23610392  0.76265391 -0.20438194 -0.04635090
## F3 -0.25613467 -0.38298261  0.30326953  0.10686285
## F4  0.48882147  0.16197363  0.42518761 -0.03038089
## F5 -0.10063221 -0.27018517 -0.54310699  0.22990381
PCAresultFood$center
##     M1     M2     M3     M4     M5     F1     F2     F3     F4     F5 
## 6.0381 5.7848 5.9471 5.6695 5.6406 5.7813 5.5639 5.3794 5.5174 5.5421
# PCAresultFood$x[,1]

plot(PCAresultFood$x[,1], PCAresultFood$x[,2], type="n")
text(PCAresultFood$x[,1], PCAresultFood$x[,2])
abline(h=seq(-2,3,1), lty=3)
abline(v=seq(-6,7,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultFood$x[,2], PCAresultFood$x[,3], type="n")
text(PCAresultFood$x[,2], PCAresultFood$x[,3])
abline(h=seq(-2,1,1), lty=3)
abline(v=seq(-2,3,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultFood$x[,3], PCAresultFood$x[,1], type="n")
text(PCAresultFood$x[,3], PCAresultFood$x[,1])
abline(h=seq(-6,6,2), lty=3)
abline(v=seq(-2,1,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)


81. 参考

このページで取り扱ったプログラムだけを抜き出して以下に列挙しておく。


# ディレクトリの移動。必須ではない。個々人の設定に応じて。
setwd("D:/home_sub3/R_Dir")  # ホームディレクトリに移動(Set Working Directory)
getwd()           # 現在のディレクトリ位置を表示
list.files()      # ファイル・ディレクトリ一覧を表示
setwd("KougiDS24")  # ディレクトリを移動
list.files()      # ファイル・ディレクトリ一覧を表示

# データの読み込み
Student24<-read.csv("StudAll24b.csv", skip=5,
                     header=TRUE, na.strings="NA")
dim(Student24)
colnames(Student24)
Student24[1:5,]

# データ切り出し
dim(Student24)
StudHW<-na.omit(Student24[,3:4]) # 欠損値を含むデータを削除
dim(StudHW)
# colnames(StudHW)
StudHW[1:5,]

StudHWC<-na.omit(Student24[,3:5]) # 欠損値を含むデータを削除
dim(StudHWC)
# colnames(StudHWC)
StudHWC[1:5,]


# 主成分分析(2変量)
PCAresultHW<-prcomp(StudHW, scale=F)  # 分散共分散行列に基づく

# 分析結果要約
summary(PCAresultHW)

PCAresultHW$rotation
PCAresultHW$center

PCAresultHW$x[1:10,]

summary(StudHW)

plot(StudHW[,2],StudHW[,1])
abline(v=seq(40,100,10), lty=3)
abline(h=seq(150,180,10), lty=3)
abline(h=169.4, lty=1)
abline(v=59.8, lty=1)

plot(PCAresultHW$x[,1], PCAresultHW$x[,2])
abline(h=seq(-20,10,5), lty=3)
abline(v=seq(-30,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

# 別表記(何も指定をしないと第1軸と第2軸を描画)
plot(PCAresultHW$x)
abline(v=seq(-30,400,10), lty=3)
abline(h=seq(-20,10,5), lty=3)
abline(v=0, lty=1)
abline(h=0, lty=1)

StudHW_PCA<-cbind(StudHW, PCAresultHW$x)  # 列側(横方向、column bind)に連結(結合)
StudHW_PCA[1:10,]
PCAresultHW$center                      # 回転の中心点
PCAresultHW$rotation                    # 各軸の重み

plot(StudHW_PCA[,3], StudHW_PCA[,4], type="n") # マークを描かない
text(StudHW_PCA[,3], StudHW_PCA[,4])           # その場所に番号を表示
abline(h=seq(-20,10,5), lty=3)
abline(v=seq(-30,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)


# 主成分分析(3変量)
PCAresultHWC<-prcomp(StudHWC, scale=F)  # 分散共分散行列に基づく

# 分析結果要約
summary(PCAresultHWC)

PCAresultHWC$rotation

PCAresultHWC$center

PCAresultHWC$x[1:10,]

# 別表記
# plot(PCAresultHWC$x)
# abline(v=seq(-20,40,10), lty=3)
# abline(h=seq(-50,10,10), lty=3)
# abline(v=0, lty=1)
# abline(h=0, lty=1)

plot(PCAresultHWC$x[,1], PCAresultHWC$x[,2])
abline(h=seq(-10,40,10), lty=3)
abline(v=seq(-30,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC$x[,2], PCAresultHWC$x[,3])
abline(h=seq(-15,10,5), lty=3)
abline(v=seq(-10,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC$x[,3], PCAresultHWC$x[,1])
abline(h=seq(-20,40,20), lty=3)
abline(v=seq(-15,10,5), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

# 番号で描画してみる。
plot(PCAresultHWC$x[,1], PCAresultHWC$x[,2], type="n") # マークを描かない
text(PCAresultHWC$x[,1], PCAresultHWC$x[,2])           # その場所に番号を表示
abline(h=seq(-10,40,10), lty=3)
abline(v=seq(-30,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC$x[,2], PCAresultHWC$x[,3], type="n") # マークを描かない
text(PCAresultHWC$x[,2], PCAresultHWC$x[,3])           # その場所に番号を表示
abline(h=seq(-15,10,5), lty=3)
abline(v=seq(-10,40,10), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC$x[,3], PCAresultHWC$x[,1], type="n") # マークを描かない
text(PCAresultHWC$x[,3], PCAresultHWC$x[,1])           # その場所に番号を表示
abline(h=seq(-20,40,20), lty=3)
abline(v=seq(-15,10,5), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)


# 主成分分析(3変量, 相関行列を用いて)
PCAresultHWC2<-prcomp(StudHWC, scale=T)  # 相関行列に基づく(scale=Tがそのことを示す)

# 分析結果要約
summary(PCAresultHWC2)

PCAresultHWC2$rotation
PCAresultHWC2$center

# PCAresultHWC2$x[1:10,]

# 別表記
# plot(PCAresultHWC2$x)
# abline(h=seq(-5,1,1), lty=3)
# abline(v=seq(-4,4,1), lty=3)
# abline(h=0, lty=1)
# abline(v=0, lty=1)

plot(PCAresultHWC2$x[,1], PCAresultHWC2$x[,2])
abline(h=seq(-5,1,1), lty=3)
abline(v=seq(-3,4,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC2$x[,2], PCAresultHWC2$x[,3])
abline(h=seq(-1,2,0.5), lty=3)
abline(v=seq(-5,1,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultHWC2$x[,3], PCAresultHWC2$x[,1])
abline(h=seq(-4,4,2), lty=3)
abline(v=seq(-1,2,0.5), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

# ディレクトリの移動。必須ではない。個々人の設定に応じて。
setwd("D:/home_sub3/R_Dir")  # ホームディレクトリに移動(Set Working Directory)
getwd()           # 現在のディレクトリ位置を表示

list.files()      # ファイル・ディレクトリ一覧を表示

setwd("KougiDS24")  # ディレクトリを移動
list.files()      # ファイル・ディレクトリ一覧を表示

# データの読み込み
Food<-read.csv("food.csv", skip=17, header=TRUE)
dim(Food)
colnames(Food)
Food[1:5,]

# 主成分分析(10変量, 相関行列を用いて)
PCAresultFood<-prcomp(Food, scale=T)  # 相関行列に基づく(scale=Tがそのことを示す)

# 分析結果要約
summary(PCAresultFood)

PCAresultFood$rotation

PCAresultFood$center

plot(PCAresultFood$x[,1], PCAresultFood$x[,2], type="n")
text(PCAresultFood$x[,1], PCAresultFood$x[,2])
abline(h=seq(-2,3,1), lty=3)
abline(v=seq(-6,7,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultFood$x[,2], PCAresultFood$x[,3], type="n")
text(PCAresultFood$x[,2], PCAresultFood$x[,3])
abline(h=seq(-2,1,1), lty=3)
abline(v=seq(-2,3,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)

plot(PCAresultFood$x[,3], PCAresultFood$x[,1], type="n")
text(PCAresultFood$x[,3], PCAresultFood$x[,1])
abline(h=seq(-6,6,2), lty=3)
abline(v=seq(-2,1,1), lty=3)
abline(h=0, lty=1)
abline(v=0, lty=1)