tam.mml: Test Analysis Modules: Marginal Maximum Likelihood Estimation

Description Usage Arguments Details Value Note References See Also Examples

View source: R/tam.mml.R

Description

Modules for psychometric test analysis demonstrated with the help of artificial example data. The package includes MML and JML estimation of uni- and multidimensional IRT (Rasch, 2PL, Generalized Partial Credit, Rating Scale, Multi Facets, Nominal Item Response) models, fit statistic computation, standard error estimation, as well as plausible value imputation and weighted likelihood estimation of ability.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
tam(resp, irtmodel="1PL", formulaA=NULL, ...)

tam.mml(resp, Y=NULL, group=NULL, irtmodel="1PL", formulaY=NULL,
    dataY=NULL, ndim=1,  pid=NULL, xsi.fixed=NULL, xsi.inits=NULL,
    beta.fixed=NULL, beta.inits=NULL, variance.fixed=NULL,
    variance.inits=NULL, est.variance=TRUE, constraint="cases", A=NULL,
    B=NULL, B.fixed=NULL,  Q=NULL,  est.slopegroups=NULL, E=NULL,
    pweights=NULL, userfct.variance=NULL,
    variance.Npars=NULL, item.elim=TRUE, verbose=TRUE, control=list() )

tam.mml.2pl(resp, Y=NULL, group=NULL, irtmodel="2PL", formulaY=NULL,
    dataY=NULL, ndim=1, pid=NULL, xsi.fixed=NULL, xsi.inits=NULL,
    beta.fixed=NULL, beta.inits=NULL, variance.fixed=NULL,
    variance.inits=NULL, est.variance=FALSE, A=NULL, B=NULL,
    B.fixed=NULL, Q=NULL, est.slopegroups=NULL, E=NULL, gamma.init=NULL,
    pweights=NULL, userfct.variance=NULL, variance.Npars=NULL,
    item.elim=TRUE, verbose=TRUE, control=list() )

tam.mml.mfr(resp, Y=NULL, group=NULL, irtmodel="1PL", formulaY=NULL,
    dataY=NULL, ndim=1, pid=NULL, xsi.fixed=NULL, xsi.setnull=NULL,
    xsi.inits=NULL, beta.fixed=NULL, beta.inits=NULL, variance.fixed=NULL,
    variance.inits=NULL, est.variance=TRUE, formulaA=~item+item:step,
    constraint="cases", A=NULL, B=NULL,  B.fixed=NULL, Q=NULL,
    facets=NULL, est.slopegroups=NULL, E=NULL,
    pweights=NULL, verbose=TRUE, control=list(), delete.red.items=TRUE )

## S3 method for class 'tam'
summary(object, file=NULL, ...)

## S3 method for class 'tam.mml'
summary(object, file=NULL, ...)

## S3 method for class 'tam'
print(x, ...)

## S3 method for class 'tam.mml'
print(x, ...)

Arguments

resp

Data frame with polytomous item responses k=0,...,K. Missing responses must be declared as NA.

Y

A matrix of covariates in latent regression. Note that the intercept is automatically included as the first predictor.

group

An optional vector of group identifiers

irtmodel

For fixed item slopes (in tam.mml) options include PCM (partial credit model), PCM2 (partial credit model with ConQuest parametrization 'item+item*step' and RSM (rating scale model; the ConQuest parametrization 'item+step').
For estimated item slopes (only available in tam.mml.2pl) options are 2PL (all slopes of item categories are estimated; Nominal Item Response Model), GPCM (generalized partial credit model in which every item gets one and only slope parameter per dimension) and 2PL.groups (subsets of items get same item slope estimates) and a design matrix E on item slopes in the generalized partial credit model (GPCM.design, see Examples). Note that item slopes can not be estimated with faceted designs using the function tam.mml.mfr. However, it is easy to use pre-specified design matrices and apply some restrictions to tam.mml.2pl (see Example 14, Model 3).

formulaY

An R formula for latent regression. Transformations of predictors in Y (included in dataY) can be easily specified, e. g. female*race or I(age^2).

dataY

An optional data frame with possible covariates Y in latent regression. This data frame is used if an R formula in formulaY is specified.

ndim

Number of dimensions (is not needed to determined by the user)

pid

An optional vector of person identifiers

xsi.fixed

A matrix with two columns for fixing ξ parameters. 1st column: index of ξ parameter, 2nd column: fixed value

xsi.setnull

A vector of strings indicating which ξ elements should be set to zero which do have entries in xsi.setnull in their labels (see Example 10a).

xsi.inits

A matrix with two columns (in the same way defined as in xsi.fixed with initial value for ξ.

beta.fixed

A matrix with three columns for fixing regression coefficients. 1st column: Index of Y value, 2nd column: dimension, 3rd column: fixed β value.
If no constraints should be imposed on β, then set beta.fixed=FALSE (see Example 2, Model 2_4).

beta.inits

A matrix (same format as in beta.fixed) with initial β values

variance.fixed

An optional matrix with three columns for fixing entries in covariance matrix: 1st column: dimension 1, 2nd column: dimension 2, 3rd column: fixed value

variance.inits

Initial covariance matrix in estimation. All matrix entries have to be specified and this matrix is NOT in the same format like variance.fixed.

est.variance

Should the covariance matrix be estimated? This argument applies to estimated item slopes in tam.mml.2pl. The default is FALSE which means that latent variables (in the first group) are standardized in 2PL estimation.

constraint

Set sum constraint for parameter identification for items or cases (applies to tam.mml and tam.mml.mfr)

A

An optional array of dimension I \times (K+1) \times N_ξ. Only ξ parameters are estimated, entries in A only correspond to the design.

B

An optional array of dimension I \times (K+1) \times D. In case of tam.mml.2pl entries of the B matrix can be estimated.

B.fixed

An optional matrix with four columns for fixing B matrix entries in 2PL estimation. 1st column: item index, 2nd column: category, 3rd column: dimension, 4th column: fixed value.

Q

An optional I \times D matrix (the Q-matrix) which specifies the loading structure of items on dimensions.

est.slopegroups

A vector of integers of length I for estimating item slope parameters of item groups. This function only applies to the generalized partial credit model
(irtmodel="2PL.groups").

E

An optional design matrix for estimating item slopes in the generalized partial credit model (irtmodel="GPCM.design")

gamma.init

Optional initial gamma parameter vector (irtmodel="GPCM.design").

pweights

An optional vector of person weights

formulaA

Design formula (only applies to tam.mml.mfr). See Example 8. It is also to possible to set all effects of a facet to zero, e.g. item*step + 0*rater (see Example 10a).

facets

A data frame with facet entries (only applies to tam.mml.mfr)

userfct.variance

Optional user customized function for variance specification (See Simulated Example 17).

variance.Npars

Number of estimated parameters of variance matrix if a userfct.variance is provided.

item.elim

Optional logical indicating whether an item with has only zero entries should be removed from the analysis. The default is TRUE.

verbose

Logical indicating whether output should be printed during iterations. This argument replaces control$progress.

control

A list of control arguments for the algorithm:

list( nodes=seq(-6,6,len=21), snodes=0, QMC=TRUE,
convD=.001,conv=.0001, convM=.0001, Msteps=4,
maxiter=1000, max.increment=1,
min.variance=.001, progress=TRUE, ridge=0,
seed=NULL, xsi.start0=0, increment.factor=1,
fac.oldxsi=0, acceleration="none", dev_crit="absolute" )
trim_increment="half" )

nodes: the discretized θ nodes for numerical integration

snodes: number of simulated θ nodes for stochastic integration. If snodes=0, numerical integration is used.

QMC: A logical indicating whether quasi Monte Carlo integration (Gonzales at al., 2006; Pan & Thompson, 2007) should be used. The default is TRUE. Quasi Monte Carlo integration is a nonstochastic integration approach but prevents the very demanding numeric integration using Gaussian quadrature. In case of QMC=FALSE, "ordinary" stochastic integration is used (see the section Integration in Details).

convD: Convergence criterion for deviance

conv: Convergence criterion for item parameters and regression coefficients

convM: Convergence criterion for item parameters within each M step

Msteps: Number of M steps for item parameter estimation. A high value of M steps could be helpful in cases of non-convergence. In tam.mml, tam.mml.2pl and tam.mml.mfr, the default is set to 4, in tam.mml.3pl it is set to 10.

maxiter: Maximum number of iterations

max.increment: Maximum increment for item parameter change for every iteration

min.variance: Minimum variance to be estimated during iterations.

progress: A logical indicating whether computation progress should be displayed at R console

ridge: A numeric value or a vector of ridge parameter(s) for the latent regression which is added to the covariance matrix Y'Y of predictors in the diagonal.

seed: An optional integer defining the simulation seed (important for reproducible results for stochastic integration)

xsi.start0: A numeric value. The value of 0 indicates that for all parameters starting values are provided. A value of 1 means that all starting values are set to zero and a value of 2 means that only starting values of item parameters (but not facet parameters) are used.

increment.factor: A value (larger than one) which defines the extent of the decrease of the maximum increment of item parameters in every iteration. The maximum increment in iteration iter is defined as max.increment*increment.factor^(-iter) where max.increment=1. Using a value larger than 1 helps to reach convergence in some non-converging analyses (see Example 12).

fac.oldxsi: An optional numeric value f between 0 and 1 which defines the weight of parameter values in previous iteration. If ξ_t denotes a parameter update in iteration t, ξ_{t-1} is the parameter value of iteration t-1, then the modified parameter value is defined as ξ_t^*=(1-f) \cdot ξ_t + f \cdot ξ_{t-1}. Especially in cases where the deviance increases, setting the parameter larger than 0 (maybe .4 or .5) is helpful in stabilizing the algorithm (see Example 15).

acceleration: String indicating whether convergence acceleration of the EM algorithm should be employed. Options are "none" (no acceleration, the default), the monotone overrelaxation method of "Yu" (Yu, 2012) and "Ramsay" for the Ramsay (1975) acceleration method.

dev_crit: Criterion for convergence in deviance. dev_crit="absolute" refers to absolute differences in successive deviance values, while dev_crit="relative" refers to relative differences.

trim_increment: Type of method for trimming parameter increments in algorithm. Possible types are "half" or ""cut".

delete.red.items

An optional logical indicating whether redundant generalized items (with no observations) should be eliminated.

object

Object of class tam or tam.mml (only for summary.tam functions)

file

A file name in which the summary output should be written (only for summary.tam functions)

...

Further arguments to be passed

x

Object of class tam or tam.mml

Details

The multidimensional item response model in TAM is described in Adams, Wilson and Wu (1997) or Adams and Wu (2007).

The data frame resp contains item responses of N persons (in rows) at I items (in columns), each item having at most K categories k=0,...,K. The item response model has D dimensions of the θ ability vector and can be written as

P( X_{pi}=k | θ_p ) \propto exp( b_{ik} θ_p + a_{ik} ξ )

The symbol \propto means that response probabilities are normalized such that ∑ _k P( X_{pi}=k | θ_p )=1 .

Item category thresholds for item i in category k are written as a linear combination a_i ξ where the vector ξ of length N_ξ contains generalized item parameters and A=( a_{ik} )_{ik}=( a_i )_{i} is a three-dimensional design array (specified in A).

The scoring vector b_{ik} contains the fixed (in tam.mml) or estimated (in tam.mml.2pl) scores of item i in category k on dimension d.

For tam.mml.2pl and irtmodel="GPCM.design", item slopes a_i can be written as a linear combination a_i=( E γ)_i of basis item slopes which is an analogue of the LLTM for item slopes (see Example 7; Embretson, 1999).

The latent regression model regresses the latent trait θ_p on covariates Y which results in

θ_p=Y β + ε_p, ε_p ~ N_D ( 0, Σ )

Where β is a N_Y times D matrix of regression coefficients for N_Y covariates in Y.

The multiple group model for groups g=1,...,G is implemented for unidimensional and multidimensional item response models. In this case, variance heterogeneity is allowed

θ_p=Y β + ε_p, ε_p ~ N( 0, σ_g^2 )

Integration: Uni- and multidimensional integrals are approximated by posing a uni- or multivariate normality assumption. The default is Gaussian quadrature with nodes defined in control$nodes. For D-dimensional IRT models, the D-dimensional cube consisting of the vector control$nodes in all dimensions is used. If the user specifies control$snodes with a value larger than zero, then Quasi-Monte Carlo integration (Pan & Thomas, 2007; Gonzales et al., 2006) with control$snodes is used (because control$QMC=TRUE is set by default). If control$QMC=FALSE is specified, then stochastic (Monte Carlo) integration is employed with control$snodes stochastic nodes.

Value

A list with following entries:

xsi

Vector of ξ parameter estimates and their corresponding standard errors

xsi.facets

Data frame of ξ parameters and corresponding constraints for multifacet models

beta

Matrix of β regression coefficient estimates

variance

Covariance matrix. In case of multiple groups, it is a vector indicating heteroscedastic variances

item

Data frame with item parameters. The column xsi.item denotes the item difficulty of polytomous items in the parametrization irtmodel="PCM2".

item_irt

IRT parameterization of item parameters

person

Matrix with person parameter estimates. EAP is the mean of the posterior distribution and SD.EAP the corresponding standard deviation

pid

Vector of person identifiers

EAP.rel

EAP reliability

post

Posterior distribution for item response pattern

rprobs

A three-dimensional array with estimated response probabilities (dimensions are items \times categories \times theta length)

itemweight

Matrix of item weights

theta

Theta grid

n.ik

Array of expected counts: theta class \times item \times category \times group

pi.k

Marginal trait distribution at grid theta

Y

Matrix of covariates

resp

Original data frame

resp.ind

Response indicator matrix

group

Group identifier

G

Number of groups

formulaY

Formula for latent regression

dataY

Data frame for latent regression

pweights

Person weights

time

Computation time

A

Design matrix A for ξ parameters

B

Fixed or estimated loading matrix

se.B

Standard errors of B parameters

nitems

Number of items

maxK

Maximum number of categories

AXsi

Estimated item intercepts a_{ik} ξ

AXsi_

Estimated item intercepts - a_{ik} ξ. Note that in summary.tam, the parameters AXsi_ are displayed.

se.AXsi

Standard errors of a_{ik} ξ parameters

nstud

Number of persons

resp.ind.list

List of response indicator vectors

hwt

Individual posterior distribution

like

Individual likelihood

ndim

Number of dimensions

xsi.fixed

Fixed ξ parameters

xsi.fixed.estimated

Matrix of estimated ξ parameters in form of xsi.fixed which can be used for parameter fixing in subsequent estimations.

B.fixed

Fixed loading parameters (only applies to tam.mml.2pl)

B.fixed.estimated

Matrix of estimated B parameters in the same format as B.fixed.

est.slopegroups

An index vector of item groups of common slope parameters (only applies to tam.mml.2pl)

E

Design matrix for estimated item slopes in the generalized partial credit model (only applies to tam.mml.2pl in case of irtmodel="GPCM.design")

basispar

Vector of γ parameters of the linear combination a_i=( E γ)_i for item slopes (only applies to tam.mml.2pl in case of irtmodel='GPCM.design')

formulaA

Design formula (only applies to tam.mml.mfr)

facets

Data frame with facet entries (only applies to tam.mml.mfr)

variance.fixed

Fixed covariance matrix

nnodes

Number of theta nodes

deviance

Final deviance

ic

Vector with information criteria

deviance.history

Deviance history in iterations

control

List of control arguments

latreg_stand

List containing standardized regression coefficients

...

Further values

Note

For more than three dimensions, quasi-Monte Carlo or stochastic integration is recommended because otherwise problems in memory allocation and large computation time will result. Choose in control a suitable value of the number of quasi Monte Carlo or stochastic nodes snodes (say, larger than 1000). See Pan and Thompson (2007) or Gonzales et al. (2006) for more details.

In faceted models (tam.mml.mfr), parameters which cannot be estimated are fixed to 99.

Several choices can be made if your model does not converge. First, the number of iterations within a M step can be increased (Msteps=10). Second, the absolute value of increments can be forced with increasing iterations (set a value higher than 1 to max.increment, maybe 1.05). Third, change in estimated parameters can be stabilized by fac.oldxsi for which a value of 0 (the default) and a value of 1 can be chosen. We recommend values between .5 and .8 if your model does not converge.

References

Adams, R. J., Wilson, M., & Wu, M. (1997). Multilevel item response models: An approach to errors in variables regression. Journal of Educational and Behavioral Statistics, 22, 47-76. doi: 10.3102/10769986022001047

Adams, R. J., & Wu, M. L. (2007). The mixed-coefficients multinomial logit model. A generalized form of the Rasch model. In M. von Davier & C. H. Carstensen (Eds.), Multivariate and mixture distribution Rasch models: Extensions and applications (pp. 55-76). New York: Springer. doi: 10.1007/978-0-387-49839-3_4

Embretson, S. E. (1999). Generating items during testing: Psychometric issues and models. Psychometrika, 64, 407-433. doi: 10.1007/BF02294564

Gonzalez, J., Tuerlinckx, F., De Boeck, P., & Cools, R. (2006). Numerical integration in logistic-normal models. Computational Statistics & Data Analysis, 51, 1535-1548. doi: 10.1016/j.csda.2006.05.003

Pan, J., & Thompson, R. (2007). Quasi-Monte Carlo estimation in generalized linear mixed models. Computational Statistics & Data Analysis, 51, 5765-5775. doi: 10.1016/j.csda.2006.10.003

Ramsay, J. O. (1975). Solving implicit equations in psychometric data analysis. Psychometrika, 40(3), 337-360. doi: 10.1007/BF02291762

Yu, Y. (2012). Monotonically overrelaxed EM algorithms. Journal of Computational and Graphical Statistics, 21(2), 518-537. doi: 10.1080/10618600.2012.672115

Wu, M. L., Adams, R. J., Wilson, M. R. & Haldane, S. (2007). ACER ConQuest Version 2.0. Mulgrave. https://shop.acer.edu.au/acer-shop/group/CON3.

See Also

See data.cqc01 for more examples which is are similar to the ones in the ConQuest manual (Wu, Adams, Wilson & Haldane, 2007).

See tam.jml for joint maximum likelihood estimation.

Standard errors are estimated by a rather crude (but quick) approximation. Use tam.se for improved standard errors.

For model comparisons see anova.tam.

See sirt::tam2mirt for converting tam objects into objects of class mirt::mirt in the mirt package.

Examples

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
#############################################################################
# EXAMPLE 1: dichotomous data
# data.sim.rasch: 2000 persons, 40 items
#############################################################################
data(data.sim.rasch)

#************************************************************
# Model 1: Rasch model (MML estimation)
mod1 <- TAM::tam.mml(resp=data.sim.rasch)
# extract item parameters
mod1$item    # item difficulties

## Not run: 
# WLE estimation
wle1 <- TAM::tam.wle( mod1 )
# item fit
fit1 <- TAM::tam.fit(mod1)
# plausible value imputation
pv1 <- TAM::tam.pv(mod1, normal.approx=TRUE, ntheta=300)
# standard errors
se1 <- TAM::tam.se( mod1 )
# summary
summary(mod1)

#-- specification with tamaan
tammodel <- "
 LAVAAN MODEL:
   F=~ I1__I40;
   F ~~ F
 ITEM TYPE:
   ALL(Rasch)
   "
mod1t <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod1t)

#************************************************************
# Model 1a: Rasch model with fixed item difficulties from 'mod1'
xsi0 <- mod1$xsi$xsi
xsi.fixed <- cbind( 1:(length(xsi0)), xsi0 )
        # define vector with fixed item difficulties
mod1a <- TAM::tam.mml( resp=data.sim.rasch, xsi.fixed=xsi.fixed )
summary(mod1a)

# Usage of the output value mod1$xsi.fixed.estimated has the right format
# as the input of xsi.fixed
mod1aa <- TAM::tam.mml( resp=data.sim.rasch, xsi.fixed=mod1$xsi.fixed.estimated )
summary(mod1b)

#************************************************************
# Model 1b: Rasch model with initial xsi parameters for items 2 (item difficulty b=-1.8),
# item 4 (b=-1.6) and item 40 (b=2)
xsi.inits <- cbind( c(2,4,40), c(-1.8,-1.6,2))
mod1b <- TAM::tam.mml( resp=data.sim.rasch, xsi.inits=xsi.inits )

#-- tamaan specification
tammodel <- "
 LAVAAN MODEL:
   F=~ I1__I40
   F ~~ F
   # Fix item difficulties. Note that item intercepts instead of difficulties
   # must be specified.
   I2 | 1.8*t1
   I4 | 1.6*t1
 ITEM TYPE:
   ALL(Rasch)
   "
mod1bt <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod1bt)

#************************************************************
# Model 1c: 1PL estimation with sum constraint on item difficulties
dat <- data.sim.rasch
# modify A design matrix to include the sum constraint
des <- TAM::designMatrices(resp=dat)
A1 <- des$A[,, - ncol(dat) ]
A1[ ncol(dat),2, ] <- 1
A1[,2,]
# estimate model
mod1c <- TAM::tam.mml( resp=dat, A=A1, beta.fixed=FALSE,
           control=list(fac.oldxsi=.1) )
summary(mod1c)

#************************************************************
# Model 1d: estimate constraint='items' using tam.mml.mfr
formulaA=~ 0 + item
mod1d <- TAM::tam.mml.mfr( resp=dat, formulaA=formulaA,
                     control=list(fac.oldxsi=.1), constraint="items")
summary(mod1d)

#************************************************************
# Model 1e: This sum constraint can also be obtained by using the argument
# constraint="items" in tam.mml
mod1e <- TAM::tam.mml( resp=data.sim.rasch, constraint="items" )
summary(mod1e)

#************************************************************
# Model 1d2: estimate constraint='items' using tam.mml.mfr
# long format response data
resp.long <- c(dat)

# pid and item facet specifications are necessary
#     Note, that we recommend the facet labels to be sortable in the same order that the
#     results are desired.
#     compare to: facets <- data.frame( "item"=rep(colnames(dat), each=nrow(dat)) )
pid <- rep(1:nrow(dat), ncol(dat))
itemnames <- paste0("I", sprintf(paste('%0', max(nchar(1:ncol(dat))), 'i', sep='' ),
                    c(1:ncol(dat)) ) )
facets   <- data.frame( "item_"=rep(itemnames, each=nrow(dat)) )
formulaA=~ 0 + item_

mod1d2 <- TAM::tam.mml.mfr( resp=resp.long, formulaA=formulaA, control=list(fac.oldxsi=.1),
                       constraint="items", facets=facets, pid=pid)
stopifnot( all(mod1d$xsi.facets$xsi==mod1d2$xsi.facets$xsi) )

## End(Not run)


#************************************************************
# Model 2: 2PL model
mod2 <- TAM::tam.mml.2pl(resp=data.sim.rasch,irtmodel="2PL")

# extract item parameters
mod2$xsi    # item difficulties
mod2$B      # item slopes

#--- tamaan specification
tammodel <- "
 LAVAAN MODEL:
   F=~ I1__I40
   F ~~ 1*F
   # item type of 2PL is the default for dichotomous data
   "
# estimate model
mod2t <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod2t)

## Not run: 
#************************************************************
# Model 2a: 2PL with fixed item difficulties and slopes from 'mod2'
xsi0 <- mod2$xsi$xsi
xsi.fixed <- cbind( 1:(length(xsi0)), xsi0 )
        # define vector with fixed item difficulties
mod2a <- TAM::tam.mml( resp=data.sim.rasch, xsi.fixed=xsi.fixed,
                 B=mod2$B # fix slopes
            )
summary(mod2a)
mod2a$B     # inspect used slope matrix

#************************************************************
# Model 3: constrained 2PL estimation
# estimate item parameters in different slope groups
# items 1-10, 21-30 group 1
# items 11-20 group 2 and items 31-40 group 3
est.slope <- rep(1,40)
est.slope[ 11:20 ] <- 2
est.slope[ 31:40 ] <- 3
mod3 <- TAM::tam.mml.2pl( resp=data.sim.rasch, irtmodel="2PL.groups",
               est.slopegroups=est.slope )
mod3$B
summary(mod3)

#--- tamaan specification (A)
tammodel <- "
 LAVAAN MODEL:
   F=~ lam1*I1__I10 + lam2*I11__I20 + lam1*I21__I30 + lam3*I31__I40;
   F ~~ 1*F
   "
# estimate model
mod3tA <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod3tA)

#--- tamaan specification (alternative B)
tammodel <- "
 LAVAAN MODEL:
   F=~ a1__a40*I1__I40;
   F ~~ 1*F
 MODEL CONSTRAINT:
   a1__a10==lam1
   a11__a20==lam2
   a21__a30==lam1
   a31__a40==lam3
   "
mod3tB <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod3tB)

#--- tamaan specification (alternative C using DO operator)
tammodel <- "
 LAVAAN MODEL:
 DO(1,10,1)
   F=~ lam1*I%
 DOEND
 DO(11,20,1)
   F=~ lam2*I%
 DOEND
 DO(21,30,1)
   F=~ lam1*I%
 DOEND
 DO(31,40,1)
   F=~ lam3*I%
 DOEND
   F ~~ 1*F
   "
# estimate model
mod3tC <- TAM::tamaan( tammodel, data.sim.rasch)
summary(mod3tC)

#############################################################################
# EXAMPLE 2: Unidimensional calibration with latent regressors
#############################################################################

# (1) simulate data
set.seed(6778)     # set simulation seed
N <- 2000          # number of persons
# latent regressors Y
Y <- cbind( stats::rnorm( N, sd=1.5), stats::rnorm(N, sd=.3 ) )
# simulate theta
theta <- stats::rnorm( N ) + .4 * Y[,1] + .2 * Y[,2]  # latent regression model
# number of items
I <- 40
p1 <- stats::plogis( outer( theta, seq( -2, 2, len=I ), "-" ) )
# simulate response matrix
resp <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
colnames(resp) <- paste("I", 1:I, sep="")

# (2) estimate model
mod2_1 <- TAM::tam.mml(resp=resp, Y=Y)
summary(mod2_1)

# (3) setting initial values for beta coefficients
#   beta_2=.20, beta_3=.35 for dimension 1
beta.inits <- cbind( c(2,3), 1, c(.2, .35 ) )
mod2_2 <- TAM::tam.mml(resp=resp, Y=Y, beta.inits=beta.inits)

# (4) fix intercept to zero and third coefficient to .3
beta.fixed <- cbind( c(1,3), 1, c(0, .3 ) )
mod2_3 <- TAM::tam.mml(resp=resp, Y=Y, beta.fixed=beta.fixed )

# (5) same model but with R regression formula for Y
dataY <- data.frame(Y)
colnames(dataY) <- c("Y1","Y2")
mod2_4 <- TAM::tam.mml(resp=resp, dataY=dataY, formulaY=~ Y1+Y2 )
summary(mod2_4)

# (6) model with interaction of regressors
mod2_5 <- TAM::tam.mml(resp=resp, dataY=dataY, formulaY=~ Y1*Y2 )
summary(mod2_5)

# (7) no constraint on regressors (removing constraint from intercept)
mod2_6 <- TAM::tam.mml(resp=resp, Y=Y, beta.fixed=FALSE )

#############################################################################
# EXAMPLE 3: Multiple group estimation
#############################################################################

# (1) simulate data
set.seed(6778)
N <- 3000
theta <- c( stats::rnorm(N/2,mean=0,sd=1.5), stats::rnorm(N/2,mean=.5,sd=1)  )
I <- 20
p1 <- stats::plogis( outer( theta, seq( -2, 2, len=I ), "-" ) )
resp <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
colnames(resp) <- paste("I", 1:I, sep="")
group <- rep(1:2, each=N/2 )

# (2) estimate model
mod3_1 <- TAM::tam.mml( resp,  group=group )
summary(mod3_1)

#############################################################################
# EXAMPLE 4: Multidimensional estimation
# with two dimensional theta's - simulate some bivariate data,
# and regressors
# 40 items: first 20 items load on dimension 1,
#           second 20 items load on dimension 2
#############################################################################

# (1) simulate some data
set.seed(6778)
library(mvtnorm)
N <- 1000
Y <- cbind( stats::rnorm( N ), stats::rnorm(N) )
theta <- mvtnorm::rmvnorm( N,mean=c(0,0), sigma=matrix( c(1,.5,.5,1), 2, 2 ))
theta[,1] <- theta[,1] + .4 * Y[,1] + .2 * Y[,2]  # latent regression model
theta[,2] <- theta[,2] + .8 * Y[,1] + .5 * Y[,2]  # latent regression model
I <- 20
p1 <- stats::plogis( outer( theta[,1], seq( -2, 2, len=I ), "-" ) )
resp1 <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
p1 <- stats::plogis( outer( theta[,2], seq( -2, 2, len=I ), "-" ) )
resp2 <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
resp <- cbind(resp1,resp2)
colnames(resp) <- paste("I", 1:(2*I), sep="")

# (2) define loading Matrix
Q <- array( 0, dim=c( 2*I, 2 ))
Q[cbind(1:(2*I), c( rep(1,I), rep(2,I) ))] <- 1

# (3) estimate models

#************************************************************
# Model 4.1: Rasch model: without regressors
mod4_1 <- TAM::tam.mml( resp=resp, Q=Q )

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    F1=~ 1*I1__I20
    F2=~ 1*I21__I40
    # Alternatively to the factor 1 one can use the item type Rasch
    F1 ~~ F1
    F2 ~~ F2
    F1 ~~ F2
    "
mod4_1t <- TAM::tamaan( tammodel, resp, control=list(maxiter=100))
summary(mod4_1t)

#************************************************************
# Model 4.1b: estimate model with sum constraint of items for each dimension
mod4_1b <- TAM::tam.mml( resp=resp, Q=Q,constraint="items")

#************************************************************
# Model 4.2: Rasch model: set covariance between dimensions to zero
variance_fixed <- cbind( 1, 2, 0 )
mod4_2 <- TAM::tam.mml( resp=resp, Q=Q, variance.fixed=variance_fixed )
summary(mod4_2)

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    F1=~ I1__I20
    F2=~ I21__I40
    F1 ~~ F1
    F2 ~~ F2
    F1 ~~ 0*F2
  ITEM TYPE:
    ALL(Rasch)
    "
mod4_2t <- TAM::tamaan( tammodel, resp)
summary(mod4_2t)

#************************************************************
# Model 4.3: 2PL model
mod4_3 <- TAM::tam.mml.2pl( resp=resp, Q=Q, irtmodel="2PL" )

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    F1=~ I1__I20
    F2=~ I21__I40
    F1 ~~ F1
    F2 ~~ F2
    F1 ~~ F2
    "
mod4_3t <- TAM::tamaan( tammodel, resp )
summary(mod4_3t)

#************************************************************
# Model 4.4: Rasch model with 2000 quasi monte carlo nodes
# -> nodes are useful for more than 3 or 4 dimensions
mod4_4 <- TAM::tam.mml( resp=resp, Q=Q, control=list(snodes=2000) )

#************************************************************
# Model 4.5: Rasch model with 2000 stochastic nodes
mod4_5 <- TAM::tam.mml( resp=resp, Q=Q,control=list(snodes=2000,QMC=FALSE))

#************************************************************
# Model 4.6: estimate two dimensional Rasch model with regressors
mod4_6 <- TAM::tam.mml( resp=resp, Y=Y, Q=Q )

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    F1=~ I1__I20
    F2=~ I21__I40
    F1 ~~ F1
    F2 ~~ F2
    F1 ~~ F2
  ITEM TYPE:
    ALL(Rasch)
    "
mod4_6t <- TAM::tamaan( tammodel, resp, Y=Y )
summary(mod4_6t)

#############################################################################
# EXAMPLE 5: 2-dimensional estimation with within item dimensionality
#############################################################################
library(mvtnorm)
# (1) simulate data
set.seed(4762)
N <- 2000 # 2000 persons
Y <- stats::rnorm( N )
theta <- mvtnorm::rmvnorm( N,mean=c(0,0), sigma=matrix( c(1,.5,.5,1), 2, 2 ))
I <- 10
# 10 items load on the first dimension
p1 <- stats::plogis( outer( theta[,1], seq( -2, 2, len=I ), "-" ) )
resp1 <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
# 10 items load on the second dimension
p1 <- stats::plogis( outer( theta[,2], seq( -2, 2, len=I ), "-" ) )
resp2 <- 1 * ( p1 > matrix( stats::runif( N*I ), nrow=N, ncol=I ) )
# 20 items load on both dimensions
p1 <- stats::plogis( outer( 0.5*theta[,1] + 1.5*theta[,2], seq(-2,2,len=2*I ), "-" ))
resp3 <- 1 * ( p1 > matrix( stats::runif( N*2*I ), nrow=N, ncol=2*I ) )
#Combine the two sets of items into one response matrix
resp <- cbind(resp1, resp2, resp3 )
colnames(resp) <- paste("I", 1:(4*I), sep="")

# (2) define loading matrix
Q <- cbind(c(rep(1,10),rep(0,10),rep(1,20)), c(rep(0,10),rep(1,10),rep(1,20)))

# (3) model: within item dimensionality and 2PL estimation
mod5 <- TAM::tam.mml.2pl(resp, Q=Q, irtmodel="2PL" )
summary(mod5)

# item difficulties
mod5$item
# item loadings
mod5$B

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    F1=~ I1__I10 + I21__I40
    F2=~ I11__I20 + I21__I40
    F1 ~~ 1*F1
    F1 ~~ F2
    F2 ~~ 1*F2
    "
mod5t <- TAM::tamaan( tammodel, resp,  control=list(maxiter=10))
summary(mod5t)

#############################################################################
# EXAMPLE 6: ordered data - Generalized partial credit model
#############################################################################
data(data.gpcm, package="TAM")

#************************************************************
# Ex6.1: Nominal response model (irtmodel="2PL")
mod6_1 <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="2PL", control=list(maxiter=200) )
mod6_1$item # item intercepts
mod6_1$B    # for every category a separate slope parameter is estimated

# reestimate the model with fixed item parameters
mod6_1a <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="2PL",
       xsi.fixed=mod6_1$xsi.fixed.estimated,  B.fixed=mod6_1$B.fixed.estimated,
       est.variance=TRUE )

# estimate the model with initial item parameters from mod6_1
mod6_1b <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="2PL",
       xsi.inits=mod6_1$xsi.fixed.estimated,  B=mod6_1$B )

#************************************************************
# Ex6.2: Generalized partial credit model
mod6_2 <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="GPCM", control=list(maxiter=200))
mod6_2$B[,2,]    # joint slope parameter for all categories

#************************************************************
# Ex6.3: some fixed entries of slope matrix B
# B: nitems x maxK x ndim
#   ( number of items x maximum number of categories x number of dimensions)
# set two constraints
B.fixed <- matrix( 0, 2, 4 )
# set second item, score of 2 (category 3), at first dimension to 2.3
B.fixed[1,] <- c(2,3,1,2.3)
# set third item, score of 1 (category 2), at first dimension to 1.4
B.fixed[2,] <- c(3,2,1,1.4)

# estimate item parameter with variance fixed (by default)
mod6_3 <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="2PL", B.fixed=B.fixed,
                 control=list( maxiter=200) )
mod6_3$B

#************************************************************
# Ex 6.4: estimate the same model, but estimate variance
mod6_4 <- TAM::tam.mml.2pl( resp=data.gpcm, irtmodel="2PL", B.fixed=B.fixed,
               est.variance=TRUE, control=list( maxiter=350) )
mod6_4$B

#************************************************************
# Ex 6.5: partial credit model
mod6_5 <- TAM::tam.mml( resp=data.gpcm,control=list( maxiter=200) )
mod6_5$B

#************************************************************
# Ex 6.6: partial credit model: Conquest parametrization 'item+item*step'
mod6_6 <- TAM::tam.mml( resp=data.gpcm, irtmodel="PCM2" )
summary(mod6_6)

# estimate mod6_6 applying the sum constraint of item difficulties
# modify design matrix of xsi paramters
A1 <- TAM::.A.PCM2(resp=data.gpcm )
A1[3,2:4,"Comfort"] <- 1:3
A1[3,2:4,"Work"] <- 1:3
A1 <- A1[,, -3] # remove Benefit xsi item parameter
# estimate model
mod6_6b <- TAM::tam.mml( resp=data.gpcm, A=A1, beta.fixed=FALSE )
summary(mod6_6b)

# estimate model with argument constraint="items"
mod6_6c <- TAM::tam.mml( resp=data.gpcm, irtmodel="PCM2", constraint="items")

# estimate mod6_6 using tam.mml.mfr
mod6_6d <- TAM::tam.mml.mfr( resp=data.gpcm, formulaA=~ 0 + item + item:step,
    control=list(fac.oldxsi=.1), constraint="items" )
summary(mod6_6d)

#************************************************************
# Ex 6.7: Rating scale model: Conquest parametrization 'item+step'
mod6_7 <- TAM::tam.mml( resp=data.gpcm, irtmodel="RSM" )
summary(mod6_7)

#************************************************************
# Ex 6.8: sum constraint on item difficulties
#         partial credit model: ConQuest parametrization 'item+item*step'
#         polytomous scored TIMMS data
#         compare to Example 16
#

data(data.timssAusTwn.scored)
dat <- data.timssAusTwn.scored[,1:11]

## > tail(sort(names(dat)),1) # constrained item
## [1] "M032761"

# modify design matrix of xsi paramters
A1 <- TAM::.A.PCM2( resp=dat )
# constrained item loads on every other main item parameter
# with opposing margin it had been loaded on its own main item parameter
A1["M032761",,setdiff(colnames(dat), "M032761")] <- -A1["M032761",,"M032761"]
# remove main item parameter for constrained item
A1 <- A1[,, setdiff(dimnames(A1)[[3]],"M032761")]

# estimate model
mod6_8a <- TAM::tam.mml( resp=dat, A=A1, beta.fixed=FALSE )
summary(mod6_8a)
# extract fixed item parameter for item M032761
## - sum(mod6_8a$xsi[setdiff(colnames(dat), "M032761"),"xsi"])

# estimate mod6_8a using tam.mml.mfr
## fixed a bug in 'tam.mml.mfr' for differing number of categories
## per item -> now a xsi vector with parameter fixings to values
## of 99 is used
mod6_8b <- TAM::tam.mml.mfr( resp=dat, formulaA=~ 0 + item + item:step,
                        control=list(fac.oldxsi=.1), constraint="items" )
summary(mod6_8b)

#************************************************************
# Ex 6.9: sum constraint on item difficulties for irtmodel="PCM"

data(data.timssAusTwn.scored)
dat <- data.timssAusTwn.scored[,2:11]
dat[ dat==9 ] <- NA

# obtain the design matrix for the PCM parametrization and
# the number of categories for each item
maxKi <- apply(dat, 2, max, na.rm=TRUE)
des <- TAM::designMatrices(resp=dat)
A1 <- des$A

# define the constrained item category and remove the respective parameter
(par <- unlist( strsplit(dimnames(A1)[[3]][dim(A1)[3]], split="_") ))
A1 <- A1[,,-dim(A1)[3]]

# the item category loads on every other item category parameter with
# opposing margin, balancing the number of categories for each item
item.id <- which(colnames(dat)==par[1])
cat.id <- maxKi[par[1]]+1
loading <- 1/rep(maxKi, maxKi)
loading <- loading [-which(names(loading)==par[1])[1]]
A1[item.id, cat.id, ] <- loading
A1[item.id,,]

# estimate model
mod6_9 <- TAM::tam.mml( resp=dat, A=A1, beta.fixed=FALSE )
summary(mod6_9)

## extract fixed item category parameter
# calculate mean for each item
ind.item.cat.pars <- sapply(colnames(dat), grep, rownames(mod6_8$xsi))
item.means <- lapply(ind.item.cat.pars, function(ii) mean(mod6_8$xsi$xsi[ii]))

# these sum up to the negative of the fixed parameter
fix.par <- -sum( unlist(item.means), na.rm=TRUE)

#************************************************************
# Ex 6.10: Generalized partial credit model with equality constraints
#          on item discriminations

data(data.gpcm)
dat <- data.gpcm

# Ex 6.10a: set all slopes of three items equal to each other
E <- matrix( 1, nrow=3, ncol=1 )
mod6_10a <- TAM::tam.mml.2pl( dat, irtmodel="GPCM.design", E=E  )
summary(mod6_10a)
mod6_10a$B[,,]

# Ex 6.10b: equal slope for first and third item
E <- matrix( 0, nrow=3, ncol=2 )
E[c(1,3),1] <- 1
E[ 2, 2 ] <- 1
mod6_10b <- TAM::tam.mml.2pl( dat, irtmodel="GPCM.design", E=E  )
summary(mod6_10b)
mod6_10b$B[,,]

#############################################################################
# EXAMPLE 7: design matrix for slopes for the generalized partial credit model
#############################################################################

# (1) simulate data from a model with a (item slope) design matrix E
set.seed(789)
I <- 42
b <- seq( -2, 2, len=I)
# create design matrix for loadings
E <- matrix( 0, I, 5 )
E[ seq(1,I,3), 1 ] <- 1
E[ seq(2,I,3), 2 ] <- 1
E[ seq(3,I,3), 3 ] <- 1
ind <- seq( 1, I, 2 ) ; E[ ind, 4 ] <- rep( c( .3, -.2 ), I )[ 1:length(ind) ]
ind <- seq( 2, I, 4 ) ; E[ ind, 5 ] <- rep( .15, I )[ 1:length(ind) ]
E

# true basis slope parameters
lambda <- c( 1, 1.2, 0.8, 1, 1.1 )
# calculate item slopes
a <- E %*% lambda

# simulate
N <- 4000
theta <- stats::rnorm( N )
aM <- outer( rep(1,N), a[,1] )
bM <- outer( rep(1,N), b )
pM <- stats::plogis( aM * ( matrix( theta, nrow=N, ncol=I  ) - bM ) )
dat <- 1 * ( pM > stats::runif( N*I ) )
colnames(dat) <- paste("I", 1:I, sep="")

# estimate model
mod7 <- TAM::tam.mml.2pl( resp=dat, irtmodel="GPCM.design", E=E )
mod7$B

# recalculate estimated basis parameters
stats::lm( mod7$B[,2,1] ~ 0+ as.matrix(E ) )
  ##   Call:
  ##   lm(formula=mod7$B[, 2, 1] ~ 0 + as.matrix(E))
  ##   Coefficients:
  ##   as.matrix(E)1  as.matrix(E)2  as.matrix(E)3  as.matrix(E)4  as.matrix(E)5
  ##          0.9904         1.1896         0.7817         0.9601         1.2132

#############################################################################
# EXAMPLE 8: Differential item functioning                                  #
#  A first example of a Multifaceted Rasch Model                            #
#  Facet is only female; 10 items are studied                               #
#############################################################################
data(data.ex08)

formulaA <- ~ item+female+item*female
# this formula is in R equivalent to 'item*female'
resp <- data.ex08[["resp"]]
facets <- as.data.frame( data.ex08[["facets"]] )

#***
# Model 8a: investigate gender DIF on all items
mod8a <- TAM::tam.mml.mfr( resp=resp, facets=facets, formulaA=formulaA )
summary(mod8a)

#***
# Model 8a 2: specification with long format response data
resp.long <- c( data.ex08[["resp"]] )
pid <- rep( 1:nrow(data.ex08[["resp"]]), ncol(data.ex08[["resp"]]) )

itemnames <- rep(colnames(data.ex08[["resp"]]), each=nrow(data.ex08[["resp"]]))
facets.long <- cbind( data.frame( "item"=itemnames ),
                 data.ex08[["facets"]][pid,,drop=F] )

mod8a_2 <- TAM::tam.mml.mfr( resp=resp.long, facets=facets.long,
                      formulaA=formulaA, pid=pid)
stopifnot( all(mod8a$xsi.facets$xsi==mod8a_2$xsi.facets$xsi) )

#***
# Model 8b: Differential bundle functioning (DBF)
#   - investigate differential item functioning in item groups

# modify pre-specified design matrix to define 'appropriate' DBF effects
formulaA <- ~ item*female
des <- TAM::designMatrices.mfr( resp=resp, facets=facets, formulaA=formulaA)
A1 <- des$A$A.3d
# item group A: items 1-5
# item group B: items 6-8
# item group C: items 9-10
A1 <- A1[,,1:13]
dimnames(A1)[[3]][ c(12,13) ] <- c("A:female1", "B:female1")
# item group A
A1[,2,12] <- 0
A1[c(1,5,7,9,11),2,12] <- -1
A1[c(1,5,7,9,11)+1,2,12] <- 1
# item group B
A1[,2,13] <- 0
A1[c(13,15,17),2,13] <- -1
A1[c(13,15,17)+1,2,13] <- 1
# item group C (define effect(A)+effect(B)+effect(C)=0)
A1[c(19,3),2,c(12,13)] <- 1
A1[c(19,3)+1,2,c(12,13)] <- -1
#   A1[,2,]   # look at modified design matrix
# estimate model
mod8b <- TAM::tam.mml( resp=des$gresp$gresp.noStep, A=A1 )
summary(mod8b)

#############################################################################
# EXAMPLE 9: Multifaceted Rasch Model
#############################################################################

data(data.sim.mfr)
data(data.sim.facets)

# two way interaction item and rater
formulaA <- ~item+item:step + item*rater
mod9a <- TAM::tam.mml.mfr( resp=data.sim.mfr, facets=data.sim.facets, formulaA=formulaA)
mod9a$xsi.facets
summary(mod9a)

# three way interaction item, female and rater
formulaA <- ~item+item:step + female*rater + female*item*step
mod9b <- TAM::tam.mml.mfr( resp=data.sim.mfr, facets=data.sim.facets, formulaA=formulaA)
summary(mod9b)

#############################################################################
# EXAMPLE 10: Model with raters.
#   Persons are arranged in multiple rows which is indicated
#   by multiple person identifiers.
#############################################################################
data(data.ex10)
dat <- data.ex10
head(dat)
  ##     pid rater I0001 I0002 I0003 I0004 I0005
  ## 1     1     1     0     1     1     0     0
  ## 451   1     2     1     1     1     1     0
  ## 901   1     3     1     1     1     0     1
  ## 452   2     2     1     1     1     0     1
  ## 902   2     3     1     1     0     1     1

facets <- dat[, "rater", drop=FALSE ] # define facet (rater)
pid <- dat$pid      # define person identifier (a person occurs multiple times)
resp <- dat[, -c(1:2) ]        # item response data
formulaA <- ~ item * rater      # formula

mod10 <- TAM::tam.mml.mfr( resp=resp, facets=facets, formulaA=formulaA, pid=dat$pid )
summary(mod10)

# estimate person parameter with WLE
wmod10 <- TAM::tam.wle( mod10 )

#--- Example 10a
# compare model containing only item
formulaA <- ~ item + rater      # pseudo formula for item
xsi.setnull <- "rater"          # set all rater effects to zero
mod10a <- TAM::tam.mml.mfr( resp=resp, facets=facets, formulaA=formulaA,
            xsi.setnull=xsi.setnull, pid=dat$pid, beta.fixed=cbind(1,1,0))
summary(mod10a)

# A shorter way for specifying this example is
formulaA <- ~ item + 0*rater        # set all rater effects to zero
mod10a1 <- TAM::tam.mml.mfr( resp=resp, facets=facets, formulaA=formulaA, pid=dat$pid )
summary(mod10a1)

# tam.mml.mfr also appropriately extends the facets data frame with pseudo facets
# if necessary
formulaA <- ~ item     # omitting the rater term
mod10a2 <- TAM::tam.mml.mfr( resp=resp, facets=facets, formulaA=formulaA, pid=dat$pid )
  ##   Item Parameters Xsi
  ##              xsi se.xsi
  ##   I0001   -1.931  0.111
  ##   I0002   -1.023  0.095
  ##   I0003   -0.089  0.089
  ##   I0004    1.015  0.094
  ##   I0005    1.918  0.110
  ##   psfPF11  0.000  0.000
  ##   psfPF12  0.000  0.000

#***
# Model 10_2: specification with long format response data
resp.long <- c(unlist( dat[, -c(1:2) ] ))

pid <- rep( dat$pid, ncol(dat[, -c(1:2) ]) )
itemnames <- rep(colnames(dat[, -c(1:2) ]), each=nrow(dat[, -c(1:2) ]))

# quick note: the 'trick' to use pid as the row index of the facet  (cf., used in Ex 8a_2)
# is not working here, since pid already occures multiple times in the original response data
facets <- cbind( data.frame("item"=itemnames),
                 dat[ rep(1:nrow(dat), ncol(dat[,-c(1:2)])), "rater",drop=F]
)

mod10_2 <- TAM::tam.mml.mfr( resp=resp.long, facets=facets, formulaA=formulaA, pid=pid)

stopifnot( all(mod10$xsi.facets$xsi==mod10_2$xsi.facets$xsi) )

#############################################################################
# EXAMPLE 11: Dichotomous data with missing and omitted responses
#############################################################################
data(data.ex11) ; dat <- data.ex11

#***
# Model 11a: Calibration (item parameter estimating) in which omitted
#            responses (code 9) are set to missing
dat1 <- dat[,-1]
dat1[ dat1==9 ] <- NA
# estimate Rasch model
mod11a <- TAM::tam.mml( resp=dat1 )
summary(mod11a)
# compute person parameters
wmod11a <- TAM::tam.wle( mod11a )

#***
# Model 11b: Scaling persons (WLE estimation) setting omitted
#            responses as incorrect and using fixed
#            item parameters form Model 11a

# set matrix with fixed item difficulties as the input
xsi1 <- mod11a$xsi    # xsi output from Model 11a
xsi.fixed <- cbind( seq(1,nrow(xsi1) ), xsi1$xsi )
# recode 9 to 0
dat2 <- dat[,-1]
dat2[ dat2==9 ] <- 0
# run Rasch model with fixed item difficulties
mod11b <- TAM::tam.mml( resp=dat2, xsi.fixed=xsi.fixed )
summary(mod11b)
# WLE estimation
wmod11b <- TAM::tam.wle( mod11b )

#############################################################################
# EXAMPLE 12: Avoiding nonconvergence using the argument increment.factor
#############################################################################
data(data.ex12)
dat <- data.ex12

# non-convergence without increment.factor
mod1 <- TAM::tam.mml.2pl( resp=data.ex12, control=list( maxiter=1000) )

# avoiding non-convergence with increment.factor=1.02
mod2 <- TAM::tam.mml.2pl( resp=data.ex12,
            control=list( maxiter=1000, increment.factor=1.02) )
summary(mod1)
summary(mod2)

#############################################################################
# EXAMPLE 13: Longitudinal data 'data.long' from sirt package
#############################################################################
library(sirt)
data(data.long, package="sirt")
dat <- data.long
  ##   > colnames(dat)
  ##    [1] "idstud" "I1T1"   "I2T1"   "I3T1"   "I4T1"   "I5T1"   "I6T1"
  ##    [8] "I3T2"   "I4T2"   "I5T2"   "I6T2"   "I7T2"   "I8T2"

## item 1 to 6 administered at T1 and items 3 to 8 at T2
## items 3 to 6 are anchor items

#***
# Model 13a: 2-dimensional Rasch model assuming invariant item difficulties

# define matrix loadings
Q <- matrix(0,12,2)
colnames(Q) <- c("T1","T2")
Q[1:6,1] <- 1       # items at T1
Q[7:12,2] <- 1      # items at T2

# assume equal item difficulty of I3T1 and I3T2, I4T1 and I4T2, ...
# create draft design matrix and modify it
A <- TAM::designMatrices(resp=data.long[,-1])$A
dimnames(A)[[1]] <- colnames(data.long)[-1]
  ##   > str(A)
  ##    num [1:12, 1:2, 1:12] 0 0 0 0 0 0 0 0 0 0 ...
  ##    - attr(*, "dimnames")=List of 3
  ##     ..$ : chr [1:12] "Item01" "Item02" "Item03" "Item04" ...
  ##     ..$ : chr [1:2] "Category0" "Category1"
  ##     ..$ : chr [1:12] "I1T1" "I2T1" "I3T1" "I4T1" ...
A1 <- A[,, c(1:6, 11:12 ) ]
dimnames(A1)[[3]] <- substring( dimnames(A1)[[3]],1,2)
A1[7,2,3] <- -1     # difficulty(I3T1)=difficulty(I3T2)
A1[8,2,4] <- -1     # I4T1=I4T2
A1[9,2,5] <- A1[10,2,6] <- -1
  ##   > A1[,2,]
  ##        I1 I2 I3 I4 I5 I6 I7 I8
  ##   I1T1 -1  0  0  0  0  0  0  0
  ##   I2T1  0 -1  0  0  0  0  0  0
  ##   I3T1  0  0 -1  0  0  0  0  0
  ##   I4T1  0  0  0 -1  0  0  0  0
  ##   I5T1  0  0  0  0 -1  0  0  0
  ##   I6T1  0  0  0  0  0 -1  0  0
  ##   I3T2  0  0 -1  0  0  0  0  0
  ##   I4T2  0  0  0 -1  0  0  0  0
  ##   I5T2  0  0  0  0 -1  0  0  0
  ##   I6T2  0  0  0  0  0 -1  0  0
  ##   I7T2  0  0  0  0  0  0 -1  0
  ##   I8T2  0  0  0  0  0  0  0 -1

# estimate model
# set intercept of second dimension (T2) to zero
beta.fixed <- cbind( 1, 2, 0 )
mod13a <- TAM::tam.mml( resp=data.long[,-1], Q=Q, A=A1, beta.fixed=beta.fixed)
summary(mod13a)

#--- tamaan specification
tammodel <- "
  LAVAAN MODEL:
    T1=~ 1*I1T1__I6T1
    T2=~ 1*I3T2__I8T2
    T1 ~~ T1
    T2 ~~ T2
    T1 ~~ T2
    # constraint on item difficulties
    I3T1 + I3T2 | b3*t1
    I4T1 + I4T2 | b4*t1
    I5T1 + I5T2 | b5*t1
    I6T1 + I6T2 | b6*t1
    "
# The constraint on item difficulties can be more efficiently written as
  ##       DO(3,6,1)
  ##         I%T1 + I%T2 | b%*t1
  ##       DOEND
# estimate model
mod13at <- TAM::tamaan( tammodel, resp=data.long,  beta.fixed=beta.fixed )
summary(mod13at)

#***
# Model 13b: invariant item difficulties with zero mean item difficulty
#           of anchor items

A <- TAM::designMatrices(resp=data.long[,-1])$A
dimnames(A)[[1]] <- colnames(data.long)[-1]
A1 <- A[,, c(1:5, 11:12 ) ]
dimnames(A1)[[3]] <- substring( dimnames(A1)[[3]],1,2)
A1[7,2,3] <- -1     # difficulty(I3T1)=difficulty(I3T2)
A1[8,2,4] <- -1     # I4T1=I4T2
A1[9,2,5] <- -1
A1[6,2,3] <- A1[6,2,4] <- A1[6,2,5] <- 1     # I6T1=-(I3T1+I4T1+I5T1)
A1[10,2,3] <- A1[10,2,4] <- A1[10,2,5] <- 1  # I6T2=-(I3T2+I4T2+I5T2)
A1[,2,]
  ##      I1 I2 I3 I4 I5 I7 I8
  ## I1T1 -1  0  0  0  0  0  0
  ## I2T1  0 -1  0  0  0  0  0
  ## I3T1  0  0 -1  0  0  0  0
  ## I4T1  0  0  0 -1  0  0  0
  ## I5T1  0  0  0  0 -1  0  0
  ## I6T1  0  0  1  1  1  0  0
  ## I3T2  0  0 -1  0  0  0  0
  ## I4T2  0  0  0 -1  0  0  0
  ## I5T2  0  0  0  0 -1  0  0
  ## I6T2  0  0  1  1  1  0  0
  ## I7T2  0  0  0  0  0 -1  0
  ## I8T2  0  0  0  0  0  0 -1

mod13b <- TAM::tam.mml( resp=data.long[,-1], Q=Q, A=A1, beta.fixed=FALSE)
summary(mod13b)

#***
# Model 13c: longitudinal polytomous data
#

# modifiy Items I1T1, I4T1, I4T2 in order to be trichotomous (codes: 0,1,2)
set.seed(42)
dat <- data.long
dat[(1:50),2] <- sample(c(0,1,2), 50, replace=TRUE)
dat[(1:50),5] <- sample(c(0,1,2), 50, replace=TRUE)
dat[(1:50),9] <- sample(c(0,1,2), 50, replace=TRUE)
  ##   > colnames(dat)
  ##    [1] "idstud" "I1T1"   "I2T1"   "I3T1"   "I4T1"   "I5T1"   "I6T1"
  ##    [8] "I3T2"   "I4T2"   "I5T2"   "I6T2"   "I7T2"   "I8T2"

## item 1 to 6 administered at T1, items 3 to 8 at T2
## items 3 to 6 are anchor items

# (1) define matrix loadings
Q <- matrix(0,12,2)
colnames(Q) <- c("T1","T2")
Q[1:6,1] <- 1       # items at T1
Q[7:12,2] <- 1      # items at T2

# (2) assume equal item difficulty of anchor items
#     create draft design matrix and modify it
A <- TAM::designMatrices(resp=dat[,-1])$A
dimnames(A)[[1]] <- colnames(dat)[-1]
  ## > str(A)
  ## num [1:12, 1:3, 1:15] 0 0 0 0 0 0 0 0 0 0 ...
  ## - attr(*, "dimnames")=List of 3
  ## ..$ : chr [1:12] "I1T1" "I2T1" "I3T1" "I4T1" ...
  ## ..$ : chr [1:3] "Category0" "Category1" "Category2"
  ## ..$ : chr [1:15] "I1T1_Cat1" "I1T1_Cat2" "I2T1_Cat1" "I3T1_Cat1" ...

# define matrix A
# Items 1 to 3 administered at T1, Items 3 to 6 are anchor items
# Item 7 to 8 administered at T2
# Item I1T1, I4T1, I4T2 are trichotomous (codes: 0,1,2)
A1 <- A[,, c(1:8, 14:15) ]
dimnames(A1)[[3]] <- gsub("T1|T2", "",  dimnames(A1)[[3]])

# Modifications are shortened compared to Model 13 a, but are still valid
A1[7,,] <- A1[3,,]  # item 7, i.e. I3T2, loads on same parameters as
                    # item 3, I3T1
A1[8,,] <- A1[4,,]  # same for item 8 and item 4
A1[9,,] <- A1[5,,]  # same for item 9 and item 5
A1[10,,] <- A1[6,,] # same for item 10 and item 6
  ## > A1[8,,]
  ##           I1_Cat1 I1_Cat2 I2_Cat1 I3_Cat1 I4_Cat1 I4_Cat2 I5_Cat1 ...
  ## Category0       0       0       0       0       0       0       0
  ## Category1       0       0       0       0      -1       0       0
  ## Category2       0       0       0       0      -1      -1       0

# (3) estimate model
#     set intercept of second dimension (T2) to zero
beta.fixed <- cbind( 1, 2, 0 )
mod13c <- TAM::tam.mml( resp=dat[,-1], Q=Q, A=A1, beta.fixed=beta.fixed,
                   irtmodel="PCM")
summary(mod13c)
wle.mod13c <- TAM::tam.wle(mod13c) # WLEs of dimension T1 and T2

#############################################################################
# EXAMPLE 14: Facet model with latent regression
#############################################################################
data( data.ex14 )
dat <- data.ex14

#***
# Model 14a: facet model
resp <- dat[, paste0("crit",1:7,sep="") ]    # item data
facets <- data.frame( "rater"=dat$rater )     # define facets
formulaA <- ~item+item*step + rater
mod14a <- TAM::tam.mml.mfr( resp, facets=facets, formulaA=formulaA, pid=dat$pid )
summary(mod14a)

#***
# Model 14b: facet model with latent regression
#   Note that dataY must correspond to rows in resp and facets which means
#   that there must be the same rows in Y for a person with multiple rows
#   in resp
dataY <- dat[, c("X1","X2") ]        # latent regressors
formulaY <- ~ X1+X2            # latent regression formula
mod14b <- TAM::tam.mml.mfr( resp, facets=facets, formulaA=formulaA,
            dataY=dataY, formulaY=formulaY, pid=dat$pid)
summary(mod14b)

#***
# Model 14c: Multi-facet model with item slope estimation
# use design matrix and modified response data from Model 1
# item-specific slopes

resp1 <- mod14a$resp      # extract response data with generalized items
A <- mod14a$A             # extract design matrix for item intercepts
colnames(resp1)

# define design matrix for slopes
E <- matrix( 0, nrow=ncol(resp1), ncol=7 )
colnames(E) <- paste0("crit",1:7)
rownames(E) <- colnames(resp1)
E[ cbind( 1:(7*7), rep(1:7,each=7) ) ] <- 1

mod14c <- TAM::tam.mml.2pl( resp=resp1, A=A, irtmodel="GPCM.design", E=E,
        control=list(maxiter=100) )
summary(mod14c)

#############################################################################
# EXAMPLE 15: Coping with nonconvergent models
#############################################################################

data(data.ex15)
data <- data.ex15
# facet model 'group*item' is of interest

#***
# Model 15a:
mod15a <- TAM::tam.mml.mfr(resp=data[,-c(1:2)],facets=data[,"group",drop=FALSE],
    formulaA=~ item + group*item, pid=data$pid )
# See output:
  ##
  ##   Iteration 47     2013-09-10 16:51:39
  ##   E Step
  ##   M Step Intercepts   |----
  ##     Deviance=75510.2868 | Deviance change: -595.0609
  ##   !!! Deviance increases!                                        !!!!
  ##   !!! Choose maybe fac.oldxsi > 0 and/or increment.factor > 1    !!!!
  ##     Maximum intercept parameter change: 0.925045
  ##     Maximum regression parameter change: 0
  ##     Variance:  0.9796  | Maximum change: 0.009226

#***
# Model 15b: Follow the suggestions of changing the default of fac.oldxsi and
#            increment.factor
mod15b <- TAM::tam.mml.mfr(resp=data[,-c(1:2)],facets=data[,"group",drop=FALSE],
            formulaA=~ group*item, pid=data$pid,
            control=list( increment.factor=1.03, fac.oldxsi=.4 ) )

#***
# Model 15c: Alternatively, just choose more iterations in M-step by "Msteps=10"
mod15c <- TAM::tam.mml.mfr(resp=data[,-c(1:2)],facets=data[,"group",drop=FALSE],
    formulaA=~ item + group*item, pid=data$pid,
    control=list(maxiter=250, Msteps=10))

#############################################################################
# EXAMPLE 16: Differential item function for polytomous items and
#             differing number of response options per item
#############################################################################

data(data.timssAusTwn.scored)
dat <- data.timssAusTwn.scored
# extract item response data
resp <- dat[, sort(grep("M", colnames(data.timssAusTwn.scored), value=TRUE)) ]
# some descriptives
psych::describe(resp)
# define facets: 'cnt' is group identifier
facets <- data.frame( "cnt"=dat$IDCNTRY)
# create design matrices
des2 <- TAM::designMatrices.mfr2( resp=resp, facets=facets,
                   formulaA=~item*step + item*cnt)
# restructured data set: pseudoitem=item x country
resp2 <- des2$gresp$gresp.noStep
# A design matrix
A <- des2$A$A.3d
    # redundant xsi parameters must be eliminated from design matrix
xsi.elim <- des2$xsi.elim
A <- A[,, - xsi.elim[,2] ]
# extract loading matrix B
B <- des2$B$B.3d
# estimate model
mod1 <- TAM::tam.mml( resp=resp2, A=A, B=B, control=list(maxiter=100) )
summary(mod1)
# The sum of all DIF parameters is set to zero. The DIF parameter for the last
# item is therefore obtained
xsi1 <- mod1$xsi
difxsi <- xsi1[ intersect( grep("cnt",rownames(xsi1)),
              grep("M03",rownames(xsi1))), ]   - colSums(difxsi)
    # this is the DIF effect of the remaining item

#############################################################################
# EXAMPLE 17: Several multidimensional and subdimension models
#############################################################################

library(mirt)
#*** (1) simulate data in mirt package
set.seed(9897)
# simulate data according to the four-dimensional Rasch model
# variances
variances <- c( 1.45, 1.74, .86, 1.48  )
# correlations
corrs <- matrix( 1, 4, 4 )
dd1 <- 1 ; dd2 <- 2 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .88
dd1 <- 1 ; dd2 <- 3 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .85
dd1 <- 1 ; dd2 <- 4 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .87
dd1 <- 2 ; dd2 <- 3 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .84
dd1 <- 2 ; dd2 <- 4 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .90
dd1 <- 3 ; dd2 <- 4 ; corrs[dd1,dd2] <- corrs[dd2,dd1] <- .90
# covariance matrix
covar <- outer( sqrt( variances), sqrt(variances) )*corrs
# item thresholds and item discriminations
d <- matrix( stats::runif(40, -2, 2 ), ncol=1 )
a <- matrix(NA, nrow=40,ncol=4)
a[1:10,1] <- a[11:20,2] <- a[21:30,3] <- a[31:40,4] <- 1
# simulate data
dat <- mirt::simdata(a=a, d=d, N=1000, itemtype="dich", sigma=covar )
# define Q-matrix for testlet and subdimension models estimated below
Q <- matrix( 0, nrow=40, ncol=5 )
colnames(Q) <- c("g", paste0( "subd", 1:4) )
Q[,1] <- 1
Q[1:10,2] <- Q[11:20,3] <- Q[21:30,4] <- Q[31:40,5] <- 1

# define maximum number of iterations and number of quasi monte carlo nodes
# maxit <- 5  ; snodes <- 300    # this specification is only for speed reasons
maxit <- 200 ; snodes <- 1500

#*****************
# Model 1: Rasch testlet model
#*****************

# define a user function for restricting the variance according to the
# Rasch testlet model
variance.fct1 <- function( variance ){
            ndim <- ncol(variance)
            variance.new <- matrix( 0, ndim, ndim )
            diag(variance.new) <- diag(variance)
            variance <- variance.new
            return(variance)
                    }
variance.Npars <- 5    # number of estimated parameters in variance matrix
# estimation using tam.mml
mod1 <- TAM::tam.mml( dat, Q=Q, userfct.variance=variance.fct1,
             variance.Npars=variance.Npars, control=list(maxiter=maxit, QMC=TRUE,
                          snodes=snodes))
summary(mod1)

#*****************
# Model 2: Testlet model with correlated testlet effects
#*****************

# specify a testlet model with general factor g and testlet effects
# u_1,u_2,u_3 and u_4. Assume that Cov(g,u_t)=0 for all t=1,2,3,4.
# Additionally, assume that \sum_t,t' Cov( u_t, u_t')=0, i.e.
# the sum of all testlet covariances is equal to zero
#=> testlet effects are uncorrelated on average.

# set Cov(g,u_t)=0 and sum of all testlet covariances equals to zero
variance.fct2 <- function( variance ){
            ndim <- ncol(variance)
            variance.new <- matrix( 0, ndim, ndim )
            diag(variance.new) <- diag(variance)
            variance.new[1,2:ndim] <- variance.new[2:ndim,1] <- 0
            # calculate average covariance between testlets
            v1 <- variance[ -1, -1] - variance.new[-1,-1]
            M1 <- sum(v1) / ( ( ndim-1)^2 - ( ndim - 1))
            v1 <- v1 - M1
            variance.new[ -1, -1 ] <- v1
            diag(variance.new) <- diag(variance)
            variance <- variance.new
            return(variance)
                    }
variance.Npars <- 1 + 4 + (4*3)/2 - 1
# estimate model in TAM
mod2 <- TAM::tam.mml( dat, Q=Q, userfct.variance=variance.fct2,
                variance.Npars=variance.Npars,
                control=list(maxiter=maxit, QMC=TRUE, snodes=snodes) )
summary(mod2)

#*****************
# Model 3: Testlet model with correlated testlet effects (different identification)
#*****************

# Testlet model like in Model 2. But now the constraint is
# \sum _t,t' Cov(u_t, u_t') + \sum_t Var(u_t)=0, i.e.
# the sum of all testlet covariances and variances is equal to zero.
variance.fct3 <- function( variance ){
            ndim <- ncol(variance)
            variance.new <- matrix( 0, ndim, ndim )
            diag(variance.new) <- diag(variance)
            variance.new[1,2:ndim] <- variance.new[2:ndim,1] <- 0
            # calculate average covariance and variance between testlets
            v1 <- variance[ -1, -1]
            M1 <- mean(v1)
            v1 <- v1 - M1
            variance.new[ -1, -1 ] <- v1
            # ensure positive definiteness of covariance matrix
            eps <- 10^(-2)
            diag(variance.new) <- diag( variance.new) + eps
            variance.new <- psych::cor.smooth( variance.new )  # smoothing in psych
            variance <- variance.new
            return(variance)
                    }
variance.Npars <- 1 + 4 + (4*3)/2 - 1
# estimate model in TAM
mod3 <- TAM::tam.mml( dat, Q=Q, userfct.variance=variance.fct3,
                variance.Npars=variance.Npars,
                control=list(maxiter=maxit, QMC=TRUE, snodes=snodes) )
summary(mod3)

#*****************
# Model 4: Rasch subdimension model
#*****************

# The Rasch subdimension model is specified according to Brandt (2008).
# The fourth testlet effect is defined as u4=- (u1+u2+u3)
# specify an alternative Q-matrix with 4 dimensions
Q2 <- Q[,-5]
Q2[31:40,2:4] <- -1

# set Cov(g,u1)=Cov(g,u2)=Cov(g,u3)=0
variance.fixed <- rbind( c(1,2,0), c(1,3,0), c(1,4,0) )
# estimate model in TAM
mod4 <- TAM::tam.mml( dat, Q=Q2,variance.fixed=variance.fixed,
                control=list(maxiter=maxit, QMC=TRUE, snodes=snodes) )
summary(mod4)

#*****************
# Model 5: Higher-order model
#*****************

# A four-dimensional model with a higher-order factor is specified.
# F_t=a_t g + eps_g
Q3 <- Q[,-1]
# define fitting function using the lavaan package and ULS estimation
N0 <- nrow(dat)         # sample size of dataset
library(lavaan)        # requires lavaan package for fitting covariance
variance.fct5 <- function( variance ){
    ndim <- ncol(variance)
    rownames(variance) <- colnames(variance) <- paste0("F",1:ndim)
    lavmodel <- paste0(
        "FHO=~", paste0( paste0( "F", 1:ndim ), collapse="+" ) )
    lavres <- lavaan::cfa( model=lavmodel, sample.cov=variance, estimator="ULS",
                       std.lv=TRUE, sample.nobs=N0)
    variance.new <- fitted(lavres)$cov
    variance <- variance.new
    # print coefficients
    cat( paste0( "\n **** Higher order loadings: ",
            paste0( paste0( round( coef(lavres)[ 1:ndim ], 3 )), collapse=" ")
                        ), "\n")
    return(variance)
                    }
variance.Npars <- 4+4
# estimate model in TAM
mod5 <- TAM::tam.mml( dat, Q=Q3, userfct.variance=variance.fct5,
                variance.Npars=variance.Npars,
                control=list(maxiter=maxit, QMC=TRUE, snodes=snodes) )
summary(mod5)

#*****************
# Model 6: Generalized Rasch subdimension model (Brandt, 2012)
#*****************

Q2 <- Q[,-5]
Q2[31:40,2:4] <- -1
# fixed covariances
variance.fixed2 <- rbind( c(1,2,0), c(1,3,0), c(1,4,0)  )
# design matrix for item loading parameters
#      items x category x dimension x xsi parameter
E <- array( 0, dim=c( 40, 2, 4, 4 ) )
E[ 1:10, 2, c(1,2), 1 ] <- 1
E[ 11:20, 2, c(1,3), 2 ] <- 1
E[ 21:30, 2, c(1,4), 3 ] <- 1
E[ 31:40, 2, 1, 4 ] <- 1
E[ 31:40, 2, 2:4, 4 ] <- -1

# constraint on slope parameters, see Brandt (2012)
gammaconstr <- function( gammaslope ){
        K <- length( gammaslope)
        g1 <- sum( gammaslope^2  )
        gammaslope.new <- sqrt(K) / sqrt(g1) * gammaslope
        return(gammaslope.new)
                    }
# estimate model
mod6 <- TAM::tam.mml.3pl( dat, E=E, Q=Q2, variance.fixed=variance.fixed2,
           skillspace="normal", userfct.gammaslope=gammaconstr, gammaslope.constr.Npars=1,
           control=list(maxiter=maxit, QMC=TRUE, snodes=snodes ) )
summary(mod6)

#############################################################################
# EXAMPLE 18: Partial credit model with dimension-specific sum constraints
#             on item difficulties
#############################################################################

data(data.Students, package="CDM")
dat <- data.Students[, c( paste0("sc",1:4), paste0("mj",1:4) ) ]
# specify dimensions in Q-matrix
Q <- matrix( 0,  nrow=8, ncol=2 )
Q[1:4,1] <- Q[5:8,2] <- 1
# partial credit model with some constraint on item parameters
mod1 <- TAM::tam.mml( dat, Q=Q, irtmodel="PCM2", constraint="items")
summary(mod1)

#############################################################################
# EXAMPLE 19: Partial credit scoring: 0.5 points for partial credit items
#             and 1 point for dichotomous items
#############################################################################

data(data.timssAusTwn.scored)
dat <- data.timssAusTwn.scored
# extrcat item response data
dat <- dat[, grep("M03", colnames(dat) ) ]

# select items with do have maximum score of 2 (polytomous items)
ind <- which( apply( dat,  2, max, na.rm=TRUE )==2 )
I <- ncol(dat)
# define Q-matrix with scoring variant
Q <- matrix( 1, nrow=I, ncol=1 )
Q[ ind, 1 ] <- .5    # score of 0.5 for polyomous items

# estimate model
mod1 <- TAM::tam.mml( dat, Q=Q, irtmodel="PCM2", control=list(nodes=seq(-10,10,len=21)))
summary(mod1)

#############################################################################
# EXAMPLE 20: Specification of loading matrix in multidimensional model
#############################################################################

data(data.gpcm)
psych::describe(data.gpcm)
resp <- data.gpcm

# define three dimensions and different loadings of item categories
# on these dimensions in B loading matrix
I <- 3  # 3 items
D <- 3  # 3 dimensions
# define loading matrix B
# 4 categories for each item (0,1,2,3)
B <- array( 0, dim=c(I,4,D) )
for (ii in 1:I){
    B[ ii, 1:4, 1 ] <- 0:3
    B[ ii, 1,2 ] <- 1
    B[ ii, 4,3 ] <- 1
            }
dimnames(B)[[1]] <- colnames(resp)
B[1,,]
  ##   > B[1,,]
  ##        [,1] [,2] [,3]
  ##   [1,]    0    1    0
  ##   [2,]    1    0    0
  ##   [3,]    2    0    0
  ##   [4,]    3    0    1
#-- test run
mod1 <- TAM::tam.mml( resp, B=B, control=list( snodes=1000, maxiter=5)  )
summary(mod1)

# Same model with TAM::tam.mml.3pl instead

dim4 <- sum(apply(B, c(1, 3), function(x) any(!(x==0))))
E1 <- array(0, dim=c(dim(B), dim4))

kkk <- 0
for (iii in seq_len(dim(E1)[1])) {
    for (jjj in seq_len(dim(E1)[3])) {
        if (any(B[iii,, jjj] !=0)) {
            kkk <- kkk + 1
            E1[iii,, jjj, kkk] <- B[iii,, jjj]
        }
    }
}
if (kkk !=dim4) stop("Something went wrong in the loop, because 'kkk !=dim4'.")

mod2 <- TAM::tam.mml.3pl(resp, E=E1, est.some.slopes=FALSE, control=list(maxiter=50))
summary(mod2)

cor(mod1$person$EAP.Dim3, mod2$person$EAP.Dim3)
cor(mod1$person$EAP.Dim2, mod2$person$EAP.Dim2)
cor(mod1$person$EAP.Dim1, mod2$person$EAP.Dim1)
cor(mod1$xsi$xsi, mod2$xsi$xsi)

#############################################################################
# EXAMPLE 21: Acceleration of EM algorithm | dichotomous data
#############################################################################

N <- 1000      # number of persons
I <- 100       # number of items
set.seed(987)
# simulate data according to the Rasch model
dat <- sirt::sim.raschtype( stats::rnorm(N), b=seq(-2,2,len=I)  )
# estimate models
mod1n <- TAM::tam.mml( resp=dat, control=list( acceleration="none") )  # no acceler.
mod1y <- TAM::tam.mml( resp=dat, control=list( acceleration="Yu") )  # Yu acceler.
mod1r <- TAM::tam.mml( resp=dat, control=list( acceleration="Ramsay") )  # Ramsay acceler.
# compare number of iterations
mod1n$iter ; mod1y$iter ; mod1r$iter
# log-likelihood values
logLik(mod1n); logLik(mod1y) ; logLik(mod1r)

#############################################################################
# EXAMPLE 22: Acceleration of EM algorithm | polytomous data
#############################################################################

data(data.gpcm)
dat <- data.gpcm

# no acceleration
mod1n <- TAM::tam.mml.2pl( resp=dat, irtmodel="GPCM",
                control=list( conv=1E-4, acceleration="none") )
# Yu acceleration
mod1y <- TAM::tam.mml.2pl( resp=dat, irtmodel="GPCM",
                control=list( conv=1E-4, acceleration="Yu") )
# Ramsay acceleration
mod1r <- TAM::tam.mml.2pl( resp=dat, irtmodel="GPCM",
                control=list( conv=1E-4, acceleration="Ramsay") )
# number of iterations
mod1n$iter ; mod1y$iter ; mod1r$iter

#############################################################################
# EXAMPLE 23: Multidimensional polytomous Rasch model in which
#             dimensions are defined by categories
#############################################################################

data(data.Students, package="CDM")
dat <- data.Students[, grep( "act", colnames(data.Students) ) ]

# define multidimensional model in which categories of item define dimensions

# * Category 0 -> loading of one on Dimension 0
# * Category 1 -> no loadings
# * Category 2 -> loading of one on Dimension 2

# extract default design matrices
res <- TAM::designMatrices( resp=dat )
A <- res$A
B0 <- 0*res$B
# create design matrix B
B <- array( 0, dim=c( dim(B0)[c(1,2) ], 2  ) )
dimnames(B)[[1]] <- dimnames(B0)[[1]]
dimnames(B)[[2]] <- dimnames(B0)[[2]]
dimnames(B)[[3]] <- c( "Dim0", "Dim2" )
B[,1,1]  <- 1
B[,3,2]  <- 1

# estimate multdimensional Rasch model
mod1 <- TAM::tam.mml( resp=dat, A=A, B=B, control=list( maxiter=100) )
summary(mod1)

# alternative definition of B
# * Category 1: negative loading on Dimension 1 and Dimension 2
B <- array( 0, dim=c( dim(B0)[c(1,2) ], 2  ) )
dimnames(B)[[1]] <- dimnames(B0)[[1]]
dimnames(B)[[2]] <- dimnames(B0)[[2]]
dimnames(B)[[3]] <- c( "Dim0", "Dim2" )
B[,1, 1]  <- 1
B[,3, 2]  <- 1
B[,2, c(1,2)]  <- -1

# estimate model
mod2 <- TAM::tam.mml( resp=dat, A=A, B=B, control=list( maxiter=100) )
summary(mod2)

#############################################################################
# EXAMPLE 24: Sum constraint on item-category parameters in partial credit model
#############################################################################

data(data.gpcm,package="TAM")
dat <- data.gpcm

# check number of categories
c1 <- TAM::tam.ctt3(dat)

#--- fit with PCM
mod1 <- TAM::tam.mml( dat )
summary(mod1, file="mod1")

#--- fit with constraint on sum of categories
#** redefine design matrix
A1 <- 0*mod1$A
A1 <- A1[,, - dim(A1)[[3]]]
str(A1)
NP <- dim(A1)[[3]]
# define item category parameters
A1[1,2,1] <- A1[1,3,2] <- A1[1,4,3] <- -1
A1[2,2,4] <- A1[2,3,5] <- A1[2,4,6] <- -1
A1[3,2,7] <- A1[3,3,8] <- -1
A1[3,4,1:8] <- 1
# check definition
A1[1,,]
A1[2,,]
A1[3,,]

#** estimate model
mod2 <- TAM::tam.mml( dat, A=A1, beta.fixed=FALSE)
summary(mod2, file="mod2")

#--- compare model fit
IRT.compareModels(mod1, mod2 )  # -> equivalent model fit

#############################################################################
# EXAMPLE 25: Different GPCM parametrizations in IRT packages
#############################################################################

library(TAM)
library(mirt)
library(ltm)

data(data.gpcm, package="TAM")
dat <- data.gpcm

#*** TAM
mod1 <- TAM::tam.mml.2pl(dat, irtmodel="GPCM")
#*** mirt
mod2 <- mirt::mirt(dat, 1, itemtype="gpcm", verbose=TRUE)
#*** ltm
mod3 <- ltm::gpcm( dat, control=list(verbose=TRUE) )
mod3b <- ltm::gpcm( dat, control=list(verbose=TRUE), IRT.param=FALSE)

#-- comparison log likelihood
logLik(mod1)
logLik(mod2)
logLik(mod3)
logLik(mod3b)

#*** intercept parametrization (like in TAM)

# TAM
mod1$B[,2,1]   # slope
mod1$AXsi      # intercepts
# mirt
coef(mod2)
# ltm
coef(mod3b, IRT.param=FALSE)[, c(4,1:3)]

#*** IRT parametrization

# TAM
mod1$AXsi / mod1$B[,2,1]
# mirt
coef(mod2, IRTpars=TRUE)
# ltm
coef(mod3)[, c(4,1:3)]

#############################################################################
# EXAMPLE 26: Differential item functioning in multdimensional models
#############################################################################

data(data.ex08, package="TAM")
formulaA <- ~ item+female+item*female
resp <- data.ex08[["resp"]]
facets <- as.data.frame(data.ex08[["facets"]])

#***  Model 8a: investigate gender DIF in undimensional model
mod8a <- TAM::tam.mml.mfr(resp=resp, facets=facets, formulaA=formulaA)
summary(mod8a)

#*** multidimensional 2PL Model
I <- 10
Q <- array(0, dim=c(I, 3))
Q[cbind(1:I, c(rep(1, 3), rep(2, 3), rep(3, 4)))] <- 1
rownames(Q) <- colnames(resp)
mod3dim2pl <- TAM::tam.mml.2pl(resp=resp, Q=Q, irtmodel="2PL",
                          control=list(snodes=2000))

#*** Combine both approaches
thisRows <- gsub("-.*", "", colnames(mod8a$resp)) #select item names

#*** uniform DIF (note irtmodel="2PL.groups" & est.slopegroups)
mod3dim2pl_udiff <- TAM::tam.mml.2pl(resp=mod8a$resp, A=mod8a$A, Q=Q[thisRows, ],
                               irtmodel="2PL.groups",
                               est.slopegroups=as.numeric(as.factor(thisRows)),
                               control=list(snodes=2000))

#*** non-uniform DIF (?); different slope parameters for each item administered to each group
mod3dim2pl_nudiff <- TAM::tam.mml.2pl(resp=mod8a$resp, A=mod8a$A, Q=Q[thisRows, ],
                                irtmodel="2PL", control=list(snodes=2000))

#*** check results
print(mod8a$xsi)
print(mod3dim2pl_udiff$xsi)
summary(mod3dim2pl_nudiff)

#*** within item dimensionality (one item loads on two dimensions)
Q2 <- Q
Q2[4,1] <- 1

# uniform DIF (note irtmodel="2PL.groups" & est.slopegroups)
mod3dim2pl_udiff2 <- TAM::tam.mml.2pl(resp=mod8a$resp, A=mod8a$A, Q=Q2[thisRows, ],
                                irtmodel="2PL.groups",
                                est.slopegroups=as.numeric(as.factor(thisRows)),
                                control=list(snodes=2000))
print(mod8a$xsi)
print(mod3dim2pl_udiff2$xsi)
print(mod3dim2pl_udiff2$xsi)

#############################################################################
# EXAMPLE 27: IRT parameterization for generalized partial credit model (GPCM) in TAM
#############################################################################

#--- read item parameters
pars <- as.numeric(miceadds::scan.vec(
"0.19029 1.25309 0.51737 -1.77046 0.94803
  0.19407 1.22680 0.34986 -1.57666 1.29726
  -0.00888 1.07093 0.31662 -1.38755 1.14809
  -0.33810 1.08205 0.48490 -1.56696 0.79547
  -0.18866 0.99587 0.37880 -1.37468 0.81114" ))
pars <- matrix( pars, nrow=5, byrow=TRUE)
beta <- pars[,1]
alpha <- pars[,5]
tau <- pars[,2:4]

#--- data simulation function for GPCM
sim_gpcm_irt_param <- function(alpha, beta, tau, N, mu=0, sigma=1)
{
    theta <- stats::rnorm(N, mean=mu, sd=sigma)
    I <- length(beta)
    K <- ncol(tau)
    dat <- matrix(0, nrow=N, ncol=I)
    colnames(dat) <- paste0("I",1:I)
    for (ii in 1:I){
        probs <- matrix(0, nrow=N, ncol=K+1)
        for (kk in 1:K){
            probs[,kk+1] <- probs[,kk] + alpha[ii]*( theta - beta[ii] - tau[ii,kk] )
        }
        probs <- exp(probs)
        probs <- probs/rowSums(probs)
        rn <- stats::runif(N)
        cumprobs <- t(apply(probs,1,cumsum))
        for (kk in 1:K){
            dat[,ii] <- dat[,ii] + ( rn > cumprobs[,kk] )
        }
    }
    return(dat)
}

#-- simulate data
N <- 20000     # number of persons
set.seed(98)
dat1 <- sim_gpcm_irt_param(alpha=alpha, beta=beta, tau=tau, N=N, mu=0, sigma=1)
head(dat1)

#* generate design matrix for IRT parameterization
A1 <- TAM::.A.PCM2( resp=dat1)

#- estimate GPCM model
mod1 <- TAM::tam.mml.2pl( resp=dat1, irtmodel="GPCM", A=A1)
summary(mod1)

# compare true and estimated slope estimates (alpha)
cbind( alpha, mod1$B[,2,] )

# compare true and estimated item difficulties (beta)
cbind( beta, mod1$xsi$xsi[1:5] / mod1$B[,2,1] )

# compare true and estimated tau parameters
cbind( tau[,-3], matrix( mod1$xsi$xsi[-c(1:5)], nrow=5, byrow=TRUE ) / mod1$B[,2,1] )

## End(Not run)

TAM documentation built on July 2, 2020, 2:40 a.m.