--- title: "Comparing forests through their proximities" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Comparing forests through their proximities} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") has_rf <- requireNamespace("randomForest", quietly = TRUE) ``` Two forests fitted on the same data with different hyperparameters represent that data differently. The question "how differently?" has no answer in the usual toolkit, because the objects being compared are $n \times n$ matrices rather than vectors of predictions. ## Two forests to compare A shallow forest and a deep one, on the same rows. The shallow one is capped at four leaves per tree, so it can only carve the data into coarse blocks; the deep one is free to isolate individual observations. ```{r setup, eval = has_rf} library(Proximum) library(randomForest) set.seed(1) rows <- sample(nrow(iris), 60) df <- iris[rows, ] set.seed(2) rf_shallow <- randomForest(Species ~ ., data = df, ntree = 300, maxnodes = 4, keep.inbag = TRUE) set.seed(3) rf_deep <- randomForest(Species ~ ., data = df, ntree = 300, keep.inbag = TRUE) px1 <- as_proximity(rf_shallow, newdata = df) px2 <- as_proximity(rf_deep, newdata = df) px1 ``` ## Correlation between two proximity matrices The Mantel test correlates the off-diagonal entries of two matrices and gets its null distribution by permuting the rows and columns of one of them together, which preserves the dependence induced by the matrix structure as an ordinary permutation of the entries would not. That is not a technicality. Each entry of a proximity matrix shares an observation with $2(n-2)$ others, so the $n(n-1)/2$ entries are a long way from independent. Referring the correlation to the usual table would reject almost whatever you fed it. ```{r mantel, eval = has_rf} mantel_test(px1, px2, n_perm = 999) ``` The result carries the number of pairs it used. On an in-bag proximity that is every pair; on an out-of-bag one it can be fewer, because a pair that was never jointly out of bag has no proximity to correlate. The test uses what is defined and reports how much that was. ```{r oob, eval = has_rf} oob <- as_proximity(rf_deep, newdata = df, type = "oob") mantel_test(px1, oob, n_perm = 999)$parameter ``` ## Alignment of representations `cka()` treats each proximity matrix as a kernel and measures the alignment of the two implied feature spaces. It is invariant to isotropic rescaling and to orthogonal transformation, which is why it has become the standard way of comparing learned representations in the deep learning literature. Applying it to forests puts ensemble explainability and representation similarity on the same footing. It is a coefficient and not a test: there is no p-value, and a high alignment is not evidence of anything by itself. `mantel_test()` is where the evidence is. ```{r cka, eval = has_rf} cka(px1, px2) rv_coefficient(px1, px2) ``` The two differ only in the centring. `cka()` double-centres each matrix first, which is what removes the mean similarity any two kernels on the same observations share whether or not they have learned the same structure. ### Out-of-bag matrices have to be repaired first An out-of-bag proximity is not positive semi-definite, and no number of trees repairs it: its entries are ratios whose denominators count only the trees where each pair was jointly out of bag, and those denominators differ across pairs. An alignment computed on it is not an alignment between kernels, so `cka()` refuses it. The refusal is not about the number leaving `[0, 1]`. Over 600 out-of-bag comparisons it never did. It is that the number is quietly too low: the negative eigenvalues subtract from the numerator, so the uncorrected alignment understated the corrected one in every one of those 600 comparisons, by 0.089 on average and by as much as 0.154. ```{r refuse, error = TRUE, eval = has_rf} cka(oob, oob) ``` `make_psd()` is the repair, and it belongs to you rather than to `cka()`, because which correction to use and what it costs is a decision about your data. ```{r repair, eval = has_rf} cka(make_psd(px2), make_psd(oob)) ``` ## Superimposing the two configurations `protest()` asks the same question geometrically. Each matrix is reduced to a configuration of `k` dimensions by classical scaling, and the two are superimposed by the best translation, rotation, reflection and rescaling. What is left over is the residual $m^2$, and the reported statistic is $r = \sqrt{1 - m^2}$, so a larger value is more agreement. ```{r protest, eval = has_rf} protest(px1, px2, n_perm = 999) ``` The null permutes the rows of the second configuration. That is the same null as permuting the second proximity matrix and scaling it again, because classical scaling commutes with relabelling, and it costs one permutation of a small matrix rather than one eigendecomposition of a large one. `k` is part of the question. The residual is not monotone in it, because both configurations are rescaled to unit sum of squares before the fit, so a further dimension changes what is being compared rather than adding to it. ```{r protest_k, eval = has_rf} vapply(c(2, 4, 6), function(k) unname(protest(px1, px2, k = k, n_perm = 99)$statistic), numeric(1)) ``` Measured over 600 replicates, on forests fitted to unrelated data the mean correlation rose by half again between two dimensions and six, so a `k` chosen after seeing the answer is a `k` chosen to flatter it. Fix it first. The level holds across the range either way: on independent data the test rejected between 0.033 and 0.050 of the time at a nominal 0.05, at every `k` tried. ## Partitioning one matrix The two functions above compare two matrices. `permanova()` takes one apart, asking how much of the structure the forest learned is explained by the response and how much by covariates the model never saw. ```{r permanova, eval = has_rf} permanova(px2, ~ Species + Sepal.Width, data = df, n_perm = 999) ``` `Species` is what the forest was trained on, so it takes most of the variation and the p-value is the smallest the permutation count allows. `Sepal.Width` is a predictor the forest did see, and what it gets here is what it adds after `Species`, not what it would take on its own. ### The order of the terms is part of the question 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. That has a consequence worth stating plainly, because it decides how you write the formula. A permutation destroys the whole matrix, the other terms' contribution included. So the observed pseudo-F of an early term is divided by a residual that a strong later term has already shrunk, while its permuted values are divided by residuals nothing has shrunk. Measured over 600 replicates on a term that explains nothing by construction, at a nominal level of 0.05: | The model | Rejection rate | |---|---| | the null term alone | 0.047 | | beside another null term | 0.048 | | before a term taking a seventh of the variation | 0.105 | | after that same term | 0.020 | `adonis2()` was measured on the same replicates and gave the same rates, so this is what sequential permutation testing is, not a defect of this implementation. The rule that follows is short: put the terms you already believe in first, and the term you are testing last. ### Undefined pairs are refused rather than dropped `mantel_test()` can drop a pair that was never jointly out of bag, because a correlation is a sum over pairs. A sum of squares is a quadratic form over the whole matrix, and classical scaling needs every distance, so `permanova()` and `protest()` refuse such a matrix instead. In the calibration study, a 25-tree out-of-bag proximity was refused in every one of 300 replicates; at 200 trees it was refused in none. ```{r refuse_undefined, error = TRUE, eval = has_rf} small <- randomForest(Species ~ ., data = df, ntree = 15, keep.inbag = TRUE) permanova(as_proximity(small, newdata = df, type = "oob"), ~ Species, data = df) ``` ## Still to come Phase F2 is complete. What is not is the motivating case underneath all of it: `e2tree` approximates a forest with a single tree, and whether that tree preserves the forest's view of the data is exactly a question about two proximity matrices. That comparison needs a `as_proximity()` method for `e2tree` objects, which does not exist yet.