The two-parameter partial credit model

Item \(i\) has \(K\) categories \(k = 0, \dots, K-1\), that is \(L = K - 1\) steps. The two-parameter partial credit model (TPPCM; Yu, 1991) gives every step \(l\) of Masters’ (1982) partial credit model its own discrimination \(a_{il}\):

\[ P(X_i = k \mid \theta) \propto \exp\Big(\sum_{l \le k} a_{il}(\theta - b_{il})\Big) = \exp\Big(\sum_{l \le k} (a_{il}\theta + c_{il})\Big), \qquad c_{il} = -a_{il} b_{il}. \]

The package labels parameters with words to avoid the clash of symbols between references and programs:

Note that d in mirt and psychotools is the intercept cumulated up to the category, not the step difficulty. The third column gives the names of the category boundary literature (Preston et al., 2011), which the DSF functions use: CBD (category boundary discrimination) \(= a_{il}\), CBL (category boundary location) \(= b_{il}\), IBD (intercept boundary difference) \(= c_{il} = -a_{il}b_{il}\). The mapping is documented once in ?tppcm-package.

The common polytomous models are special cases. The package describes them with a formula for the step discriminations, read as a two-way item × step layout on the log scale, like a two-way ANOVA:

\[ \log a_{il} = \mu + \alpha_i + \gamma_l + (\alpha\gamma)_{il}. \]

~ item:step and ~ item + step + item:step are the same model as ~ item * step. The models form a hierarchy, not a single ladder:

        ~ item * step     a_il                (interaction: items differ in shape)
              |
        ~ item + step     a_il = alpha_i gamma_l   (main effects: same shape)
          /        \
   ~ item           ~ step
   (alpha_i)        (gamma_l)
          \        /
             ~ 1          a_il = a

The interaction means that items differ in the shape of their step discriminations: the ratios \(a_{il}/a_{il'}\) differ between items. Rank-1 keeps one shape \(\gamma_l\) and lets items differ only in level \(\alpha_i\); the GPCM has a flat shape; ~ step has the same shape and level for all items. So the score test of the GPCM against ~ item + step asks whether there is a step effect, and the test of Rank-1 against ~ item * step whether there is an item × step interaction (Section 4). Rank-1 is multiplicative and cannot be written as a linear design for TAM: tppcm(dat, design = ~ item + step) stops and points to xxirt_tppcm() (Section 5). For the four linear models the log scale does not matter. It does presume \(a_{il} > 0\): Yu (1991) restricted the step discriminations to be positive; TAM does not impose it (xxirt_tppcm() does).

The same item × step language describes the locations in TAM’s and ConQuest’s facet formulas: formulaA = ~ item + step is the rating scale model (step parameters common to all items) and ~ item + item:step the partial credit model. design = ~ step is the rating-scale idea applied to the discriminations. With the latent variance fixed at 1, “one common \(a\)” (~ 1) is the same model as the usual PCM with \(a = 1\) and a free variance.

The saturated TPPCM is a reparameterization of the nominal response model (NRM; Bock, 1972): adjacent differences of the NRM category slopes are the step discriminations. The saturated model therefore needs no new estimator. The value of the TPPCM parameterization is that every restriction in the table can be written down and tested directly.

What the package does and does not do

tppcm does not estimate. Estimation stays with existing software:

• TAM: tam.mml.3pl() with the design produced by tppcm(); its slope parameters are then exactly \(a_{il}\). The PCM and GPCM can also be fitted with tam.mml() and tam.mml.2pl().
• sirt: the two things TAM cannot express, the product form (bilinear) and the bound \(a_{il} \ge 0\), through xxirt_tppcm().

tppcm adds what TAM lacks:

Every function takes the fitted TAM or xxirt object directly; the responses are read from the object.

Data

The four mathematics self-concept items of CDM::data.Students (Austrian 8th graders, \(N = 2400\)) are scored 0–3 (do not agree at all … completely agree); items sc2 and sc3 are reversed in the data set. About 2% of the responses are missing; TAM’s marginal maximum likelihood uses the observed responses as they are.

library(tppcm)   # also attaches TAM

data(data.Students, package = "CDM")
sc <- data.Students[, c("sc1", "sc2", "sc3", "sc4")]
sapply(sc, table, useNA = "ifany")
#>       sc1 sc2  sc3 sc4
#> 0     718 910 1060 469
#> 1    1081 880  728 963
#> 2     476 429  358 665
#> 3      77 128  189 240
#> <NA>   48  53   65  63

The saturated TPPCM

tppcm(dat) returns the design array E for tam.mml.3pl(); the slope parameters gammaslope are then exactly \(a_{il}\). The latent distribution has to be fixed at \(N(0,1)\) to identify the scale, so est.variance = FALSE is required.

m_tppcm <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE)
matrix(m_tppcm$gammaslope, nrow = 4, byrow = TRUE)    # rows = items, columns = steps
#>          [,1]     [,2]      [,3]
#> [1,] 2.943619 3.793098 1.2319493
#> [2,] 2.266647 1.623264 0.6768669
#> [3,] 2.577356 2.261299 1.0218987
#> [4,] 2.060255 2.634662 1.1709056

It is the same model as TAM’s NRM, only parameterized differently, so the log-likelihoods agree:

m_nrm <- tam.mml.2pl(sc, irtmodel = "2PL", verbose = FALSE)
c(TPPCM = logLik(m_tppcm), NRM = logLik(m_nrm))
#>     TPPCM       NRM 
#> -9419.699 -9419.699

When tam.mml.3pl() does not converge

On typical rating scales the saturated design can make TAM’s EM algorithm diverge. The design is cumulative, \(a_{il}\) enters the slopes of categories \(l, \dots, K_i - 1\), so the steps of an item are strongly correlated in the likelihood. TAM updates each slope with its own gradient and second derivative, without the cross-derivatives, and with 5 or more steps per item and sparse low categories this step overshoots: the deviance rises from iteration to iteration, slopes drift to the bound \(\pm 9.99\) and maxiter is exhausted. The psych::bfi agreeableness items (6 categories, few responses in the lowest ones) show it:

bfi <- psych::bfi[, paste0("A", 1:5)] - 1       # bfi is coded 1-6; tppcm needs 0-5
bfi$A1 <- 5 - bfi$A1                            # reverse-keyed item
m_bad <- tam.mml.3pl(bfi, E = tppcm(bfi), est.variance = FALSE, verbose = FALSE,
                     control = list(maxiter = 5000))
m_bad$iter                                      # 5000: not converged
range(m_bad$gammaslope)                         # -9.99 and 9.99: at the bound

Such a fit is not an estimate. The remedy is to damp the M-step with max.increment, the largest change TAM allows a parameter per iteration (default 1): with 0.1 the same fit converges in a few hundred iterations and a few seconds. For the saturated design on real data, use control = list(maxiter = 3000, max.increment = 0.1) as the standard call.

m_ok <- tam.mml.3pl(bfi, E = tppcm(bfi), est.variance = FALSE, verbose = FALSE,
                    control = list(maxiter = 3000, max.increment = 0.1))
m_nrm_bfi <- tam.mml.2pl(bfi, irtmodel = "2PL", verbose = FALSE)   # same model, converges as is
c(TPPCM = logLik(m_ok), NRM = logLik(m_nrm_bfi))

Always compare with the NRM as above: its log-likelihood is the maximum of the same model, and tam.mml.2pl() reaches it without tuning. If the damped fit stays below it, it sits in a local optimum (on these items it does, by about 60 units). Then pass the NRM fit itself to the functions of the package, or, as a suggestion, start the damped fit from the NRM solution (gammaslope = as.vector(t(get_parts(m_nrm_bfi)$a)), xsi.inits = cbind(seq_along(m_nrm_bfi$xsi$xsi), m_nrm_bfi$xsi$xsi); this needs the same number of categories in every item) and check its log-likelihood again. Quasi-Monte Carlo integration (QMC = TRUE) does not help; the problem is the step length, not the quadrature. The restricted designs of the next section have fewer and less correlated slopes and converge with TAM’s defaults.

The common designs

design gives the four linear models of the hierarchy in Section 1:

m_gpcmE <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ item), est.variance = FALSE, verbose = FALSE)
m_step  <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ step), est.variance = FALSE, verbose = FALSE)
m_pcmE  <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ 1),  est.variance = FALSE, verbose = FALSE)
attr(tppcm(sc, design = ~ step), "index")     # one parameter per column (step)
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    1    2    3
#> [3,]    1    2    3
#> [4,]    1    2    3
m_step$gammaslope                    # gamma_1, gamma_2, gamma_3
#>    step1    step2    step3 
#> 2.458218 2.437532 1.019129

The common slope of ~ 1 equals the latent standard deviation of TAM’s own PCM; the log-likelihoods differ only through numerical integration (see the box below):

m_pcm <- tam.mml(sc, verbose = FALSE)
c(common_a = unname(m_pcmE$gammaslope), tam_sd = sqrt(m_pcm$variance[1]))
#> common_a   tam_sd 
#> 2.174145 2.096062
c(tppcm = logLik(m_pcmE), tam = logLik(m_pcm))
#>     tppcm       tam 
#> -9501.691 -9508.697

Items may have different numbers of categories. The TIMSS data below mix dichotomous items with items scored 0–2; steps that do not exist are NA in the index and blank in the report:

data(data.timssAusTwn.scored, package = "TAM")
timss <- data.timssAusTwn.scored[, 1:11]
attr(tppcm(timss), "index")
#>       [,1] [,2]
#>  [1,]    1   NA
#>  [2,]    2   NA
#>  [3,]    3    4
#>  [4,]    5    6
#>  [5,]    7   NA
#>  [6,]    8   NA
#>  [7,]    9   10
#>  [8,]   11   12
#>  [9,]   13   NA
#> [10,]   14   NA
#> [11,]   15   NA
m_timss <- tam.mml.3pl(timss, E = tppcm(timss), est.variance = FALSE, verbose = FALSE)
irt_pars(m_timss, fit = FALSE, se = FALSE)
#> Model: saturated TPPCM, design ~ item * step (tam.mml.3pl) | kernel disc * (theta - diff)
#> 
#> $items
#>          disc1 disc2  diff1  diff2
#> M032166  1.095       -0.927       
#> M032721  0.610       -0.058       
#> M032757  0.732 1.671  1.768 -1.762
#> M032760A 2.190 3.566  0.740 -0.513
#> M032760B 3.232        0.358       
#> M032760C 4.242        0.537       
#> M032761  3.247 2.659  0.361  0.493
#> M032692  1.442 1.615  1.875 -0.911
#> M032626  1.946       -0.362       
#> M032595  1.939       -0.633       
#> M032673  1.587       -0.666       
#> 
#> $group
#>  group MEAN VAR
#>      1    0   1

Restrictions at the step level

More general restrictions go into index: an item × step matrix of integers, where steps with the same index share one parameter. For example, items 1 and 2 saturated and items 3 and 4 GPCM:

idx <- rbind(1:3,
             4:6,
             rep(7, 3),
             rep(8, 3))
m_mix <- tam.mml.3pl(sc, E = tppcm(sc, index = idx), est.variance = FALSE, verbose = FALSE)

To fix a step discrimination use TAM’s own gammaslope.fixed with the index of the step and the value; here step 2 of item 1 (index 2) is fixed at 3:

m_fix <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE,
                     gammaslope.fixed = cbind(2, 3))

The GPCM can also come from TAM’s own function. Likelihood-ratio tests along the hierarchy:

m_gpcm <- tam.mml.2pl(sc, irtmodel = "GPCM", verbose = FALSE)
lr <- IRT.compareModels(m_pcmE, m_gpcm, m_step, m_mix, m_tppcm)$LRtest
nested <- !paste(lr$Model1, lr$Model2) %in% c("m_gpcm m_step", "m_step m_gpcm", "m_step m_mix", "m_mix m_step")
lr[nested, ]
#>    Model1  Model2      Chi2 df            p
#> 1  m_pcmE  m_gpcm  45.01727  3 9.174822e-10
#> 2  m_pcmE  m_step  97.42024  2 0.000000e+00
#> 3  m_pcmE   m_mix 101.39636  7 0.000000e+00
#> 4  m_pcmE m_tppcm 163.98340 11 0.000000e+00
#> 6  m_gpcm   m_mix  56.37909  4 1.669798e-11
#> 7  m_gpcm m_tppcm 118.96613  8 0.000000e+00
#> 9  m_step m_tppcm  66.56317  9 7.156853e-11
#> 10  m_mix m_tppcm  62.58704  4 8.288925e-13

Only nested pairs are shown: the GPCM and the step model do not contain each other, nor do the step model and m_mix, and the likelihood ratio between them is meaningless. On these data every step is worth it: the step model beats the PCM by 97 units on 2 df (steps differ), the GPCM by 45 on 3 df (items differ), and the saturated model beats both.

Which PCM? For likelihood-ratio comparisons use tam.mml.3pl(E = tppcm(dat, design = ~ 1)) (common slope, variance 1), so that all models share one parameterization and one set of quadrature nodes. tam.mml() estimates the variance instead; the two parameterizations are nested mathematically, but with TAM’s default 21 nodes their numerical integration errors differ. Above, the two PCM log-likelihoods differ by 7 units; with more items the likelihood ratio against a slope-parameterized model can even become negative. Alternatively use more nodes, e.g. control = list(nodes = seq(-6, 6, len = 61)).

irt_pars(): a compact parameter table

Modeled on mirt’s coef(mod, IRTpars = TRUE, simplify = TRUE): a fit header, one row per item, estimates and SEs in separate tables, and the parameters of the latent distribution at the end. Any unidimensional TAM object (PCM, RSM, GPCM, NRM, tam.mml.3pl()) and any xxirt_tppcm() object can be passed.

irt_pars(m_tppcm)
#> Model: saturated TPPCM, design ~ item * step (tam.mml.3pl) | kernel disc * (theta - diff)
#> Fit: N = 2400, logLik = -9419.70, npar = 24, AIC = 18887.40, BIC = 19026.20
#>      M2(42) = 616.96, p = .000, RMSEA = .076 [.070, .081], SRMSR = .033, CFI = .913, TLI = .888
#> 
#> $items
#>     disc1 disc2 disc3  diff1 diff2 diff3
#> sc1 2.944 3.793 1.232 -0.574 0.776 2.813
#> sc2 2.267 1.623 0.677 -0.267 0.983 2.874
#> sc3 2.577 2.261 1.022 -0.058 0.909 1.802
#> sc4 2.060 2.635 1.171 -1.017 0.357 1.818
#> 
#> $se
#>     disc1 disc2 disc3 diff1 diff2 diff3
#> sc1 0.219 0.354 0.283 0.033 0.036 0.328
#> sc2 0.152 0.140 0.174 0.034 0.058 0.465
#> sc3 0.187 0.217 0.203 0.033 0.050 0.152
#> sc4 0.176 0.191 0.159 0.048 0.034 0.131
#> 
#> $group
#>  group MEAN VAR
#>      1    0   1

The two Fit: lines give

• the log-likelihood, number of parameters, AIC and BIC (identical to TAM’s IRT.IC());
• limited-information fit: M2, RMSEA2 with its 90% interval, SRMSR, CFI and TLI (see m2() below).

The report makes two things visible at once. The third step of every item has a markedly lower discrimination than the first two (the move from “rather agree” to “completely agree” separates students less sharply), and the standard errors of the third-step difficulties are large because few students reach that category of the reversed items.

IRTpars = FALSE switches to the slope–intercept form slope * theta + int:

irt_pars(m_gpcm, IRTpars = FALSE)
#> Model: GPCM, design ~ item | kernel slope * theta + int
#> Fit: N = 2400, logLik = -9479.18, npar = 16, AIC = 18990.36, BIC = 19082.90
#>      M2(50) = 845.96, p = .000, RMSEA = .081 [.077, .086], SRMSR = .032, CFI = .879, TLI = .870
#> 
#> $items
#>     slope1 slope2 slope3  int1   int2   int3
#> sc1  2.944  2.944  2.944 1.677 -2.475 -5.922
#> sc2  1.766  1.766  1.766 0.474 -1.638 -3.308
#> sc3  2.159  2.159  2.159 0.066 -1.962 -3.225
#> sc4  2.053  2.053  2.053 2.083 -0.809 -3.023
#> 
#> $se
#>     slope1 slope2 slope3  int1  int2  int3
#> sc1  0.169  0.169  0.169 0.121 0.148 0.313
#> sc2  0.093  0.093  0.093 0.073 0.097 0.175
#> sc3  0.129  0.129  0.129 0.078 0.124 0.204
#> sc4  0.108  0.108  0.108 0.117 0.080 0.158
#> 
#> $group
#>  group MEAN VAR
#>      1    0   1

Fixed parameters get no SE. TAM’s PCM fixes \(a = 1\) and estimates the variance, so \(a\) has no SE and the variance is reported with one:

irt_pars(m_pcm)
#> Model: PCM, design ~ 1 (TAM: slope 1, free variance) | kernel disc * (theta - diff)
#> Fit: N = 2400, logLik = -9508.70, npar = 13, AIC = 19043.39, BIC = 19118.58
#>      M2(53) = 833.78, p = .000, RMSEA = .078 [.074, .083], SRMSR = .046, CFI = .882, TLI = .879
#> 
#> $items
#>     disc1 disc2 disc3  diff1 diff2 diff3
#> sc1 1.000 1.000 1.000 -1.270 1.912 4.549
#> sc2 1.000 1.000 1.000 -0.584 1.873 3.811
#> sc3 1.000 1.000 1.000 -0.056 1.925 3.148
#> sc4 1.000 1.000 1.000 -2.133 0.826 3.084
#> 
#> $se
#>     disc1 disc2 disc3 diff1 diff2 diff3
#> sc1                   0.078 0.084 0.161
#> sc2                   0.075 0.087 0.139
#> sc3                   0.074 0.092 0.127
#> sc4                   0.090 0.076 0.111
#> 
#> $group
#>  group MEAN  VAR SE_VAR
#>      1    0 4.39  0.229

Standard errors

irt_pars() computes standard errors from the observed information of the whole model (Louis’ formula), so the covariance between slopes and thresholds is taken into account.

irt_pars(long = TRUE) returns a long table with confidence intervals:

head(irt_pars(m_tppcm, long = TRUE, level = 0.90))
#> saturated TPPCM, design ~ item * step (tam.mml.3pl), 90% intervals 
#>           item step  par    est    se  lower  upper
#> sc1_disc1  sc1    1 disc  2.944 0.219  2.584  3.304
#> sc1_disc2  sc1    2 disc  3.793 0.354  3.210  4.376
#> sc1_disc3  sc1    3 disc  1.232 0.283  0.767  1.697
#> sc1_diff1  sc1    1 diff -0.574 0.033 -0.629 -0.519
#> sc1_diff2  sc1    2 diff  0.776 0.036  0.717  0.835
#> sc1_diff3  sc1    3 diff  2.813 0.328  2.273  3.353

Overall fit: m2()

m2() computes the limited-information statistics of Maydeu-Olivares and Joe (2006). m2() differentiates the margins with respect to the saturated parameters and multiplies by the Jacobian from get_parts(), so it works for any step design.

m2(m_gpcm)            # full M2: all univariate and bivariate margins
#> M2 = 845.96, df = 50, p = <2e-16
#> RMSEA = 0.081 [0.077, 0.086], SRMSR = 0.032, CFI = 0.879, TLI = 0.870
#> (missing responses: margins from available cases, valid under MCAR only)
m2(m_tppcm)
#> M2 = 616.96, df = 42, p = <2e-16
#> RMSEA = 0.076 [0.070, 0.081], SRMSR = 0.033, CFI = 0.913, TLI = 0.888
#> (missing responses: margins from available cases, valid under MCAR only)
m2(m_gpcm, "C2")      # univariate margins plus bivariate cross products
#> C2 = 71.57, df = 2, p = 2.87e-16
#> RMSEA = 0.120 [0.097, 0.145], SRMSR = 0.032, CFI = 0.985, TLI = 0.954
#> (missing responses: margins from available cases, valid under MCAR only)

• "M2" uses all univariate and bivariate margins and has the most degrees of freedom. The default "auto" uses it when there are at most 2000 margins and switches to "C2" otherwise.
• "C2" (Cai & Monroe, 2014) and "M2*" reduce the margins to moments and have far fewer degrees of freedom; for a saturated TPPCM with few items they can have df ≤ 0 and cannot be computed.
• RMSEA2 \(=\sqrt{\max(M_2-df,0)/(N\,df)}\); SRMSR is the root mean square difference between observed and model-implied item-score correlations; CFI and TLI use the independence model (independent items with free category proportions) as the baseline.
• One group only. With missing data every margin is computed from the respondents who answered the items involved and the covariance matrix is adjusted for the pairwise overlap; this is valid under missing completely at random (MCAR).
• When step discriminations differ within items, the full M2 can reject the GPCM while C2 does not: C2 and M2* collapse the categories into moments of the item scores, and differences between steps are averaged away. To check step discriminations use the full M2 or the score tests of Section 4.

Information: step_info()

The Fisher information of an item equals the conditional variance of the item score on the step scale, \(A_i(X) = \sum_{l \le X} a_{il}\), and splits into the contributions of the steps:

\[ I_i(\theta) = \mathrm{Var}(A_i(X) \mid \theta) = \sum_l c_{il}(\theta), \qquad c_{il}(\theta) = a_{il} \sum_m a_{im}\{S_{\max(l,m)} - S_l S_m\},\quad S_l = P(X_i \ge l \mid \theta). \]

When all \(a_{il} \ge 0\) every step contributes non-negatively. Item and test information agree with TAM’s IRT.informationCurves(); the step contributions are new.

inf <- step_info(m_tppcm)
inf
#> Fisher information on 161 theta values in [-4, 4]
#> max: maximum of the item (test) information and where it is reached;
#> max_step*: maximum of each step's contribution (each at its own theta)
#>        max theta_at_max max_step1 max_step2 max_step3
#> sc1  4.019         0.75     2.182     3.736     0.380
#> sc2  1.757        -0.10     1.401     0.783     0.139
#> sc3  2.627         0.35     1.831     1.574     0.366
#> sc4  2.357         0.30     1.099     1.930     0.363
#> test 9.623         0.60
op <- par(mfrow = c(1, 2)); plot(inf); plot(inf, item = "sc1"); par(op)

Per-item Wald tests: are the step discriminations of this item equal?

Under the saturated model, test \(a_{i1} = a_{i2} = a_{i3}\) for every item, that is whether the item follows a GPCM:

wald_test(m_tppcm)
#> Wald tests of equal step discriminations within items (H0: the item follows the GPCM) | p_adj: Holm 
#>  item statistic df      p  p_adj
#>   sc1    39.622  2 <1e-04 <1e-04
#>   sc2    49.246  2 <1e-04 <1e-04
#>   sc3    34.257  2 <1e-04 <1e-04
#>   sc4    36.870  2 <1e-04 <1e-04

Items that are already restricted to equal steps are skipped:

wald_test(m_mix)     # items 3 and 4 are GPCM items in m_mix
#> Wald tests of equal step discriminations within items (H0: the item follows the GPCM) | p_adj: Holm 
#>  item statistic df      p  p_adj
#>   sc1    26.166  2 <1e-04 <1e-04
#>   sc2    35.920  2 <1e-04 <1e-04

Score tests from the GPCM

The Wald test needs the saturated model. A score test needs only the restricted model (here the GPCM) and still tests every direction of relaxation, and it can point to the item and the step where the restriction fails.

GPCM against the saturated TPPCM:

score_test(m_gpcm)
#> Score test (Louis information) of the fitted GPCM, design ~ item against the saturated TPPCM (~ item * step)
#>   S = 213.1316, df = 8, p = < 2.2e-16

GPCM against Rank-1: is there a step effect? This test relaxes only the step pattern \(\gamma_2, \gamma_3\) shared by all items, df = 2 (against can be any larger model of the hierarchy, e.g. ~ item or ~ step for a PCM fit):

score_test(m_gpcm, against = ~ item + step)
#> Score test (Louis information) of the fitted GPCM, design ~ item against the Rank-1 product form (~ item + step)
#>   S = 172.7664, df = 2, p = < 2.2e-16

Per-item tests relax all step discriminations of item \(i\):

item_test(m_gpcm)
#> Score tests for freeing the step discriminations of each item in the fitted GPCM, design ~ item | p_adj: Holm 
#>  item statistic df      p  p_adj
#>   sc1    33.273  2 <1e-04 <1e-04
#>   sc2    40.383  2 <1e-04 <1e-04
#>   sc3    19.790  2 <1e-04 <1e-04
#>   sc4    27.045  2 <1e-04 <1e-04

Modification indices (MI) and expected parameter changes (EPC) per step:

mi(m_gpcm)
#> Modification indices (1-df score tests) for freeing single step discriminations of the fitted GPCM, design ~ item | p_adj: Holm 
#>        par     MI df       p    p_adj    EPC
#>  sc2_disc3 31.477  1 < 1e-04  < 1e-04 -1.657
#>  sc1_disc3 28.683  1 < 1e-04  < 1e-04 -3.477
#>  sc2_disc1 22.547  1 < 1e-04  < 1e-04  0.839
#>  sc4_disc3 21.419  1 < 1e-04  < 1e-04 -1.226
#>  sc3_disc3 18.597  1 < 1e-04 0.000129 -1.556
#>  sc4_disc2 18.128  1 < 1e-04 0.000145  0.815
#>  sc1_disc2 14.961  1 0.00011 0.000659  1.089
#>  sc3_disc1  7.111  1 0.00766 0.038302  0.610
#>  sc1_disc1  2.694  1 0.10073 0.402907 -0.491
#>  sc2_disc2  2.056  1 0.15160 0.454811 -0.261
#>  sc4_disc1  0.983  1 0.32148 0.642962 -0.206
#>  sc3_disc2  0.005  1 0.94491 0.944906  0.017

Reading the EPC: if this step were freed, its discrimination is expected to be this much higher (positive) or lower (negative) than the common GPCM slope of the item. Here every third step has a large negative EPC and the first two steps positive ones: the pattern is the same in all four items, which is exactly what the product form describes, one step pattern \(\gamma\) shared by the items and an item-specific size \(\alpha_i\).

Impact is not DIF: model the group distributions first

The TIMSS 2011 mathematics items of TAM::data.timssAusTwn.scored were answered by 1050 Australian and 719 Taiwanese fourth graders; gender is also recorded. dif_test() takes the fitted model and the covariate. parm chooses the parameters (default c("disc", "diff"), or only "disc" or "diff"), by the unit of testing:

• "param": every parameter separately, e.g. each step discrimination;
• "step": the discrimination and location of every step jointly (2 df), independent of the parameterization;
• "item": the parameters of an item jointly;
• "block": all selected parameters jointly.

When the covariate is unrelated to ability, a single fit to the whole sample suffices. For gender:

cov <- data.timssAusTwn.scored[, c("IDCNTRY", "ITSEX")]
cov$country <- factor(cov$IDCNTRY, labels = c("AUS", "TWN"))
cov$sex     <- factor(cov$ITSEX,   labels = c("girl", "boy"))
sp_timss <- get_parts(m_timss)                  # compute once, reuse below
dif_test(sp_timss, cov$sex, by = "item")
#> single-group fit: a difference in ability between the groups of 'covariate' (impact) is reported as DIF in the difficulties; refit with group = covariate if that is possible (see ?dif_test)
#> Score-based DIF tests along the covariate (LM functional) | unit: item (disc + diff) | p_adj: Holm 
#>      unit npar statistic df      p p_adj
#>   M032166    2     2.243  2 0.3257 1.000
#>   M032721    2     0.163  2 0.9218 1.000
#>   M032757    4     2.531  4 0.6390 1.000
#>  M032760A    4     9.896  4 0.0422 0.464
#>  M032760B    2     2.760  2 0.2515 1.000
#>  M032760C    2     1.768  2 0.4132 1.000
#>   M032761    4     6.739  4 0.1503 1.000
#>   M032692    4     3.909  4 0.4185 1.000
#>   M032626    2     1.001  2 0.6064 1.000
#>   M032595    2     1.149  2 0.5631 1.000
#>   M032673    2     0.508  2 0.7758 1.000

Gender shows no DIF after the Holm correction (p_adj; adjust = chooses another method of p.adjust()). On a single-group fit the case-wise scores are ordered by the covariate and aggregated with a functional chosen from its type (Merkle & Zeileis, 2013): a factor gives the LM test, an ordered factor the weighted maximum LM test ("maxLMo"), a numeric covariate the Cramér–von Mises type "CvM"; functional = overrides the choice. The message is a reminder of the condition above: a single-group fit is only valid when the groups do not differ in ability, which the next example violates.

dif_test(sp_timss, cov$country, parm = "diff", by = "item")
#> single-group fit: a difference in ability between the groups of 'covariate' (impact) is reported as DIF in the difficulties; refit with group = covariate if that is possible (see ?dif_test)
#> Score-based DIF tests along the covariate (LM functional) | unit: item (diff) | p_adj: Holm 
#>      unit npar statistic df      p    p_adj
#>   M032166    1     1.671  1 0.1962 0.323111
#>   M032721    1    16.213  1 <1e-04 0.000226
#>   M032757    2     7.582  2 0.0226 0.067707
#>  M032760A    2    20.776  2 <1e-04 0.000154
#>  M032760B    1   143.165  1 <1e-04  < 1e-04
#>  M032760C    1   194.653  1 <1e-04  < 1e-04
#>   M032761    2   245.685  2 <1e-04  < 1e-04
#>   M032692    2     3.646  2 0.1616 0.323111
#>   M032626    1   202.791  1 <1e-04  < 1e-04
#>   M032595    1    61.761  1 <1e-04  < 1e-04
#>   M032673    1    40.352  1 <1e-04  < 1e-04

Almost every item is flagged, with enormous statistics. This is not DIF but impact: Taiwanese students have a much higher mean ability, and the single-group model assumes one \(N(0,1)\) distribution for everybody. The scores of all difficulty parameters then move with the covariate. The same holds for any covariate correlated with \(\theta\).

The remedy is to let the model have the group difference: fit the multigroup model with group = (latent means and variances of groups 2, …, G estimated, item parameters common). All DSF tests below start from this fit:

m_mg <- tam.mml.3pl(timss, E = tppcm(timss), group = cov$country,
                    est.variance = FALSE, verbose = FALSE)
irt_pars(m_mg, fit = FALSE)$group
#>   group     MEAN      VAR    SE_VAR    SE_MEAN
#> 1   AUS 0.000000 1.000000        NA         NA
#> 2   TWN 1.540193 1.439903 0.1573386 0.07599146

For a numeric covariate correlated with ability use the latent regression Y = of TAM in the same way (the regression coefficients are then treated as fixed; the function warns).

A DSF workflow for known groups

1. Screen every step. On the multigroup fit, dif_test(by = "step") tests the CBD and location of every step jointly. It is the Lagrange multiplier test for DIF of Glas (1998), computed as the score test of the model in which the parameters of the step are split by group (2 df for two groups). The default info = "robust" is the generalized score test of Boos (1992), which keeps its level when the model does not fit exactly:

st <- dif_test(m_mg, cov$country, by = "step")
st
#> LM tests of DIF between the 2 groups of the covariate (robust information) | unit: step (discrimination and location jointly) | p_adj: Holm 
#>        unit npar statistic df       p    p_adj
#>   M032166:1    2    80.349  2 < 1e-04  < 1e-04
#>   M032721:1    2     6.175  2 0.04561 0.182444
#>   M032757:1    2    11.594  2 0.00304 0.024297
#>   M032757:2    2    10.055  2 0.00655 0.045876
#>  M032760A:1    2     7.645  2 0.02187 0.131233
#>  M032760A:2    2    19.365  2 < 1e-04 0.000561
#>  M032760B:1    2    28.579  2 < 1e-04  < 1e-04
#>  M032760C:1    2    21.409  2 < 1e-04 0.000224
#>   M032761:1    2     5.767  2 0.05593 0.182444
#>   M032761:2    2     3.347  2 0.18755 0.375091
#>   M032692:1    2    67.169  2 < 1e-04  < 1e-04
#>   M032692:2    2    99.273  2 < 1e-04  < 1e-04
#>   M032626:1    2    65.576  2 < 1e-04  < 1e-04
#>   M032595:1    2     2.791  2 0.24765 0.375091
#>   M032673:1    2     7.123  2 0.02839 0.141961
flagged <- st$unit[st$p_adj < 0.05]
flagged
#> [1] "M032166:1"  "M032757:1"  "M032757:2"  "M032760A:2" "M032760B:1" "M032760C:1"
#> [7] "M032692:1"  "M032692:2"  "M032626:1"

Nine of the 15 boundaries differ between the countries after Holm’s correction.

2. Decompose the flagged steps. dsf() refits the multigroup model and releases one parameter of a step at a time, with all other steps as anchors and the latent means and variances estimated (lavaan’s group.partial). The location row (ibd) releases the IBD of the step; the CBD row releases the CBD with the location of the step free, so it does not depend on the origin of the scale and location DSF does not leak into it. The two LR statistics add up to the joint 2-df LR test of the step. The units of dif_test(by = "step") can be passed as they are; every test refits the model, so here we decompose two steps:

dsf(m_mg, items = c("M032166:1", "M032760C:1"), verbose = FALSE)
#> LR tests of DSF, one step parameter released at a time vs the scalar baseline (CBD and IBD equal) | reference group: AUS | p_adj: Holm
#>      item step par    AUS    TWN     LR df       p    p_adj
#>   M032166    1 cbd  1.243  1.206  0.045  1 0.83179 0.831788
#>   M032166    1 ibd  0.907 -0.520 80.471  1 < 1e-04  < 1e-04
#>  M032760C    1 cbd  4.484  2.264 17.088  1 < 1e-04 0.000107
#>  M032760C    1 ibd -3.738 -4.430  9.732  1 0.00181 0.003621

Item M032166 is harder for Taiwanese students than for Australian students of the same ability (its IBD, the intercept, is lower; LR = 80), but it discriminates equally well in both countries (CBD LR = 0.05). Item M032760C in addition discriminates less sharply in Taiwan (CBD 4.5 vs 2.3; LR = 17). dsf(m_mg, items = flagged) runs all nine (about a minute).

3. Omnibus tests: configural, metric, scalar. dsf_design() builds the multigroup TPPCM for tam.mml.3pl() with lavaan-style constraints across groups. The CBD plays the role of a factor loading and the IBD that of an intercept:

fit_mg <- function(d) tam.mml.3pl(d$resp, A = d$A, E = d$E, group = cov$country,
                                  variance.fixed = d$variance.fixed,
                                  beta.fixed = d$beta.fixed, verbose = FALSE)
configural <- fit_mg(dsf_design(timss, cov$country, group.equal = NULL))
metric     <- fit_mg(dsf_design(timss, cov$country, group.equal = "cbd"))
scalar     <- fit_mg(dsf_design(timss, cov$country))
partial    <- fit_mg(dsf_design(timss, cov$country, group.partial = flagged))
IRT.compareModels(configural, metric, scalar, partial)$LRtest
#>    Model1     Model2      Chi2 df            p
#> 1  metric configural  45.39690 14 3.516314e-05
#> 2  scalar configural 377.30862 28 0.000000e+00
#> 3 partial configural  15.71242 10 1.081659e-01
#> 4  scalar     metric 331.91172 14 0.000000e+00
#> 5  metric    partial  29.68448  4 5.674311e-06
#> 6  scalar    partial 361.59620 18 0.000000e+00

The scalar model is the multigroup fit m_mg itself (same log-likelihood). Both the CBDs (metric vs configural) and the IBDs (scalar vs metric) differ between the countries. Releasing the nine flagged steps (partial, group.partial = flagged) leaves no significant difference to the configural model (row partial–configural), while the scalar model is clearly worse than the partial one. Read only nested pairs: metric and partial are not nested.

Why the joint test is the screen

• Clean null, invariant. The joint test of a step asks whether the boundary is the same in all groups. Its null hypothesis does not depend on how the boundary is parameterized: (CBD, CBL) and (CBD, IBD) give the same statistic. This is the logic of Lord’s (1980) \(\chi^2\) test and of the Wald tests of Woods, Cai and Wang (2013), which test all parameters of an item jointly, and of Stark, Chernyshenko and Drasgow (2006), who test loadings and intercepts simultaneously. Thissen (2024) discusses grouping DIF tests by parameter or by item.
• The 1-df rows leak. With by = "param" the disc row of a step tests its CBD with its CBL held equal across groups. When only the location differs, the CBD has to absorb part of the difference. On the TIMSS data, the disc row of M032166 is highly significant, although the LR test of its CBD with the location free (step 2 above) is 0.05:

dif_test(m_mg, cov$country, parm = c("M032166_disc1", "M032166_diff1"))
#> LM tests of DIF between the 2 groups of the covariate (robust information) | unit: parameter (selected parameters) | p_adj: Holm 
#>           unit npar statistic df      p  p_adj
#>  M032166_disc1    1    49.397  1 <1e-04 <1e-04
#>  M032166_diff1    1    78.832  1 <1e-04 <1e-04

A small simulation shows the same (600 persons per group, 200 replications, DSF only in step 2 of item 2; rejection rates at \(\alpha = .05\)):

A location shift is reported as CBD DSF in one of six data sets. The 1-df CBD row is more powerful against a pure CBD difference, but its null hypothesis is not “the CBD is invariant”. Screen with the joint test and decompose with dsf().

• Which information matrix. A score test needs an information matrix, and under misspecification the choice matters (Falk & Monroe, 2018; Guastadisegni et al., 2021). In a simulation of the per-step screen (two groups, 500 replications per cell; misfit = local dependence between two items), the rejection rates of the joint 2-df tests were:

The robust test was the only one that stayed at or below the nominal level, within Monte Carlo error, in every null condition; it is slightly conservative with 300 persons per group. After size adjustment all four have the same power. info = "louis" reproduces the LR test closely (correlation about .99) without refitting and shares its slight excess in small samples; "opg" is liberal.

Why fit the TPPCM rather than the GPCM?

dif_test() can only test parameters the fitted model has. With a GPCM fit discrimination exists only at the item level, so DSF in single steps is missed or shows up as location DIF. Fit the saturated TPPCM to test steps (by = "step" requires it).

tppcmtree(): DIF trees

tppcmtree() uses partykit::mob(): every node is fitted with tam.mml.3pl() and the design of tppcm(), and the parameters are tested for stability along the covariates. parm decides which parameters drive the splits:

• parm = "disc": discriminations only (step discriminations with design = ~ item * step, item discriminations with design = ~ item);
• parm = "diff": step difficulties only;
• in the slope–intercept parameterization (param = "si") the counterparts are parm = "slope" and parm = "int";
• parm = c("disc", "diff") (default): all parameters.

The response matrix goes into a data frame as a matrix column. On the TIMSS data the tree splits by country first, which is the impact of the previous section, and then by gender within Taiwan for the difficulties:

cov$resp <- as.matrix(timss)
tr_diff <- tppcmtree(resp ~ sex + country, data = cov, parm = "diff", minsize = 300)
tr_disc <- tppcmtree(resp ~ sex + country, data = cov, parm = "disc", minsize = 300)
tr_diff
#> TPPCM tree | node design: ~ item * step | splits driven by: diff (irt parameterization)
#> Formula: resp ~ sex + country
#> [1] root
#> |   [2] country in AUS: n = 1050, -logLik = 6193.9
#> |   [3] country in TWN
#> |   |   [4] sex in girl: n = 353, -logLik = 1773.9
#> |   |   [5] sex in boy: n = 366, -logLik = 1923.0
#> 3 terminal node(s), 30 parameters per node; coef(x) gives the node parameters
tr_disc
#> TPPCM tree | node design: ~ item * step | splits driven by: disc (irt parameterization)
#> Formula: resp ~ sex + country
#> [1] root
#> |   [2] country in AUS: n = 1050, -logLik = 6193.9
#> |   [3] country in TWN
#> |   |   [4] sex in girl: n = 353, -logLik = 1773.9
#> |   |   [5] sex in boy: n = 366, -logLik = 1923.0
#> 3 terminal node(s), 30 parameters per node; coef(x) gives the node parameters

The impact problem of the previous section exists in every DIF tree based on a marginal model: a covariate correlated with ability will be split on. Read the first split of a tree on such a covariate as a difference in ability, and compare the item parameters of its daughter nodes with coef() after bringing them to a common scale; or use impact = (latent means and variances estimated in every node), or the multigroup workflow above.

The tree needs a few hundred respondents per node for the saturated TPPCM. In a node where a category of some item is no longer observed the fit stops with a message asking for a larger minsize; design = ~ item or ~ 1 needs fewer respondents.

A detail matters here: when partykit::mob() tests a subset of the parameters it decorrelates all parameters first and then picks the components to test. The estimates of \(a\) and \(b\) are correlated, so instability of \(b\) leaks into the test of \(a\). tppcmtree() uses orthogonal = TRUE by default: the scores of the tested block are residualized on the other block first (efficient scores, \(s_a - s_b \hat I_{bb}^{-1}\hat I_{ba}\)).

Testing the blocks separately matters as well: when all parameters are tested together, strong location DIF splits the sample first and step-discrimination DIF is easily buried, and testing many parameters at once is liberal in moderate samples.

The person scores used by the tree are available with get_parts(fit, "estfun") (or "estfun.si"). With param = "irt" (the default) the discrimination test means “difficulty fixed, discrimination changes”; param = "si" means “intercept fixed, slope changes”. These are different hypotheses.

data("VerbalAggression", package = "psychotools")
tr_va <- tppcmtree(resp ~ gender + anger, data = VerbalAggression,
                   design = ~ 1, parm = "diff", minsize = 50)
partykit::character_split(tr_va$node$split, tr_va$data)$name
#> [1] "gender"

Workflow

1. Fit the PCM, GPCM and step model with TAM (tppcm(dat, design = ~ 1), ~ item, ~ step) and compare them with IRT.compareModels(); the M2, RMSEA, SRMSR, CFI and TLI in the header of irt_pars() show the overall fit of each.
2. From the GPCM: score_test() for the overall picture, item_test() for the items, mi() for the steps. From the PCM, score_test(m_pcm, against = ~ item) and score_test(m_pcm, against = ~ step) test differences between items and between steps.
3. If the EPC pattern is the same across items, test the product form (score_test(m_gpcm, against = ~ item + step), then fit it with xxirt_tppcm()).
4. If only a few items deviate, free only those with tppcm(dat, index = ).
5. Report the final model with irt_pars() (with long = TRUE when needed).
6. With known groups, check DSF: fit the multigroup TPPCM (tam.mml.3pl(E = tppcm(dat), group = g, est.variance = FALSE), so that impact is modelled), screen every step with dif_test(fit, g, by = "step") (Holm-adjusted p_adj), decompose the flagged steps with dsf(fit, items = flagged) (LR tests of the CBD with the location free, and of the IBD), and test configural, metric and scalar models with dsf_design(). Do not read the 1-df disc rows of by = "param" as tests of the CBD alone. For a numeric covariate related to ability use Y =; with many covariates and an unknown grouping use tppcmtree().

Caveats

• Categories: tppcm() allows different numbers of categories per item (every category from 0 to the maximum of each item must be observed); xxirt_tppcm() and the product form need the same number of categories for all items.
• Coding: categories must be 0, 1, …; subtract 1 from 1–5 Likert data, otherwise tppcm() stops with an error.
• Supported fits: unidimensional divide-by-total models (the adjacent-category logits are linear in \(\theta\)) from TAM and xxirt, including binary items, the RSM, multiple groups, weights (pweights) and latent regression (with a warning). Models with guessing parameters, facet models (tam.mml.mfr()), JML and multidimensional models are rejected. With weights, m2() and dif_test() are unavailable; with latent regression, m2() is unavailable.
• Comparing with the PCM: for LR comparisons use the PCM from tppcm(dat, design = ~ 1) or more quadrature nodes. tam.mml() (variance estimated) and the slope-parameterized models cannot be subtracted at the default nodes.