射的データでみる雲の分散と分散の関係VbetweenとVwithin

/home/kazuo/RStudio/MCA-chap5/bwtween_within.Rmd.rmarkdown

Author
Affiliation
kazuo fujimoto

津田塾大学 数学・計算機科学研究所

Published

May 6, 2023

1 概要 第二章で理解すべきポイント

本稿は、『多重対応分析』の第二章の空間の寄与、寄与率の計算について整理したものであるが、この内容は、第３章以降の内容を理解するために必要なものである。

MCAとは、多変数表を多カテゴリ表、指標行列へと変換してそれに対するCAとして理解される。そこでは、その指標行列の残差行列Sに対するSVD（特異値分解）が処理の中心であり、それによって得られるものは、SVDされた三つの行列、U、V、そしてDである。

S = U D V^{t}

これらは、Dが変換された空間の座標軸に対応し、UとVはそれぞれ、行座標、列座標に関係する。そして、この三つの行列から得られるものが出発点であり、MCAにまつわる、様々な統計量は、この生成された空間座標に対する演算から取得されるものである。

この第２章でまとめられている平均点座標の計算、慣性（分散）の計算、そして寄与率の計算は、このためにある。

つまり、ここが理解できれば、基本的なMCAの処理コードをかけるようになる。もちろん、すべてを自前で書く必要はなく、GDAtools、FactoMineRといったパッケージが、CAやMCAに関連する諸統計量を算出し、図示するfunctionを提供してくれている。実際の分析においてはそれらを用いることになるが、そこで、実行されている処理を理解すれば、応用的な展開も可能になる。

こうした視点で見た場合、第２章は、平面（x1軸、x2軸）をまとめた値として計算しているが、本文にもあるように、実際の分析、解釈では軸ごとに行うことになるので、本稿でも、軸ごとに計算し、そののちに平面での総合を行っている。

出発点としているデータ表は、Target Example（射的データ）である。このデータは、的につくられた弾痕の分布をもとに計算している。この弾痕データは、MCAのリザルトに対応している。

また、三つの部分集合（A、B、C）に分割しているが、この分割ラベルは、追加変数のカテゴリに対応している。

MCAを分散の分解として理解する際に、この追加変数による部分集合への分割は、重要な要素になるが、他方、Rの処理においては、こうした処理を行うfunctionとして、dplyr系を理解しておくとよい。Excel　での処理と比べたときに、Rでdplyrを用いた計算のメリットは、処理の流れをスクリプトととして表現できるところにある。処理の流れは、パイプ処理で表現され、さらに、group_by によるグループごとに集約された計算が可能になっている。この第２章の計算をdplyr系のfunctionで表現することで、全体の流れはよりはっきりするはずである。

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

2 構成

  • 分析対象データを用意する。分割グループ情報も加えておく。
  • グラフを描画するのに必要な、全体の平均点座標、グループごとの平均点の座標を計算して付加しておく。
  • 全体、群内、群間の分散を計算する。
  • 寄与率を計算する。
  • 差の寄与率を計算する。

3 射的データ（Target Example）を用意する（分割名を加えてある：GP）

target <- tribble(~point,~x1,~x2,~GP,
                  #-----------------
                  "M1",0,-12,"A",
                  "M2",6,-10,"A",
                  "M3",14,-6,"C",
                  "M4",6,-2,"C",
                  "M5",12,0,"C",
                  "M6",-8,2,"B",
                  "M7",2,4,"C",
                  "M8",6,4,"C",
                  "M9",10,10,"C",
                  "M10",12,10,"C")
target
# A tibble: 10 × 4
   point    x1    x2 GP   
   <chr> <dbl> <dbl> <chr>
 1 M1        0   -12 A    
 2 M2        6   -10 A    
 3 M3       14    -6 C    
 4 M4        6    -2 C    
 5 M5       12     0 C    
 6 M6       -8     2 B    
 7 M7        2     4 C    
 8 M8        6     4 C    
 9 M9       10    10 C    
10 M10      12    10 C

3.1 分割グループの度数（wtを加えて）、x1,x2の平均点座標、群ごとの平均点座標を加える

target %>%  mutate(x1ave=mean(x1),x2ave=mean(x2)) %>%
  group_by(GP) %>% mutate(wt = n()) %>% 
                   mutate(x1Gave=mean(x1),x2Gave=mean(x2)) %>% 
  ungroup() -> .d
.d # 基本データ表
# A tibble: 10 × 9
   point    x1    x2 GP    x1ave x2ave    wt x1Gave x2Gave
   <chr> <dbl> <dbl> <chr> <dbl> <dbl> <int>  <dbl>  <dbl>
 1 M1        0   -12 A         6     0     2   3    -11   
 2 M2        6   -10 A         6     0     2   3    -11   
 3 M3       14    -6 C         6     0     7   8.86   2.86
 4 M4        6    -2 C         6     0     7   8.86   2.86
 5 M5       12     0 C         6     0     7   8.86   2.86
 6 M6       -8     2 B         6     0     1  -8      2   
 7 M7        2     4 C         6     0     7   8.86   2.86
 8 M8        6     4 C         6     0     7   8.86   2.86
 9 M9       10    10 C         6     0     7   8.86   2.86
10 M10      12    10 C         6     0     7   8.86   2.86

3.2 射的データの基本配置を図示する（pbase）

library(ggrepel)
.d %>% ggplot(aes(x=x1,y=x2)) + geom_point() + coord_fixed(ratio = 1) + geom_vline(xintercept = 0) + geom_hline(yintercept = 0) -> p
p <- p + geom_text_repel(aes(x=x1,y=x2,label=point)) 
p <- p + geom_point(aes(x=x1ave,y=x2ave)) + annotate("text",x=.d$x1ave,y=.d$x2ave+0.6,label="G") 
p <- p + geom_point(aes(x=x1Gave ,y= x2Gave)) + geom_text(aes(x=x1Gave ,y= x2Gave+0.6,label=GP))
pbase <- p

4 雲の分散を計算する

4.1 各点（M1〜M10）の平均点G（6.0）からの平方距離/まず、グラフ

4.2 Vtotal 雲全体の分散

p <- pbase + geom_segment(data = .d, aes(x = x1, y = x2, xend = x1ave, yend = x2ave),
               arrow = arrow(length = unit(0.1, "cm"), type = "closed"))
p

4.3 値を計算する

.d %>% transmute(var1=(x1 - x1ave)^2,var2= (x2 - x2ave)^2,dist=var1+var2) # p29 各点（M1〜M10）から平均点（6,0）までの平方距離
# A tibble: 10 × 3
    var1  var2  dist
   <dbl> <dbl> <dbl>
 1    36   144   180
 2     0   100   100
 3    64    36   100
 4     0     4     4
 5    36     0    36
 6   196     4   200
 7    16    16    32
 8     0    16    16
 9    16   100   116
10    36   100   136

4.4 軸ごと（x1、x2）の分散

.d %>% transmute(var1=(x1 - x1ave)^2,var2= (x2 - x2ave)^2) %>% summarize(sumvar1=sum(var1),sumvar2=sum(var2))
# A tibble: 1 × 2
  sumvar1 sumvar2
    <dbl>   <dbl>
1     400     520
.d %>% transmute(var1=(x1 - x1ave)^2,var2= (x2 - x2ave)^2) %>% summarize(sumvar1=sum(var1)/nrow(.),sumvar2=sum(var2)/nrow(.)) # x1 x2 の分散
# A tibble: 1 × 2
  sumvar1 sumvar2
    <dbl>   <dbl>
1      40      52

4.5 共分散

.d %>% transmute(cvar1=(x1 - x1ave)*(x2 - x2ave)) %>% summarize(cumvar=sum(cvar1)/nrow(.)) #　x1,x2の共分散
# A tibble: 1 × 1
  cumvar
   <dbl>
1      8

4.6 まとめて集計

.d %>% transmute(var1=(x1 - x1ave)^2,var2= (x2 - x2ave)^2) %>% summarize(sumvar1=sum(var1)/nrow(.),sumvar2=sum(var2)/nrow(.)) %>% mutate(varall=sumvar1+sumvar2,std_var=sqrt(varall))
# A tibble: 1 × 4
  sumvar1 sumvar2 varall std_var
    <dbl>   <dbl>  <dbl>   <dbl>
1      40      52     92    9.59

5 群化した部分集合の分散を計算する

5.1 群内分散（Groupごとの群平均と群内ポイントの分散合計）

p <- pbase + geom_segment(data = .d, aes(x = x1, y = x2, xend = x1Gave, yend = x2Gave,group=GP,colour=GP),
               arrow = arrow(length = unit(0.1, "cm"), type = "closed"))
p

5.2 Vwithin：群内分散の計算

群内の各点（M1〜M10）とそれぞれの群平均（A,B,C）の平方距離の総和

.d %>% group_by(GP) %>% transmute(aveM1=(x1-x1Gave)^2,abeM2=(x2-x2Gave)^2) %>% mutate(Vwithin0 = aveM1+abeM2) 
# A tibble: 10 × 4
# Groups:   GP [3]
   GP    aveM1 abeM2 Vwithin0
   <chr> <dbl> <dbl>    <dbl>
 1 A      9     1       10   
 2 A      9     1       10   
 3 C     26.4  78.4    105.  
 4 C      8.16 23.6     31.8 
 5 C      9.88  8.16    18.0 
 6 B      0     0        0   
 7 C     47.0   1.31    48.3 
 8 C      8.16  1.31     9.47
 9 C      1.31 51.0     52.3 
10 C      9.88 51.0     60.9 
.d %>% group_by(GP,wt) %>% transmute(aveM1=(x1-x1Gave)^2,abeM2=(x2-x2Gave)^2) %>% mutate(Vwithin0 = aveM1+abeM2) %>% 
  summarise(Vwithin_M=mean(Vwithin0)) %>% ungroup() %>% summarise(Vwithin=sum(wt*Vwithin_M)/10) %>% unlist -> Vwithin
`summarise()` has grouped output by 'GP'. You can override using the `.groups`
argument.
Vwithin
 Vwithin 
34.57143

5.3 Vbetween：群間分散の計算（群平均点のGから分散）

p <- pbase + geom_segment(data = .d, aes(x = x1Gave, y = x2Gave, xend = x1ave, yend = x2ave,group=GP,colour=GP),
               arrow = arrow(length = unit(0.1, "cm"), type = "closed"))
p

5.4 Vbetwenn：群間分散

各平均点（A,B,C）と重心Gの距離の重み付き平均。p31

p31の指揮では、GA2、GB2、GC^2をwtで重みをつけて平均しているが、データとしては、グループに集約されてないので、そのまま加算して要素数で割ればよい。つまり、平均すればよい。

NN <- nrow(.d)
.d %>%  group_by(GP) %>% transmute(x1G2=(x1ave-x1Gave)^2,x2G2=(x2ave-x2Gave)^2) %>% summarize(sumX1=sum(x1G2),sumX2=sum(x2G2)) %>% mutate(Vbetween0=sumX1+sumX2)
# A tibble: 3 × 4
  GP    sumX1 sumX2 Vbetween0
  <chr> <dbl> <dbl>     <dbl>
1 A      18   242        260 
2 B     196     4        200 
3 C      57.1  57.1      114.
.d %>%  group_by(GP) %>% transmute(x1G2=(x1ave-x1Gave)^2,x2G2=(x2ave-x2Gave)^2) %>% summarize(sumX1=sum(x1G2),sumX2=sum(x2G2)) %>% mutate(Vbetween0=sumX1+sumX2) %>% 
  summarise(Vbetween=sum(Vbetween0)/NN) %>% unlist -> Vbetween
Vbetween
Vbetween 
57.42857

5.5 検算 Vtotal = Vbetween + Vwithin

Vbetween + Vwithin
Vbetween 
      92

6 差の寄与率の検討

差の寄与率 “contribution of deviation”、という概念がある（Le Roux&Rouanet 1998 ）。 注目している部分集合（重み、平均点を持つ）が2つのとき、群間寄与（V_between）は、差の寄与率と等しくなるという。この差の寄与率、という考え方は、第３章にある、MCAが生成した変数空間の軸解釈に関係してそうなのだが、そこがイマイチわからない。

MCA2010=2021 では、pxx

7 参考文献

Robette N. (2023), GDAtools : Geometric Data Analysis in R, version 2.0, https://nicolas-robette.github.io/GDAtools/

Le Roux B., Rouanet H. (1998). Interpreting Axes in MCA: Method of the contributions of points and deviations, in Visualization of categorical data, ed. J. Blasius & M. Greenacre, Academic Press.