SIRUS: Stable and Interpretable RUle Set for Classification

TL;DR

SIRUS is a new classification method that produces simple, stable, and interpretable rule sets with accuracy comparable to random forests, suitable for critical decision-making contexts.

Contribution

The paper introduces SIRUS, a stable and interpretable rule-based classifier derived from random forests, addressing the interpretability and stability issues of existing models.

Findings

SIRUS achieves higher stability than existing methods.

SIRUS maintains predictive accuracy close to random forests.

The approach is validated through extensive experiments.

Abstract

State-of-the-art learning algorithms, such as random forests or neural networks, are often qualified as "black-boxes" because of the high number and complexity of operations involved in their prediction mechanism. This lack of interpretability is a strong limitation for applications involving critical decisions, typically the analysis of production processes in the manufacturing industry. In such critical contexts, models have to be interpretable, i.e., simple, stable, and predictive. To address this issue, we design SIRUS (Stable and Interpretable RUle Set), a new classification algorithm based on random forests, which takes the form of a short list of rules. While simple models are usually unstable with respect to data perturbation, SIRUS achieves a remarkable stability improvement over cutting-edge methods. Furthermore, SIRUS inherits a predictive accuracy close to random forests, combined with the simplicity of decision trees. These properties are assessed both from a theoretical and empirical point of view, through extensive numerical experiments based on our R/C++ software implementation sirus available from CRAN.

Figures

Tables

Table 1. Table 1: Description of UCI datasets

Dataset	Sample size	Total number of variables	Number of categorical variables
Authentification	1372	4	0
Breast Wisconsin	699	9	9
Credit Approval	690	15	9
Credit German	1000	20	13
Diabetes	768	8	0
Haberman	306	3	0
Heart C2	303	13	7
Heart H2	294	13	7
Heart Statlog	270	13	3
Hepatitis	155	19	0
Ionosphere	351	33	0
Kr vs Kp	3196	36	36
Liver Disorders	345	6	0
Mushrooms	8124	21	21
SECOM	1567	590	0
Sonar	208	60	0
Spambase	4601	57	0
Titanic	887	6	1
Vote	435	16	16
Wilt	4339	5	0

Table 2. Table 2: Mean model size over a 10 10 10 -fold cross-validation for UCI datasets. Results are averaged over 10 10 10 repetitions of the cross-validation.

	Decision tree	Classical rule learning	Frequent pattern mining		Tree ensemble
Dataset	CART	RIPPER	CBA	BRL	RuleFit	Node harvest	SIRUS
Authentification	21	7	7	17	49	30	13
Breast Wisconsin	7	12	24	7	24	32	24
Credit Approval	5	4	55	4	15	27	16
Credit German	18	3	69	4	33	33	20
Diabetes	13	3	17	6	26	31	8
Haberman	2	1	2	2	3	17	5
Heart C2	10	3	34	4	23	36	20
Heart H2	5	2	29	3	12	24	12
Heart Statlog	10	3	27	4	22	35	16
Hepatitis	2	2	14	2	8	14	12
Ionosphere	4	4	38	4	20	35	15
Kr vs Kp	16	15	29	9	18	13	24
Liver Disorders	15	3	2	3	19	33	17
Mushrooms	4	8	25	11	10	22	23
Sonar	6	4	33	2	32	83	19
Spambase	14	16	126	16	68	60	21
Titanic	13	4	4	3	19	23	6
Vote	2	2	25	NA	12	10	7
Wilt	9	5	3	10	31	19	24

Table 3. Table 3: Mean stability over a 10 10 10 -fold cross-validation for UCI datasets. Results are averaged over 10 10 10 repetitions of the cross-validation. Values within 10% of the maximum are displayed in bold. Algorithms are ranked with a Mann-Whitney-Wilcoxon test, the p-value with the previous performing algorithm determines the final rank (10%-level test).

	Decision tree	Classical rule learning	Frequent pattern mining		Tree ensemble
Dataset	CART	RIPPER	CBA	BRL	RuleFit	Node harvest	SIRUS
Authentification	0.41	0.36	0.87	0.86	0.48	0.59	0.81
Breast Wisconsin	0.21	0.55	0.80	0.78	0.34	0.71	0.70
Credit Approval	0.52	0.32	0.43	0.52	0.25	0.23	0.75
Credit German	0.46	0.22	0.51	0.41	0.24	0.48	0.66
Diabetes	0.29	0.21	0.46	0.73	0.39	0.45	0.81
Haberman	0.83	0.09	0.79	0.50	0.46	0.52	0.65
Heart C2	0.25	0.35	0.38	0.60	0.39	0.49	0.71
Heart H2	0.46	0.27	0.52	0.73	0.29	0.29	0.65
Heart Statlog	0.30	0.41	0.41	0.75	0.35	0.48	0.83
Hepatitis	0.26	0.16	0.24	0.34	0.26	0.49	0.68
Ionosphere	0.96	0.39	0.13	0.70	0.17	0.33	0.69
Kr vs Kp	0.71	0.74	0.84	0.80	0.19	0.27	0.87
Liver Disorders	0.23	0.10	0.91	0.50	0.24	0.31	0.58
Mushrooms	1	0.84	0.98	0.80	0.69	0.48	0.86
Sonar	0.34	0.04	0.09	0.19	0.09	0.20	0.55
Spambase	0.49	0.10	0.46	0.86	0.28	0.66	0.78
Titanic	0.55	0.42	0.69	0.88	0.37	0.36	0.76
Vote	1	0.52	0.68	NA	0.21	0.30	0.75
Wilt	0.36	0.32	0.72	0.94	0.47	0.64	0.73
Average Rank	4.2	5.9	3.3	2.8	5.6	4.3	1.9
p-values	0.07	0.33	0.33	0.08	0.05	0.98
Final Rank	4	6	2	2	6	4	1

Table 4. Table 4: Model error (1-AUC) over a 10 10 10 -fold cross-validation for UCI datasets. Results are averaged over 10 10 10 repetitions of the cross-validation. Values within 10% of the minimum are displayed in bold, random forest is put aside. Algorithms are ranked with a Mann-Whitney-Wilcoxon test, the p-value with the previous performing algorithm determines the final rank (10%-level test).

	Black box	Decision tree	Classical rule learning	Frequent pattern mining		Tree ensemble
Dataset	Random Forest	CART	RIPPER	CBA	BRL	RuleFit	Node harvest	SIRUS
Authentification	${10}^{-4}$	0.02	0.02	0.14	0.009	${9.10}^{-\mathrm{𝟒}}$	0.02	0.03
Breast Wisconsin	0.009	0.06	0.07	0.05	0.02	0.01	0.01	0.01
Credit Approval	0.07	0.14	0.15	0.14	0.11	0.08	0.07	0.09
Credit German	0.20	0.29	0.38	0.40	0.33	0.23	0.26	0.25
Diabetes	0.17	0.25	0.29	0.30	0.25	0.18	0.19	0.19
Haberman	0.31	0.48	0.39	0.50	0.43	0.37	0.34	0.35
Heart C2	0.10	0.19	0.23	0.17	0.23	0.12	0.12	0.10
Heart H2	0.11	0.23	0.24	0.24	0.16	0.11	0.11	0.12
Heart Statlog	0.10	0.20	0.21	0.17	0.22	0.12	0.12	0.10
Hepatitis	0.12	0.48	0.39	0.36	0.33	0.20	0.23	0.17
Ionosphere	0.02	0.11	0.12	0.13	0.10	0.04	0.07	0.07
Kr vs Kp	${9.10}^{-4}$	0.02	0.009	0.05	0.01	0.009	0.04	0.04
Liver Disorders	0.23	0.33	0.35	0.48	0.44	0.27	0.30	0.35
Mushrooms	0	0.007	${3.10}^{-5}$	${5.10}^{-4}$	${2.10}^{-\mathrm{𝟓}}$	${5.10}^{-4}$	0.002	${6.10}^{-4}$
Sonar	0.07	0.27	0.26	0.25	0.44	0.12	0.16	0.2
Spambase	0.01	0.11	0.08	0.12	0.05	0.02	0.04	0.07
Titanic	0.13	0.19	0.21	0.27	0.21	0.14	0.16	0.17
Vote	0.01	0.06	0.04	0.06	NA	0.02	0.02	0.02
Wilt	0.007	0.18	0.13	0.48	0.07	0.02	0.08	0.11
Average Rank		5	4.9	5.8	4.4	1.4	2.4	2.8
p-values		0.22	0.24	0.01	${6.10}^{-3}$		0.01	0.34
Final Rank		4	4	7	4	1	2	2

Table 5. Table 5: Accuracy, measured by 1-AUC on UCI datasets, for two algorithms: Breiman’s random forests and random forests with splits limited to q 𝑞 q -quantiles, for q ∈ { 2 , 5 , 10 , 20 } 𝑞 2 5 10 20 q\in\{2,5,10,20\} .

Dataset	Breiman’s RF	q=2	q=5	q=10	q=20
Authentification	0.0002	0.08	0.002	0.0005	0.0004
Diabetes	0.17	0.23	0.18	0.18	0.18
Haberman	0.32	0.35	0.30	0.32	0.30
Heart Statlog	0.10	0.10	0.10	0.10	0.10
Hepatitis	0.13	0.15	0.14	0.14	0.13
Ionosphere	0.02	0.07	0.03	0.02	0.02
Liver Disorders	0.23	0.32	0.27	0.25	0.24
Sonar	0.07	0.09	0.07	0.07	0.07
Spambase	0.01	0.14	0.03	0.02	0.01
Titanic	0.13	0.15	0.14	0.14	0.13
Wilt	0.007	0.15	0.03	0.02	0.02

Table 6. Table 6: SIRUS error (1-AUC) over a 10 10 10 -fold cross-validation (averaged over 10 10 10 repetitions) when tree depth is limited to 1 1 1 , 2 2 2 or 3 3 3 . Values within 10% of the minimum are displayed in bold, except for datasets with no significant variations.

Dataset	SIRUS - depth = 1	SIRUS - depth = 2	SIRUS - depth = 3
Authentification	0.07	0.03	0.03
Breast Wisconsin	0.01	0.01	0.01
Credit Approval	0.11	0.09	0.09
Credit German	0.25	0.25	0.26
Diabetes	0.19	0.19	0.19
Haberman	0.35	0.35	0.35
Heart C2	0.11	0.10	0.11
Heart H2	0.12	0.12	0.12
Heart Statlog	0.11	0.10	0.10
Hepatitis	0.15	0.17	0.18
Ionosphere	0.07	0.07	0.07
Kr vs Kp	0.05	0.04	0.06
Liver Disorders	0.38	0.35	0.35
Mushrooms	${3.10}^{-3}$	${6.10}^{-4}$	${3.10}^{-\mathrm{𝟒}}$
Sonar	0.19	0.2	0.2
Spambase	0.06	0.07	0.07
Titanic	0.19	0.17	0.16
Vote	0.02	0.02	0.02
Wilt	0.19	0.11	0.11

Table 7. Table 7: Mean model size over a 10 10 10 -fold cross-validation for UCI datasets (averaged over 10 10 10 repetitions).

Dataset	C5.0	PART	FOIL	SIRUS
Authentification	11	8	20	13
Breast Wisconsin	5	10	41	24
Credit Approval	9	32	40	16
Credit German	22	68	101	22
Diabetes	12	7	36	8
Haberman	2	2	4	5
Heart C2	10	20	31	20
Heart H2	4	15	29	12
Heart Statlog	10	18	28	15
Hepatitis	7	8	14	12
Ionosphere	9	6	28	15
Kr vs Kp	11	21	24	24
Liver Disorders	14	7	2	17
Mushrooms	7	9	14	23
Sonar	10	6	20	19
Spambase	29	46	73	21
Titanic	7	15	17	6
Vote	5	7	19	7
Wilt	10	8	10	24

Table 8. Table 8: Mean stability over a 10 10 10 -fold cross-validation for UCI datasets (averaged over 10 10 10 repetitions). Values within 10% of the maximum are displayed in bold.

Dataset	C5.0	PART	FOIL	SIRUS
Authentification	0.44	0.43	0.81	0.81
Breast Wisconsin	0.17	0.49	0.36	0.70
Credit Approval	0.18	0.31	0.17	0.75
Credit German	0.03	0.16	0.11	0.65
Diabetes	0.07	0.15	0.18	0.81
Haberman	0.28	0.25	0.64	0.65
Heart C2	0.09	0.15	0.16	0.71
Heart H2	0.32	0.31	0.39	0.65
Heart Statlog	0.11	0.15	0.15	0.82
Hepatitis	0.10	0.15	0.05	0.68
Ionosphere	0.24	0.13	0.07	0.69
Kr vs Kp	0.65	0.51	0.85	0.87
Liver Disorders	0.05	0.07	0.69	0.58
Mushrooms	0.79	0.78	0.93	0.86
Sonar	0.06	0.06	0.04	0.55
Spambase	0.08	0.08	0.11	0.78
Titanic	0.49	0.27	0.77	0.76
Vote	0.67	0.40	0.39	0.75
Wilt	0.34	0.37	0.48	0.73

Table 9. Table 9: Model error (1-AUC) over a 10 10 10 -fold cross-validation for UCI datasets (averaged over 10 10 10 repetitions). Values within 10% of the minimum are displayed in bold.

Dataset	C5.0	PART	FOIL	SIRUS
Authentification	0.02	0.01	0.08	0.03
Breast Wisconsin	0.06	0.07	0.08	0.01
Credit Approval	0.15	0.17	0.15	0.09
Credit German	0.37	0.36	0.41	0.25
Diabetes	0.28	0.30	0.28	0.19
Haberman	0.46	0.42	0.50	0.35
Heart C2	0.20	0.23	0.19	0.10
Heart H2	0.23	0.23	0.23	0.12
Heart Statlog	0.21	0.24	0.20	0.10
Hepatitis	0.34	0.34	0.39	0.17
Ionosphere	0.10	0.10	0.13	0.07
Kr vs Kp	0.006	0.008	0.02	0.04
Liver Disorders	0.34	0.38	0.50	0.35
Mushrooms	0.001	0	${6.10}^{-\mathrm{𝟓}}$	${6.10}^{-4}$
Sonar	0.26	0.26	0.26	0.2
Spambase	0.07	0.07	0.12	0.07
Titanic	0.20	0.20	0.25	0.17
Vote	0.04	0.05	0.05	0.02
Wilt	0.15	0.17	0.46	0.11

Table 10. Table 10: Values of S ^ M , n , p 0 subscript ^ 𝑆 𝑀 𝑛 subscript 𝑝 0 \hat{S}_{M,n,p_{0}} averaged over p 0 ∈ V ^ M , n subscript 𝑝 0 subscript ^ 𝑉 𝑀 𝑛 p_{0}\in\hat{V}_{M,n} when the stopping criterion (5.1) is used to set M 𝑀 M , for UCI datasets. Results are averaged over 10 10 10 repetitions and standard deviations are displayed in parentheses.

Dataset	Mean stability
Haberman	0.950 (0.01)
Diabetes	0.950 (0.007)
Heart Statlog	0.954 (0.007)
Liver Disorders	0.951 (0.006)
Heart C2	0.955 (0.009)
Heart H2	0.952 (0.009)
Credit German	0.950 (0.008)
Credit Approval	0.941 (0.02)
Ionosphere	0.950 (0.009)

Table 11. Table 11: Number of trees M 𝑀 M determined by the stopping criterion (5.1) for UCI datasets. Results are averaged over 10 10 10 repetitions and standard deviations are displayed in parentheses.

Dataset	Nb of trees (sd)
Haberman	10 920 (877)
Diabetes	18 830 (1538)
Heart Statlog	7840 (994)
Liver Disorders	14 650 (1242)
Heart C2	6840 (1270)
Heart H2	4220 (529)
Credit German	7940 (672)
Credit Approval	20 650 (8460)
Ionosphere	7320 (487)

Equations

$$\mathcal{P}=\{(j_k,r_k,s_k),k=1,\dots ,d\},$$

$${\hat{H}}_n(\mathcal{P})=\left\{\text{x}\in \text{mathds}{R}^p:\{\begin{array}{ll}{\text{x}}^{(j_k)}<{\hat{q}}_{n,r_k}^{(j_k)}& \text{if }{s}_k=L\\ {\text{x}}^{(j_k)}\ge {\hat{q}}_{n,r_k}^{(j_k)}& \text{if }{s}_k=R\end{array},k=1,\dots ,d\right\}.$$

$$\forall \text{x}\in \text{mathds}{R}^p,{\hat{g}}_{n,\mathcal{P}}(\text{x})=\{\begin{array}{ll}\frac{1}{N_n({\hat{H}}_n(\mathcal{P}))}{\sum }_{i=1}^n{Y}_i\text{mathds}{1}_{{\text{X}}_i\in {\hat{H}}_n(\mathcal{P})}& \text{if}\text{x}\in {\hat{H}}_n(\mathcal{P})\\ \frac{1}{n-N_n({\hat{H}}_n(\mathcal{P}))}{\sum }_{i=1}^n{Y}_i\text{mathds}{1}_{{\text{X}}_i\notin {\hat{H}}_n(\mathcal{P})}& \text{otherwise}\end{array},$$

$${\hat{g}}_{n,\mathcal{P}}(\text{x})=\{\begin{array}{ll}0.20& \text{if}{x}^{(6)}<10.5\\ 0.50& \text{if}{x}^{(6)}\ge 10.5\end{array}.$$

$${\hat{p}}_{M,n}(\mathcal{P})=\frac{1}{M}\sum _{\ell =1}^M\text{mathds}{1}_{\mathcal{P}\in T({\Theta }_{\ell },{\mathcal{D}}_n)},$$

$$p_n(\mathcal{P})=\mathbb{P}(\mathcal{P}\in T(\Theta ,{\mathcal{D}}_n)|{\mathcal{D}}_n).$$

$${\hat{\mathcal{P}}}_{M,n,p_0}=\{\mathcal{P}\in \Pi :{\hat{p}}_{M,n}(\mathcal{P})>p_0\}.$$

$$\mathbb{E}[{\hat{p}}_{M,n}(\mathcal{P})]=\text{mathds}P(\mathcal{P}\in T(\Theta ,{\mathcal{D}}_n)),$$

$${\hat{\eta }}_{M,n,p_0}(\text{x})=\frac{1}{|{\hat{\mathcal{P}}}_{M,n,p_0}|}\sum _{\mathcal{P}\in {\hat{\mathcal{P}}}_{M,n,p_0}}{\hat{g}}_{n,\mathcal{P}}(\text{x}).$$

$$\displaystyle\hat{S}_{M,n,p_{0}}=\frac{2\big{|}\hat{\mathscr{P}}_{M,n,p_{0}}\cap\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}\big{|}}{\big{|}\hat{\mathscr{P}}_{M,n,p_{0}}\big{|}+\big{|}\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}\big{|}}$$

$$\displaystyle\begin{split}L_{n}(H,\hat{q}_{n,r}^{(j)})&\stackrel{{\scriptstyle\tiny{\rm def}}}{{=}}\frac{1}{N_{n}(H)}\sum_{i=1}^{n}(Y_{i}-\overline{Y}_{H})^{2}\mathds{1}_{\textbf{X}_{i}\in H}\\ &\quad-\frac{1}{N_{n}(H)}\sum_{i=1}^{n}\big{(}Y_{i}-\overline{Y}_{H_{L}}\mathds{1}_{X_{i}^{(j)}<\hat{q}_{n,r}^{(j)}}-\overline{Y}_{H_{R}}\mathds{1}_{X_{i}^{(j)}\geq\hat{q}_{n,r}^{(j)}}\big{)}^{2}\mathds{1}_{\textbf{X}_{i}\in H},\end{split}$$

$$H_L\stackrel{\mathrm{def}}{=}\{\text{x}\in H:{\text{x}}^{(j)}<{\hat{q}}_{n,r}^{(j)}\},H_R\stackrel{\mathrm{def}}{=}\{\text{x}\in H:{\text{x}}^{(j)}\ge {\hat{q}}_{n,r}^{(j)}\},$$

$$\displaystyle\hat{q}_{n,r}^{(j)}=\inf\big{\{}x\in\mathds{R}:\frac{1}{n}\stackrel{{\scriptstyle[}}{{i}}=1]{n}{\sum}\mathds{1}_{X_{i}^{(j)}\leq x}\geq\frac{r}{q}\big{\}}.$$

$$p_n(\mathcal{P})=\mathbb{P}(\mathcal{P}\in T(\Theta ,{\mathcal{D}}_n)|{\mathcal{D}}_n),$$

$${\hat{p}}_{M,n}(\mathcal{P})=\frac{1}{M}\sum _{\ell =1}^M\text{mathds}{1}_{\mathcal{P}\in T({\Theta }_{\ell },{\mathcal{D}}_n)}.$$

$$\underset{M\to \infty }{\lim }{\hat{p}}_{M,n}(\mathcal{P})=p_n(\mathcal{P})\text{a.s.}$$

$$q_{r}^{\star(j)}=\textrm{inf}\big{\{}x\in\mathds{R}:\mathds{P}(X^{(j)}\leq x)\geq\frac{r}{q}\big{\}},$$

$$L^{\star }(H,q_r^{\star (j)})$$

$$-\text{mathds}P(X^{(j)}<q_r^{\star (j)}|\text{X}\in H)\times \mathbb{V}[Y|X^{(j)}<q_r^{\star (j)},\text{X}\in H]$$

$$-\text{mathds}P(X^{(j)}\ge q_r^{\star (j)}|\text{X}\in H)\times \mathbb{V}[Y|X^{(j)}\ge q_r^{\star (j)},\text{X}\in H].$$

$$p^{\star }(\mathcal{P})=\text{mathds}P(\mathcal{P}\in T^{\star }(\Theta )),$$

$$\underset{n\to \infty }{\lim }{\hat{p}}_{M_n,n}(\mathcal{P})=p^{\star }(\mathcal{P})\text{in probability.}$$

$$\underset{n\to \infty }{\lim }\text{mathds}P({\hat{\mathcal{P}}}_{M_n,n,p_0}={\mathcal{P}}_{p_0}^{\star })=1,$$

$$\underset{n\to \infty }{\lim }{\hat{S}}_{M_n,n,p_0}=1\text{in probability.}$$

$$\displaystyle\mathds{P}\big{(}$$

$${\hat{S}}_{M_n,n,p_0}=\frac{2\sum _{\mathcal{P}\in \Pi }\text{mathds}{1}_{{\hat{p}}_{M_n,n}(\mathcal{P})>p_0\cap {\hat{p}}_{M_n,n}^{\prime }(\mathcal{P})>p_0}}{\sum _{\mathcal{P}\in \Pi }\text{mathds}{1}_{{\hat{p}}_{M_n,n}(\mathcal{P})>p_0}+\text{mathds}{1}_{{\hat{p}}_{M_n,n}^{\prime }(\mathcal{P})>p_0}}.$$

$$\underset{n\to \infty }{\lim }\text{mathds}{1}_{{\hat{p}}_{M_n,n}(\mathcal{P})>p_0}=\text{mathds}{1}_{p^{\star }(\mathcal{P})>p_0}\text{in probability}.$$

$${\epsilon }_{M,n,p_0}=\frac{{\sum }_{\mathcal{P}\in \Pi }\Phi (M{p}_0,M,{\hat{p}}_{M,n}(\mathcal{P}))(1-\Phi (M{p}_0,M,{\hat{p}}_{M,n}(\mathcal{P})))}{{\sum }_{\mathcal{P}\in \Pi }(1-\Phi (M{p}_0,M,{\hat{p}}_{M,n}(\mathcal{P})))},$$

$$\displaystyle\underset{M}{\textrm{argmin}}\Big{\{}\frac{1}{|\hat{V}_{M,n}|}\sum_{p_{0}\in\hat{V}_{M,n}}\varepsilon_{M,n,p_{0}}<\alpha\Big{\}}.$$

$${\hat{\eta }}_{M,n,p_0}(\text{x})=\frac{1}{|{\hat{\mathcal{P}}}_{M,n,p_0}|}\sum _{\mathcal{P}\in {\hat{\mathcal{P}}}_{M,n,p_0}}{\hat{g}}_{n,\mathcal{P}}(\text{x}).$$

$$\displaystyle\ln\Big{(}\frac{\hat{\eta}_{M,n,p_{0}}(\textbf{x})}{1-\hat{\eta}_{M,n,p_{0}}(\textbf{x})}\Big{)}=\sum_{\mathscr{P}\in\hat{\mathscr{P}}_{M,n,p_{0}}}\beta_{\mathscr{P}}\hat{g}_{n,\mathscr{P}}(\textbf{x}),$$

Full text

Keywords: classification, interpretability, rules, stability, random forests.

1 Introduction

State-of-the-art learning algorithms, typically tree ensembles or neural networks, are well-known for their remarkable predictive performance. However, this high accuracy comes at the price of complex prediction mechanisms: a large number of operations are computed for a given prediction. Because of this complexity, learning algorithms are often considered as black-boxes. This lack of interpretability is a serious limitation for many applications involving critical decisions, such as healthcare, criminal justice, or industrial process optimization. This latter example is interesting to illustrate how interpretability can be essential. Indeed, in the manufacturing industry, production processes involve complex physical and chemical phenomena, whose control and efficiency are of critical importance. In practice, data is collected along the manufacturing line, describing both the production environment and its conformity. The retrieved information enables to infer a link between the manufacturing conditions and the resulting quality at the end of the line, and then to increase the process efficiency. Since the quality of the produced entities is often characterized by a pass or fail output, the problem is in fact a classification task, and state-of-the-art learning algorithms can successfully catch patterns of these complex and nonlinear physical phenomena. However, any decision impacting the production process has long-term and heavy consequences, and therefore cannot simply rely on a blind stochastic modelling. As a matter of fact, a deep physical understanding of the forces in action is required, and this makes black-box algorithms inappropriate. In a word, models have to be interpretable, i.e., provide an understanding of the internal mechanisms that build a relation between inputs and outputs, to provide insights to guide the physical analysis. This is for example typically the case in the aeronautics industry, where the manufacturing of engine parts involves sensitive casting and forging processes. Interpretable models allow us to gain knowledge on the behavior of such production processes, which can lead, for instance, to identify or fine-tune critical parameters, improve measurement and control, optimize maintenance, or deepen understanding of physical phenomena. In the following paragraphs, we deepen the discussion about the definition of interpretability to highlight the limitations of the most popular interpretable nonlinear models: decision trees and rule algorithms (Guidotti et al., 2018). Despite their high predictivity and simple structure, these methods are unstable, which is a strong operational limitation. The goal of this article is to introduce SIRUS (Stable and Interpretable RUle Set), an interpretable rule classification algorithm which considerably improves stability over state-of-the-art methods, while preserving their simple structure, accuracy, and computational complexity.

As stated in Rüping (2006), Lipton (2016), Doshi-Velez and Kim (2017), or Murdoch et al. (2019), to date, there is no agreement in statistics and machine learning communities about a rigorous definition of interpretability. There are multiple concepts behind it, many different types of methods, and a strong dependence on the area of application and the audience. Here, we focus on models intrinsically interpretable, which directly provide insights on how inputs and outputs are related, as opposed to the post-processing of black-box models. In that case, we argue that it is possible to define minimum requirements for interpretability through the triptych “simplicity, stability, and predictivity”, in line with the framework recently proposed by Yu and Kumbier (2019). Indeed, in order to grasp how inputs and outputs are related, the structure of the model has to be simple. The notion of simplicity is implied whenever interpretability is invoked (e.g., Rüping, 2006; Freitas, 2014; Letham, 2015; Letham et al., 2015; Lipton, 2016; Ribeiro et al., 2016; Murdoch et al., 2019) and essentially refers to the model size, complexity, or the number of operations performed in the prediction mechanism. Yu (2013) defines stability as another fundamental requirement for interpretability: conclusions of a statistical analysis have to be robust to small data perturbations to be meaningful. Indeed, a specific analysis is likely to be run multiple times, eventually adding a small new batch of data, and an interpretable algorithm should be insensitive to such modifications. Otherwise, unstable models provide us with a partial and arbitrary analysis of the underlying phenomena, and arouses distrust of the domain experts. Finally, if the predictive accuracy of an interpretable model is significantly lower than the one of a state-of-the-art black-box algorithm, it clearly misses strong patterns in the data and will therefore be useless, as explained in Breiman (2001b). For example, the trivial model that outputs the empirical mean of the observations for any input is simple, stable, but brings in most cases no useful information. Thus, we add a good predictivity as an essential requirement for interpretability.

Decision trees are a class of supervised learning algorithms that recursively partition the input space and make local decisions in the cells of the resulting partition. Trees can model highly nonlinear patterns while having a simple structure, and are therefore good candidates when interpretability is required. However, trees are unstable to small data perturbations (Oates and Jensen, 1997; Guidotti and Ruggieri, 2019). More precisely, as explained in Breiman (2001b): by randomly removing only $2-3\%$ of the training data, the tree structure can be quite different, which is a strong limitation to their practical use. Another class of supervised learning methods that can model nonlinear patterns while retaining a simple structure are the so-called rule models. As such, a rule is defined as a conjunction of constraints on input variables, which form a hyperrectangle in the input space where the estimated output is constant. A collection of rules is combined to form a model. Here, the term “rule” does not stand for “classification rule” but, as is traditional in the rule learning literature, to a piecewise constant estimate that simply reads “if conditions on x, then response, else default response”. Despite their simplicity and excellent predictive skills, rule algorithms are unstable and, from this point of view, share the same limitation as decision trees (Letham et al., 2015; Murdoch et al., 2019).

In line with the above, we design SIRUS in the present paper, a new rule classification algorithm which inherits an accuracy close to random forests and the simplicity of decision trees, while having a stable structure. The core aggregation principle of random forests is kept, but instead of aggregating predictions, SIRUS focuses on the probability that a given hyperrectangle (i.e., a node) is contained in a randomized tree. The nodes with the highest probability are robust to data perturbation and represent strong patterns. They are therefore selected to form a stable rule ensemble model. Here, we provide a first illustration of SIRUS with a simple and real case: the Titanic dataset (Piech, 2016). The survival status of $887$ passengers are recorded, as well as various personal characteristics: age, sex, class, number of siblings and parents aboard, and the paid fare. SIRUS outputs the following simple set of $7$ rules, which enables to grasp at a glance the main patterns to explain passenger survival:

Average survival rate $p_{s}=39\%$.

if sex is male then $p_{s}=19\%$ else $p_{s}=74\%$

if $1^{st}$ or $2^{nd}$ class then $p_{s}=56\%$ else $p_{s}=24\%$

if

$1^{st}$ or $2^{nd}$ class

& sex is female

then $p_{s}=95\%$ else $p_{s}=25\%$

if fare $<10.5$£ then $p_{s}=20\%$ else $p_{s}=50\%$

if

no parents or

children aboard

then $p_{s}=35\%$ else $p_{s}=51\%$

if

$2^{st}$ or $3^{nd}$ class

& sex is male

then $p_{s}=14\%$ else $p_{s}=64\%$

if

sex is male

& age $\geq 15$

then $p_{s}=16\%$ else $p_{s}=72\%$

To generate the prediction for a new query point x, SIRUS checks for each rule whether the conditions are satisfied to assign one of the two possible $p_{s}$ output values. Let us say for example that $x^{(sex)}$ is female, then x satisfies the condition of the first rule, which returns $p_{s}=74\%$. Next, the $7$ rule outputs are averaged to provide the predicted probability of survival for x. The model is stable: when a $10$-fold cross-validation is run to simulate data perturbation, $5$ to $6$ rules are consistent across two folds in average. The model error (1-AUC) is $0.17$, close to the $0.13$ of random forests, whereas simplicity is drastically increased: $7$ rules versus about $10^{4}$ operations for a forest prediction.

First, we review the main rule algorithms and present their mechanism principles in Section 2. Next, Section 3 is devoted to the detailed description of SIRUS. One of the main contributions of this work is the development of a software implementation, via the R/C++ package sirus (Benard and Wright, 2020) available from CRAN, based on ranger, a high-performance random forest implementation (Wright and Ziegler, 2017). In Section 4, we show that the good empirical behavior of SIRUS is theoretically understood by proving its asymptotic stability. Then, in Section 5, we illustrate the efficiency of our algorithm through numerical experiments on real datasets. Finally, Section 6 summarizes the main contributions of the article and provides directions for future research.

2 Related Work

As stated in the introduction, SIRUS has two types of competitors: decision trees and rule algorithms. More precisely, the latter can further be split into three different kinds: classical rule algorithms based on greedy heuristics, those built on top of frequent pattern mining algorithms, and those extracted from tree ensembles.

Decision trees may be the most popular competitors of SIRUS because of their simple structure. The main algorithms are CART (Breiman et al., 1984) and C5.0 (Quinlan, 1992). However, trees are unstable as we have already highlighted. A widespread method to stabilize decision trees is bagging (Breiman, 1996), in which multiple trees are grown on perturbed data and aggregated together. Random forests is an algorithm developped by Breiman (2001a) that improves over bagging by randomizing the tree construction. Predictions are stable, accuracy is increased, but the final model is unfortunately a black box. Thus, simplicity of trees is lost, and some post-treatment mechanisms are needed to understand how random forests make their decisions. Nonetheless, even if they are useful, such treatments only provide partial information and can be difficult to operationalize for critical decisions (Rudin, 2018). For example, variable importance (Breiman, 2001a, 2003a) identifies variables that have a strong impact on the output, but not which inputs values are associated to output values of interest. Similarly, local approximation methods such as LIME (Ribeiro et al., 2016) or Tolomei et al. (2017) do not provide insights on the global relation.

Rule learning originates from the influential AQ system of Michalski (1969). Many algorithms based on greedy heuristics were subsequently developped in the 1980’s and 1990’s, including Decision List (Rivest, 1987), CN2 (Clark and Niblett, 1989), FOIL (First-Order Inductive Learner, Quinlan, 1990; Quinlan and Cameron-Jones, 1995), IREP (Incremental Reduced Error Pruning, Fürnkranz and Widmer, 1994), RIPPER (Repeated Incremental Pruning to Produce Error Reduction, Cohen, 1995), PART (Partial Decision Trees, Frank and Witten, 1998), SLIPPER (Simple Learner with Iterative Pruning to Produce Error Reduction, Cohen and Singer, 1999), LRI (Leightweight Rule Induction, Weiss and Indurkhya, 2000), and ENDER (Ensemble of Decision Rules, Dembczyński et al., 2010). Since these methods are based on greedy heuristics, they are computationally fast, but similarly to decision trees, they are unstable and their accuracy is often limited.

At the end of the 1990’s a new type of rule algorithms based on frequent pattern mining is introduced with CBA (Classification Based on Association Rules, Liu et al., 1998), then extended with CPAR (Classification based on Predictive Association Rules, Yin and Han, 2003). Frequent pattern mining is originally used to identify frequent occurrences in database mining. Since the output $Y\in\{0,1\}$ is discrete and the input data can be discretized, we can generate candidate rules for classification by identifying frequent patterns associated with each output label. This exhaustive search for association rules is computationally costly (exponential with the input dimension), and efficient heuristics are used, essentially Apriori (Agrawal et al., 1993) and Eclat (Zaki et al., 1997). The rule aggregation mechanism is specific to each algorithm. More recently, BRL (Bayesian Rule List, Letham et al., 2015) uses a more sophisticated Bayesian framework for the rule aggregation than the simple approach of CBA and CPAR, while IDS (Lakkaraju et al., 2016, Interpretable Decision Sets) uses a multi-objective optimization to select interpretable rules. Finally, CORELS (Angelino et al., 2017, Certifiably Optimal RulE ListS) generates optimal rule lists for categorical data. Interestingly, these methods exhibit quite good stability properties as we will see, higher than decision trees, but on the other hand, their predictive accuracy is worse.

The last decade has seen a resurgence of rule models through powerful algorithms based on rule extraction from tree ensembles, especially RuleFit (Friedman and Popescu, 2008) and Node harvest (Meinshausen, 2010). Notice that SIRUS is also based on this principle. More specifically, RuleFit extracts all the rules of a boosted tree ensemble (Friedman and Popescu, 2003), while Node harvest is based on random forests. Then, the extracted rules are linearly combined in a sparse linear model, respectively a logistic regression with a Lasso penalty (Tibshirani, 1996) for RuleFit, and a constraint quadratic linear program for Node harvest. These two methods have a computational complexity comparable to random forests and SIRUS, since the main step of all these algorithms is to grow a tree ensemble with a large number of trees. However, both algorithms are unstable, and both output quite complex and long lists of rules. Even running RuleFit or Node harvest multiple times on the same dataset produces quite different rule lists because of the randomness in the tree ensembles—see Appendix A.1. On the other hand, SIRUS is built to have its structure converged for the given dataset, as explained later in Section 3.

To the best of our knowledge, the signed iterative random forest method (s-iRF, Kumbier et al., 2018) is the only procedure that tackles both rule learning and stability. Using random forests, s-IRF manages to extract stable signed interactions, i.e., feature interactions enriched with a thresholding behavior for each variable, lower or higher, but without specific thresholding values. Therefore, s-IRF can be difficult to operationalize since it does not provide any specific input thresholds, and thus no precise information about the influence of input variables. On the other hand, an explicit rule model identifies specific regions of interest in the input space.

3 SIRUS Algorithm

Within the general framework of supervised (binary) classification, we assume to be given an i.i.d. sample $\mathscr{D}_{n}=\{(\textbf{X}_{i},Y_{i}),i=1,\ldots,n\}$. Each $(\textbf{X}_{i},Y_{i})$ is distributed as the generic pair $(\textbf{X},Y)$ independent of $\mathscr{D}_{n}$, where $\textbf{X}=(X^{(1)},\ldots,X^{(p)})$ is a random vector taking values in $\mathds{R}^{p}$ and $Y\in\{0,1\}$ is a binary response. Throughout the document, the distribution of $(\textbf{X},Y)$ is assumed to be unknown and is denoted by $\mathds{P}_{\textbf{X},Y}$. For $\textbf{x}\in\mathds{R}^{p}$, our goal is to accurately estimate the conditional probability $\eta(\textbf{x})=\mathds{P}(Y=1|\textbf{X}=\textbf{x})$ with few simple and stable rules.

To tackle this problem, SIRUS first builds a (slightly modified) random forest. Next, each hyperrectangle of each tree of the forest is turned into a simple decision rule, and the collection of these elementary rules is ranked based on their frequency of appearance in the forest. Finally, the most significant rules are retained and are averaged together to form an ensemble model. We describe the four steps of SIRUS algorithm in the following paragraphs: the rule generation, rule selection, rule post-treatment, and the rule aggregation. This section ends with a discussion of SIRUS stability.

Rule generation.

SIRUS uses at its core the random forest method (Breiman, 2001a), slightly modified for our purpose. As in the original procedure, each single tree in the forest is grown with a greedy heuristic that recursively partitions the input space using a random variable $\Theta$. The essential difference between our approach and Breiman’s one is that, prior to all tree constructions, the empirical $q$-quantiles of the marginal distributions over the whole dataset are computed: in each node of each tree, the best split can be selected among these empirical quantiles only. This constraint is critical to stabilize the forest structure and keeps almost intact the predictive accuracy, provided $q$ is not too small (typically of the order of 10—see the experimental Subsection 5.4). Apart from this difference, the tree growing is similar to Breiman’s original procedure. The tree randomization $\Theta$ is independent of the sample and has two independent components, denoted by $\Theta^{(S)}$ and $\Theta^{(V)}$, which are respectively used for the subsampling mechanism and randomization of the split direction. Throughout the manuscript, we let $\hat{q}_{n,r}^{(j)}$ be the empirical $r$-th $q$-quantile of $\{X_{1}^{(j)},\ldots,X_{n}^{(j)}\}$, with typically $q=10$. The construction of the individual trees is summarized in Algorithm 1 below.

The main step of SIRUS is to extract rules from the modified random forest. The cornerstone of this extraction mechanism is the notion of path in a decision tree. Indeed, a path describes the sequence of splits to go from the root of the tree to a specific (inner or terminal) node. Since a hyperrectangle is associated to each node, a rule can be defined as a piecewise constant estimate with this hyperrectangle as support. Therefore, to rigorously define the rule extraction, we introduce the symbolic representation of a path in a tree. We insist that such definition is valid for both terminal leaves and inner nodes, which are all used by SIRUS. To begin, we follow the example shown in Figure 1 with a tree of depth $2$ partitioning the input space $\mathds{R}^{2}$. For instance, let us consider the node $\mathscr{P}_{6}$ defined by the sequence of two splits $\smash{X_{i}^{(2)}\geq\hat{q}_{n,4}^{(2)}}$ and $\smash{X_{i}^{(1)}\geq\hat{q}_{n,7}^{(1)}}$. The first split is symbolized by the triplet $(2,4,R)$, whose components respectively stand for the variable index $2$, the quantile index $4$, and the right side $R$ of the split. Similarly, for the second split we cut coordinate 1 at quantile index 7, and pass to the right. Thus, the path to the considered node is defined by $\mathscr{P}_{6}=\{(2,4,R),(1,7,R)\}$. Also notice that the first split already defines the path $\mathscr{P}_{2}=\{(2,4,R)\}$, associated to the right inner node at the first level of the tree. Of course, this generalizes to each path $\mathscr{P}$ of length $d$ under the symbolic compact form

[TABLE]

where, for $k\in\{1,\ldots,d\}$, the triplet $(j_{k},r_{k},s_{k})$ describes how to move from level $(k-1)$ to level $k$, with a split using the coordinate $j_{k}\in\{1,\ldots,p\}$, the index $r_{k}\in\{1,\ldots,q-1\}$ of the corresponding quantile, and a side $s_{k}=L$ if we go the the left and $s_{k}=R$ if we go to the right.

The set of all possible such paths is denoted by $\Pi$. It is important to note that $\Pi$ is in fact a deterministic (that is, non random) quantity, which only depends upon the dimension $p$ and the order $q$ of the quantiles. Of course, given a path $\mathscr{P}\in\Pi$ one can recover the hyperrectangle (i.e., the tree node) $\hat{H}_{n}(\mathscr{P})$ associated with $\mathscr{P}$ and the entire dataset $\mathscr{D}_{n}$ via the correspondence

[TABLE]

Finally, an elementary rule $\hat{g}_{n,\mathscr{P}}$ can be defined from $\hat{H}_{n}(\mathscr{P})$ as a piecewise constant estimate: $\hat{g}_{n,\mathscr{P}}(\textbf{x})$ returns the empirical probability that the output $Y$ is of class $1$ conditional on whether the query point x belongs to $\hat{H}_{n}(\mathscr{P})$ or not. Thus, the rule $\hat{g}_{n,\mathscr{P}}$ associated to the path $\mathscr{P}\in\Pi$ is formally defined by

[TABLE]

using the convention $0/0=0$, and where $N_{n}(\hat{H}_{n}(\mathscr{P}))$ is the number of observations in the node associated with $\mathscr{P}$. This formal definition can be illustrated with the Titanic dataset presented in the introduction. For the fourth rule, fare is the $6$th variable and since $\hat{q}_{n,4}^{(6)}=10.5$, the corresponding path is $\mathscr{P}=\{(6,4,L)\}$, and the associated rule is thus

[TABLE]

Finally, a $\Theta$-random tree generates a collection of paths in $\Pi$, one for each internal and terminal nodes. In the sequel, we let $T(\Theta,\mathscr{D}_{n})$ be the list of such extracted paths, a random subset of $\Pi$.

Rule selection.

Using our modified random forest algorithm, we are able to generate a large number $M$ of trees, randomized by $\Theta_{1},\ldots,\Theta_{M}$, i.i.d. copies of the generic variable $\Theta$, and then to extract a large collection of rules. Since we are interested in selecting the most important rules, i.e., those which represent strong patterns between the inputs and the output, we select rules that are shared by a large portion of trees. Such occurrence frequency is formally defined by

[TABLE]

which is the Monte-Carlo estimate of the probability that a path $\mathscr{P}$ belongs to a $\Theta$-random tree, that is

[TABLE]

As a general strategy, once the modified random forest has been built, we draw the list of all paths that appear in the forest and only retain those that occur with a frequency larger than the threshold $p_{0}\in(0,1)$, the only influential parameter of SIRUS—see Subsection 5.4 for its tuning procedure. We are thus interested in the set of the extracted paths

[TABLE]

An important feature of SIRUS algorithm is to stop the growing of the forest with an appropriate number of trees $M$. Although the right order of magnitude for $M$ is required, no fine tuning is necessary. Indeed, the uncertainty of the importance estimate $\hat{p}_{M,n}(\mathscr{P})$ of each rule decreases with $M$, whereas the computational cost linearly increases with $M$. Thus, to obtain a robust rule extraction, $M$ needs to be high enough to make the uncertainty of $\hat{p}_{M,n}(\mathscr{P})$ negligible. More precisely, $M$ is set to get the same list of selected rules $\hat{\mathscr{P}}_{M,n,p_{0}}$ when SIRUS is run multiple times on the same dataset $\mathscr{D}_{n}$. On the other hand, $M$ should be small enough to avoid useless computations. Therefore, the growing of the forest is automatically stopped when $95\%$ of the selected rules would be shared by a new run of SIRUS on $\mathscr{D}_{n}$ in average, as it is possible to derive a simple stopping criterion based on the properties of the estimates $\hat{p}_{M,n}(\mathscr{P})$—all the technical details are provided in Subsection 5.4. A random forest is usually built with around $500$ trees, as the predictive accuracy cannot be significantly increased by adding more trees. SIRUS typically grows $10$ times more trees to obtain a robust rule extraction.

Besides, we insist that the quantile discretization is critical for the rule selection. The expected value of the rule importance is

[TABLE]

but without the discretization, the list of extracted paths from a random tree $T(\Theta,\mathscr{D}_{n})$ takes values in an uncountable space when at least one component of X is a continuous random variable, and therefore the above quantity is null, making the path selection procedure unstable with respect to data perturbation.

Rule post-treatment.

By construction, there is some redundancy in the list of rules generated by the set of distinct paths $\hat{\mathscr{P}}_{M,n,p_{0}}$. The hyperrectangles associated with the paths extracted from a $\Theta$-random tree overlap, and so the corresponding rules are linearly dependent. Therefore a post-treatment to filter $\hat{\mathscr{P}}_{M,n,p_{0}}$ is needed to remove redundancy and obtain a compact rule model. The general idea is straightforward: if the rule associated with the path $\mathscr{P}\in\hat{\mathscr{P}}_{M,n,p_{0}}$ is a linear combination of rules associated with paths with a higher frequency in the forest, then $\mathscr{P}$ is removed from $\hat{\mathscr{P}}_{M,n,p_{0}}$.

To illustrate the post-treatment, let the tree of Figure 1 be the $\Theta_{1}$-random tree grown in the forest. Since the paths of the first level of the tree, $\mathscr{P}_{1}$ and $\mathscr{P}_{2}$, always occur in the same trees, we have $\hat{p}_{M,n}(\mathscr{P}_{1})=\hat{p}_{M,n}(\mathscr{P}_{2})$. If we assume these quantities to be greater than $p_{0}$, then $\mathscr{P}_{1}$ and $\mathscr{P}_{2}$ both belong to $\hat{\mathscr{P}}_{M,n,p_{0}}$. However, by construction, $\mathscr{P}_{1}$ and $\mathscr{P}_{2}$ are associated with the same rule, and we therefore enforce SIRUS to keep only $\mathscr{P}_{1}$ in $\hat{\mathscr{P}}_{M,n,p_{0}}$. Each of the paths of the second level of the tree, $\mathscr{P}_{3}$, $\mathscr{P}_{4}$, $\mathscr{P}_{5}$, and $\mathscr{P}_{6}$, can occur in many different trees, and their associated $\hat{p}_{M,n}$ are distinct (except in very specific cases). Assume for example that $\hat{p}_{M,n}(\mathscr{P}_{1})>\hat{p}_{M,n}(\mathscr{P}_{4})>\hat{p}_{M,n}(\mathscr{P}_{5})>\hat{p}_{M,n}(\mathscr{P}_{3})>\hat{p}_{M,n}(\mathscr{P}_{6})>p_{0}$. Since $\hat{g}_{n,\mathscr{P}_{3}}$ is a linear combination of $\hat{g}_{n,\mathscr{P}_{4}}$ and $\hat{g}_{n,\mathscr{P}_{1}}$, $\mathscr{P}_{3}$ is removed. Similarly $\mathscr{P}_{6}$ is redundant with $\mathscr{P}_{1}$ and $\mathscr{P}_{5}$, and it is therefore removed. Finally, among the six paths of the tree, only $\mathscr{P}_{1}$, $\mathscr{P}_{4}$, and $\mathscr{P}_{5}$ are kept in the list $\hat{\mathscr{P}}_{M,n,p_{0}}$.

Rule aggregation.

Now, the resulting small set of rules $\hat{\mathscr{P}}_{M,n,p_{0}}$ is combined to form a simple, compact, and stable rule classification model. We simply average the set of elementary rules $\{\hat{g}_{n,\mathscr{P}}:\mathscr{P}\in\hat{\mathscr{P}}_{M,n,p_{0}}\}$ that have been selected in the first steps of SIRUS. The aggregated estimate $\hat{\eta}_{M,n,p_{0}}(\textbf{x})$ of $\eta(\textbf{x})$ is thus defined by

[TABLE]

Finally, the classification procedure assigns class $1$ to an input x if the aggregated estimate $\hat{\eta}_{M,n,p_{0}}(\textbf{x})$ is above a given threshold, and class [math] otherwise. In the introduction, we presented an example of a list of $7$ rules for the Titanic dataset. In this case, for a new input x, $\hat{\eta}_{M,n,p_{0}}(\textbf{x})$ is simply the average of the output probability of survival $p_{s}$ over the $7$ selected rules.

In past works on rule ensemble models, such as RuleFit (Friedman and Popescu, 2008) and Node harvest (Meinshausen, 2010), rules are also extracted from a tree ensemble and then combined together through a regularized linear model. In our case, it happens that the parameter $p_{0}$ alone is enough to control sparsity. Indeed, in our experiments, we observe that adding such linear model in the aggregation method hardly increases the accuracy and hardly reduces the size of the final rule set, while it can significantly reduce stability, add a set of coefficients that makes the model less straightforward to interpret, and requires more intensive computations. We refer to the experiments in Appendix A.3 for a comparison between $\hat{\eta}_{M,n,p_{0}}$ defined a as simple average (3.3) versus a definition with a logistic regression.

Categorical and numerical discrete variables.

For the sake of clarity, the description of SIRUS algorithm is limited to the case of numerical continuous variables. However, SIRUS can obviously handle numerical discrete and categorical data, as it is the case for random forests. On one hand, numerical discrete variables are left untouched since the number of possible split points is already finite, and the rule definition introduced for continuous variables also applies. On the other hand, we naturally extend the rule definition for categorical variables to “if $X^{(1)}$ is category 1 or 2 then response else default response”—see the Titanic dataset example in the introduction. Originally, categorical variables are efficiently handled in trees by transformation in ordered variables. Such ordering of categories is done with respect to the output mean for each category—see Breiman et al. (1984); Friedman et al. (2001), and we follow ranger implementation. Notice that trees are likely to often cut on categorical variables with a high number of categories, as highlighted in Strobl et al. (2006). Consequently, SIRUS is likely to output irrelevant rules associated to such categorical variables. Thus, it is best to discard categorical variables with a high number of categories, or transform them by regrouping categories or using one-hot-encoding before running SIRUS. Finally, note that ordinal variables (e.g. $X^{(1)}\in\{$small, medium, big$\}$) are treated like categorical variables.

Stability.

The three main properties to assess the interpretability of SIRUS are simplicity, stability, and predictivity, as already stated. On one hand, a measure of simplicity is naturally provided by the number of rules, and predictivity is given by the missclassification rate or the AUC. On the other hand, stability requires a more thorough discussion. In the statistical learning theory, stability refers to the stability of predictions (e.g., Vapnik, 1998). In particular, Rogers and Wagner (1978), Devroye and Wagner (1979), Bousquet and Elisseeff (2002), and Poggio et al. (2004) show that stability and predictive accuracy are closely connected. In our case, we are more concerned by the stability of the internal structure of the model, and, to our knowledge, no general definition exists. So, we state the following tentative definition: a rule learning algorithm is stable if two independent estimations based on two independent samples result in two similar lists of rules. Thus, given a new sample $\mathscr{D}_{n}^{\prime}$ independent of $\mathscr{D}_{n}$, we define $\hat{p}^{\prime}_{M,n}(\mathscr{P})$ and the corresponding set of paths $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}}$ based on a modified random forest drawn with a parameter $\Theta^{\prime}$ independent of $\Theta$. Then, we measure the stability of SIRUS by the proportion of rules shared by the two sets $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}}$ and $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}}$, selected over these two runs of SIRUS on independent samples. We take advantage of a dissimilarity measure between two sets, the so-called Dice-Sorensen index, often used to assess the stability of variable selection methods (Chao et al., 2006; Zucknick et al., 2008; Boulesteix and Slawski, 2009; He and Yu, 2010; Alelyani et al., 2011). This index is defined by

[TABLE]

with the convention $0/0=1$. This is a measure of stability taking values between [math] and $1$: if the intersection between $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}}$ and $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}}$ is empty, then $\smash{\hat{S}_{M,n,p_{0}}=0}$, while if $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}=\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}}$, then $\smash{\hat{S}_{M,n,p_{0}}=1}$. Notice that it is possible to use other metrics to assess the distance between two finite sets (Zucknick et al., 2008): the Jaccard Index is another popular example. Although the stability values slightly vary with a different definition, both the asymptotic stability of SIRUS—see Section 4—and the empirical stability comparisons between algorithms—see Section 5—are insensitive to the stability metric choice.

4 Theoretical Analysis of Stability

Among the three minimum requirements for interpretability defined in Section 1, simplicity and predictivity are quite easily met for rule models (Cohen and Singer, 1999; Meinshausen, 2010; Letham et al., 2015). On the other hand, as Letham et al. (2015) recall, building a stable rule ensemble is challenging. Therefore the main goal of this section is to prove the asymptotic stability of SIRUS, i.e., provided that the sample size is large enough, SIRUS systematically outputs the same list of rules when run multiple times with independent samples. On the other hand, we also argue that existing tree-based rule algorithms are unstable by design.

In order to show the asymptotic stability of SIRUS, we first need to introduce formal definitions of the mathematical elements involved in the empirical algorithm. We additionally define the theoretical counterpart of SIRUS, an abstract procedure which is not based on the sample $\mathscr{D}_{n}$, but only on the unknown distribution $\mathds{P}_{\textbf{X},Y}$. Next, we will prove the stochastic convergence of SIRUS towards its theoretical counterpart. This means that the list of selected rules does not depend on the training data $\mathscr{D}_{n}$, but only on $\mathds{P}_{\textbf{X},Y}$, provided that the sample size is large enough. Therefore, the same list of rules is output when SIRUS is run multiple times on independent samples. This mathematical analysis highlights that the remarkable stable behavior of SIRUS in practice has theoretical groundings, and that the discretization of the cut values with the quantiles, as well as using random forests, are the cornerstones to stabilize rule models extracted from tree ensembles.

Empirical algorithm.

First, we define the empirical CART-splitting criterion used to find the optimal split at each node of each tree of the forest. In our context of binary classification where the output $Y\in\{0,1\}$, maximizing the so-called empirical CART-splitting criterion is equivalent to maximizing the criterion based on Gini impurity (see, e.g., Biau and Scornet, 2016). More precisely, at node $H$ and for a cut performed along the $j$-th coordinate at the empirical $r$-th $q$-quantile $\hat{q}_{n,r}^{(j)}$, this criterion reads

[TABLE]

where $\overline{Y}_{H}$ is the average of the $Y_{i}$’s such that $\textbf{X}_{i}\in H$, $N_{n}(H)$ is the number of data points $\textbf{X}_{i}$ falling into $H$,

[TABLE]

and for $r\in\{1,\ldots,q-1\}$ the empirical $r$-th $q$-quantile of $\{X_{1}^{(j)},\ldots,X_{n}^{(j)}\}$ is defined by

[TABLE]

Note that, for the ease of reading, (4.1) is defined for a tree built with the entire dataset $\mathscr{D}_{n}$ without resampling. As it is often the case in the theoretical analysis of random forests, we assume throughout this section that the subsampling of $a_{n}$ observations to build each tree is done without replacement to alleviate the mathematical analysis.

Recall that the rule selection is based on the probability $p_{n}(\mathscr{P})$ that a $\Theta$-random tree of the forest contains a particular path $\mathscr{P}\in\Pi$, that is,

[TABLE]

and that the Monte-Carlo estimate $\hat{p}_{M,n}(\mathscr{P})$ of $p_{n}(\mathscr{P})$ is directly computed using the random forest, and takes the form

[TABLE]

Clearly, $\hat{p}_{M,n}(\mathscr{P})$ is a good estimate of $p_{n}(\mathscr{P})$ when $M$ is large since, by the law of large numbers, conditional on $\mathscr{D}_{n}$,

[TABLE]

We also see that $\hat{p}_{M,n}(\mathscr{P})$ is unbiased since $\mathds{E}[\hat{p}_{M,n}(\mathscr{P})|\mathscr{D}_{n}]=p_{n}(\mathscr{P}).$

Theoretical algorithm.

Next, we define all theoretical counterparts of the empirical quantities involved in SIRUS, which do not depend on $\mathscr{D}_{n}$ but only on the unknown distribution $\mathds{P}_{\textbf{X},Y}$ of $(\textbf{X},Y)$. For a given integer $q\geq 2$ and $r\in\{1,\ldots,q-1\}$, the theoretical $q$-quantiles are defined by

[TABLE]

i.e., the population version of $\hat{q}_{n,r}^{(j)}$ defined in (4.2). Similarly, for a given hyperrectangle $H\subseteq\mathds{R}^{p}$, we let the theoretical CART-splitting criterion be

[TABLE]

Based on this criterion, we denote by $T^{\star}(\Theta)$ the list of all paths contained in the theoretical tree built with randomness $\Theta$, where splits are chosen to maximize the theoretical criterion $L^{\star}$ instead of the empirical one $L_{n}$, defined in (4.1). We stress again that the list $T^{\star}(\Theta)$ does not depend upon $\mathscr{D}_{n}$ but only upon the unknown distribution of $(\textbf{X},Y)$. Next, we let $p^{\star}(\mathscr{P})$ be the theoretical counterpart of $p_{n}(\mathscr{P})$, that is

[TABLE]

and finally define the theoretical set of selected paths $\mathscr{P}^{\star}_{p_{0}}$ by $\{\mathscr{P}\in\Pi:p^{\star}(\mathscr{P})>p_{0}\}$ (with the same post-treatment as for the empirical procedure—see Section 3). Notice that, in the case where multiple splits have the same value of the theoretical CART-splitting criterion, one is randomly selected.

Consistency of the path selection.

The construction of the rule ensemble model essentially relies on the path selection and on the estimates $\hat{p}_{M,n}(\mathscr{P})$, $\mathscr{P}\in\Pi$. Therefore, our theoretical analysis first focuses on the asymptotic properties of those estimates in Theorem 1. Our consistency results hold under conditions on the subsampling rate $a_{n}$ and the number of trees $M_{n}$, together with some assumptions on the distribution of the random vector X. They are given below.

1. (A1)

The subsampling rate $a_{n}$ satisfies $\lim\limits_{n\to\infty}a_{n}=\infty$ and $\lim\limits_{n\to\infty}\frac{a_{n}}{n}=0$. 2. (A2)

The number of trees $M_{n}$ satisfies $\lim\limits_{n\to\infty}M_{n}=\infty$. 3. (A3)

X has a strictly positive density $f$ with respect to the Lebesgue measure. Furthermore, for all $j\in\{1,\ldots,p\}$, the marginal density $f^{(j)}$ of $X^{(j)}$ is continuous, bounded, and strictly positive.

We can now state the consistency of the occurrence frequency of each possible path $\mathscr{P}\in\Pi$ in the modified random forest.

Theorem 1**.**

If Assumptions (A1)-(A3) are satisfied, then, for all $\mathscr{P}\in\Pi$, we have

[TABLE]

Stability.

The only source of randomness in the selection of the rules lies in the estimates $\hat{p}_{M_{n},n}(\mathscr{P})$. Since Theorem 1 states the consistency of such an estimation, the path selection consistency follows, for all threshold values $p_{0}$ that do not belong to the finite set $\mathcal{U}^{\star}=\{p^{\star}(\mathscr{P}):\mathscr{P}\in\Pi\}$ of all theoretical probabilities of appearance for each path $\mathscr{P}$. Indeed, if $p_{0}=p^{\star}(\mathscr{P})$ for some $\mathscr{P}\in\Pi$, then $\mathds{P}(\hat{p}_{M_{n},n}(\mathscr{P})>p_{0})$ does not necessarily converge to [math] and the path selection can be inconsistent. Then, we can deduce that SIRUS is asymptotically stable in the following Corollary 1.

Corollary 1**.**

Assume that Assumptions (A1)-(A3) are satisfied. Then, provided $p_{0}\in[0,1]\setminus\mathcal{U}^{\star}$, we have

[TABLE]

and then

[TABLE]

Competitors.

As we will discuss further in the experimental Section 5, CART, C5.0, RuleFit, and Node harvest are top competitors of SIRUS, which are also based on rule extraction from trees. However, these algorithms do not include a pre-processing step of discretization, which makes them unstable by design. To see this, we first adapt the definition of an extracted path without discretization as $\mathscr{P}=\{(j_{k},z_{k},s_{k}),\,k=1,\ldots,d\}$, where $z_{k}\in\mathds{R}$ is now the cutting value of the $k$-th split. For any rule algorithm, we also define $\hat{S}_{M,n}$ as the proportion of rules shared between the output rule lists over two runs with two independent samples. Note that $M=1$ for CART and C5.0, and as already mentioned, it is possible to define a rule algorithm from CART, by extracting its nodes, as in C5.0. Thus, we obtain that for any tree-based rule algorithm, $\hat{S}_{M,n}=0$ almost surely. Indeed, since the input X takes continuous values (Assumption (A3)) and decision trees can cut at the middle of two observations in all directions, the probability that a cutting value from the tree built with $\mathscr{D}_{n}$ and one from the tree built with $\mathscr{D}_{n}^{\prime}$ are equal is null.

However, recall that in the experiments, we include a pre-processing discretization step to stabilize competitors and enable fair comparisons. With this modification, they reach a value of $\hat{S}_{M,n}>0$, but still not in par with SIRUS. This shows that the high stability improvement of SIRUS does not only come from the discretization, but mainly from the rule selection procedure, based on the probability of the rule occurrence in a random tree.

Proofs.

The proof of Theorem 1 is to be found in Appendix C. It is however interesting to give a sketch of the proof here. Corollary 1 is a direct consequence of Theorem 1, the full proof follows.

Sketch of proof of Theorem 1.

The consistency is obtained by showing that $\hat{p}_{M_{n},n}(\mathscr{P})$ is asymptotically unbiased with a null variance. The result for the variance is quite straightforward since the variance of $\hat{p}_{M_{n},n}(\mathscr{P})$ can be broken into two terms: the variance generated by the Monte-Carlo randomization, which goes to [math] as the number of trees increases (Assumption (A2)), and the variance of $p_{n}(\mathscr{P})$. Following Mentch and Hooker (2016), since $p_{n}(\mathscr{P})$ is a bagged estimate it can be seen as an infinite-order U-statistic, and a classic bound on the variance of U-statistics gives that $\mathbb{V}[p_{n}(\mathscr{P})]$ converges to [math] if $\lim_{n\to\infty}\frac{a_{n}}{n}=0$, which is true by Assumption (A1). Next, proving that $\hat{p}_{M_{n},n}(\mathscr{P})$ is asymptotically unbiased requires to dive into the internal mechanisms of the random forest algorithm. To do this, we have to show that the CART-splitting criterion is consistent (Lemma 3) and asymptotically normal (Lemma 4) when cuts are limited to empirical quantiles (estimated on the same dataset) and the number of trees grows with $n$. When cuts are performed on the theoretical quantiles, the law of large numbers and the central limit theorem can be directly applied, so that the proof of Lemmas 3 and 4 boils down to showing that the difference between the empirical CART-splitting criterion evaluated at empirical and theoretical quantiles converges to [math] in probability fast enough. This is done in Lemma 2 thanks to Assumption (A3). ∎

Proof of Corollary 1.

The first result is a consequence of Theorem 1 since

[TABLE]

Next, we have

[TABLE]

Since $p_{0}\notin\mathcal{U}^{\star}$, we deduce from Theorem 1 and the continuous mapping theorem that, for all $\mathscr{P}\in\Pi$,

[TABLE]

Therefore, $\lim\limits_{n\to\infty}\hat{S}_{M_{n},n,p_{0}}=1$ in probability. ∎

5 Experiments

We begin this section by providing overall experimental settings. Next, we focus on a case study to illustrate SIRUS with an industrial process example: the semi-conductor manufacturing process SECOM data (Dua and Graff, 2017). In particular, it shows the excellent performance of SIRUS on real data in a noisy and high-dimensional setting. In Subsection 5.3, we use $19$ UCI datasets (Dua and Graff, 2017) to perform extensive comparisons between SIRUS and its main competitors. We show that SIRUS produces much more stable rule lists, while preserving a predictive accuracy and computational complexity comparable to the top competitors. Finally, in Subsection 5.4, we detail the tuning procedure of the single hyperparameter $p_{0}$, along with a thorough discussion on the design of SIRUS. In particular, the cut limitations to the quantiles and the number of constraints in the selected rules are analyzed, and we also provide the stopping criterion for the number of trees.

5.1 Experiment Description

Performance metrics.

We first introduce relevant metrics to assess the three interpretability properties in the experiments. By definition, the size (i.e., the simplicity) of the rule ensemble is the number of selected rules, i.e., $|\hat{\mathscr{P}}_{M,n,p_{0}}|$. To measure the error, 1-AUC is used and estimated by $10$-fold cross-validation (repeated $10$ times for robustness and standard deviation estimates). With respect to stability, an independent dataset is not available for real data to compute $\hat{S}_{M,n,p_{0}}$ as defined in (3.4) in the Section 3. Nonetheless, we can take advantage of the cross-validation process to compute a stability metric: the proportion of rules shared by two models built during the cross-validation, averaged over all possible pairs (Guidotti and Ruggieri, 2019).

Datasets.

We have conducted experiments on the SECOM data, as well as $19$ diverse public datasets from the UCI repository (Dua and Graff, 2017; data is described in Table 1).

These experiments aim at illustrating the good behavior of SIRUS over its competitors in various settings. To compare stability of the different methods, data is discretized using the 10-empirical quantiles for each continuous variable and the same stability metric is used for all algorithm comparisons. For simplicity and predictivity metrics, we do not apply this pre-processing step of discretization, unless the algorithm only handles categorical data.

Competitors.

For decision trees, we run both CART and C5.0, and trees are pruned to maximize their performance. Notice that, to enable simplicity and stability comparisons for CART, a list of rules is extracted from its nodes, as it is originally possible for C5.0. For rule algorithms based on greedy heuristics, we evalute RIPPER, PART, and FOIL. Next, for rule algorithms based on tree ensembles, we evaluate RuleFit and Node harvest. Note that categorical features are transformed in multiple binary variables as it is required by the two software implementations, and RuleFit is limited to rule predictors. For RuleFit, the lasso penalty is tuned by cross-validation as defined in Friedman and Popescu (2008). As advertised in Meinshausen (2010), Node harvest does not require parameter tuning by default, but it is also possible to add a regularization term to reduce the model size. We use the same tuning procedure as for SIRUS to maximize accuracy with the smallest possible model—see Subsection 5.4. Finally, for rule algorithms based on frequent pattern mining, we run the experiments for CBA and BRL. Note that we use default settings for BRL, since modifying its parameters does not significantly improve accuracy and can hurt stability. We use available R implementations: rpart (Therneau and Atkinson, 2019, CART), C50 (Kuhn and Quinlan, 2020, C5.0), RWeka (Hornik et al., 2009, RIPPER, PART), arulesCBA (Johnson and Hahsler, 2020, FOIL, CBA), pre (Fokkema, 2020, RuleFit), nodeHarvest (Meinshausen, 2015, Node harvest), and sbrl (Yang et al., 2017, BRL). We also use our R/C++ software implementation sirus (Benard and Wright, 2020) (available from CRAN), adapted from ranger, a fast random forest implementation (Wright and Ziegler, 2017). Besides, notice that for SIRUS experiments, we use the default settings of random forests well known for their excellent behavior, in particular $\texttt{mtry}=\lfloor\frac{p}{3}\rfloor$. We set $q=10$ quantiles and tune $p_{0}$ as specified in Subsection 5.4.

5.2 Case Study: Manufacturing Process Data

SIRUS is run on a real manufacturing process of semi-conductors, the SECOM dataset (Dua and Graff, 2017). Data is collected from sensors and process measurement points to monitor the production line, resulting in $590$ numeric variables. Each of the $1567$ data points represents a single production entity associated with a pass or fail output ($0/1$) for in-house line testing. As it is often the case for a production process, the dataset is unbalanced and contains $104$ fails, i.e., a failure rate $p_{f}$ of $6.6$%. We proceed to a simple pre-processing of the data: missing values (about $5$% of the total) are replaced by the median.

Figure 2 shows predictivity versus the number of rules when $p_{0}$ varies, with the optimal $p_{0}$ displayed. Notice that the relation between $p_{0}$ and the number of rules is monotone by construction, but also highly nonlinear. Therefore, we use the number of rules for the x-axis of Figure 2 to improve readability. The 1-AUC value is $0.30$ for SIRUS (for the optimal $p_{0}=0.04$), $0.29$ for Breiman’s random forests, and $0.48$ for a pruned CART tree. Thus, in that case, CART tree predicts no better than the random classifier, whereas SIRUS has a similar accuracy to random forests. The final model has $6$ rules and a stability of $0.72$, i.e., in average $4$ to $5$ rules are shared by $2$ models built in a 10-fold cross-validation process, simulating data perturbation. By comparison, Node harvest outputs $36$ rules with a value of $0.32$ for 1-AUC.

Finally, the output of SIRUS may be displayed in the simple and interpretable form of Figure 3, the output in the R console of the package sirus for the SECOM data.

Such a rule model enables to catch immediately how the most relevant variables impact failures. Among the $590$ variables, $5$ are enough to build a model as predictive as random forests, and such a selection is robust. Other rules alone may also be informative, but they do not add additional information to the model, since predictive accuracy is already minimal with the $6$ selected rules. Then, production engineers should first focus on those $6$ rules to investigate an improved setting of the production process. We insist that the stability of the output rule list is critical in practice. Indeed, the algorithm may be run multiple times during the analysis, eventually with an additional small new batch of data. The output rule list should be quite insensitive to such perturbation: domain experts are skeptical of unstable results, which are the symptoms of a partial and arbitrary modelling of the true phenomenon. SIRUS is stable, but it is not the case for decision trees or existing rule algorithms, as we show in the next subsection and illustrate in Appendix A.1.

5.3 Improvement over Competitors

Overall, we observe that SIRUS provides a high improvement of stability compared to state-of-the-art rule algorithms, while preserving the other properties. For the top competitors, experimental results are gathered in Table 2 for model size, Table 3 for stability, and Table 4 for predictive accuracy. Experiments for additional competitors are provided in Appendix A.2 in Tables 7, 8 and 9. Standard deviations are made negligible by averaging metrics over $10$ repetitions of the cross-validation and are not displayed in the tables to increase readability.

Figure 4 provides an example for the dataset “Credit German” of the dependence between predictivity and the number of rules when $p_{0}$ varies. In that case, the minimum of 1-AUC is about $0.25$ for SIRUS, $0.20$ for Breiman’s forests, and $0.29$ for CART tree. For the chosen $p_{0}$, SIRUS returns a compact set of $22$ rules and its stability is $0.66$. Figure 5 provides another example of the good practical performance of SIRUS with the “Heart Statlog” dataset. Here, the predictivity of random forests is reached with $16$ rules, with a stability of $0.83$, i.e., about $13$ rules are consistent between two different models built in a $10$-fold cross-validation. Thus, the final models are simple, quite robust to data perturbation, and have a predictive accuracy close to random forests.

We can draw the following conclusions from the experimental comparisons with competitors, displayed in Tables 2, 3, and 4. SIRUS produces more stable and predictive rule lists than decision trees, for a comparable simplicity, but at the price of a higher computational complexity since many trees are grown. SIRUS produces much more stable and shorter rule lists than RuleFit and Node harvest, for a comparable accuracy and computational complexity. Classical rule algorithms exhibit similar properties as decision trees: a smaller computational complexity, but a high instability and a reduced predictivity. Finally, algorithms based on frequent pattern mining exhibit quite good stability properties, higher than for the other types of competitors. On the other hand, their predictive accuracy is worse than decision trees. Experiments in Tables 2, 3, and 4 show that SIRUS exhibits a high stability and predictivity improvement over these methods. Besides, simplicity varies across algorithms: CBA produces much longer rule lists than SIRUS, whereas BRL generates shorter models.

5.4 SIRUS Parameters

SIRUS relies on a single tuning hyperparameter: the selection threshold $p_{0}$ involved in the definition of $\hat{\mathscr{P}}_{M,n,p_{0}}$ to filter the most important rules, which therefore controls the simplicity of the model, and consequently also its accuracy and stability. On the other hand, SIRUS is not very sensitive to the other parameters: the number of trees, the number of quantiles, and the tree depth. Therefore, they do not require fine tuning, and we simply set efficient default values as explained below.

Tuning of SIRUS.

This parameter $p_{0}$ should be set to optimize a tradeoff between the number of rules, stability, and accuracy. In practice, it is difficult to settle such a criterion, and we choose to optimize $p_{0}$ to maximize the predictive accuracy with the smallest possible set of rules. To achieve this goal, we proceed as follows. The error 1-AUC is estimated by $10$-fold cross-validation for a fine grid of $p_{0}$ values, defined such that $|\hat{\mathscr{P}}_{M,n,p_{0}}|$ varies from $1$ to $25$ rules. (We let $25$ be an arbitrary upper bound on the maximum number of rules, considering that a bigger set is not readable anymore.) The randomization introduced by the partition of the dataset in the $10$ folds of the cross-validation process has a significant impact on the variability of the size of the final model. Therefore, in order to get a robust estimation of $p_{0}$, the cross-validation is repeated multiple times (typically $10$) and results are averaged. The standard deviation of the mean of 1-AUC is computed over these repetitions for each $p_{0}$ of the grid search. We consider that all models within $2$ standard deviations of the minimum of 1-AUC are not significantly less predictive than the optimal one. Thus, among these models, the one with the smallest number of rules is selected, i.e., the optimal $p_{0}$ is shifted towards higher values to reduce the model size without decreasing predictivity—see Figures 4 and 5 for examples. This approach is very similar to the tuning procedure of the Lasso (Tibshirani, 1996).

Number of trees.

The accuracy, stability, and computational cost of SIRUS increase with the number of trees $M$. Thus, we simply design a stopping criterion to grow the minimum number of trees which ensures that accuracy and stability are higher than $95\%$ of their maximum asymptotic values with respect to $M$ and conditionally on $\mathscr{D}_{n}$. We empirically observe that the stability requirement is met for a much higher number of trees than the accuracy requirement (about $10$ times). Therefore, the stopping criterion is only based on stability. More precisely, we require that $95\%$ of the rules are identical across two runs of SIRUS on a given dataset $\mathscr{D}_{n}$ in average. Formally, the mean stability $\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]$ measures the expected proportion of rules shared by two fits of SIRUS on $\mathscr{D}_{n}$, for fixed $n$ (sample size), $p_{0}$ (threshold), and $M$ (number of trees). Thus, the stopping criterion takes the form $1-\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]<\alpha$, with typically $\alpha=0.05$.

There are two obstacles to operationalize this stopping criterion: its estimation and its dependence to $p_{0}$. We make two approximations to overcome these limitations and give empirical and theoretical evidence of their good practical behavior in Appendix B. First, Theorem 2 in Appendix B.2 provides an asymptotic equivalent with respect to $M$ of $1-\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]$, that we simply estimate by

[TABLE]

where $\Phi(Mp_{0},M,p_{n}(\mathscr{P}))$ is the cdf of a binomial distribution with parameter $p_{n}(\mathscr{P})$, $M$ trials, evaluated at $Mp_{0}$. Secondly, $\varepsilon_{M,n,p_{0}}$ depends on $p_{0}$, whose optimal value is unknown in the first step of SIRUS, when trees are grown. It turns out however that $\varepsilon_{M,n,p_{0}}$ is not very sensitive to $p_{0}$, as shown by the experiments in Appendix B.1. Consequently, our strategy is to simply average $\varepsilon_{M,n,p_{0}}$ over a set $\hat{V}_{M,n}$ of many possible values of $p_{0}$ and use the resulting average as a gauge. These values are chosen to scan all possible path sets $\hat{\mathscr{P}}_{M,n,p_{0}}$, of size ranging from $1$ to $50$ paths. When a set of $50$ paths is post-treated, its size reduces to around $25$ paths (as explained in the previous paragraph, $25$ is an arbitrarily threshold on the maximum number of rules above which a rule model is not readable anymore). In order to generate path sets of such sizes, values of $p_{0}$ are chosen halfway between two distinct consecutive $\hat{p}_{M,n}(\mathscr{P}),\mathscr{P}\in\Pi$, restricted to the highest $50$ values. Thus, in the experiments, we utilize the following criterion to stop the growing of the forest, with typically $\alpha=0.05$:

[TABLE]

Quantile discretization.

In the modified random forest grown in the first step of SIRUS, the split at each tree node is limited to the empirical $q$-quantiles of each component of X, as described in Section 3. Thus, we check that this modification alone of the forest has little impact on its accuracy. Using the R package ranger, 1-AUC is estimated for each dataset with 10-fold cross-validation for $q\in\{2,5,10,20\}$. We leave aside datasets with a majority of categorical variables, results are averaged over 10 repetitions of the cross-validation, and displayed in Table 5. Clearly, the decrease of accuracy generated by this discretization is small, and not very sensitive to $q$, provided that $q$ is not too small. Thus, $q=10$ appears to be a good default choice from the experiments. In fact, the small impact of the discretization on the forest error is not surprising: with only $p=10$ input variables, the input space is split in a fine grid of $10^{10}$ hyperrectangles for $q=10$ quantiles, providing a high flexibility to the modified random forest to identify local patterns.

Tree depth.

When SIRUS is fit using fully grown trees, the final set of rules $\hat{\mathscr{P}}_{M,n,p_{0}}$ contains almost exclusively rules made of one or two splits, and rarely of three splits. Although this may appear surprising at first glance, this phenomenon is in fact expected. Indeed, rules made of multiple splits are extracted from deeper tree levels and are thus more sensitive to data perturbation by construction. This results in much smaller values of $\hat{p}_{M,n}(\mathscr{P})$ for rules with a high number of splits, and then deletion from the final set of path through the threshold $p_{0}$: $\smash{\hat{\mathscr{P}}_{M,n,p_{0}}=\{\mathscr{P}\in\Pi:\hat{p}_{M,n}(\mathscr{P})>p_{0}\}}$. To illustrate this, let us consider the following typical example with $p=100$ input variables and $q=10$ quantiles. There are $qp=100\times 10=10^{3}$ possible splits at the root node of a tree, and then $2pq=2.10^{3}$ paths of one split. Since the left and right paths of one split at the root node are associated to the same rule, there are $qp=10^{3}$ distinct rules of one split, about $(2qp)^{2}\approx 10^{6}$ distinct rules of two splits, and about $(2qp)^{3}\approx 10^{10}$ distinct rules of three splits. Using only rules of one split is too restrictive since it generates a small model class (a thousand rules for $100$ input variables) and does not handle variable interactions. On the other hand, rules of two splits are numerous (about one million) and thus provide a large flexibility to SIRUS. More importantly, since there are $10$ billion rules of three splits, a stable selection of a few of them is clearly a difficult task, and such complex rules are naturally discarded by SIRUS.

In the software implementation sirus, the tree depth parameter max.depth is a modifiable input, set to $2$ by default to reduce the computational cost while leaving the output list of rules almost untouched as explained above. We conduct experiments where SIRUS is run with a tree depth of $1$, $2$, and $3$, and results are displayed in Table 6. Over the nineteen UCI datasets, rules of three splits appear in SIRUS rule list in only four cases, and a significant accuracy improvement over a tree depth of $2$ occurs only once, for the ‘Mushrooms’ dataset. On the other hand, for all datasets except two, SIRUS outputs rules of two constraints, and predictivity is improved over a tree depth of $1$ for half of the datasets. The Titanic example shows how the rule list is drastically simplified by limiting tree depth to $1$, lowering the insights provided by SIRUS:

Average survival rate $p_{s}=39\%$.

if sex is male then $p_{s}=19\%$ else $p_{s}=74\%$

if $1^{st}$ or $2^{nd}$ class then $p_{s}=56\%$ else $p_{s}=24\%$

This analysis of tree depth is not new. Indeed, both RuleFit (Friedman and Popescu, 2008) and Node harvest (Meinshausen, 2010) articles discuss the optimal tree depth for the rule extraction from a tree ensemble in their experiments. They both conclude that the optimal depth is $2$. Hence, the same hard limit of $2$ is used in Node harvest. RuleFit is slightly less restrictive: for each tree, its depth is randomly sampled with an exponential distribution concentrated on $2$, but allowing few trees of depth $1$, $3$, and $4$. We insist that they both reach such conclusion without considering stability issues, but only focusing on accuracy. Further considering stability properties consolidates that growing shallow trees is optimal for rule extraction from tree ensembles.

6 Conclusion

Interpretability of learning algorithms is required for applications involving critical decisions, for example the analysis of production processes in the manufacturing industry. Although interpretability does not have a precise definition, we argued that simplicity, stability, and predictivity are minimum requirements. In particular, decision trees and rule algorithms both combine a simple structure and a good accuracy for nonlinear data, and are thus considered as state-of-the-art interpretable algorithms. However, these methods are unstable with respect to data perturbation, which is a strong operational limitation. Therefore, we proposed a new rule algorithm for classification, SIRUS (Stable and Interpretable RUle Set), which takes the form of a short list of rules. We proved that SIRUS considerably improves stability over state-of-the-art algorithms, while preserving simplicity, accuracy, and computational complexity of top competitors. The principle of SIRUS is to extract rules from a random forest, based on their probability of occurrence in a random tree, and to stop the growing of the forest when the rule selection is converged. Thus, SIRUS inherits the computational complexity of random forests, and has only one tuning parameter. A software implementation, the R/C++ package sirus (Benard and Wright, 2020), is available from CRAN. Besides, we believe that the extension of SIRUS to regression is a promising future research direction: the main challenge is the construction of an appropriate rule aggregation framework to accurately estimate continuous outputs without hurting stability. Furthermore, although SIRUS has the ability to handle high-dimensional data, as illustrated with the SECOM dataset (590 inputs), specific variable selection strategies could be used to reduce the number of possible rules and then improve SIRUS performance.

Appendix A Additional Experiments

A.1 Robustness Illustration

For the SECOM dataset used in the experimental Section 5 of the article, only three rule algorithms achieve the same predictivity as random forests: RuleFit, Node harvest, and SIRUS (1-AUC of $0.30$, whereas CART and BRL are no better than the random classifier with an error of 1-AUC $=0.5$). SIRUS produces a short and stable list of $6$ rules, while RuleFit and Node harvest generate complex, long, and unstable rule lists. Rule algorithms based on tree ensembles are stochastic since they rely on the tree randomness $\Theta_{1},\ldots,\Theta_{M}$. Consequently, RuleFit and Node harvest output different rule lists when run multiple times on the same dataset. Such behavior is a strong limitation in practice, as domain experts become skeptical of the algorithm conclusions. On the other hand, SIRUS is built to have a robust rule extraction mechanism, and the same list of rules is output over multiple repetitions with the same data, as proved in Theorem 2 in the next Section.

To illustrate this, we run each algorithm twice on the SECOM dataset, and display the output models in Figure 6 for SIRUS, Figure 7 for Node harvest, and Figure 8 for RuleFit. We set the regularization parameter of Node harvest and SIRUS as explained in Subsection $5.3$ of the article, to maximize accuracy with the smallest possible model: for Node harvest $\lambda=4$, and for SIRUS $p_{0}=0.04$. RuleFit is tuned as defined in Friedman and Popescu (2008). Figures 7 and 8 show that the rule lists output by RuleFit and Node harvest are quite different across multiple runs with the exact same data, while SIRUS has the same output.

We also observe that for the same accuracy, RuleFit and Node harvest models are longer and more complex than SIRUS. In addition, rules are aggregated using weights to generate predictions. This is not the case for SIRUS, which simply averages the $6$ output rules. Finally, we can also mention that manually increasing the regularization of Node harvest, to reduce the model size to $6$ rules as in SIRUS, strongly hurts accuracy, which drops to $0.39$.

A.2 Additional Competitors

Additional experiments are provided to compare SIRUS to other competitors: C5.0 (Quinlan, 1992) (decision tree), PART (Frank and Witten, 1998), and FOIL (Quinlan and Cameron-Jones, 1995) (classical rule learning algorithms). Model size results are provided in Table 7, stability in Table 8, and error in Table 9. The stability and accuracy improvement of SIRUS is clear.

A.3 Rule Aggregation

In Section $3$ of the article, $\hat{\eta}_{M,n,p_{0}}(\textbf{x})$ ($3.3$) is a simple average of the set of rules, defined as

[TABLE]

To tackle our binary classification problem, a natural approach would be to use a logistic regression and define

[TABLE]

where the coefficients $\beta_{\mathscr{P}}$ have to be estimated. To illustrate the performance of the logistic regression (A.2), we consider again the UCI dataset, “Credit German”. We augment the previous results from Figure $4$ (in Section $5$ of the article) with the logistic regression error in Figure 9.

One can observe that the predictive accuracy is slightly improved but it comes at the price of an additional set of coefficients that can be hard to interpret (some can be negative), and an increased computational cost. Notice that categorical variables are one-hot-encoded in this example.

Appendix B Stopping Criterion for the Number of Trees $M$

We recall that the definition of the stopping criterion ($5.1$) of the forest growing is provided in Section $5$ of the main article. First, we provide three groups of experiments to show its good empirical efficiency. In the second subsection, we provide theoretical properties of the stopping criterion.

B.1 Experiments

The following experiments on the UCI datasets show the good empirical performance of the stopping criterion ($5.1$). Recall that the goal of this criterion is to determine the minimum number of trees $M$ ensuring that two independent fits of SIRUS on the same dataset result in two lists of rules with an overlap of $95$% in average. This is checked with a first batch of experiments—see next paragraph. Secondly, the stopping criterion (5.1) does not consider the optimal $p_{0}$, unknown when trees are grown in the first step of SIRUS. Then, another batch of experiments is run to show that the stability approximation $1-\varepsilon_{M,n,p_{0}}$ is quite insensitive to $p_{0}$. Finally, a last batch of experiments provides examples of the number of trees grown when SIRUS is fit. Notice that for these experiments, categorical variables are one-hot-encoded.

Experiments 1.

For each dataset, the following procedure is applied. SIRUS is run a first time using criterion (5.1) to stop the number of trees. This initial run provides the optimal number of trees $M$ as well as the set $\hat{V}_{M,n}$ of possible $p_{0}$. Then, SIRUS is fit twice independently using the precomputed number of trees $M$. For each $p_{0}\in\hat{V}_{M,n}$, the stability metric $\hat{S}_{M,n,p_{0}}$ (with $\mathscr{D}_{n}^{\prime}=\mathscr{D}_{n}$) is computed over the two resulting lists of rules. Finally $\hat{S}_{M,n,p_{0}}$ is averaged across all $p_{0}$ values in $\hat{V}_{M,n}$. This procedure is repeated $10$ times: results are averaged and presented in Table 10, with standard deviations in parentheses. Across the considered datasets, resulting values range from $0.941$ to $0.955$, and are thus close to $0.95$ as expected by construction of criterion (5.1).

Experiments 2.

The second type of experiments illustrates that $\varepsilon_{M,n,p_{0}}$ is quite insensitive to $p_{0}$ when $M$ is set with criterion (5.1). For the “Credit German” dataset, we fit SIRUS and then compute $1-\varepsilon_{M,n,p_{0}}$ for each $p_{0}\in\hat{V}_{M,n}$. Results are displayed in Figure 10. $1-\varepsilon_{M,n,p_{0}}$ ranges from $0.90$ to $1$, where the extreme values are reached for $p_{0}$ corresponding to very small number of rules, which are not of interest when $p_{0}$ is selected to maximize predictive accuracy. Thus, $1-\varepsilon_{M,n,p_{0}}$ is quite concentrated around $0.95$ when $p_{0}$ varies.

Experiments 3.

Finally, we display in Table 11 the optimal number of trees when the growing of SIRUS is stopped using criterion (5.1). It ranges from $4220$ to $20\,650$ trees. In Breiman’s forests, the number of trees above which the accuracy cannot be significantly improved is typically $10$ times lower. However SIRUS grows shallow trees, and is thus not computationally more demanding than random forests overall.

B.2 Theoretical Properties

We emphasize that growing more trees does not improve predictive accuracy or stability with respect to data perturbation for a fixed sample size $n$. Indeed, the instability of the rule selection is generated by the variance of the estimates $\hat{p}_{M,n}(\mathscr{P}),\mathscr{P}\in\Pi$. Upon noting that we have two sources of randomness—$\Theta$ and $\mathscr{D}_{n}$—, the law of total variance shows that $\mathbb{V}[\hat{p}_{M,n}(\mathscr{P})]$ can be broken down into two terms: the variance generated by the Monte Carlo randomness $\Theta$ on the one hand, and the sampling variance on the other hand. In fact, equation (C.3) in the proof of Theorem 1 below reveals that

[TABLE]

The stopping criterion (5.1) ensures that the first term becomes negligible as $M\to\infty$, so that $\mathbb{V}[\hat{p}_{M,n}(\mathscr{P})]$ reduces to the sampling variance $\mathbb{V}[p_{n}(\mathscr{P})]$, which is independent of $M$. Therefore, stability with respect to data perturbation cannot be further improved by increasing the number of trees. Additionally, the trees are only involved in the selection of the paths. For a given set of paths $\hat{\mathscr{P}}_{M,n,p_{0}}$, the construction of the final aggregated estimate $\hat{\eta}_{M,n,p_{0}}$ (see Section $3$ of the article) is independent of the forest. Thus, if further increasing the number of trees does not impact the path selection, neither it improves the predictive accuracy.

Next, Theorem 2 states that conditionally on $\mathscr{D}_{n}$ and with $\mathscr{D}_{n}^{\prime}=\mathscr{D}_{n}$, $\hat{S}_{M,n,p_{0}}$ should be close to $1$, and also provides an asymptotic approximation of $\smash{\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]}$ for large values of the number of trees $M$, which quantifies the influence of $M$ on the mean stability, conditional on $\mathscr{D}_{n}$. We let $\smash{\mathcal{U}_{n}\stackrel{{\scriptstyle\tiny{\rm def}}}{{=}}\{p_{n}(\mathscr{P}):\mathscr{P}\in\Pi\}}$ be the empirical counterpart of $\mathcal{U}^{\star}$.

Theorem 2**.**

If $p_{0}\in[0,1]\setminus\mathcal{U}_{n}$ and $\mathscr{D}^{\prime}_{n}=\mathscr{D}_{n}$, then, conditional on $\mathscr{D}_{n}$, we have

[TABLE]

In addition, for all $p_{0}<\max~{}\mathcal{U}_{n}$,

[TABLE]

where $\Phi(Mp_{0},M,p_{n}(\mathscr{P}))$ is the cdf of a binomial distribution with parameter $p_{n}(\mathscr{P})$, $M$ trials, evaluated at $Mp_{0}$, and, for all $\mathscr{P},\mathscr{P}^{\prime}\in\Pi$,

[TABLE]

and

[TABLE]

The proof of Theorem 2 is to be found in Section D. The equivalent provided in Theorem 2 is defined when the sets of rules $\hat{\mathscr{P}}_{M,n,p_{0}}$ and $\hat{\mathscr{P}}_{M,n,p_{0}}^{\prime}$ are not post-treated. It considerably simplifies the analysis of the asymptotic behavior of $\smash{\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]}$. Since the post-treatment is deterministic, this operation is not an additional source of instability. Then, if the estimation of the rule set without post-treatment is stable, it is also the case when the post-treatment is added. Finally, despite its apparent complexity, the asymptotic approximation of $1-\mathbb{E}[\hat{S}_{M,n,p_{0}}|\mathscr{D}_{n}]$ can be easily estimated, and an efficient stopping criterion for the number of trees is therefore derived in (5.1).

Appendix C Proof of Theorem 1

We recall Assumptions (A1)-(A3) and Theorem 1 for the sake of clarity.

1. (A1)

The subsampling rate $a_{n}$ satisfies $\lim\limits_{n\to\infty}a_{n}=\infty$ and $\lim\limits_{n\to\infty}\frac{a_{n}}{n}=0$. 2. (A2)

The number of trees $M_{n}$ satisfies $\lim\limits_{n\to\infty}M_{n}=\infty$. 3. (A3)

X has a strictly positive density $f$ with respect to the Lebesgue measure. Furthermore, for all $j\in\{1,\ldots,p\}$, the marginal density $f^{(j)}$ of $X^{(j)}$ is continuous, bounded, and strictly positive.

Theorem 1**.**

If Assumptions (A1)-(A3) are satisfied, then, for all $\mathscr{P}\in\Pi$, we have

[TABLE]

First, we prove Theorem 1 for a path of one split. The proof is extended for a path of two splits in the next subsection and follows the same steps. Finally, the proof can be easily extended to a path of any depth $d\in\mathbb{N}^{\star}$ by recursion.

C.1 Proof of Theorem 1 for a path of one split

We consider $\mathscr{P}_{1}=\{(j_{1},r_{1},s_{1})\}$ a path of one split, where $j_{1}\in\{1,\ldots,p\}$, $r_{1}\in\{1,\ldots,q-1\}$, and $s_{1}\in\{L,R\}$. We assume throughout that Assumptions (A1)-(A3) are satisfied.

Before proving Theorem 1, we state five lemmas (Lemma 1 to Lemma 5). Their proof can be found in the Subsection C.3. Lemma 1 is a preliminary technical result used to state both Lemmas 2 and 4 - case (b).

Lemma 1**.**

Let X be a random variable distributed on $\mathds{R}^{p}$ such that Assumptions (A1) and (A3) are satisfied. Then, for all $j\in\{1,\ldots,p\}$ and all $r\in\{1,\ldots,q-1\}$, we have

[TABLE]

and

[TABLE]

Lemma 2 is used to prove both consistency (Lemma 3) and convergence rate (Lemma 4) of the CART-splitting criterion when the root node of the tree is cut at an empirical quantile. Lemma 5 is an intermediate result to prove Theorem 1.

Lemma 2**.**

If Assumptions (A1) and (A3) are satisfied, then for all $j\in\{1,\ldots,p\}$, all $r\in\{1,\ldots,q-1\}$, and all $H\subseteq\mathds{R}^{p}$ such that $\mathbb{P}(\textbf{X}\in H,X^{(j)}<q^{\star(j)}_{r})>0$ and $\mathbb{P}(\textbf{X}\in H,X^{(j)}\geq q^{\star(j)}_{r})>0$, we have

[TABLE]

Lemma 3**.**

If Assumptions (A1) and (A3) are satisfied, then for all $j\in\{1,\ldots,p\}$, all $r\in\{1,\ldots,q-1\}$, and all $H\subseteq\mathds{R}^{p}$ such that $\mathbb{P}(\textbf{X}\in H,X^{(j)}<q^{\star(j)}_{r})>0$ and $\mathbb{P}(\textbf{X}\in H,X^{(j)}\geq q^{\star(j)}_{r})>0$, we have

[TABLE]

When splitting a node, if the theoretical CART-splitting criterion has multiple maxima, one is randomly selected. This random selection follows a discrete probability law, which is not necessarily uniform and is based on $\mathds{P}_{\textbf{X},Y}$ as specified in Definition 1. In order to derive the limit of the probability that a given split occurs in a $\Theta$-random tree in the empirical algorithm, one needs to assess the convergence rate of the empirical CART-splitting criterion when it has multiple maxima.

Lemma 4**.**

Consider that Assumptions (A1) and (A3) are satisfied. Let $\mathcal{C}_{1}\subset\{1,\ldots,p\}\times\{1,\ldots,q-1\}$ be a set of splits of cardinality $c_{1}\geq 2$, such that, for all $(j,r)\in\mathcal{C}_{1}$, $L^{\star}(\mathds{R}^{p},q_{r}^{\star(j)})\stackrel{{\scriptstyle\tiny{\rm def}}}{{=}}L^{\star}_{\mathcal{C}_{1}}$, i.e., the theoretical CART-splitting criterion is constant for all splits in $\mathcal{C}_{1}$. Let $(j_{1},r_{1})\in\mathcal{C}_{1}$ and let $\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1})}$ be a random vector where each component is the difference between the empirical CART-splitting criterion for the splits $(j,r)\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$ and $(j_{1},r_{1})$, that is

[TABLE]

(a)

If $L^{\star}_{\mathcal{C}_{1}}>0$, then we have

[TABLE]

where, for all $(j,r),(j^{\prime},r^{\prime})\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$, each element of the covariance matrix $\Sigma$ is defined by $\Sigma_{(j,r),(j^{\prime},r^{\prime})}=\emph{Cov}[Z_{j,r},Z_{j^{\prime},r^{\prime}}],$ with

[TABLE]

Besides, for all $(j,r)\in\mathcal{C}_{1}$, $\mathds{V}[Z_{j,r}]>0$.

(b)

If $L^{\star}_{\mathcal{C}_{1}}=0$, then we have

[TABLE]

where $\mathbf{V}$ is a Gaussian vector of covariance matrix $\emph{Cov}[\mathbf{Z}]$. If $\mathcal{C}_{1}$ is explicitly written $\mathcal{C}_{1}=\{(j_{k},r_{k})\}_{k=1,\ldots,c_{1}}$, $\mathbf{Z}$ is defined, for $k\in\{1,\ldots,c_{1}\}$, by

[TABLE]

where $p_{L,k}=\mathds{P}(X^{(j_{k})}<q_{r_{k}}^{\star(j_{k})})$, $p_{R,k}=\mathds{P}(X^{(j_{k})}\geq q_{r_{k}}^{\star(j_{k})})$, and $h_{\mathscr{P}_{1}}$ is a multivariate quadratic form defined as

[TABLE]

*Besides, the variance of each component of $h_{\mathscr{P}_{1}}(\mathbf{V})$ is strictly positive. *

Definition 1** (Theoretical splitting procedure).**

Let $\theta_{1}^{(V)}$ be the set of eligible variables to split the root node of a theoretical random tree. The set of best theoretical cuts at the root node is defined as

[TABLE]

If $\mathcal{C}_{1}^{\star}(\theta_{1}^{(V)})$ has multiple elements, then $(j_{1},r_{1})$ is randomly drawn with probability

[TABLE]

where $\Phi_{\theta_{1}^{(V)},(j_{1},r_{1})}$ is the cdf of the limit law defined in Lemma 4 for $\mathcal{C}_{1}=\mathcal{C}_{1}^{\star}(\theta_{1}^{(V)})$. This definition is extended for the second split in Definition 2.

Recall that the randomness in a tree can be decomposed as $\Theta=(\Theta^{(S)},\Theta^{(V)})$, where $\Theta^{(S)}$ corresponds to the subsampling and $\Theta^{(V)}$ is related to the variable selection. $\Theta^{(V)}$ takes values in the finite set $\Omega^{(V)}=\{1,\ldots,p\}^{3\times\texttt{mtry}}$.

Lemma 5**.**

If Assumptions (A1)-(A3) are satisfied, then for all $\theta^{(V)}\in\Omega^{(V)}$, we have

[TABLE]

We are now equipped to prove Theorem 1 in the case of one single split. Recall that

[TABLE]

Since $\Theta^{(V)}$ takes values in the finite set $\Omega^{(V)}$, according to Lemma 5, we have

[TABLE]

Therefore,

[TABLE]

To finish the proof, we just have to show that $\lim\limits_{n\to\infty}\mathds{V}[\hat{p}_{M_{n},n}(\mathscr{P}_{1})]=0$.

The law of total variance gives

[TABLE]

Following the approach of Mentch and Hooker (2016), $p_{n}(\mathscr{P}_{1})$ is a complete infinite order U-statistic with the kernel $\mathbb{E}[\mathds{1}_{\mathscr{P}_{1}\in T(\Theta,\mathscr{D}_{n})}|\Theta^{(S)},\mathscr{D}_{n}]$. From Hoeffding (1948),

[TABLE]

where $\xi_{a_{n},a_{n}}=\mathbb{V}[\mathbb{E}[\mathds{1}_{\mathscr{P}_{1}\in T(\Theta,\mathscr{D}_{n})}|\Theta^{(S)},\mathscr{D}_{n}]|\Theta^{(S)}]$. Since $\xi_{a_{n},a_{n}}$ is bounded and $\lim\limits_{n\to\infty}\frac{a_{n}}{n}=0$,

[TABLE]

Using equality (C.3), since $p_{n}(\mathscr{P}_{1})$ is bounded and $\lim\limits_{n\to\infty}M_{n}=\infty$,

[TABLE]

Finally,

[TABLE]

which concludes the proof.

C.2 Proof of Theorem 1 for a path of two split

The proof of Theorem 1 is extended for a path of two splits. We consider $\mathscr{P}_{1}=\{(j_{1},r_{1},s_{1})\}$ a path of one split and $\mathscr{P}_{2}=\{(j_{k},r_{k},s_{k}),k=1,2\}$ a path of two splits, where $j_{1},j_{2}\in\{1,...,p\}$, $r_{1},r_{2}\in\{1,...,q-1\}$ and $s_{1},s_{2}\in\{L,R\}$. We assume assumptions (A1)-(A3) are satisfied.

The path $\mathscr{P}_{2}=\{(j_{1},r_{1},s_{1}),(j_{2},r_{2},s_{2})\}$ can occur in trees where the split at the root node is $(j_{1},r_{1})$ and the split of one of the child node is $(j_{2},r_{2})$, and in trees where the splits are made in the reversed order, $(j_{2},r_{2})$ at the root node and $(j_{1},r_{1})$ at one of the child node. Since these two events are disjoint, \mathds{P}\big{(}\mathscr{P}_{2}\in T(\Theta,\mathscr{D}_{n})\big{)} is the sum of the probability of these two events. Without loss of generality, we will consider in the entire proof that the split at the root node is $(j_{1},r_{1})$. Lemmas 6 - 9 below extend Lemmas 2 - 5 to the case where the tree path contains two splits.

We need to introduce additional notations, first, the theoretical hyperrectangle based on a path $\mathscr{P}$ by

[TABLE]

with $d\in\{1,2\}$, the empirical counterpart of $\hat{H}_{n}(\mathscr{P})$ defined in (2.3). Furthermore, since from assumption (A3), X has a strictly positive density, then for $j\in\{1,...,p\}\setminus j_{1}$, and $r\in\{1,...,q-1\}$, \mathbb{P}\big{(}\textbf{X}\in H^{\star}(\mathscr{P}_{1}),X^{(j)}<q^{\star(j)}_{r}\big{)}>0 and \mathbb{P}\big{(}\textbf{X}\in H^{\star}(\mathscr{P}_{1}),X^{(j)}\geq q^{\star(j)}_{r}\big{)}>0. When $j=j_{1}$, the second cut is performed along the same direction as the first one. In that case, depending on the side $s_{1}$ of the first cut and the cut positions $r_{1}$ and $r$, one of the two child node can be empty with probability one. For example, the hyperrectangle associated to the path $\{(1,2,L),(1,3,R)\}$ is empty. In SIRUS, such splits are not considered to find the best cut for a node at the second level of the tree. Thus we define $\mathcal{C}_{\mathscr{P}_{1}}$ the set of possible splits for the second cut

[TABLE]

and \mathcal{C}_{\mathscr{P}_{1}}\big{(}\theta_{2}^{(V)}\big{)}=\big{\{}(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}\textrm{ s.t. }j\in\theta_{2}^{(V)}\big{\}} when the split directions are restricted to $\theta_{2}^{(V)}\subset\{1,...,p\}$.

Lemma 6**.**

If Assumptions (A1) and (A3) are satisfied, then for all $(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}$, we have

[TABLE]

Lemma 7**.**

If Assumptions (A1) and (A3) are satisfied, then for all $(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}$, we have

[TABLE]

Lemma 8**.**

Consider that Assumptions (A1) and (A3) are satisfied. Let $\mathcal{C}_{1}\subset\{1,...,p\}\times\{1,...,q-1\}$ and $\mathcal{C}_{2}\subset\mathcal{C}_{\mathscr{P}_{1}}$ be two sets of splits of cardinality $c_{1}\geq 1$ and $c_{2}\geq 2$, such that the theoretical CART-splitting criterion is constant for all splits in $\mathcal{C}_{1}$ on one hand, and in $\mathcal{C}_{2}$ on the other hand, i.e.,

[TABLE]

where $H_{1}=\mathds{R}^{p}$ and $H_{2}=H^{\star}(\mathscr{P}_{1})$. Let $(j_{1},r_{1})\in\mathcal{C}_{1}$, $(j_{2},r_{2})\in\mathcal{C}_{2}$, and let $\textbf{L}_{n,\mathscr{P}_{2}}^{(\mathcal{C}_{1},\mathcal{C}_{2})}$ a the random vector where each component is the difference between the empirical CART-splitting criterion for the splits $(j,r)\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$ and $(j_{1},r_{1})$ for the first $c_{1}-1$ components, and for the splits $(j,r)\in\mathcal{C}_{2}\setminus(j_{2},r_{2})$ and $(j_{2},r_{2})$ for the remaining $c_{2}-1$ components, that is

[TABLE]

(a)

If $L^{\star}_{\mathcal{C}_{1}}>0$ and $L^{\star}_{\mathcal{C}_{2}}>0$, then we have

[TABLE]

where for $l,l^{\prime}\in\{1,2\}$, for all $(j,r)\in\mathcal{C}_{l}\setminus(j_{l},r_{l})$, $(j^{\prime},r^{\prime})\in\mathcal{C}_{l^{\prime}}\setminus(j_{l^{\prime}},r_{l^{\prime}})$, each element of the covariance matrix $\Sigma$ is defined by $\Sigma_{(j,r,l),(j^{\prime},r^{\prime},l^{\prime})}=\textrm{Cov}[Z_{j,r,l},Z_{j^{\prime},r^{\prime},l^{\prime}}],$ with

[TABLE]

\mu_{L,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}<q_{r}^{\star(j)},\textbf{X}\in H_{l}\big{]}, \mu_{R,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}\geq q_{r}^{\star(j)},\textbf{X}\in H_{l}\big{]}. Besides, for all $l\in\{1,2\}$ and for all $(j,r)\in C_{l}$, $\mathds{V}[Z_{j,r,l}]>0$.

(b)

If $L^{\star}_{\mathcal{C}_{1}}=L^{\star}_{\mathcal{C}_{2}}=0$, then we have

[TABLE]

where $\mathbf{V}$ is a gaussian vector of covariance matrix $\textrm{Cov}[\mathbf{Z}]$. If $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are explicitly written $\mathcal{C}_{1}=\{(j_{k},r_{k})\}_{k\in J_{1}}$, and $\mathcal{C}_{2}=\{(j_{k},r_{k})\}_{k\in J_{2}}$, with $J_{1}=\{1,...,c_{1}+1\}\setminus 2$ and $J_{2}=\{2\}\cup\{c_{1}+2,...,c_{1}+c_{2}\}$, $\mathbf{Z}$ is defined, for $l\in\{1,2\}$ and $k\in J_{l}$, by

[TABLE]

where p_{L,k}=\mathds{P}\big{(}X^{(j_{k})}<q_{r_{k}}^{\star(j_{k})},\textbf{X}\in H_{l}\big{)}, p_{R,k}=\mathds{P}\big{(}X^{(j_{k})}\geq q_{r_{k}}^{\star(j_{k})},\textbf{X}\in H_{l}\big{)}, and $h_{\mathscr{P}_{2}}$ is a multivariate quadratic form defined as

[TABLE]

Besides, the variance of each component of $h_{\mathscr{P}_{2}}(\mathbf{V})$ is strictly positive.

(c)

If $L^{\star}_{\mathcal{C}_{1}}>0$ and $L^{\star}_{\mathcal{C}_{2}}=0$, then we have

[TABLE]

where $\mathbf{V}$ is a gaussian vector of covariance matrix $\textrm{Cov}[\mathbf{Z}]$, and $\mathbf{Z}$ is defined as, for $k\in J_{1}$,

[TABLE]

for $k\in J_{2}$,

[TABLE]

and $h^{\prime}_{\mathscr{P}_{2}}$ is a multivariate quadratic form defined as

[TABLE]

Besides, the variance of each component of $h^{\prime}_{\mathscr{P}_{2}}(\mathbf{V})$ is strictly positive.

(d)

$L^{\star}_{\mathcal{C}_{1}}=0$* and $L^{\star}_{\mathcal{C}_{2}}>0$. Symmetric to case (c). *

Definition 2** (Theoretical splitting procedure at children nodes).**

Let $\theta^{(V)}=(\theta_{1}^{(V)},\theta_{2}^{(V)},\cdot)\in\Omega^{(V)}$ be the sets of eligible variables to split the nodes of a theoretical random tree. The set of best theoretical cuts at the left children node along the variables in $\theta^{(V)}_{2}$ is defined as

[TABLE]

If \mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)} has multiple elements, then $(j_{2},r_{2})$ is randomly drawn with probability

[TABLE]

where \mathbb{P}\big{(}\mathscr{P}_{1}\in T^{\star}(\Theta)|\Theta^{(V)}=\theta^{(V)}\big{)} is defined from Definition 1, and $\Phi_{\mathscr{P}_{1},\theta^{(V)},(j_{2},r_{2})}$ is the cdf of the limit law defined in Lemma 8 for \mathcal{C}_{1}=\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} and \mathcal{C}_{2}=\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}.

Lemma 9**.**

If Assumptions (A1)-(A3) are satisfied, then for all $\theta^{(V)}\in\Omega^{(V)}$, we have

[TABLE]

Finally, the proof of Theorem 1 in the two-splits scenario is the same as in the single-split scenario.

C.3 Proofs of intermediate lemmas

Proof of Lemma 1.

Set $j\in\{1,...,p\}$, and $r\in\{1,...,q-1\}$. We define the marginal cumulative distribution function $F^{(j)}$ of $X^{(j)}$, F^{(j)}(x)=\mathbb{P}\big{(}X^{(j)}<x\big{)}, and $F_{n}^{(j)}$ the empirical c.d.f.

[TABLE]

We adapt an inequality from Serfling (2009) (section 2.3.2 page 75) to bound the following conditional probability for all $\varepsilon>0$

[TABLE]

Since $f$ is continuous and strictly positive, there exists three constants $c_{1},c_{2},\eta>0$ such that for all $x\in[q_{r}^{\star(j)},q_{r}^{\star(j)}+\eta]$, $c_{1}\leq f^{(j)}(x)\leq c_{2}$. Thus, for all $\varepsilon<\eta$, we have

[TABLE]

which leads to

[TABLE]

Consequently,

[TABLE]

For $n>1+\frac{1}{c_{1}\varepsilon}$, we can apply Hoeffding inequality to C.6,

[TABLE]

where $C=e^{2c_{1}\eta(1+2c_{1}\eta)}$. By definition, we have

[TABLE]

To bound the previous integral, we break it down in three parts. Since $f^{(j)}$ is bounded by $c_{2}$ on $[q_{r}^{\star(j)},q_{r}^{\star(j)}+\eta]$, for $n>1+\frac{1}{c_{1}\eta}$ we use inequality C.7 to get

[TABLE]

In the second integral, we introduce the following change of variable $u=\sqrt{2n}c_{1}\varepsilon$

[TABLE]

and therefore we can write

[TABLE]

From Assumption (A1), $\lim\limits_{n\to\infty}\frac{a_{n}}{n}=0$, and then

[TABLE]

The case \lim\limits_{n\to\infty}\sqrt{a_{n}}\mathbb{P}\big{(}\hat{q}_{n,r}^{(j)}\leq X_{1}^{(j)}<q_{r}^{\star(j)}\big{)}=0 is similar. ∎

C.3.1 Case 1: $\mathscr{P}_{1}$

Proof of Lemma 2.

Let $j\in\{1,...,p\}$, $r\in\{1,...,q-1\}$, and $H\subseteq\mathds{R}^{p}$ such that \mathbb{P}\big{(}\textbf{X}\in H,X^{(j)}<q^{\star(j)}_{r}\big{)}>0 and \mathbb{P}\big{(}\textbf{X}\in H,X^{(j)}\geq q^{\star(j)}_{r}\big{)}>0. Let

[TABLE]

that is

[TABLE]

where, for a generic hyperrectangle $H$, we define $N_{n}(H)=\sum_{i=1}^{a_{n}}\mathds{1}_{\textbf{X}_{i}\in H}$, and

[TABLE]

with the convention $\overline{Y}_{H_{L}}=0$ if $H_{L}$ is empty. The theoretical quantities $H_{L}^{\star}$ and $\overline{Y}_{H_{L}^{\star}}$ are defined similarly by replacing the empirical quantile by its population version. We define symmetrically $H_{R}$, $H_{R}^{\star}$, $\overline{Y}_{H_{R}}$, $\overline{Y}_{H_{R}^{\star}}$.

Simple calculations show that

[TABLE]

The first term in equation (C.8) can be rewritten as

[TABLE]

Since $Y_{i}\in\{0,1\}$, we have the following bound

[TABLE]

and finally

[TABLE]

where

[TABLE]

A close inspection of the terms inside the sum of (C.10) reveals that

[TABLE]

which tends to zero, according to Lemma 1. Thus, in probability,

[TABLE]

Regarding the remaining terms in inequality (C.9), by the law of large numbers, in probability,

[TABLE]

Additionally,

[TABLE]

which tends to zero, according to Lemma 1. Therefore, in probability,

[TABLE]

Since $\mathbb{P}(\textbf{X}\in H)>0$ and $\mathbb{P}(\textbf{X}\in H_{L}^{\star})>0$ by assumption, we can combine (C.11)-(C.13) to obtain, in probability,

[TABLE]

Using (C.11) and (C.14) and inequality (C.9), we obtain, in probability,

[TABLE]

Similar results can be derived for the other term in equation (C.8), which allows us to conclude that, in probability,

[TABLE]

∎

Proof of Lemma 3.

Let $j\in\{1,...,p\}$, $r\in\{1,...,q-1\}$ and $H\subseteq\mathds{R}^{p}$ such that \mathbb{P}\big{(}\textbf{X}\in H,X^{(j)}<q^{\star(j)}_{r}\big{)}>0 and \mathbb{P}\big{(}\textbf{X}\in H,X^{(j)}\geq q^{\star(j)}_{r}\big{)}>0.

[TABLE]

From the law of large number, in probability,

[TABLE]

Thus, according to Lemma 2, in probability,

[TABLE]

∎

Proof of Lemma 4.

We consider $\mathcal{C}_{1}$, a set of splits of cardinality $c_{1}\geq 2$ satisfying, for all $(j,r)\in\mathcal{C}_{1}$, L^{\star}\big{(}\mathds{R}^{p},q_{r}^{\star(j)}\big{)}\stackrel{{\scriptstyle\tiny{\rm def}}}{{=}}L^{\star}_{\mathcal{C}_{1}}. Fix $(j_{1},r_{1})\in\mathcal{C}_{1}$, we recall that

[TABLE]

Case (a): $L^{\star}_{\mathcal{C}_{1}}>0$

We first consider the following decomposition for $(j,r)\in\mathcal{C}_{1}$,

[TABLE]

where

[TABLE]

($\overline{Y}_{R}^{\star}$, $N_{n,R}^{\star}$ are defined symmetrically). Letting \mu_{L,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}<q_{r}^{\star(j)}\big{]} (and $\mu_{R,r}^{(j)}$ symmetrically), the first two terms of the last decomposition are standard variance estimates and we can write

[TABLE]

where

[TABLE]

Using the Central limit theorem, in probability,

[TABLE]

The same result holds for the second term of (C.17), and using Lemma 2 for the third term of (C.17), we get that, in probability,

[TABLE]

Finally,

[TABLE]

Using Equation (C.16), each component of $\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1})}$ writes, with $(j,r)\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$,

[TABLE]

We can apply the multivariate Central limit theorem and Slutsky’s theorem to obtain,

[TABLE]

where for all $(j,r),(j^{\prime},r^{\prime})\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$, each element of the covariance matrix $\Sigma$ is defined by $\Sigma_{(j,r),(j^{\prime},r^{\prime})}=\textrm{Cov}[Z_{j,r},Z_{j^{\prime},r^{\prime}}],$ with

[TABLE]

Since $L^{\star}_{\mathcal{C}_{1}}>0$, we have for all $(j,r)\in\mathcal{C}_{1}$, $\mu_{L,r}^{(j)}\neq\mu_{R,r}^{(j)}$. Besides, according to assumption (A3), X has a strictly positive density. Consequently, the variance of $Z_{j,r}$ is strictly positive. This concludes the first case.

Case (b): $L^{\star}_{\mathcal{C}_{1}}=0$

Fix $(j,r)\in\mathcal{C}_{1}$. Since L^{\star}\big{(}\mathds{R}^{p},q_{r}^{\star(j)}\big{)}=0, we have

[TABLE]

Then, simple calculations show that L_{a_{n}}\big{(}\mathbb{R}^{p},\hat{q}_{n,r}^{(j)}\big{)} writes

[TABLE]

where

[TABLE]

($N_{n,R}$, $\overline{Y}_{R}$ are defined similarly for the other cell). Letting p_{L,r}^{(j)}=\mathds{P}\big{(}X^{(j)}<q_{r}^{\star(j)}\big{)} and p_{R,r}^{(j)}=\mathds{P}\big{(}X^{(j)}\geq q_{r}^{\star(j)}\big{)} with $p_{L,r}^{(j)},p_{R,r}^{(j)}>0$, we have

[TABLE]

where

[TABLE]

By the law of large numbers, $\lim\limits_{n\to\infty}\frac{N_{n,L}^{\star}}{a_{n}}=p_{L,r}^{(j)}$ in probability. Using Equation (C.13) in the proof of Lemma 2, it comes that, in probability, $\lim\limits_{n\to\infty}\frac{N_{n,L}}{a_{n}}=p_{L,r}^{(j)}$, and consequently $\lim\limits_{n\to\infty}\frac{a_{n}N_{n,L}}{N_{n,L}^{\star 2}}=\frac{1}{p_{L,r}^{(j)}}$. Since $\sqrt{a_{n}}\frac{1}{a_{n}}\sum_{i=1}^{a_{n}}(Y_{i}-\mu)\mathds{1}_{X_{i}^{(j)}<q^{\star(j)}_{r}}$ converges in distribution to a normal distribution by the Central limit theorem,

[TABLE]

Furthermore, as for Equation (C.10) in the proof of Lemma 2,

[TABLE]

and

[TABLE]

According to Lemma 1, the right hand side term converges to [math]. Then, in probability, $\lim\limits_{n\to\infty}\varepsilon_{L}=0$. Additionally, $\lim\limits_{n\to\infty}\frac{a_{n}^{2}}{N_{n,L}N_{n,L}^{\star}}=\frac{1}{p_{L,r}^{(j)2}}$, and then, in probability,

[TABLE]

The second term of $a_{n}R_{L,r}^{(j)}$ writes

[TABLE]

where in probability, $\lim\limits_{n\to\infty}2\frac{N_{n,L}}{a_{n}}=p_{L,r}^{(j)}$, $\lim\limits_{n\to\infty}\sqrt{a_{n}}(\overline{Y}_{L}^{\star}-\overline{Y}_{L})=0$ according to equation C.19, and $\sqrt{a_{n}}(\overline{Y}_{L}^{\star}-\mu)$ converges to a normal random variable from the central limit theorem. By Slutsky theorem, in probability, $\lim\limits_{n\to\infty}-a_{n}\times 2\frac{N_{n,L}}{a_{n}}(\overline{Y}_{L}^{\star}-\overline{Y}_{L})(\overline{Y}_{L}^{\star}-\mu)=0$. Finally for the third term of $a_{n}R_{L,r}^{(j)}$ we also use equation C.19 to conclude that in probability

[TABLE]

Consequently,

[TABLE]

Symmetrically, we also have

[TABLE]

with $\lim\limits_{n\to\infty}a_{n}R_{R,r}^{(j)}=0$, in probability.

Each component of $\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1})}$ writes, with $(j,r)\in\mathcal{C}_{1}\setminus(j_{1},r_{1})$,

[TABLE]

We explicitly write $\mathcal{C}_{1}=\{(j_{k},r_{k})\}_{k=1,...,c_{1}}$. Then $\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1})}$ can be decomposed as

[TABLE]

where for $k\in\{1,...,c_{1}\}$,

[TABLE]

$h_{\mathscr{P}_{1}}$ is a multivariate quadratic form defined as

[TABLE]

and $R_{n,\mathscr{P}_{1},k}=R_{L,r_{k}}^{(j_{k})}+R_{R,r_{k}}^{(j_{k})}-R_{L,r_{1}}^{(j_{1})}-R_{R,r_{1}}^{(j_{1})}$.

From the multivariate central limit theorem, $\mathbf{V}_{n}\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{\mathscr{D}}{\longrightarrow}\mathbf{V}$, where $\mathbf{V}$ is a gaussian vector of covariance matrix $\textrm{Cov}[\mathbf{Z}]$, and $\mathbf{Z}$ is defined as, for $k\in\{1,...,c_{1}\}$,

[TABLE]

with the simplified notations $p_{L,k}=p_{L,r_{k}}^{(j_{k})}$ and $p_{R,k}=p_{R,r_{k}}^{(j_{k})}$.

Finally, since $\lim\limits_{n\to\infty}\mathbf{R}_{n,\mathscr{P}_{1}}=\mathbf{0}$ in probability, from Slutsky’s theorem and the continuous mapping theorem, $a_{n}\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1})}\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{\mathscr{D}}{\longrightarrow}h_{\mathscr{P}_{1}}(\mathbf{V})$. Note that, since X has a strictly positive density, each component of $h_{\mathscr{P}_{1}}(\mathbf{V})$ has a strictly positive variance. ∎

Proof of Lemma 5.

Consider a path $\mathscr{P}=(j_{1},r_{1},\cdot)$. Set $\theta^{(V)}=(\theta_{1}^{(V)},\cdot,\cdot)\in\Omega^{(V)}$, a realization of the randomization of the split direction. Recalling that the best split in a random tree is the one maximizing the CART-splitting criterion, condition on $\Theta^{(V)}=\theta^{(V)}$,

[TABLE]

We recall that, given $\theta^{(V)}$, we define the set of best theoretical cuts along the variables in $\theta^{(V)}_{1}$ as

[TABLE]

Obviously if $(j_{1},r_{1})\notin\theta_{1}^{(V)}\times\{1,...,q-1\}$, the probability to select $\mathscr{P}_{1}$ in the empirical and theoretical tree is null. In the sequel, we assume that $(j_{1},r_{1})\in\theta_{1}^{(V)}\times\{1,...,q-1\}$ and distinguish between four cases: $(j_{1},r_{1})$ is not among the best theoretical cuts \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}, is the only element in \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}, is one element of \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} with a positive value of the theoretical CART-splitting criterion, or finally, is one element of \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} that all have a null value of the theoretical CART-splitting criterion.

Case 1

We assume that (j_{1},r_{1})\notin\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}. By definition of the theoretical random forest,

[TABLE]

Let \big{(}j^{\star},r^{\star}\big{)}\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}, thus

[TABLE]

Using equation (C.20), we have:

[TABLE]

Therefore, according to Lemma 3,

[TABLE]

Case 2

We assume that \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}=\big{\{}(j_{1},r_{1})\big{\}}. By definition of the theoretical random forest,

[TABLE]

Conditional on $\Theta^{(V)}=\theta^{(V)}$,

[TABLE]

which leads to

[TABLE]

From Lemma 3, for all $j\in\theta^{(V)}_{0}$, $r\in\{1,...,q-1\}$ such that $(j,r)\neq(j_{1},r_{1})$, in probability,

[TABLE]

Using inequality (C.23) and equation (C.24), we finally obtain,

[TABLE]

Case 3

We assume that (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}, \big{|}\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}\big{|}>1, and L^{\star}\big{(}\mathbb{R}^{p},q_{r_{1}}^{\star(j_{1})}\big{)}>0. On one hand, conditional on $\Theta^{(V)}=\theta^{(V)}$,

[TABLE]

On the other hand, conditional on $\Theta^{(V)}=\theta^{(V)}$,

[TABLE]

Combining the two previous inclusions,

[TABLE]

Using the same reasoning as in Case 2, we get

[TABLE]

We define the random vector $\textbf{L}_{n,\mathscr{P}_{1}}^{(\mathcal{C}_{1}^{\star})}$ where each component is the difference between the empirical CART-splitting criterion for the splits $(j,r)\in\mathcal{C}_{1}^{\star}\setminus(j_{1},r_{1})$ and $(j_{1},r_{1})$,

[TABLE]

then

[TABLE]

From Lemma 4 (case (a)),

[TABLE]

where for all $(j,r),(j^{\prime},r^{\prime})\in\mathcal{C}_{1}^{\star}\setminus(j_{1},r_{1})$, each element of the covariance matrix $\Sigma$ is defined by

[TABLE]

with

[TABLE]

\mu_{L,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}<q_{r}^{\star(j)}\big{]}, \mu_{R,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}\geq q_{r}^{\star(j)}\big{]}, and the variance of $Z_{j,r}$ is strictly positive. If $\Phi_{\theta_{1}^{(V)},(j_{1},r_{1})}$ is the c.d.f. of the multivariate normal distribution of covariance matrix $\Sigma$, we can conclude

[TABLE]

where

[TABLE]

According to Definition 1, in the theoretical random forest, if \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} has multiple elements, $(j_{1},r_{1})$ is randomly drawn with probability

[TABLE]

that is

[TABLE]

We can notice that, in the specific case where \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} has two elements, they are both selected with equal probability $\frac{1}{2}$. For more than two elements, the weights are not necessary equal, it depends on the covariance matrix $\Sigma$.

Case 4

We assume that all candidate splits have a null value for the theoretical CART-splitting criterion, i.e. for $(j,r)\in\theta_{1}^{(V)}\times\{1,...,q-1\}$, L^{\star}\big{(}\mathbb{R}^{p},q_{r}^{\star(j)}\big{)}=0. Consequently $\mathcal{C}_{1}^{\star}(\theta_{1}^{(V)})=\theta_{1}^{(V)}\times\{1,...,q-1\}$. By definition

[TABLE]

According to Lemma 4 (case (b)),

[TABLE]

where $\mathbf{V}$ is a gaussian vector of covariance matrix $\textrm{Cov}[\mathbf{Z}]$. If \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} is explicitly written \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}=\{(j_{k},r_{k})\}_{k=1,...,c_{1}}, $\mathbf{Z}$ is defined as, for $k\in\{1,...,c_{1}\}$,

[TABLE]

p_{L,k}=\mathds{P}\big{(}X^{(j_{k})}<q_{r_{k}}^{\star(j_{k})}\big{)}, p_{R,k}=\mathds{P}\big{(}X^{(j_{k})}\geq q_{r_{k}}^{\star(j_{k})}\big{)}, and $h_{\mathscr{P}_{1}}$ is a multivariate quadratic form defined as

[TABLE]

and the variance of each component of $h_{\mathscr{P}_{1}}(\mathbf{V})$ is strictly positive. If $\Phi_{\theta_{1}^{(V)},(j_{1},r_{1})}$ is the cdf of $h_{\mathscr{P}_{1}}(\mathbf{V})$, then as in Case 3,

[TABLE]

∎

C.3.2 Case 2: $\mathscr{P}_{2}$

Proof of Lemma 6.

Let $(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}$.

[TABLE]

Since $(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}$, \mathbb{P}\big{(}\textbf{X}\in H^{\star}(\mathscr{P}_{1})|X^{(j)}<q_{r}^{\star(j)}\big{)}>0 and \mathbb{P}\big{(}\textbf{X}\in H^{\star}(\mathscr{P}_{1})|X^{(j)}\geq q_{r}^{\star(j)}\big{)}>0. Then, we can directly apply Lemma 2 to the first term of this decomposition, which shows that, in probability

[TABLE]

We expand the second term

[TABLE]

with \hat{H}_{L}=\big{\{}\textbf{x}\in\hat{H}_{n}(\mathscr{P}_{1}):x^{(j)}<\hat{q}_{n,r}^{(j)}\big{\}}, H_{L}^{\star}=\big{\{}\textbf{x}\in H^{\star}(\mathscr{P}_{1}):x^{(j)}<\hat{q}_{n,r}^{(j)}\big{\}}, and for all $H\subseteq\mathds{R}^{p}$

[TABLE]

We define symmetrically $\hat{H}_{R}$ and $H_{R}^{\star}$. We obtain

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

We first consider $\Delta_{n,1}$. Simple calculations show that

[TABLE]

We consider the case $s_{1}=L$, ($s_{1}=R$ is similar). Since Y_{i}\in\big{\{}0,1\big{\}},

[TABLE]

As in the proof of Lemma 2, according to Lemma 1, $\lim\limits_{n\to\infty}\Delta_{n,1}=0$, in probability. Since $\Delta_{n,2}$ and $\Delta_{n,3}$ are the same quantities computed on each of the two daughter nodes, we study $\Delta_{n,2}$ only.

[TABLE]

Therefore

[TABLE]

As in the proof of Lemma 2, according to Lemma 1, $\lim\limits_{n\to\infty}\Delta_{n,2}=0$, in probability, which concludes the proof, since $\Delta_{n,3}$ can be studied in the same manner. ∎

Proof of Lemma 7.

Let $(j,r)\in\mathcal{C}_{\mathscr{P}_{1}}$.

[TABLE]

According to Lemma 6, the second term in equation (C.25) converges to [math] in probability. From the law of large numbers, in probability,

[TABLE]

which concludes the proof. ∎

Proof of Lemma 8.

Similar to the case with $\mathscr{P}_{1}$ (Lemma 3), where Lemma 6 is used instead of Lemma 2. ∎

Proof of Lemma 9.

Consider a path $\mathscr{P}_{2}=\{(j_{1},r_{1},L),(j_{2},r_{2},\cdot)\}$. Set \theta^{(V)}=\big{(}\theta_{1}^{(V)},\theta_{2}^{(V)}\big{)}, a realization of the randomization of the split directions at the root node and its left child node. Then, $\theta_{1}^{(V)}$ and $\theta_{2}^{(V)}$ denote the set of eligible variables for respectively the first and second split. We also consider \mathcal{C}_{\mathscr{P}_{1}}\big{(}\theta_{2}^{(V)}\big{)}\subset\mathcal{C}_{\mathscr{P}_{1}} the set of eligible second splits.

Recalling that the best split in a random tree is the one maximizing the CART-splitting criterion, conditional on $\Theta^{(V)}=\theta^{(V)}$,

[TABLE]

Recall that \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}=\underset{(j,r)\in\theta_{1}^{(V)}\times\{1,...,q-1\}}{\textrm{argmax}}L^{\star}\big{(}\mathbb{R}^{p},q_{r}^{\star(j)}\big{)}, and similarly

[TABLE]

Obviously if $(j_{1},r_{1})\notin\theta_{1}^{(V)}\times\{1,...,q-1\}$ or $(j_{2},r_{2})\notin\mathcal{C}_{\mathscr{P}_{1}}(\theta_{2}^{(V)})$, the probability to select $\mathscr{P}_{2}$ in the empirical and theoretical tree is null. In the sequel, we assume that $(j_{1},r_{1})\in\theta_{0}^{(V)}\times\{1,...,q-1\}$ and (j_{2},r_{2})\in\mathcal{C}_{\mathscr{P}_{1}}\big{(}\theta_{2}^{(V)}\big{)} and distinguish between cases, depending on whether (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} or not and (j_{2},r_{2})\in\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)} or not, as well as the cardinality of \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} and \mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}, and whether the maximum of the theoretical CART-splitting criterion is null or not.

Case 1

We assume that (j_{1},r_{1})\notin\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}. Hence, the theoretical decision tree satisfies

[TABLE]

According to Lemma 5, we have

[TABLE]

Case 2

We assume that \big{(}j_{2},r_{2}\big{)}\notin\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}. Again, for the theoretical decision tree,

[TABLE]

Letting \big{(}j^{\star},r^{\star}\big{)}\in\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)},

[TABLE]

Therefore,

[TABLE]

Consequently, according to Lemma 7,

[TABLE]

Case 3

We assume that (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} and \mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}=\{(j_{2},r_{2})\}, i.e. $(j_{2},r_{2})$ is the unique maximum of the theoretical CART-splitting criterion for the cell $H^{\star}(\mathscr{P}_{1})$. By definition of the theoretical decision tree,

[TABLE]

Conditional on $\{\Theta^{(V)}=\theta^{(V)}\}$,

[TABLE]

Consequently,

[TABLE]

For (j,r)\in\mathcal{C}_{\mathscr{P}_{1}}\big{(}\theta_{2}^{(V)}\big{)}\setminus(j_{2},r_{2}),

[TABLE]

Thus, using inequalities (C.27) and (C.28), and according to Lemma 7,

[TABLE]

and thus, using (C.26) and (C.27),

[TABLE]

where the second inequality is a direct consequence of Lemma 5.

Case 4

For the first split, we assume (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} with L^{\star}\big{(}\mathds{R}^{p},q^{\star(j_{1})}_{r_{1}}\big{)}>0, and for the second split, \big{(}j_{2},r_{2}\big{)}\in\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)} with |\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}\big{|}>1 and L^{\star}\big{(}H^{\star}(\mathscr{P}_{1}),q^{\star(j_{2})}_{r_{2}}\big{)}>0.

On one hand, conditional on the event $\{\Theta^{(V)}=\theta^{(V)}\}$,

[TABLE]

Using equation (C.29) to find a subset and a superset of $\left\{\mathscr{P}_{2}\in T(\Theta,\mathscr{D}_{n})\right\}$, we obtain

[TABLE]

We proved in Case 3 that the limit of the last two terms of the previous inequality is zero, in probability. Therefore,

[TABLE]

We define the random vector $\textbf{L}_{n,\mathscr{P}_{2}}^{(\mathcal{C}_{1}^{\star},\mathcal{C}_{2}^{\star})}$ (we drop $\theta^{(V)}$ to lighten notations) where each component is the difference between the empirical CART-splitting criterion for the splits $(j,r)\in\mathcal{C}_{1}^{\star}\setminus(j_{1},r_{1})$ and $(j_{1},r_{1})$ for the first $|\mathcal{C}_{1}^{\star}|-1$ components, and for the splits $(j,r)\in\mathcal{C}_{2}^{\star}\setminus(j_{2},r_{2})$ and $(j_{2},r_{2})$ for the remaining $|\mathcal{C}_{2}^{\star}|-1$ components, i.e.,

[TABLE]

Then, we can write

[TABLE]

According to Lemma 8,

[TABLE]

where for $l,l^{\prime}\in\{1,2\}$, for all $(j,r)\in\mathcal{C}_{l}^{\star}\setminus(j_{l},r_{l})$, $(j^{\prime},r^{\prime})\in\mathcal{C}_{l^{\prime}}^{\star}\setminus(j_{l^{\prime}},r_{l^{\prime}})$, each element of the covariance matrix $\Sigma$ is defined by $\Sigma_{(j,r,l),(j^{\prime},r^{\prime},l^{\prime})}=\textrm{Cov}[Z_{j,r,l},Z_{j^{\prime},r^{\prime},l^{\prime}}],$ with

[TABLE]

\mu_{L,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}<q_{r}^{\star(j)},\textbf{X}\in H_{l}\big{]}, \mu_{R,r}^{(j)}=\mathbb{E}\big{[}Y|X^{(j)}\geq q_{r}^{\star(j)},\textbf{X}\in H_{l}\big{]}, and the variance of $Z_{j,r,l}$ is strictly positive.

Letting $\Phi_{\mathscr{P}_{1},\theta^{(V)},(j_{2},r_{2})}$ be the c.d.f. of the multivariate normal distribution with covariance matrix $\Sigma$, and using equalities (C.30) and (C.31),

[TABLE]

We can check that

[TABLE]

In the theoretical random forest, the first cut $(j_{1},r_{1})$ is randomly selected with probability \mathbb{P}\big{(}\mathscr{P}_{1}\in T^{\star}(\Theta)|\Theta^{(V)}=\theta^{(V)}\big{)} (see the proof of Lemma 5). For the second cut, according to Definition 2, if \mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)} has multiple elements, $(j_{2},r_{2})$ is randomly drawn with probability

[TABLE]

Since the random selection at the root node of the tree and its children nodes are independent in the theoretical algorithm, $\mathscr{P}_{2}$ is selected with probability

[TABLE]

Ultimately,

[TABLE]

Case 5

We assume that (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)} and (j_{2},r_{2})\in\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}, and that the theoretical CART-splitting criterion is null for both splits: L^{\star}\big{(}\mathds{R}^{p},q^{\star(j_{1})}_{r_{1}}\big{)}=0 and L^{\star}\big{(}H^{\star}(\mathscr{P}_{1}),q^{\star(j_{2})}_{r_{2}}\big{)}=0.

Consequently \mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}=\theta_{1}^{(V)}\times\{1,...,q-1\}, and \mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}=\mathcal{C}_{\mathscr{P}_{1}}\big{(}\theta_{2}^{(V)}\big{)}. Using the same notations defined in Case 4, we have by definition

[TABLE]

According to Lemma 8 (case (b)),

[TABLE]

where $\mathbf{V}$ is a gaussian vector of covariance matrix $\textrm{Cov}[\mathbf{Z}]$. If $\mathcal{C}_{1}^{\star}$ and $\mathcal{C}_{2}^{\star}$ are explicitly written $\mathcal{C}_{1}^{\star}=\{(j_{k},r_{k})\}_{k\in J_{1}}$, and $\mathcal{C}_{2}^{\star}=\{(j_{k},r_{k})\}_{k\in J_{2}}$, with $J_{1}=\{1,...,c_{1}+1\}\setminus 2$ and $J_{2}=\{2\}\cup\{c_{1}+2,...,c_{1}+c_{2}\}$, $\mathbf{Z}$ is defined as, for $l\in\{1,2\}$ and $k\in J_{l}$,

[TABLE]

p_{L,k}=\mathds{P}\big{(}X^{(j_{k})}<q_{r_{k}}^{\star(j_{k})},\textbf{X}\in H_{l}\big{)}, p_{R,k}=\mathds{P}\big{(}X^{(j_{k})}\geq q_{r_{k}}^{\star(j_{k})},\textbf{X}\in H_{l}\big{)}. $h_{\mathscr{P}_{2}}$ is a multivariate quadratic form defined as

[TABLE]

and the variance of each component of $h_{\mathscr{P}_{2}}(\mathbf{V})$ is strictly positive.

$\Phi_{\mathscr{P}_{1},\theta^{(V)},(j_{2},r_{2})}$ is now defined as the cdf of $h_{\mathscr{P}_{2}}(\mathbf{V})$, and the end of the proof is identical to Case 4. We conclude

[TABLE]

Case 6

We assume (j_{1},r_{1})\in\mathcal{C}_{1}^{\star}\big{(}\theta_{1}^{(V)}\big{)}, (j_{2},r_{2})\in\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)} and \big{|}\mathcal{C}_{2}^{\star}\big{(}\theta_{2}^{(V)}\big{)}\big{|}>1 as in Case 4, but either L^{\star}\big{(}\mathds{R}^{p},q^{\star(j_{1})}_{r_{1}}\big{)}=0 and L^{\star}\big{(}H^{\star}(\mathscr{P}_{1}),q^{\star(j_{2})}_{r_{2}}\big{)}>0, or L^{\star}\big{(}\mathds{R}^{p},q^{\star(j_{1})}_{r_{1}}\big{)}>0 and L^{\star}\big{(}H^{\star}(\mathscr{P}_{1}),q^{\star(j_{2})}_{r_{2}}\big{)}=0.

The same reasoning than for Cases 4 and 5 applies where the limit law of $\textbf{L}_{n,\mathscr{P}_{2}}^{(\mathcal{C}_{1}^{\star},\mathcal{C}_{2}^{\star})}$ has both gaussian and $\chi$-square components and is given by case (c) or case (d) of Lemma 8.

∎

Appendix D Proof of Theorem 2

We recall Theorem 2 for the sake of clarity.

Theorem 2**.**

If $p_{0}\in[0,1]\setminus\mathcal{U}_{n}$ and $\mathscr{D}^{\prime}_{n}=\mathscr{D}_{n}$, then, conditional on $\mathscr{D}_{n}$, we have

[TABLE]

In addition for $p_{0}<\max\mathcal{U}_{n}$,

[TABLE]

where $\Phi(Mp_{0},M,p_{n}(\mathscr{P}))$ is the cdf of a binomial distribution with parameter $p_{n}(\mathscr{P})$, $M$ trials, evaluated at $Mp_{0}$, and, for all $\mathscr{P},\mathscr{P}^{\prime}\in\Pi$,

[TABLE]

and

[TABLE]

Let $p_{0}\in[0,\max~{}\mathcal{U}_{n})\setminus\mathcal{U}_{n}$ and $\mathscr{D}_{n}^{\prime}=\mathscr{D}_{n}$. Before proving Theorem 2, we need the following two lemmas.

Lemma 10**.**

Let $F$ be the hypergeometric function. Then, for $(a,c)\in\mathbb{Z}^{2}$ and $\mathscr{P}\in\Pi$ such that $p_{n}(\mathscr{P})>p_{0}$, we have

[TABLE]

Lemma 11**.**

Let $\mathscr{P}^{\prime}\in\Pi$. For all $\mathscr{P}\in\Pi$ such that $p_{n}(\mathscr{P})>p_{0}$, we have

[TABLE]

and

[TABLE]

Symmetrically, for all $\mathscr{P}\in\Pi$ such that $p_{n}(\mathscr{P})\leq p_{0}$, we have

[TABLE]

We are now in a position to prove Theorem 2.

Proof of Theorem 2.

The first statement, identity (D.1), is proved similarly to Corollary 2, using the law of large numbers instead of Theorem 1. For the second statement, we first recall that, by definition,

[TABLE]

Taking the expectation conditional on $\mathscr{D}_{n}$ gives

[TABLE]

where $U_{M}=\underset{\mathscr{P}\in\Pi}{\sum}\mathds{1}_{\hat{p}_{M,n}(\mathscr{P})>p_{0}\cap\hat{p}_{M,n}^{\prime}(\mathscr{P})\leq p_{0}}$, $V_{M}=\underset{\mathscr{P}\in\Pi}{\sum}\mathds{1}_{\hat{p}_{M,n}(\mathscr{P})>p_{0}}$, and $V_{M}^{\prime}=\underset{\mathscr{P}\in\Pi}{\sum}\mathds{1}_{\hat{p}_{M,n}^{\prime}(\mathscr{P})>p_{0}}$. Note that

[TABLE]

Also,

[TABLE]

For all $\mathscr{P}\in\Pi$,

[TABLE]

where $\Phi$ is the cdf of the binomial distribution. As a direct consequence of Lemma 11,

[TABLE]

which yields

[TABLE]

This is the desired result. ∎

D.1 Proof of intermediate lemmas

Proof of lemma 10.

Cvitković et al. (2017) provides an asymptotic expansion of the hypergeometric function $F$ in the case where the first and third parameters goes to infinity with a constant ratio. For $a,c,z,\varepsilon\in\mathds{R}$, $b\notin\mathbb{Z}\setminus\mathbb{N}$, such that $\varepsilon>1$, and $z\varepsilon<1$, Cvitković et al. (2017) gives in the section 2.2.2 (end of page 10)

[TABLE]

We can then derive the limit of the following ratio

[TABLE]

We use D.3 in the specific case where $b=1$, $a,c\in\mathbb{Z}$, $\varepsilon=\frac{1}{1-p_{0}}>1$, $z=1-p_{n}(\mathscr{P})$ for $\mathscr{P}\in\Pi$ such that $p_{n}(\mathscr{P})>p_{0}$ (and then $z\varepsilon<1$), and $\lambda=M(1-p_{0})$, if follows that

[TABLE]

∎

Proof of lemma 11.

Fix $\mathscr{D}_{n}$. Let $\mathscr{P}^{\prime},\mathscr{P}\in\Pi$. In what follows, when there is no ambiguity, we will replace $T(\Theta,\mathscr{D}_{n})$ by $T_{n}(\Theta)$ to lighten notations.

Case 1: $p_{n}(\mathscr{P})>p_{0}$

[TABLE]

since

[TABLE]

because, conditional on $\mathscr{D}_{n}$, the events $\mathscr{P}^{{}^{\prime}}\in T_{n}(\Theta_{1}),\ldots,\mathscr{P}^{{}^{\prime}}\in T_{n}(\Theta_{M})$ are independent. We can rewrite,

[TABLE]

yielding, using equation (D.5),

[TABLE]

Besides, by definition of the correlation

[TABLE]

simple calculations show that

[TABLE]

which, together with equation (D.8) leads to,

[TABLE]

Regarding the probability in the right-hand side of equation (D.10), we have

[TABLE]

Using standard formulas, $\Phi$ can be expressed with the incomplete beta function,

[TABLE]

and the regularized beta function is related to the hypergeometric function $F$, for $a>0$, $b>0$, and $p\in[0,1]$ (Olver et al., 2010),

[TABLE]

Then, we can express the cdf of the binomial distribution using the hypergeometric function, and it follows

[TABLE]

According to Lemma 10,

[TABLE]

Consequently,

[TABLE]

and using this limiting result with equation (D.10) yields,

[TABLE]

Regarding the conditional variance,

[TABLE]

[TABLE]

where

[TABLE]

Then, we follow the same reasoning that leads to equation (D.13). We can fully expand $C_{M}$ using Bayes formula, depending whether $\mathscr{P}\in T_{n}(\Theta_{1})$ or $\mathscr{P}\in T_{n}(\Theta_{2})$. Note that, since all the trees are independent conditional on $\mathscr{D}_{n}$, we can reduce the conditioning event \big{\{}\mathscr{P}\in T_{n}(\Theta_{1}),\mathscr{P}\in T_{n}(\Theta_{2}),\hat{p}_{M,n}(\mathscr{P})\leq p_{0},\mathscr{D}_{n}\big{\}} to \big{\{}\mathscr{P}\in T_{n}(\Theta_{1}),\mathscr{P}\in T_{n}(\Theta_{2}),\mathscr{D}_{n}\big{\}}, then

[TABLE]

Conditional on $\mathscr{D}_{n}$, $T_{n}(\Theta_{1})$ and $T_{n}(\Theta_{2})$ are independent, then

[TABLE]

we can rewrite $C_{M}$

[TABLE]

where

[TABLE]

We first consider the term

[TABLE]

Equation (D.13) directly gives,

[TABLE]

On the other hand

[TABLE]

Again, as for equation (D.11), we can express the cdf of the binomial distribution using the hypergeometric function $F$

[TABLE]

and from Lemma 10,

[TABLE]

that is

[TABLE]

Using equations (D.15) and (D.17), we conclude

[TABLE]

We follow the same reasoning for $\Delta_{M,3}$, equation (D.13) gives

[TABLE]

On the other hand,

[TABLE]

From Lemma 10,

[TABLE]

And finally $\lim\limits_{M\to\infty}\Delta_{M,3}=0$. The term $\Delta_{M,2}$ can be treated in a similar way, since equation (D.13) gives

[TABLE]

Simple identity shows

[TABLE]

Taking the limit of the previous equation and using equations (D.17) and (D.19), we get

[TABLE]

Using (D.13) and (D.20), $\lim\limits_{M\to\infty}\Delta_{M,2}=0$. Since $\Delta_{M,1},\Delta_{M,2},\Delta_{M,3}\to 0$, we obtain $\lim\limits_{M\to\infty}C_{M}=0$, that is,

[TABLE]

Finally combining equations (D.14) and (D.21),

[TABLE]

Case 2: $p_{n}(\mathscr{P})\leq p_{0}$

By the law of large numbers, \lim\limits_{M\to\infty}\hat{p}_{M,n}\big{(}\mathscr{P}\big{)}=p_{n}(\mathscr{P}) in probability, and consequently \lim\limits_{M\to\infty}\mathds{P}\big{(}\hat{p}_{M,n}\big{(}\mathscr{P}\big{)}\leq p_{0}\big{)}=1. Additionally, we can simply write

[TABLE]

Again, by the law of large numbers, \lim\limits_{M\to\infty}\hat{p}_{M,n}\big{(}\mathscr{P}^{\prime}\big{)}=p_{n}(\mathscr{P}^{\prime}) in probability. Then, if $p_{n}(\mathscr{P}^{\prime})>p_{0}$, \lim\limits_{M\to\infty}\mathds{P}\big{(}\hat{p}_{M,n}\big{(}\mathscr{P}^{\prime}\big{)}>p_{0}\big{)}=1, and it follows that \lim\limits_{M\to\infty}\mathds{P}\big{(}\hat{p}_{M,n}(\mathscr{P}^{\prime})>p_{0},\hat{p}_{M,n}(\mathscr{P})\leq p_{0}|\mathscr{D}_{n}\big{)}=1. If $p_{n}(\mathscr{P}^{\prime})\leq p_{0}$, \lim\limits_{M\to\infty}\mathds{P}\big{(}\hat{p}_{M,n}\big{(}\mathscr{P}^{\prime}\big{)}>p_{0}\big{)}=0, and consequently \lim\limits_{M\to\infty}\mathds{P}\big{(}\hat{p}_{M,n}\big{(}\mathscr{P}^{{}^{\prime}}\big{)}>p_{0},\hat{p}_{M,n}(\mathscr{P})\leq p_{0}|\mathscr{D}_{n}\big{)}=0. This can be compacted under the form

[TABLE]

The proof for the case \mathds{P}\big{[}\hat{p}_{M,n}\big{(}\mathscr{P}^{{}^{\prime}}\big{)}>p_{0}|\hat{p}_{M,n}(\mathscr{P})>p_{0},\mathscr{D}_{n}\big{]} is similar. ∎