predict(fit, newdata) holds the estimated parameters (and, for biclustering,
the field memberships of the items) fixed and runs one E-step on new responses,
returning the membership probabilities, the estimated class or rank and each
respondent's log-likelihood. It covers LCA() (binary, nominal, ordinal,
rated), LRA() (binary, ordinal, rated) and Biclustering() / Ranklustering
(binary, nominal, ordinal, rated). Items are matched by name, missing responses
are left out of the likelihood, and a response in a category the fitted data
never used for that item is treated as missing with a warning. On the data the
model was fitted to it reproduces the stored membership.
To make this possible the fits now record how they coded the categories:
CategoryMap (the raw codes behind each item's renumbered categories) in the
nominal, ordinal and rated LCA() and LRA(), CategoryShift in the nominal
and ordinal Biclustering(), and LogPrior (the rank prior set by
trapezoidal) in the ordinal and rated LRA(). Fits made with earlier versions
lack these; predict() asks for a refit where it cannot rebuild them.
LCA(), the binary LRA() (GTM and isotonic) and every Biclustering()
returned the class posterior of the last E-step, which was computed from the
parameters of the previous cycle; the parameters themselves then moved once
more. The membership and the parameters returned did not belong together, and
an E-step with the returned parameters on the same data did not reproduce the
membership. Each now runs one more E-step after the loop with the final
parameters (the ordinal and rated LRA() and the SOM already did). Memberships,
class estimates and the fit quantities built on them move slightly.
When the biclustering EM stops because its objective decreased, it restored the
reference profiles of the previous cycle but kept the new field (and, for the
nominal and ordinal data, class) memberships. All are now rolled back together.
On J35S515 (binary Ranklustering, GTM, 6 x 5) the field memberships of a few
items and the fit indices change (log-likelihood -7273.06 to -7273.90). Shojima's
Mathematica reference implementation had both defects; it was corrected in step
and the reference fixtures regenerated.
Every latent class, latent rank and biclustering fit now carries four
additional quantities for model comparison (elbo, log_lik_marginal,
entropy, ICL; rated models report the nominal layer with the suffix
_nominal). Their definitions are documented in the return value of each
function and in R/00_VariationalFit.R:
elbo is the evidence lower bound of the observed-data log-likelihood under
the posterior memberships the EM produced. For a one-mode model (LCA, LRA) the
posterior is exact and this is the observed-data log-likelihood itself; for
biclustering the respondent and item posteriors are taken independent and the
value is a lower bound.log_lik_marginal is the observed-data log-likelihood of the respondent-side
mixture with the items fixed at their MAP field assignment.entropy is the posterior entropy of the memberships. elbo minus the
entropy is the expected complete-data log-likelihood.ICL is the integrated completed likelihood criterion (Biernacki, Celeux and
Govaert, 2000) in the BIC-type approximation used for the latent block model
by Keribin, Brault, Celeux and Govaert (2015), on the same scale as AIC and
BIC here (smaller is better).The default fit of Biclustering() on ordinal data changes. The EM is
deterministic but converges to a local optimum, and the start used up to now
(fields = items split by their median, block probabilities = a ladder common
to every field) stopped at one in over 80% of simulated data sets. The new
argument init chooses the start:
"auto" (default) runs from both "quantile" and "ward" and keeps the fit
with the larger ELBO."quantile": classes from equal-count bins of the respondent score."ward": classes from Ward clustering of the respondent score, which does
not assume equal class sizes."legacy": the earlier start; reproduces earlier results."auto" also adds nstart random starts (default 10; seed sets them, and
the session's random number stream is left untouched). The two data-based
starts always run to convergence, so the random starts can only raise the
ELBO; the random starts are run for 100 EM cycles and only the best 3 continue.
A random start has random classes; its fields are conf if given, and
otherwise those of the "quantile" start for half of the starts and random for
the other half. On real data the likelihood surface often had many more local
optima than in simulations, and the data-based starts alone were often beaten
by random ones; raise nstart for such data.
In "quantile" and "ward" the fields start from a Ward clustering of the
items' class-mean profiles (centred per item), and the block probabilities from
the category frequencies of each block. Both starts are invariant to the order
of the rows. In simulations they reached the level of a start at the true
partition: "quantile" under equal class sizes, "ward" under unequal ones,
and "auto" under both, at about the cost of one legacy fit. With conf the
given fields are kept, and with conf_class the given classes are the single
start. The fit records init, init_selected and init_elbo.
The isotonic fit of ordinal Ranklustering constrains the field reference
profiles, but SOACflg and WOACflg were judged on a test reference profile
built from data-weighted counts, which the constraint does not reach. A single
dip from sampling noise turned SOAC off even when every FRP was monotone, and
more often with more ranks. TRP of Biclustering() (ordinal) and
Biclustering_IRM() (ordinal) is now computed from the model, as for binary
and rated data: the sum over fields of field size times the field expected
score. A monotone FRP thus gives a monotone TRP, and the isotonic fit always
satisfies SOAC. The checks allow the feasibility tolerance of the isotonic step
(1e-6). The profile from the data is returned as TRP_empirical (mean observed
score of the class members, summed over fields); its departure from TRP
shows how far the model is from the data. TRP values change.
Up to now the polytomous biclustering models renumbered each item's
observed categories to 1..k, so on an item where nobody chose a category the
codes above it moved down. Since 1.15.0, Biclustering() and
Biclustering_IRM() on nominal, ordinal and rated data did this through
remap_category_codes(), which was meant only to make 0-based codes usable.
In these models the items of a field share the block's category
probabilities, so code q must mean the same category on every item; on an item
whose lowest category was empty the renumbering misaligned every response of
that item with the rest of its field. This is common in practice: in
simulated data most data sets had an item whose bottom category was empty.
LCA(), LRA() on ordinal and rated data, GRM())
keep renumbering each item, because a category with no responses leaves its
threshold or probability unidentified, but now warn which categories of which
items were dropped.dataFormat() codes a factor by its declared levels: an unused level keeps
its code, instead of being dropped and closing the gap. The category labels
still list the observed levels only.M2() and add_M2() on a rated LCA() fit compared the raw codes against
the renumbered profile; they now use the renumbered codes.DistractorAnalysis() on a rated biclustering fit lost a correct answer
coded 0.The log_lik of Biclustering() and Biclustering_IRM() on nominal and
ordinal data, and the AIC, BIC and CAIC built on it, change. Up to 2.1.0 these
reported sum_sj log sum_fc m_jf m_sc pi_fcq: the probability of each observed
cell under a mixture of the block probabilities weighted by that respondent's
and that item's posterior memberships. That is the likelihood of a
mixed-membership model with the memberships plugged in as free parameters, not
of the mixture model being fitted, and it is not the quantity the EM ascends.
Because the plugged-in memberships are never counted in the penalty, BIC on
this quantity favours more blocks the larger the sample, which is what a
simulation study of the grid search showed (over-selection of the number of
classes rising with the number of respondents). The reported value is now the
evidence lower bound, which the EM also monitors for convergence. The binary
implementation keeps the expected log-likelihood of the reference
implementation (Shojima, 2022) as its log_lik; the new fields are added
alongside it there, so its fit indices are unchanged.
Consequences for ordinal and nominal data:
log_lik, TestFitIndices$model_log_like, AIC, CAIC, BIC and the
chi-square based indices of Biclustering() and Biclustering_IRM() take
different values.verbose trace prints elbo instead of log_lik.dataFormat() now stores a missing cell as 0 in U (binary and rated
data) and keeps the indicator in Z; Q (polytomous data) still carries -1.
Until 2.1.0 U held -1 for a missing cell, and every routine that summed U
without masking it by Z counted each missing cell as minus one correct
response. The audit before this release found that in the binary LCA and LRA
cores (emclus(), emclus_isotonic(), somclus(), item_log_lik()), in
BNM(), LDLRA(), LDB() and BINET(), and in the quantile reference of
rated LRA(). With the new convention those sums are right by construction.
Independently of the sentinel, the following counted a missing cell as an
incorrect response or divided by the wrong count and are corrected:
ItemFit(), TestFit() and TestFitSaturated() weighted the null model by
the number of respondents instead of the number of responses, and the
benchmark model took "group size minus correct" as the incorrect count. With
missing data the chi-squares came out too large (or negative for IRT); every
model's NFI...RMSEA and chi-square based AIC/BIC on such data change.LRA(method = "SOM") did not receive Z at all; its C++ epoch now skips
missing cells in both the winner search and the update.LDB() counted the field size, not the number of observed items per field.StudentAnalysis() reported the number of items, not of responses, as NR.Biclustering() on rated data ordered the ranks by correct responses per
respondent rather than per observed response, which inverted the order when
missing cells were unevenly spread across classes.BFRP in ordinal biclustering and IRM averaged over
all cells, counting a missing cell as category 0.M2() divided the observed margins by the full sample size; with missing
data every margin of an item was deflated by its response rate and M2 was
systematically inflated (20% MCAR from a true model gave p = 0 every time).
Margins now use available cases, and caveat says so when data are missing.LRA() on rated data with minFreqRatio > 0 put missing cells into the
collapsed category.LRA(method = "GTM") flipped the direction of the ranks on almost
every cycle: the orientation test compared two boundary rows instead of
the first and last rank columns, so whether the result came out ascending or
descending depended on the number of EM cycles (the convergence tolerance).
The isotonic default was unaffected.BNM(), BNM_GA(), BNM_PBIL(), LDLRA() and LDLRA_PBIL() fitted a
different graph than the one given whenever the item labels were not in
alphabetical order: the adjacency matrix was reordered alphabetically
while the data were read by position. With labels such as Item1 ... Item15 or Q1 ... Q12 the parents of an item were attached to another
item. The matrix now stays in item-label order. LDLRA_PBIL() in addition
evaluated its candidates by position but returned them by name, so the
fitted DAG was not the one the search had scored; both now use the same
interpretation. It also ignored method, beta1 and beta2 for the final
fit and failed when filename was given.LCA() on rated data looked up the keyed category by its raw code; with
0-based or gappy codes IRP, TRP and the binary-layer fit indices came
from another category.ItemStatistics() on ordinal data with 0-based or gappy codes (such as the
bundled J15S3810) returned wrong thresholds, entropies and item-total
correlations because tabulate() skips code 0.ItemStatistics() and InterItemAnalysis() stopped on rated data
(J35S5000).GridSearch(fun = "LRA") passed ncls where LRA() expects nrank, so
every cell was fitted with two ranks and "2" was always optimal. The grid
can now also be searched on ICL, elbo and log_lik_marginal.M2() on an ordinal LRA() fit compared the observed categories by raw
code against profiles in remapped order, giving wrong values for 0-based
codes.MutualInformation() returned NaN for an empty cell (0 log 0).LCA(), Biclustering() and the EM of Biclustering_IRM() with
alpha < 1 produced negative probabilities; the MAP numerator is now
clamped and renormalised as in the ordinal M-step.LCA(), ordinal and rated LRA(), and all Biclustering() EMs ran
maxiter + 1 cycles and reported that count.LRA() stopped when a score quantile group had a single
respondent, or when every item had two categories (a dimension was dropped
by indexing). Rated LRA()'s default maxiter is now 1000, as elsewhere.GRM() built its item fit from the EAP-plugged likelihood, which is not
nested in the benchmark and gave negative chi-squares even without missing
data; it now uses the expected log-likelihood under the posterior, as
IRT() does. It also warns when the optimiser did not converge, and its
C++ score function no longer underflows on long tests.OmegaCoefficient() started from random loadings: with three items it
stopped nine times out of ten with an empty error message, and otherwise
returned a different value on every call. It starts from the first
principal component and warns instead of stopping.IRT() validates model (2, 3 or 4) instead of failing obscurely, and
stops rather than silently reusing the previous item's parameters when the
optimiser fails for an item.Biclustering() on rated data failed with conf or conf_class because it
ignored the class and field counts the nominal fit had derived from them.Biclustering_IRM() failed when no class fell below minSize
(EMt undefined), nominal Biclustering_IRM() failed when the modal
categories were not exactly 1..k, and all three IRMs failed obscurely when
every class was below minSize (they now stop at one class with a warning).mic = TRUE in ordinal biclustering reordered the reference profiles
without the memberships, so a class column no longer described the same
respondents; both are permuted together (and nothing is reordered under
conf_class). In ordinal IRM the same fix stops the EM from terminating on
the reordering cycle.LDLRA() with beta1 = beta2 = 1 stopped whenever some parent pattern was
unobserved; such cells are now excluded. An adjacency matrix without
dimnames is accepted (taken in item-label order) instead of silently losing
every edge. BNM_PBIL() and LDLRA_PBIL() no longer halve the number of
survivors with a warning on every call, run exactly maxGeneration
generations, and LDLRA_PBIL() keeps the stored fitness of its elites
aligned after sorting. LDLRA_PBIL() rejects an unknown method.dataFormat(): a numeric ID column that is unique in every row is taken as
an ID; character response columns (a missing-value symbol in a CSV column)
are converted to numeric; na may be a vector; Z may be a data frame;
factor columns keep only their observed levels, their labels stay in item
order when mixed with numeric columns, and the labels of an excluded item
are dropped with it. longdataFormat() returns one label vector per item,
so LCA() runs on its output.print() of a rated LCA() result that had been through add_M2() stopped
with "non-numeric argument to mathematical function": the M2 block carries a
character caveat field that was rounded along with the numbers.plot(type = "IRP", items = <one item>) stopped with "argument is of length
zero" for LCA, LRA and LDLRA, because a single row of the profile matrix
dropped to a vector. With several items the panel titles counted 1, 2, ...
instead of naming the items plotted; they now show the item labels.print() of a BNM() fit with an unobserved parent pattern turned the
"NaN(0/0)" marker into NA with a coercion warning.plot(type = "FCRP") and plot(type = "FCBR") failed with "figure margins
too large" at the default nc = 1 once a fit had more than a handful of
fields, which an IRM fit decides from the data. The columns now default to a
near-square layout when nc is left at 1 and more than three fields are
drawn; nr is honoured when given.plot(items = 0) is now rejected outside IRT and GRM, where 0 selects the
test-level curve; elsewhere it indexed past the first item.Biclustering_IRM() now carries elbo, log_lik_marginal, entropy
and ICL like the other biclustering fits (its assignments are hard, so the
entropy is zero).print() of these fits shows a "Variational Fit" block.LRA() on ordinal and rated data now returns log_lik at the top level, and
evaluates it (and TestFitIndices$model_log_like) at the returned estimates.
It used to be the value computed inside the last EM cycle before that
cycle's M-step, i.e. one iterate behind the profiles the object carries; the
two now describe the same state. The difference is of the order of the
convergence tolerance (about 2 log-likelihood units on J15S3810).A minor release. The estimates returned by LRA(method = "SOM") change, so the
fixes below could not ship as a patch. They came out of reading the R code
against the Mathematica routine it was ported from
(Module_LRA.wl): the port had grown a seed argument and a convergence
report that the original never had, and had moved the monotonicity sort to the
end of each epoch.
Ordinal biclustering rolled back only half of its state. When an EM sweep
failed to improve the likelihood, Biclustering() on ordinal data restored the
category reference array BCRM but left the boundary array BBRM at the
rejected values. The two are then one iteration apart, and since
test_log_lik is computed from BCRM while the shape-restricted parameter
count is read off BBRM, the reported nparam, df, AIC, BIC and CAIC
could describe a different iterate than the likelihood does. It affects
method = "R" with estimation = "isotonic" only; the other branches derive
the count from the filter matrix or the rank and field counts instead, and the
binary implementation rolls back the array it later reads. Both arrays are now
saved and restored together. In practice the rejected sweep differs from the
accepted one by around 1e-10 by the time the branch is reached, so the counted
number of distinct boundary levels -- and hence every figure above -- came out
unchanged in the cases tested; the fix removes the inconsistency rather than a
visible error.
The SOM method reused one frozen presentation order whenever seed was
given. somclus() reseeded inside the epoch loop with set.seed(seed), so
every epoch presented the respondents in exactly the same order. Online
learning relies on that order being redrawn: with it fixed, the bias of one
particular ordering is never averaged out. The original Mathematica routine
reseeds each epoch with SeedRandom[Total[uuu] + somt], and the seed = NULL
path already matched it; only the user-supplied seed was affected. The epoch
number is now added to the seed, so a run stays reproducible while the order
changes from epoch to epoch. Estimates from method = "SOM" with an explicit
seed change accordingly.
mic = TRUE sorted the rank reference matrix once per epoch instead of once
per respondent. The original applies Sort /@ refmat immediately after each
respondent is presented, which changes the matrix that decides the next
winner; deferring it to the end of the epoch is a different algorithm, not a
cheaper form of the same one. The sort -- and the conf constraint alongside
it -- now runs after every respondent. Estimates from method = "SOM" with
mic = TRUE or conf change accordingly.
Ties in the SOM winner search are resolved toward the larger rank, as in
the original (Sort[Transpose[{mlrank, clsnum}]][[-1]]). which.max() had
been picking the smaller rank.
method = "SOM" no longer warns that it may not have converged. SOM runs
an annealing schedule for maxiter epochs and has no convergence criterion --
neither does the original implementation -- so finishing the schedule is
normal termination. The warning fired on every default run and converge was
always FALSE, which left no way to tell a healthy run from a broken one.
converge is now TRUE unless BIC.check early stopping was requested and
failed to trigger within ten times maxiter. The messages on that path, and
the documentation of BIC.check, now describe it as early stopping rather
than a convergence test; its threshold is unchanged.
LRA(method = "SOM") no longer overwrites the caller's random number
stream. It reseeds once per epoch and never restored .Random.seed, so any
simulation that called it silently lost its own stream. The state is now
restored on exit.
src/som_core.cpp). One epoch of online
updates -- winner search, neighbourhood update, conf, mic and the prior
update -- now runs in compiled code, while the presentation order is still
drawn in R so that set.seed() keeps governing reproducibility. On
J15S500 with nrank = 6 and maxiter = 1000, estimation went from 5.2 s to
1.9 s with mic = FALSE, and from 103.5 s to 2.6 s with mic = TRUE (the
per-respondent sort that mic now requires would otherwise have made it
twenty times slower). The rank posterior is also computed once after the
schedule instead of every epoch, since only BIC.check reads it in between.grconvertX(), grconvertY() and rasterImage(), used by the rasterised
Array plot since 2.0.1, are now imported from graphics in NAMESPACE.
This clears the R CMD check note about undeclared global functions.A patch release: two bug fixes, no new features and no change to any estimate. Released on GitHub only (2026-08-25); these fixes reached CRAN with 2.1.0.
Respondent IDs and item labels now reach every model's output. They
survived only where the matrix a field happened to be derived from carried
dimnames, so which fields were labelled differed from model to model:
Biclustering() on ordinal data lost them completely (FieldEstimated,
FieldMembership and ClassMembership all came back as bare numbers), and
ClassEstimated was unnamed in every model, including the binary ones that
labelled everything else. The respondent axis was the more affected of the
two, which is the opposite of what the code looks like -- Students had its
row names set by hand, and nothing else did.
Labelling is now done once, in one place, immediately before each constructor
returns: respondent-indexed fields (ClassEstimated, ClassMembership,
SmoothedMembership, Students, ability, ...) take the IDs, item-indexed
fields (FieldEstimated, FieldMembership, ...) take the item labels, and the
membership matrices get column names (Class1/Rank1..., Field1...).
dataFormat() also puts the IDs on the rows of U, Q and Z; it already
put the item labels on the columns.
FieldAnalysis is deliberately excluded: its rows are sorted by correct
response rate and field, so attaching labels positionally would attach the
wrong ones. It keeps the labels it inherits, in its own order.
Only names are added -- no estimate changes.
plot(type = "Array") no longer washes rows out to white. Every cell was
drawn as a rect() with a white border, and a border cannot be thinner than
one device pixel: once the respondents outnumbered the pixels available to
them, the borders covered the fill. Which rows disappeared depended on where
the cell boundaries fell on the pixel grid, so the loss came out mottled
rather than uniform. Measured on an all-black 400x600 plot, mean luminance
rose from 0.01 at 50 rows to 0.29 at 821 rows -- roughly a third of the ink
gone at the sizes a real test data set reaches.
rasterImage()) instead of
one rect() call per cell. A 3,810 x 15 data set took 114,300 drawing calls
for the two panels and now takes two; the plot renders in 0.03 s. Cell
borders are drawn -- as grid lines, nrows + ncols of them rather than
nrows * ncols rectangles -- only when a cell is at least 6 device pixels on
both sides, which is where they can be seen without swallowing the cell.The version number was raised to 2.0.0 because the EM convergence fix below changes the estimates of every EM-based model: the same code on the same data now returns different numbers, and no longer the numbers printed in Shojima (2022). A minor bump would have understated that. The deprecation removals planned for 2.0.0 ride along, so that the break happens once.
The names deprecated in earlier releases are gone. Each carried a deprecation warning for several releases: the function renames since 1.7.0, the field-name aliases since 1.8.0.
Functions removed, with their replacements:
| Removed | Use instead |
|---|---|
| IRM() | Biclustering_IRM() |
| StrLearningGA_BNM() | BNM_GA() |
| StrLearningPBIL_BNM() | BNM_PBIL() |
| StrLearningPBIL_LDLRA() | LDLRA_PBIL() |
Field aliases removed from returned objects. The snake_case names have been present alongside them since 1.8.0, so code that already reads the new names is unaffected:
| Removed | Use instead |
|---|---|
| $Nclass | $n_class |
| $Nfield | $n_field |
| $Nrank | $n_rank |
| $N_Cycle | $n_cycle |
| $LogLik | $log_lik |
GridSearch() still accepts "LogLik" as a value for its index argument;
that is an input alias for the criterion name, not a field on the result.
Biclustering() ignored maxiter on binary data. The binary branch set
its own iteration cap to 1000 outright, so the argument was accepted and
discarded; the polytomous branches always passed it through. The default is
now 1000 to match what the function has actually been doing, so results do not
change for anyone who left the argument alone, and the help says 1000 rather
than 100.
LRA(method = "GTM") ignored maxiter. The GTM branch called the EM loop
without passing the argument on, so it always stopped at emclus()'s default
of 100 cycles no matter what the caller asked for -- while the help said the
default was 1000. The SOM and isotonic branches were unaffected. Fits that
needed more than 100 cycles returned with converge = FALSE and estimates
short of the maximum; on the packaged datasets 100 was enough, which is why it
went unnoticed. Found while building a small cross-validation dataset that
needs 113 cycles, where the log-likelihood came out 7.2e-4 away from the
reference and lands on it exactly once the argument is honoured.
Posterior memberships were normalised by subtracting the row minimum, and
polytomous biclustering silently merged fields once the sample passed roughly
700 respondents. The E-steps turn a matrix of log-likelihoods into
memberships. Doing that safely means subtracting the row maximum before
exponentiating, so every exponent is at most zero. Several E-steps subtracted
the row minimum instead, which pushes every exponent positive, and then clipped
at exp(700) to stop the overflow. Everything above the clip lands on the same
value: two fields differing by a factor of exp(400) came out equally likely.
How far a row spreads depends on what is being summed over. A class posterior sums over items, so it stays small and the clip is rarely reached. A field posterior sums over respondents, so it grows with the sample -- which is why the visible symptom was backwards, more data giving a worse answer. Measured on ordinal data with 24 items in 3 true fields: fields recovered exactly up to 500 respondents, degrading at 700, and by 1000 two of the three fields had merged and one came out empty (adjusted Rand index 1.00 falling to 0.55). The class side was almost unaffected in the same runs, which is what made it look like a field-clustering quirk rather than a normalisation bug.
Affected: Biclustering() and Ranklustering() on ordinal and nominal data
(both the class and field steps), LCA() on nominal and rated data, and
LDB(). Binary biclustering was already correct. Fixed by routing all of them
through one row_softmax() helper.
The same helper replaces the unstabilised exp(ll) / rowSums(exp(ll)) in
binary LCA()/LRA(), LDLRA(), and LRA() on ordinal and rated data.
Those subtract nothing at all, so their exponents are negative and they
underflow rather than overflow -- reachable only on very long tests, and loud
(NaN) rather than silent when reached. The returned values are unchanged
wherever the old code did not overflow or underflow.
No test in the suite caught this: the biclustering reference fixtures top out at 515 respondents, below where the clip starts to bite. Regression tests over a range of sample sizes were added.
BINET() uses the same unstabilised form but is deliberately left alone: it
pins boundary-class respondents by writing a literal 1 next to the raw
likelihoods, so rescaling the rows would change how hard that pin holds. That
the strength of the pin depends on the test length looks unintended and is
being reviewed separately.
Order-restricted biclustering returned no field assignments at all on tests longer than about 40 items. The E-step needs the log of a category probability, formed as the difference of two upper-cumulative probabilities. That difference cannot be negative in exact arithmetic, but the order-restricted M-step pools adjacent categories, so exact ties are the norm rather than the exception (11 of 30 differences on a 72-item fit), and a tie computed as a subtraction lands anywhere within a few ulp of zero.
The guard against that was + const, and const is exp(-nitems). It
therefore shrinks as the test grows -- exp(-72) is 5e-32 against rounding
errors of 1e-14 -- while the errors it has to absorb do not. Past roughly 37
items the guard stops working: log() of a negative difference gave NaN,
the NaN reached the field posterior through the softmax, and max.col()
returned NA for every item, so the fit came back with all fields empty and
converge = FALSE. The cost was visible too, since a fit that never
stabilises runs to the iteration cap: 72 items in 6 fields now converges in 5
cycles and 0.2 seconds.
Fixed by clamping the difference at zero before adding const. Only the
rounding is repaired: a genuinely tied pair still becomes log(const) exactly
as before, and every already-positive value is unchanged to the bit, so no
previously-working fit moves. Ordinal LRA() is unaffected -- its isotonic
path works with the probabilities themselves rather than a difference.
const = exp(-nitems) is doing two jobs that pull apart as tests get longer:
a modelling penalty for an impossible response, which should scale with the
test, and a numerical floor, which should not. Separating them is left for a
later release.
M2() says how much memory it needs, and refuses a size the machine cannot
hold. The margin covariance is a dense m x m matrix and m grows with the
square of the test length: 35 items over 5 categories is 0.7 GB, 40 items is
1.2 GB, 72 items is 32 GB. Allocating tens of gigabytes unannounced is not
something a user can recover from -- the operating system kills the process,
taking the session and everything in it -- so M2() and add_M2() now report
the size before starting whenever the computation needs 2 GB or more, and stop
with an explicit error when the projected peak exceeds four fifths of physical
memory. No question is asked: a prompt would have to be limited to
interactive() to keep scripts and R CMD check from waiting on input, which
would make the behaviour depend on how the code was started.
options(exametrika.m2_max_gb = ) replaces the memory-derived ceiling with an
explicit one, and Inf removes the check for anyone who wants the old
behaviour back.
M2() on a rank-deficient margin covariance is roughly 300x faster. When
the Cholesky factorisation of Xi fails, the fallback was an
eigendecomposition -- an order of magnitude more arithmetic than the
factorisation it replaces, and it needs room for the full eigenvector matrix
beside Xi itself. It is now a pivoted Cholesky, which costs the same
m^3/3 and reports the rank directly. On the largest margin set met so far
(biclustering, 40 items x 5 categories, m = 12,640, Xi at 1.2 GB) that is
16 seconds against about 80 minutes. Both are generalised inverses of Xi
and the statistic only uses the residual's component in the column space, so
the value is unchanged; the equivalence is pinned by tests.
This mattered most for biclustering, where Xi is rank-deficient at any
useful size, so the slow path was the normal path rather than the exception.
The order-restricted solver no longer inherits the EM's iteration cap.
Biclustering() and LRA() on ordinal data passed their own maxiter
straight into iso_dual_map(), so raising the EM cap from 100 to 1000 for the
convergence fix silently raised the inner solver's cap tenfold as well. Two
unrelated loops shared one knob. The solver now uses its own budget. No
results change -- it converges well inside either cap, and the fit is
identical for every EM cap from 100 to 1000 -- but the coupling was not
intended.
maxiter documentation corrected for ordinal and nominal Biclustering(),
which said 100 after the default became 1000.
LRA() on ordinal or rated data could die on a lumpy score distribution.
The EM starting values come from grouping respondents by total score at the
score quantiles. When the scores pile up -- few categories with a strong floor
or ceiling is enough -- a quantile group comes out empty, its mean is 0/0, and
the NaN propagates until sort() drops it and an array assignment fails with
a message about replacement lengths that says nothing about scores. Measured
instance: 20 items, 3 categories, 6 ranks, 300 respondents, with 66 people at
the lowest possible score and 77 at the highest.
The grouping rule is unchanged -- it is the reference implementation's -- and
only the undefined values are repaired: an empty group borrows from the nearest
group that has respondents, which is meaningful because the groups are ordered
by score. Empty groups also had to survive table(), which drops unused levels
and left the rank-by-quantile cross-tabulation too small for its own dimnames.
M2() no longer fails on a singular margin covariance. Two causes, both
routine on real data. A category nobody chose gives an exact zero in the
profile, and the covariance correction then computes 0/0; the value that
division stands for is zero. A margin with a vanishingly small probability
leaves a direction with no usable variance, which a Cholesky factorisation
rejects and a check on the diagonal alone cannot catch, the dependence being
off-diagonal. M2() now falls back to an eigendecomposition and works in the
subspace with non-negligible eigenvalues, reporting how many directions were
dropped; the degrees of freedom count the subspace actually used.
EM convergence was declared far too early, and on short tests after a single cycle. This affects every EM-based model and is the largest behavioural change in this release: estimates move, and they no longer reproduce the numbers printed in Shojima (2022). Three separate defects, all inherited faithfully from the reference Mathematica implementation:
The starting value. test_log_lik was initialised to -1 / exp(-J),
i.e. -exp(J), meant to sit below any attainable log-likelihood. For a
short test with many respondents it does not: J = 6 gives -403 against a
real log-likelihood near -3900. The first cycle then looked like a
decrease, the "likelihood went down" guard fired, and the loop exited
after one iteration -- while reporting converge = TRUE. Measured:
binary LCA() on 6 items x 1000 respondents stopped at 1 cycle;
Biclustering() on the bundled J5S1000 stopped at 1 cycle. Tests never
caught it because every Mathematica-referenced dataset has J >= 15, where
exp(J) is safely below the log-likelihood. Now -Inf, with the
comparison skipped on the first pass.
The monitored quantity. The criterion watched the expected log-posterior
Q(theta_t | theta_{t-1}) (Shojima 2022, eqs. 5.11-5.12), not the
observed-data log-likelihood. EM guarantees the latter increases; the
former conditions on a moving point and need not. On 6-item data it does
decrease at the second cycle (-3869.26 -> -3922.42), firing the guard on a
healthy iteration. LCA() (binary/nominal/rated), LRA() (isotonic and
GTM) and LDLRA() now monitor the observed-data log-likelihood, computed
by log-sum-exp. LCA.nominal()'s reported log_lik is therefore now a
proper marginal likelihood, which also makes its AIC/BIC the textbook
quantities. Biclustering keeps its original monitored quantity
deliberately -- a two-sided mixture has no single obvious observed-data
likelihood -- so only its starting value and tolerance changed.
The tolerance. A relative 1e-4 on a value of order 1e3 is a threshold
of ~0.4 nats, so EM stopped well short of convergence. J15S500 with 5
classes ran 73 cycles to an observed log-likelihood of -4121.53, where a
properly converged run reaches -4118.16 in 440. The tolerance is now
1e-8 for the EM engines, matching what the nominal engines and GRM()
already used, and maxiter defaults rise from 100 to 1000 (largest
observed requirement across the bundled data: 685 cycles, 0.25s). IRT()
uses 1e-6: each of its cycles runs one optim() per item, so extra
cycles cost real time, and the attainable accuracy is bounded by
optim()'s own reltol anyway.
Effect on the published example (LCA(J15S500, ncls = 5)): 73 cycles -> 337,
and the class reference matrix of the book's Table 5.1 changes in the second
decimal (Item 1: .519 .700 .764 .856 .860 -> .590 .662 .772 .853 .884). The
Mathematica implementation in develop/mtmk15forVer13/ was corrected in step
and agrees with R to full double precision (-3651.903821773122 against
-3651.9038217731227); the reference fixtures under
tests/testthat/fixtures/mathematica_reference/ were regenerated from it.
IRT() is substantially faster. The M-step objective is the hot spot of
the whole fit (80% of self time). It evaluated the response function twice
per call where once suffices and used ifelse() on scalars; optim() was
also called without a gradient, so L-BFGS-B approximated one by finite
differences at a cost of one extra objective evaluation per parameter per
step. An analytic gradient is now supplied. IRT(J35S515, model = 2):
1.55s -> 0.34s at the new tolerance, faster than the 0.63s the old, looser
tolerance took. Estimates are unchanged to 5e-6.
Fixed Biclustering() on ordinal data (Biclustering.ordinal) under
Ranklustering (method = "R"): the GTM neighbour-smoothing filter was
computed (smoothed_memb <- clsmemb %*% Fil) but never used in estimation —
the field E-step and the M-step both consumed the raw, unsmoothed class
membership, so ordinal Ranklustering silently reduced to plain Biclustering
plus a cosmetic class relabelling by mean score. Both steps now use the
smoothed membership, matching Biclustering.binary and LRA.ordinal (for
plain Biclustering the filter is the identity, so results there are
unchanged). Fit improves slightly (e.g. RMSEA 0.0533 -> 0.0527, AIC
7288 -> 6720 on J35S500, ncls = 5, nfld = 5). No Mathematica reference
exists for ordinal Biclustering, so the affected regression-snapshot fit
indices in test-polytomous-biclustering.R were updated (R/16_Biclustering_ordinal.R).
Fixed the GRM item information function grm_iif(). The previous
implementation deviated from Samejima's (1969) information in four
ways: the denominator used the cumulative (boundary) probability
instead of the category probability, the accumulator started at 1
instead of 0, the first category's term was missing, and a 1.702
scaling constant was applied even though GRM() estimates its
parameters on the pure logistic metric (no scaling constant).
Together these inflated the reported information severalfold
(e.g. 1.93 vs the correct 0.70 at theta = 0 for a = 1.5,
b = c(-1, -0.5, 0.5, 1)). The function now computes
I(theta) = sum_k [P'_k(theta)]^2 / P_k(theta) on the same logistic
metric as the estimation, verified against numerical differentiation
and consistent with the posterior standard deviations
(1/PSD^2 ~ TIF + 1). Affects plot(x, type = "IIF") and
plot(x, type = "TIF") for GRM objects, and downstream ggExametrika.
Fixed the ability estimates of GRM() (EAP, MAP, PSD). The
posterior over the quadrature grid was built as the sum of the item
category probabilities (W %*% L_weighted) instead of their
product, which is not a valid likelihood. Rank order was nearly
preserved (r > 0.99 with the correct estimates on J5S1000), but the
scale was strongly shrunken (EAP sd 0.29 vs the correct 0.79) and
PSD did not correspond to the test information. The posterior is
now accumulated in log space with a log-sum-exp normalization, and
missing responses are skipped explicitly (previously they were
skipped only by an R indexing accident: v[NA] <- 1 silently does
nothing). ItemFitIndices/TestFitIndices change accordingly, since
the analysis-model log-likelihood is evaluated at each examinee's
EAP.
GRM plots (type = "IRF") no longer produce NA colors for items
with more than 8 response categories; the default palette is now
recycled.
The BNM/LDLRA structure-search functions (BNM_GA(), BNM_PBIL(),
LDLRA_PBIL()) no longer rebuild the full model output per candidate.
Fitness evaluation now goes through BIC-only internal kernels
(BNM_fit_BIC() for BNM, LD_param_est(fit_only = TRUE) for LDLRA) that
skip the output-table construction, the igraph round-trips and DAG checks
(GA/PBIL candidates live on the upper triangle, so they are always simple
DAGs), and the per-candidate benchmark-model recomputation inside
TestFit() — the benchmark statistics depend only on the data and are now
computed once per run (BNM_bench_stats()). Results are unchanged and
agree with the previous implementation exactly (identical(), same seeds:
BNM_GA all three crossover types, BNM_PBIL estimate 1/4, LDLRA_PBIL).
Measured: BNM_PBIL(J35S515, population = 20, maxGeneration = 15)
14.8s -> 0.7s (~22x); LDLRA_PBIL(J15S500, ncls = 3, population = 20, maxGeneration = 20) 3.8s -> 1.1s (a further ~3.5x on top of the earlier
PIRP vectorisation, ~14x cumulative from the original 15.3s).
BNM()'s parent-response pattern (PIRP) counting is vectorised: the
per-(subject, item) binary encoding and the subject x item x pattern
one-hot array are replaced by one matrix product and one rowsum() per
item (shared helper BNM_pirp_counts(), also used by the fitness kernel).
fill_adj()'s element-wise double loop is replaced by a single indexed
assignment. Same numbers, no S x J loop, and much less memory on large
samples. Profiling after the change shows no remaining hot spot in the
GA/PBIL path (the largest self-time item is rowsum()'s internal
bookkeeping at ~16%), so C++ is not warranted here; the package-level
next bottleneck remains the emclus() E-step.
BNM_GA(crossover = 2) no longer rebuilds the two-point crossover
candidate cut-point matrix (combn(gene_length, 2) filtered) per child
per generation; it is built once per run. RNG draw order is unchanged.
LRA() on ordinal and rated data no longer fails on large tables. The
monitored log-likelihood was accumulated as sum(rankProf * log(exp(ll))) --
a round trip that is algebraically the identity but underflows on the way:
once the test is large enough for the entries of ll to fall below about
-700 (2000 respondents x 60 items x 5 categories is enough), exp() returns
zero, log(0) is -Inf, and 0 * -Inf is NaN. The convergence check then
received NA and the fit died with "missing value where TRUE/FALSE needed".
The sum is now taken on the log scale directly, which is the same quantity
without the excursion through the exponent. Found while sizing a simulation
study; it affects LRA.ordinal() and LRA.rated() under both isotonic and
GTM.
The corresponding line in each saturated-model loop adds a small constant
inside the logarithm and is deliberately left alone: removing it moves the
saturated log-likelihood substantially (where exp(ll) < const the old
expression was pinned near log(const)), and the Mathematica reference values
were produced with that behaviour in place. Changing it is a separate,
breaking decision.
The order-restricted M-step for ordinal data is now implemented in C++
(src/isotonic_core.cpp). LRA.ordinal(method = "isotonic") and
Biclustering.ordinal(estimation = "isotonic") call it through the same
internal iso_dual_map() entry point. It exploits the fact that raising one
dual multiplier theta[c, q] can only change ranks c and c+1, so a trial
value rebuilds two rows rather than the whole table.
Both root finds use the residual's value rather than only its sign: Newton's
method for the multiplier that normalises a rank, and the Illinois method for
the dual variable. The pure-R iso_dual_map_ref() does the same and the two
are checked against each other to a tolerance -- not to the bit, which a
value-driven search cannot promise: a difference in the last bits moves the
next trial point and the two paths then differ on their way to the same root.
For the multiplier, the parameterisation matters more than the root finder.
Categories with a zero count are dropped first: their probability is zero
whatever the multiplier is, but if one of them attains the smallest offset it
drags the lower end of the domain far away from the root. On the shifted
variable u = lambda + min(d) the bracket then follows in closed form,
sum_{d'=0} m <= u <= sum m, and the search runs on log u, because those two
ends can sit thirteen orders of magnitude apart on a rank holding almost no
weight -- far enough that bisection does not converge within any sane
iteration cap and a Newton step from the upper end lands past zero. The
probabilities are formed from u + d' rather than lambda + d, which would
cancel the shift back out and lose digits when min(d) is large.
The Illinois search judges its bracket width relatively. On the same low-weight ranks an endpoint can reach 1e8, where the spacing between representable doubles (~3e-8) is coarser than an absolute 1e-12 threshold: the interval stops shrinking, every interpolated point rounds onto an endpoint, and the loop has no way to exit.
Each sweep also starts its bracket from the value the previous sweep found, which is close to the answer from the second sweep on; the check of whether the constraint binds at zero is kept, being complementary slackness itself.
The sweeps stop on both halves of the KKT picture, not on the log-likelihood
alone. Near the optimum the likelihood moves quadratically while the order
violation is still shrinking geometrically, so the likelihood goes quiet
several sweeps before the restriction is actually met: on real tables it was
stopping with the returned probabilities breaching the stochastic order by
around 1e-3. Convergence now also requires the largest violation to fall
below viol_tol (1e-6), and the sweep cap rises to 1000 because the number
of sweeps grows with the table -- 6 ranks need about 100, 20 ranks about 500.
Across 60 tables from real EM runs the worst violation drops from 3.2e-3 to
9.9e-7 for a 1.2x cost, and Biclustering(J35S500, ncls = 6, nfld = 5, method = "R", estimation = "isotonic") is unchanged at 2.7s.
Measured on the dual solver alone: 111x (3 ranks x 3 categories), 289x
(12 x 7), 509x (20 x 6; 274s -> 0.54s), and a further 6.8x from the two root
finds and the warm start, on tables taken from a real EM run. At the model level, where the
remaining EM machinery is still R, Biclustering(J35S500, ncls = 6, nfld = 5, method = "R", estimation = "isotonic") drops from 265s to 11s.
The multiplier exists to make each rank's probabilities sum to one, so the row sums measure how well it is found. Across 60 expected-count tables captured from real EM runs the largest row-sum error is 9.6e-15, and every table converges. For binary data the solution must coincide with weighted PAVA, and it does.
This matters because the ordinal isotonic path was the slowest thing in the
package and effectively ruled out large simulation studies. The pure-R
implementation is kept as iso_dual_map_ref() (internal, unexported) purely
so the test suite can check the two against each other and so the algorithm
stays readable.
Spell-checking became usable again: tools/spell_check.R skips the
Japanese vignette (morphemes are permanent hunspell noise) and NEWS.md
(history is not edited), inst/WORDLIST was brought up to date (136 ->
228 words), and release_bullets() reminds every release to keep it so.
Zero findings is now the normal state, so the next finding is a real typo.
Restructured the test suite around a two-tier policy: CRAN runs tiny
Mathematica-verified fixtures plus a smoke test per user-facing function,
while full-size reference comparisons, reproducibility runs and heavy
regressions carry skip_on_cran() and run on every push in CI
(GitHub Actions, three platforms, NOT_CRAN=true pinned) and locally.
The CRAN-side suite now takes ~17 s sequentially on an M-series laptop
(was 38 s), with no test file above 1.5 s, keeping the Windows check
comfortably inside CRAN's 10-minute limit. Enabled
Config/testthat/parallel: true. No user-facing behaviour changes.
Renamed the internal helper iso_surv() in R/00_isotonic_CORE.R to
iso_upper_cum(). It computes the upper-cumulative (boundary) probabilities
S_cb = P(>= category b+1); the old name suggested a survival function,
which is not the concept being used here.
New: M2(), the limited-information goodness-of-fit statistic of
Maydeu-Olivares & Joe (2005, 2006), for LCA.nominal and LCA.rated. It tests
whether the model reproduces the item-pair cross tables. The reference point is
the saturated model of the first- and second-order margins rather than of the
full response-pattern table, so no benchmark log-likelihood is needed — which is
exactly what nominal data lacks. The analogy that usually lands: an SEM
chi-square does not test the whole multivariate distribution, only whether the
model reproduces the covariance matrix; M2() is the categorical counterpart,
with cross tables in place of the covariance matrix.
Two things worth knowing. The degrees of freedom are m - rank(Delta), not
m - n_param: with the class proportions fixed at 1/C, the second-order
margins see the class deviations only through their Gram matrix, which is
invariant to rotations of the class space, so the Jacobian loses
(ncls - 1)(ncls - 2) / 2 of its rank and the parameters are not identified
from bivariate margins once ncls >= 3 (the margins themselves still are, so
the statistic is well defined). And the statistic is only interpretable for a
maximum likelihood fit — not for filter-based (GTM) estimation, which is a
regularisation rather than an MLE, nor for order-restricted estimation, whose
limiting distribution is a mixture of chi-squares.
Cost is dominated by the Cholesky factorisation of a dense m x m matrix:
20 items with 5 categories is m = 3120 and well under a second; 50 items is
m = 19800, 2.9 GB and around 20 seconds.
LCA.nominal() and LCA.rated() now keep Q and Z in the returned object,
as Biclustering() already did, so M2() can be computed after the fact.
New: add_M2() and a second set of fit indices. fit <- add_M2(fit)
computes M2 for the model and for the independence baseline and attaches
TestFitIndicesM2: NFI, RFI, IFI, TLI, CFI and RMSEA built entirely from
margin-based chi-squares. print() then shows the response-pattern indices and
the margin-based ones as two labelled blocks, and takes fit_indices = "both" | "pattern" | "margin" to show one alone.
They are kept apart on purpose. The response-pattern chi-square is a likelihood
ratio against the saturated model; M2 is a quadratic form in margin
residuals. Combining one with the other inside a single index — a ratio of a
pattern chi-square to a margin chi-square, say — has no defensible reading
(Shojima, personal communication, 2026-07-26). Within the margin block both
the model and the baseline statistic are M2, so the indices are internally
consistent.
For the same reason the margin block carries no AIC/BIC/CAIC. Those need a
log-likelihood: chi_sq - 2 df works on the pattern side only because that
chi-square is a likelihood ratio, and M2 is not. The likelihood-based
information criteria these models already report are the ones to use.
This is a separate call rather than part of the fit because the cost grows with the square of the item count; for nominal data it is also the only route to an absolute fit statistic at all, since no saturated model can be formed.
M2() also accepts fits from LRA() on ordinal data and from
Biclustering() on ordinal or nominal data, so the three ways of treating the
same polytomous data — order-restricted, filter-smoothed, unordered — can be
put on one scale. Biclustering shares one category profile across the items of
a field, so its Jacobian columns are shared too, and its parameter count is
fields x classes x (categories - 1) rather than items x classes x (categories
Only a maximum likelihood fit makes the reference distribution hold, so for filter-based (GTM) and order-restricted estimation the statistic is descriptive: recorded, not tested.
add_M2() covers LRA() on ordinal data and Biclustering() on ordinal or
nominal data as well, so any of them can carry both sets of fit indices. For
ordinal data the response-pattern indices are complete — a saturated model
exists — so two full sets stand side by side, which is precisely where keeping
them apart matters.
Where the reference distribution does not hold, the margin block says so, and names the reason that applies rather than issuing a generic warning: the field partition is taken as given (any biclustering), the order restriction can bind so the limit is a mixture of chi-squares (isotonic), filter smoothing is a regularisation rather than a maximum likelihood estimator (GTM). A nominal LCA gets no such line, being an honest maximum likelihood fit. The numbers are still worth reading as description; the p value is not a test.
M2() and add_M2() take gc = TRUE, which releases the workspace before
returning. The margin covariance and its factor are the largest objects the
package allocates and R otherwise keeps the block. Interactive use wants this;
a loop over many fits should pass FALSE, where the collection costs time and
buys nothing. With verbose = TRUE and a large problem, add_M2() now says
which of its two statistics it is working on.
New: LCA() supports nominal data (LCA.nominal). The model is a finite
mixture of product-multinomial distributions — one free category distribution
per item within each latent class, with no ordering imposed on the classes or
on the categories. Category counts may differ across items. Estimation runs
through the shared EM engine emclus_nominal(); the new alpha argument sets
a Dirichlet prior on the category profiles (default 1, i.e. the plain
multinomial MLE). The returned object carries ICRP (Item Category Reference
Profile: one row per item-category, one column per class), LCD, CMD and
Students, and prints and plots (type = "ICRP", "LCD", "CMP"). The
"ICRP" plot uses grouped bars rather than lines across classes, because
neither axis is ordered.
As in Biclustering.nominal, no benchmark (saturated) model is fitted — with
many items and categories nearly every response pattern is unique — so only
AIC, BIC and CAIC are reported and the chi-square based indices are NA. Note
that these follow the -2 log L + k convention, unlike the chi-square based
AIC/BIC/CAIC on the binary path; they are not comparable across response
types. Class sizes are not free parameters in this formulation (every class
carries the same implicit prior, as in emclus() for binary data), so the
parameter count is ncls * sum(ncat - 1).
New: LCA() supports rated data (LCA.rated), i.e. multiple-choice items
with an answer key. Estimation is the nominal one — LCA.nominal() is called
internally — and the key is then used to recover the quantities that need a
correct answer: IRP[j, c] is the model-implied probability of the keyed
category rho[j, CA[j] | c], and TRP is its weighted item sum. Fit is
reported in two layers: TestFitIndices from correct/incorrect responses
(chi-square based, comparable with binary LCA) and TestFitIndicesNominal
from the internal nominal fit (AIC/BIC/CAIC). Full category probabilities
remain in ICRP for distractor analysis. Unlike Biclustering.rated, the
classes are not sorted by correct response rate: latent classes carry no
order, and sorting them would imply one.
Ordered rating data passed to LCA() is now analysed with the nominal model
after a message, instead of stopping: latent classes carry no order, so the
category ordering has nothing to attach to. Use LRA() when the ordering
should be respected. Symmetrically, nominal data passed to LRA() is handed
to LCA.nominal() with a message (a rank ordering is not defined without
ordered categories); an nrank argument is carried over as ncls.
LCA() is now an S3 generic dispatching on response type, matching
LRA() and Biclustering(). The former function body is unchanged and
now lives in LCA.binary(); LCA.default() formats raw input with
dataFormat() and routes it. Binary results are bit-identical to 1.15.0 and
non-binary input still raises the same error. This is the wiring for the
forthcoming LCA.nominal().
Biclustering() on binary data (Biclustering.binary) gains an
order-restricted estimation of the Field Reference Profiles under
Ranklustering (method = "R"), selected by the new estimation argument:
"isotonic" (now the default) or "GTM" (the original filter-based
smoothing of Shojima 2012). Under "isotonic" the rank ordering — each
field's reference profile monotone non-decreasing across ranks — is imposed
directly in the M-step by a weighted PAVA along the class axis (weights are
the per-cell expected counts), reusing pava_up() from
R/00_isotonic_CORE.R; for the default flat prior this is the exact
order-restricted MLE (Ayer et al. 1955), and the crude mic sort is skipped
as unnecessary. On typical data it attains a higher likelihood than GTM
(e.g. -6941 vs -7273 on J35S515, ncls = 6, nfld = 5). The
shape-restricted degrees of freedom (total PAVA block count across fields)
feed the fit indices. estimation is ignored for plain Biclustering
(method = "B"), whose classes are unordered (R/07_Biclustering.R).
Biclustering() on ordinal data (Biclustering.ordinal) gains the same
estimation argument for Ranklustering (method = "R"): "isotonic" (now
the default) or "GTM". Under "isotonic" the rank ordering is imposed
as a genuine stochastic-order restriction — each field's upper-cumulative
(boundary) probabilities are monotone non-decreasing across ranks at every
category threshold — solved per field in the M-step by the Fenchel-dual
algorithm (El Barmi & Dykstra 1994) via iso_dual_map() in
R/00_isotonic_CORE.R, replacing the independent per-cell boundary MLE and
the crude mic class relabelling. Unlike the binary case, isotonic does not
maximise the likelihood relative to GTM: the (now filter-smoothed) GTM path
leaves the field profiles unconstrained and so fits better by violating the
order (e.g. only 1 of 5 fields is monotone on J35S500), whereas isotonic
enforces the order in every field at a modest likelihood cost — this is the
intended trade and the point of the method. Shape-restricted df = distinct
free boundary levels per field. estimation is ignored for plain
Biclustering (method = "B") (R/16_Biclustering_ordinal.R).
LRA() on binary data (LRA.binary) gains an order-restricted estimation
method, method = "isotonic", which is now the default (previously
"GTM"). Rather than the GTM filter, the rank ordering — item reference
profiles monotone non-decreasing across ranks — is imposed directly in the
M-step by a weighted PAVA down each item column; for the default flat prior
(beta1 = beta2 = 1) this is the exact order-restricted MLE (Ayer et al.
1955). On typical data it attains a higher likelihood than GTM/SOM. The
shape-restricted degrees of freedom (per-item PAVA block count) feed the
item/test fit indices. method = "GTM" and method = "SOM" are unchanged
and still selectable. New shared internal core R/00_isotonic_CORE.R
(pava_up(), emclus_isotonic()); ItemFit() now accepts a per-item
nparam vector in addition to a scalar.
LRA() on ordinal data (LRA.ordinal) gains an order-restricted estimation
method, method = "isotonic", now the default (previously the GTM
filter). The stochastic-order restriction (item category boundaries monotone
across ranks) is solved exactly in the M-step by the Fenchel-dual algorithm
(El Barmi & Dykstra 1994) in R/00_isotonic_CORE.R, rather than by filter
smoothing; a alpha argument sets the Dirichlet concentration (the
polytomous analogue of binary beta1/beta2, default 1 = ML). It attains a
higher likelihood than GTM on typical data, handles mixed category counts
per item, and reports a shape-restricted degrees of freedom (distinct
boundary levels per item). method = "GTM" is unchanged and still selectable.
LRA() on ordinal data (LRA.ordinal) now supports items with
differing numbers of response categories. Previously the ordinal
method assumed every item had the same number of categories and
raised an error otherwise (directing users to Biclustering.ordinal),
because its saturation/restricted reference matrices were laid out on
a fixed nitems * max(ncat) stride. The reference matrices, the
one-hot design (uuMat), and the null-model term now use a ragged
sum(ncat) layout addressed through the design1 index map,
so mixed category counts are handled directly. On uniform-category
data the results are numerically identical to before. ICBR/ICRP
now carry one row per (item, category) pair, so their row counts vary
by item under mixed categories (R/12_LRA_ordinal.R).
dataFormat()'s id argument, and longdataFormat()'s Sid/Qid/Resp/w
arguments, now accept a column name (character string) in addition to a
column number. Previously only numeric column indices were accepted;
dataFormat() even raised an explicit error when a character value was
passed. Column names can now be supplied directly, e.g.
dataFormat(data, id = "StudentID") or
longdataFormat(data, Sid = "Sid", Qid = "Qid", Resp = "Resp"). Specifying
an unknown or ambiguous (duplicated) column name now raises a clear error
listing the available columns (R/01_dataFormat.R).An API audit across all exported model functions surfaced a number of gratuitous inconsistencies in argument names, argument order, and defaults. They are unified in this release; the changes below alter default behavior or (rarely) positional-argument calls, so please review if you relied on the old defaults.
The verbose argument now defaults to FALSE in every function that has
one. Previously IRT(), GRM(), LCA(), GridSearch(), all
Biclustering() methods, all Biclustering_IRM() methods, BNM_GA(),
and BNM_PBIL() defaulted to TRUE while LRA(), LDLRA(), LDB(),
BINET(), LDLRA_PBIL(), Glasso(), and chatterjee_matrix() defaulted
to FALSE. Estimation results are unaffected; pass verbose = TRUE to
restore progress messages. Several roxygen pages that claimed
"default is TRUE" for functions whose default was already FALSE were
corrected as well.
Biclustering_IRM()'s iteration-limit argument is renamed from
max_iter to maxiter, matching Biclustering(), LCA(), LRA(), and
the rest of the package. Because the IRM methods accept ..., a
maxiter = argument passed under the old assumption was silently
swallowed with no effect — an easy trap when switching between
Biclustering() and Biclustering_IRM(). The old name is still honored
(with a deprecation warning) and will be removed in v2.0.0.
Biclustering_IRM() for ordinal data now defaults to mic = FALSE and
EM_limit = 20, matching LRA()/Biclustering() (ordinal) and the
other Biclustering_IRM() methods respectively. It previously defaulted
to mic = TRUE and EM_limit = 100 with no documented rationale. Pass
mic = TRUE explicitly if you want monotonically increasing Field
Reference Profiles.
BNM(), LDLRA(), LDB(), BINET(), BNM_GA(), BNM_PBIL(),
LDLRA_PBIL(), and Biclustering_IRM() (binary) now declare their
missing-data arguments in the order na, Z, w — the order used by
CTT(), IRT(), GRM(), LCA(), LRA(), and Biclustering() — where
they previously used Z, w, na. Only calls that passed these three
arguments positionally are affected; named usage (Z = ...) is
unchanged.
Biclustering() for rated data now exposes conf_class in its signature
and documentation. It was already forwarded to the internal nominal
engine via ..., but never appeared on the help page. Also documented
that the rated method defaults to method = "R" (Ranklustering, classes
sorted by correct response rate) while binary/ordinal default to
method = "B".
LDLRA()'s beta1/beta2 defaults remain 2 (the other network models
use 1); the help page now notes this follows the original Mathematica
implementation of LDLRA rather than being an accident.
longdataFormat() falsely reported "Duplicated IDs found" whenever a
student answered more than one item. The duplicate check ran on the raw
student-ID column, which is expected to repeat once per item in long
format; it now checks for duplicated (student, item) pairs instead, so
ordinary long-format data (one row per student-item response) no longer
triggers a spurious error. A genuine duplicate — the same student
answering the same item twice — is still caught (R/01_dataFormat.R).
IRT(), Biclustering(), BNM(), LDLRA(), LDB(), and BINET() silently
treated missing responses as incorrect when computing item-total correlation
or correct response rate (crr()). In each case, a response matrix that
had already been stripped of its exametrika class (with missing values
recoded to 0 via tmp$U * tmp$Z) was passed to ItemTotalCorr(),
ItemThreshold(), or crr(). Those functions re-run dataFormat()
internally when given an unclassed matrix, but without the original missing
-value mask, so missing responses were counted as incorrect answers. This
affected IRT()'s initial parameter values (rho/tau, which only seed
the EM algorithm — final estimates were unaffected), Biclustering()'s
FieldAnalysis$CRR column, and the $crr field returned by BNM(),
LDLRA(), LDB(), and BINET() (including their BNM_GA()/BNM_PBIL()/
LDLRA_PBIL() wrappers). Fixed by passing the already-formatted object
through instead of the stripped matrix (R/04C_ParameterEstimation.R,
R/07_Biclustering.R, R/08A_BNM.R, R/09_LDLRA.R, R/10_LDB.R,
R/11_BINET.R).
CTT(), BNM(), LDLRA(), LDB(), and BINET() crashed when passed raw
(unformatted) matrix/data.frame input instead of a pre-built
exametrika object, even though this is documented as supported. The
response-type check in each function read U$response.type from the raw
input argument instead of from the dataFormat()-formatted object, which
is NULL for raw input and made the check error out before model fitting
ever started. CTT() had a more severe variant of the same bug: its
inherits(U, "exametrika") branch was inverted, so raw input skipped
dataFormat() entirely and crashed immediately. Fixed to consistently
check the formatted object (R/03_CTT.R, R/08A_BNM.R, R/09_LDLRA.R,
R/10_LDB.R, R/11_BINET.R).
BINET() could crash with a cryptic 'dimnames' length error (via NaN
values silently produced by log()) on data with missing responses. Two
compounding issues were found. First, Ccj (the per-class correct-response
count used throughout the RISP estimation) was computed as
t(clsmemb) %*% tmp$U, which does not mask out missing cells (encoded as
-1 in tmp$U) the way the corresponding incorrect-response count Fcj
already did (t(clsmemb) %*% (tmp$Z * (1 - tmp$U))); missing responses
therefore contributed -1 to the correct-response count instead of 0,
corrupting the conditional correct-response-rate estimates. Second, the
smoothing constant that keeps those conditional rates away from 0/0 when
a class-by-field cell has no non-missing observations was hardcoded to
gamp <- 1 (a flat Beta(1,1) prior, under which 0/0 is undefined) with
no way for the caller to change it, unlike the equivalent beta1/beta2
arguments already exposed by BNM(), LDLRA(), and LDB(). Fixed the
Ccj computation to mask missing cells via tmp$Z, added beta1 = 1 /
beta2 = 1 arguments to BINET() so the prior can be strengthened when
needed, and replaced the eventual downstream crash with an immediate,
actionable error when a class-by-field cell is still undefined under the
requested smoothing (R/11_BINET.R).
0-indexed or non-contiguous polytomous category codes were silently
dropped in Biclustering()/Biclustering_IRM() for nominal and ordinal
data, with no error or warning. Biclustering.nominal(),
Biclustering.ordinal(), Biclustering_IRM.nominal(), and
Biclustering_IRM.ordinal() (and, via delegation, the rated/rated_IRM
variants) built a one-hot response array by indexing directly with each
item's raw category code; a raw code of 0 is a no-op under R's
array-indexing rules, so every response in category 0 vanished from the
likelihood and classification with no diagnostic. On the package's own
J15S3810 sample data (categories coded 0..3), this silently discarded
3,854 observed responses. GRM() already remapped each item's categories
to contiguous 1-based indices to avoid exactly this; that remap is now
shared (remap_category_codes(), R/00_BiclucterUtils.R) and applied in
all four functions.
Biclustering_IRM.nominal()/Biclustering_IRM.ordinal() corrupted
results whenever the data had missing responses. The one-hot response
array was built by writing directly at Uq[s, j, tmp$Q[s, j]] <- 1
without checking for missing cells; a missing cell's -1 sentinel, used
as a raw array index, sets every category except the first to 1 under
R's negative-indexing rules instead of leaving the cell at 0. This fed
through into the Gibbs sampler's small-class-adjustment step and the final
class/field assignments. With 8% missing data injected into J20S600,
7% of final class assignments differed from the correct (missing-masked)
result. Fixed by building the one-hot array with the same
missing-value-masked construction already used in the binary/nominal/
ordinal Biclustering() path.
LDB() crashed with "object 'conf_mat' not found" when conf was
passed as a matrix or data.frame (only the vector form of conf worked).
This is the same class of bug already fixed in Biclustering.ordinal()/
Biclustering.nominal() on 2026-04-27, which LDB() was missed by at the
time (R/10_LDB.R).
LDB()'s beta1/beta2 smoothing arguments had their roles swapped
relative to BNM()/LDLRA()/BINET(). All four functions compute a
posterior-mode correct-response rate of the form
(count + α - 1) / (N + α + β - 2) under a Beta(α, β) prior; BNM(),
LD_param_est() (used by LDLRA()), and BINET() consistently use
beta1 as α (the "success" pseudo-count, matching the count in the
numerator), but LDB() used beta2. At the package default
beta1 = beta2 = 1 this made no numerical difference, but a caller
following BINET()'s own advice to set an asymmetric prior (e.g.
beta1 = 2, beta2 = 2 to avoid a degenerate 0/0 cell) would see the
opposite effect in LDB() versus its sibling functions. Traced to the
original Mathematica Chapter 10 module
(develop/mtmk15forVer13/mod/Module_LDB.nb), which uses beta2 here —
this is not an R porting bug but an inconsistency inherited from the
reference implementation. Changed LDB() to use beta1, matching the
mathematically conventional Beta(α, β) parameterization and the rest of
the package (R/10_LDB.R).
LDLRA(..., method = "C") (local dependence latent Class model)
mislabeled its output as the Rank model. A hardcoded model <- 2
right before building the output tables silently overrode the model
value that had correctly been set from method earlier in the function,
so result$model, result$msg, and the "Rank"/"Class" column headers in
Estimation_table/CCRR_table always read "Rank" regardless of the
requested method. Numerical estimates were unaffected — only the labeling
(R/09_LDLRA.R).
BNM()'s and BINET()'s graph-validity checks (acyclicFLG,
connectedFLG) never actually detected cycles or disconnection. Both
used adj^i in a loop intending to count length-i walks, but ^ on a
matrix in R is elementwise, not matrix power, so the loop only ever
re-tested self-loops and "any edge present." BINET()'s
stop("Your graph is not a DAG") guard was consequently a no-op for any
real cycle of length ≥ 2, and BNM() silently fit cyclic "Bayesian
networks" with no warning. (BINET()'s connectedFLG was further
computing off the student response matrix, not a graph object, so it
could not have tested connectivity at all.) Replaced both with
igraph::is_dag() / igraph::is_connected(mode = "weak") on the actual
adjacency matrix being validated (R/08A_BNM.R, R/11_BINET.R).
Biclustering_IRM.binary()'s small-class-reduction loop tracked
bestclass/class, which resolve to base R's class() function rather
than any local variable, so the intended "roll back to the best-BIC
class assignment" step silently did nothing. In practice this had no
effect on the final returned model (the class/field assignments are
recomputed from scratch after the loop regardless), but the code read as
doing something it did not; fixed to track the actual cls vector
(R/07_IRM.R).
CCRR() returned JCRR()'s result when called on non-exametrika
input — a copy-paste bug in CCRR.default()'s fallback branch called
JCRR(U) instead of recursing into CCRR(U) (R/02_TestItemFunctions.R).
MutualInformation()'s base argument was ignored for binary data;
MutualInformation.binary() always computed log base 2 regardless of the
requested base, unlike the ordinal/nominal/rated path which already
respected it (R/02_TestItemFunctions.R).
JSR() and CSR() crashed instead of falling back cleanly when given
binary data. Both message "using X instead" and call the binary
equivalent, but never return()ed its result, so execution fell through
into polytomous-only code (NCOL(U$Q), U$Z == 0 on nonexistent U$Q)
and crashed (R/02_TestItemFunctions.R).
ItemTotalCorr.ordinal() did not mask missing responses before summing
total scores. total <- rowSums(U$Q) included the -1 missing
sentinel directly in the sum, so any student with even one missing item
got a wrong total score, corrupting their item-total correlation with
every item, not just the missing one. Fixed to mask by Z before
summing, matching ItemThreshold.ordinal()'s existing pattern
(R/02_TestItemFunctions.R).
ITBiserial(), ScoreReport(), and ItemReport() crashed when passed
raw (unformatted) data — each called dataFormat() but stored the
result in a variable that was never used, leaving the original
unformatted argument in place for the subsequent $response.type check
and computation (R/02_TestItemFunctions.R, R/02_QitemFunctions.R).
dataFormat() silently skipped value validation for ordinal/rated/
nominal data. The branch was guarded by
response.type == "polytgomous" (a typo — no response.type value is
ever literally "polytomous", since the actual values are "ordinal"/
"rated"/"nominal"), so the "must be non-negative integers" check
never ran for any polytomous type. Fixed to check
response.type %in% c("ordinal", "rated") (deliberately excluding
"nominal", whose category codes are arbitrary labels rather than
ordered counts and may legitimately be negative) (R/01_dataFormat.R).
longdataFormat() misaligned ID/response-matrix rows whenever
numeric student or item IDs were not an exact 1..N sequence (e.g.
student IDs 10, 20, 30). Numeric IDs were used directly as row/column
indices instead of being remapped to a dense index, so nrow(result$U)
could be far larger than the actual number of students, with most rows
entirely empty. Both Sid/Qid are now always remapped through
as.factor() to a dense, correctly-ordered 1-based index, regardless of
whether the raw IDs are numeric or character (R/01_dataFormat.R).
longdataFormat()'s w (item weight) argument was misused as a row
index. w_vec[unique(Qid_num)] indexed the long-format weight column
by item-ID values rather than by the row where each item first
appears, so weights were essentially always wrong unless item IDs
happened to be an exact 1..N sequence starting from the first item
column. Fixed to look up each item's weight via its first occurrence in
the long-format data (R/01_dataFormat.R).
print(GridSearch(...)) always crashed. print.exametrika()'s
internal switch() had no case for the "GridSearch" class, and a
trailing comma after the last real case created an empty unnamed
fallback branch that errors when evaluated. GridSearch() is a
top-level function whose result is auto-printed at the console in normal
interactive use, so this reliably crashed on first use. Added a proper
GridSearch print case summarizing the fit-index grid and optimal
setting (R/00_exametrikaPrint.R).
plot(x) without an explicit type argument crashed with a confusing
"the condition has length > 1" error instead of the intended
"The 'type' argument must be specified." message. The missing(type)
check ran after type had already been evaluated (forcing its
multi-element default) for an unrelated uses_layout computation, so by
the time missing() was checked, type was no longer missing in the
technical sense needed for the check to matter — and the multi-element
default itself broke the preceding if first. Reordered so the
missing(type) check runs first (R/00_exametrikaPlot.R).
plot(grm_result, type = "ICC") failed even though the identical call
works for IRT models. GRM's plot-type alias handling for "ICC" was
implemented but never reachable, because valid_types$GRM in the
upstream validation list omitted "ICC" (only "IIC"/"TIC" were
listed). Added "ICC" to valid_types$GRM (R/00_exametrikaPlot.R).
print(bnm_result, digits = ...)/print(ldlra_result, digits = ...)
ignored the requested digits for the "Conditional Correct Response
Rate" table. A typo (x$CCRR_tabje instead of x$CCRR_table, in two
places) wrote the formatted string into a nonexistent list field instead
of the one actually printed immediately afterward. Also propagated the
existing digits argument to two related tables (IRPIndex/FRPIndex
in the LDLRA/LDB print sections) that were missing it while equivalent
tables elsewhere already had it (R/00_exametrikaPrint.R).
print(InterItemAnalysis(ordinal_data)) never displayed the
Conditional Selection Ratio — it printed x$JSR (Joint Selection
Ratio, already shown just above) a second time instead of x$CSR
(R/00_exametrikaPrint.R).
GridSearch() could report the wrong optimal (ncls, nfld) when
multiple grid cells tied for the best fit index. which(ret == ..., arr.ind = TRUE) returns one row per tied match; optimal_idx[1]/
optimal_idx[2] read the matrix's column-major-flattened first two
values instead of the first tie's row/col pair, so with ≥ 2 ties the
reported optimum could combine the row of one tie with the column of
another — a combination that was not actually the best (or even
necessarily tied). Fixed to read optimal_idx[1, "row"]/
optimal_idx[1, "col"] (R/00_GridSearch.R).
GridSearch(..., max_ncls = 1) silently tested ncls = 2 anyway.
2:max_ncls is R's descending-sequence trap: with max_ncls = 1 it
evaluates to c(2, 1) rather than an empty range, so a value the caller
explicitly excluded was tested regardless. Added an explicit
max_ncls/max_nfld < 2 guard that stops with a clear message instead
(R/00_GridSearch.R).
Glasso()'s edge_tol argument was documented as controlling edge
detection but did not affect model selection. The EBIC computation
used internally to choose lambda_opt had its own hardcoded 1e-6
edge-count threshold, so changing edge_tol only affected the
reported n_edge, not which lambda was actually selected. Threaded
edge_tol through to compute_EBIC_glasso() so it consistently governs
both (R/22_GlassoUnit.R).
dataFormat() could silently miss items with no real variance whenever
the item also had missing responses, letting them through into IRT(),
which then crashed with the uninformative "object 'result' not found"
(the EM optimizer's initial parameter, seeded from an item-total
correlation of NaN, made every retry of optim() fail). The existing
zero-variance check computed sd() on the raw response column, which
still contains the -1 missing-value sentinel; an item that only a
handful of students answered, all identically, therefore looked
non-constant (-1 mixed with the real value) even though it truly had
zero variance among actual respondents. sd.check is now computed on
each item's valid (Z == 1) responses only. Additionally, items with no
variance among valid responses are now excluded from the data (with a
message identifying them) instead of merely warned about and left in,
matching the original Mathematica dataformat[], which drops such items
before any analysis; an item with CA specified (rated data) has its
CA entry dropped along with it so the two stay aligned. All-missing
items (no valid responses at all) are unaffected and continue to be kept
with a separate warning (R/01_dataFormat.R).
DistractorAnalysis() counted raw category codes against columns
1..maxQ of its frequency tables. For rated data coded from 0 (or with
gaps in the codes), every category-0 response silently vanished from the
frequency/proportion tables and the remaining categories were shifted by
one column — inconsistent with the fitted model itself, whose category
probabilities are computed on codes remapped to 1..K per item since the
category-code fix earlier in this release. The stored Q and correct
answers are now remapped the same way before tabulation
(R/21_DistractorAnalysis.R).
A conf given as a data.frame was wrongly rejected ("The conf matrix
should only contain 0s and 1s") by Biclustering() and LDB(), because
the 0/1 validation ran before coercion and %in% compares whole columns
on a data.frame. The input is now coerced with as.matrix() first
(R/07_Biclustering.R).
GridSearch()'s help-page example called a nonexistent function
grid_serch(data_matrix, ...); it now calls GridSearch(J35S515, ...)
(R/00_GridSearch.R).
LDLRA() now stops with an informative error when a
rank-by-parent-pattern cell has no (smoothed) observations under
beta1 = beta2 = 1 (0/0), instead of silently propagating NaN through
the log-likelihood, matching the diagnostic added to BINET() earlier in
this release (R/09_LDLRA.R).
LRA() for rated data has a minFreqRatio argument that is documented to
pool response categories rarer than nobs * minFreqRatio into a single
bucket, but the feature was doubly broken: the collapse test compared each
category's frequency against the correct-answer code (catfreq != CA,
almost always true), so no category was ever pooled; and had it fired, the
report assembly indexed the full category-label set and would have stopped
on a length mismatch. The condition now reads "keep a category if it is
frequent enough or it is the correct answer" (the correct answer must
never be pooled), and the Item-Category Reference Profile labels the
pooled bucket "CatX", using only the labels of the categories actually
estimated. The default minFreqRatio = 0 keeps every category and is
unaffected (R/13_LRA_rated.R).
Fixed user-visible typos in printed output: "Dimensionality Analyeis" →
"Dimensionality Analysis", "Cummurative Percentage" → "Cumulative
Percentage", and "Conditonal" → "Conditional" (in "Conditional Correct
Response Ratio" and "Conditional Selection Ratio" headings)
(R/00_print_ctt_irt.R).
LLtheta_mat()/EAP_PSD() (R/04B_AbilityEstimation.R) and
the corresponding dead print.exametrika() "IRT_EAP_PSD" case.
Neither function was exported or called anywhere in the package; IRT()
computes ability EAP/PSD inline instead. Also removed a handful of other
now-dead assignments left over from the same tmp$U * tmp$Z
class-stripping pattern fixed above, and one duplicated-then-discarded
array/loop in LDLRA()'s internal LD_param_est(), none of which were
ever read after being written (R/06_LRA.R, R/08C_BNM_GA.R,
R/09B_LDLRA_GA.R, R/09_LDLRA.R, R/11_BINET.R).Follow-up cleanup after the audit above, consolidating logic that had drifted out of sync across near-identical call sites (the same root cause behind several of the bugs fixed above):
R/00_EMclus.R, R/04C_ParameterEstimation.R,
R/07_Biclustering.R, R/08C_BNM_GA.R, R/09B_LDLRA_GA.R,
R/15_Biclustering_nominal.R, and R/16_Biclustering_ordinal.R used a
carriage-return (\r) "overwrite in place" style, which reads as one
long run-on line when logged to a file or captured non-interactively
rather than the intended per-iteration progress. Switched to \n.beta_posterior_mode() (R/08A_BNM.R) for the
(count + beta1 - 1) / (total + beta1 + beta2 - 2) formula duplicated
across BNM(), LD_param_est() (used by LDLRA()), and BINET() — the
same formula whose numerator/denominator role LDB() had swapped, above.build_conf_mat() (R/07_Biclustering.R) for the conf
(vector/matrix/data.frame) parsing and validation duplicated across
Biclustering(), Biclustering.nominal(), Biclustering.ordinal(), and
LDB() — the same parsing whose missing as.matrix(conf) step caused
the LDB() crash fixed above.stop_if_all_grid_failed(), report_failed_grid_settings(),
and select_optimal_grid_index() (R/00_GridSearch.R) for the
all-cells-failed check, failed-settings warning, and optimal-index
selection duplicated between GridSearch()'s Biclustering (2-D) and
LCA/LRA (1-D) branches — the same duplication in which both the tie-break
bug and the max_ncls < 2 bug above were independently present.tests/testthat/test-regression-1150.R) and smoke tests for
AlphaCoefficient(), OmegaCoefficient(), and BiserialCorrelation(),
which previously had no direct tests
(tests/testthat/test-smoke-reliability.R).R/00_BiclucterUtils.R to R/00_BiclusterUtils.R (filename
typo) and fixed a "Bicluter" comment typo; removed a stray .png file
accidentally committed to the repository root in 2024, and broadened the
.Rbuildignore rule for Rplots.pdf so a leftover
tests/testthat/Rplots.pdf can no longer slip into the CRAN tarball.plot.exametrika() now forwards graphical parameters supplied via ...
to every plot type (R Journal review request). Previously the ...
argument documented on the plot.exametrika() help page was captured but
never passed on: the internal dispatch functions (plot_irt_model(),
plot_grm_model(), plot_common_profiles(), the polytomous Biclustering
plotters, the array plot, and the network plotters) did not accept or
propagate it, so standard graphical parameters were silently ignored. For
example, plot(result.LRA, type = "IRP", items = 1:4, nc = 2, nr = 2, las = 2, pch = 16) changed neither the axis-label orientation (las)
nor the plotting symbol (pch).
User-supplied parameters are now forwarded consistently to the underlying
base R plotting calls (plot, barplot, image, lines, curve) for
all plot types, including manually drawn axes so that las,
cex.axis, and similar parameters take effect on every axis. Standard
parameters such as pch, las, cex, col, lty, and lwd work as
expected, and user values take precedence over the package defaults, so
xlab, ylab, main, etc. can be overridden. This is achieved without
adding any dependency: the package remains lightweight and base-graphics
only.
Internally this is implemented via small helpers (merge_plot_dots(),
call_plot(), draw_curve()) that merge the package defaults with the
user's ... (user wins) and dispatch through do.call(); the IRT/GRM
curves are drawn by grid evaluation through draw_curve() instead of
curve() so that the non-standard evaluation of curve()'s first
argument does not block parameter forwarding.
IRT Test Information Function plot title typo. The type = "TIF"
plot for IRT models had main = "Test Informaiton Function"; it now reads
"Test Information Function" (R/00_plot_irt.R).
User-facing message typos and missing-word fixes. Four cosmetic but
user-visible string fixes following a family-wide audit (exametrika,
ggExametrika, shinyExametrika) prompted by the Clusterd -> Clustered
rename. No behaviour change; messages now read correctly.
R/04C_ParameterEstimation.R: IRT slope warning had a missing space
after the period ("... exceeds 10.Please exercise caution ..."
-> "... exceeds 10. Please exercise caution ...").R/02_TestItemFunctions.R: CCRR.nominal() info message was missing
the verb "Using" ("... binary data only. Conditional Selection Rate for your polytomous data instead." -> "... binary data only. Using Conditional Selection Rate for your polytomous data instead."), now
matching the parallel messages in CCRR.ordinal/CCRR.rated/CCRR.binary.R/09_LDLRA.R and R/10_LDB.R: max-parents warning had a comma
where a period belongs ("... is N, Please check." -> "... is N. Please check.").plot.exametrika() Biclustering panel: typo "Clusterd" -> "Clustered".
Fixed a long-standing misspelling in the base-R Biclustering array plot
(R/00_plot_biclustering.R): the panel title "Clusterd Data" is now
"Clustered Data", the inline comment ## Clusterd Plot is corrected to
## Clustered Plot, and the internal variable clusterd_data is renamed
to clustered_data. No user-facing API change. The downstream
ggExametrika::plotArray_gg() had inherited and propagated the typo into
its argument names (Clusterd, Clusterd_lines, Clusterd_lines_color);
that is being corrected in ggExametrika 1.1.1 as a coordinated release.
Binary IRM: crash on real data with missing values. Biclustering_IRM
on a binary exametrika data object containing missing cells
(tmp$Z[s, j] == 0) crashed with "確率ベクトル中にNAがあります"
(NA in probability vector) inside rmultinom(), after log() produced
NaN in the Gibbs sampler. The cause was that tmp$U retains the
-1 missingness marker from dataFormat() and the IRM aggregated it
directly (tmp$U %*% fld01), producing negative counts in CcfPlus
whose log was NaN. The pre-existing line U <- tmp$U * tmp$Z was a
dead variable. R/07_IRM.R now assigns the masked matrix back as
tmp$U <- tmp$U * tmp$Z, so all downstream tmp$U references see
zero in missing cells. Reproduced on HCI (651×20, 1.5% missing) and
SAT12 (600×32 binarised, 0.4% missing).
Graphical Lasso: graceful handling of numerical divergence. When the
inner block coordinate descent (BCD) of Glasso() produced non-finite
values — typically with high-dimensional, near-singular polychoric
correlation matrices (N <= p regime) and very small lambdas — the
previous implementation crashed in R's if (diff < eps) test with
"missing value where TRUE/FALSE needed", discarding all of the
intermediate work along the lambda grid.
glasso_one() now detects non-finite values in the BCD update of
beta, in the inner convergence statistic diff, and in the working
covariance W, and returns converged = FALSE instead of crashing.Glasso() now reacts to converged = FALSE (or a non-finite EBIC)
by emitting a warning() that names the offending lambda and the
best lambda found so far, and breaks the lambda search, returning
the best solution encountered up to that point. If divergence happens
on the very first lambda the call still errors out, since no solution
exists to return.path element of the result preserves NA for any lambda values
skipped by the early break, making the breakpoint visible to the user.Observed concretely on a 98-item / 143-respondent ordinal dataset
(polychoric condition number ~1.4e5): the old Glasso() crashed at
lambda = 0.0096 after ~25 minutes of work; the new implementation
returns the same lambda_opt = 0.1003 solution with a clear warning
and a path whose two trailing entries are NA.
Resubmission of the 1.13.0 release. The 1.13.0 submission was rejected
by the CRAN auto-check service for a single Overall checktime NOTE on
r-devel-windows-x86_64 (11 min vs. the 10 min limit). The package itself
passed cleanly on both Windows and Debian (Status: OK / OK).
R CMD check comfortably under the 10-minute Windows limit. The
same tests continue to run locally and on R-hub / win-devel via the
NOT_CRAN environment variable that testthat sets.
test-grm.R: the J15S3810 regression (a 15-item / 3,810-respondent
fit, ~60s) and the nitems >= 8 underflow regression (~8s) are
skipped on CRAN. Default-coverage GRM tests on small data continue
to run on CRAN.test-irm.R: the J35S515 shared-fixture Gibbs run (~23s) and the
two reproducibility tests that re-run Gibbs (~22s combined) are
skipped on CRAN. The Default seed is 123 structural check still
runs on CRAN.This release also incorporates the development changes from 1.12.0–1.12.2 (never published to CRAN). Their entries are retained below for traceability.
Graphical Lasso (Glasso): New function for sparse precision matrix
estimation from ordinal item response data. The polychoric correlation
matrix is computed internally and the optimal regularization parameter
is selected by Extended Bayesian Information Criterion (EBIC; Foygel
and Drton 2010). Implements block coordinate descent (Friedman, Hastie,
Tibshirani 2008; Algorithm 17.2 of Hastie, Tibshirani, Friedman 2009)
with cyclical coordinate descent for the inner lasso step. Warm-starting
across the lambda grid accelerates the search.
Internal helpers glasso_one() (single-lambda solver) and
compute_EBIC_glasso() (EBIC computation) are available but not
exported.
Returns a list with theta (selected precision matrix), lambda_opt
(selected lambda), ebic_opt, n_edge, and path (data frame of
lambda, ebic, n_edge over the search grid).
print.exametrika Glasso method: Added a Glasso branch to the
shared print.exametrika() dispatcher to summarize the estimated
model (optimal lambda, EBIC value, edge count, precision matrix).
Chatterjee's xi correlation (chatterjee_xi, xi_stable,
chatterjee_matrix): New family of functions implementing
Chatterjee's (2021) rank-based correlation coefficient.
chatterjee_xi() computes the single-shot value with random
tie-breaking. xi_stable() averages B replications (default
B = 1000) to stabilize against tie-induced variability and returns
a list with the mean, standard deviation, standard error, and B.
chatterjee_matrix() produces the p x p asymmetric pairwise xi
matrix from ordinal data with pairwise-complete handling of
missing values; the asymmetry between xi(j, k) and xi(k, j)
enables direction detection in graphical-model construction.
Glasso() and chatterjee_matrix() in
backticks (e.g., `[0, 1]`, `[j, k]`, `Q[i, j]`) so
that roxygen2 no longer rewrites them as \link{} cross-references.
Resolves R CMD check WARNINGs about missing link targets.Glasso() ... documented: Added @param ... to Glasso() so
that R CMD check's "Undocumented arguments" WARNING is cleared.compute_EBIC_glasso() partial-match fix: Replaced
determinant(Theta, log = TRUE) with the full argument name
logarithm = TRUE, removing the partial-argument-match NOTE.dataFormat(response.type = "rated") no longer errors when a CA
category is unobserved: previously, if no respondent chose item
j's correct-answer category (e.g., everyone got it wrong on a hard
item, or the sample was small enough that the CA category did not
appear by chance), dataFormat() aborted with
"CA for item j is not a valid response category". This rejected
legitimate test data where an item is uniformly difficult or where
the sample size is small. The check has been relaxed to a
warning(); the affected items are encoded as all-incorrect (U[, j] == 0), which is the correct downstream behavior. This applies to
both the matrix-input path (dataFormat) and the long-format path
(longdataFormat). Mistyped CAs (e.g., CA[j] = 8 when categories
are 1–7) still surface as warnings, so the diagnostic value is
preserved.Biclustering.nominal() and Biclustering.ordinal() no longer abort
with "missing value where TRUE/FALSE needed" when extreme grid
configurations (e.g., very small ncls combined with very large
nfld) leave fields or classes empty during EM. The two
log-likelihood checks inside the EM loop — the relative convergence
test and the monotonicity guard — now treat a non-finite
test_log_lik as a non-converged exit instead of throwing. Affected
cells are returned with converge = FALSE so that GridSearch()
skips them automatically when comparing criteria. This also resolves
the same crash in Biclustering.rated(), which calls
Biclustering.nominal() internally.Class-side Confirmatory Biclustering: Biclustering() now accepts
a conf_class argument that fixes class memberships during EM. Like
conf (which fixes field memberships), conf_class accepts either a
vector of class labels (one per respondent) or a 0/1 membership matrix
(respondents x classes). When supplied, the class-side E-step is
overridden every iteration. For Ranklustering (method = "R"), the
neighbour-smoothing step is skipped (smoothed_memb <- clsmemb) since
smoothing pre-fixed labels would defeat the purpose of fixing them.
Available for binary, ordinal, nominal, and rated data; can be
combined with conf to fix both fields and classes simultaneously.
Note: rated re-orders classes by correct rate after estimation, so
the output class labels may not match the input labels (the
individual-to-class mapping is preserved up to relabeling).
Number of model parameters (nparam) is intentionally not adjusted
when conf or conf_class is in effect: only PiFR / BCRM cell
probabilities count as parameters in this implementation, and
membership matrices are latent posteriors, not parameters.
Confirmatory ordinal Biclustering: field membership now stays fixed
during EM: Biclustering.ordinal() was overwriting fldmemb with the
E-step estimate every iteration, so the conf/conf_mat argument only
affected the initial value. Aligned the implementation with
Biclustering.binary() by re-applying fldmemb <- conf_mat immediately
after the field-side E-step. With this fix result$FieldEstimated
matches the user-supplied assignment exactly.
Confirmatory nominal Biclustering: conf argument is now actually
honored: Biclustering.nominal() validated length(conf) against
NCOL(U), but U is an exametrika list object so NCOL(U) always
returned 1. The check rejected every well-formed conf vector with
"conf vector size does NOT match with data.", making confirmatory
nominal Biclustering unreachable since v1.10.0. Replaced with
NCOL(U$Q) (and the matrix-form NROW(conf) check, plus the
conf_mat allocation, in the same way).
Confirmatory ordinal Biclustering: conf length check actually
validates input: same NCOL(U) issue as nominal, but in
Biclustering.ordinal(). Fixed by switching to NCOL(U$Q).
Confirmatory binary Biclustering: conf size checks now reference the
formatted response matrix explicitly: Biclustering.binary() happened
to work because U is rebound to tmp$U * tmp$Z before the conf block,
so NCOL(U) returned the item count. Switched to NCOL(tmp$U) for
consistency with the ordinal/nominal fixes and to make the intent
obvious to readers.
Confirmatory Biclustering: matrix-form conf is no longer silently
dropped: when conf was supplied as a membership matrix
(items x fields) instead of a vector, all three implementations
(binary, ordinal, nominal) validated the matrix but never copied
it into conf_mat, so the subsequent nfld <- NCOL(conf_mat) failed
with object 'conf_mat' not found. Added conf_mat <- as.matrix(conf)
in the matrix branch.
Biclustering.ordinal() now honors the maxiter argument: The
inner EM iteration cap was hardcoded to 100 (maxemt <- 100) and
ignored the user-supplied maxiter. This mirrors the bug that was
fixed for Biclustering.nominal() in 1.11.0. Callers that pass a
larger maxiter (e.g. 2000 in Monte Carlo studies) now actually
use that ceiling.
GRM() fit indices no longer return NaN for moderate or larger
item counts: The benchmark model used
const <- exp(-nitems * 100) as a log-domain epsilon inside
sum(cat_counts * log(cat_probs + const)). For about 8 or more items
this expression underflows to exactly 0 in IEEE-754 double precision,
so log(0) = -Inf propagates through 0 * -Inf = NaN and poisons
every downstream index (model_Chi_sq, NFI, CFI, RMSEA, AIC, CAIC,
BIC). The benchmark and null loops now skip zero-count categories
explicitly, removing the need for an additive epsilon.
GRM() degrees of freedom were computed incorrectly: The model
df was set to n_pattern * (ncat - 1) + 1 and the null df to
n_pattern * (ncat - 1), which is neither (# bench params) - (# model params) nor (# bench params) - (# null params). The inflated
df then clamped CFI, TLI, IFI, and RMSEA to zero whenever
chi^2 < df (and the previous df_A > df_B inequality was the wrong
direction to begin with). The convention now matches
Biclustering.ordinal():
bench_nparam_j = n_pattern * (ncat_j - 1)null_nparam_j = ncat_j - 1model_nparam_j = ncat_j (one slope plus ncat_j - 1 thresholds)df_A_j = bench_nparam_j - model_nparam_jdf_B_j = bench_nparam_j - null_nparam_jGRM() now accepts any integer-coded ordinal responses, not just
1..K: Category counts were derived from apply(dat, 2, max), which
assumes responses are already 1-indexed. Data coded from 0 (e.g.
0..3) or with gaps (e.g. 1, 2, 4) undercounted ncat[j] by one or
more and then indexed grm_prob()[resp] / v[resp] out of range,
producing either an outright invalid subscript type 'list' error
(on J15S3810-style data) or silently truncated threshold columns.
Each item's responses are now remapped to contiguous 1..K codes via
match() against the sorted unique values on entry, with ncat[j]
derived from the mapped levels; both the R model-fit blocks and the
C++ log-likelihood receive 1-based input regardless of the user's
coding.
GridSearch() now tolerates per-cell fit errors: Previously, a
single Biclustering() (or LCA() / LRA()) call that raised an
error at a grid corner (for example, empty-cluster edge cases at
large ncls/nfld with small nobs/nitems) would propagate out
of GridSearch() and abort the entire grid. The call to the
underlying analysis function is now wrapped in tryCatch, and
errors are handled the same way as non-convergence: the cell is
marked NA in the index matrix and the (ncls, nfld) pair is
recorded in failed_settings. GridSearch() still raises only
when all grid cells fail, preserving the existing "all-failed"
error.
C++ implementation of the IRM Gibbs sampler core: The collapsed
Gibbs sampler shared by Biclustering_IRM.nominal(),
Biclustering_IRM.ordinal(), and Biclustering_IRM.rated() is now
implemented in C++ via Rcpp (src/irm_gibbs_core.cpp). The R
reference implementation in R/00_IRM_Gibbs_CORE.R is preserved
behind irm_gibbs_core(..., use_cpp = FALSE) for cross-checking.
Numerical reproducibility (revised 2026-05-07): RNG calls are
routed through R-level sample.int() and rmultinom() (via
Rcpp::Function) so the unif_rand() consumption order matches
base R exactly. The C-level entry points (Rcpp::sample,
R::rmultinom) consume RNG differently from base R and were
explicitly avoided.
The C++ and R reference paths are bit-identical for the first
Gibbs iterations on the bundled test datasets and on real
polytomous data (verified through iteration 2 on J143S32 ordinal),
but diverge from iteration 3 onward on real data due to
sub-LSB floating-point ordering differences accumulating through
the lmvbeta() calls in the CRP likelihood. Once the divergence
flips a single rmultinom() outcome, the two chains follow
different sample paths.
Empirically the two paths still target the same posterior:
on a 50-seed comparison (M = 66 ordinal items, N = 143), the
marginal distributions of n_field and n_class are
statistically indistinguishable between paths (Wilcoxon
p = 0.56 / 0.81). Use either path for inference; do not rely
on a single seed reproducing the same (n_field, n_class)
across the two paths.
The parity tests in tests/testthat/test-irm-gibbs-cpp.R
cover only short runs (<= 10 iterations) and therefore did not
catch this. A longer-iteration test on real-data-like
configurations is on the roadmap.
Wall-clock: roughly 4x speedup on the inner Gibbs loop on the
bundled test datasets (J20S600 nominal: 82 to 20 ms/iter;
J35S500 ordinal: 76 to 19 ms/iter). The end-to-end
Biclustering_IRM() call also benefits proportionally on larger
iteration counts (e.g. simulation runs).
No new dependencies: Rcpp is already in LinkingTo:; the
implementation uses only Rcpp::Function and <vector>/<cmath>.
Vectorized EM hot path in Biclustering.ordinal(): The per-iteration
cost of the ordinal EM loop has been reduced by rewriting five hot spots
in base R without introducing new package dependencies.
apply(Ufcq_prior, c(1, 2), function(x) rev(cumsum(rev(x))))
with an in-place reverse cumulative sum loop over the category axis.
The previous form dispatched nfld * ncls R-level calls per EM
iteration; the new form performs maxQ - 1 vectorized array
additions.apply(X, 1, min) and apply(X, 1, max) with
do.call(pmin.int, as.data.frame(X)) /
do.call(pmax.int, as.data.frame(X)). These compute the row-wise
reduction in a small number of C-level calls rather than one
R-level call per row. apply(X, 1, which.max) was similarly
replaced by max.col(X, ties.method = "first").log(BBRM[,,q] - BBRM[,,q+1] + const) once per EM
iteration instead of independently in the class-side and field-side
E-steps for every q.Z * Uq[,,q] once per fit (call it ZU) via
column-major recycling (Uq * as.vector(Z)), replacing about eight
separate elementwise products per EM iteration and reusing the
same array in the post-EM fit-index blocks. The replicate(maxQ, Z)
allocation that produced Zrep has been removed as a byproduct.apply(X, c(1, 2), sum) and apply(X, c(2, 3), sum) on
3-D arrays with rowSums(X, dims = 2) / colSums(X, dims = 1),
and apply(M, 2, sum) on matrices with colSums(M).Output is bit-identical to 1.11.0 on every tested configuration
(max(abs(old - new)) == 0 for BCRM, BBRM, ClassMembership,
FieldMembership, TestFitIndices, and the full EM trajectory).
Single-fit wall-clock on typical sizes (ncls, nfld up to 20 by 15 on
J35S500) improves by roughly 1.2 to 1.3 times; the per-EM-iteration
speedup compounds across GridSearch() grids.
Vectorized 3-D apply reductions in Biclustering.nominal(): The
same rowSums/colSums(X, dims = ...) substitutions were applied to
the nominal M-step and null-model blocks, and Zrep was eliminated
in favor of the ZU precomputation. Output is bit-identical to
1.11.0.
Vectorized Uq one-hot encoding in Biclustering.ordinal() and
Biclustering.nominal(): The construction of the one-hot response
array Uq[i, j, tmp$Q[i,j]] = 1 previously used a nobs * nitems
nested R loop. It now writes all ones in one C-level matrix-index
assignment, restricted to non-missing cells (tmp$Z == 1). Single-fit
wall-clock on the validation matrix improves by 1.2-5.3x over 1.11.0
(up to 5x on small cases where the loop overhead dominated). Output
is bit-identical to 1.11.0 at every downstream location (the missing
entries of the old Uq were never read because every consumer
applies the tmp$Z mask).
Biclustering.ordinal() now returns SmoothedMembership: The
smoothed (Ranklustering-filtered) class membership matrix is now part
of the return value, mirroring Biclustering.binary(). This fills a
long-standing gap (only the binary form previously exposed it) and
makes Ranklustering with conf_class introspectable: callers can
verify that the smoothing step is correctly skipped when class labels
are pre-fixed (SmoothedMembership equals ClassMembership in that
case).print() smoke test for Biclustering_IRM.rated() results
(test-irm-rated.R).nrank = 4 LRA.rated regression test for
DistractorAnalysis() to cover varying rank counts
(test-distractor.R).Imports or Depends. The new C++ Gibbs
sampler uses Rcpp, which was already declared in LinkingTo:.conf_class argument on
Biclustering() is additive and defaults to NULL, preserving the
previous behavior for all existing callers.Biclustering_IRM.rated(): Rated IRM (Infinite Relational Model for multiple-choice data): Performs IRM-based biclustering for rated data (items with correct answers and multiple response categories). Internally calls Biclustering_IRM.nominal() for estimation via Dirichlet-Multinomial collapsed Gibbs sampler, then post-processes the results: classes are sorted by correct response rate (ranklustering), and two layers of fit indices are reported — binary (item correct/incorrect, with full benchmark model including CFI/RMSEA) and nominal (category-level, AIC/BIC/CAIC only). Returns TestFitIndices (binary, default) and TestFitIndices_nominal (nominal) in the output.
DistractorAnalysis(): Distractor analysis for rated models: S3 generic function for analyzing response category distributions by rank/class. Computes observed frequency tables, proportion tables, chi-square tests against chance level, and Cramer's V effect sizes for each item-by-rank cell. Works with LRA.rated and Biclustering.rated / Biclustering_IRM.rated results. For Biclustering models, output is grouped by field. Supports items and ranks filtering in both print() and plot() methods.
J21S300: New rated sample dataset (21 items, 300 students, 4 response categories) for testing rated Biclustering and IRM models. Smaller and faster than J35S5000 (35 items, 5000 students).Fixed maxiter parameter ignored in Biclustering.nominal(): The maxiter parameter was accepted but internally overridden by a hardcoded maxemt <- 100. Now correctly uses maxiter as the iteration limit.
Fixed spurious "No ID column detected" message in Biclustering.rated() and Biclustering_IRM.rated(): The crr() function was called with a raw matrix (tmp$Z * tmp$U) instead of an exametrika object, triggering an unnecessary dataFormat() call. Now creates a temporary binary-typed exametrika object to avoid the message.
Biclustering.rated(): Rated (multiple-choice) Biclustering: Performs biclustering for rated data (items with correct answers and multiple response categories). Internally calls Biclustering.nominal() for estimation, then post-processes the results: classes are sorted by correct response rate (ranklustering), and two layers of fit indices are reported — binary (item correct/incorrect, with full benchmark model) and nominal (category-level, AIC/BIC/CAIC only). The binary layer uses item-level correct response rates per class without field constraints (nparam = nitems * ncls). Supports both Biclustering (method = "B") and Ranklustering (method = "R") modes. Returns TestFitIndices (binary) and TestFitIndices_nominal (nominal) in the output. Added 16 tests for rated Biclustering using J35S5000.ItemFRP to quasiFRP in Biclustering.rated() output: The item-level correct response rate matrix (J x C) is computed without field constraints because the field assignments come from nominal analysis and have different meaning in the binary context. The name quasiFRP (quasi Field Reference Profile) clarifies this distinction from the standard field-constrained FRP.subscript out of bounds in nominal/ordinal Biclustering with large nfld: The initial field assignment ceiling(1:nitems / (nitems / nfld)) could exceed nfld due to floating-point rounding. Added pmin(..., nfld) clamping, consistent with the binary Biclustering fix already in place.dataFormat() call in Biclustering.rated(): The function was unnecessarily re-calling dataFormat() on already-formatted data, causing repeated "No ID column detected" messages during GridSearch(). Now reuses the existing exametrika object directly.x$U (binary correctness matrix) was used instead of x$Q (polytomous responses). Now correctly uses x$Q for polytomous models.invalid graphics state error: par(mfrow=...) and layout() conflicted when plotting FCRP/FCBR. Now skips par(mfrow=...) for plot types that use layout() internally.id parameter validation in dataFormat(): When a vector (e.g., id = tmp$ID) was passed instead of a column number, the error message was cryptic (the condition has length > 1). Now validates that id is a single integer and provides a clear error message.ell_B = 0) and chi-square based indices (NFI, RFI, IFI, TLI, CFI, RMSEA) meaningless. These indices are now reported as NA for nominal data. Information criteria (AIC, BIC, CAIC) are computed directly from the model log-likelihood and remain available. Also fixed a length(benchGroup) → length(unique(benchGroup)) bug in Biclustering.nominal() (the unique() call was missing).nfld. Affects Biclustering() (binary/ordinal/nominal), Biclustering_IRM() (binary/ordinal/nominal).field to fld in Biclustering_IRM.binary(): Internal variable name changed from field to fld for consistency with other Biclustering models. No change to the public API (FieldEstimated output name is unchanged).@format field incorrectly stated "nominal responses" but the dataset is ordinal (response.type=ordinal). Also fixed "A ordinal" to "An ordinal" in @description.eval=FALSE to computationally intensive code chunks in vignettes (GridSearch, IRM, LDLRA, LDLRA_PBIL, LDB, BINET, ordinal/nominal Biclustering, ordinal/rated LRA). The Japanese guide (guide-ja) now uses eval=FALSE for all model estimation chunks to avoid duplicating computation from English vignettes. Full rendered output is available on the pkgdown site.log(n + 1) but should be log(n) + 1 per Bozdogan (1987, Psychometrika, 52(3), p.358, Proposition 2, Eq.44). The original Mathematica implementation had this error (Log[nobs + 1]), and the R port inherited it. Both the R version (R/00_ModelFitModule.R) and the Mathematica version (develop/mtmk15forVer13/mod/Module_ModelFit.nb) have been corrected. The numerical difference is approximately 1 (constant), but the corrected formula now matches the published definition: CAIC(k) = -2 log L + k * (log(n) + 1). This affects all models that compute fit indices: IRT, LCA, LRA (binary/ordinal/rated), Biclustering (binary/ordinal/nominal), IRM, BNM, LDLRA, LDB, BINET, and GRM. All 873 tests pass with the corrected formula.\donttest to \dontrun: GRM's multi-panel plot examples caused "invalid graphics state" errors in non-interactive environments (pkgdown, CI). Changed to \dontrun to prevent build failures.CA parameter: When CA (correct answer vector) was provided but response.type was not explicitly specified, dataFormat() incorrectly classified the data as "ordinal" instead of "rated". This caused $U (binary scoring matrix) to be NULL. The auto-detection now checks for CA before ordinal/nominal detection: binary → rated (if CA provided) → ordinal → nominal.id parameter not working for column 1: Changed id default from 1 to NULL. Previously, id = 1 (both default and explicit) triggered auto-detection heuristics, making it impossible to explicitly specify the first column as the ID column. Now id = NULL triggers auto-detection, and any numeric value (including 1) forces that column to be used as ID.1:N) as ID regardless of whether they also look like valid response values. Previously, 1:10 was misclassified as response data because looks_like_response_data() returned TRUE.U matrix: When the response matrix contained missing values (-1), the binary scoring matrix U incorrectly scored them as 0 (incorrect) instead of -1 (missing). Now U[i,j] = -1 when Q[i,j] = -1.drop=FALSE missing in item exclusion: When items were excluded due to invalid variance (e.g., containing Inf), the column subsetting response.matrix[, mask] dropped matrix dimensions if only one item remained, causing a crash. Added drop = FALSE.dataFormat() now reports via message() when it detects problematic data:
g_list / adj_list Input Path Fixg_csv variable undefined error when using g_list or adj_list input: BINET() crashed with an undefined variable error at the graph construction step when the DAG was specified via g_list or adj_list parameters. The internal variable g_csv (used to build the integrated graph object all_g) was only defined in the adj_file code path. Added logic to reconstruct g_csv from adj_list for the g_list and adj_list input paths, ensuring all_g is correctly built with Field edge attributes regardless of input method.g_list / adj_list length validation: The length check for g_list and adj_list incorrectly compared against ncls (number of classes) instead of nfld (number of fields). In BINET, each element of adj_list represents the DAG structure at a specific field, so the list length should equal nfld. Also fixed g_list type check from g_list[[1]] (always checking the first element) to g_list[[j]] (checking each element).... to Biclustering_IRM.binary() signature: The S3 generic Biclustering_IRM(U, ...) requires all methods to include ... in their formal arguments. The .binary method was missing it, causing an R CMD check WARNING. Added ... to match the generic and the .nominal/.ordinal methods.msg Field Fixmsg field to correctly distinguish Rank vs Class models: Binary IRM and Ordinal IRM are Ranklustering models (ordered latent classes), so msg = "Rank" is correct. Nominal IRM has no Ranklustering concept (unordered latent classes), so msg = "Class" is correct. Previously all three methods had inconsistent values; now Biclustering_IRM.binary and Biclustering_IRM.ordinal return msg = "Rank", while Biclustering_IRM.nominal returns msg = "Class". Also updated internal column names of cls01 and FRP from "Class" to "Rank" for binary/ordinal IRM. The shared Gibbs core (irm_gibbs_core()) output column names were also updated.LRD AliasLRD (Latent Rank Distribution) alias to Biclustering_IRM.binary() return value: The binary IRM was the only Biclustering model that did not include LRD in its return value (it only had LCD). Other Biclustering models (binary Biclustering, nominal IRM, ordinal IRM) all return both LRD and LCD. Added LRD = clsdist as an alias for LCD to ensure consistency across all Biclustering-family models. This improves interoperability with downstream packages (ggExametrika) that look up LRD when msg == "Rank".Biclustering_IRM() seed default back to 123: Ensures reproducibility by default.t(apply()) Dimension Drop Fixt(apply(log_S, 1, irm_log_to_prob)) dimension drop in Nominal/Ordinal IRM: When the Gibbs sampler converged to ncls=1 or nfld=1, apply() on a single-column matrix returned a vector instead of a matrix, causing t() to produce a 1-row matrix instead of an N-row matrix. This led to incorrect class/field membership assignment. Replaced all 6 instances of t(apply(mat, 1, fun)) with explicit for-loops in both Biclustering_IRM.nominal() (3 locations: EM E-step, final class membership, final field membership) and Biclustering_IRM.ordinal() (3 locations: same). Consistent with the project's apply() caution policy.Biclustering_IRM() Gibbs sampler: Changed exp(ptab - min(ptab)) to exp(ptab - max(ptab)) for numerical stability. The previous implementation subtracted the minimum log-probability, which could cause overflow (exp(large positive) = Inf) when the range of log-probabilities was large, resulting in NaN after normalization. Subtracting the maximum ensures the largest value becomes exp(0) = 1 and all others underflow harmlessly to near-zero values.Biclustering.nominal() using stale log-likelihood in model fit indices: When the EM algorithm exited due to log-likelihood decrease, BCRM was correctly reverted to the previous iteration's value, but test_log_lik was not. The model fit section recalculated the correct log-likelihood as testell, but chi_A, model_log_like, log_lik, and LogLik all referenced the stale test_log_lik instead of testell. Also removed a duplicated model fit code block (copy-paste error).J20S400 response type from nominal to binary: The J20S400.rda dataset was incorrectly stored with response.type = "nominal" despite being a binary (0/1) dataset with -99 as missing values. Missing values were not properly handled (Z was all 1s, Q contained raw -99 values). Regenerated from the original CSV source (develop/sampleData/J20S400.csv) with dataFormat(..., na = -99), now correctly typed as binary with 86 missing values properly masked in Z.Biclustering_IRM.nominal() for nominal/polytomous data: Extends the Infinite Relational Model (IRM) from binary to nominal scale data using a Dirichlet-Multinomial collapsed Gibbs sampler. The Chinese Restaurant Process (CRP) automatically determines the optimal number of classes and fields. After the Gibbs sampling phase, small classes are consolidated and refined with an EM algorithm. The Dirichlet prior concentration parameter alpha controls smoothing of category probabilities.Biclustering_IRM.ordinal() for ordinal/polytomous data: Extends IRM to ordinal scale data. Shares the same Dirichlet-Multinomial collapsed Gibbs sampler as nominal IRM (Phase 1), then applies ordinal-specific EM refinement (Phase 2) with cumulative normalization to enforce monotonic category boundaries. The mic parameter (default TRUE) enforces monotone increasing class ordering by expected score sum. Reports BFRP (Bicluster Field Reference Profile), FRPIndex, TRP, and Strongly/Weakly Ordinal Alignment Conditions (SOACflg/WOACflg).Biclustering_IRM() is now an S3 generic: Dispatches to Biclustering_IRM.binary (existing binary IRM), Biclustering_IRM.nominal (new), and Biclustering_IRM.ordinal (new). Raw data is automatically formatted and dispatched based on response.type.irm_gibbs_core() (in R/00_IRM_Gibbs_CORE.R), used by both nominal and ordinal IRM. Helper functions (irm_calc_Ufcq, irm_lmvbeta, irm_log_to_prob, irm_bic_calc, etc.) are also shared.U_fcq, computing only the contribution of the target student/item rather than recalculating the full array at each step.LCA() and LRA() now support a conf parameter for confirmatory analysis: The conf argument accepts an IRP matrix (ncls/nrank x testlength) where non-NA values are held fixed throughout EM estimation and NA values are freely estimated. This enables test equating scenarios where anchor items retain their known IRPs while new items are calibrated against them.conf has column names, items are matched by label rather than by position. This allows conf to contain a subset of items (anchor items only) or items in a different order than the data. Items in the data but not in conf are automatically set to freely estimated. Unknown labels in conf produce an informative error.somclus() internal function: The Self-Organizing Maps estimation code previously inlined in LRA.binary() (~100 lines) has been extracted into a standalone internal function somclus() in 00_EMclus.R. This parallels emclus() (GTM/EM) and returns the same structure. Also fixed a pre-existing typo (h_cout → h_count) and another (oldsBIC → oldBIC).tidyverse and readxl for reading Excel-based Mathematica reference data. Replaced with 24 modern test files using base R read.csv() and test_path() to load CSV fixtures from tests/testthat/fixtures/mathematica_reference/. The new test suite has zero external package dependencies beyond testthat.readxl/tidyverse from test dependencies: DESCRIPTION Suggests no longer requires any packages beyond knitr, rmarkdown, and testthat.test-irm-nominal.R): 18 test blocks for Biclustering_IRM.nominal() using J20S600 data — basic execution, dimensions, FRP validity, membership, fit indices, backward compatibility, seed reproducibility, alpha validation.test-irm-ordinal.R): 25 test blocks for Biclustering_IRM.ordinal() using J35S500 data — basic execution, dimensions, FRP validity, expected scores, TRP, BFRP, FRPIndex, SOAC/WOAC flags, mic parameter, fit indices, backward compatibility, seed reproducibility, alpha validation.stem(nrs(dataFormat(J15S500))) to demonstrate score distribution visualization.^tests$ and ^inst$ entries that were incorrectly excluding tests and vignettes from the built package.index Parameter FixesGridSearch() now accepts common aliases for fit indices. "loglik", "log_lik", "LogLik", and "LL" are mapped to "model_log_like". "Chi_sq" and "chi_sq" are mapped to "model_Chi_sq". Previously, using these aliases caused silent NULL extraction and eventual errors.GridSearch(), before the computationally expensive grid search loop runs. The error message lists all valid options including available aliases.model_log_like is now correctly treated as a maximization target (larger log-likelihood = better fit). Previously, it was incorrectly placed in the minimization group, which would have selected the worst-fitting model.tolerance = 1e-2 in expect_equal) to absolute difference comparison (threshold 0.005). Q3 residual correlations include values near zero where relative tolerance comparisons are unreliable.layout and plot.new from graphics to NAMESPACE. These functions are used by the legend strip layout helpers for polytomous Biclustering plots.apply(U$Q, 2, unique) returning matrix instead of list: When all items have the same number of response categories (e.g., all 5-point Likert), apply() returns a matrix rather than a list. The subsequent lapply() then iterates over individual elements instead of per-column vectors, causing ncat to be a vector of 1s and crashing the algorithm. Replaced with lapply(seq_len(nitems), function(j) sort(unique(U$Q[, j]))) which always returns a proper list. Same fix applied to catfreq999 computation. Both LRA.ordinal and LRA.rated were affected.apply(tmp$Q, 2, table) returning matrix instead of list: Same class of bug as the LRA.ordinal/LRA.rated fix above. In GRM(), the null model log-likelihood computation used apply(tmp$Q, 2, table) to compute per-item category frequencies. When all items have the same number of response categories (e.g., all 5-point Likert), apply() returns a matrix instead of a list, causing response_list[[j]] to extract a single number rather than the full frequency table. This produced incorrect null_log_like values and consequently wrong chi-square statistics, RMSEA, TLI, and CFI in ItemFitIndices. Replaced with lapply(seq_len(nitems), function(j) table(tmp$Q[, j])) which always returns a proper list.LRA.ordinal() now raises an informative error when items have different numbers of response categories (e.g., some items with 3 categories and others with 5). The internal matrix algebra uses fixed-stride indexing that assumes uniform category counts. The error message suggests alternatives (LRA.rated, Biclustering.ordinal) that support mixed category counts via list-based designs.plot(lca_result, type = "FRP") passed validation but failed at runtime because the $FRP field does not exist in LCA/LRA return values. Now properly rejected with an informative error message at the validation stage.stat parameter supporting "mean" (default), "median", and "mode".style parameter supporting "line" (default) and "bar" (stacked bar chart).stat parameter.FRPIndex (Field Reference Profile indices) including location parameters (Alpha, Beta), slope parameters (A, B), and monotonicity indices (Gamma, C).plot.exametrika()stat: Controls the summary statistic for FRP and RRV plots ("mean", "median", "mode").style: Controls the display style for FCRP plots ("line", "bar").#808080) to distinguish from white (incorrect) and black (correct). Polytomous data uses black (#000000) to distinguish from the category color palette.print() now displays FRPIndex (Field Reference Profile Indices) with a note that the values are based on normalized expected scores (E[score]-1)/(maxQ-1).?Biclustering, including detailed definitions and polytomous adaptation logic.Systematic unification of return value structures across all analysis functions for consistency and interoperability.
n_class (retaining Nclass for backward compatibility)n_rank, n_field (retaining Nrank, Nfield)n_class, n_field (retaining Nclass, Nfield)n_class, n_field, n_cycle (retaining Nclass, Nfield, N_Cycle)n_class, n_field, n_cycle (retaining Nclass, Nfield, N_Cycle)n_cycle, N_Cycle (retaining em_cycle, EM_Cycle as IRM-specific aliases)log_lik Added to All Functionslog_lik at the top level of the return object, matching TestFitIndices$model_log_like.ModelFit class, matching the format used by IRT/LCA/LRA and other functions. Previously these used bare calcFitIndices() output (9 fields, no class). Added fields: model_log_like, bench_log_like, null_log_like, model_Chi_sq, null_Chi_sq, model_df, null_df.ModelFit class.TestFitIndices added as primary name for multigroup fit indices (previously only MG_FitIndices). MG_FitIndices retained as backward-compatible alias. SM_FitIndices (saturated model) remains unchanged.Estimate column (most probable class assignment) to the Students matrix.Estimate column (most probable class assignment) to the Students matrix.FRPIndex (Field Reference Profile indices) for consistency with Biclustering.binary and LDB.TRP from matrix(1×ncls) to numeric vector, consistent with all other functions.class = c("exametrika", "GridSearch") to return value for method dispatch support.msg field assignment: Fixed msg <- "Class" to msg = "Class" inside structure() call (line 150 of 05_LCA.R). The <- operator was being interpreted as a standalone assignment rather than a named list element, causing the msg field name to be empty.plot(r, type="RMP", students=1)). Added drop = FALSE to prevent matrix-to-vector coercion when extracting a single row from the Students matrix.TestFitIndices log-likelihood fields: Fixed null_log_like which incorrectly stored the saturated model log-likelihood (log_lik_satu) instead of the null model log-likelihood (sum(null_itemll)). Added missing bench_log_like field containing the saturated model log-likelihood (sum(satu_itemll2)). Note: the chi-square-based fit indices (NFI, CFI, TLI, RMSEA, AIC, BIC, etc.) were always computed correctly; only the stored log-likelihood labels were affected.TestFitIndices and ItemFitIndices to the standard 16-field structure with ModelFit class (c("exametrika", "ModelFit")), matching all other analysis functions. Added model_log_like, bench_log_like, null_log_like to ItemFitIndices. Removed ScoreRankCorr / RankQuantCorr from TestFitIndices (already available at the top level of the return object).test-12OLR.R and test-13NLR.R to handle the new ModelFit class (add unclass() before as.data.frame()) and adjusted column/index references to match the unified 16-field structure.layout() with a thin dedicated row (height ratio 0.2) to reduce visual clutter in data panels. Added setup_legend_layout() and draw_legend_strip() internal helper functions.P(Q>=1)=1.0 reference line at the top of FCBR plots for visual completeness.dataFormat() now correctly identifies factor-type ID columns before converting factors to numeric. Previously, the factor-to-numeric conversion occurred before ID detection, causing factor ID columns with many levels (>=20) to trigger a "Too many categories" error instead of being recognized as IDs.is_response_data()) that was defined but never called.obj$Q instead of obj$U for test length calculationpmax(Ufcq_prior, 1e-10) to prevent division by zero and NaN errorsThis release improves naming consistency across the package while maintaining full backward compatibility through a deprecation path.
All analysis functions now return results with snake_case field names for better consistency:
n_class - Number of latent classes (replaces Nclass)n_field - Number of latent fields (replaces Nfield)n_rank - Number of latent ranks (replaces Nrank)n_cycle - Number of EM iterations (replaces N_Cycle)log_lik - Log-likelihood value (replaces LogLik)The old field names continue to work for backward compatibility:
Nclass, Nfield, Nrank, N_Cycle, LogLik - Deprecated but functional# Old code (deprecated but still works)
result <- LCA(data, ncls = 3)
n_classes <- result$Nclass
# New code (recommended)
result <- LCA(data, ncls = 3)
n_classes <- result$n_class
Progress messages now display properly in R Markdown documents:
verbose parameter (default: TRUE)
ncls = 2: nfld=2 nfld=3 nfld=4 ...verbose = FALSE to suppress all progress messagesiter 1: match=0 nfld=15 ncls=30Adjusting classes: BIC=-99592.5 ncls=21 (min size < 20)verbose = TRUE to verbose = FALSE (consistent with binary/rated versions)Saturation Model - iter 1: log_lik=-1.234567Restricted Model - iter 1: log_lik=-0.987654verbose = TRUE to verbose = FALSE (consistent with binary/ordinal versions)All progress messages now use proper line breaks instead of carriage returns, ensuring clean output in R Markdown/knitr documents and web documentation.
sort(unique(...)) to ensure consistent ordering: 0 (white/incorrect) → 1 (black/correct)testEll → test_log_lik)log_lik terminologyconverge variable to all EM-based functions to indicate algorithm convergence status
New function LRA() now supports ordinal response data
Bug fixes and improvements
Standardized terminology: unified the usage of "class" and "rank" throughout the package
Various minor bug fixes