tppcm recipes

This document shows every function of tppcm in the form Problem – Solution – Discussion – See also. Some examples are deliberately wrong to show where the package stops you and what the message looks like. The concepts and the mathematics are in vignette("tppcm-tutorial").

The examples use data sets that ship with other packages:

Data Content
CDM::data.Students Austrian 8th graders (\(N = 2400\)): mathematics self-concept sc1–sc4 and mathematics enjoyment mj1–mj4, scored 0–3, about 2% missing; covariates female, urban
TAM::data.timssAusTwn.scored TIMSS 2011 mathematics, 1050 Australian and 719 Taiwanese students, 11 items of which 4 are scored 0–2 and 7 are dichotomous; covariates country and gender
psychotools::VerbalAggression De Boeck and Wilson (2004): 316 students, 24 items scored 0–2; covariates gender, anger
mirt::Science, TAM::data.gpcm the Euro-Barometer science attitude items (\(N = 392\)), scored 1–4 in Science and 0–3 in data.gpcm 
library(tppcm)   # also attaches TAM
data(data.Students, package = "CDM")
sc <- data.Students[, c("sc1", "sc2", "sc3", "sc4")]

A. Data and estimation

A1. Getting raw questionnaire data into shape

Problem: the questionnaire is scored 1–4, there is an ID column, one variable is a factor, and one item has a category nobody chose. What happens if it goes in as is?

Solution: let tppcm() check it. It reports one problem at a time:

data(Science, package = "mirt")
head(Science, 3)
#>   Comfort Work Future Benefit
#> 1       4    4      3       2
#> 2       3    3      3       3
#> 3       3    2      2       3
tppcm(Science)                                   # scored 1-4
#> Error:
#> ! categories must be coded 0, 1, ..., K-1, but the lowest observed value is 1; recode, e.g. dat - 1
sci <- Science - 1
tppcm(cbind(id = seq_len(nrow(sci)), sci))       # an ID column slipped in
#> Error:
#> ! item(s) id have more than 20 categories; is this an ID or another non-item column?
s_fac <- sci; s_fac$Comfort <- factor(s_fac$Comfort)
tppcm(s_fac)                                     # a factor column
#> Error:
#> ! column(s) Comfort are not numeric; code responses as numbers 0, 1, ..., K-1
s_one <- sci; s_one$Work <- 0
tppcm(s_one)                                     # an item with one category
#> Error:
#> ! item(s) Work have a single category; drop them
s_gap <- sci; s_gap$Work[s_gap$Work == 2] <- 3
tppcm(s_gap)                                     # nobody in category 2 of Work
#> Error:
#> ! unobserved middle categories for item(s) Work; collapse or recode the categories

Collapse the empty category (3 of Work moves down to 2) and check again:

s_gap$Work[s_gap$Work == 3] <- 2
attr(tppcm(s_gap), "index")                      # Work now has 2 steps
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    4    5   NA
#> [3,]    6    7    8
#> [4,]    9   10   11

Discussion: categories must be coded 0, 1, …, and every category from 0 to the maximum of each item must be observed, otherwise the corresponding step cannot be estimated. Items may have different numbers of categories (A4). Missing responses need no treatment; TAM’s marginal maximum likelihood uses what was answered.

See also: ?tppcm, A4.

A2. Fitting the saturated TPPCM

Problem: every step of every item should have its own discrimination.

Solution: tppcm() builds the design array for TAM::tam.mml.3pl(). The latent variance must be fixed:

fit_sat <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE)
fit_sat$gammaslope[1:6]               # the step discriminations
#> sc1_disc1 sc1_disc2 sc1_disc3 sc2_disc1 sc2_disc2 sc2_disc3 
#> 2.9437627 3.7930345 1.2318079 2.2666576 1.6232439 0.6768411

Forgetting est.variance = FALSE leaves the scale unidentified; the other functions warn:

fit_bad <- tam.mml.3pl(sc, E = tppcm(sc), verbose = FALSE)
p_bad <- irt_pars(fit_bad, fit = FALSE)
#> Warning: the latent variance was estimated together with free slopes, so the scale is not
#> identified; refit with tam.mml.3pl(..., est.variance = FALSE)

Discussion: the saturated TPPCM is the same model as TAM’s NRM (tam.mml.2pl(irtmodel = "2PL")) with the same log-likelihood; tppcm() makes the parameters step discriminations directly and allows restrictions (A3).

See also: ?tppcm, vignette("tppcm-tutorial") Section 2.

A3. GPCM, step model, PCM and custom restrictions

Problem: fit something more parsimonious than the saturated model, or free only some items.

Solution: design = for the common models, a formula for the step discriminations read as an item × step layout on the log scale, \(\log a_{il} = \mu + \alpha_i + \gamma_l + (\alpha\gamma)_{il}\):

design = Model \(a_{il}\) Shape of the steps
~ item * step (default) saturated TPPCM \(a_{il}\) item-specific (interaction)
~ item + step Rank-1 (A5; not linear, xxirt_tppcm()) \(\alpha_i\gamma_l\) same shape, item-specific level
~ item GPCM \(\alpha_i\) flat
~ step step model \(\gamma_l\) same shape and level for all items
~ 1 PCM \(a\) flat, the same for all items

index = gives any other equality restriction (same number = same parameter) and takes precedence over design; TAM’s gammaslope.fixed fixes a step:

fit_gpcm <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ item), est.variance = FALSE, verbose = FALSE)
fit_step <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ step), est.variance = FALSE, verbose = FALSE)
fit_pcm  <- tam.mml.3pl(sc, E = tppcm(sc, design = ~ 1),  est.variance = FALSE, verbose = FALSE)
dimnames(tppcm(sc, design = ~ step))[[4]]                  # parameter names
#> [1] "step1" "step2" "step3"
idx <- rbind(1:3, 4:6, rep(7, 3), rep(8, 3))      # only sc1 and sc2 saturated
fit_mix <- tam.mml.3pl(sc, E = tppcm(sc, index = idx), est.variance = FALSE, verbose = FALSE)
fit_fix <- tam.mml.3pl(sc, E = tppcm(sc), est.variance = FALSE, verbose = FALSE,
                       gammaslope.fixed = cbind(2, 3))                  # sc1 step 2 fixed at 3

Discussion: the numbers in index need not be consecutive; they are renumbered. ~ item:step and ~ item + step + item:step are the same as ~ item * step. The interaction term means that items differ in the shape of their step discriminations (the ratios \(a_{il}/a_{il'}\)); ~ step is the rating-scale idea of TAM’s formulaA = ~ item + step applied to the discriminations. Rank-1 is multiplicative, so tppcm() refuses it:

tppcm(sc, design = ~ item + step)
#> Error:
#> ! design = ~ item + step (Rank-1, a_il = alpha_i * gamma_l) is multiplicative and cannot be written as a linear design for TAM; fit it with xxirt_tppcm(dat, design = ~ item + step)

See also: B1 (report), C3 (model comparison).

A4. Items with different numbers of categories

Problem: the test mixes dichotomous and partial-credit items, or the top category of some items was collapsed.

Solution: nothing special. The TIMSS items have 2 or 3 categories:

data(data.timssAusTwn.scored, package = "TAM")
timss <- data.timssAusTwn.scored[, 1:11]
fit_timss <- tam.mml.3pl(timss, E = tppcm(timss), est.variance = FALSE, verbose = FALSE)
irt_pars(fit_timss, se = FALSE)
#> Model: saturated TPPCM, design ~ item * step (tam.mml.3pl) | kernel disc * (theta - diff)
#> Fit: N = 1769, logLik = -10418.73, npar = 30, AIC = 20897.46, BIC = 21061.81
#>      M2(86) = 343.09, p = .000, RMSEA = .041 [.037, .046], SRMSR = .035, CFI = .985, TLI = .982
#> 
#> $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.243        0.537       
#> M032761  3.247 2.659  0.361  0.493
#> M032692  1.442 1.615  1.875 -0.910
#> M032626  1.946       -0.362       
#> M032595  1.939       -0.633       
#> M032673  1.587       -0.666       
#> 
#> $group
#>  group MEAN VAR
#>      1    0   1

Discussion: steps that do not exist are blank in the report and NA in index. Rank-1 (xxirt_tppcm(dat, design = ~ item + step)) and score_test(against = ~ item + step) need the same number of categories in every item.

See also: A1, ?tppcm.

A5. Rank-1 (product form) and the lower bound of the discriminations

Problem: estimate \(a_{il} = \alpha_i \gamma_l\), or keep discriminations from going negative. TAM cannot do either.

Solution: xxirt_tppcm() fits with sirt:

fit_r1 <- xxirt_tppcm(as.matrix(sc), design = ~ item + step)
irt_pars(fit_r1, restricted = TRUE)              # alpha_i and gamma_l
#> Parameters of the fitted model (Rank-1 product form, design ~ item + step), 95% intervals 
#>             est    se lower upper
#> sc1_alpha 3.215 0.209 2.804 3.625
#> sc2_alpha 1.947 0.115 1.722 2.172
#> sc3_alpha 2.451 0.156 2.145 2.758
#> sc4_alpha 2.431 0.158 2.121 2.741
#> gamma2    0.993 0.058 0.879 1.107
#> gamma3    0.419 0.049 0.324 0.514
dim(xxirt_vcov(fit_r1))                          # covariance matrix of the xxirt fit
#> [1] 18 18

An estimate at the bound 0 is marked in the report. The third step of Work in TAM::data.gpcm is such a case: TAM gives \(a = -0.04\) with a meaningless difficulty, xxirt gives \(a = 0\):

data(data.gpcm, package = "TAM")
fit_g3 <- tam.mml.3pl(data.gpcm, E = tppcm(data.gpcm), est.variance = FALSE, verbose = FALSE)
irt_pars(fit_g3, fit = FALSE)$items
#>             disc1     disc2       disc3     diff1     diff2    diff3
#> Comfort 1.4764438 0.6181333  2.24445913 -1.092584 3.8290106 2.200135
#> Work    0.7811138 0.3291185 -0.03993933 -2.040273 2.4080559       NA
#> Benefit 1.2492408 1.4054680  3.77027046 -1.243889 0.7109339 1.682515
irt_pars(xxirt_tppcm(as.matrix(data.gpcm)), fit = FALSE)
#> Model: saturated TPPCM, design ~ item * step | kernel disc * (theta - diff)
#>
#> $items
#>         disc1 disc2 disc3  diff1 diff2 diff3
#> Comfort 1.486 0.619 2.252 -1.089 3.824 2.200
#> Work    0.781 0.319 0.000 -2.041 2.473
#> Benefit 1.232 1.397 3.845 -1.252 0.712 1.678
#>
#> $se
#>         disc1 disc2 disc3 diff1 diff2 diff3
#> Comfort 0.560 0.358 1.244 0.236 1.898 0.386
#> Work    0.254 0.181       0.535 1.346
#> Benefit 0.597 0.782 4.838 0.308 0.230 0.259
#>
#> On the bound 0 (no SE; tests involving them are not valid): Work_disc3
#>
#> $group
#>  group MEAN VAR
#>      1    0   1

Discussion: a parameter on the bound has no normal-theory SE, and LR and score tests involving it no longer follow a \(\chi^2\) distribution. A discrimination near 0 makes the step difficulty \(b = d/a\) unstable (here \(-27\) with SE 205); report such items in the slope–intercept form (irt_pars(IRTpars = FALSE)), whose intercepts are well behaved. A negative discrimination usually means the scoring direction of that step is doubtful. xxirt_tppcm() uses TAM’s quadrature nodes, so its log-likelihood is comparable with TAM’s (C3).

See also: ?xxirt_tppcm, ?xxirt_vcov, C2.

B. Reports and components

B1. One parameter table

Problem: a compact table like mirt’s coef(, IRTpars = TRUE, simplify = TRUE).

Solution:

irt_pars(fit_sat)
#> 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.95, 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

Other forms:

irt_pars(fit_sat, IRTpars = FALSE, fit = FALSE, se = FALSE)   # slope-intercept: slope * theta + int
#> Model: saturated TPPCM, design ~ item * step (tam.mml.3pl) | kernel slope * theta + int
#> 
#> $items
#>     slope1 slope2 slope3  int1   int2   int3
#> sc1  2.944  3.793  1.232 1.689 -2.943 -3.465
#> sc2  2.267  1.623  0.677 0.604 -1.596 -1.945
#> sc3  2.577  2.261  1.022 0.150 -2.054 -1.841
#> sc4  2.060  2.635  1.171 2.095 -0.940 -2.128
#> 
#> $group
#>  group MEAN VAR
#>      1    0   1
head(irt_pars(fit_sat, long = TRUE, level = 0.90))            # long format with intervals
#> 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
irt_pars(fit_gpcm, restricted = TRUE)                         # parameters of the fitted model itself
#> Parameters of the fitted model (GPCM, design ~ item (tam.mml.3pl)), 95% intervals 
#>            est    se lower upper
#> sc1_disc 2.941 0.168 2.611 3.271
#> sc2_disc 1.767 0.093 1.584 1.950
#> sc3_disc 2.160 0.129 1.908 2.412
#> sc4_disc 2.053 0.108 1.841 2.264

For a paper, the pieces are plain data frames:

p <- irt_pars(fit_sat, fit = FALSE)
knitr::kable(cbind(p$items[, 1:3], p$se[, 1:3]), digits = 2,
             col.names = c(paste0("a", 1:3), paste0("SE(a", 1:3, ")")))
a1 a2 a3 SE(a1) SE(a2) SE(a3)
sc1 2.94 3.79 1.23 0.22 0.35 0.28
sc2 2.27 1.62 0.68 0.15 0.14 0.17
sc3 2.58 2.26 1.02 0.19 0.22 0.20
sc4 2.06 2.63 1.17 0.18 0.19 0.16

Discussion: the Fit: header has the log-likelihood, number of parameters, AIC and BIC (as TAM’s IRT.IC()) and the M2 family (C4). diff is Masters’ step difficulty, equal to TAM’s beta + tau and mirt’s b with IRTpars = TRUE.

See also: ?irt_pars, B2, C4.

B2. Standard errors and confidence intervals

Problem: standard errors and confidence intervals for the step discriminations.

Solution:

st <- irt_pars(fit_sat, long = TRUE, level = 0.90)
head(st[st$par == "disc", ])
#> 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
#> sc2_disc1  sc2    1 disc 2.267 0.152 2.017 2.517
#> sc2_disc2  sc2    2 disc 1.623 0.140 1.392 1.854
#> sc2_disc3  sc2    3 disc 0.677 0.174 0.391 0.963

Discussion: irt_pars() uses the observed information of the whole model (Louis’ formula), including the covariance between slopes and thresholds. level = sets the confidence level.

See also: vignette("tppcm-tutorial") Section 3.

B3. Getting at the components

Problem: further computations need the information matrix, the person scores, the Jacobian.

Solution: get_parts(), used like lme4::getME() or lavaan::lavInspect():

get_parts(fit_gpcm, "est")
#> sc1_disc sc2_disc sc3_disc sc4_disc 
#> 2.941187 1.766827 2.159744 2.052681
get_parts(fit_gpcm, c("loglik", "npar", "nobs"))
#> $loglik
#> [1] -9479.182
#> 
#> $npar
#> [1] 16
#> 
#> $nobs
#> [1] 2400
dim(get_parts(fit_gpcm, "information"))
#> [1] 24 24
head(get_parts(fit_gpcm, "estfun")[, 1:4], 3)
#>         sc1_disc   sc2_disc   sc3_disc    sc4_disc
#> [1,]  0.02977269  0.1433597 0.10612541 -0.01434813
#> [2,] -0.05007352 -0.2427045 0.11333258  0.13069880
#> [3,]  0.11521202  0.1773845 0.06919201  0.14183058

When several things are to be done with one model, take the whole bundle once and pass it on:

sp <- get_parts(fit_gpcm)
sp
#> Components of a fitted model: GPCM, design ~ item (tam.mml.3pl) | 4 items, 4 categories, 2400 persons
#> get_parts(x, what = ) with what in:
#>   disc, diff, int, est, se, vcov, par, partable, vcov.sat, se.sat, information,
#>   information.opg, gradient, scores, estfun, estfun.si, jacobian, curvature,
#>   posterior, nodes, prior, probs, data, nobs, loglik, npar, model, formula,
#>   groups, design
logLik(sp); nobs(sp)
#> 'log Lik.' -9479.182 (df=16)
#> [1] 2400
get_parts(fit_gpcm, "hessian")      # no such component
#> Error:
#> ! unknown component 'hessian'; see ?get_parts

Discussion: printing the bundle lists all available what. The internal order of the saturated parameters is \(a_{i1..L}, d_{i1..L}\) per item with \(d = a \cdot \text{diff}\).

See also: ?get_parts.

B4. Information functions

Problem: where on the trait does each item, each step, give most information?

Solution:

inf <- step_info(fit_sat)
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.622         0.60
op <- par(mfrow = c(1, 2)); plot(inf); plot(inf, item = "sc1"); par(op)

Discussion: the item information \(I_i(\theta) = \mathrm{Var}(\sum_{l \le X} a_{il} \mid \theta)\) is split into step contributions that add up to it. Item and test information agree with TAM’s IRT.informationCurves(); the step contributions are new.

See also: ?step_info.

C. Tests

C1. Which items have unequal step discriminations

Problem: under the saturated model, test \(a_{i1} = a_{i2} = a_{i3}\) item by item.

Solution:

wald_test(fit_sat)
#> 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.630  2 <1e-04 <1e-04
#>   sc2    49.249  2 <1e-04 <1e-04
#>   sc3    34.260  2 <1e-04 <1e-04
#>   sc4    36.874  2 <1e-04 <1e-04
wald_test(fit_gpcm)          # in a GPCM the steps are equal by construction
#> no item has free step discriminations; fit the saturated TPPCM, e.g. E = tppcm(dat)
#> Wald tests of equal step discriminations within items (H0: the item follows the GPCM) | p_adj: Holm 
#> (nothing to test)

Discussion: the Wald test needs the saturated model (or at least the item in question free). When only the GPCM has been fitted, use the score tests of C2, which do not need the larger model.

See also: C2.

C2. Is the GPCM enough

Problem: only the GPCM was fitted. Should the step discriminations be freed, and where?

Solution:

score_test(fit_gpcm)                    # GPCM against the saturated TPPCM
#> Score test (Louis information) of the fitted GPCM, design ~ item (tam.mml.3pl) against the saturated TPPCM (~ item * step)
#>   S = 212.9262, df = 8, p = < 2.2e-16
score_test(fit_gpcm, against = ~ item + step)           # GPCM against the product form
#> Score test (Louis information) of the fitted GPCM, design ~ item (tam.mml.3pl) against the Rank-1 product form (~ item + step)
#>   S = 172.6925, df = 2, p = < 2.2e-16
item_test(fit_gpcm)                     # which items
#> Score tests for freeing the step discriminations of each item in the fitted GPCM, design ~ item (tam.mml.3pl) | p_adj: Holm 
#>  item statistic df      p  p_adj
#>   sc1    33.159  2 <1e-04 <1e-04
#>   sc2    40.406  2 <1e-04 <1e-04
#>   sc3    19.816  2 <1e-04 <1e-04
#>   sc4    27.033  2 <1e-04 <1e-04
head(mi(fit_gpcm), 5)                   # which steps (modification indices and EPCs)
#> Modification indices (1-df score tests) for freeing single step discriminations of the fitted GPCM, design ~ item (tam.mml.3pl) | p_adj: Holm 
#>        par     MI df      p    p_adj    EPC
#>  sc2_disc3 31.493  1 <1e-04  < 1e-04 -1.658
#>  sc1_disc3 28.572  1 <1e-04  < 1e-04 -3.464
#>  sc2_disc1 22.559  1 <1e-04  < 1e-04  0.840
#>  sc4_disc3 21.408  1 <1e-04  < 1e-04 -1.225
#>  sc3_disc3 18.620  1 <1e-04 0.000128 -1.558
score_test(fit_sat)                     # already saturated
#> Error:
#> ! the fitted model (saturated TPPCM, design ~ item * step (tam.mml.3pl)) already contains the model ~ item * step; nothing to test
score_test(fit_timss, against = ~ item + step)          # TIMSS items have 1 or 2 steps
#> Error:
#> ! Rank-1 (against = ~ item + step) needs the same number of categories in every item

Discussion: a score test uses only the estimates of the restricted model. against is a design formula of a larger model (~ item * step, the default; ~ item + step; ~ item; ~ step), or D = gives arbitrary directions. Against ~ item + step the GPCM is tested for a step effect, Rank-1 against ~ item * step for an item × step interaction. The sign of an EPC says whether the step should be steeper (positive) or flatter (negative) than the common GPCM slope; here every third step is flatter, which Rank-1 captures with one \(\gamma_3\). The model label GPCM, design ~ item (tam.mml.3pl) means a GPCM fitted through tppcm(dat, design = ~ item); a fit from tam.mml.2pl(irtmodel = "GPCM") is labeled GPCM, design ~ item and gives the same results; get_parts(fit, "formula") returns the design.

See also: C1, C3, vignette("tppcm-tutorial") Section 4.

C3. Likelihood-ratio tests of nested models

Problem: compare PCM, GPCM, step model, Rank-1 and the saturated model.

Solution: fit all models in the slope parameterization (the PCM with tppcm(dat, design = ~ 1)) and hand them to TAM:

lr <- IRT.compareModels(fit_pcm, fit_gpcm, fit_step, fit_sat)$LRtest
lr
#>     Model1   Model2      Chi2 df            p
#> 1  fit_pcm fit_gpcm  45.01679  3 9.176983e-10
#> 2  fit_pcm fit_step  97.42346  2 0.000000e+00
#> 3  fit_pcm  fit_sat 163.98336 11 0.000000e+00
#> 4 fit_step fit_gpcm -52.40667  1 1.000000e+00
#> 5 fit_gpcm  fit_sat 118.96657  8 0.000000e+00
#> 6 fit_step  fit_sat  66.55990  9 7.167356e-11

The \(\chi^2\) in row 4 is negative: TAM treats the model with more parameters as the larger one, but the step model (3 parameters) and the GPCM (4) do not contain each other, so that row is meaningless. Keep the nested pairs only:

lr[-4, ]
#>     Model1   Model2      Chi2 df            p
#> 1  fit_pcm fit_gpcm  45.01679  3 9.176983e-10
#> 2  fit_pcm fit_step  97.42346  2 0.000000e+00
#> 3  fit_pcm  fit_sat 163.98336 11 0.000000e+00
#> 5 fit_gpcm  fit_sat 118.96657  8 0.000000e+00
#> 6 fit_step  fit_sat  66.55990  9 7.167356e-11

Mixing TAM and xxirt objects (the product form):

suppressWarnings(IRT.compareModels(fit_gpcm, fit_r1, fit_sat))$LRtest
#>     Model1  Model2      Chi2 df            p
#> 1 fit_gpcm  fit_r1  91.90638  2 0.0000000000
#> 2 fit_gpcm fit_sat 118.96657  8 0.0000000000
#> 3   fit_r1 fit_sat  27.06019  6 0.0001410955

Discussion:

  • IRT.compareModels() compares all pairs and does not check nesting. The nesting is: ~ 1 ⊂ ~ item, ~ 1 ⊂ ~ step ⊂ ~ item + step ⊂ ~ item * step; ~ item ⊂ ~ item + step; ~ item and ~ step do not contain each other.
  • Do not compare TAM’s own PCM from tam.mml() (variance estimated) with the slope-parameterized models by LR: the two are nested mathematically, but their numerical integration errors differ at TAM’s default 21 nodes; the LR can even be negative.
  • With mixed TAM and xxirt objects, look only at LRtest; TAM warns when computing the GHP column for the xxirt object, and that column is wrong.

See also: A3, A5.

C4. Overall fit

Problem: M2, RMSEA, SRMSR, CFI and TLI in addition to \(\chi^2\) comparisons.

Solution:

m2(fit_gpcm)                 # full M2
#> M2 = 845.85, 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(fit_sat)
#> M2 = 616.95, 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(fit_gpcm, "C2")           # C2, margins reduced to moments
#> 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(fit_timss)                # missing data and unequal categories are fine
#> M2 = 343.09, df = 86, p = <2e-16
#> RMSEA = 0.041 [0.037, 0.046], SRMSR = 0.035, CFI = 0.985, TLI = 0.982
w <- runif(nrow(sc), 0.5, 2)
fit_w <- tam.mml.2pl(sc, irtmodel = "GPCM", pweights = w, verbose = FALSE)
m2(fit_w)                    # not available with weights
#> Error:
#> ! m2() does not support person weights

Discussion: C2 and M2* collapse the categories into moments of the item scores, so differences between steps are averaged away; to check step discriminations use the full M2 or the score tests of C2. With missing data the margins use all available pairs, which assumes MCAR.

See also: ?m2, B1.

D. Differential item and step functioning

D1. Known groups: which step has DIF

Problem: do boys and girls differ in the discrimination or difficulty of some step?

Solution: when the covariate is unrelated to ability, fit the saturated model once to the whole sample, then dif_test():

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(fit_timss)                 # reused below
d_disc <- dif_test(sp_timss, cov$sex, parm = "disc")   # every step discrimination
#> 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)
d_disc[d_disc$p < 0.10, ]
#> Score-based DIF tests along the covariate (LM functional) | unit: parameter (disc) | p_adj: Holm 
#>              unit npar statistic df      p p_adj
#> 5  M032760A_disc1    1     2.929  1 0.0870     1
#> 10  M032761_disc2    1     3.018  1 0.0823     1
dif_test(sp_timss, cov$sex, 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.565  1 0.21099 1.0000
#>   M032721    1     0.154  1 0.69499 1.0000
#>   M032757    2     1.908  2 0.38526 1.0000
#>  M032760A    2     9.637  2 0.00808 0.0889
#>  M032760B    1     0.282  1 0.59566 1.0000
#>  M032760C    1     0.130  1 0.71855 1.0000
#>   M032761    2     3.914  2 0.14130 1.0000
#>   M032692    2     1.298  2 0.52244 1.0000
#>   M032626    1     0.992  1 0.31918 1.0000
#>   M032595    1     0.901  1 0.34261 1.0000
#>   M032673    1     0.158  1 0.69139 1.0000

Discussion: every unit is residualized on all other parameters first (efficient score), so difficulty DIF is not mistaken for discrimination DIF and vice versa. Nothing survives the Holm correction here. With a GPCM fit, discrimination exists only at the item level, and step DIF is either invisible or misreported as difficulty DIF. The message about impact is a reminder that a single-group fit is valid only if the groups do not differ in ability; otherwise fit the groups (D2) and use the step-level workflow of D5.

See also: ?dif_test, D2, D3, D5.

D3. Numeric and ordered covariates

Problem: the covariate is trait anger (numeric) or an ordered grouping.

Solution: dif_test() picks the functional from the covariate type:

data("VerbalAggression", package = "psychotools")
va <- VerbalAggression
fit_va <- tam.mml.3pl(va$resp, E = tppcm(va$resp, design = ~ item), est.variance = FALSE, verbose = FALSE)
sp_va <- get_parts(fit_va)
d_anger <- dif_test(sp_va, va$anger, parm = "diff", by = "item")        # numeric: Cramer-von Mises
head(d_anger[order(d_anger$p), ], 5)
#> Score-based DIF tests along the covariate (CvM functional) | unit: item (diff) | p_adj: Holm 
#>           unit npar statistic df       p p_adj
#> 6    S1DoShout    2     1.500 NA 0.00500 0.120
#> 9  S2WantScold    2     1.163 NA 0.00635 0.146
#> 13 S3WantCurse    2     1.150 NA 0.00683 0.150
#> 11 S2WantShout    2     1.069 NA 0.00976 0.205
#> 2    S1DoCurse    2     0.967 NA 0.01724 0.345
anger3 <- cut(va$anger, quantile(va$anger, 0:3 / 3), include.lowest = TRUE, ordered_result = TRUE)
dif_test(sp_va, anger3, parm = "disc", by = "block")                     # ordered: weighted max LM
#> 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 (maxLMo functional) | unit: all disc jointly | p_adj: Holm 
#>  unit npar statistic df     p p_adj
#>   all   24     1.992 NA 0.878 0.878
dif_test(sp_va, va$anger, parm = "disc", by = "block", functional = "maxLM")
#> Score-based DIF tests along the covariate (maxLM functional) | unit: all disc jointly | p_adj: Holm 
#>  unit npar statistic df     p p_adj
#>   all   24    47.832 NA 0.068 0.068
dif_test(sp_va, replace(va$anger, 1, NA))                                # covariate with NA
#> Error:
#> ! 'covariate' has missing values; drop those persons and refit
dif_test(sp_va, va$anger[-1])                                            # wrong length
#> Error:
#> ! 'covariate' must have one value per person of the fit (316), not 315

Discussion: functional = can be "LM", "maxLMo", "DM", "CvM" or "maxLM". maxLM also estimates the location of the change point. The covariate may not have missing values; drop those persons and refit. With 24 items and 316 students the GPCM is the sensible node model here; the saturated TPPCM would have 96 parameters.

See also: D1, D4.

D4. Unknown grouping: a DIF tree

Problem: several covariates, and the data should find the grouping.

Solution: tppcmtree(), looking at difficulties and discriminations separately. With the verbal aggression data:

tr_diff <- tppcmtree(resp ~ gender + anger, data = va, design = ~ 1, parm = "diff", minsize = 50)
partykit::width(tr_diff)
#> [1] 2
partykit::character_split(tr_diff$node$split, tr_diff$data)$name
#> [1] "gender"
tr_disc <- tppcmtree(resp ~ gender + anger, data = va, design = ~ item, parm = "disc", minsize = 50)
partykit::width(tr_disc)
#> [1] 1
tppcmtree(resp ~ gender + anger, data = va, design = ~ item, parm = "diff", minsize = 50)
#> Error:
#> ! not all categories observed for item(s) S3DoShout, S4DoShout in a node; increase minsize (e.g. minsize = 100) or use design = ~ 1

Discussion: the difficulty tree splits on gender; the discrimination tree does not split. The failed call shows the limiting factor: in some node a category of an item is no longer observed, and the node model cannot be estimated. A larger minsize or a less demanding node model (design = ~ 1) helps. coef(tree) gives the parameters per node.

See also: ?tppcmtree, D1, D2.

D5. Which step functions differently between known groups (DSF)

Problem: the groups are known (here the countries of D2). Which category boundaries differ, and is it their discrimination (CBD) or their location (IBD/CBL)?

Solution: on the multigroup fit, screen every step with the joint 2-df score test, decompose flagged steps with likelihood ratio tests, and compare configural, metric and scalar models:

st <- dif_test(fit_mg, cov$country, by = "step")       # 1. screen (Holm-adjusted p_adj)
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"
dsf(fit_mg, items = "M032692", verbose = FALSE)        # 2. decompose one item (refits)
#> 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
#>  M032692    1 cbd  0.638  1.244   5.455  1 0.01952 0.0195
#>  M032692    1 ibd -3.733 -2.269  80.896  1 < 1e-04 <1e-04
#>  M032692    2 cbd  0.648  1.382   6.816  1 0.00903 0.0181
#>  M032692    2 ibd  0.114  2.019 129.658  1 < 1e-04 <1e-04
d_part <- dsf_design(timss, cov$country, group.partial = list(ibd = "M032692"))
d_part                                                 # 3. a partial-invariance model
#> Multigroup TPPCM design: 11 items x 2 groups (reference group: AUS)
#>   group.equal: cbd, ibd -> scalar (CBD and IBD equal); within groups: design ~ item * step
#>   group.partial cbd: none
#>   group.partial ibd: M032692:1, M032692:2
#>   15 CBD and 17 IBD parameters; free latent means: groups 2; variances: groups 2
#> Fit with TAM::tam.mml.3pl(d$resp, A = d$A, E = d$E, group = <group>,
#>         variance.fixed = d$variance.fixed, beta.fixed = d$beta.fixed)
m_part <- tam.mml.3pl(d_part$resp, A = d_part$A, E = d_part$E, group = cov$country,
                      variance.fixed = d_part$variance.fixed,
                      beta.fixed = d_part$beta.fixed, verbose = FALSE)
IRT.compareModels(fit_mg, m_part)$LRtest
#>   Model1 Model2     Chi2 df p
#> 1 fit_mg m_part 131.1685  2 0
dif_test(fit_mg, cov$country, parm = "disc", by = "step")   # a step is tested jointly
#> Error:
#> ! by = "step" tests the discrimination and location of a step jointly; use parm = c("disc", "diff")
dsf(fit_timss)                                              # needs a multigroup fit
#> Error:
#> ! fit has a single group; refit with tam.mml.3pl(..., group = g, est.variance = FALSE)

Discussion:

  • dif_test(by = "step") tests the CBD and location of every step jointly (Glas, 1998; info = "robust", the generalized score test of Boos, 1992, is the default). Its null hypothesis does not depend on the parameterization. Do not screen with the 1-df disc rows of by = "param": they hold the location of the step equal across groups, so location DSF leaks into them (vignette("tppcm-tutorial"), Section 6).
  • dsf() releases one parameter at a time against the scalar baseline (group.equal = c("cbd", "ibd"), lavaan’s loadings and intercepts) with all other steps as anchors. The CBD row frees the location of the step as well, so it does not depend on the origin of the scale. Both steps of M032692 differ mainly in their IBD (LR 81 and 130); their CBDs differ less (LR 5.5 and 6.8). Each test refits the model, so pass only the flagged steps: dsf(fit_mg, items = flagged).
  • dsf_design() gives the configural (group.equal = NULL), metric ("cbd") and scalar models and partial models (group.partial), with the identification of the latent means and variances; fit them with tam.mml.3pl() as above. The scalar model has the log-likelihood of fit_mg.

See also: ?dif_test, ?dsf, ?dsf_design, D2, vignette("tppcm-tutorial") Section 6.

E. Simulation, unsupported models, weights and groups

E1. Simulating data and estimating power

Problem: simulate TPPCM data to estimate the power of a test.

Solution: sim_tppcm(); trailing NA in a row means fewer steps for that item:

disc <- rbind(c(1, 1.5, 0.7), c(1, 1.5, 0.7), c(1, 1, NA))
diff <- rbind(c(-1, 0, 1),    c(-1, 0, 1),    c(-0.5, 0.5, NA))
head(sim_tppcm(5, disc, diff))
#>   I1 I2 I3
#> 1  2  1  1
#> 2  3  0  2
#> 3  2  1  2
#> 4  2  3  0
#> 5  3  0  1

A small power simulation (score test GPCM against the saturated model, \(N = 300\)):

set.seed(3)
A <- outer(c(1, 1, 1.2, 0.8, 1), c(1, 1.5, 0.7)); B <- matrix(c(-1, 0, 1), 5, 3, byrow = TRUE)
p <- replicate(10, {
  d <- sim_tppcm(300, A, B)
  score_test(tam.mml.2pl(d, irtmodel = "GPCM", verbose = FALSE))$p.value
})
mean(p < 0.05)
#> [1] 0.7
sim_tppcm(10, disc, diff[, 1:2])                 # shapes differ
#> Error:
#> ! disc and diff must have the same dimensions (items x steps)

Discussion: use many more replications in practice. theta = sets the person parameters, e.g. for a non-normal latent distribution.

See also: ?sim_tppcm.

E2. Models the package does not handle

Problem: which fitted objects are refused?

Solution: the package refuses and says why:

set.seed(4)
th <- rnorm(1000)
bin <- sapply(1:6, function(i) rbinom(1000, 1, 0.2 + 0.8 * plogis(th - (i - 3.5) / 2)))
irt_pars(tam.mml.3pl(bin, guess = rep(0.2, 6), est.variance = FALSE, verbose = FALSE))   # guessing
#> Error:
#> ! the fitted item response functions are not of the TPPCM (divide-by-total) form (e.g. guessing parameters); this model is not supported
irt_pars(tam.mml.2pl(sc, irtmodel = "GPCM", Q = cbind(c(1, 1, 0, 0), c(0, 0, 1, 1)),
                     verbose = FALSE))                                                   # two dimensions
#> Error:
#> ! only unidimensional models are supported
irt_pars(sc)                                                                             # not a fit
#> Error:
#> ! expected a fitted model from TAM (tam.mml(), tam.mml.2pl(), tam.mml.3pl()) or xxirt_tppcm(), not an object of class 'data.frame'

Discussion: the package handles unidimensional divide-by-total models (adjacent-category logits linear in \(\theta\)), including binary items, the RSM, multiple groups, weights and latent regression. The form is checked from the response probabilities, so other models that violate it are stopped as well. Facet models (tam.mml.mfr()) and JML are not supported.

See also: E3, ?get_parts.

E3. Weights, multiple groups and latent regression

Problem: which functions work with sampling weights, several groups or a latent regression?

Solution:

irt_pars(fit_w, fit = FALSE)$se[1:2, ]                                            # weights enter the SEs
#>          disc1      disc2      disc3      diff1     diff2      diff3
#> sc1 0.15182098 0.15182098 0.15182098 0.02941799 0.0318009 0.05891535
#> sc2 0.08690086 0.08690086 0.08690086 0.03378903 0.0411399 0.06523290
fit_w1 <- tam.mml.2pl(sc, irtmodel = "GPCM", pweights = w / mean(w), verbose = FALSE)
irt_pars(fit_w1, fit = FALSE)$se[1:2, ]                                           # after normalizing
#>          disc1      disc2      disc3      diff1      diff2      diff3
#> sc1 0.16973002 0.16973002 0.16973002 0.03288818 0.03555218 0.06586509
#> sc2 0.09715182 0.09715182 0.09715182 0.03777483 0.04599283 0.07292788
irt_pars(fit_mg, fit = FALSE)$group                                               # group means and variances
#>   group     MEAN      VAR    SE_VAR    SE_MEAN
#> 1   AUS 0.000000 1.000000        NA         NA
#> 2   TWN 1.540195 1.439849 0.1573329 0.07598995
fit_lr <- tam.mml.3pl(sc, E = tppcm(sc), Y = data.Students[, c("female", "urban")],
                      est.variance = FALSE, verbose = FALSE)
p_lr <- irt_pars(fit_lr, fit = FALSE)                                             # warns
#> Warning: latent regression: the regression coefficients are treated as fixed in the
#> standard errors and tests

Discussion: weights are treated as frequency weights, so their sum is the effective sample size; with a mean weight of about 1.25 the SEs above are smaller than without weights. Normalize sampling weights to sum to the number of persons, e.g. w / mean(w); even then the SEs are model-based, not design-based (sandwich) SEs.

Situation Available Not available
weights parameters, SEs, Wald, score tests, MI m2(), dif_test()
multiple groups parameters, SEs (the group means and variances are part of the parameter vector), all tests, dif_test(), dsf(), dsf_design() m2()
latent regression parameters, SEs (regression coefficients fixed), tests, dif_test() m2()

See also: ?get_parts, D2, D5.