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.speMCA0 <- speMCA(.d0[,3:6])
res.speMCA0 %>% flip.mca(dim = c(1,2)) -> res.speMCA

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(Age) %>% summarize_at(.vars = 1:3,.funs = mean) -> Age_mean # sup=Ageの平均点座標/個体空間
Age_mean
## # A tibble: 6 × 4
##   Age     dim.1   dim.2      dim.3
##   <fct>   <dbl>   <dbl>      <dbl>
## 1 18-24  0.589  -0.332   0.0145   
## 2 25-34  0.272  -0.191  -0.0143   
## 3 35-44  0.0890 -0.0533  0.0522   
## 4 45-54 -0.0538 -0.0697 -0.0465   
## 5 55-64 -0.367   0.101  -0.0130   
## 6 65+   -0.281   0.359  -0.0000491

AGe(年齢)の度数を計算

ind_coord_sup %>% as_tibble() %>% group_by(Age) %>% count(Age) -> Age_weight
Age_weight
## # A tibble: 6 × 2
## # Groups:   Age [6]
##   Age       n
##   <fct> <int>
## 1 18-24    93
## 2 25-34   248
## 3 35-44   258
## 4 45-54   191
## 5 55-64   183
## 6 65+     242

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

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(Age_mean,by = "Age")  %>%  
  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(5,10:12) %>% group_by(Age) %>% 
  summarize_at(.vars = 1:3,.funs = sum) %>% left_join(Age_weight,by="Age") %>% mutate(var1=d1_2/n,var2=d2_2/n,var3=d3_2/n) %>% select(-c(2:4)) -> dim_var
dim_var
## # A tibble: 6 × 5
##   Age       n  var1  var2  var3
##   <fct> <int> <dbl> <dbl> <dbl>
## 1 18-24    93 0.192 0.395 0.258
## 2 25-34   248 0.308 0.322 0.293
## 3 35-44   258 0.337 0.288 0.341
## 4 45-54   191 0.360 0.318 0.312
## 5 55-64   183 0.312 0.246 0.409
## 6 65+     242 0.340 0.314 0.308

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.321    0.307    0.324

Vbetween 群間分散（原点-年齢間分散）

ind_coord_sup %>% as_tibble() %>% group_by(Age) %>% summarize_at(.vars = 1:3,.funs = mean) %>% left_join(Age_weight,by="Age") %>% 
  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_Age
Vbetween_Age
## # A tibble: 1 × 3
##   between1 between2 between3
##      <dbl>    <dbl>    <dbl>
## 1   0.0798   0.0444  0.00100

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

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

軸の慣性 = 固有値

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

Vbetween を軸ごとの慣性で割る

Vbetween_Age %>% as.matrix() %*% diag(1/eig1_3)
##           [,1]      [,2]        [,3]
## [1,] 0.1992784 0.1264439 0.003084088