Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
40,840 characters · 12 sections · 24 citation commands
Does Residuals-on-Residuals Regression Produce Representative Estimates of Causal Effects?
Double Machine Learning (DML) bickel1993efficient,robins1994estimation,Chernozhukov2018-fl has become the standard method for estimating causal effects in large, high-dimensional datasets under conditional ignorability, which stipulates that the treatment is as good as randomly-assigned given observed covariates Imbens2004-ir. To strengthen this assumption, researchers in a wide variety of fields use the DML method to condition on rich covariates without also making restrictive functional form assumptions jares2025policy, chernozhukov2020causal,holtz2020interdependence.
DML encompasses a broad class of estimators. In this paper, we focus on the Partially Linear Model (PLM), which relates an outcome $Y_i$ to a (continuous- or discrete-valued) treatment $T_i$ conditional on pretreatment covariates $X_i$ as follows: \[ Y_i = \theta T_i + g(X_i) + e_i \quad \text{and} \quad T_i = h(X_i) + u_i. \] The PLM imposes minimal assumptions on how the treatment and outcome relate to covariates. It also motivates an intuitive two-step DML estimator of $\theta$, the residuals-on-residuals regression (RORR). RORR, a direct extension of the Frisch-Waugh-Lovell theorem, first “partials out” the effect of $X_i$ using flexible machine learning methods, then forms the residuals \[ \widetilde{Y}_i = Y_i - \widehat{g}(X_i) \quad \text{and} \quad \widetilde{T}_i = T_i - \widehat{h}(X_i), \] and lastly regresses $\widetilde{Y}_i$ on $\widetilde{T}_i$ to obtain the RORR estimate of $\theta$, which we will denote by $\hat{\theta}$ Robinson1988-dm.
When the treatment effect $\theta$ is the same for all units in the population, it is also the Average Treatment Effect (ATE) for binary treatments, the Average Causal Derivative (ACD) for continuous treatments, and the Average Incremental Effect (AIE) for integer-valued treatments. Alternatively, we study the probability limit (plim) and interpretation of $\hat{\theta}$ when treatment effects are heterogeneous. Under such heterogeneity, $\hat{\theta}$ converges to a conditional variance-weighted average of causal effects, which puts greater weight on units whose treatment values are less predictable. For example, when $X_i$ is discrete, the RORR estimand places the most weight on the treatment effects in the strata where the treatment is most variable Angrist1998-ok, Aronow2016-nn, Sloczynski2022-kg. When treatments are many-valued (for example, continuous), RORR may be subject to more nuanced biases documented in this paper.
We contribute a general analysis of RORR and its potential biases with binary and many-valued treatments. We demonstrate the empirical relevance of these biases both in a stylized numerical example and with real-world data from Netflix. We propose simple alternative estimators that coarsen the treatment into bins, establish the favorable theoretical properties of an estimator based on Augmented Inverse Propensity Weighting (AIPW) Cattaneo2010-oc, and apply this estimator to empirical data. The proposed Coarsened AIPW estimator has been used for causal decision-making at Netflix for several years lal2024doesregressionproducerepresentative, netflix2024round2 and is the default methodology used in the company's internal observational causal inference platform.
RORR is popular for many reasons. Most notably, it enables the use of highly flexible modern machine learning estimators for the nuisances in the PLM, thus freeing researchers from making strong assumptions about those parameters. If the PLM is correctly specified, it is also statistically efficient and achieves parametric convergence rates, even when the nuisances are estimated at slower rates. Attesting to its wide applicability and popularity, recent applications of RORR can be found in economics BaiardiNaghi2024, ecology FinkEtAl2023, and public health WeiEtAl2024.
In industry settings---characterized by large datasets, short timelines, and stakeholders with diverse technical backgrounds---RORR is also commonly used for practical reasons. For example, the final regression step is computationally efficient, as only a small number of statistics is needed to compute $\hat{\theta}$. Furthermore, the recipe of (1) removing variation explainable by pretreatment covariates and (2) estimating the effect of the remaining exogenous variation in $T_i$ on $Y_i$ is intuitive and easy to explain to non-experts. In its simplest form, RORR is none other than linear regression.
Unfortunately, this simplicity comes at a cost. The interpretability of $\hat{\theta}$ as an estimate of the “average” treatment effect (whether the ATE, ACD, or AIE) depends on the assumption of a homogeneous treatment effect in the PLM, which may not hold in applications. In this section, we discuss the interpretation of $\hat{\theta}$ under two violations of this assumption: binary treatments with heterogeneous effects across individual units and many-valued treatments with continuous dose-response functions.
We begin by studying the RORR estimand for binary treatments with heterogeneous treatment effects. Letting $T_i \in \{0, 1\}$ denote the binary treatment, we consider the model: \[ Y_i = \theta_i T_i + g(X_i) + e_i \quad \text{and} \quad T_i = h(X_i) + u_i. \] where $\theta_i$ is an individual treatment effect. We assume conditional ignorability of treatment given $X_i$, which in turn implies that the errors are conditionally exogenous and uncorrelated: $E[e_i|X_i] = 0$, $E[u_i|X_i] = 0$, and $E[e_i u_i|X_i] = 0$. Thus, $g(X_i) = E[Y_i - \theta_i T_i | X_i]$ and $h(X_i) = E[T_i | X_i]$. Conditional ignorability also implies $\theta_i$ is conditionally independent of $T_i$ given $X_i$, $\theta_i \protect\mathpalette{\protect\independenT}{\perp} T_i | X_i$. Lastly, we assume non-degeneracy of the treatment distribution, $E[(T_i - h(X_i))^2] > 0$.
We consider the plim of the OLS regression of $Y_i - g(X_i)$ on $T_i - h(X_i)$ with iid observations $(Y_i, T_i, X_i)$, $i = 1, \ldots N$. To be clear, $g$ and $h$ must be estimated in practice. We assume consistent estimators for these and focus on the (true) limiting $g$ and $h$ to clarify that the RORR plim is biased relative to the ATE even when the researcher has access to sufficient covariates and consistently estimates the nuisance parameters.\footnote{Note that researchers often use extremely flexible function classes for $g$ and $h$ (e.g., gradient boosted trees or deep neural networks) that are able to approximate the true nuisance functions arbitrarily closely.}
First, observe that
Using the fact that $T_i$ is binary and applying the law of iterated expectations, we rewrite the above as:
where Equation ((ref)) follows from $\theta_i$ being conditionally independent of $T_i$. This demonstrates the well-known result that, with a binary treatment, linear regression converges to a conditional variance-weighted average of individual treatment effects Angrist1998-ok,Aronow2016-nn,Sloczynski2022-kg.
For an intuitive restatement, denote the conditional variance weights by $\omega_i := \frac{(T_i - h(X_i))^2}{E[(T_i - h(X_i))^2]}$ and note that $E[\omega_i] = 1$ by construction. Then the bias of the RORR plim (which we will denote by $\tilde{\theta}$) with respect to the ATE can be written as:
In other words, the RORR bias for binary treatments when treatment effects are heterogeneous is equal to the covariance of individual treatment effects with the normalized residual variance of the treatment. This covariance will not equal zero except in special cases (e.g., the treatment is assigned uniformly at random) and therefore $\tilde{\theta} \neq E[\theta_i]$ in general.\footnote{A corollary is that ranking treatments based on their PLM coefficient is not the same as ranking them based on their ATEs lal2024doesregressionproducerepresentative.}
We now turn our attention to many-valued (e.g., continuous or integer-valued) treatments with continuous dose-response functions. Although past research has studied the effect of treatment effect heterogeneity on the interpretation of linear treatment effect estimators, it has primarily done so in the context of binary treatments and/or linear treatment effects Aronow2016-nn. However, in many applications, treatments are continuous and/or have nonlinear effects on the outcome (for example, they may have diminishing returns). Such nonlinearity is an important form of treatment effect heterogeneity Yitzhaki1996-io, Angrist1999-sp. Here, we present a novel bias decomposition for RORR with many-valued treatments that emphasizes its differences and similarities with the binary treatment case.
Specifically, we study the model:
where $f$ is a twice continuously differentiable function. As before, we assume conditional ignorability, consistent estimators for $g$ and $h$, non-degenerate treatment values, and iid observations.
Under these assumptions, the RORR estimate converges in probability to:
Since $h(X_i)$ is a constant given $X_i$ and applying the law of iterated expectations, we can rewrite the above as:
By the mean value theorem and continuity of $f$, there exists a $T_i^*$ between $T_i$ and $h(X_i)$ for all $X_i$ such that:
showing that, as in the binary treatment setting, $\hat{\theta}$ also converges to a conditional variance-weighted average of causal effects.\footnote{Not coincidentally, this representation of the RORR estimand closely resembles the representation of the Wald estimand with a binary instrument and continuous endogenous treatment as a first-stage effect-weighted average of derivatives at the mean values $T_i^*$ Angrist1999-sp.} However, unlike in the binary treatment case, the quantity being averaged cannot be interpreted as the causal effect of increasing the treatment in the population represented by the sample. This is because the mean value $T_i^*$ is not the actual treatment dose received by $i$, but a convex combination of the received treatment $T_i$ and its conditional mean. As such, $T_i^*$ may not be an observed treatment level. If $T_i$ is not continuous, it may not even be a realizable treatment value.
Proposition (ref) establishes the restrictive conditions under which $\tilde{\theta}$ converges to the ACD.
In other words, $A$ is bounded by a term that depends on the curvature of $f$, such that $\tilde{\theta}$ will be closer to the conditional variance-weighted average of $f'$ when $f$ is close to affine. The bias of the conditional variance-weighted average of $f'$ is eliminated when there is no treatment heterogeneity, which holds trivially if $f$ is affine. Therefore, $A$ and $B$ vanish when $f$ is affine (and the PLM is therefore correctly specified). However, if $f$ is not affine, then the biases do not vanish except in contrived cases (e.g., the biases cancel exactly).
To help build intuition, this section presents a stylized numerical example. Replication code for all figures and tables in this section can be found at \url{https://github.com/winston-chou/rorr}. Although we make simplifying assumptions to facilitate closed-form analysis, our choices are also intended to reflect qualitative aspects of real-world data. In particular, we assume that:
Let $X_i$ be a Categorical variable that takes on values $j = 1, \ldots, J$ with probabilities $\pi_1, \ldots, \pi_J$ and $T_i$ be conditionally Poisson given $X_i$ with parameters $\lambda_1, \ldots, \lambda_J$. Let $f(T_i) = \log(T_i + 1)$. This allows us to derive the following analytical expression for the conditional expected derivative of $Y_i$ with respect to $T_i$ given $X_i$ (see Appendix (ref)):
The ACD is then just $\sum_j \pi_j \frac{1 - \exp(-\lambda_j)}{\lambda_j}$.
We can also derive the RORR plim analytically as:
where, as before, $T_i^*$ is a point between $T_i$ and $\lambda_j$. Note the two biases relative to the ACD. First, rather than evaluate $f'$ at $T_i$, we evaluate it at $T_i^*$. Second, we also weight each $f'(T_i^*)$ by its normalized conditional variance.
Figure (ref) illustrates the resulting bias by simulating this data-generating process. First, we plot $f(T_i) = \log(T_i + 1)$ in top panel of Figure (ref), as well as tangent lines with slopes equal to $E[f'(T_i)]$ in blue and to $E[\omega_i f'(T_i^*)]$ in red, where $T_i^*$ is the “effective” treatment analyzed by RORR. The key takeaway is that RORR targets a quantity other (and smaller) than the ACD.\footnote{An analogous result in the welfare economics literature is that OLS up-weights the slopes of higher-income groups in regressions of consumption on income, leading to attenuation Yitzhaki1996-io.} The subsequent panels give intuition for this result: After weighting by $\omega_i$ and transforming $T_i$ to $T_i^*$, the effective treatment distribution is much more right-skewed than the observed treatment distribution. This means that we tend to evaluate the slope of $f$ at higher values of $T_i$. This leads to negative bias because $f''(T_i) < 0$.
In Table (ref), we report the estimated empirical RORR from simulations at varying sample sizes. For comparison, we also report the empirical ACD (calculated as the sample mean of $\frac{1}{T_i + 1}$) and the true ACD computed using ((ref)). Note that, because $T_i$ is integer-valued in this example, the more appropriate causal estimand is the Average Incremental Effect (AIE), defined as:
where $p$ is the mass function of $T_i$. However, because the RORR plim is a weighted average of derivatives, we focus on the ACD in Table (ref) and propose a consistent estimator of the AIE in Section (ref). As Table (ref) shows, the RORR plim is negatively biased for the ACD. This is because it places more weight on the derivative of the dose-response curve at larger values of the treatment, where the dose-response curve tends to be flatter.
Given how common right-skewed treatments and diminishing dose-response curves are in practice, our analysis suggests that, as a rule of thumb, the RORR estimate will tend to have a downward bias relative to the ACD. Below, we propose a consistent estimator for the ACD and provide additional evidence for this rule using real-world data from Netflix.
We have shown that residuals-on-residuals regression (RORR) with many-valued treatments targets a conditional variance-weighted average of derivatives, which does not in general equal the Average Causal Derivative (ACD). A natural remedy is to coarsen the treatment into bins, thus approximating the dose-response function by a step function. In this section, we formalize two such estimators: Coarsened RORR and Coarsened AIPW. Both proceed by partitioning the support of the treatment into $K$ disjoint intervals, estimating treatment effects between bins, then aggregating these to form a coarsened estimate of the ACD.
Assume $T_i$ is distributed on a compact interval $[\underline{t}, \overline{t}]$ of length $C$. Let $\{S_1,\dots,S_K\}$ be an evenly-spaced partition of this interval with $S_k=[t_k,t_{k+1})$, $k = 1, K - 1$. Denote by $\ell = t_{k+1} - t_k$ and $\overline{t}_k = \frac{t_{k+1} + t_k}{2}$ the length and midpoint of each bin, respectively. Define bin indicators $D_{ik} = \mathbf 1\{T_i \in S_k\}$ and probabilities $p_k(X_i)=\Pr(T_i\in S_k\mid X_i)$. Let $m_k(X_i)=\mathbb E[Y_i\mid T_i\in S_k, X_i]$.
We first consider the Coarsened RORR of the residualized outcome $\widetilde Y_i$ on the residualized bin indicators $\widetilde D_{ik} = D_{ik} - p_k(X_i)$ for $k > 1$ (omitting the first bin to avoid multicollinearity). The Coarsened RORR estimate of the ACD is given by:
where $\hat{\beta}_j$ is the estimated regression coefficient corresponding to $\widetilde{D}_{ij}$ and the weights $w_k$ are proportional to the fraction of treatment values in $S_k$: \[ w_k =
\]
For ease of exposition, we will focus on the plim of $\hat{\psi}_{CR}$ in the simple case where $K = 2$ (i.e., the treatment is binarized by choosing a cutpoint $c$ and setting $D_{i2} = 0$ if $T_i < c$ and 1 otherwise). Then $\hat{\psi}_{CR}$ is given by the estimated regression coefficient on $D_{i2}$, which converges to
where $\upsilon_i := \frac{p_2(X_i)(1 - p_2(X_i))}{E[p_2(X_i)(1 - p_2(X_i))]}$. There are two weights that distort $\beta_2$ from the ACD. First, $\upsilon_i$ corresponds to the usual conditional variance weight: strata of $X_i$ whose treatment values are more evenly distributed about $c$ will receive greater weight. The second weight is more subtle. A first-order Taylor expansion of $E[f(T_i) | T_i > c, X_i] - E[f(T_i) | T_i < c, X_i]$ yields $f'(c) \delta(x)$, where $\delta(x) := E[T_i | T_i > c, X_i = x] - E[T_i | T_i < c, X_i = x]$ is a measure of the different $T_i$ tends to be from $c$ within a given stratum.
Thus, the Coarsened RORR plim in this setting is interpretable as a weighted average of derivatives at the cutoff $c$. These weights can lead to counterintuitive biases. In particular, strata of $X_i$ for which the cutoff $c$ is less similar to the average treatment values on either side of $c$ are given more weight when estimating $f'(c)$.
A common benchmark for estimating the ACD with continuous treatments is the Generalized Propensity Score (GPS) Imbens2000-sk. However, GPS requires estimating the conditional density of the treatment, which can suffer from slow rates and instability without parametric assumptions kennedy2017non.
As an alternative to both RORR and GPS, we propose a simple coarsened Augmented Inverse Propensity Weighting (AIPW) estimator, which uses AIPW estimates of counterfactual means as building blocks Cattaneo2010-oc. As before, this estimator proceeds by first partitioning the support of $T_i$ into $K$ disjoint bins of length $\ell$. For each bin $S_k$, we estimate the bin-level propensity score $p_k(X_i)$, for example by fitting a multiclass classification model. Denote this estimate by $\hat{p}_k(X_i)$. We also estimate a flexible outcome regression for $m_k(X_i)$, which we denote by $\hat{m}_k(X_i)$.
Next, we form the usual AIPW estimator for the marginal counterfactual mean in bin $S_k$: \[ \hat{\psi}_k := \frac{1}{N}\sum_{i=1}^N \left(\left[\frac{\mathbf{1}(T_i \in S_k)}{\hat{p}_k(X_i)}(Y_i - \hat{m}_k(X_i))\right] + \hat{m}_k(X_i)\right). \] This estimator has the double-robustness property, meaning that if either $\hat{p}_k$ or $\hat{m}_k$ is consistently estimated, then $\hat{\psi}_k$ is also consistent for an average potential outcome within the bin $S_k$ kennedy2024semiparametric. Note that the plim of $\hat{\psi}_k$ has a somewhat subtle interpretation due to averaging over the bin. Under strong ignorability, it can be interpreted as the population-average potential outcome if units are treated according to the conditional distribution of $T_i$ given $X_i$ and $T_i \in S_k$.
The Coarsened AIPW estimate of the ACD is given by: \[ \hat{\psi} := \sum_{k=1}^{K-1} w_k \left(\frac{\hat{\psi}_{k + 1} - \hat{\psi}_k}{\overline{t}_{k + 1} - \overline{t}_k}\right). \] Heuristically, $\hat{\psi}$ approximates $f'$ using a piecewise linear function and computes its weighted average using the empirical distribution of the lower segment.
Proposition (ref) proves the consistency of this estimator under regularity conditions. For the proof, we introduce the potential outcomes notation:
$\Delta_k$ has a subtle interpretation: it is the “effect” of moving from segment $S_k$ to $S_{k+1}$ when units are treated according to the conditional treatment distribution in each segment. It can also be interpreted as a deterministic ATE. Given continuity of $f$, $p$, and positivity, the mean value theorem for integrals states that there exists a $\tilde{t}_k \in S_k$ such that \[ \int f(t) \frac{p(t | x) 1(t \in S_k)}{\Pr(t \in S_k | x)} dt = f(\tilde{t}_k). \] Note that $\tilde{t}_k$ depends implicitly on $x$. Thus, we can rewrite the above as:
showing that $\Delta_k$ is the ATE of fixing $T_i$ to some $\tilde{t}_{k+1} \in S_{k+1}$ relative to $\tilde{t}_k \in S_k$.
Table (ref) shows the results of applying this estimator to the simulated data from Section (ref). Because the treatment is integer-valued, we can simply set the bins to each observed treatment value. As the table shows, this yields a consistent estimator for the AIE.\footnote{Note that the AIE is less than the ACD because $f(t + 1) - f(t) = \log\left(1 + \frac{1}{t+ 1}\right) \le \frac{1}{t + 1}$ for all $t \ge 0$.}
Our proof of the consistency of the coarsened AIPW estimator relies on conditional ignorability, positivity of the conditional density $p(t|x)$ and continuous differentiability of the dose-response function $f$ and conditional density. Although the validity of these assumptions should be assessed on a case-by-case basis, one advantage of the coarsened AIPW estimator is that it enables a suite of useful diagnostics austin2015moving. For example, a best practice for strengthening the conditional ignorability assumption, which we implement in our empirical application, is to show balance on pretreatment covariates before and after weighting by the estimated propensity scores.
In our experience, a common way in which these assumptions can be violated is if the treatment is extremely sparse and/or skewed in regions of the covariate space. Such sparsity can be diagnosed by checking the distribution of estimated propensity scores. Remedies for sparsity include trimming and/or grouping extreme treatment values petersen2012diagnosing. Indeed, in Appendix (ref), we provide a justification for choosing a fairly small number of bins (on the order of $N^{1/7}$, for example, five to 10 for a dataset with 1 million observations). Plotting the treatment distribution and estimated dose-response curve, as we do in our empirical application, can also help diagnose violations of positivity and continuity.
We now demonstrate the empirical relevance of our theoretical analysis using real-world data from Netflix. Although we are limited in what we can share for confidentiality reasons, the main thrust of this section is to show that the theoretical biases discussed above can (and, in our experience, often do) appear in real-world data.
In this particular application, we sought to understand how the use of a feature, which we will call Feature A, affects future visits to Netflix. To answer this question, we drew a random sample of 2,971,128 members and counted the number of times they used Feature A over a 28 day window. We then divided this number by the member's count of visits to Netflix over the same period to define our continuous treatment, Feature A Usage Rate. Next, we defined our outcome as the count of each member's visits to Netflix in the next 28 day window. As covariates, we included the count of times each member used Feature A and the count of times each member visited Netflix in the seven, 14, and 28 days preceding the treatment period.\footnote{We complemented these six covariates with an additional 25 covariates; most of these measured the usage of other Netflix features in the 28 days preceding to the treatment period.}
We divided our dataset into roughly equal-sized training, validation, and test datasets consisting of $\approx$980,000 units each. To estimate the nuisance parameters in the PLM, we fit gradient boosted regression trees to the treatment and outcome variables observed in the training dataset, using the validation dataset to tune the number of boosting rounds. Lastly, we regressed the outcome residuals on the treatment residuals in the test dataset to obtain the RORR treatment estimate, which is shown in Table (ref). As the table shows, the RORR estimate of the effect of Feature A on subsequent visits is small, negative, and statistically significant.\footnote{Note that treatment effects are reported after standardizing the treatment and outcome by their respective standard deviations.} This finding contradicted our prior belief that Feature A would increase visits to Netflix.
Our coarsened AIPW estimator helps provide intuition for this puzzling result. To fit this estimator, we first coarsened the treatment into five bins and then fit a multiclass classifier using gradient boosting to the resulting bins. We assigned zero values (i.e., no usage of Feature A during the treatment period) to the first bin and then divided the remaining non-zero values into quartiles. We reused the RORR outcome regression.
Figure (ref) is a standard diagnostic that plots the difference in the standardized pretreatment value of the outcome in each bin and bin 1 before and after inverse propensity score weighting. As the figure shows, IPW significantly reduces pretreatment differences in the outcome variable, making the bins more comparable to each other and strengthening the credibility of conditional ignorability.
We plot the main results in Figure (ref), whose panels show, from top to bottom, the counterfactual mean of the post-treatment outcome in each treatment bin; the estimated treatment effect associated with incrementing each bin; and lastly the proportion of the dataset in each bin, which is clearly concentrated in the zero-usage bin.
As Figure (ref) shows, AIPW estimates a large positive treatment effect of moving from the zero-usage bin (bin 1) to the next bin (bin 2). Moreover, because Feature A usage is zero-inflated, bin 1 is the most representative bin. Therefore, as shown in Table (ref), the coarsened AIPW estimate is positive, statistically significant, and substantially larger in magnitude than the RORR estimate. This discrepancy arises because the coarsened AIPW estimator explicitly weights the treatment effects to be representative of the treatment distribution, whereas RORR up-weights units with higher values of the treatment, where the dose-response curve is downward-sloping. The discrepancy is highly relevant for decision making: Although RORR indicates that Feature A has a negative treatment effect on the outcome, the AIPW results show that increasing Feature A usage would have a positive effect for the vast majority of members. Indeed, all nonzero Feature A usage bins have a higher conditional means than the zero-usage bin, indicating that any Feature A usage is preferable to none.
Although DML estimators are becoming increasingly popular in both academic and commercial research, researchers must, as ever, carefully evaluate their suitability for specific applications. Focusing on the residuals-on-residuals regression (RORR), this paper studies the interpretation of RORR when treatment effects are heterogeneous. We show that, for many-valued treatments, RORR converges to a conditional variance-weighted average of causal derivatives, with the added complication that these derivatives are evaluated on a “pseudo-treatment” distribution that differs from the treatment distribution seen in the data. As our empirical application shows, the subtle biases of RORR relative to the average treatment effect can have significant consequences for decision-making. To address these biases, we propose a coarsened AIPW estimator and demonstrate that this yields more representative estimates of causal effects.
For valuable suggestions, we thank Peter Hull, Yi Zhang, and participants in the Workshop on Causal Inference and Machine Learning in Practice at KDD '25 in Toronto, CA. For building and maintaining the observational causal inference platform that utilizes the methods described in this paper, we thank Adrien Alexandre, Colin Gray, and Dan Zylberglejd.
\balance