# Systematic Review of Newton-Schulz Iterations with Unified Factorizations :
Integration in the Richardson Method and Application to Robust Failure
Detection in Electrical Networks

Alexander Stotsky
Department of Computer Science and Engineering
Chalmers University of Technology
Gothenburg SE - 412 96, Sweden
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Systematic overview of Newton-Schulz and Durand iterations with convergence
analysis and factorizations is presented in the chronological sequence in
unified framework. Practical recommendations for the choice of the order and
factorizations of the algorithms and integration into Richardson iteration are
given. The simplest combination of Newton-Schulz and Richardson iteration is
applied to the parameter estimation problem associated with the failure
detection via evaluation of the frequency content of the signals in electrical
network. The detection is performed on real data for which the software
failure was simulated, which resulted in the rank deficient information
matrix. Robust preconditioning for rank deficient matrices is proposed and the
efficiency of the approach is demonstrated by simulations via comparison with
standard LU decomposition method.

##  1 Introduction

Parameter estimation problems appear often in adaptive control, [1],[2] system
identification,[3], signal processing, [3] and in many other areas. Suppose
that the process is described by the following equation:

$y_{k}=\varphi_{k}^{T}\theta_{*}$ (1)

where $\varphi_{k}$ is the regressor vector and $\theta_{*}$ is the vector of
unknown parameters to be estimated, $k=1,2,...$
Introduction of the model in the form

$\hat{y}_{k}=\varphi_{k}^{T}\theta_{k}$ (2)

and minimization of the following function associated with the mismatch
between the model and the process

$E_{k}=\sum_{p=1}^{p=k}(\hat{y}_{p}-y_{p})^{2}$ (3)

yields the following restricted linear system of equations :

$\displaystyle A_{k}\theta_{k}$ $\displaystyle=b_{k}$ (4) $\displaystyle
A_{k}$ $\displaystyle=\sum_{p=1}^{p=k}\varphi_{p}\leavevmode\nobreak\
\varphi_{p}^{T},\leavevmode\nobreak\
b_{k}=\sum_{p=1}^{p=k}\varphi_{p}\leavevmode\nobreak\ y_{p}$ (5)

which should be solved with respect to $\theta_{k}$. Accurate and
computationally efficient solution of (4) is the most significant part of the
estimation problem. The solution may benefit from the properties of the matrix
$A_{k}$, which depend on the type of regressor and the step number associated
with the window size. The matrix $A_{k}$ is SPD (Symmetric and Positive
Definite) matrix in many cases. For example, the matrix $A$ is SPD for the
systems with harmonic regressor, [1] \- [5]. In general case the equation (4)
can be multiplied by $A^{T}_{k}$ (for any invertible matrix $A_{k}$) that
transforms the system to the SPD case with the Gram matrix, [6]. Notice that
the Gram matrix is often ill-conditioned matrix due to squaring of the
condition number [6].
The matrix $A$ is SPD ill-conditioned (a) information matrix for systems with
harmonic regressor and short window sizes, [7],[8] (b) Gram matrix due to the
squaring of the condition number in the least squares method, [3], [6], (c)
mass matrix in finite element method, [9] and mass matrix (lumped mass matrix)
for mechanical systems with singular perturbations, [10], (d) state matrix for
systems of linear equations, [11] and in many other applications. Moreover,
many processes in practice are singulary perturbed (have modes with different
time-scales), [12] and considered as stiff and ill-conditioned systems, which
potentially extends application areas of this model. In other words the model
(4) with SPD (and possibly ill-conditioned matrix $A_{k}$) covers almost all
the cases in the parameter estimation applications.
Ill-conditioning implies robustness problems (sensitivity to numerical
calculations) and imposes additional requirements on the accuracy of the
solution of (4) (which are especially pronounced in finite-digit calculations)
since small changes in $b_{k}$ due to measurement, truncation, accumulation,
rounding and other errors result in significant changes in $\theta_{k}$. In
addition, ill-conditioning implies slow convergence of the iterative
procedures.
The accuracy requirements is the main motivation for application of the
Richardson iteration which is driven by the residual error and has filtering
(averaging) property. The residual error is smoothed and remains bounded for a
sufficiently large step number in this iteration, where the bound depends on
inaccuracies, providing the best possible solution in finite digit
calculations.
Therefore convergence rate improvement and reduction of the computational
complexity of the Richardson iteration is one of the most important problems
in the area. One of the ways to convergence rate improvement is introduction
of the Newton-Schulz and Durand algorithms in the Richardson iteration, [8],
[13] \- [16].
This paper starts with short systematic overview (presented in unified
tutorial fashion) of Newton-Schulz and Durand iterations in Section 2. Special
attention is paid to practical introduction to power series factorizations for
reduction of the computational complexity in Section 3. The simplest
integration of Newton-Schulz algorithms in the Richardson iteration is
presented in Section 4 with application to practical problem in Section 5. The
problem is associated with the failure detection via estimation of the
frequency content of the signals in electrical system in the presence of
missing data. The detection is performed on real signals with simulated
failure in the software. The efficiency of the approach is demonstrated in
Section 5. The paper ends with brief conclusions in Section 6.

##  2 Overview of Iterative Matrix Inversion Algorithms: Mini Tutorial

Brief overview of iterative matrix inversion algorithms with convergence
analysis is presented in this Section. The overview is presented in the
systematic way (from simple algorithms to more sophisticated ones) in the
chronological sequence with convergence analysis and starts with the second
order Newton-Schulz algorithms in Section 2.1 followed by the description of
high order Newton-Schulz algorithms in Section 2.2. Computationally efficient
Durand iterations are described in Section 2.3, and finally the combination
for convergence rate improvement of high order Newton-Schulz and Durand
algorithms is presented in Section 2.4.

###  2.1 Second Order Newton-Schulz Iterations

Consider the following iterations, [17] \- [20]:

$G_{k}=G_{k-1}+F_{k-1}G_{k-1},\leavevmode\nobreak\ \leavevmode\nobreak\
F_{k}=I-G_{k}A,$ (6)

where $I$ is the identity matrix, $k=0,1,2...$ . Evaluation of the error in
the first step $F_{1}$ yields :

$\displaystyle F_{1}$ $\displaystyle=I-(G_{0}+F_{0}\leavevmode\nobreak\
G_{0})\leavevmode\nobreak\ A=I-(I+F_{0})\leavevmode\nobreak\
G_{0}\leavevmode\nobreak\ A$ $\displaystyle=I-(I+F_{0})\leavevmode\nobreak\
(I-F_{0})=I-I^{2}+F_{0}^{2}=F_{0}^{2}$ $\displaystyle F_{2}$
$\displaystyle=F_{1}^{2}=F_{0}^{4}$

which gives the error model:

$F_{k}=F_{0}^{2^{k}}$ (7)

Algorithm (6) which is also known as Newton-Schulz-Hotelling matrix iteration,
[18] estimates the inverse of the matrix $A$ via minimization of the
estimation error $F_{k}$ provided that the spectral radius of initial matrix

$\rho(F_{0})=\rho(I-G_{0}A)<1$ (8)

is less than one, where $G_{0}$ is preconditioner (initial guess of $A^{-1}$).
The proof of convergence is presented in [17] \- [20], where the algorithm (6)
is written in the following form: $G_{k}=G_{k-1}+G_{k-1}F_{k-1}$, where
$F_{k}=I-AG_{k}$.
The order of the iterative algorithm can be increased for convergence rate
improvement by introduction of the power series in the algorithm (6).

###  2.2 High Order Newton-Schulz Iterations

High order algorithms (which are also known as hyperpower iterative
algorithms, which are widely used for calculation of generalized inverses) can
be presented in the following form for $n=2,3,...$, [21] \- [23] :

$\displaystyle G_{k}$
$\displaystyle=\sum_{j=0}^{n-1}F^{j}_{k-1}G_{k-1}=(I+F_{k-1}+F^{2}_{k-1}+...+F^{n-1}_{k-1})\leavevmode\nobreak\
G_{k-1}$ (9)

or in the form presented in [16], [22] :

$\displaystyle G_{k}$
$\displaystyle=\sum_{j=0}^{n-1}F^{j}_{k-1}G_{k-1}=G_{k-1}+F_{k-1}\leavevmode\nobreak\
\sum_{j=0}^{n-2}F^{j}_{k-1}G_{k-1}$ (10)

The error $F_{k}$ in (9) can be written in the following form:

$\displaystyle F_{k}$
$\displaystyle=I-G_{k}A=I-(I+F_{k-1}+F^{2}_{k-1}+...+F^{n-1}_{k-1})G_{k-1}A$
$\displaystyle=I-(I+F_{k-1}+F^{2}_{k-1}+...+F^{n-1}_{k-1})(I-F_{k-1})=F_{k-1}^{n}$
$\displaystyle F_{k}$ $\displaystyle=F_{0}^{n^{k}}$ (11)

Notice that the algorithm (9) is written in following form in [21] \- [23]
(and it is widely used in this form):

$\displaystyle G_{k}$
$\displaystyle=G_{k-1}\sum_{j=0}^{n-1}F^{j}_{k-1}=G_{k-1}(I+F_{k-1}+F^{2}_{k-1}+...+F^{n-1}_{k-1})$
(12)

with $F_{k}=I-AG_{k}$. Notice that the error $F_{k}=I-G_{k}A$ is better
aligned with Richardson iteration, see Section 4 and therefore is considered
further in this paper.

###  2.3 Durand Iterations

Computationally efficient matrix inversion method described by Durand, [24],
[25] can be derived from the relation associated with the initial error in
Newton-Schulz approach, $A^{-1}=F_{0}A^{-1}+G_{0}$ by substitution of $G_{k}$
instead of $A^{-1}$ as follows :

$\displaystyle G_{k}$ $\displaystyle=F_{0}\leavevmode\nobreak\ G_{k-1}+G_{0}$
(13) $\displaystyle G_{k}-A^{-1}$ $\displaystyle=F_{0}\leavevmode\nobreak\
(G_{k-1}-A^{-1})$ (14) $\displaystyle F_{k}$
$\displaystyle=I-G_{k}A=F^{k+1}_{0}$ (15)

The error models (14), (15) are valid for the algorithm (13). Higher order
Durand iterations are discussed in [8], [25].
Notice that the convergence rate of Newton-Schulz Iterations is significantly
higher than the convergence rate of Durand iterations. However, Durand
iteration (13) requires only one matrix product in each step.
Notice also that the convergence performance can be enhanced by introduction
of the initial power series as preconditioning $G_{0}$ (which is calculated
only once) in Newton-Schulz and Durand iterations, [8], [26]. The order of
initial power series can be chosen sufficiently high, which essentially
improves the convergence rate. Computational complexity of initial power
series can be reduced via factorizations, see Section 3 for details.
Convergence rate can also be improved via combination of two Newton-Schulz
iterations, which also cover Durand iterations as special case, see next
Section 2.4.

###  2.4 Convergence Rate Improvement via Combination of High Order Newton-
Schulz Iterations

Consider the following combination of two Newton-Schulz iterations associated
with two-step iterative method, [27] :

$\displaystyle\Gamma_{k}$ $\displaystyle=I-L_{k-1}A$ (16) $\displaystyle
L_{k}$ $\displaystyle=\\{\sum_{j=0}^{n-1}\Gamma^{j}_{k}\\}\leavevmode\nobreak\
L_{k-1}$ (17) $\displaystyle\Gamma^{n}_{k}$ $\displaystyle=I-L_{k}A$ (18)
$\displaystyle G_{k}$
$\displaystyle=\underbrace{L_{k}}_{\begin{subarray}{l}\text{Newton-Schulz }\\\
\text{Iteration
}\end{subarray}}+\underbrace{\Gamma^{n}_{k}}_{\begin{subarray}{l}\text{High
Order }\\\ \text{Convergence}\\\
\text{Accelerator}\end{subarray}}\leavevmode\nobreak\ \leavevmode\nobreak\
\underbrace{\\{\sum_{j=0}^{n-1}F^{j}_{k-1}\\}\leavevmode\nobreak\
G_{k-1}}_{\begin{subarray}{l}\text{Newton-Schulz }\\\ \text{Iteration
}\end{subarray}}$ (19) $\displaystyle F_{k}$
$\displaystyle=I-G_{k}A=\Gamma^{n}_{k}\leavevmode\nobreak\ F^{n}_{k-1}$ (20)
$\displaystyle F_{k}$ $\displaystyle=F_{0}^{\displaystyle\leavevmode\nobreak\
k\leavevmode\nobreak\ n^{k+1}+\leavevmode\nobreak\ n^{k}},\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ n>1$ (21)

where $\displaystyle L_{0}=\sum_{j=0}^{n-1}\Gamma^{j}_{0}\leavevmode\nobreak\
S^{-1}$, $\Gamma_{0}=F_{0}=I-S^{-1}A$, where $\displaystyle G_{0}=S^{-1}$
satisfies (8). Notice that Durand iteration (13) can be derived from (16) -
(19) with $n=1$. The error model (20) can be proved via multiplication of (19)
by $A$ and explicit evaluation of the error $F_{k}$ (similar to Section 2.2,
see also [21], [27]) as follows:

$\displaystyle F_{k}$
$\displaystyle=I-G_{k}A=\underbrace{I-L_{k}A}_{\Gamma^{n}_{k}}-\Gamma^{n}_{k}\sum_{j=0}^{n-1}F^{j}_{k-1}\underbrace{G_{k-1}A}_{I-F_{k-1}}=\Gamma^{n}_{k}\leavevmode\nobreak\
F^{n}_{k-1}$ (22)

The error model (21) is obtained via explicit evaluation of (22). The
algorithm (16) - (21) has two independent and equally complex computational
parts: the first part is associated with calculations of $\Gamma_{k}$, $L_{k}$
and $\Gamma^{n}_{k}$, where $\Gamma_{k}=\Gamma^{n}_{k-1}$, whereas the term
$\displaystyle\sum_{j=0}^{n-1}F^{j}_{k-1}G_{k-1}$ is calculated in the second
part. This algorithm is ideally suited for parallel implementation.

Extended version of this survey is presented in Stotsky A., Recursive Versus
Nonrecursive Richardson Algorithms: Systematic Overview, Unified Frameworks
and Application to Electric Grid Power Quality Monitoring, Automatika, vol.
63, N2, pp. 328-337, 2022, https://doi.org/10.1080/00051144.2022.2039989,
[33].

Further developments can be found in Stotsky A.,Simultaneous Frequency and
Amplitude Estimation for Grid Quality Monitoring : New Partitioning with
Memory Based Newton-Schulz Corrections, IFAC PapersOnLine 55-9, pp. 42-47,
2022, https://www.sciencedirect.com/science/article/pii/S2405896322003937
,[34].

##  3 Reduction of Computational Complexity via Unified Factorizations

The following high order Newton-Schulz power series can be factorized for
composite orders $n$ as follows, [27] :

$\displaystyle G_{k}$
$\displaystyle=\\{\sum_{j=0}^{n-1}F_{k-1}^{j}\\}\leavevmode\nobreak\
G_{k-1}=\\{\sum_{j=0}^{w-1}F_{k-1}^{(p+1)j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{p}F_{k-1}^{d}\\}\leavevmode\nobreak\ G_{k-1}$ (23)
$\displaystyle F_{k-1}^{p+1}$
$\displaystyle=I-\\{\sum_{d=0}^{p}F_{k-1}^{d}\\}\leavevmode\nobreak\
G_{k-1}\leavevmode\nobreak\ A$
$\displaystyle=I-\\{\sum_{d=0}^{p}F_{k-1}^{d}\\}\leavevmode\nobreak\
\\{I-F_{k-1}\\},\leavevmode\nobreak\ F_{k-1}=I-G_{k-1}A$ (24) $\displaystyle
n$ $\displaystyle=w\leavevmode\nobreak\ (p+1),\leavevmode\nobreak\
\leavevmode\nobreak\ p=0,1,2,...,\leavevmode\nobreak\ w=1,2,3,...$ (25)

Taking into account that every prime number is followed after the composite
number the factorization for the prime orders $n+1$ can be performed as
follows:

$\displaystyle G_{k}$
$\displaystyle=(I+F_{k-1}\\{\sum_{j=0}^{w-1}F_{k-1}^{(p+1)j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{p}F_{k-1}^{d}\\})\leavevmode\nobreak\ G_{k-1}$ (26)
$\displaystyle=\\{\sum_{j=0}^{(n+1)-1}F_{k-1}^{j}\\}\leavevmode\nobreak\
G_{k-1}$

Similar factorizations can be found in [28]. The following efficiency index
introduced in [29] is widely used for quantification of the performance of
different factorizations:

$EI=\displaystyle n^{1/n_{p}}$ (27)

where $n$ is the order of the algorithm and $n_{p}$ is the number of matrix
products per each cycle. For illustration of the unified factorization
approach few examples with increasing EI are presented below.

###  3.1 Case of $n=8=w\leavevmode\nobreak\ (p+1)$, $p=3$, $w=2$

$\displaystyle G_{k}$
$\displaystyle=\\{\sum_{j=0}^{1}F_{k-1}^{4j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{3}F_{k-1}^{d}\\}\leavevmode\nobreak\ G_{k-1}$ (28)
$\displaystyle=(I+F_{k-1}^{4})(I+F_{k-1}^{2})(I+F_{k-1})\leavevmode\nobreak\
G_{k-1}$

which requires $6$ matrix products with $EI=1.4142$ .

###  3.2 Case of $n=9=w\leavevmode\nobreak\ (p+1)$, $p=2$, $w=3$

$\displaystyle G_{k}$
$\displaystyle=\\{\sum_{j=0}^{2}F_{k-1}^{3j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{2}F_{k-1}^{d}\\}\leavevmode\nobreak\ G_{k-1}$
$\displaystyle=(I+F_{k-1}^{3}+F_{k-1}^{6})(I+F_{k-1}+F_{k-1}^{2})\leavevmode\nobreak\
G_{k-1}$
$\displaystyle=(I+(I+F_{k-1}^{4})(I+F_{k-1}^{2})(F_{k-1}+F_{k-1}^{2}))\leavevmode\nobreak\
G_{k-1}$

which requires $6$ matrix products with $EI=1.4422$ .

###  3.3 Case of $n=10=w\leavevmode\nobreak\ (p+1)$, $p=1$, $w=5$

$\displaystyle G_{k}$
$\displaystyle=\\{\sum_{j=0}^{4}F_{k-1}^{2j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{1}F_{k-1}^{d}\\}\leavevmode\nobreak\ G_{k-1}$
$\displaystyle=(I+(F_{k-1}^{2}+F_{k-1}^{4})(I+F_{k-1}^{4}))(I+F_{k-1}))\leavevmode\nobreak\
G_{k-1}$

which requires $6$ matrix products with $EI=1.4678$ .

###  3.4 Case of $n=11=w\leavevmode\nobreak\ (p+1)+1$, $p=1$, $w=5$

$\displaystyle G_{k}$
$\displaystyle=(I+F_{k-1}\\{\sum_{j=0}^{4}F_{k-1}^{2j}\\}\leavevmode\nobreak\
\\{\sum_{d=0}^{1}F_{k-1}^{d}\\})\leavevmode\nobreak\ G_{k-1}$ (29)
$\displaystyle=(I+F_{k-1}(I+(F_{k-1}^{2}+F_{k-1}^{4})(I+F_{k-1}^{4}))(I+F_{k-1}))\leavevmode\nobreak\
G_{k-1}$
$\displaystyle=(I+(I+(F_{k-1}^{2}+F_{k-1}^{4})(I+F_{k-1}^{4}))(F_{k-1}+F^{2}_{k-1}))\leavevmode\nobreak\
G_{k-1}$

which requires $6$ matrix products, $EI=1.4913$, [30].
Notice that the idea of factorization is associated with Newton-Schulz method.
Therefore Newton-Schulz algorithms of low orders (order two or three) being
iterated for a number of steps can be applied instead of factorized Newton-
Schulz iterations of higher orders for the sake of robustness and efficiency.
Notice that three steps of the second order Newton-Schulz iteration are
equivalent to one step of eighth order iteration with $6$ matrix products, see
(28). However, application of the eleventh order algorithm, see (29) and [30]
requires also $6$ mmm and provides faster convergence than three steps of the
second order Newton-Schulz iteration. This algorithm has also the highest
efficiency index, $EI=1.4913$ among the algorithms presented above. Therefore
the proper choice of the order and factorization that is made for each
particular application should represent the trade-off between the robustness
and convergence rate.

##  4 Richardson Iterations

Many practical problems can be associated with restricted linear system of
equations:

$A\theta_{*}=b$ (30)

which should be solved for $\theta_{*}$. Richardson iteration is the most
promising solution to this problem which provides the best possible accuracy
in finite digit calculations. Despite the fact that the inverse of matrix $A$
is not required in order to find $\theta_{*}$ the approximate inverse
essentially improves convergence rate of the Richardson iteration, [8], [13]
\- [16]. The simplest combination of matrix inversion techniques and
Richardson algorithm can be presented as follows :

$\displaystyle\theta_{k}$
$\displaystyle=\theta_{k-1}-G_{k}\leavevmode\nobreak\
(A\theta_{k-1}-b)=F_{k-1}\theta_{k-1}+G_{k}b$ (31)

where $\theta_{k}$ is the estimate of $\theta_{*}$, $G_{k}$ is associated with
the estimate of the inverse of $A$ and $F_{k}$ is the inversion error, see
Section 2. Iteration (31) can be implemented using matrix vector products
only.
Notice that definition of the inversion error in the form $F_{k}=I-G_{k}A$ is
more convenient than the following definition $F_{k}=I-AG_{k}$ since the error
$F_{k}=I-G_{k}A$ is better aligned with Richardson iteration providing the
same convergence conditions. The convergence analysis of Richardson iteration
with the error $F_{k}=I-AG_{k}$ can be found for example in [14].
Any matrix inversion algorithms presented in Section 2 with factorizations
described in Section 3 can be applied in (31). Error models for different
choices of $G_{k}$ are presented in [8]. Since $G_{k}$ is getting closer to
$A^{-1}$ in each step the convergence rate of (31) is essentially improved.
Updating of $G_{k}$ can even be stopped after several steps for reduction of
the computational burden.

(a)

(b)

Figure 1:  The voltages which correspond to the sag and swell events in phases
’a’,’b’ and ’c’ are plotted with green, red and blue lines respectively in
Figure 1(a). The voltages are associated with the failure event with ID
$2911$. The data is obtained from the DOE/EPRI National Database of Power
System Events, [31]. The cause of this event is unknown. Amplitudes of the
first harmonic (with the corresponding colours) which show delectability of
the sag and swell events for the signals in Figure 1(a) are plotted in Figure
1(b). The detection is performed with rank deficient information matrix.

Figure 2:  The condition number of the matrix $A_{rk}$ in (40) as a function
of the regularization parameter $\beta>0$ is plotted. For a small $\beta$ the
matrix $A_{rk}$ is ill-conditioned matrix for which the spectral radius
$\displaystyle\rho(I-S_{k}^{-1}A_{rk})$, where $\displaystyle
S_{k}^{-1}=I/\alpha_{k}$ with $\alpha_{k}=\|A_{rk}\|_{\infty}/2+\varepsilon$
is very close to one.

##  5 Estimation of the Frequency Content of Electrical Signals for Robust
Failure Detection

Suppose that the measured electrical signal (voltage or current) $y_{k}$ can
be presented in the following form

$y_{k}=\varphi_{k}^{T}\vartheta+\xi_{k}$ (32)

where $\vartheta$ is the vector of unknown constant parameters and
$\varphi_{k}$ is the harmonic regressor presented in the following form:

$\displaystyle\varphi_{k}^{T}=[\cos(q_{0}k)\leavevmode\nobreak\
\sin(q_{0}k)\leavevmode\nobreak\ \cos(2q_{0}k)$
$\displaystyle\leavevmode\nobreak\ \sin(2q_{0}k)\leavevmode\nobreak\
...\leavevmode\nobreak\ \cos(mq_{0}k)\leavevmode\nobreak\ \sin(mq_{0}k)]$ (33)

where $q_{0}$ is the fundamental frequency of network (for example $q_{0}=50$
Hertz or $q_{0}=60$ Hertz), $m$ is the number of harmonics, and $\xi_{k}$ is a
zero mean white Gaussian noise, $k=1,2,...$ is the step number.
The model of the signal (32) with adjustable parameters $\theta_{*k}$ is
presented in the following form:

$\hat{y}_{k}=\varphi_{k}^{T}\theta_{*k}$ (34)

The signal $y_{k}$ is approximated by the model $\hat{y}_{k}$ in the least
squares sense in each step $k$ of moving window of a size $s$.
The estimation algorithm is based on minimization of the following error
$E_{k}$;

$E_{k}=\sum_{p=k-(s-1)}^{p=k}(\hat{y}_{p}-y_{p})^{2}$ (35)

for a fixed step $k$, where $k\geq s$.
The least squares solution for estimation of the parameter vector
$\theta_{*k}$, which should be calculated in each step can be written as
follows:

$\displaystyle A_{k}\theta_{*k}$ $\displaystyle=$ $\displaystyle b_{k}$ (36)
$\displaystyle A_{k}$ $\displaystyle=$
$\displaystyle\sum_{p=k-(s-1)}^{p=k}\varphi_{p}\leavevmode\nobreak\
\varphi_{p}^{T}$ (37) $\displaystyle b_{k}$ $\displaystyle=$
$\displaystyle\sum_{p=k-(s-1)}^{p=k}\varphi_{p}\leavevmode\nobreak\ y_{p}$
(38)

###  5.1 Missing Input Data & Rank Deficiency of Information Matrix

The matrix $A_{k}$ is symmetric and positive definite information matrix for
systems with harmonic regressor for a sufficiently large window size $s$.
However, the software, functional or system failures may occur, which usually
result in declaration of zero values in the information matrix $A_{k}$. Such
failures may result in rank deficiency of the information matrix $A_{k}$ in
each step. Nevertheless, the methods described above can be still applied in
this case.

###  5.2 Fault Detection on Real Data

Power quality monitoring and detection of sag and swell events is considered
in the presence of system failures. Monitoring is performed via estimation of
the frequency content of the voltage signals in three phase system. Real data
obtained from the DOE/EPRI National Database of Power System Events, [31] is
used for power quality monitoring. The event (Event ID $2911$) which resulted
in voltage sag and swell was considered. The cause of this event is unknown
and the event is considered as temporary fault. Measured voltage signals which
correspond to the sag and swell in phases ’a’,’b’ and ’c’ are plotted with
green, red and blue lines respectively in Figure 1(a).
It is assumed that the voltage signals can be described by the model (32) with
harmonic regressor (33) and fundamental frequency of $60$ Hertz with four
higher harmonics and unknown parameter vector $\vartheta$. The model of the
signals is given by (34) with the vector of the parameters $\theta_{*k}$ which
estimates unknown vector $\vartheta$. The estimation problem is reduced to the
restricted linear system of equations (36) which should be solved w.r.t.
$\theta_{*k}$ in each step. The $10\leavevmode\nobreak\
\times\leavevmode\nobreak\ 10$ matrix $A_{k}$ is the SPD matrix in each step
(which is verified by simulations).
Suppose that a functional or system failure occurred, which resulted in three
zero columns (columns 3,4 and 5) of the information matrix $A_{k}$ in each
step.
In addition, suppose that three components of the vector $b_{k}$ (the
components 3,4 and 5) are zeros in each step due to similar failure. Such
failures result in rank deficient information matrix $A_{k}$ in each step. The
rank of the information matrix is equal to seven in this case. Nevertheless,
the methods described above can be applied in this case.

###  5.3 Pre-Conditioning for Rank Deficient Information Matrix

In order to apply the methods described above the matrix $S_{k}^{-1}$ should
be chosen such that the spectral radius of the matrix $\displaystyle
I-S_{k}^{-1}A_{k}$ is less than one in each step,
$\displaystyle\rho(I-S_{k}^{-1}A_{k})<1$. For SPD information matrices the
pre-conditioner can be chosen as $\displaystyle S_{k}^{-1}=I/\alpha_{k}$ with
$\alpha_{k}=\|A_{k}\|_{\infty}/2+\varepsilon$, where $\|\cdot\|_{\infty}$ is
the maximum row sum matrix norm, $\varepsilon>0$, [7], [32]. In the case of
software or system failures the information matrix may become rank deficient
and special preconditioning should be applied. Consider the following
transformation of the system (36) and regularization of the matrix $A_{k}$:

$\displaystyle A_{rk}\theta_{*k}$ $\displaystyle=A^{T}_{k}b_{k}$ (39)
$\displaystyle A_{rk}$ $\displaystyle=\beta I+A^{T}_{k}A_{k}$ (40)

where $\beta$ is a small positive number, and the matrix $A_{rk}$ is SPD
matrix for which the preconditioning described above can be applied. The
choice of the parameter $\beta$ represents the trade-off between high
condition number (for very small $\beta$) and the performance of
regularization and hence the performance of the event detection (for large
$\beta$). The condition number of $A_{rk}$ as a function of the regularization
parameter $\beta$ is shown in Figure 2. For a small $\beta$ which allows high
performance detection the matrix $A_{rk}$ becomes ill-conditioned matrix
(which is sensitive to numerical operations) for which the spectral radius
$\displaystyle\rho(I-S_{k}^{-1}A_{rk})$ is very close to one. Application of
the stepwise splitting method, where the matrix $A$ is split recursively in
each step provides robust solution for ill-conditioned case, [32].
Notice that the model (39), (40) provides robust solution for parameter
estimation problems in the case of insufficient excitation and singular
information matrices and can be applied instead of classical self-tuning
adaptive and estimation algorithms in adaptive control and system
identification, [1] \- [3].

###  5.4 Simulation Results

Amplitudes of the first harmonic of the measured signals recovered by the
Richardson parameter estimation algorithm (31) with the second order matrix
inversion method (6) are presented in Figure 1(b). The Figure shows that the
sag and swell events are detectable using the approach described above even in
the presence of system failures, which result in rank deficiency of
information matrix. The estimation accuracy was evaluated as the average sum
of squares of the parameter errors for Richardson algorithm with respect to
parameter vector calculated without missing data. The parameter vector
calculated with Richardson algorithm is quite close to the parameter vector
calculated without missing data.
In addition, LU decomposition method was also applied for parameter
calculation in (36) with regularized matrix $A_{rk}$, which is close to
singular matrix (extremely ill-conditioned) for sufficiently small $\beta$. LU
decomposition method is very sensitive to numerical errors for extremely ill-
conditioned matrices, whereas Richardson algorithm delivers robust solution.
Therefore Richardson algorithms are recommended for robust detection of the
failure events in the presence of missing data.

##  6 CONCLUSIONS

The practical guide for solution of the parameter estimation problem described
by restricted linear system of equations was presented in the Richardson
framework, which provides the best possible solution in finite digit
calculations. Newton-Schulz and Durand iterations with unified factorizations,
which enhance the performance of the Richardson method are presented in the
form of tutorial (in historical framework) with the convergence analysis.
The combination of Richardson method with second order Newton-Schulz approach
was tested in the failure detection problem on real signals from the
electrical network. In addition, the software or system failure was simulated
which resulted in rank deficiency of the information matrix. New
preconditioning method was proposed for rank deficient matrices in the problem
of estimation of the frequency contents of electrical signals. It was shown
that the Richardson method is applicable in the failure detection problems
with rank deficient information matrices.

## References

  * [1] Fomin V., Fradkov A. and Yakubovich V., Adaptive Control of Dynamic Objects, Nauka, Moscow, 1981 (in Russian).
  * [2] Ioannou P. and Sun J., Robust Adaptive Control, Prentice Hall, Upper Saddle River, NJ, 1996.
  * [3] Ljung L., System Identification: Theory for the User, Prentice-Hall, Upper Saddle River, NJ, 1999.
  * [4] Bayard, D., A General Theory of Linear Time-Invariant Adaptive Feedforward Systems with Harmonic Regressors. IEEE Trans. Autom. Control vol. 45, N 11, pp. 1983-1996, 2000.
  * [5] Stotsky A., Recursive Trigonometric Interpolation Algorithms, Journal of Systems and Control Engineering, vol. 224, N 1, pp. 65-77, 2010.
  * [6] Björck Å., Numerical Methods for Least Squares Problems, SIAM, First edition, April 1, 1996.
  * [7] Stotsky A., Accuracy Improvement in Least-Squares Estimation with Harmonic Regressor: New Preconditioning and Correction Methods, 54-th CDC, Dec. 15-18, Osaka, Japan, 2015, pp. 4035-4040.
  * [8] Stotsky A., Unified Frameworks for High Order Newton-Schulz and Richardson Iterations: A Computationally Efficient Toolkit for Convergence Rate Improvement, Journal of Applied Mathematics and Computing, vol. 60, N 1 - 2, pp. 605-623, 2019.
  * [9] Sticko S., Towards Higher Order Immersed Finite Elements for the Wave Equation, IT Licentiate theses, 2016-008, Uppsala University, Department of Information Technology, 2016.
  * [10] Moustafa K., Real-Time Identification of Ill-Conditioned Structures, Journal of Offshore Mechanics and Arctic Engineering, August, vol. 112, pp. 181 - 186, 1990.
  * [11] Hasan A., Kerrigan E. , Constantinides G., Solving a Positive Definite System of Linear Equations via the Matrix Exponential, 50-th IEEE CDC-ECC, Orlando, FL, USA, December 12-15, 2011, pp. 2299-2304.
  * [12] Kokotovic P., Applications of Singular Perturbation Techniques to Control Problems, SIAM Review, vol. 26, N 4, pp. 501-550, 1984.
  * [13] Dubois D., Greenbaum A. and Rodrigue G., Approximating the Inverse of a Matrix for Use in Iterative Algorithms on Vector Processors, Computing, 22, pp. 257-268, 1979.
  * [14] Srivastava S., Stanimirovic P., Katsikis V. and Gupta D., A Family of Iterative Methods with Accelerated Convergence for Restricted Linear System of Equations, Mediterr. J. Math., vol. 14-222, pp. 1-26, 2017.
  * [15] Liu X., Du W., Yu Y. and Yonghui Qin Y., Efficient Iterative Method for Solving the General Restricted Linear Equation, Journal of Applied Mathematics and Physics, N 6, pp. 418-428, 2018.
  * [16] Stotsky A., Combined High-Order Algorithms in Robust Least-Squares Estimation with Harmonic Regressor and Strictly Diagonally Dominant Information Matrix, Proc. IMechE Part I: Journal of Systems and Control Engineering, vol. 229, N 2, pp. 184-190, 2015.
  * [17] Schulz G., Iterative Berechnung Der Reziproken Matrix, Zeitschrift für Angewandte Mathematik und Mechanik, vol. 13, pp. 57-59, 1933.
  * [18] Hotelling H., Some New Methods in Matrix Calculation, Annals of Math. Stat., vol. 14, pp. 1-34, 1943.
  * [19] Demidovich B., Maron I., Basics of Numerical Mathematics, Moscow, Fizmatgiz, 1963, (in Russian).
  * [20] Ben-Israel A. , A Note on an Iterative Method for Generalized Inversion of Matrices, Math. Comput. vol. 20, pp. 439-440, 1966.
  * [21] Isaacson E. and Keller H., Analysis of Numerical Methods, John Wiley & Sons, New York, 1966.
  * [22] Petryshyn W., On Generalized Inverses and on the Uniform Convergence of $(I-\beta K)^{n}$ with Application to Iterative Methods, J. Math. Anal. Appl., vol. 18, pp. 417-439, 1967.
  * [23] Stickel E., On a Class of High Order Methods for Inverting Matrices, ZAMM Z. Angew. Math. Mech. 67, pp. 331-386, 1987.
  * [24] Durand E., Solutions Numeriques des Equations Algebriques, vol. 2., Paris, Masson, 1972.
  * [25] Brezinski C. Variations on Richardson’s Method and Acceleration, in: Numerical Analysis, A Numerical Analysis Conference in Honour of Jean Meinguet, Bull. Soc. Math. Belgium, 1996, pp. 33-44.
  * [26] O’Leary, D., Yet Another Polynomial Preconditioner for the Conjugate Gradient Algorithm, Linear Algebra Appl. vol. 154-156, pp.377-388, 1991.
  * [27] Stotsky A. , Efficient Iterative Solvers in the Least Squares Method, Proc. of the 21-st IFAC World Congress, Berlin, Germany, July 12-17,2020, pp. 901-906, see also arXiv:2008.11480v1 [math.OC], 26 Aug 2020.
  * [28] Soleymani, F., Stanimirovic, P., Haghani, F., On Hyperpower Family of Iterations for Computing Outer Inverses Possessing High Efficiencies, Linear Algebra Appl., vol. 484,pp. 477-495, 2015.
  * [29] Ehrmann H., Konstruktion und Durchführung von Iterationsverfahren höherer Ordnung. Arch. Ration. Mech. Anal., N. 4, pp. 65-88,1959.
  * [30] Stanimirovic P., Kumar A. and Katsikis V., Further Efficient Hyperpower Iterative methods for the Computation of Generalized Inverses $A^{2}_{T,S}$, RACSAM, vol.113, pp. 3323-3339, 2019.
  * [31] DOE/EPRI National Database Repository of Power System Events.
  * [32] Stotsky A., Grid Frequency Estimation Using Multiple Model with Harmonic Regressor: Robustness Enhancement with Stepwise Splitting Method, IFAC PapersOnLine 50-1, pp. 12817 - 12822, 2017.
  * [33] Stotsky A. (2022). Recursive Versus Nonrecursive Richardson Algorithms: Systematic Overview, Unified Frameworks and Application to Electric Grid Power Quality Monitoring, Automatika, vol. 63, N2, pp. 328-337, https://doi.org/10.1080/00051144.2022.2039989
  * [34] Stotsky A.,Simultaneous Frequency and Amplitude Estimation for Grid Quality Monitoring : New Partitioning with Memory Based Newton-Schulz Corrections, IFAC PapersOnLine 55-9, pp. 42-47, 2022, https://www.sciencedirect.com/science/article/pii/S2405896322003937


# mSHAP: SHAP Values for Two-Part Models

Spencer Matthews  Brian Hartman

###### Abstract

Two-part models are important to and used throughout insurance and actuarial
science. Since insurance is required for registering a car, obtaining a
mortgage, and participating in certain businesses, it is especially important
that the models which price insurance policies are fair and non-
discriminatory. Black box models can make it very difficult to know which
covariates are influencing the results. SHAP values enable interpretation of
various black box models, but little progress has been made in two-part
models. In this paper, we propose mSHAP (or multiplicative SHAP), a method for
computing SHAP values of two-part models using the SHAP values of the
individual models. This method will allow for the predictions of two-part
models to be explained at an individual observation level. After developing
mSHAP, we perform an in-depth simulation study. Although the kernelSHAP
algorithm is also capable of computing approximate SHAP values for a two-part
model, a comparison with our method demonstrates that mSHAP is exponentially
faster. Ultimately, we apply mSHAP to a two-part ratemaking model for personal
auto property damage insurance coverage. Additionally, an R package (mshap) is
available to easily implement the method in a wide variety of applications.

###### keywords:

Explainability , Machine Learning , Ratemaking
Declarations of interest: none

††journal: Computational Statistics & Data Analysis

[inst1]organization=Department of Statistics, Brigham Young
University,city=Provo, postcode=84602, state=UT, country=Utah County

##  1 Introduction

One of the most popular families of machine learning models are tree-based
algorithms, which use the concept of many decision trees working together to
create more generalized predictions [9]. Current implementations include
random forests, gradient boosted forests, and others. These models are very
good at learning relationships and have proven highly accurate in diverse
areas. Already, many aspects of life are affected by these algorithms as they
have been implemented in business, technology, and more.

As these methods become more abundant, it is crucial that explanations of
model output are easily available. The exact definition of “explanation” is a
subject of debate, and Lipton [8] argues that the word is often used in a very
unscientific way due to the confusion over its meaning. In this paper, we will
regard an explainable system as what Doran et al. [3] refer to as a
comprehensible system, or one that “allow[s] the user to relate properties of
the inputs to their output.”

Explainable models are important not only because some industries require
them, but also because understanding the why behind the output is essential to
avoiding possible pitfalls. Understanding the reasoning behind model output
allows for recognition of model bias, and increased security against harmful
models being put into production. When implemented well, machine learning
models can be more accurate than traditional models. However, more accurate
model families can be less explainable, simply because of the nature of these
algorithms. Generally, as predictive performance increases so too does model
complexity, decreasing the ability to understand the effects of inputs on the
output [5].

In this paper, we propose a methodology for explaining two-part models, which
expands on the already prevalent TreeSHAP algorithm [9]. This methodology,
called mSHAP, will allow the output of two models to be multiplied together
while maintaining explainability of the resulting prediction.

In order to understand the implications and advantages of mSHAP, we will first
revisit the existing SHAP-based methods and discuss where issues arise in the
context of two-part models. We will then propose a context in which SHAP
values for two existing models can be combined to explain a two-part model.
Although this framework is robust, it does leave a part (which we call
$\alpha$) of the ultimate prediction that must be distributed back into the
contributions of the variables. To this end, we run a simulation across
different methods of distributing $\alpha$ and score the methods in comparison
to kernelSHAP, an existing method for estimating explanations of any type of
model. Included in this simulation is a comparison of mSHAP to kernelSHAP,
which highlights the advantages of mSHAP. Having scored these methods, we
select the best one, and apply the process of mSHAP on an auto insurance
dataset.

##  2 Motivation

The initial idea for this methodology came due to the problem of machine
learning in auto insurance ratemaking (or pricing). Actuaries are tasked with
taking historical data and using it to set current rates for insured
consumers. Given the sensitive nature of the data and the potential impact it
has to bias rates for different types of people, there are strict regulations
on the models. The outputs of these models must be explainable so that
regulators in the insurance industry can be sure that the rates are not
unfairly discriminatory.

Many actuaries use a two-part model to set rates, where the first part
predicts how many claims a policyholder will have (the claim frequency) and
the second part predicts the average cost of an individual claim [7, 4].
Multiplying the two outputs of these models predicts the total cost of a given
policyholder.

Two-part models are more difficult to explain than standard models, but the
complexity increases when the two models themselves are not traditional
generalized linear models. Given this difficulty and the strict requirements
of the regulators, machine learning models are not often used in actuarial
ratemaking. Despite the lack of current industry use, machine learning models
such as tree-based algorithms could improve the accuracy of ratemaking models
[1]. Since the data that actuaries work with is typically tabular, tree-based
algorithms are a good fit for predicting on the data. In recent years there
have been many advances in explaining tree based machine learning algorithms,
which could lead to greater adaptation in the field. One of the most important
is the SHAP value.

###  2.1 SHAP Values and Current Implementations

SHAP values originate in the field of economics, where they are used to
explain player contributions in cooperative game theory. Proposed by Shapley
[12], they predict what each player brings to a game. This idea was ported
into the world of machine learning by Lundberg and Lee [10]. The basic
algorithm calculates the contribution of a variable to the prediction for
every possible ordering of variables, then averages those contributions. This
becomes computationally impractical very quickly, but Lundberg and Lee created
a modified algorithm that approximates these SHAP values.

A couple years later Lundberg et al. [9] published a new paper detailing a
method called TreeSHAP. This method is a rapid way of computing exact SHAP
values for any tree-based machine learning model. The fixed structure of trees
in a tree-based model allows shortcuts to be taken in the computation of SHAP
values, which greatly speeds up the process. With this improvement, it becomes
feasible to explain millions of predictions from tree based machine learning
algorithms. These local explanations can then be combined to create an
understanding of the whole model.

###  2.2 Properties of SHAP Values

There are three essential properties of SHAP values: local
accuracy/efficiency, consistency/monotonicity, and missingness [10]. These
three properties are satisfied by the equation used to calculate SHAP values,
as implemented by Lundberg and Lee. While we focus on the local accuracy
property for the rest of this section, we note that since mSHAP is built on
top of treeSHAP, it automatically incorporates the consistency/monotonicity
and missingness properties.

####  2.2.1 Local Accuracy in Implementation

The most important of the above mentioned properties in the context of mSHAP
is the property of local accuracy/efficiency. In the context of machine
learning, this property says that the contributions of the variables should
add up to the difference between the prediction and the average prediction of
the model. The average prediction can be thought of as the model bias term,
which is what the model will predict, on average, across all inputs (assuming
representative training data). Note that in the TreeSHAP algorithm the average
prediction of the model is computed as the mean of all predictions for the
training data set. The SHAP values are then computed to explain deviance from
the average prediction.

Thus, given an arbitrary model $Y$ with prediction $\hat{y}$ based on two
predictors, $x_{1}$ and $x_{2}$, we can represent the mean prediction with
$\mu_{Y}$ and the SHAP values for the two covariates as $s_{x_{1}}$ and
$s_{x_{2}}$. Based on the property of local accuracy, we know that
$\hat{y}=\mu_{Y}+s_{x_{1}}+s_{x_{2}}$.

This principle applies to models with any number of predictors and is very
desireable in explainable machine learning [2].

####  2.2.2 The Problem of Local Accuracy

Since it is so important that the SHAP values add up to the model output, any
attempt at explaining two-part model output from the SHAP values of the
individual parts must maintain this property. However, multiplying the output
of two models blends the contributions from different variables, making it
unclear what contributions should be given to what variables. The idea of
combining models, and using the SHAP values of the individual models to get
the SHAP values for the combined model has been implemented before. In a
related github issue, Scott Lundberg assures that averaging model output is
compatible with averaging SHAP values, as long as the SHAP values (and model
output) are in their untransformed state [13]. Even though averaging SHAP
values for each variable works when averaging model outputs, the same
principle does not apply when multiplying model outputs.

When considered, this is apparent. In the most simple of cases, we see that if
we have two models that both predict some outcome based on two covariates
$x_{1}$ and $x_{2}$, we can average their results and likely get a better
prediction. We will call these models $A$ and $B$, respectively. For a given
observation, model $A$ predicts $\hat{a}$ and model $B$ predicts $\hat{b}$.
When run through a SHAP explainer, we can break down these predictions even
further. Since SHAP values are additive, we know that
$\hat{a}=\mu_{A}+s_{x_{1}a}+s_{x_{2}a}$ and
$\hat{b}=\mu_{B}+s_{x_{1}b}+s_{x_{2}b}$. It follows that

$\displaystyle\text{avg}(\hat{a},\hat{b})$
$\displaystyle=\frac{\hat{a}+\hat{b}}{2}$
$\displaystyle=\frac{\mu_{A}+s_{x_{1}a}+s_{x_{2}a}+\mu_{B}+s_{x_{1}b}+s_{x_{2}b}}{2}$
$\displaystyle=\frac{\mu_{A}+\mu_{B}}{2}+\frac{s_{x_{1}a}+s_{x_{1}b}}{2}+\frac{s_{x_{2}a}+s_{x_{2}b}}{2}$
$\displaystyle=\text{avg}(\mu_{A},\mu_{B})+\text{avg}(s_{x_{1}a},s_{x_{1}b})+\text{avg}(s_{x_{2}a},s_{x_{2}b}).$

This means that we can find the contribution to the overall model from $x_{1}$
by averaging $s_{x_{1}a}$ and $s_{x_{1}b}$, and likewise for the contribution
to the overall model from $x_{2}$.

However, if we for some reason wished to stack our models such that the two
outputs ($\hat{a}$ and $\hat{b}$) were multiplied, we run into a problem. This
occurs because, despite the longings of all algebra students,

$\hat{a}\hat{b}=(\mu_{A}+s_{x_{1}a}+s_{x_{2}a})(\mu_{B}+s_{x_{1}b}+s_{x_{2}b})\neq\mu_{A}\mu_{B}+s_{x_{1}a}s_{x_{1}b}+s_{x_{2}a}s_{x_{2}b}.$

Instead, we end up with

$\displaystyle\hat{a}\hat{b}$
$\displaystyle=(\mu_{A}+s_{x_{1}a}+s_{x_{2}a})(\mu_{B}+s_{x_{1}b}+s_{x_{2}b})$
$\displaystyle=\mu_{A}\mu_{B}+\mu_{A}s_{x_{1}b}+\mu_{A}s_{x_{2}b}+s_{x_{1}a}\mu_{B}+s_{x_{1}a}s_{x_{1}b}+s_{x_{1}a}s_{x_{2}b}+s_{x_{2}a}\mu_{B}+s_{x_{2}a}s_{x_{1}b}+s_{x_{2}a}s_{x_{2}b}.$

Even in this simple case, it is difficult to assign a single contribution to
our two different variables when presented with the SHAP values of the two
original models. This problem grows even more difficult with the addition of
other explanatory features. mSHAP is the methodology developed to solve this
problem.

##  3 The Math behind Multiplying SHAP Values

In a two-part model, the output of one model is multiplied by the output of a
second model to obtain the response. The principal driver behind mSHAP is the
explanation of these sorts of models, and it requires that the SHAP values be
multiplied together in some way to get a final SHAP value for the output. The
mathematics behind mSHAP are explained here in the general case, for any given
number of predictors with a training set of arbitrary size. Although an exact
solution for the SHAP values of a two-part model is still out of reach, this
method proves very accurate in its results.

###  3.1 Definitions

Consider three different models, $f,g,$ and $h$, and a single input (training)
matrix $A$. We will let the number of columns and rows in $A$ be arbitrary. In
other words, let $A$ be an $n\times p$ matrix where each column is a covariate
and each row is an observation. Also let $A_{i}$ denote the $i$th observation
(row) of $A$. Furthermore, define $h$ to be the product of $f$ and $g$, so
$h(A_{i})=f(A_{i})\cdot g(A_{i})$.

Recall that the sum of the SHAP values for each covariate and the average
model output must add up to the model prediction. For simplicity in
presentation, we will define $f(A_{i})=\hat{x_{i}}$, $g(A_{i})=\hat{y_{i}}$,
and $h(A_{i})=\hat{z_{i}}$ and the contribution of the $j$th predictor to
$x_{i}$ as $s_{x_{i}j}$. With these considerations in place, we can define the
output space of our three models on the training data set.

For model $f$:

$\displaystyle\hat{x_{1}}=$ $\displaystyle
s_{x_{1}1}+s_{x_{1}2}+s_{x_{1}3}+\ldots+s_{x_{1}p}+\mu_{f}$
$\displaystyle\hat{x_{2}}=$ $\displaystyle
s_{x_{2}1}+s_{x_{2}2}+s_{x_{2}3}+\ldots+s_{x_{2}p}+\mu_{f}$
$\displaystyle\hat{x_{3}}=$ $\displaystyle
s_{x_{3}1}+s_{x_{3}2}+s_{x_{3}3}+\ldots+s_{x_{3}p}+\mu_{f}$
$\displaystyle\vdots$ $\displaystyle\hat{x_{n}}=$ $\displaystyle
s_{x_{n}1}+s_{x_{n}2}+s_{x_{n}3}+\ldots+s_{x_{n}p}+\mu_{f}$

For model $g$:

$\displaystyle\hat{y_{1}}=$ $\displaystyle
s_{y_{1}1}+s_{y_{1}2}+s_{y_{1}3}+\ldots+s_{y_{1}p}+\mu_{g}$
$\displaystyle\hat{y_{2}}=$ $\displaystyle
s_{y_{2}1}+s_{y_{2}2}+s_{y_{2}3}+\ldots+s_{y_{2}p}+\mu_{g}$
$\displaystyle\hat{y_{3}}=$ $\displaystyle
s_{y_{3}1}+s_{y_{3}2}+s_{y_{3}3}+\ldots+s_{y_{3}p}+\mu_{g}$
$\displaystyle\vdots$ $\displaystyle\hat{y_{n}}=$ $\displaystyle
s_{y_{n}1}+s_{y_{n}2}+s_{y_{n}3}+\ldots+s_{y_{n}p}+\mu_{g}$

And for model $h$:

$\displaystyle\hat{z_{1}}=$ $\displaystyle
s_{z_{1}1}+s_{z_{1}2}+s_{z_{1}3}+\ldots+s_{z_{1}p}+\mu_{h}$
$\displaystyle\hat{z_{2}}=$ $\displaystyle
s_{z_{2}1}+s_{z_{2}2}+s_{z_{2}3}+\ldots+s_{z_{2}p}+\mu_{h}$
$\displaystyle\hat{z_{3}}=$ $\displaystyle
s_{z_{3}1}+s_{z_{3}2}+s_{z_{3}3}+\ldots+s_{z_{3}p}+\mu_{h}$
$\displaystyle\vdots$ $\displaystyle\hat{z_{n}}=$ $\displaystyle
s_{z_{n}1}+s_{z_{n}2}+s_{z_{n}3}+\ldots+s_{z_{n}p}+\mu_{h}$

Furthermore, given our training data $A$, we can extract the values of
$\mu_{f},\mu_{g},$ and $\mu_{h}$. As explained above, these are the average
value of the model predictions on the training set.

$\displaystyle\mu_{f}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\hat{x_{i}}=\frac{\hat{x_{1}}+\hat{x_{2}}+\hat{x_{3}}+\ldots+\hat{x_{n}}}{n}$
$\displaystyle\mu_{g}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\hat{y_{i}}=\frac{\hat{y_{1}}+\hat{y_{2}}+\hat{y_{3}}+\ldots+\hat{y_{n}}}{n}$
$\displaystyle\mu_{h}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\hat{z_{i}}=\frac{\hat{z_{1}}+\hat{z_{2}}+\hat{z_{3}}+\ldots+\hat{z_{n}}}{n}$

In practice, it is necessary to be able to pull $\mu_{h}$ out of
$\hat{x_{i}}\hat{y_{i}}$. When implemented, it is important to note that
$\mu_{f}\mu_{g}=\mu_{f}\mu_{g}-\mu_{h}+\mu_{h}$. Since every expansion of SHAP
values from $\hat{x_{i}}\hat{x_{j}}$ contains $\mu_{f}\mu_{g}$, we substitute
in $\mu_{f}\mu_{g}-\mu_{h}+\mu_{h}$, where $\mu_{h}$ is essential and
$\mu_{f}\mu_{g}-\mu_{h}$ becomes a term that we label $\alpha$ and distribute
among all the SHAP values. A more formalized definition of $\alpha$ is given
in B.

###  3.2 Obtaining $z_{i}$’s SHAP Values

We now derive the individual SHAP values for each variable as it pertains to
the prediction of model $h$. Again, we will let this output be an arbitrary
$\hat{z_{i}}$. Recall that

$\hat{z_{i}}=\hat{x_{i}}\hat{y_{i}}=(s_{x_{i}1}+s_{x_{i}2}+s_{x_{i}3}+\ldots+s_{x_{i}p}+\mu_{f})(s_{y_{i}1}+s_{y_{i}2}+s_{y_{i}3}+\ldots+s_{y_{i}p}+\mu_{g}).$

Using a tabular form for visual simplicity, we obtain the following expansion:

$\begin{array}[]{c|ccccccccccc}&s_{x_{i}1}&+&s_{x_{i}2}&+&s_{x_{i}3}&+&\ldots&+&s_{x_{i}p}&+&\mu_{f}\\\\[6.0pt]
\hline\cr\\\
s_{y_{i}1}&s_{x_{i}1}s_{y_{i}1}&&s_{x_{i}2}s_{y_{i}1}&&s_{x_{i}3}s_{y_{i}1}&&\ldots&&s_{x_{i}p}s_{y_{i}1}&&\mu_{f}s_{y_{i}1}\\\
+\\\
s_{y_{i}2}&s_{x_{i}1}s_{y_{i}2}&&s_{x_{i}2}s_{y_{i}2}&&s_{x_{i}3}s_{y_{i}2}&&\ldots&&s_{x_{i}p}s_{y_{i}2}&&\mu_{f}s_{y_{i}2}\\\
+\\\
s_{y_{i}3}&s_{x_{i}1}s_{y_{i}3}&&s_{x_{i}2}s_{y_{i}3}&&s_{x_{i}3}s_{y_{i}3}&&\ldots&&s_{x_{i}p}s_{y_{i}3}&&\mu_{f}s_{y_{i}3}\\\
+\\\ \vdots&\vdots&&\vdots&&\vdots&&\ddots&&\vdots&&\vdots\\\ +\\\
s_{y_{i}n}&s_{x_{i}1}s_{y_{i}p}&&s_{x_{i}2}s_{y_{i}p}&&s_{x_{i}3}s_{y_{i}p}&&\ldots&&s_{x_{i}p}s_{y_{i}p}&&\mu_{f}s_{y_{i}p}\\\
+\\\
\mu_{g}&s_{x_{i}1}\mu_{g}&&s_{x_{i}2}\mu_{g}&&s_{x_{i}3}\mu_{g}&&\ldots&&s_{x_{i}p}\mu_{g}&&\mu_{f}\mu_{g}\end{array}$

We break these terms into the SHAP values for each variable, 1 through $p$,
for $\hat{z_{i}}$. Our approach breaks $s_{z_{i}j}$ into two parts, which we
call $s_{z_{i}j}^{\prime}$ and $\alpha_{ij}$. Though the method of obtaining
$\alpha_{i}$ can take on several forms, $s_{z_{i}j}^{\prime}$ is always as
follows (where $j$ refers to the $j$th covariate):

$\displaystyle s_{z_{i}j}^{\prime}$
$\displaystyle=\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+s_{x_{i}j}s_{y_{i}j}+\sum^{p}_{a=1}(\frac{s_{x_{i}j}s_{y_{i}a}}{2}I(a\neq
j))+\sum^{p}_{a=1}(\frac{s_{y_{i}j}s_{x_{i}a}}{2}I(a\neq j))$
$\displaystyle=\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})$

In words, and with the aid of the table above, this can be described as the
sum of the $i$th row and $i$th column, where every term is divided by two
except the terms with $\mu_{f}$ and $\mu_{g}$. When applied to each variable,
this can be written as

$\hat{z_{i}}=\sum^{p}_{j=1}\left[\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]+\mu_{f}\mu_{g}$

and applying the breakdown we derived above, while simplifying as well,

$\hat{z_{i}}=\left(\sum^{p}_{j=1}s_{z_{i}j}^{\prime}\right)+\alpha+\mu_{h}.$

For a proof that this formula and the subsequent distribution of $\alpha$
maintains the local accuracy property of SHAP values, see B.1.

###  3.3 Methods for distributing $\alpha$

We now arrive at the aforementioned point of deciding how to distribute
$\alpha$ into each $s_{z_{i}j}$. There are four ways which we tested for
distributing $\alpha$, the first being simple uniform distribution, and the
others being variations of weighting based on the value of
$s_{z_{i}j}^{\prime}$. Note that all four of these methods maintain the local
accuracy property of SHAP values.

####  3.3.1 Uniformly Distributed

The simplest way of distributing $\alpha$ between all the $s_{z_{i}j}$’s is to
divide it evenly. In this case, our resulting equation for each variable’s
SHAP value would be

$s_{z_{i}j}=s_{z_{i}j}^{\prime}+\frac{\alpha}{p}.$

This method could prove a strong baseline.

####  3.3.2 Raw Weights

The computation of this method is made easier by recalling that

$\sum^{p}_{j=1}s_{z_{i}j}^{\prime}=\hat{z_{i}}-\mu_{f}\mu_{g},$

which allows us to use $\hat{z_{i}}-\mu_{f}\mu_{g}$ as the whole upon which we
base our weighting. When applied, this method defines each SHAP value as

$s_{z_{i}j}=s_{z_{i}j}^{\prime}+\frac{s_{z_{i}j}^{\prime}}{\hat{z_{i}}-\mu_{f}\mu_{g}}(\alpha).$

####  3.3.3 Absolute Weights

This method differs from that of the raw weights, in that instead of summing
the $s_{z_{i}j}^{\prime}$’s we sum their absolute values. The weight for each
SHAP value is calculated with

$s_{z_{i}j}=s_{z_{i}j}^{\prime}+\frac{|s_{z_{i}j}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha).$

####  3.3.4 Squared Weights

Finally, instead of working with the absolute values, we could work with
squares. Similar to the equation above, the SHAP values under this method are
computed by

$s_{z_{i}j}=s_{z_{i}j}^{\prime}+\frac{(s_{z_{i}j}^{\prime})^{2}}{\sum^{p}_{k=1}(s_{z_{i}k}^{\prime})^{2}}(\alpha).$

##  4 Simulation Study for Distributing $\alpha$

To test the differences between these methods of distributing $\alpha$, we
simulated various multiplicative models based on known equations, and compared
the results of our multiplicative method with the output from kernelSHAP.
KernelSHAP is an existing generalized method for estimating the contributions
based on any prediction function. However, it is extremely computationally
expensive when compared with TreeSHAP. When training on millions of rows with
many variables, it becomes unrealistic to use kernelSHAP for computing the
SHAP values.

###  4.1 Scoring the Methods

Several factors were considered in scoring, including the mean absolute error
of the SHAP values, the directions of the SHAP values, and the rank (in
magnitude) of the SHAP values for each variable. The score needed to be a
singular way to asses how close the method gets to the kernelSHAP estimates.
Even though kernelSHAP is an estimate and not necessarily the truth, we used
it as a benchmark in the different parts of our score. This allowed us to
compare the new variations of the mSHAP method to existing methods for the
computation of SHAP values.

For ease of notation, if we define the SHAP value we are estimating as
$s_{z_{i}j}$, then we can define its counterpart as computed by kernelSHAP, as
$k_{z_{i}j}$.

####  4.1.1 General Equation for Scoring

In the end, an equation was formed to create a raw “score” based on the
direction of the SHAP value, the relative value of the SHAP value, and the
rank (importance) of the SHAP value in comparison to kernelSHAP. The score
ranges from 0 to 3 (with 3 being the best possible score), and is defined by

$\displaystyle\beta(s_{z_{i}j},k_{z_{i}j}|\theta_{1},\theta_{2})$
$\displaystyle=\lambda_{1}(s_{z_{i}j},k_{z_{i}j}|\theta_{1})+\lambda_{2}(s_{z_{i}j},k_{z_{i}j}|\theta_{2})+\lambda_{3}(s_{z_{i}j},k_{z_{i}j})$

where

$\displaystyle\lambda_{1}(s_{z_{i}j},k_{z_{i}j}|\theta_{1})$
$\displaystyle=\begin{cases}1&\quad s_{z_{i}j}k_{z_{i}j}>0\\\
\text{min}(1,\frac{1+\theta_{1}}{|s_{z_{i}j}|+|k_{z_{i}j}|+\theta_{1}})&\quad\text{otherwise}\\\
\end{cases}$ $\displaystyle\lambda_{2}(s_{z_{i}j},k_{z_{i}j}|\theta_{2})$
$\displaystyle=\text{min}(1,\frac{1+\theta_{2}}{|s_{z_{i}j}-k_{z_{i}j}|+1})$
$\displaystyle\lambda_{3}(s_{z_{i}j},k_{z_{i}j})$
$\displaystyle=\frac{1}{|\text{imp}(s_{z_{i}j})-\text{imp}(k_{z_{i}j})|+1}$

and imp$(s_{z_{i}j})$ is the importance of that SHAP value relative to the
other contributions in the observation (where importance is determined by
absolute value).

In this function (and as will be described in the following section),
$\lambda_{1}$ is the contribution from the signs of the SHAP values,
$\lambda_{2}$ is the contribution from the relative value of the SHAP values,
and $\lambda_{3}$ is the contribution from the relative ranking (importance)
of the SHAP values.

####  4.1.2 Lambda Functions

To gain some intuition about the $\lambda$ functions and the impact of
$\theta_{1}$ and $\theta_{2}$ we depict them in Figure 1.

For $\lambda_{1}$, which measures whether the two SHAP values are the same
sign, any values in the first and third quadrants return a perfect score of 1,
since the two values have the same sign. It also allows for some wiggle room
with $\theta_{1}$, by allowing anything within the lines
$k_{z_{i}j}=s_{z_{i}j}+\theta_{1}$ and $k_{z_{i}j}=s_{z_{i}j}-\theta_{1}$ to
be 1. Beyond those boundaries, the scores gradually decrease.

The function $\lambda_{2}$, which compares the values, also creates boundary
lines for the perfect score of 1 at $k_{z_{i}j}=s_{z_{i}j}+\theta_{2}$ and
$k_{z_{i}j}=s_{z_{i}j}-\theta_{2}$. In other words, as long as the difference
between $s_{z_{i}j}$ and $k_{z_{i}j}$ is less than $\theta_{2}$, the function
will return 1. Beyond that and the value begins to decrease.

(a) $\lambda_{1}$

(b) $\lambda_{2}$

Figure 1: Heat maps for the $\lambda$ functions.

Out of the three $\lambda_{3}$, the rank measure, is the easiest to
understand. In a given observation, each SHAP value is given a rank (between 1
and $p$, inclusive) based on its absolute value. These ranks are then
compared, and the closer they are together, the higher the score, with a
perfect score of 1 being obtained if the two rankings are the same.

###  4.2 Simulation Study

As mentioned above, we simulated various multiplicative models based on known
equations, and compared the results of our multiplicative method with the
output from kernelSHAP in order to test the model.

Specifically, we used three variables, $x_{1},x_{2},x_{3}$ in a variety of
response equations $y_{1},y_{2}$ to create models for $y_{1}$ and $y_{2}$ and
then multiply their outputs together. Using the multiplied output and the
covariates, we were able to use kernelSHAP to compute an estimate of the SHAP
values. We could then compare this estimate to the result from our
multiplicative method, as described above, with different ways of distributing
$\alpha$ applied.

More details on the simulation can be found in C.

For testing, we used 100 samples in each iteration for faster computation,
which allowed us to simulate over 2500 scenarios. Specifically, we worked with
all possible combinations of the following values.

Variable | Possible Values
---|---
$y_{1}$ | $x_{1}+x_{2}+x_{3}$
| $2*x_{1}+2*x_{2}+3*x_{3}$
$y_{2}$ | $x_{1}+x_{2}+x_{3}$
| $2*x_{1}+2*x_{2}+3*x_{3}$
| $x_{1}*x_{2}*x_{3}$
| $x_{1}^{2}*x_{2}^{3}*x_{3}^{4}$
| $(x_{1}+x_{2})/(x_{1}+x_{2}+x_{3})$
| $x_{1}*x_{2}/(x_{1}+x_{1}*x_{2}+x_{1}^{2}*x_{3}^{2})$
$\theta_{1}$ | 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5,
| 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5
$\theta_{2}$ | 1, 6, 11, 16, 21, 26, 31, 36, 41, 46

For each combination of values in the above table, we distributed $\alpha$ in
each of the four ways mentioned in Section 3.3. The resulting table,
therefore, had results for each model and each way of distributing $\alpha$.
In general, we averaged across all rows of the same method to obtain the
scores that were compared to each other.

It should also be noted that in our examples, our covariates were distributed
as follows:

$\displaystyle x_{1}$ $\displaystyle\sim\text{Uniform}[-10,10]$ $\displaystyle
x_{2}$ $\displaystyle\sim\text{Uniform}[0,20]$ $\displaystyle x_{3}$
$\displaystyle\sim\text{Uniform}[-5,-1].$

###  4.3 Results of the Simulation

In general, the multiplicative SHAP method did very well, when compared to the
kernelSHAP output. Since kernelSHAP is an estimation as well, it is hard to
determine exactly how well the multiplicative SHAP method does, but we will
summarize some statistics here.

####  4.3.1 Distributing $\alpha$

After trying the aforementioned four methods for distributing $\alpha$ into
the SHAP values, we came to the conclusion that the weighted by absolute value
method was the best. This came by way of the score as well as other metrics.
Details can be seen in the following table (please note that all values are
averaged across all 2520 simulations).

Method | Score |  Direction Score |  Relative Value Score | Rank Score |  Pct Same Sign |  Pct Same Rank
---|---|---|---|---|---|---
Weighted by Absolute Value | 2.27 | 0.869 | 0.594 | 0.802 | 84.8% | 62.5%
Weighted by Squared Value | 2.21 | 0.841 | 0.579 | 0.792 | 81.8% | 60.8%
Uniformly Distributed | 2.20 | 0.858 | 0.563 | 0.783 | 83.7% | 59.4%
Weighted by Raw Value | 1.99 | 0.727 | 0.494 | 0.768 | 71.4% | 56.2%

####  4.3.2 Impact of $\theta_{1}$ and $\theta_{2}$

We plotted the effects of the different values for $\theta_{1}$ and
$\theta_{2}$ on the overall score, based on type of method of distribution.

Figure 2: How $\theta_{1}$ impacts overall score, on average Figure 3: How
$\theta_{2}$ impacts overall score, on average

As can be seen, in Figures 2 and 3, changing the value of these two parameters
has a similar impact across all scoring methods.

####  4.3.3 Computational Time

The most dramatic benefit of mSHAP over kernelSHAP is the computational
efficiency of mSHAP.

(a) Fixed $n$

(b) Fixed Variables

Figure 4: Computational time of kernelSHAP and mSHAP

In Figure 4, we are able to see the comparison in run time between kernelSHAP
and the mSHAP method (including the individual treeSHAP value calculations).
Both increasing the number of variables and the number of samples causes the
time of kernelSHAP to grow greatly, while the multiplicative method remains
fairly constant. In these trials, the number of background samples was fixed
at 100 for kernelSHAP.

A case study can show the importance of this. In the auto insurance data set,
there are 5,000,000 rows in the test set, with 46 variables. For the sake of
simplicity, let’s assume that we use 45 of those variables, and that 100
background samples is enough to compute accurate SHAP values. In reality, it
would need many more background samples, but that only accentuates the point,
as a large quantity of background samples slows kernelSHAP drastically.
KernelSHAP computes SHAP values for 45 variables at a rate of about 2.268
seconds per observation on a personal laptop. In order to compute the SHAP
values for the entire test set, then, one would need about 131 days of
continuous compute time.

In contrast, our multiplicative method, using treeSHAP on two tree-based
models, computes SHAP values at a rate of about .00175 seconds per observation
for a model with 45 variables. To compute the SHAP values for the entire test
set using this method would take a little less than three hours of continuous
compute time.

###  4.4 Final Equation for mSHAP

Based on the results of the simulation, we determine that the best way of
distributing $\alpha$ is the method of weighting by absolute values (as
described above). Recall from above that in this method,

$s_{z_{i}j}=s_{z_{i}j}^{\prime}+\frac{|s_{z_{i}j}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha)$

and that $s_{z_{i}j}^{\prime}$ refers to an initial mSHAP value, before the
correction introduced by $\alpha$. It is calculated as

$s_{z_{i}j}^{\prime}=\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a}).$

Thus, the final equation for the mSHAP value of the $j$th predictor on the
$i$th observation can be written as

$s_{z_{i}j}=\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]+\frac{|s_{z_{i}j}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha).$

For a complete proof that local accuracy holds with this equation, see B.1.

##  5 Case Study

In order to prove the efficacy of mSHAP it is necessary to put it into
practice. We obtained an insurance dataset including over 20 million auto
insurance policies for a large insurance provider in the United States. Using
this data, we created a two-part model that predicts the expected property
damage cost of each policy. Both parts of this model consist of tree-based
methods, specifically random forests. After creating this model, we used the
shap python library to explain the predictions of each individual part on a
sample of 50,000 observations from our test set. We then applied the final
mSHAP method as described above, to obtain explanations for the overall model
and used the mshap R package to visualize some of the results.

###  5.1 Model Creation

As mentioned above, the model is a two-part model for predicting the expected
cost of the policy. The first part of the model predicts the frequency of the
claims. It is a random forest that predicts the probability of each of 4
possible outcomes (a multinomial model). In our data set, there existed
policies with up to 7 claims, but we chose the classes of 0, 1, 2, and 3 and
bundled everything over 3 into the 3 class. The data was heavily imbalanced,
so we used a combination of upsampling the minority classes (1, 2, 3 claims)
and downsampling the majority class (0 claims) to get a more balanced training
data set. This allowed the model to use the information to predict meaningful
probabilities instead of always assigning a very high probability to 0 claims.

The second part is a random forest which predicts the severity component of
the two-part model, or the expected cost per claim.

Once these models were created, we could calculate the expected value (or in
this case, the expected cost) of a policy in the following manner. If we let
$\hat{P_{i}}(a)$ denote the predicted probability of the for the $i$th policy
of the $a$th class and $\hat{y_{i}}$ be the predicted severity of the policy,
then

$EV=\hat{y_{i}}\left(0\hat{P_{i}}(0)+1\hat{P_{i}}(1)+2\hat{P_{i}}(2)+3\hat{P_{i}}(3)\right)$

The final two-part model was used to predict the expected cost of 50,000
policies from the test data set. For more specific details about the model and
how it was tuned, see E.

###  5.2 Model Explanation

After creating the two-part model and getting final predictions for the
expected cost of the claims, we were able to apply mSHAP to explain the final
model predictions. Before doing this, we computed the SHAP values on the
individual models so that we have the necessary data to apply the mSHAP method
for explaining two-part models. Summary plots for the five different sets of
SHAP values (one for severity, and one for each class of the frequency model)
can be created. In Figure 5 we depict the SHAP values for one of the frequency
classes from the frequency model and the SHAP values for the severity model.

(a) Summary plot of the frequency model’s SHAP values for the 0 claim class

(b) Summary plot of the severity model’s SHAP values

Figure 5: Example summary plots of SHAP values from the individual model parts

After computing these SHAP values, we applied the mSHAP method detailed in
this paper. Note that after applying mSHAP, the expected value formula above
is simply a linear combination, and we are able to perform that same linear
combination on the SHAP values before (or after) applying mSHAP. This process
left us with a single mSHAP value for each variable in every row of our test
set, and an overall expected value across the training set. The summary plot
of those final mSHAP values can be seen in Figure 6, and an example of an
observation plot is shown in Figure 7.

Figure 6: Summary plot of the two-part model’s mSHAP values. Figure 7:
Observation plot from the two-part model’s mSHAP values. This plot shows how
mSHAP can be used to explain a single observation.

The beauty of the mSHAP method is it allows for a two-part model to be
explained in the same ways that tree-based models can be easily explained with
SHAP values. As can be seen in the plots, general trends across variables can
be established, as well as specific policies dissected to see the individual
motivators behind each prediction. The ability of mSHAP to explain these types
of models opens the door to using two-part models that are both powerful and
explainable.

##  6 Conclusion

In this paper we developed mSHAP, a method for calculating SHAP values in two-
part models. Our method is much quicker than kernelSHAP and may lead to
increased adoption of machine learning (especially tree-based) components of
two-part models. This is particularly true in insurance where issues of
fairness and discrimination are very important.

##  Appendix A Shapley Values

In this section we briefly discuss the math behind Shapley values. Please note
that this section leans heavily upon the explanations and formulas as given in
Lundberg and Lee [10]. A motivated reader will find further information
regarding Shapley values in that paper.

Shapley values are a class of what is known as additive feature attribution
methods. These methods are defined as methods that have an “explanation model
that is a linear function of binary variables:”

$g(z^{\prime})=\phi_{0}+\sum^{M}_{i=1}\phi_{i}z^{\prime}_{i}$

where $M$ is the number of input features, $\phi_{0},\phi_{i}\in\mathbb{R}$
and $z^{\prime}_{i}\in\\{0,1\\}^{M}$. Essentially, every prediction of the
model (which we will denote $f(x)$), can be obtained by assigning some
contribution to each of the variables.

The Shapley values have three desirable properties, as mentioned above, and
the formal definitions for these properties are given here.

Local Accuracy. Local accuracy requires that the outputs of our model $f(x)$
and the outputs of the additive feature attribution method to be equal. In
symbols, this means that

$f(x)=g(x^{\prime})=\phi_{0}+\sum^{M}_{i=1}\phi_{i}x^{\prime}_{i}.$

Missingness. A second property is missingness. Simply stated, any variable
that has a value of 0 requires its corresponding contribution to the output to
be zero. In other words,

$x^{\prime}_{i}=0\Rightarrow\phi_{i}=0$

Consistency. The third property is consistency, which assures that if the
model changes so that an input’s contribution increases or stays the same, the
attribution of that input should not decrease. If we let
$f_{x}(z^{\prime})=f(h_{x}(z^{\prime}))$ and $z^{\prime}/i$ denote setting
$z_{i}^{\prime}=0$, then for any two models $f$ and $f^{\prime}$, if

$f^{\prime}_{x}(z^{\prime})-f^{\prime}_{x}(z^{\prime}/i)\geq
f_{x}(z^{\prime})-f_{x}(z^{\prime}/i)$

for all inputs $z^{\prime}\in\\{0,1\\}^{M}$, then
$\phi_{i}(f^{\prime},x)\geq\phi_{i}(f,x)$.

The theorem proposed by Lundberg and Lee is as follows: Only one possible
explanation model follows the above definition and the three given properties.

$\phi_{i}(f,x)=\sum_{z^{\prime}\subseteq
x^{\prime}}\frac{|z^{\prime}|!(M-|z^{\prime}|-1)!}{M!}[f_{x}(z^{\prime})-f_{x}(z^{\prime}/i)]$

where $|z^{\prime}|$ is the number of non-zero entries in $z^{\prime}$ and
$z^{\prime}\subseteq x^{\prime}$ represents all $z^{\prime}$ vectors where the
non-zero entries are asubset of the non-zero entries in $x^{\prime}$.

##  Appendix B The Relationship Between $\mu_{f},\mu_{g}$ and $\mu_{h}$

Recall from above that

$\mu_{h}=\frac{1}{n}\sum_{i=1}^{n}\hat{z_{i}}=\frac{\hat{z_{1}}+\hat{z_{2}}+\hat{z_{3}}+\ldots+\hat{z_{n}}}{n},$

and that we defined model $h$ as the product of models $f$ and $g$. Thus, any
$\hat{z_{i}}$ is equivalent to $\hat{x_{i}}\hat{y_{i}}$.

Taking the equation above and substituting $\hat{x_{i}}\hat{y_{i}}$ for every
$\hat{z_{i}}$ we see that

$\mu_{h}=\frac{\hat{x_{1}}\hat{y_{1}}+\hat{x_{2}}\hat{y_{2}}+\hat{x_{3}}\hat{y_{3}}+\ldots+\hat{x_{n}}\hat{y_{n}}}{n}.$

Whenever we multiply $\hat{x_{i}}$ and $\hat{y_{i}}$ to get $\hat{z_{i}}$, it
is inevitable that we end up with the term $\mu_{f}\mu_{g}$ in the resulting
expansion. We will take this term and split it into two parts: $\mu_{h}$ and
$\alpha$, some correction that must be added in to the other SHAP values.
Start with the expansion of $\mu_{f}\mu_{g}$,

$\mu_{f}\mu_{g}=(\frac{1}{n}\sum_{i=1}^{n}\hat{x_{i}})\cdot(\frac{1}{n}\sum_{i=1}^{n}\hat{y_{i}})$

which can be written in tabular form for ease of explanation.

$\begin{array}[]{c|ccccccccc}&\frac{\hat{x_{1}}}{n}&+&\frac{\hat{x_{2}}}{n}&+&\frac{\hat{x_{3}}}{n}&+&\ldots&+&\frac{\hat{x_{n}}}{n}\\\\[6.0pt]
\hline\cr\\\
\frac{\hat{y_{1}}}{n}&\frac{\hat{x_{1}}\hat{y_{1}}}{n^{2}}&&\frac{\hat{x_{2}}\hat{y_{1}}}{n^{2}}&&\frac{\hat{x_{3}}\hat{y_{1}}}{n^{2}}&&\ldots&&\frac{\hat{x_{n}}\hat{y_{1}}}{n^{2}}\\\
+&&&&&&&&&\\\
\frac{\hat{y_{2}}}{n}&\frac{\hat{x_{1}}\hat{y_{2}}}{n^{2}}&&\frac{\hat{x_{2}}\hat{y_{2}}}{n^{2}}&&\frac{\hat{x_{3}}\hat{y_{2}}}{n^{2}}&&\ldots&&\frac{\hat{x_{n}}\hat{y_{2}}}{n^{2}}\\\
+&&&&&&&&&\\\
\frac{\hat{y_{3}}}{n}&\frac{\hat{x_{1}}\hat{y_{3}}}{n^{2}}&&\frac{\hat{x_{2}}\hat{y_{3}}}{n^{2}}&&\frac{\hat{x_{3}}\hat{y_{3}}}{n^{2}}&&\ldots&&\frac{\hat{x_{n}}\hat{y_{3}}}{n^{2}}\\\
+&&&&&&&&&\\\ \vdots&\vdots&&\vdots&&\vdots&&\ddots&&\vdots\\\ +&&&&&&&&&\\\
\frac{\hat{y_{n}}}{n}&\frac{\hat{x_{1}}\hat{y_{n}}}{n^{2}}&&\frac{\hat{x_{2}}\hat{y_{n}}}{n^{2}}&&\frac{\hat{x_{3}}\hat{y_{n}}}{n^{2}}&&\ldots&&\frac{\hat{x_{n}}\hat{y_{n}}}{n^{2}}\\\
\end{array}$

Note that along the diagonal are the terms that may be of interest to us,
specifically

$\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n^{2}}=\frac{\mu_{h}}{n}.$

By multiplying both sides by $n$, we see that

$n\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n^{2}}=\mu_{h}.$

And since we already have one, we can simply add $n-1$ and subtract $n-1$
summands to pull out the desired $\mu_{h}$. This can be summarized as follows

$\displaystyle\mu_{f}\mu_{g}$
$\displaystyle=\sum^{n}_{i=1}\sum^{n}_{j=1}(\frac{\hat{x_{i}}\hat{y_{j}}}{n^{2}}I(i\neq
j))-(n-1)\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n^{2}}+\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n}$
$\displaystyle=\sum^{n}_{i=1}\sum^{n}_{j=1}(\frac{\hat{x_{i}}\hat{y_{j}}}{n^{2}}I(i\neq
j))-(n-1)\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n^{2}}+\mu_{h}$
$\displaystyle=\alpha+\mu_{h}$

where
$\alpha=\sum^{n}_{i=1}\sum^{n}_{j=1}(\frac{\hat{x_{i}}\hat{y_{j}}}{n^{2}}I(i\neq
j))-(n-1)\sum^{n}_{i=1}\frac{\hat{x_{i}}\hat{y_{i}}}{n^{2}}=\mu_{f}\mu_{g}-\mu_{h}$.
This becomes a critical element in our substitutions in later steps.

###  B.1 Proof of Local Accuracy

If we define $\hat{z_{i}}$ as the prediction of our model, $h$ for the $i$th
observation, $\mu_{h}$ as the average model prediction across our training
set, and $s_{z_{i}j}$ as the contribution of the $j$th variable to the $i$th
observation’s prediction, we can define local accuracy as

$\hat{z_{i}}=\mu_{h}+\sum^{p}_{j=1}s_{z_{i}j}.$

In this section we will prove that this equation holds for our chosen
definition of $s_{z_{i}j}$.

Remember that based on our initial definition,
$\hat{z_{i}}=\hat{x_{i}}\hat{y_{i}}$, and recall from above the final equation
for the mSHAP values:

$s_{z_{i}j}=\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]+\frac{|s_{z_{i}j}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha).$

We see that

$\displaystyle\mu_{h}+\sum^{p}_{j=1}s_{z_{i}j}$
$\displaystyle=\mu_{h}+\sum^{p}_{j=1}\left(\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]+\frac{|s_{z_{i}j}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha)\right)$
$\displaystyle=\mu_{h}+\left(\mu_{f}s_{y_{i}1}+s_{x_{i}1}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}1}s_{y_{i}a}+s_{y_{i}1}s_{x_{i}a})\right]+\frac{|s_{z_{i}1}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha)\right)+\ldots$
$\displaystyle\quad\quad\quad\ldots+\left(\mu_{f}s_{y_{i}p}+s_{x_{i}p}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}p}s_{y_{i}a}+s_{y_{i}p}s_{x_{i}a})\right]+\frac{|s_{z_{i}p}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha)\right)$
$\displaystyle=\mu_{h}+\frac{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}{\sum^{p}_{k=1}|s_{z_{i}k}^{\prime}|}(\alpha)+\sum^{p}_{j=1}\left(\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]\right)$
$\displaystyle=\mu_{h}+\alpha+\sum^{p}_{j=1}\left(\mu_{f}s_{y_{i}j}+s_{x_{i}j}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a})\right]\right)$

At this point we recall the definition given in Section 3.1 that
$\mu_{f}\mu_{g}-\mu_{h}=\alpha$. With a simple manipulation, we see that
$\mu_{h}+\alpha=\mu_{f}\mu_{g}$. Thus,

$\displaystyle=\mu_{f}\mu_{g}+\left(\mu_{f}s_{y_{i}1}+s_{x_{i}1}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}1}s_{y_{i}a}+s_{y_{i}1}s_{x_{i}a})\right]\right)+\ldots$
$\displaystyle\quad\quad\quad\ldots+\left(\mu_{f}s_{y_{i}p}+s_{x_{i}p}\mu_{g}+\frac{1}{2}\left[\sum^{p}_{a=1}(s_{x_{i}p}s_{y_{i}a}+s_{y_{i}p}s_{x_{i}a})\right]\right)$
$\displaystyle=\mu_{f}\mu_{g}+\sum^{p}_{j=1}\mu_{f}s_{y_{i}j}+\sum^{p}_{j=1}s_{x_{i}j}\mu_{g}+\frac{1}{2}\sum^{p}_{j=1}\sum^{p}_{a=1}(s_{x_{i}j}s_{y_{i}a}+s_{y_{i}j}s_{x_{i}a}).$

We can expand this further to give us

$\displaystyle=\mu_{f}\mu_{g}+\sum^{p}_{j=1}\mu_{f}s_{y_{i}j}+\sum^{p}_{j=1}s_{x_{i}j}\mu_{g}+\frac{1}{2}(s_{x_{i}1}s_{y_{i}1}+s_{y_{i}1}s_{x_{i}1}+s_{x_{i}1}s_{y_{i}2}+s_{y_{i}1}s_{x_{i}2}+\ldots+s_{x_{i}1}s_{y_{i}p}$
$\displaystyle\quad\quad\quad+s_{y_{i}1}s_{x_{i}p}+s_{x_{i}2}s_{y_{i}1}+s_{y_{i}2}s_{x_{i}1}+s_{x_{i}2}s_{y_{i}2}+s_{y_{i}2}s_{x_{i}2}+\ldots+s_{x_{i}2}s_{y_{i}p}+s_{y_{i}2}s_{x_{i}p}$
$\displaystyle\quad\quad\quad+\ldots+\ldots$
$\displaystyle\quad\quad\quad+s_{x_{i}p}s_{y_{i}1}+s_{y_{i}p}s_{x_{i}1}+s_{x_{i}p}s_{y_{i}2}+s_{y_{i}p}s_{x_{i}2}+\ldots+s_{x_{i}p}s_{y_{i}p}+s_{y_{i}p}s_{x_{i}p})$
$\displaystyle=\mu_{f}\mu_{g}+\sum^{p}_{j=1}\mu_{f}s_{y_{i}j}+\sum^{p}_{j=1}s_{x_{i}j}\mu_{g}+\frac{1}{2}(2s_{x_{i}1}s_{y_{i}1}+2s_{x_{i}1}s_{y_{i}2}+2s_{x_{i}1}s_{y_{i}3}+\ldots+2s_{x_{i}1}s_{y_{i}p}$
$\displaystyle\quad\quad\quad+2s_{x_{i}2}s_{y_{i}1}+2s_{x_{i}2}s_{y_{i}2}+2s_{x_{i}2}s_{y_{i}3}+\ldots+2s_{x_{i}2}s_{y_{i}p}$
$\displaystyle\quad\quad\quad+\ldots+\ldots$
$\displaystyle\quad\quad\quad+2s_{x_{i}p}s_{y_{i}1}+2s_{x_{i}p}s_{y_{i}2}+2s_{x_{i}p}s_{y_{i}3}+\ldots+2s_{x_{i}p}s_{y_{i}p}$
$\displaystyle=(\mu_{f}+s_{x_{i}1}+s_{x_{i}2}+\ldots+s_{x_{i}p})(\mu_{g}+s_{y_{i}1}+s_{y_{i}2}+\ldots+s_{y_{i}p}).$

Since the original SHAP values have the local accuracy property, we know that

$(\mu_{f}+s_{x_{i}1}+s_{x_{i}2}+\ldots+s_{x_{i}p})(\mu_{g}+s_{y_{i}1}+s_{y_{i}2}+\ldots+s_{y_{i}p})=\hat{x_{i}}\hat{y_{i}}$

Which in turn is equal to $\hat{z_{i}}$. We see that
$\hat{z_{i}}=\mu_{h}+\sum^{p}_{j=1}s_{z_{i}j}$ and that the local accuracy
property holds for mSHAP.

##  Appendix C The Simulation

###  C.1 Simulation Process

The basic flow for the simulation involved creating a data frame with all our
desired combinations of $y_{1}$, $y_{1}$, $\theta_{1}$ and $\theta_{2}$ and
then mapping through the following steps for each row:

  1. 1.

Using randomly distributed data as the covariates, create the response
variables by evaluating $y_{1}$ and $y_{2}$ and then multiplying them
together.

  2. 2.

Create two gradient boosted forests, one to predict $y_{1}$ and the other to
predict $y_{2}$, based on the covariates.

  3. 3.

Multiply the model predictions together and run kernelSHAP to approximate
explanations for the final model output.

  4. 4.

Use TreeSHAP to obtain exact explanations for the predictions of $y_{1}$ and
$y_{2}$.

  5. 5.

Multiply the TreeSHAP values together, using the method described in Section 3
to calculate mSHAP values for each variable.

  6. 6.

Distribute $\alpha$ into the subsequent mSHAP values, in each of the four
proposed ways.

  7. 7.

Compare the mSHAP values to the kernelSHAP values, using the scoring metrics
described in Section 4.1.

  8. 8.

Record the resulting scores in a data frame.

As previously mentioned, final scores were calculated by taking the average
across all variables and all combinations of the inputs. The code used to
perform the simulation can be found in the github repo at
https://github.com/srmatth/CAS, inside the mSHAP directory.

###  C.2 Additional Simulations

Since the initial simulation only used data with three explanatory variables,
we have completed additional simulations with different numbers of variables.
The goal of this is to ascertain that the weighted by absolute value is the
best method no matter the number of variables.

Our additional simulations used between 10 and 50 covariates across over 250
combinations of $y_{1}$, $y_{2}$, $\theta_{1}$, and $\theta_{2}$. For these
simulations all of our covariates were distributed uniformly between -1 and 1.
After performing the simulation, we saw that the absolute value method of
weighting alpha is again the best (but just barely) based on overall score and
in other metrics as well.

Method | Score |  Direction Score |  Relative Value Score | Rank Score |  Pct Same Sign |  Pct Same Rank
---|---|---|---|---|---|---
Weighted by Absolute Value | 2.13 | 0.884 | 0.770 | 0.480 | 74.4% | 24.9%
Uniformly Distributed | 2.13 | 0.890 | 0.766 | 0.470 | 75.0% | 23.7%
Weighted by Squared Value | 2.12 | 0.880 | 0.768 | 0.475 | 73.9% | 24.3%
Weighted by Raw Value | 2.00 | 0.780 | 0.753 | 0.468 | 63.5% | 23.2%

Due to these results, we are assured that the absolute weighting method of
distributing $\alpha$ is the best based on our chosen metrics, across
different numbers of covariates. It can be seen in Figure 8 that the general
score decreases as we add more variables. However, this is consistent with
what we see when we compare TreeSHAP (exact) to kernelSHAP (on singular
models, not two-part models), as demonstrated in Figure 9.

Figure 8: How the number of covariates impacts overall score, on average for
mSHAP compared to kernelSHAP Figure 9: How the number of covariates impacts
overall score, on average, for TreeSHAP compared to kernelSHAP

##  Appendix D The Data

The data used to create the model is a Property Damage data set, which is not
available publicly but can be obtained through the Casualty Actuarial Society.

##  Appendix E The Model

Both the severity model and the frequency model were tuned in R using a h2o
backend [6]. Tuning parameters are given in Table 1 and model metrics are
given in Table 2. Note that all model metrics were computed on the test (hold-
out) subset of data. These tuning results were then used to create the final
model in python using scikit-learn [11]. Note that scikit-learn was used to
create the models because multinomial predictions do not have SHAP support in
H2O as of the time of writing.

Tuning Parameter | Severity Model | Frequency Model
---|---|---
ntrees | 200 | 100
max_depth | 30 | 20
mtries | 20 | 20
min_split_improvement | .0001 | .001
sample_rate | 0.632 | 0.632
Table 1: Tuning parameters for the frequency and severity models Model | MAE | MSE | Logloss
---|---|---|---
Severity Model | 2,832 | 16,359,170 | NA
Frequency Model | NA | 0.074 | 0.427
Two-Part Model | 683 | 830,351 | NA
Table 2: Model metrics for all models

##  Appendix F Code Availability

The code used to tune the model (as well as additional code focused on working
with the CAS datasets) can be found at this github link:
https://github.com/srmatth/CAS.

mSHAP has been developed into an R package as well. The R package can be
downloaded from CRAN, with the R code

⬇

install.packages(”mshap”)

or get the development version from https://github.com/srmatth/mshap by
running

⬇

devtools::install_github(”srmatth/mshap”)

in R.

The mSHAP package repository (https://github.com/srmatth/mshap) also contains
all the code and data used to generate the plots in this paper, as well as the
code used to run the various simulations mentioned. It can be found in the
inst/paper directory under the main directory of the package. Be aware that
installing the package by following the steps above will not download the code
used in this paper, it must be obtained from the github repository.

Acknowledgments

  * 1.

Brigham Young University Department of Statistics Computing Cluster

  * 2.

Brian Fanin and the Casualty Actuarial Society for providing the data

  * 3.

Isabelle Matthews for proof-reading

Funding: This paper was funded by an individual grant from the Casualty
Actuarial Society.

## References

  * Akinyemi and Leiser [2020] Akinyemi, K., Leiser, B., 2020\.  The use of advanced predictive analytics for rate making in insurance .
  * Arrieta et al. [2020] Arrieta, A.B., Díaz-Rodríguez, N., Del Ser, J., Bennetot, A., Tabik, S., Barbado, A., García, S., Gil-López, S., Molina, D., Benjamins, R., et al., 2020\.  Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai.  Information Fusion 58, 82–115.
  * Doran et al. [2017] Doran, D., Schulz, S., Besold, T.R., 2017.  What does explainable ai really mean? a new conceptualization of perspectives.  arXiv preprint arXiv:1710.00794 .
  * Frees and Sun [2010] Frees, E.W., Sun, Y., 2010.  Household life insurance demand: A multivariate two-part model.  North American Actuarial Journal 14, 338–354.
  * Gunning [2017] Gunning, D., 2017.  Explainable artificial intelligence (xai).  Defense Advanced Research Projects Agency (DARPA), nd Web 2, 2.
  * H2O.ai [2021] H2O.ai, 2021.  h2o R Package.  URL: https://github.com/h2oai/h2o-3. version 3.34.0.1.
  * Heras et al. [2018] Heras, A., Moreno, I., Vilar-Zanón, J.L., 2018.  An application of two-stage quantile regression to insurance ratemaking.  Scandinavian Actuarial Journal 2018, 753–769.
  * Lipton [2018] Lipton, Z.C., 2018.  The mythos of model interpretability: In machine learning, the concept of interpretability is both important and slippery.  Queue 16, 31–57.
  * Lundberg et al. [2020] Lundberg, S.M., Erion, G., Chen, H., DeGrave, A., Prutkin, J.M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., Lee, S.I., 2020\.  From local explanations to global understanding with explainable ai for trees.  Nature machine intelligence 2, 2522–5839.
  * Lundberg and Lee [2017] Lundberg, S.M., Lee, S.I., 2017\.  A unified approach to interpreting model predictions, in: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (Eds.), Advances in Neural Information Processing Systems 30. Curran Associates, Inc., pp. 4765–4774.  URL: http://papers.nips.cc/paper/7062-a-unified-approach-to-interpreting-model-predictions.pdf.
  * Pedregosa et al. [2011] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Shapley [1953] Shapley, L.S., 1953.  A value for n-person games.  Contributions to the Theory of Games 2, 307–317.
  * slundberg [2020] slundberg, 2020.  SHAP Values for Ensemble of XGBoost Models.  https://github.com/slundberg/shap/issues/112.  Online; accessed 7 April 2021; First Post June 2018; Last post July 2020.