A large fit is expensive, and the questions asked of it afterwards are many: does this item behave differently for girls and boys, does that one drift with age, does school type predict ability once the rest is accounted for? Refitting for each question does not scale. Score-based tests answer them from one fit: each person’s contribution to the score of the marginal likelihood says which way, and how strongly, that person’s responses pull on every parameter, and a covariate that orders those pulls is evidence that the parameter depends on it.
scoreContributions() computes the contributions, and
scoreTest() turns them into tests: parametric score (LM)
tests with a one-step estimate of the effect, and nonparametric
fluctuation tests that need no functional form.
The covariate z raises ability (impact), and makes items
1 to 4 easier for persons with high z (uniform DIF).
set.seed(1)
N <- 1500; J <- 20
a <- rlnorm(J, 0, .2); b <- rnorm(J)
z <- runif(N)
theta <- rnorm(N) + (z > .5)
B <- matrix(b, N, J, byrow = TRUE)
B[, 1:4] <- B[, 1:4] - 1.2 * (z - .5)
P <- plogis(matrix(a, N, J, byrow = TRUE) * theta - B)
dat <- data.frame(id = rep(1:N, J), Item = rep(1:J, each = N), Scale = 1,
score = rbinom(N * J, 1, P))
dat$z <- z[dat$id]Score contributions have mean zero, for every person and
conditionally on the person’s covariates, when the model is right. A
covariate that is related to ability is therefore part of the model:
left out, its effect on ability shows up as DIF in every item. So it
goes into personPreds, in a form flexible enough to carry
the relation – here a spline, since the impact is a step that a linear
term would not capture.
Bz <- splines::ns(z, df = 3)
dat[c("z1", "z2", "z3")] <- Bz[dat$id, ]
fit <- fitIRT(dat, pl = 2, personPreds = c("z1", "z2", "z3"),
laplaceRefine = "agq", cores = 2, verbose = 0)laplaceRefine = "agq" finishes the fit on the quadrature
marginal, the same one the score contributions are computed from, so
their sums vanish at the estimates. It matters most when each person
answers few items.
Putting the covariate in the fit is the cleanest route, but not a
requirement: given the fit, scoreTest() adds a covariate
the ability model lacks after the fact (see below).
sc <- scoreContributions(fit)
cov <- data.frame(id = rownames(fit$pars$Ability))
cov$z <- z[as.integer(cov$id)]
dif <- scoreTest(sc, cov, parameters = "B", form = "linear")
dif
#> Score tests of B: 20 tests, rule gh9^1
#> z: 9 of 20 at p < .05 (7 after FDR); smallest p 4.59e-13
#> covariate item parameters form n LM df p
#> <char> <char> <char> <char> <int> <num> <int> <num>
#> 1: z 2 B linear 1500 52.375007 1 4.585224e-13
#> 2: z 3 B linear 1500 21.727548 1 3.142461e-06
#> 3: z 1 B linear 1500 13.452135 1 2.447275e-04
#> 4: z 4 B linear 1500 9.848323 1 1.699867e-03
#> 5: z 8 B linear 1500 9.019691 1 2.670866e-03
#> 6: z 11 B linear 1500 8.274742 1 4.020037e-03
#> 7: z 14 B linear 1500 8.027334 1 4.607660e-03
#> 8: z 16 B linear 1500 4.256911 1 3.909093e-02
#> 9: z 17 B linear 1500 3.932277 1 4.736777e-02
#> 10: z 19 B linear 1500 1.846425 1 1.741997e-01
#> info_ratio ability_adjusted p_fdr p_holm
#> <num> <lgcl> <num> <num>
#> 1: 0.9530758 FALSE 9.170447e-12 9.170447e-12
#> 2: 0.9659637 FALSE 3.142461e-05 5.970677e-05
#> 3: 0.9720477 FALSE 1.631516e-03 4.405095e-03
#> 4: 1.1159979 FALSE 8.499334e-03 2.889774e-02
#> 5: 0.8930107 FALSE 1.068346e-02 4.273385e-02
#> 6: 1.0674953 FALSE 1.316474e-02 6.030055e-02
#> 7: 1.0232933 FALSE 1.316474e-02 6.450725e-02
#> 8: 1.0327843 FALSE 9.772733e-02 5.081821e-01
#> 9: 1.0420700 FALSE 1.052617e-01 5.684132e-01
#> 10: 1.0128440 FALSE 3.483994e-01 1.000000e+00The coefficients are one-step estimates of each item’s
linear DIF slope on the intercept scale (positive is harder), with
standard errors: what a refit with BitemPreds = "z" and
itemSpecificBetas = TRUE would estimate.
head(dif$coefficients[order(item)], 6)
#> covariate item parameter term estimate se z
#> <char> <char> <char> <char> <num> <num> <num>
#> 1: z 1 B linear -0.8255522 0.2250862 -3.6677153
#> 2: z 10 B linear 0.2239522 0.2146193 1.0434857
#> 3: z 11 B linear 0.7140795 0.2482386 2.8765852
#> 4: z 12 B linear 0.1881709 0.2290706 0.8214537
#> 5: z 13 B linear 0.2486350 0.2143140 1.1601438
#> 6: z 14 B linear 0.5973658 0.2108408 2.8332551Uniform DIF that every item shares cannot be told apart from an effect on ability, and the ability model absorbs it; DIF is measured relative to it.
The first screen flags the four DIF items, and also a few clean ones
with positive slopes. That is contamination, not noise: four items
easier for high z pull those persons’ ability posteriors
up, and the ability coefficient on z with them, so the
clean items look harder for them. purify = TRUE anchors on
the unflagged items – every person’s posterior and the ability
coefficients are recomputed from those alone, the item parameters held –
and retests every item against that anchor, until the flagged set stops
changing. It needs the fit rather than its scores, but no refit.
pur <- scoreTest(fit, cov, parameters = "B", purify = TRUE, adjust = "none")
pur$tests[order(p)][1:8, .(item, p, p_fdr, anchor, stage)]
#> item p p_fdr anchor stage
#> <char> <num> <num> <lgcl> <int>
#> 1: 2 1.698414e-15 3.396828e-14 FALSE 4
#> 2: 3 9.273267e-08 9.273267e-07 FALSE 4
#> 3: 1 5.504112e-06 3.669408e-05 FALSE 4
#> 4: 4 2.505248e-04 1.252624e-03 FALSE 4
#> 5: 8 2.306478e-02 9.225912e-02 TRUE 4
#> 6: 14 4.522655e-02 1.507552e-01 TRUE 4
#> 7: 15 1.066258e-01 3.046451e-01 TRUE 4
#> 8: 11 2.032501e-01 5.051868e-01 TRUE 4(adjust = "none" because z is already in
the ability model, as z1 to z3; see the next
section.)
A large fit is made once, and the covariates to screen keep coming.
Given the fit, scoreTest() adds a covariate missing from
the ability model after the fact: a natural spline of it (or indicators,
for a factor) enters the ability regression, its coefficients and the
ability SD are estimated with the item parameters held, and the scores
are recomputed under that model before testing. Here the fit leaves
z out entirely.
fit0 <- fitIRT(dat, pl = 2, laplaceRefine = "agq", cores = 2, verbose = 0)
plain <- suppressWarnings(scoreTest(fit0, cov, parameters = "B", adjust = "none"))
adjusted <- scoreTest(fit0, cov, parameters = "B")
flag <- function(t) sort(as.integer(t$tests$item[t$tests$p_fdr < .05]))
list(unadjusted = flag(plain), adjusted = flag(adjusted), refit = flag(dif))
#> $unadjusted
#> [1] 1 2 3 4
#>
#> $adjusted
#> [1] 1 2 3 4 8 11 14
#>
#> $refit
#> [1] 1 2 3 4 8 11 14With the adjustment the screen is the one on the fit that had
z in its ability model, item for item – including items 8,
11 and 14, the contamination the previous section removes. Without it
the test is not valid: unmodelled impact makes every item look easier
for high z. Here that happens to cancel the contamination,
which makes the clean items look harder, so the unadjusted screen looks
tidier than it is; in simulations with impact and no DIF, ignoring
impact flagged 40 per cent of the items. Purification works on the
adjusted model too:
A covariate given per response (age at each response, say) enters per response, so ability may also change with it within person.
The same contributions test non-uniform DIF
(parameters = "A"), both jointly
(c("B", "A")), the asymptotes of a 3PL or 4PL
("C", "D"), and other functional forms:
form = "spline" for smooth nonlinear relations,
"factor" for groups, and "fluctuation" for the
nonparametric double-max, Cramer-von Mises and maxLM tests, whose
p-values are simulated from the process’s exact covariance (so they hold
even though, with impact, information does not accrue evenly along
z).
fl <- scoreTest(sc, cov, parameters = "B", form = "fluctuation", items = 1:6, seed = 1)
fl$tests[, .(item, DM, p_DM, p_CvM, p_maxLM)]
#> item DM p_DM p_CvM p_maxLM
#> <char> <num> <num> <num> <num>
#> 1: 1 1.8376292 0.00059994 0.00049995 0.00119988
#> 2: 2 3.0885369 0.00009999 0.00009999 0.00009999
#> 3: 3 2.1821347 0.00009999 0.00009999 0.00009999
#> 4: 4 1.3834992 0.01979802 0.00309969 0.01499850
#> 5: 5 0.6622962 0.61438561 0.80519481 0.48651349
#> 6: 6 0.4790536 0.90809191 0.88611389 0.94105894parameters = "Ability" asks whether a covariate predicts
ability beyond what the ability model already holds. w here
is noise, and z is already in the model through its spline,
so neither should be flagged.
cov$w <- rnorm(nrow(cov))
scoreTest(sc, cov[c("id", "w")], parameters = "Ability")$tests
#> covariate item parameters form n LM df p
#> <char> <char> <char> <char> <int> <num> <int> <num>
#> 1: w Ability Ability linear 1500 0.07797347 1 0.7800627
#> p_fdr p_holm
#> <num> <num>
#> 1: 0.7800627 0.7800627Covariates can also be given per response – one row per row of the
data passed to fitIRT() – for example age at the time of
each response, or the item’s position in the test. Then
by = "all" tests one coefficient shared by every item, and
a vector of item groups tests one per group.
pos <- ave(dat$id, dat$id, FUN = function(x) sample(seq_along(x)))
scoreTest(sc, data.frame(position = pos), parameters = "B", by = "all")$tests
#> covariate group parameters form n LM df p
#> <char> <char> <char> <char> <int> <num> <int> <num>
#> 1: position all B linear 30000 4.098849 1 0.04291241
#> info_ratio ability_adjusted p_fdr p_holm
#> <num> <lgcl> <num> <num>
#> 1: 1.00822 FALSE 0.04291241 0.04291241