De Boeck, P., & Partchev, I. (2012). IRTrees: Tree-based item response models of the GLMM family. Journal of Statistical Software, 48, 1-28.
library(irtrees)Estimation: Multidimensional Item Response Model
Data Transformation
library(mirt)## Loading required package: stats4## Loading required package: latticemapping <- cbind(c(0, 1, 1), c(NA, 0, 1))
head(VerbAgg3[ , -(1:2)])## S1WantCurse S1WantScold S1WantShout S2WantCurse S2WantScold S2WantShout
## [1,] 1 1 1 1 1 1
## [2,] 1 1 1 1 1 1
## [3,] 2 2 2 2 1 2
## [4,] 2 2 2 2 2 2
## [5,] 2 1 2 2 1 1
## [6,] 3 3 1 3 1 1
## S3WantCurse S3WantScold S3WantShout S4wantCurse S4WantScold S4WantShout
## [1,] 1 1 2 3 1 1
## [2,] 1 1 1 1 1 1
## [3,] 2 1 1 1 1 1
## [4,] 2 1 1 1 1 1
## [5,] 2 1 1 2 1 1
## [6,] 2 2 2 1 1 1
## S1DoCurse S1DoScold S1DoShout S2DoCurse S2DoScold S2DoShout S3DoCurse
## [1,] 2 1 2 2 1 1 2
## [2,] 1 1 1 1 1 1 2
## [3,] 1 2 2 1 1 2 1
## [4,] 2 2 2 3 2 2 2
## [5,] 2 2 1 2 1 1 2
## [6,] 3 1 1 3 1 1 3
## S3DoScold S3DoShout S4DoCurse S4DoScold S4DoShout
## [1,] 1 1 3 3 3
## [2,] 1 1 1 1 1
## [3,] 1 1 2 1 1
## [4,] 1 1 1 1 1
## [5,] 1 1 2 1 1
## [6,] 1 1 3 3 1VerbAgg3T_wide <- WtoW_single.tree(VerbAgg3[ , -(1:2)], mapping)
head(VerbAgg3T_wide)## S1WantCurse:node1 S1WantScold:node1 S1WantShout:node1 S2WantCurse:node1
## 1 0 0 0 0
## 2 0 0 0 0
## 3 1 1 1 1
## 4 1 1 1 1
## 5 1 0 1 1
## 6 1 1 0 1
## S2WantScold:node1 S2WantShout:node1 S3WantCurse:node1 S3WantScold:node1
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 1 1 0
## 4 1 1 1 0
## 5 0 0 1 0
## 6 0 0 1 1
## S3WantShout:node1 S4wantCurse:node1 S4WantScold:node1 S4WantShout:node1
## 1 1 1 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 1 0 0
## 6 1 0 0 0
## S1DoCurse:node1 S1DoScold:node1 S1DoShout:node1 S2DoCurse:node1
## 1 1 0 1 1
## 2 0 0 0 0
## 3 0 1 1 0
## 4 1 1 1 1
## 5 1 1 0 1
## 6 1 0 0 1
## S2DoScold:node1 S2DoShout:node1 S3DoCurse:node1 S3DoScold:node1
## 1 0 0 1 0
## 2 0 0 1 0
## 3 0 1 0 0
## 4 1 1 1 0
## 5 0 0 1 0
## 6 0 0 1 0
## S3DoShout:node1 S4DoCurse:node1 S4DoScold:node1 S4DoShout:node1
## 1 0 1 1 1
## 2 0 0 0 0
## 3 0 1 0 0
## 4 0 0 0 0
## 5 0 1 0 0
## 6 0 1 1 0
## S1WantCurse:node2 S1WantScold:node2 S1WantShout:node2 S2WantCurse:node2
## 1 NA NA NA NA
## 2 NA NA NA NA
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 NA 0 0
## 6 1 1 NA 1
## S2WantScold:node2 S2WantShout:node2 S3WantCurse:node2 S3WantScold:node2
## 1 NA NA NA NA
## 2 NA NA NA NA
## 3 NA 0 0 NA
## 4 0 0 0 NA
## 5 NA NA 0 NA
## 6 NA NA 0 0
## S3WantShout:node2 S4wantCurse:node2 S4WantScold:node2 S4WantShout:node2
## 1 0 1 NA NA
## 2 NA NA NA NA
## 3 NA NA NA NA
## 4 NA NA NA NA
## 5 NA 0 NA NA
## 6 0 NA NA NA
## S1DoCurse:node2 S1DoScold:node2 S1DoShout:node2 S2DoCurse:node2
## 1 0 NA 0 0
## 2 NA NA NA NA
## 3 NA 0 0 NA
## 4 0 0 0 1
## 5 0 0 NA 0
## 6 1 NA NA 1
## S2DoScold:node2 S2DoShout:node2 S3DoCurse:node2 S3DoScold:node2
## 1 NA NA 0 NA
## 2 NA NA 0 NA
## 3 NA 0 NA NA
## 4 0 0 0 NA
## 5 NA NA 0 NA
## 6 NA NA 1 NA
## S3DoShout:node2 S4DoCurse:node2 S4DoScold:node2 S4DoShout:node2
## 1 NA 1 1 1
## 2 NA NA NA NA
## 3 NA 0 NA NA
## 4 NA NA NA NA
## 5 NA 0 NA NA
## 6 NA 1 1 NAUnidimensional model
IRTree.mod1 = 'degree = 1-48'
tree.mirt1 = mirt(data = VerbAgg3T_wide,
model = IRTree.mod1,
itemtype = 'Rasch')##
Iteration: 1, Log-Lik: -7478.058, Max-Change: 5.05521
Iteration: 2, Log-Lik: -6393.893, Max-Change: 0.37432
Iteration: 3, Log-Lik: -6369.233, Max-Change: 0.15296
Iteration: 4, Log-Lik: -6360.713, Max-Change: 0.08112
Iteration: 5, Log-Lik: -6355.947, Max-Change: 0.05219
Iteration: 6, Log-Lik: -6352.507, Max-Change: 0.03911
Iteration: 7, Log-Lik: -6340.422, Max-Change: 0.01685
Iteration: 8, Log-Lik: -6340.070, Max-Change: 0.01418
Iteration: 9, Log-Lik: -6339.765, Max-Change: 0.01175
Iteration: 10, Log-Lik: -6338.677, Max-Change: 0.01384
Iteration: 11, Log-Lik: -6338.613, Max-Change: 0.00720
Iteration: 12, Log-Lik: -6338.576, Max-Change: 0.00500
Iteration: 13, Log-Lik: -6338.499, Max-Change: 0.00651
Iteration: 14, Log-Lik: -6338.485, Max-Change: 0.00213
Iteration: 15, Log-Lik: -6338.478, Max-Change: 0.00187
Iteration: 16, Log-Lik: -6338.459, Max-Change: 0.00368
Iteration: 17, Log-Lik: -6338.456, Max-Change: 0.00102
Iteration: 18, Log-Lik: -6338.455, Max-Change: 0.00064
Iteration: 19, Log-Lik: -6338.453, Max-Change: 0.00106
Iteration: 20, Log-Lik: -6338.453, Max-Change: 0.00034
Iteration: 21, Log-Lik: -6338.452, Max-Change: 0.00023
Iteration: 22, Log-Lik: -6338.452, Max-Change: 0.00048
Iteration: 23, Log-Lik: -6338.452, Max-Change: 0.00012
Iteration: 24, Log-Lik: -6338.452, Max-Change: 0.00007tree.mirt1##
## Call:
## mirt(data = VerbAgg3T_wide, model = IRTree.mod1, itemtype = "Rasch")
##
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 24 EM iterations.
## mirt version: 1.37.1
## M-step optimizer: nlminb
## EM acceleration: Ramsay
## Number of rectangular quadrature: 61
## Latent density type: Gaussian
##
## Log-likelihood = -6338.452
## Estimated parameters: 49
## AIC = 12774.9
## BIC = 12958.94; SABIC = 12803.52coef(tree.mirt1, simplify = TRUE)## $items
## a1 d g u
## S1WantCurse:node1 1 1.125 0 1
## S1WantScold:node1 1 0.510 0 1
## S1WantShout:node1 1 0.059 0 1
## S2WantCurse:node1 1 1.625 0 1
## S2WantScold:node1 1 0.643 0 1
## S2WantShout:node1 1 -0.004 0 1
## S3WantCurse:node1 1 0.477 0 1
## S3WantScold:node1 1 -0.649 0 1
## S3WantShout:node1 1 -1.423 0 1
## S4wantCurse:node1 1 0.994 0 1
## S4WantScold:node1 1 -0.339 0 1
## S4WantShout:node1 1 -0.979 0 1
## S1DoCurse:node1 1 1.125 0 1
## S1DoScold:node1 1 0.347 0 1
## S1DoShout:node1 1 -0.820 0 1
## S2DoCurse:node1 1 0.798 0 1
## S2DoScold:node1 1 -0.068 0 1
## S2DoShout:node1 1 -1.382 0 1
## S3DoCurse:node1 1 -0.211 0 1
## S3DoScold:node1 1 -1.403 0 1
## S3DoShout:node1 1 -2.765 0 1
## S4DoCurse:node1 1 0.643 0 1
## S4DoScold:node1 1 -0.371 0 1
## S4DoShout:node1 1 -1.859 0 1
## S1WantCurse:node2 1 0.090 0 1
## S1WantScold:node2 1 -0.181 0 1
## S1WantShout:node2 1 -1.037 0 1
## S2WantCurse:node2 1 0.023 0 1
## S2WantScold:node2 1 -0.252 0 1
## S2WantShout:node2 1 -0.646 0 1
## S3WantCurse:node2 1 -1.041 0 1
## S3WantScold:node2 1 -2.036 0 1
## S3WantShout:node2 1 -2.600 0 1
## S4wantCurse:node2 1 -0.702 0 1
## S4WantScold:node2 1 -1.307 0 1
## S4WantShout:node2 1 -1.356 0 1
## S1DoCurse:node2 1 -0.254 0 1
## S1DoScold:node2 1 -0.713 0 1
## S1DoShout:node2 1 -1.349 0 1
## S2DoCurse:node2 1 -0.268 0 1
## S2DoScold:node2 1 -1.078 0 1
## S2DoShout:node2 1 -1.778 0 1
## S3DoCurse:node2 1 -1.814 0 1
## S3DoScold:node2 1 -2.442 0 1
## S3DoShout:node2 1 -3.306 0 1
## S4DoCurse:node2 1 -0.855 0 1
## S4DoScold:node2 1 -1.465 0 1
## S4DoShout:node2 1 -2.196 0 1
##
## $means
## degree
## 0
##
## $cov
## degree
## degree 1.329Multidimension model
IRTree.mod2 = 'admission = 1-24
affirmation = 25-48
COV = admission*affirmation'
tree.mirt2 = mirt(data = VerbAgg3T_wide,
model = IRTree.mod2,
itemtype = 'Rasch')##
Iteration: 1, Log-Lik: -6674.100, Max-Change: 3.94435
Iteration: 2, Log-Lik: -6230.576, Max-Change: 0.63205
Iteration: 3, Log-Lik: -6192.488, Max-Change: 0.25989
Iteration: 4, Log-Lik: -6175.944, Max-Change: 0.16037
Iteration: 5, Log-Lik: -6164.881, Max-Change: 0.12241
Iteration: 6, Log-Lik: -6156.280, Max-Change: 0.12810
Iteration: 7, Log-Lik: -6149.474, Max-Change: 0.11313
Iteration: 8, Log-Lik: -6144.235, Max-Change: 0.09087
Iteration: 9, Log-Lik: -6140.366, Max-Change: 0.06848
Iteration: 10, Log-Lik: -6133.246, Max-Change: 0.10034
Iteration: 11, Log-Lik: -6132.428, Max-Change: 0.04052
Iteration: 12, Log-Lik: -6132.243, Max-Change: 0.01662
Iteration: 13, Log-Lik: -6132.163, Max-Change: 0.01133
Iteration: 14, Log-Lik: -6132.132, Max-Change: 0.00771
Iteration: 15, Log-Lik: -6132.112, Max-Change: 0.00503
Iteration: 16, Log-Lik: -6132.043, Max-Change: 0.00621
Iteration: 17, Log-Lik: -6132.049, Max-Change: 0.00273
Iteration: 18, Log-Lik: -6132.051, Max-Change: 0.00156
Iteration: 19, Log-Lik: -6132.046, Max-Change: 0.00263
Iteration: 20, Log-Lik: -6132.050, Max-Change: 0.00096
Iteration: 21, Log-Lik: -6132.051, Max-Change: 0.00071
Iteration: 22, Log-Lik: -6132.050, Max-Change: 0.00153
Iteration: 23, Log-Lik: -6132.053, Max-Change: 0.00050
Iteration: 24, Log-Lik: -6132.054, Max-Change: 0.00026
Iteration: 25, Log-Lik: -6132.054, Max-Change: 0.00044
Iteration: 26, Log-Lik: -6132.055, Max-Change: 0.00016
Iteration: 27, Log-Lik: -6132.055, Max-Change: 0.00013
Iteration: 28, Log-Lik: -6132.055, Max-Change: 0.00025
Iteration: 29, Log-Lik: -6132.055, Max-Change: 0.00009tree.mirt2##
## Call:
## mirt(data = VerbAgg3T_wide, model = IRTree.mod2, itemtype = "Rasch")
##
## Full-information item factor analysis with 2 factor(s).
## Converged within 1e-04 tolerance after 29 EM iterations.
## mirt version: 1.37.1
## M-step optimizer: nlminb
## EM acceleration: Ramsay
## Number of rectangular quadrature: 31
## Latent density type: Gaussian
##
## Log-likelihood = -6132.056
## Estimated parameters: 51
## AIC = 12366.11
## BIC = 12557.65; SABIC = 12395.9coef(tree.mirt2, simplify = TRUE)## $items
## a1 a2 d g u
## S1WantCurse:node1 1 0 1.225 0 1
## S1WantScold:node1 1 0 0.567 0 1
## S1WantShout:node1 1 0 0.082 0 1
## S2WantCurse:node1 1 0 1.753 0 1
## S2WantScold:node1 1 0 0.711 0 1
## S2WantShout:node1 1 0 0.014 0 1
## S3WantCurse:node1 1 0 0.532 0 1
## S3WantScold:node1 1 0 -0.685 0 1
## S3WantShout:node1 1 0 -1.528 0 1
## S4wantCurse:node1 1 0 1.085 0 1
## S4WantScold:node1 1 0 -0.348 0 1
## S4WantShout:node1 1 0 -1.043 0 1
## S1DoCurse:node1 1 0 1.225 0 1
## S1DoScold:node1 1 0 0.392 0 1
## S1DoShout:node1 1 0 -0.870 0 1
## S2DoCurse:node1 1 0 0.876 0 1
## S2DoScold:node1 1 0 -0.055 0 1
## S2DoShout:node1 1 0 -1.483 0 1
## S3DoCurse:node1 1 0 -0.209 0 1
## S3DoScold:node1 1 0 -1.505 0 1
## S3DoShout:node1 1 0 -2.985 0 1
## S4DoCurse:node1 1 0 0.711 0 1
## S4DoScold:node1 1 0 -0.383 0 1
## S4DoShout:node1 1 0 -2.002 0 1
## S1WantCurse:node2 0 1 0.379 0 1
## S1WantScold:node2 0 1 0.131 0 1
## S1WantShout:node2 0 1 -0.765 0 1
## S2WantCurse:node2 0 1 0.233 0 1
## S2WantScold:node2 0 1 0.019 0 1
## S2WantShout:node2 0 1 -0.388 0 1
## S3WantCurse:node2 0 1 -0.916 0 1
## S3WantScold:node2 0 1 -1.776 0 1
## S3WantShout:node2 0 1 -2.474 0 1
## S4wantCurse:node2 0 1 -0.540 0 1
## S4WantScold:node2 0 1 -1.075 0 1
## S4WantShout:node2 0 1 -1.089 0 1
## S1DoCurse:node2 0 1 -0.103 0 1
## S1DoScold:node2 0 1 -0.470 0 1
## S1DoShout:node2 0 1 -1.107 0 1
## S2DoCurse:node2 0 1 -0.150 0 1
## S2DoScold:node2 0 1 -0.875 0 1
## S2DoShout:node2 0 1 -1.493 0 1
## S3DoCurse:node2 0 1 -1.747 0 1
## S3DoScold:node2 0 1 -2.291 0 1
## S3DoShout:node2 0 1 -2.926 0 1
## S4DoCurse:node2 0 1 -0.775 0 1
## S4DoScold:node2 0 1 -1.308 0 1
## S4DoShout:node2 0 1 -1.797 0 1
##
## $means
## admission affirmation
## 0 0
##
## $cov
## admission affirmation
## admission 1.999 0.777
## affirmation 0.777 1.889Save fixed effects and random effects
d = coef(tree.mirt2, simplify = TRUE)$item[,'d']
fscores = fscores(tree.mirt2)Estimation: Generalized Linear Mixed Model
Data Transformation
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:mirt':
##
## fixefhead(VerbAgg3[ , -(1:2)])## S1WantCurse S1WantScold S1WantShout S2WantCurse S2WantScold S2WantShout
## [1,] 1 1 1 1 1 1
## [2,] 1 1 1 1 1 1
## [3,] 2 2 2 2 1 2
## [4,] 2 2 2 2 2 2
## [5,] 2 1 2 2 1 1
## [6,] 3 3 1 3 1 1
## S3WantCurse S3WantScold S3WantShout S4wantCurse S4WantScold S4WantShout
## [1,] 1 1 2 3 1 1
## [2,] 1 1 1 1 1 1
## [3,] 2 1 1 1 1 1
## [4,] 2 1 1 1 1 1
## [5,] 2 1 1 2 1 1
## [6,] 2 2 2 1 1 1
## S1DoCurse S1DoScold S1DoShout S2DoCurse S2DoScold S2DoShout S3DoCurse
## [1,] 2 1 2 2 1 1 2
## [2,] 1 1 1 1 1 1 2
## [3,] 1 2 2 1 1 2 1
## [4,] 2 2 2 3 2 2 2
## [5,] 2 2 1 2 1 1 2
## [6,] 3 1 1 3 1 1 3
## S3DoScold S3DoShout S4DoCurse S4DoScold S4DoShout
## [1,] 1 1 3 3 3
## [2,] 1 1 1 1 1
## [3,] 1 1 2 1 1
## [4,] 1 1 1 1 1
## [5,] 1 1 2 1 1
## [6,] 1 1 3 3 1VerbAgg3T <- dendrify(VerbAgg3[ , -(1:2)], mapping)
head(VerbAgg3T)## value item person node sub
## 1 0 i01 p001 node1 i01:node1
## 2 0 i01 p002 node1 i01:node1
## 3 1 i01 p003 node1 i01:node1
## 4 1 i01 p004 node1 i01:node1
## 5 1 i01 p005 node1 i01:node1
## 6 1 i01 p006 node1 i01:node1Unidimensional model for linear response trees
onedim <- glmer(value ~ 0 + item:node + (1 | person),
family = binomial, data = VerbAgg3T)## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0603612 (tol = 0.002, component 1)The model did not converge, find the optimization method that makes the model converge. (Note: This part is very time-consuming, you can just skip this.)
Model Convergence Check
library(optimx)
optimx_options <- c("L-BFGS-B", "nlminb", "nlm", "bobyqa", "nmkb", "hjkb")
for(i in 1:length(optimx_options)){
print(paste0("Testing optimx: ", optimx_options[i]))
onedim <- glmer(value ~ 0 + item:node + (1 | person),
family = binomial,
data = VerbAgg3T,
control = glmerControl(optimizer = "optimx",
optCtrl = list(method = optimx_options[i])))
if(is.null(onedim@optinfo$conv$lme4$messages)){
print(paste0("One of the optimx options, ", optimx_options[i],", worked!"))
print(onedim)
break
}
}## [1] "Testing optimx: L-BFGS-B"## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00281375 (tol = 0.002, component 1)## [1] "Testing optimx: nlminb"## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.## [1] "One of the optimx options, nlminb, worked!"
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: value ~ 0 + item:node + (1 | person)
## Data: VerbAgg3T
## AIC BIC logLik deviance df.resid
## 12778.183 13137.021 -6340.091 12680.183 11146
## Random effects:
## Groups Name Std.Dev.
## person (Intercept) 1.152
## Number of obs: 11195, groups: person, 316
## Fixed Effects:
## itemi01:nodenode1 itemi02:nodenode1 itemi03:nodenode1 itemi04:nodenode1
## 1.125634 0.509815 0.058646 1.625445
## itemi05:nodenode1 itemi06:nodenode1 itemi07:nodenode1 itemi08:nodenode1
## 0.643434 -0.004791 0.476854 -0.650117
## itemi09:nodenode1 itemi10:nodenode1 itemi11:nodenode1 itemi12:nodenode1
## -1.423994 0.994664 -0.339080 -0.979583
## itemi13:nodenode1 itemi14:nodenode1 itemi15:nodenode1 itemi16:nodenode1
## 1.125634 0.346456 -0.820467 0.797933
## itemi17:nodenode1 itemi18:nodenode1 itemi19:nodenode1 itemi20:nodenode1
## -0.068197 -1.382510 -0.211144 -1.403170
## itemi21:nodenode1 itemi22:nodenode1 itemi23:nodenode1 itemi24:nodenode1
## -2.764554 0.643434 -0.371269 -1.859269
## itemi01:nodenode2 itemi02:nodenode2 itemi03:nodenode2 itemi04:nodenode2
## 0.090175 -0.180945 -1.037889 0.022392
## itemi05:nodenode2 itemi06:nodenode2 itemi07:nodenode2 itemi08:nodenode2
## -0.251934 -0.646173 -1.041544 -2.036989
## itemi09:nodenode2 itemi10:nodenode2 itemi11:nodenode2 itemi12:nodenode2
## -2.600337 -0.702944 -1.307700 -1.356975
## itemi13:nodenode2 itemi14:nodenode2 itemi15:nodenode2 itemi16:nodenode2
## -0.254287 -0.713273 -1.350166 -0.268103
## itemi17:nodenode2 itemi18:nodenode2 itemi19:nodenode2 itemi20:nodenode2
## -1.078291 -1.778981 -1.815014 -2.443082
## itemi21:nodenode2 itemi22:nodenode2 itemi23:nodenode2 itemi24:nodenode2
## -3.306843 -0.855228 -1.465708 -2.196665
## optimizer (optimx) convergence code: 0 (OK) ; 1 optimizer warnings; 0 lme4 warningsThe ‘nlimnb’ optimization method was selected.
Final Model
onedim <- glmer(value ~ 0 + item:node + (1 | person),
family = binomial, data = VerbAgg3T,
control = glmerControl(optimizer = "optimx",
optCtrl = list(method = 'nlminb')))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.summary(onedim)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: value ~ 0 + item:node + (1 | person)
## Data: VerbAgg3T
## Control: glmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
##
## AIC BIC logLik deviance df.resid
## 12778.2 13137.0 -6340.1 12680.2 11146
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0383 -0.6902 -0.2646 0.7087 12.3788
##
## Random effects:
## Groups Name Variance Std.Dev.
## person (Intercept) 1.327 1.152
## Number of obs: 11195, groups: person, 316
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## itemi01:nodenode1 1.125634 0.152798 7.367 1.75e-13 ***
## itemi02:nodenode1 0.509815 0.143845 3.544 0.000394 ***
## itemi03:nodenode1 0.058646 0.141505 0.414 0.678549
## itemi04:nodenode1 1.625445 0.165658 9.812 < 2e-16 ***
## itemi05:nodenode1 0.643434 0.145199 4.431 9.36e-06 ***
## itemi06:nodenode1 -0.004791 0.141449 -0.034 0.972982
## itemi07:nodenode1 0.476854 0.143555 3.322 0.000895 ***
## itemi08:nodenode1 -0.650117 0.144707 -4.493 7.03e-06 ***
## itemi09:nodenode1 -1.423994 0.158730 -8.971 < 2e-16 ***
## itemi10:nodenode1 0.994664 0.150297 6.618 3.64e-11 ***
## itemi11:nodenode1 -0.339080 0.142254 -2.384 0.017143 *
## itemi12:nodenode1 -0.979583 0.149218 -6.565 5.21e-11 ***
## itemi13:nodenode1 1.125634 0.152797 7.367 1.75e-13 ***
## itemi14:nodenode1 0.346456 0.142596 2.430 0.015114 *
## itemi15:nodenode1 -0.820467 0.146793 -5.589 2.28e-08 ***
## itemi16:nodenode1 0.797933 0.147166 5.422 5.89e-08 ***
## itemi17:nodenode1 -0.068197 0.141447 -0.482 0.629711
## itemi18:nodenode1 -1.382510 0.157660 -8.769 < 2e-16 ***
## itemi19:nodenode1 -0.211144 0.141723 -1.490 0.136267
## itemi20:nodenode1 -1.403170 0.158174 -8.871 < 2e-16 ***
## itemi21:nodenode1 -2.764554 0.219124 -12.616 < 2e-16 ***
## itemi22:nodenode1 0.643434 0.145200 4.431 9.36e-06 ***
## itemi23:nodenode1 -0.371269 0.142431 -2.607 0.009143 **
## itemi24:nodenode1 -1.859269 0.172467 -10.780 < 2e-16 ***
## itemi01:nodenode2 0.090175 0.162569 0.555 0.579109
## itemi02:nodenode2 -0.180945 0.173525 -1.043 0.297058
## itemi03:nodenode2 -1.037889 0.190122 -5.459 4.79e-08 ***
## itemi04:nodenode2 0.022392 0.155246 0.144 0.885313
## itemi05:nodenode2 -0.251934 0.169962 -1.482 0.138262
## itemi06:nodenode2 -0.646173 0.187698 -3.443 0.000576 ***
## itemi07:nodenode2 -1.041544 0.180550 -5.769 7.99e-09 ***
## itemi08:nodenode2 -2.036989 0.246678 -8.258 < 2e-16 ***
## itemi09:nodenode2 -2.600337 0.342643 -7.589 3.22e-14 ***
## itemi10:nodenode2 -0.702944 0.165347 -4.251 2.12e-05 ***
## itemi11:nodenode2 -1.307700 0.206323 -6.338 2.33e-10 ***
## itemi12:nodenode2 -1.356975 0.245687 -5.523 3.33e-08 ***
## itemi13:nodenode2 -0.254287 0.159865 -1.591 0.111692
## itemi14:nodenode2 -0.713273 0.176387 -4.044 5.26e-05 ***
## itemi15:nodenode2 -1.350166 0.230327 -5.862 4.57e-09 ***
## itemi16:nodenode2 -0.268103 0.165626 -1.619 0.105508
## itemi17:nodenode2 -1.078291 0.191514 -5.630 1.80e-08 ***
## itemi18:nodenode2 -1.778981 0.275793 -6.450 1.12e-10 ***
## itemi19:nodenode2 -1.815014 0.218907 -8.291 < 2e-16 ***
## itemi20:nodenode2 -2.443082 0.317403 -7.697 1.39e-14 ***
## itemi21:nodenode2 -3.306843 0.625544 -5.286 1.25e-07 ***
## itemi22:nodenode2 -0.855228 0.171073 -4.999 5.76e-07 ***
## itemi23:nodenode2 -1.465708 0.212534 -6.896 5.34e-12 ***
## itemi24:nodenode2 -2.196665 0.352779 -6.227 4.76e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 48 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (optimx) convergence code: 0 (OK)
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.head(ranef(onedim)$person)## (Intercept)
## p001 -0.0820955
## p002 -2.2716979
## p003 -0.6405218
## p004 -0.1311874
## p005 -0.6605571
## p006 0.5415838Multidimensional model for linear response trees
twodim <- glmer(value ~ 0 + item:node + (0 + node | person),
family = binomial, data = VerbAgg3T,
control = glmerControl(optimizer = "optimx",
optCtrl = list(method = 'nlminb')))
summary(twodim)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: value ~ 0 + item:node + (0 + node | person)
## Data: VerbAgg3T
## Control: glmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
##
## AIC BIC logLik deviance df.resid
## 12377.9 12751.3 -6137.9 12275.9 11144
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4803 -0.6289 -0.2239 0.6480 15.5848
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## person nodenode1 1.918 1.385
## nodenode2 1.821 1.349 0.39
## Number of obs: 11195, groups: person, 316
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## itemi01:nodenode1 1.22222 0.16259 7.517 5.59e-14 ***
## itemi02:nodenode1 0.56602 0.15400 3.676 0.000237 ***
## itemi03:nodenode1 0.08154 0.15192 0.537 0.591456
## itemi04:nodenode1 1.75012 0.17514 9.993 < 2e-16 ***
## itemi05:nodenode1 0.70893 0.15527 4.566 4.98e-06 ***
## itemi06:nodenode1 0.01319 0.15190 0.087 0.930819
## itemi07:nodenode1 0.53072 0.15374 3.452 0.000556 ***
## itemi08:nodenode1 -0.68466 0.15557 -4.401 1.08e-05 ***
## itemi09:nodenode1 -1.52558 0.17031 -8.958 < 2e-16 ***
## itemi10:nodenode1 1.08319 0.16016 6.763 1.35e-11 ***
## itemi11:nodenode1 -0.34777 0.15291 -2.274 0.022943 *
## itemi12:nodenode1 -1.04235 0.16033 -6.501 7.97e-11 ***
## itemi13:nodenode1 1.22222 0.16258 7.517 5.58e-14 ***
## itemi14:nodenode1 0.39094 0.15285 2.558 0.010538 *
## itemi15:nodenode1 -0.86951 0.15778 -5.511 3.57e-08 ***
## itemi16:nodenode1 0.87380 0.15715 5.560 2.69e-08 ***
## itemi17:nodenode1 -0.05518 0.15193 -0.363 0.716467
## itemi18:nodenode1 -1.48046 0.16920 -8.750 < 2e-16 ***
## itemi19:nodenode1 -0.20948 0.15229 -1.376 0.168971
## itemi20:nodenode1 -1.50293 0.16975 -8.854 < 2e-16 ***
## itemi21:nodenode1 -2.97858 0.23304 -12.781 < 2e-16 ***
## itemi22:nodenode1 0.70893 0.15527 4.566 4.98e-06 ***
## itemi23:nodenode1 -0.38260 0.15310 -2.499 0.012453 *
## itemi24:nodenode1 -1.99894 0.18467 -10.824 < 2e-16 ***
## itemi01:nodenode2 0.37923 0.17596 2.155 0.031144 *
## itemi02:nodenode2 0.13281 0.18752 0.708 0.478787
## itemi03:nodenode2 -0.76372 0.20414 -3.741 0.000183 ***
## itemi04:nodenode2 0.23378 0.16735 1.397 0.162427
## itemi05:nodenode2 0.01996 0.18453 0.108 0.913848
## itemi06:nodenode2 -0.38551 0.20367 -1.893 0.058385 .
## itemi07:nodenode2 -0.91474 0.19324 -4.734 2.20e-06 ***
## itemi08:nodenode2 -1.77152 0.26231 -6.753 1.44e-11 ***
## itemi09:nodenode2 -2.46670 0.36472 -6.763 1.35e-11 ***
## itemi10:nodenode2 -0.53976 0.17741 -3.042 0.002346 **
## itemi11:nodenode2 -1.07187 0.22295 -4.808 1.53e-06 ***
## itemi12:nodenode2 -1.08619 0.26004 -4.177 2.95e-05 ***
## itemi13:nodenode2 -0.10062 0.17393 -0.578 0.562928
## itemi14:nodenode2 -0.46716 0.19214 -2.431 0.015042 *
## itemi15:nodenode2 -1.10319 0.24630 -4.479 7.49e-06 ***
## itemi16:nodenode2 -0.14681 0.18014 -0.815 0.415081
## itemi17:nodenode2 -0.87163 0.20783 -4.194 2.74e-05 ***
## itemi18:nodenode2 -1.48876 0.29502 -5.046 4.50e-07 ***
## itemi19:nodenode2 -1.74260 0.23488 -7.419 1.18e-13 ***
## itemi20:nodenode2 -2.28354 0.33932 -6.730 1.70e-11 ***
## itemi21:nodenode2 -2.91340 0.67215 -4.334 1.46e-05 ***
## itemi22:nodenode2 -0.77186 0.18594 -4.151 3.31e-05 ***
## itemi23:nodenode2 -1.30416 0.22977 -5.676 1.38e-08 ***
## itemi24:nodenode2 -1.79234 0.37381 -4.795 1.63e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 48 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itfixef = summary(twodim)$coef[, 1, drop = FALSE]
fixef## Estimate
## itemi01:nodenode1 1.22222226
## itemi02:nodenode1 0.56601507
## itemi03:nodenode1 0.08153819
## itemi04:nodenode1 1.75012156
## itemi05:nodenode1 0.70893028
## itemi06:nodenode1 0.01318707
## itemi07:nodenode1 0.53072477
## itemi08:nodenode1 -0.68466448
## itemi09:nodenode1 -1.52557943
## itemi10:nodenode1 1.08319021
## itemi11:nodenode1 -0.34776969
## itemi12:nodenode1 -1.04234670
## itemi13:nodenode1 1.22222226
## itemi14:nodenode1 0.39093537
## itemi15:nodenode1 -0.86951483
## itemi16:nodenode1 0.87380107
## itemi17:nodenode1 -0.05517909
## itemi18:nodenode1 -1.48045581
## itemi19:nodenode1 -0.20948239
## itemi20:nodenode1 -1.50292781
## itemi21:nodenode1 -2.97857579
## itemi22:nodenode1 0.70893029
## itemi23:nodenode1 -0.38259936
## itemi24:nodenode1 -1.99894232
## itemi01:nodenode2 0.37922572
## itemi02:nodenode2 0.13280952
## itemi03:nodenode2 -0.76372490
## itemi04:nodenode2 0.23377898
## itemi05:nodenode2 0.01996324
## itemi06:nodenode2 -0.38551339
## itemi07:nodenode2 -0.91474072
## itemi08:nodenode2 -1.77152492
## itemi09:nodenode2 -2.46670291
## itemi10:nodenode2 -0.53975786
## itemi11:nodenode2 -1.07187062
## itemi12:nodenode2 -1.08618730
## itemi13:nodenode2 -0.10061840
## itemi14:nodenode2 -0.46715929
## itemi15:nodenode2 -1.10319060
## itemi16:nodenode2 -0.14680834
## itemi17:nodenode2 -0.87163360
## itemi18:nodenode2 -1.48875873
## itemi19:nodenode2 -1.74260170
## itemi20:nodenode2 -2.28353590
## itemi21:nodenode2 -2.91339825
## itemi22:nodenode2 -0.77185878
## itemi23:nodenode2 -1.30416130
## itemi24:nodenode2 -1.79233904head(ranef(twodim)$person)## nodenode1 nodenode2
## p001 -0.43165219 0.592751
## p002 -2.51991853 -1.051173
## p003 -0.35202598 -1.756556
## p004 0.41824318 -1.401137
## p005 -0.35274937 -1.771928
## p006 -0.01440857 1.448797ranef = ranef(twodim)$personRelationships between the estimates from MIRT and GLMM
par(mfrow=c(1,3),mgp=c(2,1,0),mar=c(3,3,1,1))
plot(fixef, d) # item:node
plot(ranef[,1], fscores[,1]) # latent trait 1
plot(ranef[,2], fscores[,2]) # latent trait 2