第11回（1月14日） Task Check and Weekly Assignment
☆課題 psychパッケージに含まれるサンプルデータ，msq(動機づけ状態調査)から， ADAC(活性ー不活性度形容詞チェックリスト)部分を抜き出し，カテゴリカル因子分析を行い， 探索的な因子構造を確認した上で構造方程式モデリングで確認しなさい。 ADAC部分の抜き出しについては，下部のサンプルコードを用いてよい。 また，msqデータに含まれる他の変数を用いても良い。
回答者の反応が「正答・誤答」や「はい・いいえ」のような二値で得られた場合，一般的な(ピアソンの積率)相関係数を求めて分析をすることはできません。
そうしたデータの場合は，新しいテスト理論ともいわれる項目反応理論(Item Response Theory/Latent Trait Model)を用いて分析することになります。
パッケージはltmを使います。
# install.packages('ltm')
library(ltm)
## Loading required package: MASS
## Loading required package: msm
## Loading required package: mvtnorm
## Loading required package: polycor
## Loading required package: sfsmisc
パッケージに含まれるサンプルデータ，LSATを使って次のコードを実行してみましょう
data(LSAT)
result.1 <- rasch(LSAT)
result.1
##
## Call:
## rasch(data = LSAT)
##
## Coefficients:
## Dffclt.Item 1 Dffclt.Item 2 Dffclt.Item 3 Dffclt.Item 4 Dffclt.Item 5
## -3.615 -1.322 -0.318 -1.730 -2.780
## Dscrmn
## 0.755
##
## Log.Lik: -2467
plot(result.1)
分析結果は数値で見ることも可能ですが，plotを見た方がわかりやすいかもしれません。plot(result.1)で示されているのは，項目特性曲線と呼ばれるもので，横軸が回答者の潜在能力（因子得点），縦軸が正答（あるいは「はい」と答える）率を示しています。
項目によってはこの確率を表す関数が右や左にずれていますが，右にずれるとより正答しにくい（「はい」と答えにくい）困難な項目であることになります。この曲線の位置を決めるのが困難度母数と呼ばれるもので，数値出力の際にDifficult.XXXXで示された数字です。
この曲線を変換し，どの値のときにもっともその項目の判別力が高くなるかを示したのが項目情報曲線です。
plot(result.1, type = "IIC")
各項目の項目情報曲線を累積し，テスト全体(ここでは5項目からなるテスト)での判別力を示すのがテスト情報曲線です。
plot(result.1, type = "IIC", items = 0)
項目の特徴を示すのに，困難度母数だけでなく，曲線の傾きに該当する識別力(Discriminant parameter)を加えるモデルもあります。
result.2 <- ltm(LSAT ~ z1)
result.2
##
## Call:
## ltm(formula = LSAT ~ z1)
##
## Coefficients:
## Dffclt Dscrmn
## Item 1 -3.360 0.825
## Item 2 -1.370 0.723
## Item 3 -0.280 0.890
## Item 4 -1.866 0.689
## Item 5 -3.124 0.657
##
## Log.Lik: -2467
plot(result.2)
ここで示されているDscrmnが識別力母数と呼ばれるもので，プロットを見るとカーブの傾きが項目毎に異なっていることがわかります。二つの母数を推定するので，二母数モデルともよばれます。
二母数モデルでも，情報曲線を描くことが可能です。
plot(result.2, type = "IIC")
plot(result.1, type = "IIC", items = 0)
さらに「あて推量母数Guessing parameter」を加えた三母数モデルというのもあります。
result.3 <- tpm(LSAT)
result.3
##
## Call:
## tpm(data = LSAT)
##
## Coefficients:
## Gussng Dffclt Dscrmn
## Item 1 0.037 -3.296 0.829
## Item 2 0.078 -1.145 0.760
## Item 3 0.012 -0.249 0.902
## Item 4 0.035 -1.766 0.701
## Item 5 0.053 -2.990 0.666
##
## Log.Lik: -2467
plot(result.3)
plot(result.3, type = "IIC")
plot(result.3, type = "IIC", items = 0)
数値出力のGussngで表示されているのがそれで，曲線全体が少し上にあがっているのが確認できます。
二値反応ではなく，数段階の反応で回答が得られる場合もあります。複数の順序的反応に対応した項目反応理論のモデルが，段階反応モデルです。
ltmパッケージに含まれるサンプルデータセット，Scienceの一部を使ってみてみましょう。この調査データは，四段階での反応が得られています。
data(Science)
head(Science)
## Comfort Environment Work Future
## 1 strongly agree strongly agree strongly agree agree
## 2 agree strongly agree agree agree
## 3 agree disagree disagree disagree
## 4 agree agree disagree disagree
## 5 agree strongly disagree strongly agree strongly agree
## 6 strongly agree agree strongly agree agree
## Technology Industry Benefit
## 1 strongly agree agree disagree
## 2 agree agree agree
## 3 strongly agree strongly agree agree
## 4 strongly agree strongly agree agree
## 5 disagree agree strongly disagree
## 6 agree strongly agree agree
result.g <- grm(Science[c(1, 3, 4, 7)])
result.g
##
## Call:
## grm(data = Science[c(1, 3, 4, 7)])
##
## Coefficients:
## Extrmt1 Extrmt2 Extrmt3 Dscrmn
## Comfort -4.672 -2.536 1.408 1.041
## Work -2.385 -0.735 1.849 1.226
## Future -2.281 -0.965 0.856 2.299
## Benefit -3.060 -0.906 1.543 1.094
##
## Log.Lik: -1609
plot(result.g)
プロットで表示されるのは，アイテムカテゴリ反応特性曲線Item Response Category Characteristiv Curveと呼ばれる者で，各段階への反応確率が示されています。また，段階反応モデルでも情報曲線を描くことが可能です。
plot(result.g, type = "IIC")
plot(result.g, type = "IIC", items = 0)
ところで，なぜデータの一部しか使わなかったのか気になった方がいるかもしれません。全てのデータを使うとどうなるのでしょうか。
result.g2 <- grm(Science)
result.g2
##
## Call:
## grm(data = Science)
##
## Coefficients:
## Extrmt1 Extrmt2 Extrmt3 Dscrmn
## Comfort -10.768 -5.645 3.097 0.411
## Environment -2.154 -0.790 0.627 1.570
## Work 32.102 9.261 -24.402 -0.074
## Future -30.602 -11.806 10.455 0.108
## Technology -2.462 -0.885 0.642 1.650
## Industry -2.870 -1.529 0.286 1.642
## Benefit -21.232 -5.982 10.297 0.136
##
## Log.Lik: -2998
plot(result.g2)
識別母数に負の数があるものが現れました。また，ICCの形も極端になだらかなおかしなものが含まれています。
これは項目反応理論が一因子(想定する潜在能力は一種類)を仮定したモデルだからであり，Scienceデータは実は複数の因子からなるデータセットなのです。
このことを確認するために，四段階という段階的反応(順序尺度水準)に対応した因子分析，カテゴリカル因子分析をして確認してみます。
下準備として，元のデータセットが因子型で入っているので，これを数値型に戻します。
summary(Science)
## Comfort Environment Work
## strongly disagree: 5 strongly disagree: 29 strongly disagree: 33
## disagree : 32 disagree : 90 disagree : 98
## agree :266 agree :145 agree :206
## strongly agree : 89 strongly agree :128 strongly agree : 55
## Future Technology Industry
## strongly disagree: 14 strongly disagree: 18 strongly disagree: 10
## disagree : 72 disagree : 91 disagree : 47
## agree :210 agree :157 agree :173
## strongly agree : 96 strongly agree :126 strongly agree :162
## Benefit
## strongly disagree: 21
## disagree :100
## agree :193
## strongly agree : 78
str(Science)
## 'data.frame': 392 obs. of 7 variables:
## $ Comfort : Factor w/ 4 levels "strongly disagree",..: 4 3 3 3 3 4 3 3 3 4 ...
## $ Environment: Factor w/ 4 levels "strongly disagree",..: 4 4 2 3 1 3 2 2 3 3 ...
## $ Work : Factor w/ 4 levels "strongly disagree",..: 4 3 2 2 4 4 2 2 3 3 ...
## $ Future : Factor w/ 4 levels "strongly disagree",..: 3 3 2 2 4 3 3 3 4 3 ...
## $ Technology : Factor w/ 4 levels "strongly disagree",..: 4 3 4 4 2 3 4 3 4 3 ...
## $ Industry : Factor w/ 4 levels "strongly disagree",..: 3 3 4 4 3 4 4 4 4 3 ...
## $ Benefit : Factor w/ 4 levels "strongly disagree",..: 2 3 3 3 1 3 4 4 2 3 ...
Science.c <- as.data.frame(data.matrix(Science))
str(Science.c)
## 'data.frame': 392 obs. of 7 variables:
## $ Comfort : int 4 3 3 3 3 4 3 3 3 4 ...
## $ Environment: int 4 4 2 3 1 3 2 2 3 3 ...
## $ Work : int 4 3 2 2 4 4 2 2 3 3 ...
## $ Future : int 3 3 2 2 4 3 3 3 4 3 ...
## $ Technology : int 4 3 4 4 2 3 4 3 4 3 ...
## $ Industry : int 3 3 4 4 3 4 4 4 4 3 ...
## $ Benefit : int 2 3 3 3 1 3 4 4 2 3 ...
psychパッケージのfa関数やfa.parallel関数には，カテゴリカルに対応したfa.poly関数やfa.parallel.poly関数が用意されています。これを使って分析すると，カテゴリカル因子分析を行うことになります。
ちなみに，polyは量的変数のピアソンの積率相関係数に対する，順序尺度水準の相関係数であるポリコリック相関係数の略です。
library(psych)
##
## Attaching package: 'psych'
##
## 以下のオブジェクトはマスクされています (from 'package:ltm') :
##
## factor.scores
##
## 以下のオブジェクトはマスクされています (from 'package:polycor') :
##
## polyserial
fa.parallel.poly(Science.c)
##
##
##
##
## See the graphic output for a description of the results
## Parallel analysis suggests that the number of factors = 3 and the number of components = 2
## Call: fa.parallel.poly(x = Science.c)
## Parallel analysis suggests that the number of factors = 3 and the number of components = 2
##
## Eigen Values of
## Original factors Simulated data Original components simulated data
## 1 1.43 0.64 2.08 1.20
## 2 0.94 0.15 1.88 1.13
## 3 0.18 0.09 0.87 1.07
fa.result <- fa.poly(Science.c, 2, rotate = "promax")
print(fa.result, sort = T, digit = 3, cut = 0.3)
## Factor Analysis using method = minres
## Call: fa.poly(x = Science.c, nfactors = 2, rotate = "promax")
## Standardized loadings (pattern matrix) based upon correlation matrix
## item MR1 MR2 h2 u2 com
## Future 4 0.757 0.569 0.431 1.21
## Work 3 0.558 0.318 0.682 1.02
## Benefit 7 0.529 0.282 0.718 1.13
## Comfort 1 0.495 0.284 0.716 1.00
## Technology 5 0.700 0.486 0.514 1.03
## Environment 2 0.660 0.432 0.568 1.07
## Industry 6 0.631 0.423 0.577 1.00
##
## MR1 MR2
## SS loadings 1.430 1.366
## Proportion Var 0.204 0.195
## Cumulative Var 0.204 0.399
## Proportion Explained 0.511 0.489
## Cumulative Proportion 0.511 1.000
##
## With factor correlations of
## MR1 MR2
## MR1 1.000 0.084
## MR2 0.084 1.000
##
## Mean item complexity = 1.1
## Test of the hypothesis that 2 factors are sufficient.
##
## The degrees of freedom for the null model are 21 and the objective function was 1.248 with Chi Square of 484.1
## The degrees of freedom for the model are 8 and the objective function was 0.075
##
## The root mean square of the residuals (RMSR) is 0.041
## The df corrected root mean square of the residuals is 0.094
##
## The harmonic number of observations is 392 with the empirical chi square 27.92 with prob < 0.000491
## The total number of observations was 392 with MLE Chi Square = 28.99 with prob < 0.000319
##
## Tucker Lewis Index of factoring reliability = 0.8806
## RMSEA index = 0.0826 and the 90 % confidence intervals are 0.0511 0.1149
## BIC = -18.78
## Fit based upon off diagonal values = 0.975
## Measures of factor score adequacy
## MR1 MR2
## Correlation of scores with factors 0.847 0.843
## Multiple R square of scores with factors 0.717 0.711
## Minimum correlation of possible factor scores 0.434 0.422
探索的因子分析ではなく，確認的因子分析(=構造方程式モデリング)をする場合は，推定の際に順序尺度水準で測定された変数があることをオプションで指定します。
library(lavaan)
## This is lavaan 0.5-15
## lavaan is BETA software! Please report any bugs.
model <- "
fa1 =~ Future+Work+Benefit+Comfort
fa2 =~ Technology+Environment+Industry
fa1 ~~ 0 * fa2
"
result.cfa <- cfa(model, Science.c, ordered = c("Future", "Work", "Benefit",
"Comfort", "Technology", "Environment", "Industry"))
summary(result.cfa, standardized = T, fit.measures = T)
## lavaan (0.5-15) converged normally after 20 iterations
##
## Number of observations 392
##
## Estimator DWLS Robust
## Minimum Function Test Statistic 47.563 41.383
## Degrees of freedom 14 14
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.254
## Shift parameter 3.459
## for simple second-order correction (Mplus variant)
##
## Model test baseline model:
##
## Minimum Function Test Statistic 683.313 530.404
## Degrees of freedom 21 21
## P-value 0.000 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.949 0.946
## Tucker-Lewis Index (TLI) 0.924 0.919
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.078 0.071
## 90 Percent Confidence Interval 0.055 0.103 0.046 0.096
## P-value RMSEA <= 0.05 0.027 0.078
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Robust.sem
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Latent variables:
## fa1 =~
## Future 1.000 0.777 0.777
## Work 0.677 0.080 8.455 0.000 0.526 0.526
## Benefit 0.683 0.084 8.161 0.000 0.531 0.531
## Comfort 0.663 0.087 7.649 0.000 0.515 0.515
## fa2 =~
## Technology 1.000 0.701 0.701
## Environment 0.945 0.109 8.657 0.000 0.662 0.662
## Industry 0.886 0.099 8.928 0.000 0.621 0.621
##
## Covariances:
## fa1 ~~
## fa2 0.000 0.000 0.000
##
## Intercepts:
## fa1 0.000 0.000 0.000
## fa2 0.000 0.000 0.000
##
## Thresholds:
## Future|t1 -1.803 0.119 -15.091 0.000 -1.803 -1.803
## Future|t2 -0.774 0.071 -10.937 0.000 -0.774 -0.774
## Future|t3 0.691 0.069 9.981 0.000 0.691 0.691
## Work|t1 -1.377 0.091 -15.155 0.000 -1.377 -1.377
## Work|t2 -0.428 0.066 -6.536 0.000 -0.428 -0.428
## Work|t3 1.079 0.079 13.693 0.000 1.079 1.079
## Benefit|t1 -1.611 0.105 -15.415 0.000 -1.611 -1.611
## Benefit|t2 -0.500 0.066 -7.531 0.000 -0.500 -0.500
## Benefit|t3 0.845 0.072 11.685 0.000 0.845 0.845
## Comfort|t1 -2.234 0.172 -12.960 0.000 -2.234 -2.234
## Comfort|t2 -1.314 0.088 -14.952 0.000 -1.314 -1.314
## Comfort|t3 0.749 0.070 10.652 0.000 0.749 0.749
## Technology|t1 -1.686 0.110 -15.343 0.000 -1.686 -1.686
## Technology|t2 -0.589 0.068 -8.715 0.000 -0.589 -0.589
## Technology|t3 0.464 0.066 7.034 0.000 0.464 0.464
## Environmnt|t1 -1.447 0.094 -15.311 0.000 -1.447 -1.447
## Environmnt|t2 -0.514 0.067 -7.729 0.000 -0.514 -0.514
## Environmnt|t3 0.450 0.066 6.835 0.000 0.450 0.450
## Industry|t1 -1.951 0.134 -14.547 0.000 -1.951 -1.951
## Industry|t2 -1.056 0.078 -13.531 0.000 -1.056 -1.056
## Industry|t3 0.219 0.064 3.428 0.001 0.219 0.219
##
## Variances:
## Future 0.396 0.396 0.396
## Work 0.723 0.723 0.723
## Benefit 0.718 0.718 0.718
## Comfort 0.734 0.734 0.734
## Technology 0.509 0.509 0.509
## Environment 0.562 0.562 0.562
## Industry 0.615 0.615 0.615
## fa1 0.604 0.074 1.000 1.000
## fa2 0.491 0.067 1.000 1.000
このように，反応段階に応じたモデリングが可能ですから，実際に自分で研究を計画する際は回答者が反応しやすい段階，回答方法を考えて，より丁寧にデータを取るように心がけましょう。
追伸） 本日の課題で用いるデータセットを作るには，次のコードを使用してください。
data(msq)
msq.sub <- subset(msq, select = c("active", "energetic", "vigorous", "wakeful",
"full.of.pep", "lively", "sleepy", "tired", "drowsy"))