主成分分析
佐藤圭一
2023-09-15
1 例題
事例：①年齢、②学歴、③財産、④職業、⑤性別、⑥人種による不公不公平感についてきいたデータがあるとする。「総合的な不公平感」に関する合成変数を作成したい。
1.1 模擬データの作成
dat <- data.frame(
age_un = c(4,5,2,3,3,2,2,2,4,3,1,1,2,5,2,4,1,3,3,2),
edu_un = c(4,2,1,1,2,0,2,1,3,1,1,1,2,4,2,3,2,2,1,1),
wealth_un = c(3,4,2,1,3,1,1,1,2,2,1,1,1,4,1,2,2,2,4,2),
job_un = c(3,4,3,2,3,4,2,2,4,3,3,3,3,4,4,3,5,2,3,2),
sex_un = c(5,2,1,1,2,1,2,1,4,1,1,1,2,4,2,4,2,2,1,1),
race_un = c(5,5,4,2,5,2,1,2,4,4,1,2,1,5,2,4,3,4,5,4)
)
# IDをrownamesとしてつける
rownames(dat) <- 1:201.2 記述統計
ひとまず記述統計を見てみる（良い習慣として）
library(psych)
describe(dat) # 記述統計## vars n mean sd median trimmed mad min max range skew kurtosis
## age_un 1 20 2.70 1.22 2.5 2.62 0.74 1 5 4 0.38 -0.93
## edu_un 2 20 1.80 1.06 2.0 1.69 1.48 0 4 4 0.63 -0.40
## wealth_un 3 20 2.00 1.08 2.0 1.88 1.48 1 4 3 0.72 -0.84
## job_un 4 20 3.10 0.85 3.0 3.06 1.48 2 5 3 0.31 -0.76
## sex_un 5 20 2.00 1.26 2.0 1.81 1.48 1 5 4 1.06 -0.23
## race_un 6 20 3.25 1.48 4.0 3.31 1.48 1 5 4 -0.22 -1.57
## se
## age_un 0.27
## edu_un 0.24
## wealth_un 0.24
## job_un 0.19
## sex_un 0.28
## race_un 0.33pairs(dat) # 変数ペアのプロット
1.3 主成分分析
library(psych)
# 主成分分析
out.pr <- principal(dat, # データ
nfactors = 6, # 主成分の数（ここでは変数の数）
rotate = "none", # 軸の回転をしないオプションの場合（教科書の例）
scores = T) # 主成分得点を計算する
# 結果
out.pr # 結果の表示(全体)## Principal Components Analysis
## Call: principal(r = dat, nfactors = 6, rotate = "none", scores = T)
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 PC3 PC4 PC5 PC6 h2 u2 com
## age_un 0.88 -0.15 -0.20 0.41 0.01 0.02 1 2.2e-16 1.6
## edu_un 0.83 0.47 -0.20 -0.14 -0.11 0.10 1 7.8e-16 1.9
## wealth_un 0.83 -0.48 0.16 -0.10 -0.22 -0.06 1 1.1e-16 1.9
## job_un 0.37 0.30 0.88 0.08 0.02 0.01 1 0.0e+00 1.6
## sex_un 0.82 0.52 -0.18 -0.06 0.10 -0.11 1 5.6e-16 1.9
## race_un 0.82 -0.50 0.05 -0.17 0.22 0.04 1 7.8e-16 1.9
##
## PC1 PC2 PC3 PC4 PC5 PC6
## SS loadings 3.62 1.09 0.91 0.23 0.12 0.03
## Proportion Var 0.60 0.18 0.15 0.04 0.02 0.00
## Cumulative Var 0.60 0.78 0.94 0.98 1.00 1.00
## Proportion Explained 0.60 0.18 0.15 0.04 0.02 0.00
## Cumulative Proportion 0.60 0.78 0.94 0.98 1.00 1.00
##
## Mean item complexity = 1.8
## Test of the hypothesis that 6 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0
## with the empirical chi square 0 with prob < NA
##
## Fit based upon off diagonal values = 1# 主成分得点(ここでは最初の6行のみ)
head(out.pr$scores) ## PC1 PC2 PC3 PC4 PC5 PC6
## 1 1.7433800 0.93033103 -1.06106792 -1.2808148 0.504323548 -1.4705168
## 2 1.3001708 -1.23695172 0.94452437 1.8774056 -1.084327432 0.2807845
## 3 -0.3905668 -0.89996827 0.37195763 -0.7663059 0.915187048 0.6417905
## 4 -0.8276661 -0.30932016 -1.17849121 1.6602658 0.008434758 0.7030033
## 5 0.5705358 -0.93242858 0.02033546 -0.9871340 0.256657312 0.5119055
## 6 -1.0064446 0.03897121 1.47096778 1.5804463 1.301555270 -2.1973111# スクリープロット
plot(out.pr$values, type = "o") 
# 結果のまとめ
out <- rbind(out.pr$Structure, # 主成分得点
eigen = out.pr$values, # 固有値
prop = out.pr$values/sum(out.pr$values), # 寄与率
cum.prop = cumsum(out.pr$values/sum(out.pr$values))) # 累積寄与率
out## PC1 PC2 PC3 PC4 PC5 PC6
## age_un 0.8783263 -0.1474025 -0.20336218 0.40615398 0.01128868 0.019253825
## edu_un 0.8318908 0.4738901 -0.20370388 -0.13966698 -0.11100724 0.100305448
## wealth_un 0.8285209 -0.4763983 0.15974659 -0.09648038 -0.22006402 -0.057812142
## job_un 0.3669104 0.3011185 0.87629681 0.07778280 0.02401096 0.013475069
## sex_un 0.8173710 0.5239074 -0.18161090 -0.05757682 0.09636075 -0.108824329
## race_un 0.8178419 -0.4998052 0.05214383 -0.17373744 0.21665073 0.038577083
## eigen 3.6215301 1.0882112 0.91196802 0.23332634 0.11767768 0.027286640
## prop 0.6035884 0.1813685 0.15199467 0.03888772 0.01961295 0.004547773
## cum.prop 0.6035884 0.7849569 0.93695156 0.97583928 0.99545223 1.000000000# 結果を保存する
write.csv(out, file = "out.csv")2 単純加算
主成分分析を行う代わりに、変数の単純加算を行う場合がある（今回の場合、un = age_un + edu_un + … + race_unなど。）。主成分負荷量をすべて１と置いた主成分得点を計算しているものとして、特殊な主成分分析と位置づけることもできる。単純加算が可能かどうかを判定するためには、クロンバックのαがしばしばもちいられる。
2.1 クロンバックのα
library(psych)
alpha(dat,
check.keys = T) # 逆転項目をチェックし、自動的に反転する。##
## Reliability analysis
## Call: alpha(x = dat, check.keys = T)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.86 0.86 0.93 0.5 5.9 0.047 2.5 0.9 0.44
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.74 0.86 0.94
## Duhachek 0.77 0.86 0.95
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## age_un 0.81 0.81 0.91 0.46 4.2 0.066 0.067 0.41
## edu_un 0.82 0.82 0.87 0.47 4.4 0.057 0.067 0.41
## wealth_un 0.82 0.82 0.87 0.47 4.4 0.060 0.070 0.43
## job_un 0.89 0.89 0.96 0.63 8.5 0.042 0.039 0.65
## sex_un 0.83 0.82 0.86 0.47 4.5 0.054 0.067 0.45
## race_un 0.83 0.82 0.90 0.48 4.6 0.056 0.066 0.43
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## age_un 20 0.86 0.85 0.81 0.79 2.7 1.22
## edu_un 20 0.82 0.82 0.84 0.73 1.8 1.06
## wealth_un 20 0.83 0.82 0.82 0.75 2.0 1.08
## job_un 20 0.41 0.47 0.33 0.27 3.1 0.85
## sex_un 20 0.81 0.81 0.83 0.70 2.0 1.26
## race_un 20 0.83 0.80 0.79 0.70 3.2 1.48
##
## Non missing response frequency for each item
## 0 1 2 3 4 5 miss
## age_un 0.00 0.15 0.35 0.25 0.15 0.10 0
## edu_un 0.05 0.40 0.35 0.10 0.10 0.00 0
## wealth_un 0.00 0.40 0.35 0.10 0.15 0.00 0
## job_un 0.00 0.00 0.25 0.45 0.25 0.05 0
## sex_un 0.00 0.45 0.35 0.00 0.15 0.05 0
## race_un 0.00 0.15 0.25 0.05 0.30 0.25 0# 逆転項目があるため、警告が出る。# 合計得点により合成変数を作成する場合
un_syn <-
(dat$age_un + dat$edu_un + dat$wealth_un +
dat$job_un + dat$sex_un + dat$race_un)/6
# こちらと、上記の第一主成分は強い相関がある
plot(x = un_syn,
y = out.pr$scores[,1])
cor.test(un_syn, out.pr$scores[,1]) # r = .99の相関##
## Pearson's product-moment correlation
##
## data: un_syn and out.pr$scores[, 1]
## t = 69.409, df = 18, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9951866 0.9992796
## sample estimates:
## cor
## 0.9981371参考：1, 3, 6番目の変数のみでαを計算してみたい場合
以下のように使用する変数の列をc(… )で選択する。なお同様の選択の方法は、主成分分析で特定の変数のみを選択したい場合にも使える）
alpha(dat[ ,c(1,3,6)],
check.keys = T)2.2 参考：変数のペアごとに相関をプロットする
主成分分析や単純加算などを検討するうえで、変数のペアごとに相関係数を出して検討したい場合がある。以下のスクリプトは、相関係数の絶対値（つまり正負の情報は表示しない）とヒストグラムを同時にプロットする便利な関数である。
# 以下の関数を先に読み込む
## ヒストグラムを計算する関数
panel.hist <- function(x, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)}
## 変数間の相関係数の絶対値を計算する関数
## with size proportional to the correlations.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
# 変数のペアごとのヒストグラムと相関の同時プロット
pairs(dat,
diag.panel = panel.hist,
lower.panel = panel.smooth,
upper.panel = panel.cor,
gap=0, row1attop = FALSE)
3 PCAの結果を図示する
実際の論文においてしばしば使用される主成分分析の用い方のフローに基づいたスクリプトを示す。具体的には、はじめに主成分分析を行い、結果を解釈したうえで、主成分得点を保存し、各ケースをプロットする。
主成分分析
# 主成分分析
library(psych)
pr.out <- principal(dat,
nfactors = 2, # 二つの主成分に要約
scores = T,
rotate = "none")第二主成分の負荷量は、「仕事」「教育」「性別」には正、「人種」「収入」には負がかかっている。 日本の文脈においては、前者の不平等が特に共有される一方、後者は不平等がないと理解されてきた。 よって第二主成分を伝統的日本型不平等への認識を示すものと考える。なお、第一主成分は総合不平等認識度とする。
主成分得点の保存
# 主成分得点を別のオブジェクトに保存する
pr.score <- data.frame(pr.out$scores)
colnames(pr.score) <- c("general.inequality","Japanese.inequality")
# ケースをプロットしてみる
library(car)
plot(x = pr.score$general.inequality,
y = pr.score$Japanese.inequality,
pch = 16,
xlab = "総合的不平等認知",
ylab = "日本の典型的不平等への認知")
abline(v = 0, col = "red") # 見やすいように平均値（０）に線を引く
abline(h = 0, col = "red") # 見やすいように平均値（０）に線を引く
pointLabel(
x = pr.score$general.inequality,
y = pr.score$Japanese.inequality,
labels=rownames(dat)) # 図にIDもつけてみる。 
総合的不平等を平均以上に認識している層で見た場合、より高く総合的不平等を認識している層は日本の典型的不平等への認知も高い。
4 ポリコック相関を用いた主成分分析
主成分分析は比例尺度の相関行列に対して用いることが想定されている手法である。順序尺度のデータに対してはポリコック相関(2値)、テトラコック相関（3値）に基づいた主成分分析を行うことが考えられる。
4.1 順序尺度に基づく模擬データの作成
# 年齢に関する不平等感に関する回答で平均以上を1、以下を0とする二値データを作成する
age_un.bin <- ifelse(dat$age_un >= mean(dat$age_un), yes = 1, no = 0)
# 教育に関する不平等感に関する回答で1と2を1, 3を2, 4と5を3に分ける
edu_un.3cat <- cut(dat$edu_un,
breaks = c(-Inf, 2.5, 3.5, Inf),
labels = c("low","middle","high"),
right = F, ordered_result = T)
# 新しいデータセットをつくる
dat2 <- data.frame(age_un.bin = age_un.bin,
edu_un.3cat = as.numeric(edu_un.3cat), # 数値にすることに注意
dat[,3:6])4.2 通常の相関、ポリコック、テトラコック相関のミックス
cor.result <- mixedCor(data = dat2,
c = 3:6, # 連続変数の行番号
p = 2, # 順序変数の行番号
d = 1) # 二値変数の行番号cor.result$rho # 相関の結果## age_un.bin edu_un.3cat wealth_un job_un sex_un race_un
## age_un.bin 1.0000000 0.9936381 0.8153113 0.0000000 0.5984493 0.8877645
## edu_un.3cat 0.9936381 1.0000000 0.6280214 0.3171803 1.0000000 0.6458404
## wealth_un 0.8153113 0.6280214 1.0000000 0.2869585 0.3892495 0.8908902
## job_un 0.0000000 0.3171803 0.2869585 1.0000000 0.2948839 0.1874756
## sex_un 0.5984493 1.0000000 0.3892495 0.2948839 1.0000000 0.4238404
## race_un 0.8877645 0.6458404 0.8908902 0.1874756 0.4238404 1.0000000# 主成分分析
pr.out2 <- principal(
r = cor.result$rho, # インプットには相関マトリクスも取れるため、これを利用
nfactors = 2,
rotate = "none")
# 結果
pr.out2 ## Principal Components Analysis
## Call: principal(r = cor.result$rho, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2 com
## age_un.bin 0.96 -0.30 1.01 -0.008 1.2
## edu_un.3cat 0.97 0.19 0.97 0.027 1.1
## wealth_un 0.85 -0.22 0.77 0.235 1.1
## job_un 0.31 0.80 0.73 0.266 1.3
## sex_un 0.77 0.40 0.76 0.241 1.5
## race_un 0.87 -0.31 0.85 0.152 1.2
##
## PC1 PC2
## SS loadings 4.01 1.07
## Proportion Var 0.67 0.18
## Cumulative Var 0.67 0.85
## Proportion Explained 0.79 0.21
## Cumulative Proportion 0.79 1.00
##
## Mean item complexity = 1.2
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.14
##
## Fit based upon off diagonal values = 0.95参考までにもとの連続変数を用いた主成分分析の結果を、比較のために再掲する。
pr.out## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = "none", scores = T)
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2 com
## age_un 0.88 -0.15 0.79 0.207 1.1
## edu_un 0.83 0.47 0.92 0.083 1.6
## wealth_un 0.83 -0.48 0.91 0.087 1.6
## job_un 0.37 0.30 0.23 0.775 1.9
## sex_un 0.82 0.52 0.94 0.057 1.7
## race_un 0.82 -0.50 0.92 0.081 1.7
##
## PC1 PC2
## SS loadings 3.62 1.09
## Proportion Var 0.60 0.18
## Cumulative Var 0.60 0.78
## Proportion Explained 0.77 0.23
## Cumulative Proportion 0.77 1.00
##
## Mean item complexity = 1.6
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.09
## with the empirical chi square 4.67 with prob < 0.32
##
## Fit based upon off diagonal values = 0.97参考文献
永吉希久子，2016，『行動科学の統計学：社会調査のデータ分析』共立出版．