The streaming layer, which the previous release described as declared and deliberately not built. Everything here is additive: nothing that ran under 1.0.0 runs differently.
proximity_stream() builds an object that holds what a
proximity matrix is made of rather than the matrix.
mantel_test(), cka() and
rv_coefficient() take it and manufacture the matrix a block
of rows at a time, accumulating the statistic and discarding each block,
so the object is never allocated. The streamed answer is the dense
answer: over the cells of inst/simulations/streaming-cost.R
the two agreed to within 1.8e-13 on the statistic, on every permuted
statistic, and on the count of usable pairs, in-bag and out-of-bag
alike.mantel_test(), cka() and
rv_coefficient() gain a block_size argument,
used only when the arguments are streams. The default divides a 64 MB
budget by the sample size.README and vignette("large-n") said
mantel_test() would gain a streaming argument.
It cannot have one. mantel_test() takes a
proximity, which is a matrix that has already been built,
and a matrix that has already been built cannot be traversed instead of
built. The saving exists only at the moment the object is constructed,
so it belongs to the type of the input rather than to a flag on the
function, and the feature arrived as a constructor.print() shows both figures side by side rather
than leaving it to be assumed.method = "spearman" errors on the streaming path. A
rank is a statement about every other pair, so it cannot be accumulated
from blocks that have been discarded.make_psd(), because the repair is
an eigendecomposition of the object the stream exists in order not to
build.ensemble_inbag() collects the two engines’ spellings of
the bootstrap counts, which three call sites now need. It replaces two
copies of the same branch in as_proximity.randomForest()
and as_proximity.ranger().First release. The four phases of the roadmap are complete and the
interface is settled: as_proximity() was the last name to
move, and it moved so that this one would not have to. Breaking changes
from here go through a deprecation cycle, which is what the
stable badge is a promise about.
The proximity() generic is renamed
as_proximity(). e2tree, the package most
likely to be attached alongside this one, exports a
proximity() of its own, and two generics of one name mask
each other in whichever order the user happened to call
library(). The as_ prefix also says what the
function does: it coerces a fitted ensemble into an object of class
proximity, which keeps its name. Nothing was ever released
under the old one, so there is no deprecated alias.
as_proximity() gains a method for
e2tree fits. An e2tree is one tree, so its
proximity is an indicator rather than a proportion, and it is still a
Gram matrix over the leaf indicators: positive semi-definite, of rank
the number of leaves. This closes a loop. e2tree() takes as
its D argument the dissimilarity this package induces, and
reading its result back as a proximity makes the explanation comparable
with the ensemble it explains, on the same footing, through
mantel_test(). That is the question e2tree
exists to answer and the one thing it cannot ask of itself.
loans is a new dataset: 2,000 synthetic consumer
loans over four origination years, with the outcome, ten predictors, and
one protected attribute the model is not allowed to see. It is generated
rather than collected, because the public credit panels that can be
redistributed have had exactly that attribute stripped out, and the
fairness question the case study asks cannot be asked without one.
inst/data-raw/loans.R is the specification as code: the
default rate is 16.4 per cent and rises from 13.8 per cent in the 2019
vintage to 19.9 in 2022, utilisation carries three times the weight in
the last vintage that it carries in the first, and
applicant_group shifts the median income by 8,902 and the
median score by 29 points while affecting the default probability
nowhere.
vignette("credit-scoring-case") is no longer a
placeholder. Its three questions now have measured answers: the forest
isolates a segment of about sixty borrowers that defaults at close to
nine in ten against a book rate of one in six, and that segment is the
score-by-utilisation interaction rather than either margin; four refits
of one window agree at a Mantel correlation of 0.991 while the early and
late windows applied to the same borrowers agree at 0.86; and the
protected attribute the model never received accounts for 0.4 per cent
of the variation in the proximity after the outcome, at
p = 0.002.
The README is rewritten as the package’s illustrative page, with every output in it produced by the code above it.
autoplot() is implemented, and there are five
methods rather than one. A proximity object gets
"heatmap", "mds" and "network"; a
proximity_sparse gets the network at the threshold it was
built with; a proximity_nystrom gets its configuration,
read off the stored factor; and the two objects of the stability layer
get the views that were promised of the wrong object. Reaching for
autoplot() on any class in the package now gives a picture
rather than “no applicable method”.
Two of the four views this page promised could not be drawn from
the object they were promised of. type = "stability" was to
show the agreement between replicates against the number of trees: a
proximity object has one number of trees and no replicates,
so that view is now autoplot() of a
stability() or an n_trees_required() result,
which are the objects that hold the numbers. Its bootstrap band went for
the reason it went from stability(), that the pairwise
agreements are dependent. And type = "mds" was to be
“coloured by class and by out-of-bag error”, neither of which the object
receives; colour is now the caller’s to fill, and a wrong
length is refused.
n_trees_required() returns an integer of class
proximity_trees, which buys it a print() and
the plot of its search path. The class does not survive arithmetic:
measured, R keeps the attributes of the first operand, so without
Ops.proximity_trees the number required + 100
would carry the path and eps of a search that
stopped at required and would print an answer that search
never gave.
n_trees_required() documented an integer and
returned one only when the target was missed. 2L^k is a
double in R however integer its operands, so the grid the answer is
taken off was a double vector and so was every answer read from it;
NA_integer_ was the only integer it ever returned.
tree_grid() now builds the grid as integers.
seriation and igraph are suggested, not
required. The heatmap needs the first and the network the second, each
refused by name when absent; the other three views need neither. The
refusal is tested with the package masked rather than assumed absent,
since both are installed wherever this is developed.
Both randomised views put back the random stream they found, so
drawing a plot cannot move the permutation tests around it. The network
was the expected case, through its force-directed layout and its
community search. The heatmap was not: the seriation was written up as
deterministic and measured otherwise, because a proximity is
k/B and so takes at most B + 1 distinct values
however many pairs it has. Over the eighteen cells of
inst/simulations/view-cost.R, the dissimilarity took 43
distinct values across 19,900 pairs at 50 trees and 370 across 319,600
at 500; the ordering drew from the stream in seventeen of the eighteen
and came back different under a second seed in seven. With the ties
separated by a jitter below 1/B it came back the same in
all eighteen, though it still drew from the stream. The first version of
this entry claimed it stopped drawing too, which was one ad hoc check
rather than the study.
inst/simulations/view-cost.R is new, and measures
what the heatmap costs. Median of three draws on forests of 200 trees:
at n = 800 the seriation takes 0.08 s, drawing takes 0.25 s and the plot
object is 10.2 MB; at n = 3200, 0.46 s, 3.03 s and 156.7 MB. The
seriation is not the expensive part, which is worth knowing because it
looks like it should be. The drawing is, and it grows as
n^2 with the matrix.
reject_unused() now takes the advice at the end of
its message from the caller, and autoplot() uses it:
color for colour would otherwise have drawn an
uncoloured plot and said nothing about why.
nystrom() is implemented: the approximation
P ~ P[, m] P[m, m]^-1 P[m, ] from m landmark
observations, stored in factored form as an n by
r matrix so that the n by n
object is never built, on the way in or on the way out. At n = 10,000
with 500 landmarks that is 40 MB instead of 760. The landmark block is
inverted through a pseudo-inverse with a rank cut rather than
solve(), because it is singular as soon as two landmarks
reach the same leaves in every tree.
The entries of the approximation move much more than the geometry
does, and the geometry is what the object is for. Measured over 40
replicates per cell in
inst/simulations/scalability-stability.R, at n = 400 with
200 landmarks the relative Frobenius error of the entries was 0.207
while the configuration embedding() returns was 0.049 out,
a factor of four. At n = 800 with 200 landmarks, 0.315 against 0.102. On
data whose response is independent of every predictor the same figures
run from 0.87 to 0.46: there is no low-rank structure to find and no
number of landmarks invents one.
Stratifying the landmark sample is worth less than it sounds. On
a minority class holding 2.4 per cent of the sample, a simple sample of
twenty landmarks drew none of it in 62.5 per cent of draws and
reconstructed those rows 4 per cent worse than a stratified sample did,
winning in 80 per cent of draws. By eighty landmarks the two were
indistinguishable, and the error over the whole matrix barely moved in
either case. strata earns its place at small landmark
budgets and rare classes, and nowhere else.
sparsify() is implemented, returning an object of
class proximity_sparse rather than a
proximity. The class is separate for two reasons, both
measured. Attaching the S3 class to a Matrix object
succeeds and destroys its S4 class along with every method defined on
it, so the sparse matrix has to live in a slot rather than carry the
class. And the saving does not survive being used: across every sample
size tried, sqrt(1 - P) had between 99.5 and 99.9 per cent
of its entries non-zero and its doubly centred form had 100 per cent,
exactly. Every inference function refuses a
proximity_sparse and names the reason.
The documentation of sparsify() used to say that
most pairs never share a leaf and that the thresholded matrix is
typically very sparse. Measured on forests of 500 trees, the exact zeros
are 14.3 per cent of pairs at n = 200 and 58.5 per cent at n = 1600
in-bag. The claim was wrong at the sizes users work at and is replaced
by the table.
embedding() is new, and exported: the classical
scaling configuration of a proximity, with a method that reads it off a
Nystrom factor at O(n r^2) instead of O(n^3).
Classical scaling of sqrt(1 - P) and the principal
components of P read as a kernel are the same computation
whenever the diagonal is one, which is the identity the factored method
rests on. protest() now goes through it and accepts a
proximity_nystrom directly; its agreement with
vegan::protest() on r and m^2 to
seven decimals is unchanged by the refactor.
stability() is implemented. Its documented
description said it refits the ensemble while its signature only ever
received a list of already-computed replicates; the signature won, and
the description was rewritten. The word bootstrap went with it: the
R(R-1)/2 pairwise agreements are dependent, since every
replicate enters R-1 of them, so the reported interval is
the percentile interval of those agreements and is not a sampling
distribution of anything.
The two agreement statistics do not measure the same thing, by
enough that the choice matters. Over 20 replications of four replicates
each at n = 300: on forests refitted from the same ensemble, 0.978 for
Mantel and 0.991 for the centred kernel alignment; on forests fitted to
different predictors with the same response, 0.043 and 0.959; on forests
with nothing in common at all, 0.000 and 0.670. The alignment has a
floor near two thirds and cannot separate replicates of one ensemble
from forests that share only a response. "mantel" is the
default and the recommendation is in ?stability.
n_trees_required() is implemented, and both its
criterion and its default changed because the documented ones could not
be met. CV(B) = 0.01 with max_trees = 2000
needs on the order of fifty thousand trees. The criterion is now the
mean between-replicate standard deviation over pairs divided by the mean
proximity, rather than a coefficient of variation formed pair by pair:
8.2 per cent of pairs are zero in every replicate at every block size
measured, so the per-pair form is zero over zero and the answer would
depend on which pairs were discarded.
The default eps is 0.15, which is the tightest value
tried that was reached in all eight cells of four generating processes
by two engines, at a median of 300 trees. It is measured rather than
chosen. The level of the criterion depends on the data far more than on
the ensemble: a factor of six between the cleanest process and pure
noise at every block size, and almost nothing between
randomForest and ranger on the same
data.
CV(B) falls as a power of B, and the
exponent is not one half. Fitted across the eight cells the slope was
between -0.55 and -0.60, mean -0.56, with an R^2 of at
least 0.98 everywhere. The projection reported when no block size
reaches the target fits that slope rather than assuming it; assuming one
half would overstate the requirement by half again over one
decade.
Replicates for n_trees_required() come from disjoint
blocks of the fitted ensemble’s terminal node matrix, not from
refitting. Against four real refits, four blocks of one 400-tree forest
gave a median coefficient of variation of 0.4576 against 0.4574, so the
cheaper route measures the same thing at the cost of one fit instead of
a grid of them.
mantel_test() is implemented: the correlation
between the off-diagonal entries of two proximity matrices, Pearson or
Spearman, against a null built by permuting the rows and columns of one
of them together. Partial Mantel is supported through pxz.
Measured over 300 replicates per cell in
inst/simulations/inference-calibration.R, the rejection
rate on forests fitted to independent data was 0.053 in-bag and 0.057
out-of-bag against a nominal 0.05, and the power on forests fitted to
the same data was 1.000. Undefined pairs are used where they are defined
and counted in the result, which is what vegan::mantel()
cannot do and the reason this package has its own implementation. On
matrices that vegan can take, the two agree exactly on the
statistic.
cka() and rv_coefficient() are
implemented: the normalised Frobenius inner product of the two matrices,
double-centred for cka() and as supplied for
rv_coefficient(). Both refuse an out-of-bag proximity and
name make_psd(), because an alignment between kernels needs
two kernels. The refusal is worth its inconvenience: over 600 out-of-bag
comparisons the uncorrected alignment stayed inside [0, 1]
every time and was below the corrected one every time, by 0.089 on
average and by as much as 0.154.
permanova() is implemented: the McArdle-Anderson
partition of the induced dissimilarity across the terms of a one-sided
formula, with a permutation test on each. The sums of squares are
sequential and every term is tested against the residual of the full
model, which is what vegan::adonis2(by = "terms") does; the
two were checked to agree on Df, SumOfSqs,
R2 and F exactly, and their p-values
correlated at 1.000 over 300 replicates. Measured in
inst/simulations/permanova-protest-calibration.R over 600
replicates, the power against the response the forest was trained on was
1.000 in-bag and out-of-bag alike.
Where a term sits in the formula changes the level it is tested
at, and by enough to matter. On a term that explains nothing by
construction, at a nominal 0.05 with a Monte Carlo standard error near
0.009: 0.047 alone, 0.048 beside another null term, 0.105 placed before
a term taking a seventh of the variation, and 0.020 placed after it.
adonis2() gave the same rates on the same replicates, so
this is a property of sequential permutation testing and not of this
implementation. The rule it implies is in ?permanova: put
the terms you already believe in first and the term you are testing
last.
The out-of-bag Gower matrix is indefinite, so a term’s sum of
squares could come out negative and its R-squared could leave
[0, 1]. Over 1800 out-of-bag values neither happened, and
the range observed was 0.008 to 0.942. The risk is real in principle and
did not appear in practice, which is the same shape as the finding
behind the cka() refusal.
protest() is implemented: classical scaling of each
induced dissimilarity to k dimensions, then the symmetric
Procrustes superimposition, reported as with the PROTEST permutation
test. It agrees with vegan::protest() on the statistic and
the residual to seven decimals. Over 600 replicates the rejection rate
on forests fitted to independent data was between 0.033 and 0.050
against a nominal 0.05, at every k tried and on both
definitions of the proximity, and the power on forests fitted to the
same data was 1.000.
The Procrustes statistic is not monotone in k, which
I had assumed it was. Both configurations are rescaled to unit sum of
squares at each k, so a further dimension changes what is
compared rather than adding to it. Over 600 replicates the mean
correlation on unrelated forests rose from 0.103 at two dimensions to
0.173 at six, and on forests fitted to the same data it fell between two
and six in 78% of replicates in-bag and 98% out-of-bag.
permanova() and protest() refuse a
proximity carrying undefined pairs, rather than dropping them as
mantel_test() does. A sum of squares is a quadratic form
over the whole matrix and classical scaling needs every distance, so
there is no honest per-pair rule available to either. In the calibration
study a 25-tree out-of-bag proximity was refused in every one of 300
replicates and a 200-tree one in none.
Both gain a transform argument, "sqrt"
or "linear", passed to as_dissimilarity().
Neither could be reached otherwise: they take a proximity, and handing
them a dissimilarity would transform it twice.
as_proximity() gains a method for
ranger fits, with the same in-bag and out-of-bag
definitions. The leaf co-occurrence engine now normalises each tree’s
leaf labels, since engines disagree on where they start counting. On
iris, the in-bag proximities of a 300-tree
randomForest and a 300-tree ranger forest
correlate at 0.998.
make_psd() projects an indefinite proximity onto the
cone of positive semi-definite matrices, by "clip",
"flip" or "shift". The out-of-bag proximity is
not a kernel, so centered kernel alignment, the RV coefficient and
kernel PCA all need this first. The correction is recorded on the object
and shown by print().
as_dissimilarity(), make_psd() and
double_centre() reached the new storage classes by coercion
or by dispatch and would have densified them without saying so. All
three now refuse, through the same guard the inference functions
use.as_proximity() raises on an unknown argument
said that sparsify() and nystrom() were
scheduled for a later release. They are in this one, and the message now
points at them.ranger forest was kept at fitting time
was duplicated between as_proximity.ranger() and the node
extractor. It lives in the extractor, which is the only place that needs
it and the one nystrom() and
n_trees_required() also go through.The ... of as_proximity() no longer
swallows arguments the method does not have.
vignettes/large-n.Rmd sketches a scalability layer with
as_proximity(rf, newdata = df, sparse = TRUE, threshold = 0.05),
and against the current code that call returned a dense matrix as though
the request had been honoured. It was the one path in the package that
answered wrongly rather than failing.
DESCRIPTION no longer claims the inference,
stability, scalability and visualisation layers in the present
indicative. It is the text CRAN and every package index display, and
four of those layers did not exist. Two now do.
as_proximity() no longer relabels the matrix stored
by randomForest. randomForest() declares
oob.prox = proximity, so a fit made with
proximity = TRUE carries the out-of-bag matrix;
the previous code returned it labelled "inbag" whenever
type = "inbag" was requested. The type is now recovered
from object$call, and the fit is refused when the call does
not settle the question or when the stored type is not the one asked
for.
print.proximity() counted missing matrix cells and
reported them as pairs, so it printed twice the number
summary() did. It now counts pairs.
as_proximity() computes the leaf co-occurrence as a
sparse Gram matrix, P = Z Z'/B, instead of looping over
trees. On a forest of 300 trees fitted to 1500 observations, the in-bag
proximity went from 13.0 s to 0.16 s. The out-of-bag case is the
elementwise quotient of two such products.
proximity_from_nodes_reference() keeps the old loop as the
yardstick the new engine is tested against.double_centre() subtracts row and column means directly
rather than forming the centring matrix and multiplying, which drops the
cost from O(n^3) to O(n^2) and preserves the
dimnames of the input.These entries record how the package was built up on the way to this release. The phases they name as scheduled are all complete in 1.0.0.
as_proximity() generic with a method for
randomForest fits, supporting the in-bag and out-of-bag
definitions. Pairs never jointly out-of-bag are NA rather
than 0.as.matrix(), as.dist(),
print() and summary() methods for the
proximity class.as_dissimilarity(), is_euclidean() and
double_centre() for the transform layer.summary() records that the in-bag proximity is positive
semi-definite while the out-of-bag proximity is not, so the induced
dissimilarity is Euclidean only in the in-bag case.