Rで因子分析と分散共分散分析をする。
- 因子分析は、潜在要因をみつける方法。
- 主成分分析と回帰分析と異なり観測変数xを説明する潜在変数をみつける方法（下のx1とx2を見つける）
y1=a*x1+a*x2
y2=a*x1+a*x2
y3=a*x1+a*x2
y4=a*x1+a*x2
-
- 心理学等では、「瞬発筋力x1」や「持久力x2」や「器用さx3」など観測できない構成概念を扱うことが重要なためこの方法が発達した。
- ただし、係数の決め方や変数の決め方に不定性（同じパラメータでおなじ説明力）があるので、軸の回転をおこなうことで、説明しやすい解や理解しやすい解をさぐる。
- 各潜在要因が直交しない解がだすことができるのも因子分析の特徴
- 自由に色々選べるので調整がききやすい
- 代表的な自動回転法 promax回転（非直行）やvarimax回転（直行）。構造を目指す（出来る限り観測変数が一つの説明の要因で説明できるようにする回転。潜在要因と観測要因の関係が明確になるので説明しやすい）。
- 要因の説明力をたかめると、個別データをバラバラにする野の違いもある（主成分と因子分析）
-
- Rでは、fa関数（psychパッケージ）
参考：Rで因子分析 商用ソフトで実行できない因子分析のあれこれ
http://www.slideshare.net/simizu706/r-42283141
-
-
- Step1 下処理と標準化 scale(X)関数．na.omit(X), 正規分布っぽくする．欠損値除去
- Step2 因子数の自動推定(vss, parallel関数)
-
install.packages("psych")
library(psych)
#最低値がでる。BICがおすすめらしい
map <- vss(rr2, fm="minres")
#最大値がでる。BICがおすすめらしい
para <- fa.parallel(rr2)
-
-
-
- ほかには主成分析→おりまがり+-1を調べる plot(eigen(cor(X))$values)
- Step3 因子分析
-
-
#例：8要因でoblimin回転（promaxの進化版?）
v<-fa(rr3,nfactor=8,rotate = "oblimin")
#例：8要因でvarimax回転（直行）
v2<-fa(rr3,nfactor=8,rotate = "varimax")
#例：8要因でpromax回転（非直行）
v3<-fa(rr3,nfactor=8,rotate = "promax")
-
-
- Step4: 結果の取り出し（各潜在要因のうち因子スコアが大きい観測量をとりだし）。
-
for(i in 1:8){
cat(labels*12[i],"\n")
s<-(v2$loadings[,i][rev(order(abs(v2$loadings[,i])))][1:10])
print(s)
#print(paste(names(s),collapse=", "))
#print(paste(sprintf("%.3f",s),collapse=", "))
}
-
-
- RMSEA(root mean square error of approximation) 0.05 OK 0.1> No
- 共通因子スコアと独自因子スコア
-
- 分散共分散分析（SEM）
- (1)パッケージのインストール
#install.pakcages("lavaan")
#install.packages("semPlot")
#install.packages("knitr")
#install.packages("data.table")#library(semPlot)
library(semPlot)
library(lavaan)
-
- （２）練習1
- 例１:線形回帰
- 構造方程式（xがyを説明する場合）はy~x
- 観測方程式（潜在要因fが観測量xを説明する場合）はy=~xとかく
- 相関なしは潜在要因f1とf2が無相関はf1~0*f2とする
- 例１:線形回帰
- （２）練習1
#通常の回帰
x<-rnorm(1000,mean=0,sd=2)
y<-3*x+rnorm(1000,mean=0,sd=0.1)model1<-"y~x"
X<-cbind(x,y)
getCov(X)
#データ直接入力モード
fit<-sem(model=model1,data=X,estimator="ML")
#サマリー
summary(fit)
#サマリー、精度評価指標つき
summary(fit,fit.measure=TRUE)
#解析結果のグラフ化
semPaths(fit,whatLabels="stand",optimizeLatRes=T)
#解析結果のとりだし
parameterEstimates(object=fit)
#解析結果のとりだし標準化
standardizedSolution(object=fit)
#体重gの列を除く(体重(kg)と相関係数1の関係)
rr2b<-(rr2[,-which(colnames(rr2)=="体重(g)")])
#分散共分散行列
r_d<-cov(rr2b)
r_d3<-r_d
#データのラベルを取得
n1<-labels(rr2b)2
#すべてのデータラベルの観測方程式の（y1=~長距離走; y1=~ 短距離走; y1~=ソフトボールなげなど）を作る。
model_test<-paste("y1=~",n1,sep="")
#すべてのデータラベルの観測方程式の（y2=~長距離走; y2=~ 短距離走; y2~=ソフトボールなげなど）を作る。
model_testb<-paste("y2=~",n1,sep="")
#潜在要因y1とy2は無相関という条件をいれる
model_testc<-'y1~0*y2'
#モデルをまとめる
model_test2<-paste(c(model_test,model_testb,model_testc),collapse="\n")
#データ数の取得
data_no<-length(rr2b[,1])
#SEMを実行（分散共分散行列モード、データ数の指定が必要）
fit<-sem(model=model_test2,sample.cov=r_d3,estimator="ML",sample.nobs=data_no)
#semPaths(fit,whatLabels="stand",optimizeLatRes=T,nCharNodes=12)
#parameterEstimates(object=fit)
#標準化回帰係数の取り出し
standardizedSolution(object=fit)
-
- 実際の解析
- とりあえず、探索的因子分析（因子分析の結果）を利用して観測方程式を作ってみる。
- 潜在要因をさらに朝ごはんや学校成績で説明するモデル（構造方程式）
- 実際の解析
長距離走<-(持久力),(瞬発力)
短距離走<-(瞬発力)
ソフトボール投げ<-(瞬発力)
…
朝ごはん->(持久力)
学校成績->(持久力)
運動習慣->(持久力)
朝ごはん->(瞬発力)
学校成績->(瞬発力)
運動習慣->(瞬発力)
朝ごはん->(筋力)
学校成績->(筋力)
運動習慣->(筋力)
*1:v3$loadings
