RMdimMartinLof: Martin-Löf Test of Unidimensionality
In easyRasch2: Psychometric Analysis with Rasch Measurement Theory
RMdimMartinLof R Documentation
Martin-Löf Test of Unidimensionality
Description
Likelihood-ratio test of unidimensionality against an a priori specified
multidimensional alternative, generalised to polytomous Rasch / partial
credit models (Christensen, Bjorner, Kreiner, & Petersen, 2002). The
p-value is obtained by parametric-bootstrap (Monte Carlo) sampling under
the unidimensional null, following Christensen & Kreiner (2007), because
the asymptotic chi-square approximation is biased toward conservatism for
realistic sample sizes – especially with polytomous items, where the
degrees of freedom can be very large.
Usage
RMdimMartinLof(
data,
partition,
iterations = 1000L,
stopping = c("none", "sequential"),
h = 50L,
alpha = 0.05,
parallel = TRUE,
n_cores = NULL,
verbose = FALSE,
seed = NULL
)
Arguments
data
A data.frame or matrix of item responses (0-based,
non-negative integers). Complete cases only: rows with any NA on the
items named in partition are dropped. Missingness on items outside
partition is ignored.
partition
The hypothesised partition of items into subscales. One of:
a list of column-name or column-index vectors, e.g.
list(c("I1","I2","I3"), c("I4","I5","I6"));
a vector of length ncol(data) indicating each item's subscale
(factor, character, or integer), e.g. c(1,1,1,2,2,2).
Each subscale must contain at least two items. Subscales must not
overlap; items not assigned to any subscale are dropped with a warning.
iterations
Integer. Maximum number of Monte Carlo iterations
(default 1000).
stopping
Character. "none" (default) runs all iterations.
"sequential" uses Besag & Clifford's (1991) sequential rule: stop as
soon as h simulated statistics have exceeded the observed value. The
sequential strategy substantially reduces compute time when H0 holds
but cannot be parallelised.
h
Integer. Sequential-stopping count threshold (default 50).
Ignored when stopping = "none".
alpha
Numeric in (0, 1). Nominal significance level used only for
the rejected flag in the result; default 0.05.
parallel
Logical. Use parallel processing via mirai (default
TRUE). Ignored when stopping = "sequential".
n_cores
Integer or NULL. Number of parallel workers. When NULL,
getOption("mc.cores") is checked first; if neither is set, falls back
to sequential with a warning.
verbose
Logical. Show a progress bar (default FALSE).
seed
Integer or NULL. Random seed for reproducibility. Items are
processed internally in a fixed (alphabetical) order, so the same seed
reproduces the same p-value regardless of how the data's columns are
arranged and regardless of the parallel setting. See
easyRasch2-reproducibility for what this guarantees and how it
interacts with parallel.
Details
This is not a routine screening tool. The test requires an a priori
partition of items into subscales; using it post-hoc on, e.g., the
partition suggested by RMdimResidualPCA()'s PC1 sign would inflate the
Type-I error rate. Both source papers state this explicitly.
Test statistic. With items partitioned into D subscales, total score
t and subscores (t_1, \ldots, t_D) (Christensen et al. 2002,
eq. 22):
T = 2\Bigl[\sum_{t_1, \ldots, t_D}
n_{t_1, \ldots, t_D}\log(n_{t_1, \ldots, t_D}/N)
- \sum_t n_t\log(n_t/N)
- \ell_C(\hat{\epsilon}) + \sum_d \ell_C(\hat{\epsilon}^{(d)})\Bigr]
where \ell_C is the conditional log-likelihood and the
\hat{\epsilon}^{(d)} are CML estimates on the d-th subscale alone.
CML fits use psychotools::raschmodel() (RM) or psychotools::pcmodel()
(PCM) for speed.
Monte Carlo sampling under H0. Following Christensen & Kreiner (2007):
(a) sample N total scores from the empirical score distribution
n_t/N; (b) for each sampled score, sample an item-response vector
from the conditional distribution p(x \mid t, \hat{\epsilon}) given
by eq. 4 of the paper. Step (b) uses the recursive \gamma-function
algorithm (Andersen, 1995, eq. 15.22 and 15.27) for dichotomous and
polytomous items alike: items are filled in from the last backwards, each
conditional on the score still to be distributed over the items before it.
The nested \gamma^{(1)}, \ldots, \gamma^{(n)} depend only on the
item parameters, which are fixed across the run, so the recursion is
evaluated once rather than per person.
The faster dichotomous shortcut described by Christensen & Kreiner (2007,
p. 23), choosing s items one at a time with probabilities
\phi_i / \sum_{\mathrm{remaining}} \phi, is deliberately not used:
successive sampling of that kind follows a Wallenius-type distribution,
whereas the Rasch conditional distribution weights a set of items by the
product of their easinesses. See the note on validation below.
Iterations that fail (e.g., simulated dataset has an empty category for an
item) are silently dropped.
Item parameters are estimated once on the observed data and held fixed
across MC iterations. Christensen & Kreiner (2007) use the extended
likelihood function (Tjur, 1982) with the empirical score distribution as
a non-parametric estimate of the latent distribution, so no distributional
assumption about \theta is needed.
Missing data. Complete cases only. The statistic is built from the
joint table of subscores, so a respondent needs a defined subscore on
every subscale; there is no partial-missingness path, and the SAS macro
this implementation was validated against has the same requirement. Rows
with NA are dropped after items outside partition have been
removed, so missingness on items the test does not use costs no
respondents. sample_n reports the respondents analysed and
sample_n_total the raw input rows.
Validation. The statistic, the conditional log-likelihoods, the
conditional sampler and the expected counts and residuals of
RMdimMartinLofResiduals have been checked against the pml
SAS macro (Christensen, 2004), kindly shared by Karl Bang Christensen.
The statistic and the residual table agree with the macro to numerical
precision. The sampler follows the macro's recursion rather than the
dichotomous shortcut published in Christensen & Kreiner (2007).
Value
A list with components:
T_obsObserved Martin-Löf likelihood-ratio statistic.
p_valueMonte Carlo p-value with (n_exceed + 1) / (n + 1)
correction. The attainable p-values are k / (n + 1) for
k = 1, ..., n + 1, so the p-value's resolution is limited by the
number of iterations: with 100 iterations the smallest attainable
value is 1/101 = 0.0099. A reported p-value equal to the floor
(see p_value_floor) means no simulated statistic reached the
observed one and should be read as "p < floor" – the true p-value
may be much smaller; increase iterations for finer resolution.
p_value_floorThe smallest attainable p-value,
1 / (actual_iterations + 1).
actual_iterationsNumber of successful MC iterations
completed.
rejectedLogical: is p_value < alpha?
partitionNormalised partition (list of integer indices).
n_subscalesNumber of subscales.
is_polytomousWhether a PCM was fitted.
sample_nNumber of complete cases analysed, counted over the
items in partition only.
sample_n_totalNumber of respondents in the raw input data,
before the complete-case filter.
sample_has_naLogical. Whether the partitioned items contained
any missing values.
n_itemsNumber of items.
stoppingThe stopping strategy used.
hThe sequential-stopping count, or NA for stopping = "none".
T_repNumeric vector of successful MC test statistics.
wle_scoresdata.frame with one row per person and one column
per subscale (subscale_1_wle, ..., subscale_D_wle), giving Warm's
Weighted Likelihood Estimate of theta from a CML fit on each subscale
alone. Persons whose subscore equals the minimum or maximum on a
subscale produce non-finite WLEs (Inf / -Inf) and are excluded
from wle_correlation pairwise.
wle_correlationdata.frame of pairwise Pearson correlations
between subscale WLEs, with columns subscale_a, subscale_b, r,
ci_lower, ci_upper (95% CI from stats::cor.test), p_value,
and n (number of persons with finite WLEs on both subscales). One
row per pair; for D = 2, a single row. Useful as an effect-size
companion to p_value – a rejected test with r near 1 indicates a
small effect; r clearly below 1 indicates substantive
multidimensionality.
References
Andersen, E. B. (1995). Polytomous Rasch models and their estimation. In
G. H. Fischer & I. W. Molenaar (Eds.), Rasch models: Foundations, recent
developments, and applications (pp. 271-291). Springer-Verlag.
Besag, J., & Clifford, P. (1991). Sequential Monte Carlo p-values.
Biometrika, 78(2), 301-304. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/78.2.301")}
Christensen, K. B. (2004). pml: A SAS macro for testing
unidimensionality in polytomous Rasch models (Technical note). National
Institute of Occupational Health, Denmark, and Department of
Biostatistics, University of Copenhagen.
Christensen, K. B., Bjorner, J. B., Kreiner, S., & Petersen, J. H. (2002).
Testing unidimensionality in polytomous Rasch models. Psychometrika,
67(4), 563-574. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02295132")}
Christensen, K. B., & Kreiner, S. (2007). A Monte Carlo approach to
unidimensionality testing in polytomous Rasch models. Applied
Psychological Measurement, 31(1), 20-30.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/0146621605286204")}
See Also
RMdimResidualPCA, RMdimResidualPCACutoff
Examples
set.seed(1)
# Build 2-dimensional polytomous data: 4 items per subscale, 5 categories
n <- 400
theta1 <- rnorm(n)
theta2 <- 0.6 * theta1 + sqrt(1 - 0.6^2) * rnorm(n)
make_pcm <- function(theta, n_items, taus) {
sapply(seq_len(n_items), function(j) {
# ... toy simulation here
sample(0:4, n, replace = TRUE)
})
}
dat <- cbind(make_pcm(theta1, 4, NULL), make_pcm(theta2, 4, NULL))
colnames(dat) <- paste0("I", 1:8)
# Few iterations for a fast example; use 1000+ in real analyses
RMdimMartinLof(dat,
partition = list(c("I1","I2","I3","I4"),
c("I5","I6","I7","I8")),
iterations = 100, parallel = FALSE, seed = 1)
# Sequential stopping: stop as soon as h = 25 simulated statistics exceed
# the observed one (cuts compute time under H0).
RMdimMartinLof(dat,
partition = c(1,1,1,1,2,2,2,2),
iterations = 200, stopping = "sequential", h = 25,
seed = 1)
easyRasch2 documentation built on Sept. 13, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| RMdimMartinLof | R Documentation |
Martin-Löf Test of Unidimensionality
Description
Likelihood-ratio test of unidimensionality against an a priori specified multidimensional alternative, generalised to polytomous Rasch / partial credit models (Christensen, Bjorner, Kreiner, & Petersen, 2002). The p-value is obtained by parametric-bootstrap (Monte Carlo) sampling under the unidimensional null, following Christensen & Kreiner (2007), because the asymptotic chi-square approximation is biased toward conservatism for realistic sample sizes – especially with polytomous items, where the degrees of freedom can be very large.
Usage
RMdimMartinLof(
data,
partition,
iterations = 1000L,
stopping = c("none", "sequential"),
h = 50L,
alpha = 0.05,
parallel = TRUE,
n_cores = NULL,
verbose = FALSE,
seed = NULL
)
Arguments
data |
A data.frame or matrix of item responses (0-based,
non-negative integers). Complete cases only: rows with any |
partition |
The hypothesised partition of items into subscales. One of:
|
iterations |
Integer. Maximum number of Monte Carlo iterations
(default |
stopping |
Character. |
h |
Integer. Sequential-stopping count threshold (default |
alpha |
Numeric in (0, 1). Nominal significance level used only for
the |
parallel |
Logical. Use parallel processing via |
n_cores |
Integer or |
verbose |
Logical. Show a progress bar (default |
seed |
Integer or |
Details
This is not a routine screening tool. The test requires an a priori
partition of items into subscales; using it post-hoc on, e.g., the
partition suggested by RMdimResidualPCA()'s PC1 sign would inflate the
Type-I error rate. Both source papers state this explicitly.
Test statistic. With items partitioned into D subscales, total score
t and subscores (t_1, \ldots, t_D) (Christensen et al. 2002,
eq. 22):
T = 2\Bigl[\sum_{t_1, \ldots, t_D}
n_{t_1, \ldots, t_D}\log(n_{t_1, \ldots, t_D}/N)
- \sum_t n_t\log(n_t/N)
- \ell_C(\hat{\epsilon}) + \sum_d \ell_C(\hat{\epsilon}^{(d)})\Bigr]
where \ell_C is the conditional log-likelihood and the
\hat{\epsilon}^{(d)} are CML estimates on the d-th subscale alone.
CML fits use psychotools::raschmodel() (RM) or psychotools::pcmodel()
(PCM) for speed.
Monte Carlo sampling under H0. Following Christensen & Kreiner (2007):
(a) sample N total scores from the empirical score distribution
n_t/N; (b) for each sampled score, sample an item-response vector
from the conditional distribution p(x \mid t, \hat{\epsilon}) given
by eq. 4 of the paper. Step (b) uses the recursive \gamma-function
algorithm (Andersen, 1995, eq. 15.22 and 15.27) for dichotomous and
polytomous items alike: items are filled in from the last backwards, each
conditional on the score still to be distributed over the items before it.
The nested \gamma^{(1)}, \ldots, \gamma^{(n)} depend only on the
item parameters, which are fixed across the run, so the recursion is
evaluated once rather than per person.
The faster dichotomous shortcut described by Christensen & Kreiner (2007,
p. 23), choosing s items one at a time with probabilities
\phi_i / \sum_{\mathrm{remaining}} \phi, is deliberately not used:
successive sampling of that kind follows a Wallenius-type distribution,
whereas the Rasch conditional distribution weights a set of items by the
product of their easinesses. See the note on validation below.
Iterations that fail (e.g., simulated dataset has an empty category for an item) are silently dropped.
Item parameters are estimated once on the observed data and held fixed
across MC iterations. Christensen & Kreiner (2007) use the extended
likelihood function (Tjur, 1982) with the empirical score distribution as
a non-parametric estimate of the latent distribution, so no distributional
assumption about \theta is needed.
Missing data. Complete cases only. The statistic is built from the
joint table of subscores, so a respondent needs a defined subscore on
every subscale; there is no partial-missingness path, and the SAS macro
this implementation was validated against has the same requirement. Rows
with NA are dropped after items outside partition have been
removed, so missingness on items the test does not use costs no
respondents. sample_n reports the respondents analysed and
sample_n_total the raw input rows.
Validation. The statistic, the conditional log-likelihoods, the
conditional sampler and the expected counts and residuals of
RMdimMartinLofResiduals have been checked against the pml
SAS macro (Christensen, 2004), kindly shared by Karl Bang Christensen.
The statistic and the residual table agree with the macro to numerical
precision. The sampler follows the macro's recursion rather than the
dichotomous shortcut published in Christensen & Kreiner (2007).
Value
A list with components:
T_obsObserved Martin-Löf likelihood-ratio statistic.
p_valueMonte Carlo p-value with
(n_exceed + 1) / (n + 1)correction. The attainable p-values arek / (n + 1)fork = 1, ..., n + 1, so the p-value's resolution is limited by the number of iterations: with 100 iterations the smallest attainable value is1/101 = 0.0099. A reported p-value equal to the floor (seep_value_floor) means no simulated statistic reached the observed one and should be read as "p < floor" – the true p-value may be much smaller; increaseiterationsfor finer resolution.p_value_floorThe smallest attainable p-value,
1 / (actual_iterations + 1).actual_iterationsNumber of successful MC iterations completed.
rejectedLogical: is
p_value < alpha?partitionNormalised partition (list of integer indices).
n_subscalesNumber of subscales.
is_polytomousWhether a PCM was fitted.
sample_nNumber of complete cases analysed, counted over the items in
partitiononly.sample_n_totalNumber of respondents in the raw input data, before the complete-case filter.
sample_has_naLogical. Whether the partitioned items contained any missing values.
n_itemsNumber of items.
stoppingThe stopping strategy used.
hThe sequential-stopping count, or
NAforstopping = "none".T_repNumeric vector of successful MC test statistics.
wle_scoresdata.frame with one row per person and one column per subscale (
subscale_1_wle, ...,subscale_D_wle), giving Warm's Weighted Likelihood Estimate of theta from a CML fit on each subscale alone. Persons whose subscore equals the minimum or maximum on a subscale produce non-finite WLEs (Inf/-Inf) and are excluded fromwle_correlationpairwise.wle_correlationdata.frame of pairwise Pearson correlations between subscale WLEs, with columns
subscale_a,subscale_b,r,ci_lower,ci_upper(95% CI fromstats::cor.test),p_value, andn(number of persons with finite WLEs on both subscales). One row per pair; for D = 2, a single row. Useful as an effect-size companion top_value– a rejected test withrnear 1 indicates a small effect;rclearly below 1 indicates substantive multidimensionality.
References
Andersen, E. B. (1995). Polytomous Rasch models and their estimation. In G. H. Fischer & I. W. Molenaar (Eds.), Rasch models: Foundations, recent developments, and applications (pp. 271-291). Springer-Verlag.
Besag, J., & Clifford, P. (1991). Sequential Monte Carlo p-values. Biometrika, 78(2), 301-304. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/78.2.301")}
Christensen, K. B. (2004). pml: A SAS macro for testing unidimensionality in polytomous Rasch models (Technical note). National Institute of Occupational Health, Denmark, and Department of Biostatistics, University of Copenhagen.
Christensen, K. B., Bjorner, J. B., Kreiner, S., & Petersen, J. H. (2002). Testing unidimensionality in polytomous Rasch models. Psychometrika, 67(4), 563-574. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02295132")}
Christensen, K. B., & Kreiner, S. (2007). A Monte Carlo approach to unidimensionality testing in polytomous Rasch models. Applied Psychological Measurement, 31(1), 20-30. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/0146621605286204")}
See Also
RMdimResidualPCA, RMdimResidualPCACutoff
Examples
set.seed(1)
# Build 2-dimensional polytomous data: 4 items per subscale, 5 categories
n <- 400
theta1 <- rnorm(n)
theta2 <- 0.6 * theta1 + sqrt(1 - 0.6^2) * rnorm(n)
make_pcm <- function(theta, n_items, taus) {
sapply(seq_len(n_items), function(j) {
# ... toy simulation here
sample(0:4, n, replace = TRUE)
})
}
dat <- cbind(make_pcm(theta1, 4, NULL), make_pcm(theta2, 4, NULL))
colnames(dat) <- paste0("I", 1:8)
# Few iterations for a fast example; use 1000+ in real analyses
RMdimMartinLof(dat,
partition = list(c("I1","I2","I3","I4"),
c("I5","I6","I7","I8")),
iterations = 100, parallel = FALSE, seed = 1)
# Sequential stopping: stop as soon as h = 25 simulated statistics exceed
# the observed one (cuts compute time under H0).
RMdimMartinLof(dat,
partition = c(1,1,1,1,2,2,2,2),
iterations = 200, stopping = "sequential", h = 25,
seed = 1)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.