検証的因子分析 Last modified: Oct 07, 2008
目的
検証的因子分析（プロマックス回転）を行う
使用法
cfa(r, n, loc, lim=0.999, init.val=0.5)
引数
r 相関係数行列
n サンプルサイズ
loc 因子負荷量または設定情報（使用例を参照のこと）
lim=0.999 推定パラメータの値域 (-lim, lim)
init.val=0.5 推定パラメータの初期値
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/cfa.R", encoding="euc-jp")
# 検証的因子分析（確認的因子分析）
cfa <- function(r, # 相関係数行列
n, # サンプルサイズ
loc, # 因子負荷量または設定情報
lim=0.999, # 推定パラメータの値域 (-lim, lim)
init.val=0.5) { # 推定パラメータの初期値
tr <- function(x) sum(diag(x)) # トレース
get.chi.sq <- function(par, scaler=TRUE) # optim 用の関数
{
loadings[loadings != 0] <- par[1:p1] # 因子負荷量行列
fac.cor <- diag(nfac) # 因子間相関行列
fac.cor[lower.tri(fac.cor)] <- par[-(1:p1)]
fac.cor <- fac.cor+t(fac.cor)
diag(fac.cor) <- 1
rho <- loadings%*%fac.cor%*%t(loadings) # 母相関係数行列
diag(rho) <- 1
temp <- det(rho) # 行列式が負の値になるときの対策
retval <- tr(solve(rho)%*%r)+ # カイ二乗値の一部（定数部を除く）
(if (temp <= 0) 100 else log(temp)) # 100 は根拠のない仮の値
if (scaler) {
return(retval) # optim への戻り値
}
else {
return(list(retval=retval, # 収束後の結果諸々
loadings=loadings, fac.cor=fac.cor, rho=rho))
}
}
p <- nrow(r) # 変数の個数
if (is.matrix(loc) && nrow(loc) == p) { # 因子負荷量として与えられるとき
loadings <- loc # 変数が因子に含まれるところを 1 とする
nfac <- ncol(loadings)
}
else {
nfac <- max(loc) # 因子の個数
loadings <- matrix(0, p, nfac) # 因子負荷量行列
loadings[cbind(1:p, loc)] <- 1 # 変数が因子に含まれるときに 1 とする
}
p1 <- sum(loadings != 0) # 推定すべき因子負荷量の個数
if (length(init.val) == 1) { # 推定すべき因子負荷量の初期値
init.val <- rep(init.val, p1) # 全て同じにする
}
par <- c(init.val, rep(0, nfac*(nfac-1)/2)) # パラメータ初期値（因子負荷量行列と因子間相関行列）
df <- p*(p-1)/2-length(par) # 自由度
ans <- optim(par, get.chi.sq, method="L-BFGS-B", # 最適化（-1 〜 1 の制約付きで）
lower=rep(-lim, length(par)), upper=rep(lim, length(par)))
if (ans$convergence) {
stop(paste("convergence =", ans$convergence, # 何らかの原因で収束しなかったら，理由を明示し停止
"\nmessage =", ans$message))
}
ans2 <- get.chi.sq(ans$par, scaler=FALSE) # 最適パラメータのときの結果
rho <- ans2$rho # 母相関係数行列
chisq <- (ans2$retval-log(det(r))-p)*(n-1) # 定数部を含めたカイ二乗値
P <- pchisq(chisq, df, lower.tail=FALSE) # P 値
z1 <- solve(rho)%*%(r-rho)
z2 <- solve(rho)%*%r
GFI <- 1-tr(z1%*%z1)/tr(z2%*%z2) # GFI
return(list(loadings=ans2$loadings, fac.cor=ans2$fac.cor,
chisq=chisq, P=P, df=df, GFI=GFI,
AGFI=1-p*(p+1)*(1-GFI)/2/df, # AGFI
SRMR=sqrt(sum((rho-r)^2)/p/(p+1)), # SRMR
RMSEA=sqrt((chisq-df)/df/(n-1)) )) # RMSEA
}
使用例
r <- structure(c( # 相関係数行列
1, 0.412, 0.521, 0.358, 0.334, 0.346,
0.412, 1, 0.495, 0.499, 0.293, 0.248,
0.521, 0.495, 1, 0.525, 0.364, 0.323,
0.358, 0.499, 0.525, 1, 0.607, 0.517,
0.334, 0.293, 0.364, 0.607, 1, 0.506,
0.346, 0.248, 0.323, 0.517, 0.506, 1),
.Dim = c(6L, 6L),
.Dimnames = list(c("代数", "幾何", "解析", "英語", "国語", "古文"),
c("代数", "幾何", "解析", "英語", "国語", "古文")))
n <- 100 # サンプルサイズ
loc <- c(2,2,2,1,1,1) # 変数ごとにどの因子に含まれるか因子番号
cfa(r, n, loc) # 変数1〜変数3は第2因子，変数4〜変数6は第1因子に含まれる
loc2 <- matrix(c( # 直接因子パターンを指定するとき
0,1, # 1 番目の変数は第 2 因子
0,1, # 2 番目の変数は第 2 因子
0,1, # 3 番目の変数は第 2 因子
1,1, # 4 番目の変数は第 1 因子と同時に第 2 因子にも含まれるとする
1,0, # 5 番目の変数は第 1 因子
1,0), # 6 番目の変数は第 1 因子
byrow=TRUE, ncol=2)
cfa(r, n, loc2)
出力結果例
> ans <- factanal(covmat=r, factors=2, rotation="promax") # 探索的因子分析を行う
> print(ans, cutoff=0)
Call:
factanal(factors = 2, covmat = r, rotation = "promax")
Uniquenesses:
代数 幾何 解析 英語 国語 古文
0.609 0.577 0.363 0.334 0.393 0.575
Loadings:
Factor1 Factor2
代数 0.031 0.605
幾何 0.032 0.629
解析 -0.026 0.815
英語 0.629 0.251 2 因子構造が見られるが，「英語」は両方の因子に含まれるかも。
国語 0.826 -0.074
古文 0.653 -0.001
Factor1 Factor2
SS loadings 1.506 1.494
Proportion Var 0.251 0.249
Cumulative Var 0.251 0.500
The degrees of freedom for the model is 4 and the fit was 0.0726
> cfa(r, n, loc) 変数1〜変数3は第2因子，変数4〜変数6は第1因子に含まれるとして検証的因子分析を行う
$loadings 推定された因子負荷量行列
[,1] [,2]
[1,] 0.0000000 0.6418457
[2,] 0.0000000 0.6523311
[3,] 0.0000000 0.7816607
[4,] 0.8629436 0.0000000
[5,] 0.7094359 0.0000000
[6,] 0.6222413 0.0000000
$fac.cor 因子間相関係数行列
[,1] [,2]
[1,] 1.0000000 0.7402071
[2,] 0.7402071 1.0000000
$chisq 適合度検定のカイ二乗値
[1] 10.08075
$P P 値
[1] 0.2594043
$df 自由度
[1] 8
$GFI 以下は各種の適合度指標
[1] 0.9667234
$AGFI
[1] 0.912649
$SRMR
[1] 0.03562844
$RMSEA
[1] 0.05125626
> cfa(r, n, loc2) 変数1〜変数3は第2因子，変数4は両方に，変数5，変数6の変数は第1因子に含まれるとして検証的因子分析を行う
$loadings
[,1] [,2]
[1,] 0.0000000 0.6294310
[2,] 0.0000000 0.6615699
[3,] 0.0000000 0.7826791
[4,] 0.6074606 0.2870043
[5,] 0.7675569 0.0000000
[6,] 0.6591889 0.0000000
$fac.cor
[,1] [,2]
[1,] 1.0000000 0.6306548
[2,] 0.6306548 1.0000000
$chisq
[1] 7.387895
$P
[1] 0.3896392
$df
[1] 7
$GFI
[1] 0.9765228
$AGFI
[1] 0.9295685
$SRMR
[1] 0.02943410
$RMSEA
[1] 0.02365867
直前のページへ戻る E-mail to Shigenobu AOKI
