myscratchpad

いつの間にか、psychパッケージの因子分析のオプションが激増していた。William Revelle先生に感謝。


# データの準備
dat <- read.table("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1261838567", header = T)
# UCLAのページから拝借。感謝。
## 0をNAで置換する。この記事のコメントで教えていただいた。感謝
dat2 <- dat
dat2[dat2==0] <- NA
dat3 <- na.omit(dat2) # NAを除外する


# 記述統計
## 変数名を保存しておく
orname <- data.frame(itnumber=names(dat), itname=c("INSTRUC WELL PREPARED","INSTRUC SCHOLARLY GRASP","INSTRUCTOR CONFIDENCE","INSTRUCTOR FOCUS LECTURES","INSTRUCTOR USES CLEAR RELEVANT EXAMPLES","INSTRUCTOR SENSITIVE TO STUDENTS","INSTRUCTOR ALLOWS ME TO ASK QUESTIONS","INSTRUCTOR IS ACCESSIBLE TO STUDENTS OUTSIDE CLASS","INSTRUCTOR AWARE OF STUDENTS UNDERSTANDING","I AM SATISFIED WITH STUDENT PERFORMANCE EVALUATION","COMPARED TO OTHER INSTRUCTORS, THIS INSTRUCTOR IS","COMPARED TO OTHER COURSES THIS COURSE WAS"))
dscrb <- data.frame(ITEM=orname$itname, MEAN=mean(dat3), SD=sd(dat3), N=sapply(dat3, length))
dscrb

# スクリープロットと平行分析
library(psych)
fa.parallel(dat3)

# psychパッケージ、fa関数で最尤解の因子分析。プロマックス回転
library(psych)
factMLres <- fa(dat3, 3 ,rotate="promax", fm="ml", scores=T) # 因子得点を出すためにはscores=Tとしておく
print(factMLres, cut=0, sort=T, digits=5) # factanal関数と違って、fa関数はきちんと並べ替えされる

## factanal関数と比較してみる
fctanlres <- factanal(data=dat3, factors=3, rotation="promax", scores="regression")
print(fctanlres, cut=0.00001, digits=5)

## いずれにせよ自由度33のカイ二乗値=135.35, p < .001なので、この因子モデルは適合しない。因子数を増やして再度分析する
factMLres4 <- fa(dat3, 4 ,rotate="promax", fm="ml", scores=T)
factMLres5 <- fa(dat3, 5 ,rotate="promax", fm="ml", scores=T)
factMLres6 <- fa(dat3, 6 ,rotate="promax", fm="ml", scores=T)
## カイ二乗値とp値だけとりだす
factMLres$STATISTIC
factMLres$PVAL
factMLres4$STATISTIC
factMLres4$PVAL
factMLres5$STATISTIC
factMLres5$PVAL
factMLres6$STATISTIC
factMLres6$PVAL
## 5因子解が適合した (ただし、カイ二乗値による結果が絶対ではない。おそらく3因子解が適切だろう)
print(factMLres5, cut=0, sort=T, digits=5)
## 因子得点
factMLres5$scores



# psychパッケージのfa関数はfactanal関数で出力されなかったものも出てくるので便利。負荷量は一致するが、どうも平方和の計算がなんか違う。因子の順番を変えて計算しているようだ？

spssでは標準で出力される (らしい) Kaiser-Meyer-Olkinの標本妥当性と、Bartlettの球面性検定の関数を見つけたのでメモしておく

# Kaiser-Meyer-Olkinの標本妥当性
青木先生のページから。感謝。
kmo(dat) # datは数値のみのデータフレームか行列

全く同じものがこのページにもあった
kmo.test <- function(df)
{
        cor.sq = cor(df)^2
    cor.sumsq = (sum(cor.sq)-dim(cor.sq)[1])/2
    library(corpcor)
    pcor.sq = cor2pcor(cor(df))^2
    pcor.sumsq = (sum(pcor.sq)-dim(pcor.sq)[1])/2
    kmo = cor.sumsq/(cor.sumsq+pcor.sumsq)
    return(kmo)
}


# Bartlettの球面性
# psychパッケージは偉大すぎる
library(psych)
cortest.bartlettdat)
#cortest.bartlett(cor(dat)) # 相関行列を与えること、とマニュアルにはあるが、spssの結果と一致するのはローデータを与えたときのようだ

主因子法 (principal factor) と主軸法 (principal axis) って同じなのかな？

dat <- read.table("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1261662305", header = T)
# 小塩先生のページを参考にさせていただいた。感謝
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/da09_01.html

# 変数名が長いので別にしておこう
ornames <- names(dat)
dat2 <- dat
names(dat2) <- c(paste("A_0", 1:9, sep = ""), paste("A_", 10:23, sep = ""))

# まずは記述統計
mres <- mean(dat2)
sdres <- sd(dat2)
resdes <- data.frame(name=names(dat2), mean= mres, sd = sdres, mm1sd=mres-sdres, mp1sd=mres+sdres)
resdes
## 以下でもいい
#summary(dat2)
#sd(dat2)
#library(psych)
#describe(dat2)

# 固有値を調べプロットする
evres <- eigen(cor(dat2))
evres$value # 固有値を表示
plot(evres$value, type="b") # プロット。いわゆるスクリープロット
# psychパッケージのfa.parallelが便利
fa.parallel(dat2)

# psychパッケージのfactor.pa関数により主因子解の因子分析
factpares <- factor.pa(dat2, 3 ,rotate="promax")
print(factpares, cut=0, sort=T, digits=5)

# 小塩先生のページどおり、一部の項目を除外して因子分析
dat3 <- dat2[c(-19, -8, -23)] # A19, A08, A23を除外
fa.parallel(dat3)
factpares2 <- factor.pa(dat3, 3, rotate="promax")
print(factpares2, cut=0, sort=T, digits=5)

微妙に小塩先生のspssの結果と一致しないがまあいいだろう
因子分析なんてそんなものだ、ということにしよう

## 以下と同じ分析のはずだが、なぜか第2主軸の符号が逆になってる。たぶん上のでいいんだろう
fares <- fa(dat3, 3, rotate="promax", fm="pa")
print(fares, cut=0, sort=T, digits=5)

自分用に手順を書いておく

#1. まず何はともあれ実行しておく。因子数は理論に従う
result <- factanal(dat, factors=3, rotation="varimax", scores="regression")

#2. その上で、他の因子数モデルを実行し、カイ二乗値を比較する
result3 <- factanal(..factors=3,..)
result4 <- factanal(..factors=4,..)
...
## カイ二乗値が有意水準を超えたところが適切な因子数

#3. 適切な因子数で分析し、結果を見る
print(result, cut=0.001)
## cut=0.001とすることで、負荷量で0.1未満の数値も表示される

共通性
1-result%uniqueness
## 1から独自性を引く

寄与率
resultの中の
Proportion Var
Cumulative Var # 累積寄与率

因子負荷量
print(sort.loadings(result), cut=0.001)
## 青木先生の並べ替え関数を使わせていただく。感謝

## プロマックス回転時の因子間相関
fcor <- factanal(dat,factors=3,rotation="none")
pro <- promax(loadings(fcor),m=3)
solve(t(pro$rotmat)%*%pro$rotmat)

Rjipwiki: 　Rで共分散構造分析・構造方程式モデル　
psychパッケージのガイド (pdf) 　
semパッケージのマニュアル　
J. Fox (2006), Structural Equation Modeling With the sem Package in R. Structural Equation Modeling, 13:465-486.　　
Quick-Rの因子分析のページ　
James H. Steigerの確認的因子分析ガイド (from statpower)