This exercise requires mirt for estimating multidimensional IRT models in Exercises 1-6 and eRm for Exercise 7.

For the first 5 exercises we will continue with the Workplace 1990 Workplace Industrial Relation Survey. The data are available as data(WIRS,package="ltm"). In the previous exercise we did a lot of parametric and nonparametric modelling of this difficult data set. We found that unidimensional parametric IRT does not fit the data. We saw that some assumptions are violated: There are intersecting IRF, there are non-monotonic IRF (Item 1, perhaps Item 2). Nonparametric IRT gave us a clear indication that there was something wrong. We also saw however, that the nonparametric model fit is not yet satisfactory. This follows when looking at the itemfit of the spline model, the credibility curves and the distribution on the latent trait.

  1. One assumption we have not yet checked is whether it makes sense to assume a single latent trait. Use plot(ksm,plottype="PCA") to investigate a two-dimensional representation of the item space. What do you find? Use aisp(WIRS,search="ga") from mokken to partition the data into sub-scales that fit the Monotone Homogeneity Model. What do you find?
library(mirt)
## Loading required package: stats4
## Loading required package: lattice
library(KernSmoothIRT)
library(mokken)
## Loading required package: poLCA
## Loading required package: scatterplot3d
## Loading required package: MASS
data(WIRS,package="ltm")
plot(ksm,plottype="PCA")

aisp(WIRS,search="ga")
##        0.3
## Item 1   0
## Item 2   2
## Item 3   1
## Item 4   2
## Item 5   1
## Item 6   1
  1. We saw that the items seem not to lie on a unidimensional scale but there is some indication that a two-dimensional parametric MHM model would fit Item 2-6. Use mirt() to fit two-dimensional 2PL, 3PL, 4PL models to Items 2-6 (if there are convergence issues either set technical=list(NCYCLES=2000) or method="MHRM" or combine both). Assess fit with statistics and plots. Which model is best and which would you choose? Look at the summary() to see the factor loadings.
subWIRS<-WIRS[,2:6]
td2pl<-mirt(subWIRS,2,itemtype="2PL", verbose=FALSE)
## EM cycles terminated after 500 iterations.
td3pl<-mirt(subWIRS,2,itemtype="3PL", verbose=FALSE)
## EM cycles terminated after 500 iterations.
td4pl<-mirt(subWIRS,2,itemtype="4PL", verbose=FALSE)
## EM cycles terminated after 500 iterations.
td2pl
## 
## Call:
## mirt(data = subWIRS, model = 2, itemtype = "2PL", verbose = FALSE)
## 
## Full-information item factor analysis with 2 factor(s).
## FAILED TO CONVERGE within 1e-04 tolerance after 500 EM iterations.
## mirt version: 1.24.5 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 31
## 
## Log-likelihood = -2748
## Estimated parameters: 14 
## AIC = 5523; AICc = 5524
## BIC = 5592; SABIC = 5548
## G2 (17) = 37, p = 0.0034
## RMSEA = 0.034, CFI = NaN, TLI = NaN
td3pl #Best
## 
## Call:
## mirt(data = subWIRS, model = 2, itemtype = "3PL", verbose = FALSE)
## 
## Full-information item factor analysis with 2 factor(s).
## FAILED TO CONVERGE within 1e-04 tolerance after 500 EM iterations.
## mirt version: 1.24.5 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 31
## 
## Log-likelihood = -2740
## Estimated parameters: 19 
## AIC = 5518; AICc = 5519
## BIC = 5611; SABIC = 5551
## G2 (12) = 21.38, p = 0.0451
## RMSEA = 0.028, CFI = NaN, TLI = NaN
td4pl
## 
## Call:
## mirt(data = subWIRS, model = 2, itemtype = "4PL", verbose = FALSE)
## 
## Full-information item factor analysis with 2 factor(s).
## FAILED TO CONVERGE within 1e-04 tolerance after 500 EM iterations.
## mirt version: 1.24.5 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 31
## 
## Log-likelihood = -2739
## Estimated parameters: 24 
## AIC = 5527; AICc = 5528
## BIC = 5645; SABIC = 5568
## G2 (7) = 20.26, p = 0.005
## RMSEA = 0.043, CFI = NaN, TLI = NaN
anova(td3pl,td4pl)
## 
## Model 1: mirt(data = subWIRS, model = 2, itemtype = "3PL", verbose = FALSE)
## Model 2: mirt(data = subWIRS, model = 2, itemtype = "4PL", verbose = FALSE)
##    AIC AICc SABIC  BIC logLik   X2  df     p
## 1 5518 5519  5551 5611  -2740  NaN NaN   NaN
## 2 5527 5528  5568 5645  -2739 1.12   5 0.952
anova(td2pl,td3pl)
## 
## Model 1: mirt(data = subWIRS, model = 2, itemtype = "2PL", verbose = FALSE)
## Model 2: mirt(data = subWIRS, model = 2, itemtype = "3PL", verbose = FALSE)
##    AIC AICc SABIC  BIC logLik    X2  df     p
## 1 5523 5524  5548 5592  -2748   NaN NaN   NaN
## 2 5518 5519  5551 5611  -2740 15.62   5 0.008
plot(td3pl) #Best model

itemplot(td3pl,1) #Scale 1

itemplot(td3pl,2) #Scale 2

itemplot(td3pl,3) #Scale 1

itemplot(td3pl,4) #Scale 2

itemplot(td3pl,5) #Scale 2

plot(td3pl,type="infocontour")

summary(td3pl)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##             F1      F2    h2
## Item 2  0.0227  0.9573 0.901
## Item 3 -0.7899 -0.1118 0.570
## Item 4 -0.6802  0.4509 0.898
## Item 5 -0.7694 -0.0379 0.571
## Item 6 -0.6626 -0.0993 0.399
## 
## Rotated SS loadings:  2.118 1.144 
## 
## Factor correlations: 
## 
##        F1     F2
## F1  1.000 -0.378
## F2 -0.378  1.000
  1. What to do with Item 1? We can use the results from the Mokken analysis and specify a confirmatory model with ideal point or nominal specification for Item 1 (Dimension 1) and a two dimensional 3PL for the other items (Items 2,4 being Dimension 1; Items 3,5,6 being Dimension 2). Compare the fit of both and inspect the IRFs.
cmodel <- mirt.model('F1 = 1
                      F2 = 2,4
                      F3 = 3,5,6
                      COV = F1*F2*F3')

cmod.3PLi <- mirt(WIRS, cmodel, itemtype=c("ideal",rep("3PL",5)), method = 'MHRM', verbose=FALSE)
## MHRM terminated after 2000 iterations.
cmod.3PLi
## 
## Call:
## mirt(data = WIRS, model = cmodel, itemtype = c("ideal", rep("3PL", 
##     5)), method = "MHRM", verbose = FALSE)
## 
## Full-information item factor analysis with 3 factor(s).
## FAILED TO CONVERGE within 0.001 tolerance after 2000 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = 0, SE = 0
## Estimated parameters: 20 
## AIC = 40; AICc = 40.85
## BIC = 138.3; SABIC = 74.73
## G2 (43) = 0, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
M2(cmod.3PLi)
##          M2 df         p  RMSEA RMSEA_5 RMSEA_95    TLI    CFI
## stats 13.95  1 0.0001878 0.1136 0.06593   0.1697 0.6473 0.9765
cmod.3PLn <- mirt(WIRS, cmodel, itemtype=c("nominal",rep("3PL",5)), method = 'MHRM',technical=list(NCYCLES=10000), verbose=FALSE)
cmod.3PLn
## 
## Call:
## mirt(data = WIRS, model = cmodel, itemtype = c("nominal", rep("3PL", 
##     5)), method = "MHRM", verbose = FALSE, technical = list(NCYCLES = 10000))
## 
## Full-information item factor analysis with 3 factor(s).
## Converged within 0.001 tolerance after 9027 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = -3339, SE = 0.017
## Estimated parameters: 20 
## AIC = 6719; AICc = 6720
## BIC = 6817; SABIC = 6754
## G2 (43) = 153.3, p = 0
## RMSEA = 0.051, CFI = NaN, TLI = NaN
M2(cmod.3PLn)
##          M2 df         p  RMSEA RMSEA_5 RMSEA_95    TLI    CFI
## stats 17.61  1 2.717e-05 0.1286 0.08047   0.1842 0.5478 0.9699
anova(cmod.3PLi,cmod.3PLn) #nominal is a bit better
## 
## Model 1: mirt(data = WIRS, model = cmodel, itemtype = c("nominal", rep("3PL", 
##     5)), method = "MHRM", verbose = FALSE, technical = list(NCYCLES = 10000))
## Model 2: mirt(data = WIRS, model = cmodel, itemtype = c("ideal", rep("3PL", 
##     5)), method = "MHRM", verbose = FALSE)
##    AIC    AICc   SABIC    BIC logLik   X2  df   p
## 1 6719 6719.67 6753.55 6817.1  -3339  NaN NaN NaN
## 2   40   40.85   74.73  138.3      0 6679   0   0
  1. Both models from before do not fit very well (but the nominal specification a bit better). One reason could be that there is a nonlinear effect on the latent traits. Indeed, when looking at the IRFs this could explain the non-monotonic effect revealed before. Fit a confirmatory 2D model with nominal specification for Item 1 and 2PL specification for all other items, where each item can load on each latent trait and include a multiplicative effect of the two latent dimensions (Hint: see the help file of mirt.model() how to do this). Fit a compensatory and a partially compensatory version. Which model fits best? Compare this to an exploratory 3D 3PL. Explore the overall best model.
##multiplicative effect of the latent traits (F1*F2)=1-6
##free covariance between the two factors
cmodel2 <- mirt.model('F1 = 1-6
                       F2 = 1-6
                       (F1*F2) =  1-6   
                       COV = F1*F2')    


#nominal specifiction for Item 1, compensatory 2D
cmod.2PLqn <- mirt(WIRS, cmodel2, itemtype=c("nominal",rep("2PL",5)), method = 'MHRM', verbose=FALSE)
cmod.2PLqn
## 
## Call:
## mirt(data = WIRS, model = cmodel2, itemtype = c("nominal", rep("2PL", 
##     5)), method = "MHRM", verbose = FALSE)
## 
## Full-information item factor analysis with 2 factor(s).
## Converged within 0.001 tolerance after 271 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = -3293, SE = 0.017
## Estimated parameters: 25 
## AIC = 6636; AICc = 6637
## BIC = 6758; SABIC = 6679
## G2 (38) = 60.08, p = 0.0127
## RMSEA = 0.024, CFI = NaN, TLI = NaN
summary(cmod.2PLqn)
##            F1      F2 (F1*F2)    h2
## Item 1  0.365 -0.0605 -0.9102 0.965
## Item 2 -0.112  0.6138  0.6728 0.842
## Item 3  0.711  0.1430 -0.0663 0.530
## Item 4  0.484 -0.0586  0.4954 0.483
## Item 5  0.673  0.3036 -0.1176 0.559
## Item 6  0.597  0.2089 -0.1404 0.420
## 
## SS loadings:  1.695 0.54 1.564 
## Proportion Var:  0.283 0.09 0.261 
## 
## Factor correlations: 
## 
##       F1    F2
## F1 1.000 0.386
## F2 0.386 1.000
#nominal specifictionfor Item 1, partially-compensatory 2D
cmod.2PLqnpc <- mirt(WIRS, cmodel2, itemtype=c("nominal",rep("PC2PL",5)), method = 'MHRM', verbose=FALSE)
## MHRM terminated after 2000 iterations.
cmod.2PLqnpc
## 
## Call:
## mirt(data = WIRS, model = cmodel2, itemtype = c("nominal", rep("PC2PL", 
##     5)), method = "MHRM", verbose = FALSE)
## 
## Full-information item factor analysis with 2 factor(s).
## FAILED TO CONVERGE within 0.001 tolerance after 2000 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = 0, SE = 0
## Estimated parameters: 35 
## AIC = 70; AICc = 72.6
## BIC = 241.9; SABIC = 130.8
## G2 (28) = 0, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
summary(cmod.2PLqnpc)
##             F1     F2 (F1*F2)    h2
## Item 1 -0.5816 0.1902  -0.734 0.914
## Item 2  0.2708 0.6196   0.709 0.960
## Item 3 -0.1792 0.0575  -0.975 0.986
## Item 4  0.0364 0.2434   0.889 0.851
## Item 5  0.2411 0.5901   0.542 0.700
## Item 6 -0.0994 0.5375  -0.630 0.695
## 
## SS loadings:  0.513 1.12 3.474 
## Proportion Var:  0.086 0.187 0.579 
## 
## Factor correlations: 
## 
##      F1   F2
## F1 1.00 0.22
## F2 0.22 1.00
anova(cmod.2PLqnpc,cmod.2PLqn) #compensatory model is best
## 
## Model 1: mirt(data = WIRS, model = cmodel2, itemtype = c("nominal", rep("2PL", 
##     5)), method = "MHRM", verbose = FALSE)
## Model 2: mirt(data = WIRS, model = cmodel2, itemtype = c("nominal", rep("PC2PL", 
##     5)), method = "MHRM", verbose = FALSE)
##    AIC   AICc  SABIC    BIC logLik   X2  df   p
## 1 6636 6636.9 6678.9 6758.3  -3293  NaN NaN NaN
## 2   70   72.6  130.8  241.9      0 6586  10   0
#3 dim 3PL 
thd3pl<-mirt(WIRS,3,itemtype="3PL",technical=list(NCYCLES=10000), verbose=FALSE)
anova(thd3pl,cmod.2PLqn) #compensatory model with nonlinear latent trait is best
## 
## Model 1: mirt(data = WIRS, model = cmodel2, itemtype = c("nominal", rep("2PL", 
##     5)), method = "MHRM", verbose = FALSE)
## Model 2: mirt(data = WIRS, model = 3, itemtype = "3PL", verbose = FALSE, 
##     technical = list(NCYCLES = 10000))
##    AIC AICc SABIC  BIC logLik    X2  df   p
## 1 6636 6637  6679 6758  -3293   NaN NaN NaN
## 2 6663 6664  6710 6796  -3304 -23.4   2   1
#Overall best model exploartion
summary(cmod.2PLqn)
##            F1      F2 (F1*F2)    h2
## Item 1  0.365 -0.0605 -0.9102 0.965
## Item 2 -0.112  0.6138  0.6728 0.842
## Item 3  0.711  0.1430 -0.0663 0.530
## Item 4  0.484 -0.0586  0.4954 0.483
## Item 5  0.673  0.3036 -0.1176 0.559
## Item 6  0.597  0.2089 -0.1404 0.420
## 
## SS loadings:  1.695 0.54 1.564 
## Proportion Var:  0.283 0.09 0.261 
## 
## Factor correlations: 
## 
##       F1    F2
## F1 1.000 0.386
## F2 0.386 1.000
coef(cmod.2PLqn)
## $`Item 1`
##        a1     a2     a3 ak0 ak1 d0     d1
## par 3.328 -0.552 -8.303   0   1  0 -0.088
## 
## $`Item 2`
##         a1    a2    a3      d g u
## par -0.481 2.629 2.882 -0.079 0 1
## 
## $`Item 3`
##        a1    a2     a3      d g u
## par 1.765 0.355 -0.165 -1.478 0 1
## 
## $`Item 4`
##        a1     a2    a3      d g u
## par 1.147 -0.139 1.173 -1.941 0 1
## 
## $`Item 5`
##        a1    a2     a3      d g u
## par 1.724 0.778 -0.301 -1.011 0 1
## 
## $`Item 6`
##        a1    a2     a3     d g u
## par 1.335 0.467 -0.314 -2.33 0 1
## 
## $GroupPars
##     MEAN_1 MEAN_2 COV_11 COV_21 COV_22
## par      0      0      1  0.386      1
plot(cmod.2PLqn)

plot(cmod.2PLqn,type="infocontour")

itemplot(cmod.2PLqn,1)

itemplot(cmod.2PLqn,2)

itemplot(cmod.2PLqn,3)

itemplot(cmod.2PLqn,4)

itemplot(cmod.2PLqn,5)

itemplot(cmod.2PLqn,6)

  1. We now turn to polytomous MIRT and mixed MIRT. For this, load the data file familiness.rda. Inspect it. You’ll find data on 11 Likert items (6 points scale) for three dimensions of family involvement that characterize family businesses (“Ownership, Management, Control”= OMC, “Transgenerational Orientation”=TGO, “Family Business Involvement”= FBI) and two explanatory variables (the age of the firm and the size of the firm). Fit an exploratory 1-dimensional GPCM and a 3-dimensional GPCM. Is the 3D solution better? Fit a confirmatory model with no item cross-loading to other dimensions (so OMC only loads on OMC) etc. Would you say that the hypothesized model is worse than the exploratory model?
load("familiness.rda")

modpcm1<-mirt(familiness[,1:11],model=1,itemtype="gpcm", verbose=FALSE)
modpcm1
## 
## Call:
## mirt(data = familiness[, 1:11], model = 1, itemtype = "gpcm", 
##     verbose = FALSE)
## 
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 66 EM iterations.
## mirt version: 1.24.5 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 61
## 
## Log-likelihood = -6987
## Estimated parameters: 66 
## AIC = 14106; AICc = 14126
## BIC = 14385; SABIC = 14175
## G2 (362796989) = 7822, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
itemfit(modpcm1)
##    item   S_X2 df.S_X2 p.S_X2
## 1  OMC1  77.92      58  0.042
## 2  OMC2  73.05      67  0.286
## 3  OMC3 121.96      95  0.033
## 4  OMC4  69.66      59  0.161
## 5  TGO1  83.20      68  0.101
## 6  TGO2  93.84      72  0.043
## 7  TGO3  91.60      66  0.020
## 8  FBI1  43.68      53  0.816
## 9  FBI2  66.82      74  0.711
## 10 FBI5  71.18      77  0.665
## 11 FBI6  71.35      72  0.500
modpcm3<-mirt(familiness[,1:11],model=3,itemtype="gpcm",method="MHRM", verbose=FALSE)
modpcm3
## 
## Call:
## mirt(data = familiness[, 1:11], model = 3, itemtype = "gpcm", 
##     method = "MHRM", verbose = FALSE)
## 
## Full-information item factor analysis with 3 factor(s).
## Converged within 0.001 tolerance after 841 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = -6636, SE = 0.064
## Estimated parameters: 85 
## AIC = 13443; AICc = 13478
## BIC = 13802; SABIC = 13532
## G2 (362796970) = 7121, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
itemfit(modpcm3)
##    item   S_X2 df.S_X2 p.S_X2
## 1  OMC1  73.89      54  0.037
## 2  OMC2  78.40      63  0.091
## 3  OMC3 125.78      97  0.026
## 4  OMC4  55.18      52  0.355
## 5  TGO1  87.91      73  0.112
## 6  TGO2  79.82      69  0.175
## 7  TGO3  77.84      68  0.194
## 8  FBI1  39.39      47  0.777
## 9  FBI2  76.65      77  0.490
## 10 FBI5  71.35      74  0.566
## 11 FBI6  67.96      68  0.478
summary(modpcm3)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##             F1       F2       F3     h2
## OMC1 -0.001947  0.78598 -0.05067 0.6570
## OMC2 -0.000115  0.81498  0.05608 0.6271
## OMC3 -0.003437  0.14918 -0.04256 0.0303
## OMC4 -0.031786  0.56873 -0.06419 0.3800
## TGO1 -0.141011  0.00472 -0.54641 0.4005
## TGO2  0.011564  0.06454 -0.73359 0.5748
## TGO3  0.016506 -0.02604 -0.89550 0.7675
## FBI1 -0.625132  0.06902 -0.03662 0.4636
## FBI2 -0.380411 -0.13128 -0.08676 0.1459
## FBI5 -0.797607  0.00292  0.01583 0.6256
## FBI6 -0.606838 -0.02203  0.00872 0.3508
## 
## Rotated SS loadings:  1.561 1.655 1.659 
## 
## Factor correlations: 
## 
##        F1     F2     F3
## F1  1.000 -0.475  0.514
## F2 -0.475  1.000 -0.441
## F3  0.514 -0.441  1.000
cmod <- mirt.model('OMC = 1-4
                    TGO= 5,6,7
                    FBI = 8-11
                    COV = OMC*TGO*FBI')

cmodpcm3 <- mirt(familiness[,1:11], cmod, itemtype="gpcm", method = 'MHRM', verbose=FALSE)
cmodpcm3 #cmod3pcm and modpcm3 are similar, modpcm3 slightly better cmodpcm3 much easier though
## 
## Call:
## mirt(data = familiness[, 1:11], model = cmod, itemtype = "gpcm", 
##     method = "MHRM", verbose = FALSE)
## 
## Full-information item factor analysis with 3 factor(s).
## Converged within 0.001 tolerance after 952 MHRM iterations.
## mirt version: 1.24.5 
## M-step optimizer: NR1 
## 
## Log-likelihood = -6664, SE = 0.063
## Estimated parameters: 69 
## AIC = 13467; AICc = 13489
## BIC = 13759; SABIC = 13539
## G2 (362796986) = 7177, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
summary(cmodpcm3)
##        OMC   TGO   FBI     h2
## OMC1 0.816 0.000 0.000 0.6657
## OMC2 0.732 0.000 0.000 0.5357
## OMC3 0.179 0.000 0.000 0.0319
## OMC4 0.627 0.000 0.000 0.3927
## TGO1 0.000 0.627 0.000 0.3934
## TGO2 0.000 0.749 0.000 0.5604
## TGO3 0.000 0.863 0.000 0.7455
## FBI1 0.000 0.000 0.688 0.4737
## FBI2 0.000 0.000 0.366 0.1337
## FBI5 0.000 0.000 0.791 0.6262
## FBI6 0.000 0.000 0.577 0.3333
## 
## SS loadings:  1.626 1.699 1.567 
## Proportion Var:  0.148 0.154 0.142 
## 
## Factor correlations: 
## 
##       OMC   TGO   FBI
## OMC 1.000 0.506 0.430
## TGO 0.506 1.000 0.537
## FBI 0.430 0.537 1.000
summary(modpcm3)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##             F1       F2       F3     h2
## OMC1 -0.001947  0.78598 -0.05067 0.6570
## OMC2 -0.000115  0.81498  0.05608 0.6271
## OMC3 -0.003437  0.14918 -0.04256 0.0303
## OMC4 -0.031786  0.56873 -0.06419 0.3800
## TGO1 -0.141011  0.00472 -0.54641 0.4005
## TGO2  0.011564  0.06454 -0.73359 0.5748
## TGO3  0.016506 -0.02604 -0.89550 0.7675
## FBI1 -0.625132  0.06902 -0.03662 0.4636
## FBI2 -0.380411 -0.13128 -0.08676 0.1459
## FBI5 -0.797607  0.00292  0.01583 0.6256
## FBI6 -0.606838 -0.02203  0.00872 0.3508
## 
## Rotated SS loadings:  1.561 1.655 1.659 
## 
## Factor correlations: 
## 
##        F1     F2     F3
## F1  1.000 -0.475  0.514
## F2 -0.475  1.000 -0.441
## F3  0.514 -0.441  1.000
itemfit(cmodpcm3)
##    item   S_X2 df.S_X2 p.S_X2
## 1  OMC1  69.91      48  0.021
## 2  OMC2  77.21      59  0.056
## 3  OMC3 126.00      99  0.035
## 4  OMC4  59.54      51  0.193
## 5  TGO1  95.70      76  0.063
## 6  TGO2  80.56      69  0.161
## 7  TGO3  80.24      63  0.070
## 8  FBI1  39.80      50  0.849
## 9  FBI2  69.89      73  0.582
## 10 FBI5  71.04      73  0.543
## 11 FBI6  69.55      68  0.425
anova(modpcm3,cmodpcm3)
## 
## Model 1: mirt(data = familiness[, 1:11], model = cmod, itemtype = "gpcm", 
##     method = "MHRM", verbose = FALSE)
## Model 2: mirt(data = familiness[, 1:11], model = 3, itemtype = "gpcm", 
##     method = "MHRM", verbose = FALSE)
##     AIC  AICc SABIC   BIC logLik    X2  df   p
## 1 13467 13489 13539 13759  -6664   NaN NaN NaN
## 2 13443 13478 13532 13802  -6636 55.87  16   0
  1. For the familiness data, include a fixed effect for firm Age and Size in a latent regression model on the confirmatory model structure (i.e., decomposing the latent trait, lr.fixed). Inspect the coefficients and the model fit and interpret it. Compare the explanatory model to the confirmatory model from 5).
covs<-familiness[,c(12,13)] #covariates
covs[,1]<-scale(covs[,1])  #scale covariates for numerical stability

ecmod3 <- mixedmirt(familiness[,1:11], covdata=covs, model=cmod, lr.fixed = ~ 1 + Age + Size,  itemtype = 'gpcm', verbose=FALSE)
summary(ecmod3)
## 
## Call:
## mixedmirt(data = familiness[, 1:11], covdata = covs, model = cmod, 
##     itemtype = "gpcm", lr.fixed = ~1 + Age + Size, verbose = FALSE)
## 
## 
## --------------
## RANDOM EFFECT COVARIANCE(S):
## Correlations on upper diagonal
## 
## $Theta
##       OMC   TGO   FBI
## OMC 1.000 0.517 0.433
## TGO 0.517 1.000 0.526
## FBI 0.433 0.526 1.000
## 
## --------------
## LATENT REGRESSION FIXED EFFECTS:
## 
##               OMC   TGO   FBI
## (Intercept) 0.000 0.000 0.000
## Age         0.028 0.082 0.037
## Size        0.018 0.197 0.065
## 
##             Std.Error_OMC Std.Error_TGO Std.Error_FBI z_OMC z_TGO z_FBI
## (Intercept)            NA            NA            NA    NA    NA    NA
## Age                 0.041         0.041         0.036 0.692 2.017 1.003
## Size                  NaN           NaN         0.022   NaN   NaN 2.994
coef(ecmod3)
## $OMC1
##            a1 a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2     d3     d4
## par     2.253  0  0   0   1   2   3   4   5  0 3.045 5.469  9.171 11.693
## CI_2.5  1.310 NA NA  NA  NA  NA  NA  NA  NA NA 1.231 2.315  4.860  6.341
## CI_97.5 3.197 NA NA  NA  NA  NA  NA  NA  NA NA 4.858 8.622 13.482 17.044
##             d5
## par     13.781
## CI_2.5   8.017
## CI_97.5 19.545
## 
## $OMC2
##            a1 a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2     d3     d4
## par     1.831  0  0   0   1   2   3   4   5  0 2.345 5.892  8.056 10.158
## CI_2.5  1.262 NA NA  NA  NA  NA  NA  NA  NA NA 1.082 4.165  5.815  7.524
## CI_97.5 2.400 NA NA  NA  NA  NA  NA  NA  NA NA 3.608 7.618 10.296 12.792
##             d5
## par     11.685
## CI_2.5   9.042
## CI_97.5 14.327
## 
## $OMC3
##            a1 a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0     d1     d2     d3    d4
## par     0.306  0  0   0   1   2   3   4   5  0 -0.436  0.167  0.217 0.783
## CI_2.5  0.219 NA NA  NA  NA  NA  NA  NA  NA NA -0.836 -0.202 -0.160 0.450
## CI_97.5 0.394 NA NA  NA  NA  NA  NA  NA  NA NA -0.036  0.535  0.595 1.116
##            d5
## par     0.886
## CI_2.5  0.588
## CI_97.5 1.184
## 
## $OMC4
##            a1 a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3     d4
## par     1.394  0  0   0   1   2   3   4   5  0 1.877 4.023 6.721  8.208
## CI_2.5  0.915 NA NA  NA  NA  NA  NA  NA  NA NA 0.725 2.453 4.793  5.982
## CI_97.5 1.872 NA NA  NA  NA  NA  NA  NA  NA NA 3.028 5.594 8.648 10.434
##             d5
## par      9.707
## CI_2.5   7.531
## CI_97.5 11.883
## 
## $TGO1
##         a1    a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3    d4
## par      0 1.363  0   0   1   2   3   4   5  0 1.269 3.006 3.692 3.508
## CI_2.5  NA 1.075 NA  NA  NA  NA  NA  NA  NA NA 0.634 2.293 2.901 2.693
## CI_97.5 NA 1.650 NA  NA  NA  NA  NA  NA  NA NA 1.904 3.719 4.483 4.322
##            d5
## par     2.293
## CI_2.5  1.456
## CI_97.5 3.129
## 
## $TGO2
##         a1    a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3    d4
## par      0 1.866  0   0   1   2   3   4   5  0 1.658 3.852 4.611 4.513
## CI_2.5  NA 1.411 NA  NA  NA  NA  NA  NA  NA NA 0.850 2.764 3.285 3.058
## CI_97.5 NA 2.320 NA  NA  NA  NA  NA  NA  NA NA 2.466 4.940 5.937 5.967
##            d5
## par     4.074
## CI_2.5  2.672
## CI_97.5 5.477
## 
## $TGO3
##         a1    a2 a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3    d4
## par      0 2.846  0   0   1   2   3   4   5  0 2.301 5.346 6.602 6.457
## CI_2.5  NA 1.876 NA  NA  NA  NA  NA  NA  NA NA 0.924 3.138 3.823 3.477
## CI_97.5 NA 3.815 NA  NA  NA  NA  NA  NA  NA NA 3.679 7.553 9.380 9.438
##            d5
## par     5.021
## CI_2.5  2.351
## CI_97.5 7.691
## 
## $FBI1
##         a1 a2    a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2     d3     d4
## par      0  0 1.587   0   1   2   3   4   5  0 3.297 7.341  9.902 11.873
## CI_2.5  NA NA 1.210  NA  NA  NA  NA  NA  NA NA 1.508 5.215  7.365  9.061
## CI_97.5 NA NA 1.965  NA  NA  NA  NA  NA  NA NA 5.087 9.466 12.439 14.685
##            d5
## par     12.46
## CI_2.5   9.62
## CI_97.5 15.30
## 
## $FBI2
##         a1 a2    a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3    d4
## par      0  0 0.660   0   1   2   3   4   5  0 1.156 2.904 4.026 4.219
## CI_2.5  NA NA 0.504  NA  NA  NA  NA  NA  NA NA 0.334 2.096 3.170 3.341
## CI_97.5 NA NA 0.816  NA  NA  NA  NA  NA  NA NA 1.979 3.712 4.883 5.098
##            d5
## par     3.658
## CI_2.5  2.803
## CI_97.5 4.513
## 
## $FBI5
##         a1 a2    a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2     d3     d4
## par      0  0 2.274   0   1   2   3   4   5  0 3.417 7.234  8.897 10.062
## CI_2.5  NA NA 1.606  NA  NA  NA  NA  NA  NA NA 2.058 5.160  6.259  7.121
## CI_97.5 NA NA 2.941  NA  NA  NA  NA  NA  NA NA 4.776 9.308 11.534 13.003
##             d5
## par      9.370
## CI_2.5   6.526
## CI_97.5 12.214
## 
## $FBI6
##         a1 a2    a3 ak0 ak1 ak2 ak3 ak4 ak5 d0    d1    d2    d3    d4
## par      0  0 1.207   0   1   2   3   4   5  0 1.425 4.425 5.754 5.985
## CI_2.5  NA NA 0.921  NA  NA  NA  NA  NA  NA NA 0.445 3.339 4.468 4.595
## CI_97.5 NA NA 1.493  NA  NA  NA  NA  NA  NA NA 2.404 5.510 7.039 7.375
##            d5
## par     5.733
## CI_2.5  4.391
## CI_97.5 7.075
## 
## $GroupPars
##         MEAN_1 MEAN_2 MEAN_3 COV_11 COV_21 COV_31 COV_22 COV_32 COV_33
## par          0      0      0      1  0.517  0.433      1  0.526      1
## CI_2.5      NA     NA     NA     NA  0.413  0.331     NA  0.439     NA
## CI_97.5     NA     NA     NA     NA  0.621  0.536     NA  0.612     NA
## 
## $lr.betas
##         OMC_(Intercept) OMC_Age OMC_Size TGO_(Intercept) TGO_Age TGO_Size
## par                   0   0.028    0.018               0   0.082    0.197
## CI_2.5               NA  -0.052      NaN              NA   0.002      NaN
## CI_97.5              NA   0.109      NaN              NA   0.161      NaN
##         FBI_(Intercept) FBI_Age FBI_Size
## par                   0   0.037    0.065
## CI_2.5               NA  -0.035    0.023
## CI_97.5              NA   0.108    0.108
anova(cmodpcm3, ecmod3)
## 
## Model 1: mirt(data = familiness[, 1:11], model = cmod, itemtype = "gpcm", 
##     method = "MHRM", verbose = FALSE)
## Model 2: mixedmirt(data = familiness[, 1:11], covdata = covs, model = cmod, 
##     itemtype = "gpcm", lr.fixed = ~1 + Age + Size, verbose = FALSE)
##     AIC  AICc SABIC   BIC logLik    X2  df     p
## 1 13467 13489 13539 13759  -6664   NaN NaN   NaN
## 2 13462 13489 13541 13779  -6656 16.47   6 0.011
  1. Load the AMC data object. AMC 2 is an exam on accounting at WU. We look at the 12 dichotomously scored items. The students were assessed on two different points in time with a break in between. The first round was asking multiple choice format, the second asked the same questions but in an open format. It was of interest to see whether multiple choice is different than open format. Fit an LLRA model for the AMCit with two measurement timepoints. This quantifies the change over time (trend) for each item as its own dimension. Assess the model and plot it. What do the results mean in light of the design setup? What is problematic in the setup? From the AMCcov data frame, select HS the type of high school (AHS is standard high school, HAK is a high school with focus on commerce, other is the rest). Include it as a groups covariate in the analysis. Assess the model in comparison to the model without a group and plot it. What do we learn from this analysis?
library(eRm)
## 
## Attaching package: 'eRm'
## The following objects are masked from 'package:mirt':
## 
##     itemfit, personfit
load("AMC.rda")

# no groups
modllra0 <- LLRA(AMCit,mpoints=2)
## Reference group:  CG
summary(modllra0)
## 
## Results of LLRA via LPCM estimation: 
## 
## Call:  LLRA(X = AMCit, mpoints = 2) 
## 
## Conditional log-likelihood: -859.6 
## Number of iterations: 14 
## Number of parameters: 12 
## 
## Estimated parameters with 0.95 CI:
##              Estimate Std.Error lower.CI upper.CI
## trend.I1.t2    -2.244     0.219   -2.674   -1.815
## trend.I2.t2    -1.905     0.253   -2.401   -1.410
## trend.I3.t2    -1.792     0.230   -2.243   -1.340
## trend.I4.t2     1.478     0.217    1.052    1.904
## trend.I5.t2    -1.616     0.197   -2.001   -1.230
## trend.I6.t2    -1.826     0.220   -2.257   -1.395
## trend.I7.t2    -2.205     0.272   -2.738   -1.671
## trend.I8.t2    -3.531     0.383   -4.282   -2.779
## trend.I9.t2    -0.883     0.181   -1.239   -0.528
## trend.I10.t2   -3.124     0.308   -3.727   -2.520
## trend.I11.t2    0.912     0.176    0.566    1.258
## trend.I12.t2    0.418     0.143    0.137    0.698
## 
## Reference Group:  CG
plotTR(modllra0) #Items 4, 11 ,12 are relatively easier in open question the others more difficult; confounding of response effect with trend though. We have no comparison group at t2/t1 of the other format. 

modllra1 <- LLRA(AMCit,mpoints=2,groups=AMCcov$HS)
## Reference group:  AHS
summary(modllra1,gamma=0.95)
## 
## Results of LLRA via LPCM estimation: 
## 
## Call:  LLRA(X = AMCit, mpoints = 2, groups = AMCcov$HS) 
## 
## Conditional log-likelihood: -839.8 
## Number of iterations: 39 
## Number of parameters: 36 
## 
## Estimated parameters with 0.95 CI:
##              Estimate Std.Error lower.CI upper.CI
## other.I1.t2    -0.194     0.798   -1.758    1.369
## HAK.I1.t2       1.276     0.470    0.354    2.198
## other.I2.t2    -0.519     0.806   -2.100    1.061
## HAK.I2.t2       0.192     0.585   -0.953    1.338
## other.I3.t2     1.159     0.587    0.008    2.310
## HAK.I3.t2       0.871     0.549   -0.205    1.947
## other.I4.t2     0.086     0.625   -1.138    1.311
## HAK.I4.t2      -0.025     0.492   -0.990    0.940
## other.I5.t2     1.159     0.496    0.187    2.131
## HAK.I5.t2       1.041     0.476    0.109    1.974
## other.I6.t2     1.028     0.536   -0.023    2.078
## HAK.I6.t2       0.790     0.554   -0.297    1.877
## other.I7.t2     0.454     0.721   -0.959    1.868
## HAK.I7.t2       0.272     0.647   -0.995    1.539
## other.I8.t2    -0.110     1.131   -2.328    2.107
## HAK.I8.t2       0.606     0.883   -1.125    2.337
## other.I9.t2     1.127     0.460    0.226    2.028
## HAK.I9.t2       0.078     0.500   -0.903    1.058
## other.I10.t2    0.328     0.823   -1.285    1.941
## HAK.I10.t2      0.077     0.819   -1.529    1.682
## other.I11.t2    0.305     0.522   -0.717    1.328
## HAK.I11.t2      0.069     0.412   -0.738    0.876
## other.I12.t2   -0.636     0.365   -1.352    0.080
## HAK.I12.t2     -0.826     0.365   -1.542   -0.110
## trend.I1.t2    -2.639     0.327   -3.281   -1.998
## trend.I2.t2    -1.879     0.324   -2.513   -1.244
## trend.I3.t2    -2.257     0.350   -2.944   -1.571
## trend.I4.t2     1.472     0.296    0.891    2.053
## trend.I5.t2    -2.140     0.305   -2.738   -1.542
## trend.I6.t2    -2.217     0.317   -2.839   -1.595
## trend.I7.t2    -2.351     0.370   -3.077   -1.626
## trend.I8.t2    -3.651     0.506   -4.643   -2.658
## trend.I9.t2    -1.127     0.240   -1.597   -0.657
## trend.I10.t2   -3.190     0.386   -3.946   -2.434
## trend.I11.t2    0.847     0.230    0.396    1.298
## trend.I12.t2    0.731     0.195    0.349    1.113
## 
## Reference Group:  AHS
plotTR(modllra1) #picture similar to before

plotGR(modllra1) #Items 1, 3, 5, 6 are easier for students who were in HAK in open question formats (For item 1, 5 its significant at 0.05); Item 12 is more difficult for HAK students in open format (significant)

Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.