Random forests methodologies for :
Libraries we use :
As a mention, we use our own implementation of LDA and PLS from (Friedman, Hastie, and Tibshirani 2001, vol. 181, 114), PLS is optimized for univariate, see 5.1. For linear algebra optimization purposes on large reftables, the Linux version of binaries (standalone and python wheel) are statically linked with Intel’s Math Kernel Library, in order to leverage multicore and SIMD extensions on modern cpus.
There is one set of binaries, which contains a Macos/Linux/Windows (x64 only) binary for each platform. There are available within the “Releases” tab, under “Assets” section (unfold it to see the list).
This is pure command line binary, and they are no prerequisites or library dependencies in order to run it. Just download them and launch them from your terminal software of choice. The usual caveats with command line executable apply there : if you’re not proficient with the command line interface of your platform, please learn some basics or ask someone who might help you in those matters.
The standalone is part of a specialized Population Genetics graphical interface DIYABC-RF, presented in MER (Molecular Ecology Resources, Special Issue), (Collin et al. 2021).
pip install pyabcranger
- ABC Random Forest - Model choice or parameter estimation command line options
Usage:
../build/abcranger [OPTION...]
-h, --header arg Header file (default: headerRF.txt)
-r, --reftable arg Reftable file (default: reftableRF.bin)
-b, --statobs arg Statobs file (default: statobsRF.txt)
-o, --output arg Prefix output (modelchoice_out or estimparam_out by
default)
-n, --nref arg Number of samples, 0 means all (default: 0)
-m, --minnodesize arg Minimal node size. 0 means 1 for classification or
5 for regression (default: 0)
-t, --ntree arg Number of trees (default: 500)
-j, --threads arg Number of threads, 0 means all (default: 0)
-s, --seed arg Seed, generated by default (default: 0)
-c, --noisecolumns arg Number of noise columns (default: 5)
--nolinear Disable LDA for model choice or PLS for parameter
estimation
--plsmaxvar arg Percentage of maximum explained Y-variance for
retaining pls axis (default: 0.9)
--chosenscen arg Chosen scenario (mandatory for parameter
estimation)
--noob arg number of oob testing samples (mandatory for
parameter estimation)
--parameter arg name of the parameter of interest (mandatory for
parameter estimation)
-g, --groups arg Groups of models
--help Print help
--chosenscen
, --parameter
and --noob
, parameter estimation mode is selected.Example :
abcranger -t 10000 -j 8
Header, reftable and statobs files should be in the current directory.
With the option -g
(or --groups
), you may “group” your models in several groups splitted . For example if you have six models, labeled from 1 to 6 `-g “1,2,3;4,5,6”
Four files are created :
modelchoice_out.ooberror
: OOB Error rate vs number of trees (line number is the number of trees)modelchoice_out.importance
: variables importance (sorted)modelchoice_out.predictions
: votes, prediction and posterior error ratemodelchoice_out.confusion
: OOB Confusion matrix of the classifierWhen specifying the parameter (option --parameter
), one may specify simple composite parameters as division, addition or multiplication of two existing parameters. like t/N
or T1+T2
.
The --plsmaxvar
option (defaulting at 0.90) fixes the number of selected pls axes so that we get at least the specified percentage of maximum explained variance of the output. The explained variance of the output of the $m$ first axes is defined by the R-squared of the output:
$$Yvar^m = \frac{\sum{i=1}^{N}{(\hat{y}^{m}{i}-\bar{y})^2}}{\sum{i=1}^{N}{(y{i}-\hat{y})^2}}$$
where $\hat{y}^{m}$ is the output $Y$ scored by the pls for the $m$th component. So, only the $n_{comp}$ first axis are kept, and :
$$n_{comp} = \underset{Yvar^m \leq{} 0.90*Yvar^M, }{\mathop{\text{argmax}}}$$
Note that if you specify 0 as --plsmaxvar
, an “elbow” heuristic is activiated where the following condition is tested for every computed axis :
$$\frac{Yvar^{k+1}+Yvar^{k}}{2} \geq 0.99(N-k)\left(Yvar^{k+1}-Yvar^ {k}\right)$$
If this condition is true for a windows of previous axes, sized to 10% of the total possible axis, then we stop the PLS axis computation.
In practice, we find this $n{heur}$ close enough to the previous $n{comp}$ for 99%, but it isn’t guaranteed.
noob
parameterThe median global/local statistics and confidence intervals (global) measures for parameter estimation need a number of OOB samples (--noob
) to be reliable (typlially 30% of the size of the dataset is sufficient). Be aware than computing the whole set (i.e. assigning --noob
the same than for --nref
) for weights predictions (Raynal et al. 2018) could be very costly, memory and cpu-wise, if your dataset is large in number of samples, so it could be adviseable to compute them for only choose a subset of size noob
.
Example (working with the dataset in test/data
) :
abcranger -t 1000 -j 8 --parameter ra --chosenscen 1 --noob 50
Header, reftable and statobs files should be in the current directory.
Five files (or seven if pls activated) are created :
estimparam_out.ooberror
: OOB MSE rate vs number of trees (line number is the number of trees)estimparam_out.importance
: variables importance (sorted)estimparam_out.predictions
: expectation, variance and 0.05, 0.5, 0.95 quantile for predictionestimparam_out.predweights
: csv of the value/weights pairs of the prediction (for density plot)estimparam_out.oobstats
: various statistics on oob (MSE, NMSE, NMAE etc.)if pls enabled :
estimparam_out.plsvar
: variance explained by number of componentsestimparam_out.plsweights
: variable weight in the first component (sorted by absolute value)Comment When there isn’t any missing data, stages $2.1$ and $2.2$ could be replaced by $w{k}=\frac{X{k-1}^{T} y{k-1}}{\left|X{k-1}^{T} y{k-1}\right|}$ and $2.3$ by $t{k}=X{k-1}w{k}$
To get $W$ so that $T=XW$ we compute : $$\mathbf{W}=\mathbf{W}^{}\left(\widetilde{\mathbf{P}} \mathbf{W}^{\}\right)^{-1}$$
where $\widetilde{\mathbf{P}}{K \times p}=\mathbf{t}\left[p{1}, \ldots, p{K}\right]$ and $\mathbf{W}^{*}{p \times K} = [w_1, \ldots, w_K]$
This have been the subject of a proceedings in JOBIM 2020, PDF and video (in french), (Collin et al. 2020).
[^1]: The term “online” there and in the code has not the usual meaning it has, as coined in “online machine learning”. We still need the entire training data set at once. Our implementation is an “online” one not by the sequential order of the input data, but by the sequential order of computation of the trees in random forests, sequentially computed and then discarded.
[^2]: We only use the C++ Core of ranger, which is under MIT License, same as ours.