RDAforest 2.11.0

Fix: projection of out-of-bag samples in mtrySelJack() and ordinationJackknife()

Each jackknife replicate builds an ordination from its in-bag samples and places all samples on it. Earlier versions did this with predict(ord, Y, type = "sp"), which treats the columns of the distance matrix as species and so projects the distances themselves, not the double-centred squared distances that principal coordinates are made of. The projected samples did not keep the shape of the data: even in-bag samples, which should land exactly on their own scores, came out with each axis stretched or shrunk by a different factor (0.86 to 1.26 over the first eight axes on the wolves), and with covariates the axes also lost up to 7% of their correlation with the true scores. With Euclidean test data, where the true geometry is known, distances among projected samples were off by up to 17% (in-bag) and 26% (out-of-bag) of the largest distance.

Samples are now placed by Gower's (1968) add-a-point formula, which puts the in-bag samples exactly on their own scores and the out-of-bag samples where the same geometry puts them (Euclidean test data: all distances reproduced exactly). With covariates, each sample is adjusted for its own covariate values using the in-bag regression, as Condition() does for the in-bag samples. The same projection is available as the new exported function projectOrdination(), for placing new individuals on the axes of a model built without them, and, without vegan, as pcoa_project(): a base-R principal coordinates analysis of a subset of samples with all samples projected (optionally conditioned on covariates), in the units of the distances.

On the wolves (same seeds as before):

• mtrySelJack() selected MeanT_DryQ, T_Season and MaxT_WarmM; Prec_DryQ, previously just above the cutoff (proportion positive 0.52 at prop.positive.cutoff = 0.5), fell to 0.19 and is now dropped as a proxy of T_Season (r = -0.74).
• ordinationJackknife() with the same four predictors: predicted genetic PCs agree with the old ones at Procrustes r = 0.95 (per-axis r 0.78 to 1.00), and genetic offsets at r = 0.97.
• The lone wolf of the example sits in the same place (its five nearest model wolves are High Arctic either way), but the old projection had shrunk it towards the centre, to 83% of its distance.

In simulations of the mtry test (40 data sets with a driver, a proxy correlated with it at r = 0.95, and noise), the proxy was rejected at least as often as before (kept in 17% of data sets, against 23%).

ordinationJackknife(): replicates keep their own ordination (align)

Each replicate builds its ordination from its in-bag samples and projects the out-of-bag samples onto it. Earlier versions then rotated it onto the ordination of all samples over all axes before fitting anything. Because the replicate's ordination reproduces the in-bag samples' distances exactly, that rotation put the in-bag samples back almost exactly where the full ordination has them (on the wolves within 0.4-1.1% of the variance), so the replicates differed mainly in where the out-of-bag samples landed, not in their axes.

New argument align:

• "predictions" (default): the replicate keeps its own axes for everything it fits - forests, gradient forest, importances, turnover curves (which are in predictor units and are averaged without alignment) - and only its predicted gPC scores for newX are rotated onto the full ordination's top top.pcs axes (Procrustes rotation and scaling fitted on the samples) before being averaged. ojPredict() applies the same rotation to kept models.
• "top.pcs": the replicate's sample scores are rotated onto the full top top.pcs axes before fitting.
• "all": rotation over all axes before fitting, as before.

In simulations (12 landscapes; linear, U-shaped and threshold responses; IBD or barrier; 3 true and 3 random predictors; prediction at 1000 unsampled locations), accuracy of predicted adaptation was the same under all three (R2 0.660, 0.661, 0.663), and so was the ranking of true over random predictors. They differed in how much importances vary between replicates (coefficient of variation for the true predictors 0.057, 0.023 and 0.036 for all, top.pcs, predictions). On the wolves, "predictions" gave importances with CVs of 0.09-0.11 (0.05 under "all"), predicted gPCs that agree with "all" at Procrustes r = 0.99, and genetic offsets at r = 0.99.

Raising oob from 0.2 to 1/3 did not change prediction accuracy in the simulations and raised importance variation modestly (wolves, "predictions": CVs 0.10-0.15).

mtrySelJack() already used each replicate's own axes. It now also returns all.importances, the per-replicate importances behind its medians. In the simulations they varied more between replicates than ordinationJackknife() importances (CV 0.052 against 0.036), partly because they come from the higher of its two mtry values.

Novelty: which environments are too unlike the samples to predict

ordinationJackknife() and ojPredict() now decide what to predict with two distance-based tests, and an environment in newX is predicted only if it passes both. The range rule of earlier versions (extra, extra.npred) is still available, unchanged.

res <- ordinationJackknife(Y, X, newX = list(present = envc, future = envf))  # sigma = 4, nn.max = 3
res$future$novelty      # data frame: sigma and nn for every location (map them)
res$future$limits       # the limits used
res$future$goodrows     # predicted: passed every test in use

Both tests measure distance in the space of the model's predictors, counting each predictor by its median importance in the fit and taking their correlations into account. The predictors are standardised and decorrelated by symmetric (ZCA) whitening, which keeps each whitened component tied to its own predictor, and each component is weighted by its predictor's importance. (Rescaling the raw predictors by importance would do nothing: the Mahalanobis distance is unchanged by rescaling.)

• sigma (default 4): distance from the centre of the sampled environments, converted to the equivalent number of standard deviations of a single normal variable (chi-square approximation with the effective number of predictors as degrees of freedom), so sigma = 4 means the same whatever the number of predictors. With equal importances it is the ordinary Mahalanobis distance. It catches values beyond the sampled range and combinations of values the samples never showed.
• nn.max (default 3): mean distance to the nn.k (default 5) nearest distinct sampled environments, in units of their typical spacing (the median of the same quantity among the sampled environments, each leaving itself out). It catches what the centre-based distance misses when the samples are lopsided, and gaps within the sampled range. This is the idea behind the "area of applicability" of Meyer & Pebesma (2021), with the predictors also decorrelated; their data-driven limit (upper whisker of the samples' own values) is available as nn.max = "auto".
• extra, extra.npred: giving either switches to the range rule of earlier versions, with results identical to 2.10 (checked on simulated data and on the wolves). The distance tests then default to off; name sigma or nn.max as well to apply them too.
• Each limit may be Inf (test off) and may take one value per environment when newX is a list. Environments within the limits are predicted with each predictor clamped to its sampled range, as before. The sampled environments themselves (newX = NULL) are always predicted.
• A message says, per environment, how many locations were left out and by which test.
• New exported function envNovelty(X, newX, importance, nn.k) returns both measures, for mapping or for use with predict_rf() / predict_gf().
• goodrows is now always logical (it was an index vector when newX was omitted).

On the wolves (four selected predictors, same seed and forests as the 2.10 run):

Every location beyond sigma = 4 was also beyond nn.max = 3. The nearest- sample test adds the US south, which is hotter in summer than anywhere a wolf was sampled (the samples sit at the cool end of that gradient, so the distance from the centre alone was lenient there), north-eastern Greenland, and a few patches along the Pacific and Alaskan coasts and in the interior. Predictions at locations kept by both the new and the old rules are identical.

RDAforest 2.10.0

lfmmClean(): removing neutral genetic structure

New function implementing the procedure that came out best in simulation: fit a ridge latent factor mixed model (LFMM; Caye et al. 2019) to the principal coordinates of the genetic distances, with the leading four environmental PCs plus their squares and pairwise products as explanatory variables, choose the number of latent factors K by the broken-stick rule on the residual eigenvalues, and subtract the fitted latent structure.

cl <- lfmmClean(cordist, env)          # all candidate predictors
cl                                     # K, variance removed
mm <- mtrySelJack(cl$dist, env, ...)   # no covariates needed
oj <- ordinationJackknife(cl$dist, env[, mm$goodvars], ...)

In simulations (200 samples; isolation by distance with or without a barrier; linear, U-shaped, threshold, plateau and asymmetric responses, additive or interacting), it raised the accuracy of predicted adaptation at unsampled locations from R2 = 0.49 to 0.65 (from 0.52 to 0.63 in a second set with threshold-type responses), was best or near-best in every setting, and was never worse than no correction. Subtracting the fitted structure beat conditioning the ordination on the latent factors; quadratic terms beat linear ones by 0.12 and matched spline terms.

When the leading genetic structure is itself explained by the environmental terms, the broken-stick rule finds K = 0 and the distances are returned unchanged, with a message. This is the case for the wolves example, where the environmental terms explain 59-93% of the three leading genetic axes.

The ridge LFMM estimator is ported from the lfmm package (GPL-3); its authors, Kevin Caye and Olivier Francois, are credited as contributors.

RDAforest 2.9.10

Changes since 2.8.1.

ordinationJackknife(): several environments in one run

newX may now be a named list of environments, for example newX = list(present = envc, future = envf), with extra and extra.npred given once or once per environment (extra = c(0.1, 0.2)). Each replicate fits its models, scores every environment, and lets the models go, so all environments are scored through the same replicates while the ensemble is never held in memory. The result is a list of ordinary results under the same names (res$present, res$future); unnamed elements are called newX1, newX2, ..., and a NULL element stands for X. With a single data frame, or with newX omitted (which scores X), the result is exactly as before, bit for bit.

This matters for genetic offsets. An offset is the difference between two prediction sets, and if the two come from separate runs that difference also carries the sampling noise of two independent sets of forests; on the wolves that noise was 20-50% of the offset itself. Scoring both environments in one run removes it at no extra cost: fitting dominates the running time, so two environments take about as long as one.

Keeping the models: keep.models and ojPredict()

keep.models = TRUE (default FALSE) also returns the fitted models of every replicate in $models (once, next to the environments, when newX is a list), and the new ojPredict() scores further environments with them, without refitting:

res <- ordinationJackknife(Y, X, newX = envc, extra = 0.1, keep.models = TRUE)
fut <- ojPredict(res, newX = envf, extra = 0.2, ncores = 4)
gen_offset_oj(res, fut, envc[, 1:2], envf[, 1:2])

ojPredict() returns exactly what ordinationJackknife() would have returned had those environments been in its newX, with the same ensemble.id, so its predictions can be set against the original run's in gen_offset_oj(). Like ordinationJackknife(), it takes one environment, a named list of them, or none (which scores X), with extra and extra.npred once or per environment, and ncores: on the wolves, scoring the future grid through 25 replicates took 42 s on two cores instead of 74 s, with identical results.

Only what prediction reads is kept: each replicate's random forests without their training-only parts (node sums of squares and improvements, out-of-bag predictions, the response), and its gradient forest as the turnover curve of every predictor. For the wolves at the defaults that is 1.2 GB, against about 2.7 GB for the full models; a message reports the size above a gigabyte. Predictions are identical whether models are kept or not.

Every result carries ensemble.id, a token of the run it came from, and gen_offset_oj() warns when its two arguments come from different runs. Objects without a token pass unchecked.

Parallel replicates: ncores

mtrySelJack() and ordinationJackknife() take ncores (default 1). The replicates are independent, so with ncores > 1 they run that many at a time; on two cores a wolves mtrySelJack() run went from 52 s to 30 s, and the gain grows with the number of cores up to nreps.

• Forking on Linux and macOS, socket clusters on Windows; options(RDAforest.parallel = "psock") forces sockets where forking misbehaves.
• Reproducibility: with ncores = 1 nothing changes. With ncores > 1 each replicate draws from its own L'Ecuyer-CMRG stream derived from the session's stream, so set.seed() still fixes the result, and the result is identical for any number of cores above one, and between forking and sockets; it differs from the ncores = 1 result only by the luck of the draw. The session's own generator is left as it was.
• Memory: each worker needs its own memory for fitting, roughly 0.8 GB for the wolves at the defaults, so choose ncores by RAM as well as by processor count. A worker that runs out is killed by the operating system, and the function stops with "a worker was probably killed, most often for lack of memory. Try fewer cores." ordinationJackknife() dispatches ncores replicates at a time and sums their predictions as they return, so beyond the workers its memory stays flat.

mtrySelJack()

Few predictors. The test compares importances at two mtry values. With five or fewer predictors the default values coincided (3 and 3), so the test compared a forest with itself and its outcome was a coin toss: in simulations with four predictors, true drivers were kept only about half the time at prop.positive.cutoff = 0.5. The defaults are now:

Explicit values are honoured but must satisfy 1 <= mintry < maxtry <= number of predictors. The values used are returned in $mtry.

Dropped predictors correlated with retained ones. The test treats a predictor that loses importance at higher mtry as a proxy of a correlated, more important one. Two predictors that are both real drivers but correlated compete the same way, and one of them can be dropped (in simulations: 10-20% of the time at r = 0.8, 25-40% at r = 0.9). Every dropped predictor that correlates with a retained one at |r| >= report.cor (default 0.7) is now listed in a message and returned in $dropped.pairs (dropped, kept, r, importance, prop.positive, reason). report.cor = NULL switches this off. Selection itself is unchanged. The predictor correlations and the threshold are returned too ($cor.X, $report.cor).

Reselect() now recomputes and reports dropped.pairs for the new selection (it used to leave the old list, which could then name as dropped a predictor that was now kept). It takes report.cor (by default the value used in mtrySelJack()), and X = for objects made by earlier versions, which do not store the correlations. It also applies the importance cutoff exactly as mtrySelJack() does (>=; it was >) and returns goodvars in the same order, so that with the original cutoffs it reproduces the original selection exactly.

Fixes

• A distance matrix without row names, such as 1 - cor(t(genotypes)), made mtrySelJack() and ordinationJackknife() stop with vegan's "'newdata' does not have named rows matching one or more of the original rows": each replicate's projection matches samples by name. Unnamed matrices are now given names 1..n.
• A single predictor, or top.pcs = 1, failed in predict_rf(), predict_gf(), mtrySelJack() and ordinationJackknife(), because a one-column subset collapsed to a vector; and ordinationJackknife()'s default mtry of 2 exceeded the number of predictors. Both fixed.

Distance matrices only

mtrySelJack() and ordinationJackknife() now accept only a distance matrix as Y (square and symmetric, or a dist object). Raw data matrices are refused with an error suggesting as.matrix(dist(Y)), which gives the same ordination as an RDA of Y. This also removes a failure: a raw matrix together with covariates stopped with vegan's "no 'wa' scores available (yet) in partial RDA". Results for distance input are unchanged.

RDAforest 2.8.1

extra.npred in ordinationJackknife()

extra says how far a predictor in newX may reach beyond the range the model was fitted on before its row is dropped. It now applies only to the first extra.npred predictors, counted in the order they were supplied in X; the rest may extend without limit, so no row is ever dropped on their account.

ordinationJackknife(Y, X, newX = rasters, extra = 0.1, extra.npred = 2)

With five predictors, that lets predictors 1 and 2 reach a tenth of their span past the modelled range and no further, while predictors 3 to 5 are unrestricted. Useful when some predictors should be extrapolated cautiously (climate, say) while others legitimately cover ground the samples did not (coordinates, depth).

The default is extra.npred = ncol(X) which, with extra = 0, is exactly the old behaviour: no extension anywhere.

Fixes

• Automatic site detection no longer overrides an explicitly named site.repeats. Detection supplies the sites; the mode is whatever the caller asked for, and only a caller who left site.repeats alone gets the "genetic.median" default. Affected rdaf_randomForest(), rdaf_gradientForest(), makeGF(), makeGF_simple(), predict_rf(), mtrySelJack() and ordinationJackknife(). When a mode is named, the message now says which one was used instead of naming the genetic median.

RDAforest 2.8.0

site.repeats: three ways to handle repeated samples within a site

The oob.blocks / oob.block.holdout pair introduced in 2.7.0 is replaced by sites plus site.repeats. Pass the site membership of each sample as sites, and choose what is done about the repeats:

sites = NULL (the default) means no repeated sampling and leaves everything as it was, so code that does not use the option is unaffected.

"genetic.median"

Each site is replaced by a single "false individual": not one of the real samples, but the point with the smallest sum of Euclidean distances to all of them, found by Weiszfeld's iteration in the space of all the principal coordinates being analysed. The iteration starts at the centroid and repeatedly replaces the estimate by the mean of the points weighted by the reciprocal of their distances to it, until it stops moving; if it lands exactly on a data point that point is returned.

This is not the medoid (which is a real sample, and so carries that individual's own noise) and not the centroid (which minimises squared distances and is pulled around by outliers). The forest is then fitted to one row per site, so there is no repeated sampling left to exploit.

New exported helpers: geometric_median() for a single set of points, genetic_medians() for the per-site collapse, and detect_sites() for the automatic check described below.

"blocked.resampling"

What oob.blocks did in 2.7.0. Each tree withholds a random block.holdout fraction of whole sites (default 1/3), then draws its training sample from the rows of the sites that remain, in exactly the ordinary way: with replacement by default, honouring replace and sampsize (the latter scaled to the size of the pool). The out-of-bag set is the withheld sites and nothing else.

"direct.use"

The ordinary random forest, bit for bit: what RDAforest did before 2.7.0.

Automatic detection

sites no longer has to be supplied. Left NULL, RDAforest inspects the predictors before fitting. Coordinate columns are recognised by name (lon/lat, x/y, easting/northing and similar) or given with coords; with no coordinates identified, a site is a set of rows whose predictor values are identical throughout.

• No shared locations — the individuals are used directly, silently.
• Shared location, same environment — treated as repeated sampling and handled according to site.repeats. With site.repeats left at its default the message is "Data appear to contain multiple samples per site. Such samples will be represented as their genetic median; see option site.repeats for alternatives."; naming a mode overrides the default, and the message then says which mode was used. Detection supplies the sites, never the mode.
• Shared location, different environment — "some samples share a spatial location but have different environmental entries, check" is warned, and the individuals are used directly, leaving the decision to you.

Telling the second case from the third needs the coordinates to be identifiable. With no recognisable column names and no coords, a site can only be a set of identical rows, so a sample that shares a location but differs in one environmental value counts as a site of its own and nothing is warned about.

detect_sites() runs the check on its own and returns what it found. auto.sites = FALSE switches it off; passing sites overrides it. The check runs once per top-level call, so a jackknife does not repeat the message on every replicate.

Why it matters

In a simulation of 40 sites with 10 individuals each, two of four site-level predictors truly driving the genetics at R² = 0.6:

Under "direct.use" a predictor that merely tags an unusual site looks like a driver, and R² overstates a relationship whose truth is 0.6. The two honest settings agree with each other and now understate it, because they are working from 40 effective observations rather than 400.

Where the option lives

rdaf_randomForest(), rdaf_gradientForest(), makeGF(), makeGF_simple(), predict_rf(), predict_gf(), mtrySelJack() and ordinationJackknife().

site_repeats_summary() reports what a given setting implies (how many sites, how many rows will be fitted, and for blocked resampling how large the eligible pool is per tree) before you commit to a long run. It replaces oob_blocks_summary().

spatial_blocks() no longer fails when a re-split lands on a block too small to halve (hclust needs at least two objects), which could happen with a small min.size and many blocks.

As before, the oob argument of ordinationJackknife() and mtrySelJack(), which governs how many samples are withheld when the ordination is rebuilt in each replicate, remains a plain random draw over rows. site.repeats governs the random forest fitting only.

Deprecations

oob.blocks and oob.block.holdout are still accepted by rdaf_randomForest() and rdaf_gradientForest(), with a warning; they map to sites with site.repeats = "blocked.resampling" and to block.holdout. They will be removed in a future version. oob_blocks_summary() is gone; use site_repeats_summary().

RDAforest 2.7.0

Blocked out-of-bag sampling

Random forest fitting can now hold out whole pre-specified blocks of samples rather than whichever individual rows the bootstrap happens to miss.

Pass a categorical vector, one entry per sample, as oob.blocks. oob.block.holdout sets the fraction of blocks withheld from each tree (default 1/3, rounded to a whole number of blocks).

Blocking changes only which samples are eligible; nothing else. Each tree first withholds a random fraction of the blocks, then draws its training sample from the rows of the blocks that remain, in exactly the ordinary way: with replacement by default, so about 1 − 1/e of the eligible rows are drawn, and honouring replace and sampsize (the latter scaled to the size of the pool). The out-of-bag set is the withheld blocks and nothing else — a retained row that the draw happened to miss is not scored, because its block-mates are in the training sample.

# 12 geographic blocks of at least 4 samples each
blocks <- spatial_blocks(latlon, nblocks = 12, min.size = 4)
oob_blocks_summary(blocks)

oj <- ordinationJackknife(Y = cordist, X = env, newX = envc,
                          covariates = latlon.gcd, nreps = 25,
                          oob.blocks = blocks)

Why it matters: when samples come in groups (several individuals per site, siblings, sequencing batches, repeated years), a row the bootstrap missed almost always has a near-duplicate of itself in the training set. Out-of-bag R-squared and variable importance are then optimistic, and predictors that merely tag the group look important. Blocked sampling scores every tree on groups it has never seen.

The option is available on rdaf_randomForest(), rdaf_gradientForest(), makeGF(), makeGF_simple(), predict_rf(), predict_gf(), mtrySelJack() and ordinationJackknife(). It defaults to NULL, which is the ordinary bootstrap, so existing scripts are unaffected.

Note: the oob argument of ordinationJackknife() and mtrySelJack(), which governs how many samples are withheld when the ordination is rebuilt in each replicate, is still a plain random draw over rows. oob.blocks changes the random forest fitting only.

New functions

• spatial_blocks() builds a blocking vector by clustering samples on great circle distances (Ward's method), merging and re-splitting until every block holds at least min.size samples.
• oob_blocks_summary() reports what a blocking scheme implies: number of blocks, block sizes, blocks withheld per tree, and the resulting range of eligible-pool sizes.

The package is now self-contained

extendedForest and gradientForest are no longer dependencies. The functions RDAforest used from them are incorporated into the package, together with the C and Fortran random forest engine, which is what made blocked sampling possible.

Incorporated functions are renamed with an rdaf_ prefix so that RDAforest can be loaded alongside the original packages without masking anything:

Fitted objects have class rdaf_randomForest / rdaf_gradientForest, and print(), plot(), predict() and density() dispatch on them as before.

Everything else is unchanged. With the same random seed and no oob.blocks, the incorporated code reproduces extendedForest and gradientForest exactly: the package's test suite asserts identical() on the forest structure, MSE, R-squared, out-of-bag predictions, in-bag matrix, importances, cumimp() curves and predict() output.

The parts of gradientForest that RDAforest never used — the combinedGradientForest family — were not carried over.

Consequences for installation

RDAforest now contains compiled code, so installing from source needs a C compiler: Rtools on Windows, the Xcode command line tools on macOS, r-base-dev on Linux.

No Fortran compiler is required. The classification tree builder that extendedForest inherits from Breiman and Cutler (rfsub.f: buildtree, findbestsplit, movedata and their helpers) has been translated to C in src/rfsub.c, and gfortran is no longer part of the build. The translation is deliberately literal, keeping the original column-major layout, control flow and — critically — the exact order in which random numbers are drawn, including two tie-breaking quirks that look like bugs but shift every subsequent draw. The test suite pins this down: with the same seed, classification forests are identical() to the Fortran ones across fourteen configurations covering multiple classes, mtry, nodesize, sampling with and without replacement, class weights and cutoffs, stratified sampling, proximity, local importance, conditional permutation, categorical predictors (both the exhaustive and the sampled-split branches) and a deliberately tie-heavy dataset.

If you call the old functions directly

Scripts that only use RDAforest's own functions need no changes. A script that called gradientForest() or randomForest() directly should either add the rdaf_ prefix, or keep loading the original packages, which still works.

Documentation

Every incorporated function now has its own help page, and the new blocked sampling option is documented on each function that accepts it. ?RDAforest gives an overview.