検証的因子分析     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

Made with Macintosh