library(tidyverse)
library(GDAtools)
library(showtext)
showtext_auto(TRUE)

データ読み込み

load("taste_d_J.rda")
.d_J
## # A tibble: 1,253 × 9
##       ID Isup   TV            Film               Art   Eat   Gender Age   Income
##    <int> <fct>  <fct>         <fct>              <fct> <fct> <fct>  <fct> <fct> 
##  1     1 Active TV-メロドラマ 映画-アクション    芸術… 外食… 女性   55-64 £20-29
##  2     2 Active TV-メロドラマ 映画-ホラー        芸術… 外食… 女性   45-54 <£9   
##  3     3 Active TV-自然       映画-アクション    芸術… 外食… 女性   55-64 <£9   
##  4     4 Active TV-メロドラマ 映画-時代劇        芸術… 外食… 女性   65+   £10-19
##  5     5 Active TV-コメディー 映画-ホラー        芸術… 外食… 女性   35-44 £10-19
##  6     6 Active TV-コメディー 映画-ホラー        芸術… 外食… 女性   18-24 <£9   
##  7     7 Active TV-ニュース   映画-アクション    芸術… 外食… 女性   25-34 £10-19
##  8     8 Active TV-ニュース   映画-ドキュメンタ… 芸術… 外食… 男性   65+   £10-19
##  9     9 Active TV-メロドラマ 映画-時代劇        芸術… 外食… 女性   65+   <£9   
## 10    10 Active TV-ニュース   映画-アクション    芸術… 外食… 女性   65+   £10-19
## # ℹ 1,243 more rows

Isup == “Active” の1215行を切り出す

.d_J %>% filter(Isup=="Active") -> .d0

MCAを実行

res.speMCA <- speMCA(.d0[,3:6])

resultから分散を計算する

個体、座標を取り出す

res.speMCA$ind$coord -> coord_ind

追加要素を接続する

cbind(coord_ind[,1:3],.d0[,7:9]) -> ind_coord_sup

分散の計算

個体空間における性別平均点を計算

ind_coord_sup %>% as_tibble() %>% group_by(Gender) %>% summarize_at(.vars = 1:3,.funs = mean) -> Gender_mean # sup=Genderの平均点座標/個体空間
Gender_mean
## # A tibble: 2 × 4
##   Gender   dim.1  dim.2  dim.3
##   <fct>    <dbl>  <dbl>  <dbl>
## 1 男性    0.112   0.158  0.300
## 2 女性   -0.0822 -0.115 -0.219

性別の度数を計算

ind_coord_sup %>% as_tibble() %>% group_by(Gender) %>% count(Gender) -> Gender_wait
Gender_wait
## # A tibble: 2 × 2
## # Groups:   Gender [2]
##   Gender     n
##   <fct>  <int>
## 1 男性     513
## 2 女性     702

全体雲の平均点がゼロであることを確認

ind_coord_sup %>% as_tibble()  %>% summarize_at(.vars = 1:3,.funs = mean) #  これはゼロとうこと。
## # A tibble: 1 × 3
##      dim.1    dim.2     dim.3
##      <dbl>    <dbl>     <dbl>
## 1 2.13e-16 1.14e-16 -1.08e-16

性別/軸ごとの分散を計算

ind_coord_sup %>% as_tibble() %>% left_join(Gender_mean,by = "Gender")  %>%  
  mutate(d1_2=(dim.1.x - dim.1.y)^2,d2_2=(dim.2.x - dim.2.y)^2,d3_2=(dim.3.x - dim.3.y)^2,) %>% 
  select(4,10:12) %>% group_by(Gender) %>% 
  summarize_at(.vars = 1:3,.funs = sum) %>% left_join(Gender_wait,by="Gender") %>% mutate(var1=d1_2/n,var2=d2_2/n,var3=d3_2/n) %>% select(-c(2:4)) -> dim_var
dim_var
## # A tibble: 2 × 5
##   Gender     n  var1  var2  var3
##   <fct>  <int> <dbl> <dbl> <dbl>
## 1 男性     513 0.292 0.253 0.257
## 2 女性     702 0.464 0.392 0.261

Vwithin 郡内分散（性別内分散）

性別ごとの軸内分散の加重平均

dim_var %>% summarise(within_1 = weighted.mean(x=var1,w=n),within_2 = weighted.mean(x=var2,w=n),within_3 = weighted.mean(x=var3,w=n))
## # A tibble: 1 × 3
##   within_1 within_2 within_3
##      <dbl>    <dbl>    <dbl>
## 1    0.391    0.333    0.259

Vbetween 群間分散（原点-性別間分散）

ind_coord_sup %>% as_tibble() %>% group_by(Gender) %>% summarize_at(.vars = 1:3,.funs = mean) %>% left_join(Gender_wait,by="Gender") %>% 
  mutate(Gdim.1 = (dim.1^2),Gdim.2 = (dim.2^2),Gdim.3 = (dim.3^2)) %>% 
  summarize(between1=(weighted.mean(x=Gdim.1,w=n)),between2=(weighted.mean(x=Gdim.2,w=n)),between3=(weighted.mean(x=Gdim.3,w=n))) -> Vbetween_gender
Vbetween_gender
## # A tibble: 1 × 3
##   between1 between2 between3
##      <dbl>    <dbl>    <dbl>
## 1  0.00924   0.0182   0.0657

η2 相関比（軸ごとの慣性に対する、between分散の比率）

この値がおおきければ、そこんは差異がある。

軸の慣性 = 固有値

res.speMCA$eig$eigen[1:3] -> eig1_3

Vbetween を軸ごとの慣性で割る

Vbetween_gender %>% as.matrix() %*% diag(1/eig1_3)
##            [,1]       [,2]      [,3]
## [1,] 0.02307637 0.05178419 0.2020366