R/imputeCellGLoc.R
imputeCellGLoc.RdEstimates the mean structure \(B\), scatter \(\Sigma\) and cell weights \(W\) of the model \(x_i | u_i \sim N(B' u_i, \Sigma)\), where \(u_i\) is a design row built from the categorical variables. Detection uses the conditional residual of a cell given the other continuous cells in its row and the row's categorical pattern.
a data.frame with continuous and categorical columns.
one-sided formula for the categorical mean structure.
~ . (default) is main effects over all categorical columns,
~ .^2 adds interactions, ~ 1 is intercept only.
"soft" for redescending weights in [0, 1],
"binary" for the penalised cellwise MCD objective.
maximum number of outer iterations. 0 returns the
starting fit; under the categorical EM that is the first start's, and the
second start (see categorical) does not run. The default is 200
rather than 50 because the relaxed iteration needed up to 51 steps across
the sweep in .gloc_damp; converged fits leave the loop early, so the
cap costs nothing. Failing to
converge within maxit warns, and the warning says whether the
limit was merely exhausted or the iteration is cycling. Note that
maxit changes the trajectory and not only the stopping
point, because the fallback in .gloc_damp fires on a stall counter
and what is left of the budget then decides whether the restarted run
finishes: truncating a non-converged fit early is not the same answer
sooner, and measured deviations are non-monotone in maxit.
convergence tolerance, applied to the scaled change in the
fitted means, to the relative change in the scatter, and to the cell
weights' fixed-point residual \(\max|f(W) - W|\). The weight step is
divided by damp before the comparison, because the step itself is
damp times that residual: without the division a relaxed run would
stop at a proportionally looser residual than an unrelaxed one, and
eps would not mean the same thing at two relaxation factors. The
mean and scatter terms are not divided, since neither is exactly
damp times a fixed-point residual; the weight term is in practice
the binding one.
minimum fraction of unflagged cells per column (binary corner).
The robust start of the soft corner runs cellMCD at its own tolerance,
.gloc_start_alpha = 0.5, and ignores this argument.
tuning constant of the Tukey bisquare (soft corner).
Inf disables downweighting.
a cell is conditioned on only when its weight exceeds
this, so that a downweighted peer is treated as absent rather than as
evidence. The threshold is applied over a narrow band rather than at a
point – a peer whose weight is within .gloc_peer_band of it is
conditioned on with its information discounted, which is what keeps the
weight map continuous; see .gloc_cond_resid. Cells outside the band
are unaffected, so the meaning of this argument is unchanged.
The default 0.5 is the conventional 1% flagging rule. Raising
it discards more peers. A threshold of 0 does not disable peer
filtering: the damped weight update multiplies a weight by
\(1 - d\) each time the bisquare sends it to zero, so a contaminated
weight decays geometrically towards zero without ever attaining it
(measured around 4e-6 at convergence, with no cell exactly 0), and a
zero threshold readmits those cells at full influence. Use a negative
value to condition on every finite peer, which is useful only for
demonstrating what the filtering buys. The same rule selects the cells a
missing cell is imputed from (since 7.4.1).
half-width of the band around peer_w_min over which
a peer fades out of the conditioning set instead of leaving it at a step;
see .gloc_peer_band for the value and .gloc_cond_resid for
the construction. 0 restores the hard cut, which is useful only for
demonstrating what the band buys: it makes the weight map discontinuous
again, and with it the limit cycles that made converged = FALSE the
usual outcome of a design = ~ . fit. It does not reproduce a
release from before 7.4.0, because the relaxation floor moved at the same
time and the two interact: measured on 260 pooled fits, the hard cut
converges 208 times at the current floor of 0.5 against 245 at the old
0.25. A low floor helps under the hard cut and merely costs iterations
under the band.
Setting it changes the estimator and not only the iteration. Measured over
520 paired fits against peer_band = 0 on the same data, the scatter
moves by up to about 20% either way (the largest relative Frobenius
changes seen were -19.9% and +20.7%, and 22% of fits move by more than
1%). The mean effect is near zero – 0.8401 against 0.8389 in relative
error against a known truth – the direction is not predictable from
anything but the design arm itself, and detection is unaffected (F1 0.4837
against 0.4842). The band is a fix for convergence, not for accuracy, and
is not claimed to improve the estimate.
relaxation factor for the soft corner's weight update.
NULL (default) runs the adaptive schedule of
.gloc_damp_start / .gloc_damp_shrink / .gloc_damp:
it starts unrelaxed, strengthens only when the weight change stops
falling, and falls back once to the cold start at the floor if it is
still cycling there. A number in (0, 1] pins the factor instead, with no
backoff and no fallback; damp = 0.25 is the relaxation factor used
by releases before 7.4.0 and is what the tests compare the schedule
against, the floor having since risen to 0.5 because the band removed the
cycling the low floor was there to suppress.
It is not bit-identical to those releases, because they also
compared the undivided weight step against eps and so stopped
four times earlier than eps asked; see eps. Ignored for
weights = "binary", which takes \(W\) from
cellWise::cellMCD and never relaxes it.
convergence tolerance of the EM inside
cellWise::cwLocScat, the scatter step. That step is about 95% of an
outer iteration, and cwLocScat's own default of 1e-12 is seven
orders of magnitude tighter than eps, so it refines digits this
function immediately discards. See .gloc_scatter_soft.
print progress.
starting values for the soft corner. "robust" (default
since 7.4.1) fits each continuous column on the categorical design alone
along robustbase::lmrob's M-S path (an L1 fit followed by an M-step
at the L1 residual scale) and takes the starting flags from
cellWise::cellMCD on those residuals at .gloc_start_alpha;
see .gloc_start_robust. "classical" starts with every
observed cell at weight 1 and the mean structure by least squares. A
redescending weight function started from a non-robust fit can settle on
a masked solution, which is why the default changed. Ignored for
weights = "binary", whose first step already calls
cellWise::cellMCD. Under categorical = "em", "robust"
can run the fit from two starts and keep the better fixed point; the
paragraph "Two starts" of categorical says when.
start selects the fixed point, not only the path to it: the
iteration can have several fixed points, and the starting values choose
among them. On clean data (n = 200, six continuous and six categorical
variables, 20% missing) at eps = 1e-8, the two starts reached
different fixed points in 7 of 10 fits with design = ~ . (relative
scatter difference 0.009 to 0.032, 2 to 14 cells flagged differently,
unchanged as the tolerance tightens) and in 4 of 9 converged fits with
design = ~ 1. Which fixed point is reached depends on the starting
mean and the starting weights together, and on how the weights are built:
in the 7 differing ~ . fits, the start's hard flags with a
least-squares starting mean reached the robust start's fixed point in 4;
soft starting weights on the robust residuals reached it in 6 when computed
from conditional residuals and in none when computed from marginal ones.
The conditional construction used the start's hard flags twice: as the
cell weights of the scatter it conditioned with, and as the set of peers
it conditioned on.
(Corrected: an earlier version said the hard starting flags were not the
cause, on the evidence of one construction of soft starting weights.) Under
contamination the classical start can mask: with 20% of cells shifted by 10
its scatter error was 6.45 against 0.14 for the robust start.
With "classical" the estimates B, Sigma and W
reproduce VIM 7.4.0; the imputations do not, because since 7.4.1 a missing
cell is imputed from unflagged cells only. (Corrected: this page said
"classical" "reproduces VIM 7.4.0 exactly", which stopped being true
for $imputed under both starts.)
how missing categorical cells are treated. "em"
(default since 7.5.0) imputes them under the model itself. Each
categorical variable gets a multinomial logistic regression on the other
categorical variables. A row with a missing cell enters the estimation
once per candidate level, weighted by that level's posterior probability,
and the posterior combines the prior with the density of the row's
unflagged continuous cells (the peer rule of detection and imputation). In
the soft corner each level is judged with its own cell weights, and flagging a
cell under a level costs a penalty (since 7.5.1; see "Per-level detection" in
Details). The
fit is an EM-type algorithm for a pseudo-likelihood, because conditional
models define no joint distribution, and no monotonicity is claimed. A row
missing several categorical cells is handled by mean-field sweeps, an
approximation; cat_multi_missing is the share of rows with two or more
missing categorical cells, and a row with more level combinations than the cap
takes the most frequent level instead of being swept, with a warning. In the
binary corner Sigma comes from cellWise::cellMCD on one residual
row per observation, so the level uncertainty enters B but not
Sigma; with maxit = 0 no iteration runs, and Sigma then
comes from the weighted pseudo-rows in both corners. Observed categorical cells
are trusted; see cat_prob_observed. "level" keeps the 7.4.x
behaviour: a missing categorical value becomes an extra level of the design and
stays NA in imputed. Without a missing categorical cell the two
give the same B, Sigma, W, imputed and
criterion, bit for bit.
Two starts. Under "em", when a categorical cell is missing,
weights = "soft", start = "robust", the design has a
categorical term and maxit >= 1, the fit runs through the same
iteration from two starts and returns one of the two fixed points. The first
start is the robust start on the expected design rows, fitted on the rows
whose categorical cells are all observed. The second is the robust start of
"level", fitted on every row with the missing value as a level, with
that level's coefficients dropped; where the extra level changes the coding
of the others, as with an ordered factor's polynomial contrasts, its fitted
means are carried over by least squares instead. Both begin with the
posteriors at the priors. The first start sees fewer rows, and in a heavily
contaminated column it can miss part of the contamination, which the
iteration then keeps masked. In the investigation behind this rule (block
contamination at eps = 0.20 with a shift of 6, n = 200, six continuous
and six categorical variables), the fit ended at the good fixed point in 37
of 40 datasets with both starts against 30 of 40 with the first alone.
The second start's fixed point is returned only if that run converged and
the first did not, or both converged and the second has the smaller
binary-corner objective: minus twice the Gaussian log-likelihood of the
retained residual cells \(x_{ij} - (B' \bar u_i)_j\) (weight at least 1/2),
plus \(\lambda_j = \chi^2_{1;0.99} + \log 2\pi + \log c_j\) for each flagged
cell, where \(c_j = 1 / (S^{-1})_{jj}\) and \(S\) is the scatter of the
first start's cellWise::cellMCD call, the same for both runs. Ties,
and two runs that did not converge, keep the first start. The objective tells
a masked fixed point from an unmasked one, not two nearly equivalent ones.
The choice is discontinuous in the data, as for a multi-start MCD: a small
change in the data can switch the fixed point returned. Such a fit takes
about twice as long. If the first start's cellMCD call fails, or its
scatter gives no valid penalty (it is singular, or its inverse has a
non-positive diagonal entry), only the first start runs. An error in the
second start does not stop the fit: it warns and returns the first start's
fit. The warnings, iterations, converged and criterion
are those of the run returned, and the warnings of a returned second start
read "cellGLoc: (second start) ..."; em_starts reports both
runs.
"em" needs cellWise once a categorical cell is missing: the
scatter is then taken over the pseudo-rows, whose cell weights carry each
level's posterior, and only cellWise::cwLocScat turns those into a
case weight linearly. The weighted pairwise fallback of
.gloc_scatter_soft would square them, so it is refused rather than
taken here; install cellWise or use "level".
With "em" the cat_prob_observed diagnostic fits the multinomial
priors on every call, including a call on data with no missing categorical
cell, where the estimation itself needs no priors at all.
a list with B, Sigma, W, U,
imputed, converged, iterations, criterion,
cat_posterior, cat_prob_observed, cat_multi_missing,
cat_priors, em_starts and cat_weights.
W is the \(n x p\) matrix of cell weights of the continuous
variables, 0 on missing cells. Under per-level detection ("em", soft
corner, a missing categorical cell; see Details) a row with a missing
categorical cell reports the posterior mixture
\(W_i = \sum_k r_k W_k\) of its candidates' weights under the final E-step;
a row without one, and a row above the combination cap, reports its own.
imputed is data with its missing cells filled: a continuous
cell by its conditional expectation given the unflagged cells of its row, a
categorical cell (under "em") by its posterior mode. Under per-level
detection a missing continuous cell of a row that also misses a
categorical cell is imputed instead by the posterior mixture
\(\sum_k r_k E[x_{mis} | \textrm{cells clean under } k, u_k]\) over that
row's candidates: each candidate contributes the conditional expectation at
its own design row, conditioned on the cells that are clean under it, and the
candidates are mixed by their posterior probabilities \(r_k\). That is the
exact posterior mean of the cell. Every other row, a row above the
combination cap included, is imputed at its expected design row with the
returned W, as in 7.5.0.
cat_weights is NULL unless per-level detection ran. Then it is a
list with one entry per candidate of the rows below the combination cap that
miss a categorical cell, in the order of the candidate table: row, the
data row (integer); levels, a data frame of the candidate's values of
every categorical variable, observed values included – its columns are the
candidate table's factors, on the level sets the fit observed, so a
character or logical column of data comes back as a
factor; prob, the
candidate's posterior probability, the product of the marginal posteriors when
several cells are missing; and W, the candidates' cell weights, a
matrix with one row per candidate and one column per continuous variable.
W[i, ] of the fit equals colSums(prob * W) over row i's
entries.
Under categorical = "em", U holds each row's expected design
row, so U %*% B are the fitted means. cat_posterior is a
list with one matrix per categorical variable that has a missing cell: its
rows are those cells' rows (row names are row indices), its columns the
levels, its entries posterior probabilities. imputed holds the
posterior mode (ties: the first level). For multiple imputation, draw the
level from these probabilities first and the continuous cells given it.
cat_prob_observed is an \(n x k\) matrix, NA where the cell is
missing: for every observed categorical cell, the posterior probability of its
own level computed as if the cell were missing, without a refit. A small value
points at a miscoded cell; the estimation does not use it. It is NA in
every column of a row above the combination cap as well, since such a row has
no pseudo-rows to compute it from. Under per-level detection it is not on the
same scale as cat_posterior: the diagnostic keeps 7.5.0's un-penalised
Gaussian density, evaluated at the returned (mixture) W, while the fit's
own posteriors now come from the cellwise-penalised per-candidate score
(.gloc_cat_score). Per-level scoring of observed cells is deliberately
deferred, so compare a value of this matrix with other values of it and with
the same fit's flagging rate, not with an entry of cat_posterior.
cat_multi_missing is the share of rows with two or more missing
categorical cells. cat_priors holds the final multinomial prior models,
so that level posteriors for other rows can be computed under this fit (for
example in bootstrap multiple imputation). All four are NULL under
categorical = "level", and under "em" when the data have no
categorical column.
em_starts is NULL unless the fit ran from two starts (see
categorical). Then it is a list: chosen, "complete" for
the first start or "na_level" for the second; objective,
converged and iterations, each a vector named complete
and na_level; and lambda, the penalty per continuous variable
in the objective. When the second start failed, its three entries are
NA and error holds its error message. iterations,
converged and criterion of the fit are those of the chosen run.
criterion is the named vector (means, scatter, weights,
scatter_spread, categorical): the three stopping residuals as of the last
iteration, compared against eps, plus the elementwise spread of
\(\Sigma\) over the last .gloc_stall_iters iterations relative to its
largest variance. scatter_spread is reported only when converged
is FALSE, and is NA otherwise: on a converged fit the window
still holds the last steps of the approach, so a value there would say nothing
about stability. When the fit did not converge it says how much the returned
scatter depends on where maxit happened to stop – a settled cycle
reports its amplitude, a run still drifting reports the drift – and a caller
can test it instead of parsing a warning string. The fifth entry,
categorical, is the largest change of a categorical posterior
probability in the last iteration; it is 0 without the EM, NA when no
iteration ran, and part of the stopping rule with the EM. Under per-level
detection weights is taken over the candidates: the largest change of a
candidate's cell weight times the candidate's posterior probability (1 for a
row that is a single candidate or lies above the combination cap), divided by
the relaxation factor. The adaptive relaxation and the stall detector use the
same number, so a candidate the posterior has all but ruled out cannot hold up
convergence.
It does not tell those two apart, and it is not an error estimate. On 10
non-converged weights = "binary" fits it was positive for all 10,
while raising maxit left 8 of them bit-identical and moved 2 (by
0.008 and 0.096 relative, the latter matching its reported spread
exactly). That is the intended behaviour – all 10 return a scatter that
depends on the stopping point – but the rank correlation with what
raising maxit actually does is only 0.16, so read it as "this
answer is not settled", never as "it would move by this much".
Every entry is NA at maxit = 0.
With design = ~ 1 the estimator reduces to the published continuous-only
estimators: to the cellwise MCD of Raymaekers and Rousseeuw (2024) with
weights = "binary", and to the cellwise weighted maximum likelihood
estimator of Rousseeuw (2026) with weights = "soft".
The soft corner corrects the weighted scatter for its downweighting.
Without a correction the bisquare downweighting deflates the scatter by
about 21% at the default tuning, which would inflate the standardised
residuals and make the estimator over-flag. The correction is per column
and depends on \(\Sigma\), not a single constant: the weights act on
conditional residuals, so they shrink only the unpredictable part of
each cell, and they leave its covariances with the other cells unbiased to
first order. Since 7.5.1 the correction therefore adds the missing part of
each conditional variance to the diagonal and keeps the off-diagonals; see
.gloc_correct_scatter. It vanishes at psi_c = Inf. It is a
first-order correction: on clean data with five variables at correlation
0.7 (n = 10000) it returns conditional variances 4.5% too large and a block
correlation of 0.682, where the symmetric rescaling used until 7.5.0 returned
conditional variances 13% too small and a block correlation of 0.742.
The estimator is the triple \((B, \Sigma, W)\) and iteration stops only
when all three have settled:
\(\max|U(B - B_{old})| / \sqrt{\max \mathrm{diag}(\Sigma)} < eps\),
\(\max|\Sigma - \Sigma_{old}| / \max \mathrm{diag}(\Sigma) < eps\) and
\(\max|W - W_{old}| / damp < eps\). The first term is a change in the
fitted values measured in units of the scatter, not a change in a
coefficient measured against a location, so the criterion is invariant to
shifting the data; the second is relative, so it is invariant to rescaling
it; the third is the weight map's fixed-point residual, relaxation divided
out. converged = TRUE therefore means that nothing this function
returns is still moving, which is the only reading of convergence a
detection method can defend: \(W\) is not a nuisance quantity here, it is
the flag set the caller uses.
The scatter condition is a formal guard rather than an active one. Across every fit measured so far it has never been the binding constraint: at the iteration where the other two conditions are first met it stands at most a seventh of its own tolerance. It is there because \(\Sigma\) is returned and was the one returned quantity the rule did not test, and because it costs nothing; it is not there because it was observed to catch anything.
The soft corner's weight update is relaxed (see .gloc_damp) to keep
the fixed-point iteration inside its contraction radius, and the factor is
adapted rather than fixed; see damp. Relaxation alone does not make
the criterion attainable: conditioning on a peer only when its weight
exceeds a threshold makes the weight map discontinuous, and a discontinuous
self-map of \([0,1]^{n \times p}\) need not have a fixed point at all, so
before 7.4.0 the iteration could be asked to reach a state that did not
exist – and in the design = ~ . arm it usually was. The peer band
(see .gloc_peer_band) removes that, and the two work together rather
than one of them doing the work: over 260 fits of the
categorical-mean-structure arm, converged is 131 under the 7.3.1
schedule and hard cut, 181 with the relaxation floor at 0.5 and the scatter
condition but the cut still hard, and 245 with the band as well. About 40%
of the gain is the floor, the rest the band.
Convergence is not universal and should not be assumed. On that 260-fit grid it is 245, or 94%, with the failures concentrated at correlations of 0.6 and above combined with 20% of cells contaminated.
All convergence figures in this section were measured before 7.4.1, with
the classical cold start (start = "classical"); the robust start has
not been measured on that grid.
When the design does not use a factor's own contrasts attribute – the
factor has an unused level, or, under categorical = "level", its missing
values become a level of their own – the fit says so once for that factor
(since 7.5.1; see .gloc_contrast_lost), and base R's own "contrasts
dropped" warning is muffled at every design site here; the fitted means,
\(\Sigma\), \(W\) and the imputations are the same under any full-rank
coding, so only B's rows change.
Continuous columns that are integer in data stay
integer in $imputed; their conditional expectations are
rounded.
A missing continuous cell is imputed by its conditional expectation given
the unflagged cells in its row, selected by the same peer rule detection
uses (peer_w_min with the peer band), in both weight corners. Until
7.4.0 the imputation conditioned on every observed cell, flagged ones
included, so a grossly contaminated cell was carried into the imputation of
its row-mates; B, Sigma and W were not affected. See
.gloc_impute.
Per-level detection (since 7.5.1). Under categorical = "em"
in the soft corner, a row with a missing categorical cell is judged once per
candidate level, or per level combination when several cells are missing. Each
candidate carries its own cell weights: its cells are compared with its own
fitted mean and conditioned on the cells that are clean under it. The E-step
scores a level by the cellwise-penalised likelihood of the binary corner, the
Gaussian log-density of the cells retained under the level minus
\(\lambda_j / 2\) for each cell flagged under it, with
\(\lambda_j = \chi^2_{1;0.99} + \log 2\pi + \log c_j\),
\(c_j = 1 / (\Sigma^{-1})_{jj}\) from the current scatter, and the peer band
in between (see .gloc_cat_score). B and Sigma are updated
with each candidate's weights times its posterior probability, and W
reports the posterior mixture of the candidates' weights, which
cat_weights returns. Why: until 7.5.0 detection in such a row ran at its
expected design row, where a cell that decides the level looks outlying
whenever the true level is far from the prior mean. That cell was flagged, the
E-step dropped it, and the posterior could not move, a self-locking fixed
point. On the datasets of the categorical arms of the 7.5.0 simulation at
eps = 0 (n = 200, six continuous and six categorical variables, 10
replicates per pattern), the decisive cell was flagged in 17.6-19.8% of the
rows missing the factor, and the imputed level was right in 0.809-0.861 of those
rows, against 0.965-0.986 for a classifier that uses the fit's own \(B\) and
\(\Sigma\), the true prior probabilities and every observed continuous cell. A
contaminated cell is flagged under every level and pays the same penalty under
each, so it still does not steer the level. The binary corner keeps one weight
row per observation, since cellWise::cellMCD returns one: its detection
stays at the expected design row, and the lock can occur there. The soft corner
is the default. The missing continuous cells of such a row are imputed by the
posterior mixture of its candidates' conditional expectations, each taken at
that candidate's design row and conditioned on the cells that are clean under
it (.gloc_impute_mix); see imputed under Value. (Corrected: until
7.5.0 they were imputed once, at the expected design row with the row's weights,
which was the exact mixture expectation only while every level shared one
weight row.)
Raymaekers, J. and Rousseeuw, P. J. (2024). The cellwise minimum covariance determinant estimator. JASA 119(548), 2610-2621. doi:10.1080/01621459.2023.2267777
Rousseeuw, P. J. (2026). Analyzing cellwise weighted data. Econometrics and Statistics 38, 31-41. doi:10.1016/j.ecosta.2023.01.007