1 因子分析の準備

dat<-read.csv("data/data_ch9-1.csv",header=T,row.names=1)
head(dat)

	item1	item2	item3	item4	item5	item6	item7	item8	item9
s001	3	3	3	3	3	3	3	3	3
s002	2	1	4	4	3	4	4	4	2
s003	5	5	4	5	3	3	3	5	4
s004	5	5	5	5	2	5	5	5	5
s005	5	5	3	5	4	5	4	5	4
s006	5	4	2	5	2	3	5	5	3

1.1 サンプルサイズ

KMO（Kaiser-Meyer-Olkin factor adequacy）サンプリングの適切性指標（0.7以上で適切）

KMO(dat)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = dat)
## Overall MSA =  0.76
## MSA for each item = 
## item1 item2 item3 item4 item5 item6 item7 item8 item9 
##  0.79  0.76  0.61  0.86  0.83  0.86  0.70  0.76  0.74

1.2 観測変数間の相関

Bartlettの球面性検定(共分散行列が球形性を持たないという帰無仮説の検定)
有意であれば変数間に相関があり因子分析を行うには妥当(帰無仮説は無相間)

注意. Bartlett の球形検定と Mauchly の球形検定の違い.

Bartlettは相関行列の無相関性の検定.
Mauchlyは被験者内要因の分散が等しい検定.
（被験者間計画の分散分析でいう分散の等質性の仮定に相当する） http://riseki.php.xdomain.jp/index.php?ANOVA君%2F球面性検定の出力

相関が弱ければ共通因子がないことになり、因子分析は意味をなさなくなる。（相関0.3以上の項目がなければ因子分析には適さない）.
相関が0.9以上のものは多重共線性（同一内容測定）が考えられ、どちらかの変数を除外する。

options(digits=2)　#桁数指定
cor(dat)

##       item1 item2  item3 item4 item5 item6 item7 item8  item9
## item1  1.00  0.64  0.169  0.37 0.212  0.33 0.228  0.42  0.377
## item2  0.64  1.00  0.050  0.33 0.160  0.34 0.153  0.45  0.472
## item3  0.17  0.05  1.000  0.22 0.048  0.17 0.469  0.12 -0.054
## item4  0.37  0.33  0.221  1.00 0.230  0.20 0.274  0.50  0.346
## item5  0.21  0.16  0.048  0.23 1.000  0.13 0.047  0.22  0.275
## item6  0.33  0.34  0.165  0.20 0.134  1.00 0.174  0.16  0.154
## item7  0.23  0.15  0.469  0.27 0.047  0.17 1.000  0.24  0.070
## item8  0.42  0.45  0.123  0.50 0.224  0.16 0.239  1.00  0.644
## item9  0.38  0.47 -0.054  0.35 0.275  0.15 0.070  0.64  1.000

#n=nrow(dat);n
cortest.bartlett(cor(dat),n=nrow(dat))

## $chisq
## [1] 477
## 
## $p.value
## [1] 2.4e-78
## 
## $df
## [1] 36

pairs.panels(dat,lm=TRUE,density=FALSE)

帰無仮説「観測変数は無相関」は棄却される。

2 探索的因子分析

2.1 因子数の決定

r.eigen<-eigen(cor(dat))
r.eigen

## eigen() decomposition
## $values
## [1] 3.25 1.45 1.00 0.89 0.65 0.60 0.51 0.34 0.30
## 
## $vectors
##        [,1]   [,2]  [,3]   [,4]  [,5]   [,6]    [,7]   [,8]    [,9]
##  [1,] -0.42 -0.025  0.33 -0.056 -0.46 -0.275  0.0032  0.587  0.2953
##  [2,] -0.41 -0.167  0.37 -0.166 -0.33 -0.070  0.0230 -0.634 -0.3495
##  [3,] -0.16  0.660 -0.12  0.056 -0.14  0.027 -0.6962 -0.126  0.0351
##  [4,] -0.37  0.087 -0.29 -0.024  0.44 -0.714  0.1390 -0.169  0.1432
##  [5,] -0.22 -0.145 -0.30  0.870 -0.25  0.040  0.0864 -0.034 -0.1063
##  [6,] -0.26  0.143  0.63  0.338  0.59  0.224 -0.0041  0.075 -0.0073
##  [7,] -0.23  0.579 -0.15 -0.125 -0.12  0.332  0.6652 -0.039  0.0574
##  [8,] -0.42 -0.143 -0.32 -0.242  0.19  0.212 -0.1359  0.402 -0.6143
##  [9,] -0.38 -0.362 -0.22 -0.141  0.10  0.449 -0.1644 -0.196  0.6139

print(r.eigen$values,digit=2)#固有値

## [1] 3.25 1.45 1.00 0.89 0.65 0.60 0.51 0.34 0.30

2.2 スクリープロット

plot(r.eigen$values,type="b")
abline(h=1)

2.3 MAP基準の確認

VSS(dat, n = nrow(dat), fm = "ml")

## 
## Very Simple Structure
## Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm, 
##     n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
## VSS complexity 1 achieves a maximimum of 0.69  with  2  factors
## VSS complexity 2 achieves a maximimum of 0.81  with  4  factors
## 
## The Velicer MAP achieves a minimum of 0.04  with  1  factors 
## BIC achieves a minimum of  -55  with  3  factors
## Sample Size adjusted BIC achieves a minimum of  -17  with  3  factors
## 
## Statistics by number of factors 
##   vss1 vss2   map dof   chisq    prob sqresid  fit RMSEA   BIC SABIC complex
## 1 0.64 0.00 0.043  27 1.4e+02 5.8e-17     5.6 0.64  0.14  -4.4  81.1     1.0
## 2 0.69 0.76 0.056  19 6.8e+01 1.9e-07     3.7 0.76  0.11 -32.5  27.7     1.3
## 3 0.57 0.78 0.074  12 8.4e+00 7.5e-01     2.9 0.82  0.00 -55.2 -17.2     1.4
## 4 0.60 0.81 0.130   6 1.1e+00 9.8e-01     2.0 0.87  0.00 -30.6 -11.6     1.5
## 5 0.50 0.69 0.180   1 3.6e-02 8.5e-01     2.3 0.85  0.00  -5.3  -2.1     1.7
## 6 0.50 0.70 0.267  -3 8.5e-07      NA     2.0 0.87    NA    NA    NA     2.0
## 7 0.43 0.64 0.471  -6 7.2e-10      NA     1.9 0.88    NA    NA    NA     2.0
## 8 0.44 0.64 1.000  -8 0.0e+00      NA     1.7 0.89    NA    NA    NA     2.1
## 9 0.53 0.75    NA  -9 2.1e+01      NA     3.2 0.80    NA    NA    NA     1.6
##    eChisq    SRMR  eCRMS  eBIC
## 1 1.6e+02 1.0e-01 0.1211  15.3
## 2 4.9e+01 5.8e-02 0.0802 -51.8
## 3 6.4e+00 2.1e-02 0.0364 -57.2
## 4 5.7e-01 6.3e-03 0.0154 -31.2
## 5 2.2e-02 1.2e-03 0.0074  -5.3
## 6 5.1e-07 5.9e-06     NA    NA
## 7 8.2e-10 2.4e-07     NA    NA
## 8 1.6e-14 1.1e-09     NA    NA
## 9 1.4e+01 3.1e-02     NA    NA

MAP 基準では因子数 1、BIC 基準では因子数 3 が提案されています。 MAP 基準は少なめの因子数を、BIC基準や平行分析は多めの因子数を提案する性質があります。

平行分析や MAP/BIC 基準を利用することが現在は推奨されており、カイザー基準やスクリーテストにはあまり頼らない方が良いと考えられます。

2.4 平行分析

fa.parallel(cor(dat), fm = "ml", n.obs = nrow(dat), n.iter = 100)

## Parallel analysis suggests that the number of factors =  3  and the number of components =  2

平行分析は因子=3、および主成分components=2を提案する。
（平行分析は、因子分析にも対応しているため、両方の結果を提案してくれる）

2.5 最尤法・オブリミン回転

r.fa <- fa(dat, nfactors = 3, fm = "ml", rotate = "oblimin")

## Loading required namespace: GPArotation

print(r.fa, sort = TRUE, digits = 2)

## Factor Analysis using method =  ml
## Call: fa(r = dat, nfactors = 3, rotate = "oblimin", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##       item   ML1   ML3   ML2   h2   u2 com
## item8    8  0.89 -0.04  0.08 0.77 0.23 1.0
## item9    9  0.70  0.13 -0.18 0.60 0.40 1.2
## item4    4  0.43  0.11  0.25 0.36 0.64 1.8
## item5    5  0.24  0.09  0.01 0.09 0.91 1.3
## item2    2  0.03  0.84 -0.08 0.71 0.29 1.0
## item1    1  0.01  0.75  0.11 0.61 0.39 1.0
## item6    6 -0.11  0.46  0.16 0.21 0.79 1.4
## item3    3 -0.03  0.01  0.73 0.53 0.47 1.0
## item7    7  0.10  0.05  0.62 0.43 0.57 1.1
## 
##                        ML1  ML3  ML2
## SS loadings           1.63 1.59 1.09
## Proportion Var        0.18 0.18 0.12
## Cumulative Var        0.18 0.36 0.48
## Proportion Explained  0.38 0.37 0.25
## Cumulative Proportion 0.38 0.75 1.00
## 
##  With factor correlations of 
##      ML1  ML3  ML2
## ML1 1.00 0.61 0.15
## ML3 0.61 1.00 0.19
## ML2 0.15 0.19 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 3 factors are sufficient.
## 
## The degrees of freedom for the null model are  36  and the objective function was  2.4 with Chi Square of  477
## The degrees of freedom for the model are 12  and the objective function was  0.04 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  0.04 
## 
## The harmonic number of observations is  200 with the empirical chi square  6.4  with prob <  0.9 
## The total number of observations was  200  with Likelihood Chi Square =  8.4  with prob <  0.75 
## 
## Tucker Lewis Index of factoring reliability =  1
## RMSEA index =  0  and the 90 % confidence intervals are  0 0.052
## BIC =  -55
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1  ML3  ML2
## Correlation of (regression) scores with factors   0.92 0.91 0.82
## Multiple R square of scores with factors          0.85 0.83 0.68
## Minimum correlation of possible factor scores     0.70 0.65 0.36

print( r.fa$loadings, digits = 2, cutoff = 0.3 , sort = TRUE)

## 
## Loadings:
##       ML1   ML3   ML2  
## item8  0.89            
## item9  0.70            
## item1        0.75      
## item2        0.84      
## item3              0.73
## item7              0.62
## item4  0.43            
## item5                  
## item6        0.46      
## 
##                 ML1  ML3  ML2
## SS loadings    1.55 1.50 1.06
## Proportion Var 0.17 0.17 0.12
## Cumulative Var 0.17 0.34 0.46

biplot(r.fa)

fa.diagram( r.fa )

3 確認的因子分析

3.1 lavaanパッケージの読み込み

library("lavaan")

## This is lavaan 0.6-15
## lavaan is FREE software! Please report any bugs.

## 
## Attaching package: 'lavaan'

## The following object is masked from 'package:psych':
## 
##     cor2cov

3.2 3因子構造のモデルを記述

model.1 <- '
  LV.1 =~ item8 + item9 + item4
  LV.2 =~ item3 + item7
  LV.3 =~ item2 + item1 + item6 '

3.3 確認的因子分析

fit.1 <- cfa(model.1, data = dat, estimator = "ML")
summary(fit.1, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6.15 ended normally after 38 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        19
## 
##   Number of observations                           200
## 
## Model Test User Model:
##                                                       
##   Test statistic                                38.362
##   Degrees of freedom                                17
##   P-value (Chi-square)                           0.002
## 
## Model Test Baseline Model:
## 
##   Test statistic                               465.296
##   Degrees of freedom                                28
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.951
##   Tucker-Lewis Index (TLI)                       0.920
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -2148.404
##   Loglikelihood unrestricted model (H1)      -2129.223
##                                                       
##   Akaike (AIC)                                4334.809
##   Bayesian (BIC)                              4397.477
##   Sample-size adjusted Bayesian (SABIC)       4337.283
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.079
##   90 Percent confidence interval - lower         0.046
##   90 Percent confidence interval - upper         0.113
##   P-value H_0: RMSEA <= 0.050                    0.072
##   P-value H_0: RMSEA >= 0.080                    0.519
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.061
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   LV.1 =~                                                               
##     item8             1.000                               0.765    0.868
##     item9             0.950    0.103    9.244    0.000    0.727    0.728
##     item4             0.881    0.119    7.430    0.000    0.674    0.565
##   LV.2 =~                                                               
##     item3             1.000                               0.564    0.491
##     item7             1.719    0.819    2.099    0.036    0.969    0.955
##   LV.3 =~                                                               
##     item2             1.000                               0.913    0.819
##     item1             0.937    0.105    8.938    0.000    0.855    0.789
##     item6             0.477    0.091    5.234    0.000    0.436    0.408
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   LV.1 ~~                                                               
##     LV.2              0.114    0.064    1.781    0.075    0.264    0.264
##     LV.3              0.459    0.075    6.148    0.000    0.657    0.657
##   LV.2 ~~                                                               
##     LV.3              0.131    0.075    1.749    0.080    0.255    0.255
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .item8             0.192    0.052    3.725    0.000    0.192    0.247
##    .item9             0.469    0.065    7.204    0.000    0.469    0.470
##    .item4             0.970    0.108    9.005    0.000    0.970    0.681
##    .item3             0.998    0.177    5.631    0.000    0.998    0.758
##    .item7             0.092    0.433    0.211    0.833    0.092    0.089
##    .item2             0.409    0.087    4.714    0.000    0.409    0.329
##    .item1             0.444    0.081    5.498    0.000    0.444    0.378
##    .item6             0.949    0.099    9.544    0.000    0.949    0.833
##     LV.1              0.586    0.089    6.555    0.000    1.000    1.000
##     LV.2              0.318    0.170    1.873    0.061    1.000    1.000
##     LV.3              0.834    0.140    5.950    0.000    1.000    1.000

3.4 適合度指標の検討

fitMeasures(fit.1)

##                  npar                  fmin                 chisq 
##               1.9e+01               9.6e-02               3.8e+01 
##                    df                pvalue        baseline.chisq 
##               1.7e+01               2.0e-03               4.7e+02 
##           baseline.df       baseline.pvalue                   cfi 
##               2.8e+01               0.0e+00               9.5e-01 
##                   tli                  nnfi                   rfi 
##               9.2e-01               9.2e-01               8.6e-01 
##                   nfi                  pnfi                   ifi 
##               9.2e-01               5.6e-01               9.5e-01 
##                   rni                  logl     unrestricted.logl 
##               9.5e-01              -2.1e+03              -2.1e+03 
##                   aic                   bic                ntotal 
##               4.3e+03               4.4e+03               2.0e+02 
##                  bic2                 rmsea        rmsea.ci.lower 
##               4.3e+03               7.9e-02               4.6e-02 
##        rmsea.ci.upper        rmsea.ci.level          rmsea.pvalue 
##               1.1e-01               9.0e-01               7.2e-02 
##        rmsea.close.h0 rmsea.notclose.pvalue     rmsea.notclose.h0 
##               5.0e-02               5.2e-01               8.0e-02 
##                   rmr            rmr_nomean                  srmr 
##               7.3e-02               7.3e-02               6.1e-02 
##          srmr_bentler   srmr_bentler_nomean                  crmr 
##               6.1e-02               6.1e-02               6.9e-02 
##           crmr_nomean            srmr_mplus     srmr_mplus_nomean 
##               6.9e-02               6.1e-02               6.1e-02 
##                 cn_05                 cn_01                   gfi 
##               1.4e+02               1.8e+02               9.5e-01 
##                  agfi                  pgfi                   mfi 
##               8.9e-01               4.5e-01               9.5e-01 
##                  ecvi 
##               3.8e-01

3.5 2因子構造のモデルを記述

model.2 <- '
  LV.1 =~ item8 + item9 + item4
  LV.2 =~ item2 + item1 + item6 '

3.6 確認的因子分析

fit.2 <- cfa(model.2, data = dat, estimator = "ML")
summary(fit.2, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6.15 ended normally after 26 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        13
## 
##   Number of observations                           200
## 
## Model Test User Model:
##                                                       
##   Test statistic                                13.786
##   Degrees of freedom                                 8
##   P-value (Chi-square)                           0.088
## 
## Model Test Baseline Model:
## 
##   Test statistic                               377.790
##   Degrees of freedom                                15
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.984
##   Tucker-Lewis Index (TLI)                       0.970
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1581.788
##   Loglikelihood unrestricted model (H1)      -1574.895
##                                                       
##   Akaike (AIC)                                3189.576
##   Bayesian (BIC)                              3232.455
##   Sample-size adjusted Bayesian (SABIC)       3191.269
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.060
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.112
##   P-value H_0: RMSEA <= 0.050                    0.327
##   P-value H_0: RMSEA >= 0.080                    0.305
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.036
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   LV.1 =~                                                               
##     item8             1.000                               0.756    0.857
##     item9             0.980    0.106    9.261    0.000    0.741    0.742
##     item4             0.881    0.121    7.296    0.000    0.666    0.558
##   LV.2 =~                                                               
##     item2             1.000                               0.929    0.834
##     item1             0.905    0.102    8.846    0.000    0.842    0.776
##     item6             0.463    0.089    5.169    0.000    0.430    0.403
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   LV.1 ~~                                                               
##     LV.2              0.465    0.075    6.197    0.000    0.662    0.662
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .item8             0.207    0.052    3.987    0.000    0.207    0.266
##    .item9             0.449    0.065    6.887    0.000    0.449    0.450
##    .item4             0.982    0.109    9.021    0.000    0.982    0.689
##    .item2             0.379    0.089    4.273    0.000    0.379    0.305
##    .item1             0.468    0.081    5.762    0.000    0.468    0.398
##    .item6             0.954    0.100    9.564    0.000    0.954    0.838
##     LV.1              0.571    0.089    6.434    0.000    1.000    1.000
##     LV.2              0.864    0.143    6.041    0.000    1.000    1.000

3.7 model1とmodel2の比較

anova(fit.1, fit.2)

## Warning in lavTestLRT(object = object, ..., model.names = NAMES): lavaan
## WARNING: some models are based on a different set of observed variables

	Df	AIC	BIC	Chisq	Chisq diff	RMSEA	Df diff	Pr(>Chisq)
fit.2	8	3190	3232	14	NA	NA	NA	NA
fit.1	17	4335	4397	38	25	0.09	9	0

3.8 因子得点の計算

head(lavPredict(fit.2))

##       LV.1  LV.2
## [1,] -0.80 -0.92
## [2,] -0.62 -1.88
## [3,]  0.72  0.83
## [4,]  0.96  1.04
## [5,]  0.75  0.98
## [6,]  0.45  0.37

3.9 尺度得点の計算

#str(dat)
LV.1 <- rowMeans(dat[, c(8, 9, 4)])#item8,9,4の平均
LV.2 <- rowMeans(dat[, c(2, 1, 6)])#item2,1,6の平均
head(cbind(LV.1, LV.2))

##      LV.1 LV.2
## s001  3.0  3.0
## s002  3.3  2.3
## s003  4.7  4.3
## s004  5.0  5.0
## s005  4.7  5.0
## s006  4.3  4.0

par(mfrow = c(1, 2))
plot(lavPredict(fit.2))
plot(cbind(LV.1, LV.2))

4 項目反応理論

dat.A<-read.csv("data/data_ch9-2.csv",header=T,row.names=1)
head(dat.A)

	CI1	CI2	CI3	CI4	CI5	CI6	CI7	CI8	CI9	CI10	CI11	CI12	CI13	CI14	CI15	CI16	CI17	CI18	CI19	CI20	A1	A2	A3	A4	A5	A6	A7	A8	A9	A10	A11	A12	A13	A14	A15	A16	A17	A18	A19	A20	A21	A22	A23	A24	A25	A26	A27	A28	A29	A30
S001	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	0	1	0	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1	1
S002	1	1	1	1	1	1	1	1	1	0	0	1	0	0	1	0	0	1	0	1	1	0	0	0	0	0	1	0	0	1	1	1	0	0	1	1	0	1	1	0	1	1	0	0	0	1	1	1	1	1
S003	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	0
S004	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1	0	0	1	0	1	1	0	1	0	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1
S005	0	1	1	0	1	1	1	1	1	1	1	0	0	1	1	0	0	1	0	1	1	0	1	0	1	1	1	0	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1
S006	1	1	0	1	1	1	0	1	1	1	0	1	0	0	1	1	0	1	0	1	1	0	0	0	1	0	1	0	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1

dat.B<-read.csv("data/data_ch9-3.csv",header=T,row.names=1)
head(dat.B)

	CI1	CI2	CI3	CI4	CI5	CI6	CI7	CI8	CI9	CI10	CI11	CI12	CI13	CI14	CI15	CI16	CI17	CI18	CI20	B1	B2	B3	B4	B5	B6	B7	B8	B9	B10	B11	B12	B13	B14	B15	B16	B17	B18	B19	B20	B21	B22	B23	B24	B25	B26	B27	B29	B30
S171	1	1	1	1	1	1	1	1	0	1	0	1	1	1	1	1	0	1	1	1	0	1	1	1	1	1	1	1	1	1	1	0	1	1	1	0	1	1	1	1	1	0	1	1	0	0	0	0
S172	1	1	0	1	1	1	0	1	1	0	1	1	1	1	1	1	0	1	1	1	0	0	0	1	0	0	1	1	1	1	1	1	1	1	1	1	0	1	1	1	0	0	0	1	0	0	0	0
S173	1	1	0	1	1	1	0	1	1	0	0	0	1	1	1	0	0	1	1	1	1	0	1	0	1	1	0	1	0	1	1	1	0	1	1	1	1	0	1	0	0	0	0	0	0	0	0	0
S174	1	1	1	1	1	1	1	1	1	1	1	0	1	1	1	1	0	1	1	1	1	1	1	0	0	0	1	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	0	1	1	0	0	0
S175	1	1	1	1	1	1	1	1	1	0	1	1	0	0	0	1	1	1	1	1	0	1	1	1	0	0	1	1	1	1	1	0	1	1	1	0	1	1	1	1	1	0	1	1	1	1	1	1
S176	0	1	0	1	1	1	0	1	1	0	0	1	0	1	1	0	0	1	1	0	0	0	0	0	0	1	1	1	0	1	1	0	1	1	1	0	0	0	1	1	1	0	1	0	1	1	0	0

4.1 1パラメターモデルで分析

test.A <- rasch(dat.A)
test.B <- rasch(dat.B)
#test.A
#test.B

# 2つの図を1つにまとめて表示する設定
par(mfrow = c(1, 2))
# 項目特性曲線の作成
plot(test.A, type = "ICC", items = c(23, 24, 29))
# 項目情報量曲線の作成
plot(test.A, type = "IIC", items = c(23, 24, 29))

4.2 能力推定値の計算

theta.A <- factor.scores.rasch(test.A, resp.pattern = dat.A)
theta.B <- factor.scores.rasch(test.B, resp.pattern = dat.B)
A.theta <- theta.A$score.dat$z1
B.theta <- theta.B$score.dat$z1

4.3 計算結果の冒頭を確認

head(A.theta)

## [1]  1.217 -0.980  1.411  0.857  0.223 -0.064

head(B.theta)

## [1]  0.31 -0.20 -0.59  0.41  0.41 -0.66

4.4 図9.6

par(family = "HiraKakuProN-W3") #日本語フォントの指定
raw　<-　rowSums(dat.A[, -1])
head(cbind(raw, A.theta))

##      raw A.theta
## S001  43   1.217
## S002  28  -0.980
## S003  44   1.411
## S004  41   0.857
## S005  38   0.223
## S006  35  -0.064

par(ps = 20, mai = c(1, 1, 1, 1), mfrow = c(1, 2))
hist(A.theta, main = "", xlab = "能力推定値", ylab = "度数")
abline(v = mean(A.theta), lty = 2)
plot(raw, A.theta, pch = 1, xlab = "素点", ylab = "能力推定値", cex = 1.5)
abline(lm(A.theta ~ raw), col = "black", lty = 2)

4.5 等化

4.5.1 共通項目の難易度の平均値

(dffcltA.mean <- mean(test.A$coefficients[1:20, 1]))

## [1] 1.3

(dffcltB.mean <- mean(test.B$coefficients[1:20, 1]))

## [1] 0.95

4.5.2 変換式の作成

Intercept <- dffcltA.mean - dffcltB.mean

4.5.3 テストB受験者の能力推定値をテストAに合わせる処理

B.adjusted <- B.theta + Intercept
# 等化前と等化後のテストＢ受験者の能力推定値の変化を箱ひげ図で可視化

par(family = "HiraKakuProN-W3") #日本語フォントの指定
boxplot(A.theta, B.theta, B.adjusted, names = c("テストA", "テストB（等化前）", "テストB（等化後）"), main = NA, xlab = NA, ylab = "能力推定値", col = "grey")

4.6 段階反応モデル

dat.C <- read.csv("data/data_ch9-4.csv",header=T,row.names=1)
head(dat.C)

	Q1	Q2	Q3	Q4	Q5	Q6	Q7	Q8	Q9	Q10	Q11	Q12	Q13	Q14	Q15	Q16	Q17	Q18
S001	1	1	1	1	1	1	2	1	1	1	1	1	1	1	1	2	1	2
S002	2	2	2	2	1	2	2	1	2	2	2	2	3	2	3	2	1	2
S003	NA	NA	NA	NA	NA	NA	2	3	2	3	2	2	2	2	2	3	2	2
S004	2	3	2	2	1	2	NA	NA	NA	2	2	2	2	2	2	2	2	1
S005	2	1	2	3	1	3	2	1	1	2	2	2	2	2	2	0	0	0
S006	0	0	0	2	1	2	2	2	2	3	2	2	1	1	2	2	2	2

4.6.1 パッケージの読み込み

## Loading required package: stats4

## Loading required package: lattice

## 
## Attaching package: 'mirt'

## The following object is masked from 'package:ltm':
## 
##     Science

4.6.2 段階反応モデルで分析

test.C <- mirt(data = dat.C, model = 1, itemtype = "graded")

## 
Iteration: 1, Log-Lik: -1845.020, Max-Change: 2.41028
Iteration: 2, Log-Lik: -1783.668, Max-Change: 0.42607
Iteration: 3, Log-Lik: -1778.501, Max-Change: 0.17289
Iteration: 4, Log-Lik: -1776.516, Max-Change: 0.10180
Iteration: 5, Log-Lik: -1775.439, Max-Change: 0.07488
Iteration: 6, Log-Lik: -1774.792, Max-Change: 0.05800
Iteration: 7, Log-Lik: -1773.827, Max-Change: 0.02622
Iteration: 8, Log-Lik: -1773.714, Max-Change: 0.02055
Iteration: 9, Log-Lik: -1773.633, Max-Change: 0.01707
Iteration: 10, Log-Lik: -1773.490, Max-Change: 0.02325
Iteration: 11, Log-Lik: -1773.458, Max-Change: 0.01041
Iteration: 12, Log-Lik: -1773.440, Max-Change: 0.01317
Iteration: 13, Log-Lik: -1773.406, Max-Change: 0.00982
Iteration: 14, Log-Lik: -1773.397, Max-Change: 0.00623
Iteration: 15, Log-Lik: -1773.391, Max-Change: 0.00546
Iteration: 16, Log-Lik: -1773.383, Max-Change: 0.00689
Iteration: 17, Log-Lik: -1773.379, Max-Change: 0.00432
Iteration: 18, Log-Lik: -1773.377, Max-Change: 0.00303
Iteration: 19, Log-Lik: -1773.373, Max-Change: 0.00283
Iteration: 20, Log-Lik: -1773.372, Max-Change: 0.00274
Iteration: 21, Log-Lik: -1773.371, Max-Change: 0.00233
Iteration: 22, Log-Lik: -1773.370, Max-Change: 0.00157
Iteration: 23, Log-Lik: -1773.370, Max-Change: 0.00134
Iteration: 24, Log-Lik: -1773.369, Max-Change: 0.00222
Iteration: 25, Log-Lik: -1773.369, Max-Change: 0.00225
Iteration: 26, Log-Lik: -1773.368, Max-Change: 0.00107
Iteration: 27, Log-Lik: -1773.368, Max-Change: 0.00090
Iteration: 28, Log-Lik: -1773.368, Max-Change: 0.00208
Iteration: 29, Log-Lik: -1773.367, Max-Change: 0.00036
Iteration: 30, Log-Lik: -1773.367, Max-Change: 0.00020
Iteration: 31, Log-Lik: -1773.367, Max-Change: 0.00017
Iteration: 32, Log-Lik: -1773.367, Max-Change: 0.00079
Iteration: 33, Log-Lik: -1773.367, Max-Change: 0.00073
Iteration: 34, Log-Lik: -1773.367, Max-Change: 0.00013
Iteration: 35, Log-Lik: -1773.367, Max-Change: 0.00060
Iteration: 36, Log-Lik: -1773.367, Max-Change: 0.00020
Iteration: 37, Log-Lik: -1773.367, Max-Change: 0.00012
Iteration: 38, Log-Lik: -1773.367, Max-Change: 0.00054
Iteration: 39, Log-Lik: -1773.367, Max-Change: 0.00023
Iteration: 40, Log-Lik: -1773.367, Max-Change: 0.00012
Iteration: 41, Log-Lik: -1773.367, Max-Change: 0.00047
Iteration: 42, Log-Lik: -1773.367, Max-Change: 0.00033
Iteration: 43, Log-Lik: -1773.367, Max-Change: 0.00013
Iteration: 44, Log-Lik: -1773.367, Max-Change: 0.00008

4.6.3 能力推定値の計算

C.theta <- fscores(test.C)
C.theta

##             F1
##   [1,] -1.4292
##   [2,]  0.4212
##   [3,]  1.0227
##   [4,]  0.5310
##   [5,] -0.0790
##   [6,]  0.0083
##   [7,] -1.5416
##   [8,] -0.5366
##   [9,] -1.2110
##  [10,]  0.6645
##  [11,]  0.7080
##  [12,] -0.7024
##  [13,]  0.1022
##  [14,] -0.8055
##  [15,] -0.4268
##  [16,] -0.5757
##  [17,]  0.8797
##  [18,]  0.7214
##  [19,]  1.3447
##  [20,] -1.4733
##  [21,]  0.2118
##  [22,]  0.0967
##  [23,] -0.0069
##  [24,]  0.8669
##  [25,]  0.2725
##  [26,]  0.8798
##  [27,] -1.2056
##  [28,]  0.6687
##  [29,] -1.9780
##  [30,] -0.8400
##  [31,]  0.1953
##  [32,] -0.2778
##  [33,] -0.7077
##  [34,] -0.0617
##  [35,] -0.2004
##  [36,] -1.5622
##  [37,]  0.1735
##  [38,] -1.5454
##  [39,]  0.9616
##  [40,] -0.9544
##  [41,] -2.1809
##  [42,]  1.4465
##  [43,] -0.3944
##  [44,] -2.1097
##  [45,]  0.1969
##  [46,] -0.2604
##  [47,]  0.9203
##  [48,] -0.4965
##  [49,] -0.0478
##  [50,] -1.0390
##  [51,] -0.2361
##  [52,] -0.9947
##  [53,]  0.5227
##  [54,] -1.1788
##  [55,] -2.4309
##  [56,] -0.6462
##  [57,] -2.2711
##  [58,] -0.4693
##  [59,]  0.4170
##  [60,] -0.5319
##  [61,] -0.4861
##  [62,] -0.3668
##  [63,] -0.6595
##  [64,] -0.8457
##  [65,] -0.4757
##  [66,]  0.5896
##  [67,]  0.1182
##  [68,] -1.5386
##  [69,]  0.6905
##  [70,]  0.6390
##  [71,]  0.1547
##  [72,] -0.3926
##  [73,] -0.2049
##  [74,] -0.6998
##  [75,]  0.8377
##  [76,]  2.2189
##  [77,]  0.2006
##  [78,]  0.6335
##  [79,]  0.4060
##  [80,] -1.4073
##  [81,]  0.0896
##  [82,] -0.4436
##  [83,]  0.8118
##  [84,]  1.4186
##  [85,]  0.4793
##  [86,]  0.3356
##  [87,]  1.2159
##  [88,]  0.7306
##  [89,]  0.5174
##  [90,]  0.3321
##  [91,] -0.8510
##  [92,] -1.1131
##  [93,]  1.2468
##  [94,]  1.3477
##  [95,]  0.0756
##  [96,]  1.3897
##  [97,]  1.7998
##  [98,]  0.2048
##  [99,]  1.2631
## [100,]  1.7926
## [101,]  0.0483
## [102,]  1.5431
## [103,]  1.2934
## [104,]  1.1067
## [105,]  1.0287
## [106,]  0.8986
## [107,]  0.5714
## [108,]  1.0786
## [109,]  0.3197
## [110,] -0.1719
## [111,]  1.0295
## [112,] -0.7104
## [113,]  0.6658
## [114,] -0.3107
## [115,] -0.5365
## [116,] -0.0738
## [117,] -1.0721
## [118,]  0.3533
## [119,] -0.3671
## [120,] -0.3069
## [121,]  0.5159
## [122,]  0.1661

4.6.4 項目特性曲線の作成

plot(test.C, type = 'trace', which.items = c(6, 10))

9章_因子分析