Shapley Value Confidence Intervals for Attributing Variance Explained

The coefficient of determination, the R 2 , is often used to measure the variance explained by an affine combination of multiple explanatory covariates. An attribution of this explanatory contribution to each of the individual covariates is often sought in order to draw inference regarding the importance of each covariate with respect to the response phenomenon. A recent method for ascertaining such an attribution is via the game theoretic Shapley value decomposition of the coefficient of determination. Such a decomposition has the desirable efficiency, monotonicity, and equal treatment properties. Under a weak assumption that the joint distribution is pseudo-elliptical, we obtain the asymptotic normality of the Shapley values. We then utilize this result in order to construct confidence intervals and hypothesis tests for Shapley values. Monte Carlo studies regarding our results are provided. We found that our asymptotic confidence intervals required less computational time to competing bootstrap methods and are able to exhibit improved coverage, especially on small samples. In an expository application to Australian real estate price modeling, we employ Shapley value confidence intervals to identify significant differences between the explanatory contributions of covariates, between models, which otherwise share approximately the same R 2 value. These different models are based on real estate data from the same periods in 2019 and 2020, the latter covering the early stages of the arrival of the novel coronavirus, COVID-19.


INTRODUCTION
When conducting statistical estimation and computation, the assumption of randomness of data necessitates that we address not only the problem of point estimation, but also variability quantification. In Ref. 1, variability for the coefficient of determination Shapley values were quantified via the use of bootstrap confidence intervals (CIs). Combined with the usual computational intensiveness of bootstrap resampling (see, e.g., Refs. 2 and 3; Ch. 12)), the combinatory nature of the computation of Eq. 5 (notice that |P| d!) compounds the time complexity of such a method, which is already of order O(2 d ). In this article, we seek to provide an asymptotic method for computing CIs for the Shapley values.
Our approach uses the joint asymptotic normality result of the elements in a correlation matrix, under an elliptical assumption, via [4], combined with asymptotic normality results concerning the determinants of a correlation matrix, of Refs. 5 and 6. Using these results, we derive the asymptotic joint distribution for the R 2 (Z n ) Shapley values, which allows us to construct CIs for each of the values and their contrasts. We assess the finite sample properties of our constructions via a comprehensive Monte Carlo study and demonstrate the use of our CIs via applications to real estate price data.
The remainder of the article proceeds as follows. In Section 3, we present our main results regarding the asymptotic distribution of the coefficient of determination Shapley values, and their CI constructions. In Section 4, we present a comprehensive Monte Carlo study of our CI construction method. In Section 5, we demonstrate how our results can be applied to real estate price data. Conclusions are lastly drawn in Section 6. All data and scripts to reproduce the analyses are available at https://github. com/ex2o/shapley_confidence.

The Coefficient of Determination
The multiple linear regression model (MLR) is among the most commonly applied tools for statistics inference; see Refs. 7 and 8; Part I) for thorough introductions to MLR models. In the usual MLR setting, one observes an independent and identically distributed (IID) sample of data pairs Z u , . . . , n}, and d, n ∈ N. The MLR model is then defined via the linear relationship where β u (β 0 , . . . , β n ) ∈ R d+1 , and X u i (X i1 , . . . , X id ). We shall also write Z u i (Z i0 , Z i1 , . . . , Z id ), when it is convenient to do so. Here, (·) u is the transposition operator.
The usual nomenclature is to call the Y i and X i elements of each pair, the response (or dependent) variable and the explanatory (or covariate) vector, respectively. Here, the jth element of X i : X ij (j ∈ [d]), is referred to as the jth explanatory variable (or the jth covariate). We may put the covariate and response pairs into a dataset Z n {Z i } n i 1 . Let R jk (Z n ) denote the sample correlation coefficient for each j, k ∈ {0}∪ [d]. Here Z j n −1 n i 1 Z ij is the sample mean of variable j. Write U4{0} ∪ [d] be a nonempty subset, where U {u 1 , . . . , u |U| }, where |U| is the cardinality of U. We refer to the matrix of correlations between the variables in U as A common inferential task is to determine the degree to which the response can be explained by the covariate vector, in totality. The usual device for addressing this question is via the coefficient of determination (or squared coefficient of multiple correlation), which is defined as where |C| is the matrix determinant of C (cf.

The Shapley Value
A refinement to the question that is addressed by the R 2 (Z n ) coefficient, is that of eliciting the contribution of each of the covariates to the total value of R 2 (Z n ).
In the past, this question has partially been resolved via the use of partial correlation coefficients (see, e.g., Ref. 8 The Shapley decomposition is a game-theoretic method for decomposing the contribution to the value of a utility function in the context of cooperative games. Let π u (π 1 , . . . , π d ) be a permutation of the set be the elements of [d] that appear before π j when [d] is permuted by π. We may define R 2 S j (π) (Z n ) and R 2 {j}∪ S j (π) (Z n ) in a similar manner to Eq. 3, using the generic definition for nonempty subsets S4[d], and R 2 { } (Z n ) 0 for the empty set. Treating the coefficient of determination as a utility function, we may conduct a Shapley partition of the R 2 (Z n ) coefficient by computing the jth Shapley value, for each of the d covariates, defined by where P is the set of all possible permutations of [d].

Uniqueness of the Shapley Value
Compared to other decompositions of the coefficient of determination, such as those considered in Refs. 13 and 14, the Shapley values, obtained from the partitioning above, have the favorable axiomatic properties that were well exposed in Ref.
Monotonicity: For pairs of samples Z m and Z n , of sizes m, n ∈ N, Equal treatment: If covariates j, k ∈ [d] are substitutes in the sense that for each π ∈ P such that k ∉ S j (π) and j ∉ S k (π), then the Shapley decomposition is such that V j (Z n ) V k (Z n ).
We note that equal treatment is also often referred to as symmetry in the literature. The uniqueness of the Shapley decomposition in exhibiting the three described properties is often used as the justification for its application. Furthermore, there are numerous sets of axiomatic properties that lead to the Shapley value decomposition as a solution (see, e.g., Ref. 16

(6)
Remark 1. We note that the elliptical distribution assumption above can be replaced by a broader pseudo-elliptical assumption, as per Refs. 21 and 22. This is a wide class of distributions that includes some that may not be symmetric. Due to the complicated construction of the class, we refer the interested reader to the source material for its definition.
Remark 2. We may state a similar result that replaces the elliptical assumption by a fourth moments existence assumption instead, using Proposition 2 of Ref. 4. In order to make practical use of such an assumption, we require the estimation of (d + 1)!/[(d − 3)!4!] fourth order moments instead of a single kurtosis term κ. Such a result may be useful when the number of fourth order moments is small, but become infeasible rapidly, as d increases.
Let ν4{0}∪[d], where ν {v 1 , . . . , v |ν| }. Define C n (ν) in the same manner as Eq. 2, and let The following theorem is adapted from a result of Ref. 5 (also appearing as Theorem 1 in Ref. 23). Our result expands upon the original theorem, to allow for inference regarding elliptically distributed data, and not just normally distributed data. We further fix some typographical matters that appear in both Refs. 5 and 23.
where acov(ζ gh , ζ jk ) is as per Eq. 6, and r ν (j, k) (j, k ∈ ν) is defined similarly. Proof. The result is due to an application of the delta method (see, e.g., Ref. 24; Thm. 3.1)) and the fact that for any matrix R, the derivative of its determinant is z|R|/zR |R|R −u [25]; Sec. 17.45). Notice that we use the unconstrained case of the determinant derivative, since we sum over each pair of [25]; Sec. 17.45). Using this fact, we may write Eq. 7 in the alternative, and more computationally efficient form

The Coefficient of Determination
Let plim denote convergence in probability, so that for any sequence {X n }, and any random variable X, the statement for every ε > 0, can be written as plim n → ∞ X n X. Now, recall definition Eq. 4, and further let We adapt from and expand upon [23]; Thm. 2) in the following result. This result also fixes typographical errors that appear in the original theorem, as well as in Ref. 5.
Lemma 3. Assume the same conditions as in Eq. 1. Then, the normalize coefficient of determination λ(S) n 1/2 (R 2 S (Z n ) − ρ 2 S ) (where S and T are nonempty subsets of [d]) converges to a jointly normal distribution, with asymptotic mean and covariance elements 0 and Proof. We apply the delta method again, using the functional form Eq. 4, and using the fact that Remark 4. When S T , Eq. 8 yields the usual form for the asymptotic variance of R 2

A Shapley Value Confidence Interval
where Using this functional form Eq. 9, we may apply the delta method once more, in order to derive the following joint asymptotic normal distribution result regarding the Shapley values Remark 5. The form Eq. 9 is a useful computational trick that reduces the computational time of form Eq. 5 and results in more efficient computations for fixed d. It is unclear whether other formulations such as that of Ref. 26 converge to a jointly normal distribution, with asymptotic mean and covariance elements 0 and Using the result above, we may apply the delta method again in order to construct asymptotic CIs or hypothesis tests regarding any continuous function of the d Shapley values for the coefficient of determination. Of particular interest is the asymptotic CI for each of the individual Shapley values and the hypothesis test for the difference between two Shapley values.
The asymptotic 100(1 − α)% CI for the jth expected Shapley value v j has the usual form where avar(ξ j ) acov(ξ j , ξ j ) denotes the asymptotic variance of ξ j and Φ −1 is the inverse cumulative distribution function of the standard normal distribution. The z-statistic for the test of the null and alternative hypotheses where Δ n has an asymptotic standard normal distribution. Remark 6. In practice, we do not know the necessary elements κ and ρ jk (j, k ∈ {0}∪ [d]; j ≠ k) that are required in order to specify the asymptotic covariance terms in Eqs 1-3 and Eq. 1. However, by Slutsky's theorem, we have the usual result that any acov (or avar) term can be replaced by the estimator acov n (or avar n ), which replaces κ by the estimator of Ref. 19 where Σ jk (Z n ) is the sample covariance between the jth and kth dimension of Z (j, k ∈ {0}∪ [d]). For example, the estimated test statistic for the hypotheses Eq. 10, retains the property of having an asymptotically standard normal distribution. The confidence interval calculation is summarized in Algorithm 1.
Algorithm 1: A Shapley value confidence interval 1. Input: dataset Z n , significance level α, and feature index j 2. Compute the estimate V j (Z n ) of the Shapley value v j , via Eq. 9 3. Compute a jj , b jj , c jj , d jj from Theorem 1, taking care to replace ρ jk in Lemmas 1-3 by the estimate R jk (Z n ) (see Remark 6) 4. Compute avar n (ξ j ) a jj + b jj − c jj − d jj 5. Output: a 100(1 − α)% confidence interval for v j is then given by avar n ξ j n ⎞ ⎟ ⎟ ⎠

MONTE CARLO STUDIES AND BENCHMARKS
In each of the following three Monte Carlo studies, we simulate a large number N of random samples Z (i) n , i ∈ [N], of size n, from a chosen distribution D. For each sample, we apply Eq. 1 to calculate an asymptotic 95% CI for the first Shapley value v 1 , producing a set of N observed intervals Here, the population Shapley value v 1 has the form: where R 2 S is defined by replacing C n in Eq. 4 by the known population correlation matrix cor(Z), as determined by the chosen distribution D. In Studies A and B, this population correlation matrix is the (d + 1) × (d + 1) matrix Σ, with diagonal elements equal to 1 and off-diagonal elements equal to a constant correlation c ∈ [0, 1). That is, where J d+1 denotes a (d + 1) × (d + 1) matrix with all entries equal to 1. Note that, for fixed variance, larger magnitudes of c map to larger regression coefficients. Thus, these simulations can be viewed as assigning equal regression coefficients to each covariates that are increasing for increasing values of c.
In Study C, we are concerned with covariance matrices with off-diagonal elements deviating from c. We aim to capture a case where the off-diagonal elements of cor(Z) are not uniform, and where some may be negative. This is achieved by sampling Z (i) n from a multivariate normal distribution D i , with a symmetric positive definite covariance matrix Σ i that is sampled at random from a Wishart distribution with scale matrix Σ; see Section 4.3. Accordingly, the population Shapley value v (i) 1 is unique to sample i, and thus we adjust the coverage estimator covr v 1 ([L, U]) by replacing v 1 on the RHS of Eq. 12 by v (i) 1 . To accompany our estimates of covr v 1 ([L, U]), we also provide Clopper-Pearson CIs [27] for the coverage probability. We also report the average CI widths, and middle 95% percentile intervals for the widths. For comparison, we estimate coverage probabilities of non-parametric bootstrap confidence intervals in each of the three studies. To obtain the bootstrap CIs, we set some large number N b and take random resamples R (r) n , r ∈ [N b ], of size n, with replacement, from Z (i) n . From these resamples, we calculate the set of estimated Shapley values The ith 95% bootstrap CI is then taken as the middle 95% percentile interval of L, and the coverage is estimated as in Eq. 12.

MC Study A
Here, we choose D N d+1 (0, Σ), so that each sample Z (i) n , for i ∈ [N], is drawn from a multivariate normal distribution, with covariance Σ given in Eq. 13.
The simulation results in Figure 1 (and in Supplementary Figure S12) show very similar coverage and width performance between the two assessed CIs for moderate and high correlations c > 0.3. For lower correlations c ≤ 0.2, coverage convergence appears to be slower in n than the bootstrap CI for large sample sizes (n ≥ 100). The opposite trend seems to hold for small sample sizes (n ≤ 50), see the discussion under MC Study C. Also, for the highest correlation c 0.99, coverage performance of the asymptotic CI is overall slightly better than the bootstrap CI.
For all sample sizes n and correlations c, coverage and width performances are similar to MC Study A (see Supplementary  Figures S13 and S14). Of particular interest, in both MC Studies A and B (but not in MC Study C), we observe that for c 0, the estimated coverage probability of the asymptotic CI is almost equal to 1, for all sample sizes greater than 10, while the corresponding bootstrap CIs have estimated coverage equal to 0 (Figure 2 left). Despite this, the average CI widths, though large under small samples, are somewhat smaller than those for bootstrap (Figure 2

MC Study C
Here, we set D i N d+1 (0, Σ i ), so that the sample Z (i) n , for i ∈ N, is drawn from a multivariate normal distribution, with covariance matrix Σ i realized from a Wishart distribution W d+1 (Σ, ]) with scale matrix Σ and ν degrees of freedom. This set up is different from Studies A and B in that the distributions D i , and therefore the population Shapley values v (i) 1 , are allowed to differ between samples.
The distribution W d+1 (Σ, ]) can be understood as the distribution of the sample covariance matrix of a sample of size ] + 1 from the distribution N d+1 (0, Σ) (cf. Ref. 28). This implies that each covariance matrix Σ i can have non-uniform and negative off-diagonal elements, with variability between the off-diagonal elements increasing as ν decreases. For this study, we set ] 100.
Aside for the case c 0, coverage and width statistics are again similar to MC Studies A and B, for all n and c (see Supplementary  Material). Interestingly, in all three Studies, for small sample sizes (n ≤ 50), coverage is often higher than for the bootstrap CI, with slightly smaller average widths, as seen in Figure 3 (and in Supplementary Figures S12, S13). For the c 0 case, the observed behavior differs from MC Studies A and B, with bootstrap performing comparatively well for large sample sizes (Figures 4, 5).

Computational Benchmarks
From Figure 6, we see that the memory usage (left) and mean execution time (right) for the bootstrap CIs are both higher than that for the asymptotic CIs, and that the ratio increases with sample size given the benchmarking parameters in Table 1. As n increases, asymptotic CIs become increasingly efficient, compared to the bootstrap CIs. On the other hand, as d increases, with n fixed, we expect an increase in the relative efficiency of the bootstrap, since the complexity of calculating acov(ξ j , ξ k ) in Eq. 1 grows faster in d than the complexity of the bootstrap procedure.

Summary of Results and Recommendations for Use
The most relevant observation is that the widths all appear to shrink at a rate that is comparable to the expected n √ convergence dictated by asymptotic normality. The remaining observations below are focused on computational time, sample size, and correlation, in comparison with bootstrap intervals. Below follows a summary of the general tendencies and observations from our results, and recommendations regarding when to use the asymptotic CIs.
• For all correlations c in all three Studies, the estimated coverage probability of asymptotic intervals is above 0.85 for all sample sizes n ≥ 10. • For small correlations, in particular c ≤ 0.1 and sample sizes 10 ≤ n ≤ 100, the lower bound of the confidence interval for the coverage of the asymptotic CIs never drops below 0.91. • For c 0 and n ≥ 15, the lower bound of the asymptotic CI for coverage never drops below 0.95 in Studies A and B, while in Study C the lower bound is at least 0.88. • For sample sizes 15 ≤ n ≤ 50, the coverage of the asymptotic CI tends to be higher when c is closer to the boundaries of [0, 1], as shown in Figure 7. • The average asymptotic CI width is lower when c is nearer to the boundaries of [0, 1], see Figure 8.
We now make some general observations which apply to all three Studies. As sample size increases, the estimated coverage initially increases rapidly, as can be seen, for example, in the left column of Figure 3. For small sample sizes between n 5 and n 50, the asymptotic CIs typically outperform bootstrap CIs,    especially when c lies farther from 0.5; there is a clear drop in coverage as c approaches 0.5, for small samples, as can be seen in Figure 7. In many cases, the estimated coverage is above 0.9 for n ≥ 10. However, empirical coverage does not appear to be an increasing function of sample size in general. On the top row in the left column of Figure 1, we observe one example of a clear and extended dip in coverage for n in [100, 1000]. This gives rise to the general observation that the asymptotic intervals have preferable coverage statistics over bootstrap for small samples, but not for a certain range of large samples, depending on c. We further observe that, for all n, there is a general increase in the average CI width as c approaches 0.5 from either direction, as in Figure 8. In all Studies, over all sample sizes and correlations, the bootstrap CI average widths were smaller than the asymptotic CI widths by at most 0.0289, and vice versa by at most 0.0667. In general, the asymptotic intervals display favourable widths, though less so near c 0.5.

Recommendations for Use
Based on these observations, we recommend using asymptotic CIs over bootstrap CIs under the following conditions: (i) Computational time is relevant (e.g., estimating a large number of Shapley values). (ii) The sample size is small (e.g., n ≤ 50). (iii) The correlation between explanatory variables and the response variable is expected to be beyond ± 0.2 from 0.5, or when this is where the highest precision is desired.
We note that our observations are made from an incomplete albeit comprehensive set of simulation scenarios. There are of course an infinite number of combinations of simulation cases and thus we cannot guarantee that our observation applies to all possible DGPs.  Frontiers in Applied Mathematics and Statistics | www.frontiersin.org December 2020 | Volume 6 | Article 587199

Application: Real Estate and COVID-19
For an interesting application of our methods, in this section we identify significant changes in the behavior of the local real estate market in Melbourne, Australia, between the period 1 February 2019 to 1 April 2019 and the period 1 February 2020 to 1 April 2020. In 2020, this corresponds to an early period of growing public concern regarding the novel coronavirus COVID-19. We obtain the Shapley decomposition of the coefficient of multiple correlation R 2 between observed house prices and a number of property features. We also find significant differences in behavior between real estate near and far from the Central Business District (CBD), where near is defined to be within 25 km of Melbourne CBD, and far is defined as non-near (see Figure 9).
Note that the nature of this investigation is exploratory and expository; our conclusions are not intended to be taken as evidence for the purposes of policy or decision making. On 1 February the Australian government announced a temporary ban on foreign arrivals from mainland China, and by 1 April a number of social distancing measures were in place. We scraped real estate data from the AUHousePrices website (https://www.auhouseprices.com/), to obtain a data set of 13, 291 (clean) house sales between 1 January and 18 July in each of 2019 and 2020. We then reduced this date range to capture only the spike in sales observed between 1 February and 1 April (see Figure 10), giving a remaining sample size of 5,110, which was partitioned into the four subgroups in Table 2. By capturing this spike in sales, we hope to detect a significant change in the Shapley values between the same period in 2019 and 2020. Within each of the four subgroups we perform a Yeo-Johnson transformation [29] to reduce any violation of the assumption of joint pseudo-ellipticity. We decompose R 2 among the covariates: distance to CBD (CBD); images used in advertisement (images); property land size (land); distance to nearest school (school); distance to nearest station (station); and number of bedrooms + bathrooms + car ports (room); along with the response variable, house sale price (price). We expect the room covariate to act as a proxy for house size. Thus we decompose R 2 for the linear model, price β 0 + β CBD, images, land, school, station, room T , Fitted to each of the four subgroups, we obtain R 2 0.37 for model Eq. 14, for each subgroup except for the near (2020) subgroup, for which R 2 0.39. The resulting Shapley values and 95% confidence intervals are listed in Table 3, and shown graphically in Figure 11. Also in Figure 11, the corresponding bootstrap CIs are given for completeness, and are almost identical to the asymptotic CIs. From these results we make the following observations regarding attributions of the total variability in house prices explained by model Eq. 14: (i) In both 2019 and 2020, the attribution for distance to CBD was significantly higher for house sales near to the CBD, compared to house sales farther from the CBD. Correspondingly, the attribution for roominess was significantly lower for house sales near to the CBD. (ii) Among sales that were near to the CBD, distances to the CBD received significantly greater attribution than both land size and roominess in 2019. Unlike the case for sales that were far from the CBD, these differences remained significant in 2020. (iii) Distances to stations and schools, as well as images used in advertising, have apparently had little overall impact. In all four subgroups, the attribution for distance to the nearest school is not significantly different from 0. However, distance to a station does receive significantly more attribution among houses that are near to the CBD, compared to those farther away. (iv) Interestingly, while not a significant result, the number of images used in advertising did appear to receive greater attribution among house sales that were far from the CBD, compared to those near to it. (v) Among sales that were far from the CBD, land size and roominess both received significantly more attribution than distance to the CBD, in 2019. However, this difference vanished in 2020, with distances to the CBD apparently gaining more attribution, while roominess and land size apparently lost some attribution, in a relative leveling out of these three covariates. Item (i) is perhaps unsurprizing: distances were less relevant far from the city, where price variability was influenced more by roominess and land size. Indeed, we can assume we are less likely to find small and expensive houses far from the CBD. However, the authors find Item (v) interesting: near the city, the behavior did not change significantly during the 2020 period. However, far from the city, the behavior did change significantly, moving toward the near-city behavior. Distance to the city became more important for explaining variability in price, while land size and roominess both became less important, compared with their 2019 estimates. Our initial guess at an explanation was that near-city buyers, with nearcity preferences, were temporarily motivated by urgency to buy farther out. However, according to Table 2, the observed ratio of near-city buyers to non-near buyers actually increased in this period, from 1.26 in 2019 to 1.61 in 2020. We will not take this expository analysis further here, but we hope that the interested reader is motivated to take it further, and to this end we have made the data, as well as R and Julia scripts, available at https:// github.com/ex2o/shapley_confidence.

DISCUSSION
In this work, we focus on regression Shapley values as a means for attributing goodness of fit. There has also been much recent interest in the machine learning literature in using Shapley values to attribute predictions. These approaches, such as SHAP [30], TreeExplainer [31], and Quantitative Input Influence [32] are focused on providing local explanations that depend on the behavior of a fitted model. For more on the differences between prediction, estimation and attribution, see the recent exposition of Ref. 33.
Research into the attribution of goodness of fit using the Shapley value, predates research into attributing predictions. As described in Section 2, regression Shapley values were first introduced in Ref. 10 Our work is the first to determine the asymptotic distribution of the regression Shapley values. In Section 3, we show that under an elliptical (or pseudo-elliptical) joint distribution assumption, the Shapley value decomposition of R 2 (Z n ) is asymptotically normal. This involved correcting a number of errors that appeared in the expositions of Lemmas 2 and 3 in prior sources, and stating those results in a notation conforming to the language of cooperative games and set-valued functions. Our result regarding the asymptotic distribution leads immediately to Shapley value confidence intervals and hypothesis tests.
In Section 4, we use Monte Carlo simulations to examine the coverage and width statistics of these asymptotic CIs. We compare the coverage and width statistics for different data generating processes, correlations and sample sizes. In Section 4.5 we provide recommendations for when asymptotic CIs should be used rather than bootstrap CIs. We find that the asymptotic CIs have estimated coverage probabilities of at least 0.85 across all studies, are preferable over the bootstrap CIs for small sample sizes (n ≤ 50), and are often (although not always) favourable for large sample sizes. We discuss the generality of our findings that asymptotic CIs are computationally much more efficient than bootstrap CIs.
In Section 5, we derive asymptotic CIs to Melbourne house prices during a period of altered consumer behavior in the initial stages of the arrival of COVID-19 in Australia. Using the CIs, we identify significant changes in model behavior between 2019 and 2020, and attribute these changes among features, highlighting an interesting direction of future research into the period. Implementations of our methods as well as data sets are openly available for use at https://github.com/ex2o/shapley_ confidence.
We aim to use our developments to derive the asymptotic distributions of variance inflation factors and their generalizations [34], as well as the closely related Owen values decomposition of the coefficient of determination for exogenously grouped sets of features [1]. We will also investigate calculating Shapley values capable of uncovering non-linear dependencies and their confidence intervals in future work.

AUTHOR'S NOTE
This manuscript has been released as a pre-print on arXiv [35].

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DF Provided Monte Carlo simulations, theoretical contributions, real estate data scraping and analysis, figures, computational benchmarks and general writing. IS Provided real estate data analysis, theoretical contributions, general writing and proofreading. HN Provided theoretical proofs, theoretical contributions, general writing and proofreading.