matrixdist: An R Package for Statistical Analysis of Matrix Distributions

Recently, a regression model for frequency modelling based on DPH distributions was introduced in Bladt and Yslas, 2023c as a generalisation of mixture-of-experts (MoE) models. The main idea is to incorporate the covariate information on the vector of initial probabilities, which play the role of expert weights. The resulting distributional model is flexible, as it possesses the property of denseness on more general regression models satisfying some mild technical conditions. Next, we present the fundamental concepts of this regression model. Define the mapping

where $\Delta^{p-1}=\{(\alpha_{1},\dots,\alpha_{p})\in\mathbb{R}^{p}\mid\sum_{k}\alpha_{k}=1,\;\alpha_{k}\geq 0\;\forall k\}$ is the standard $(p-1)$ simplex. Then, for any given vector of covariate information $\mathbf{X}=(X_{1},\dots,X_{h})\in\mathbb{R}^{h}$, the initial probabilities of the underlying MC of a DPH distribution are taken to be $\mathbb{P}(Z_{0}=k)=\alpha_{k}(\mathbf{x}):=(\bm{\alpha}(\mathbf{x}))_{k}$, $k=1,\dots,p$. Thus, we say that $N\mid\mathbf{X}\sim\mbox{DPH}(\bm{\alpha}(\mathbf{X}),\mathbf{S})$ follows the DPH-MoE specification. For $D=\mathbb{R}^{h}$, we have a convenient parametrisation of initial probabilities, namely the softmax parametrisation given by

where $\bm{\gamma}_{k}\in\bar{\mathbb{R}}^{h}$, and $\bm{\gamma}=(\bm{\gamma}_{1}^{\top},\dots,\bm{\gamma}_{p}^{\top})^{\top}\in\bar{\mathbb{R}}^{p\times h}$. More details about this model as well as precise denseness considerations can be found in Bladt and Yslas, 2023c . We also refer to Bladt and Yslas, 2023b for an extension of this regression model to continuous phase-type distributions.

At least two regression models based on IPH distributions exist in the literature. The first model, known as the proportional intensities (PI) model, was introduced in Albrecher et al., 2022a . The authors, inspired by the proportional hazards model (cf. Cox,, 1972), formulated a way of regressing on the intensity function $\lambda$ for different covariate information. More specifically, the main idea of the PI model is that the inhomogeneous intensity $\lambda(\cdot;\bm{\theta})$, which is assumed to be fully determined by a vector of parameters $\bm{\theta}$, is proportionally affected for different covariate information $\mathbf{X}\in\mathbb{R}^{h}$ according to a positive real-valued and measurable function $m(\cdot)$. In other words, the model assumes that

where $\bm{\gamma}$ is a $h$-dimensional column vector containing regression coefficients. With this assumption, $Y\mid\mathbf{X}\sim\mbox{IPH}(\bm{\alpha},\mathbf{S},\lambda(\cdot;\mathbf{X},\bm{\theta},\bm{\gamma}))$ with corresponding density and cumulative functions

and

Typically, the measurable function $m()$ is taken to be $m(x)=\exp(x)$, which is the default implementation in matrixdist. In this framework, when $p=1$, one recovers the proportional hazard model. However, for $p>1$, the implied hazard functions between different subgroups can deviate from proportionality in the distribution body but are asymptotically proportional in the tail. For more details on the PI model, we refer the reader to Albrecher et al., 2022a and Bladt, (2022).

The second regression model, called the phase-type mixture-of-experts (PH-MoE) model, was introduced in Bladt and Yslas, 2023b as a generalisation of mixture-of-experts models to the PH distributions framework and follows the same rationale as the DPH-MoE specification described above (cf. Bladt and Yslas, 2023b for further details).

Let $N\sim\mbox{DPH}_{p}(\bm{\alpha},\mathbf{S})$ be a DPH distributed random variable with underlying Markov chain $(Z_{n})_{n\in\mathbb{N}_{0}}$. Let $\mathbf{r}_{j}=(r_{j}(1),\dots,r_{j}(p))^{\top}$ be column vectors of dimension $p$ taking values in $\mathbb{N}_{0}^{p}$, $j=1,\dots,d$, and let $\mathbf{R}=(\mathbf{r}_{1},\dots,\mathbf{r}_{d})$ be a $p\times d$ matrix, called the reward matrix. We define

for all $j=1,\dots,d$. Then, we say that the random vector $\mathbf{N}=(N^{(1)},\dots,N^{(d)})$ has a multivariate discrete phase-type distribution of the MDPH* type, and we write $\mathbf{N}\sim\mbox{MDPH*}(\bm{\alpha},\mathbf{S},\mathbf{R})$. This class of distributions was introduced in Campillo-Navarro, (2018), and the construction principle can be seen as an extension of the so-called transformation via rewards for univariate DPH distributions. In this contribution, the author showed that elements in this class possess explicit expressions for joint probability-generating function, joint moment-generating function, and joint moments. Moreover, the class is dense in the class of distributions with support in $\mathbb{N}^{d}$, and margins are DPH distributed. However, general closed-form expressions for the joint density and distribution functions are not available, and the estimation of this general class is still an open question, limiting their use in practice. Nonetheless, there are subclasses of MDPH* distributions with explicit expressions of these two functionals that preserve the denseness property of the MDPH* class. The matrixdist package provides implementations for two of these subclasses, which we describe next.

Let $\mathbf{N}=(N^{(1)},N^{(2)})\sim\mbox{MDPH*}(\bm{\alpha},\mathbf{S},\mathbf{R})$ with

Here, $\mathbf{S}_{11}$ and $\mathbf{S}_{22}$ are sub-transition matrices of dimensions $p_{1}$ and $p_{2}$ ($p_{1}+p_{2}=p$), respectively, $\mathbf{S}_{12}$ is a $p_{1}\times p_{2}$ matrix satisfying $\mathbf{S}_{11}\mathbf{e}+\mathbf{S}_{12}\mathbf{e}=\mathbf{e}$, and $\bm{\pi}$ is a $p_{1}$ dimensional vector of initial probabilities. Then, one can show that the joint density of this random vector is given by

where $\mathbf{s}_{2}=\mathbf{e}-\mathbf{S}_{22}\mathbf{e}$. In this case, we say that $\mathbf{N}$ is bivariate discrete phase-type distributed of the feed-forward type, and we write $\mathbf{N}\sim\mbox{fMDPH}(\bm{\pi},\mathbf{S}_{11},\mathbf{S}_{12},\mathbf{S}_{22})$. In particular, the marginals are parametrised by $N^{(1)}\sim\mbox{DPH}(\bm{\pi},\mathbf{S}_{11})$ and $N^{(2)}\sim\mbox{DPH}(\bm{\pi}\left(\mathbf{I}-\mathbf{S}_{11}\right)^{-1}\mathbf{S}_{12},\mathbf{S}_{22})$. For more details on this class of bivariate DPH distributions, we refer the reader to Bladt and Yslas, 2023c . Although the above construction principle can be extended to higher dimensions, the implementation of the methods associated with this class becomes more challenging. Fortunately, another subclass of MDPH* distributions introduced in Bladt and Yslas, 2023c can easily be tracked and implemented in higher dimensions, as we describe next.

Let $(Z_{n}^{(j)})_{n\in\mathbb{N}_{0}}$, $j=1,\dots,d$, be Markov chains on a common state space with $p$ transient states. All chains are assumed to start in the same state and then evolve independently until respective absorptions. Formally,

If $N^{(j)}$, $j=1,\dots,d$, are the univariate DPH distributed absorption times of $(Z_{n}^{(j)})_{n\in\mathbb{N}_{0}}$, we say that the vector $\mathbf{N}=(N^{(1)},\ldots,N^{(d)})$ has a multivariate discrete phase-type distribution of the mDPH type, and we write

where $\bm{\alpha}$ is the common vector of initial probabilities and $\widetilde{S}=\{\mathbf{S}_{1},\dots,\mathbf{S}_{d}\}$, with $\mathbf{S}_{j}$ the sub-transition matrix associated with $(Z_{n}^{(j)})_{n\in\mathbb{N}_{0}}$, $j=1,\dots,d$. Explicit expressions for several functionals of interest of this class can easily be obtained by conditioning on the starting value of the Markov chains. For instance, the joint density function of $\mathbf{N}\sim\mbox{mDPH}(\bm{\alpha},\widetilde{S})$ is given by

We consider the AutoBi dataset from the insuranceData package, which contains information about automobile bodily injury claims. The dataset comprises eight variables: one for the economic loss (which we aim to describe), six covariates, and one variable to identify each claim. We start by providing general descriptive statistics for the variable of interest (LOSS).

# Load the data
data("AutoBi", package = "insuranceData")
claim <- AutoBi
# Economic losses resulting from an automobile accident
loss <- claim$LOSS
# Summary statistics of losses
summary(loss)

#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#>    0.005    0.640    2.331    5.954    3.995 1067.697

# 10 largest values of losses
tail(sort(loss), n = 10)

#>  [1]   82.000   96.007  114.604  150.000  162.047  188.720  193.000  222.405
#>  [9]  273.604 1067.697

We observe that the losses range from 0.005 to 1067.697, with 75% of the data smaller than 3.995. The presence of large values, compared to the body of the distribution, serves as an initial indication of a heavy tail. To further confirm this observation, we create a Hill plot using the hill() function from the evir package.

library(evir)
hill(loss)

The plot exhibits a flat behaviour for the largest order statistics, suggesting a Pareto-type tail (or regularly varying). Therefore, we opt for an IPH model to more accurately capture the tail of the losses. More specifically, we use a matrix-Pareto distribution. The fitting process starts by creating an initial IPH distribution, which provides the starting parameters for the EM algorithm. This can be accomplished by employing the iph() method in matrixdist. For our analysis, we create an IPH object with random parameters of a general structure, a dimension 6, and a Pareto inhomogeneity function.

set.seed(22)
# Create an IPH object
x <- iph(gfun = "pareto", structure = "general", dimension = 6)

With an initial IPH distribution in place, we can employ the fit() method to fit the data. It is important to note that fit() requires specifying the number of steps for the EM algorithm. In this instance, we use 1500 steps, as the changes in the loglikelihood become negligible with this number (the same criteria is employed to determine the number of iterations in the subsequent illustration). The method outputs an IPH object containing the fitted parameters. Further steps could be made on the output object if required.

# Fitting procedure
z <- fit(x, y = loss, stepsEM = 1500)

Figures 5.2 and 5.3 display a QQ plot for the fitted distribution and a histogram of the (log) data against the fitted density, respectively. These plots indicate that the fitted distribution serves as a good approximation of the data. In particular, the sample and fitted quantiles align closely with the identity line, especially within the distribution’s body. Although there is a slight deviation from the identity line for quantiles in the tail, they remain relatively close to the blue line. The histogram further supports the notion that the fitted IPH distribution is a satisfactory approximation.

We now shift our focus to the regression framework. In particular, we employ the PI specification to capture the influence of covariate information on economic losses. The initial step involves removing any incomplete covariate data. Next, we need to split the covariates and the response variable into two different objects. In our case, we investigate the effect of the variables ATTORNEY, CLMSEX and CLMAGE. The variable ATTORNEY indicates whether the claimant was represented by an attorney or not, CLMSEX is the claimant’s sex, and CLMAGE corresponds to their age.

# Covariate information, excluding missing entries
cov <- na.omit(claim)
# Separation of the information
X <- cov[, c("ATTORNEY", "CLMSEX", "CLMAGE")]
loss <- cov[, "LOSS"]

Now, we can employ the reg() method to perform the PI regression. This function works similarly to the fit() method, with the main difference being that estimation of the regression parameters is included in each iteration of the EM algorithm (see Albrecher et al., 2022a for details). Along with the data, the reg() method requires an initial IPH distribution. In this case, we employ the previously fitted matrix-Pareto distribution without covariates.

# PI regression
z.reg <- reg(x = z, y = loss, X = X, stepsEM = 250)

The resulting object (z.reg) contains all matrix parameters along with the regression coefficients. This information can be easily accessed via the coef() function as follows.

# Fitted parameters
z.reg.par <- coef(z.reg)
z.reg.par

#> $alpha
#> [1] 2.356094e-04 3.257751e-06 5.920828e-02 3.755407e-04 4.561202e-01
#> [6] 4.840571e-01
#>
#> $S
#>               [,1]          [,2]        [,3]          [,4]          [,5]
#> [1,] -3.412707e+01  7.592078e-06  1.42774265  1.533919e-01  6.975181e-06
#> [2,]  2.145191e+01 -3.473845e+01  0.38024154  1.288647e+01  3.162116e-05
#> [3,]  3.090083e-01  6.287658e-01 -1.83234510  4.643455e-01  3.670351e-01
#> [4,]  5.666353e-01  1.166288e-04  7.77223071 -2.987094e+01  2.335434e-04
#> [5,]  1.432727e-04  3.578715e+01  0.49172745  2.247409e-04 -3.628114e+01
#> [6,]  4.561651e-06  1.243609e-03  0.08311972  3.771681e-05  4.707313e-02
#>               [,6]
#> [1,]  3.037507e+01
#> [2,]  2.799554e-05
#> [3,]  3.352413e-02
#> [4,]  2.066459e+01
#> [5,]  1.118600e-03
#> [6,] -6.641175e+00
#>
#> $beta
#> [1] 33.10073
#>
#> $B
#>
#>  1.15139949 -0.10888748 -0.01368395
#>
#> $C
#> numeric(0)

The above information can be used, for example, to create customised loss distributions for individual claimants based on their distinctive characteristics.

Take, for instance, a new claimant: a 60-year-old male, unrepresented by an attorney. The following code leverages our model to estimate his loss distribution.

# Covariate information of the new claimant
claimant <- data.frame(ATTORNEY = 2, CLMSEX = 1, CLMAGE = 60)
# Multiplicative effect in the sub-intensity matrix
prop <- exp(sum(z.reg.par$B * claimant))
# IPH distribution for claimant’s loss
claimant.iph <- iph(
  alpha = z.reg.par$alpha, S = z.reg.par$S * prop,
  gfun = "pareto", gfun_pars = z.reg.par$beta
)

If we wish to visualize this claimant’s survival function, we can plot it as follows.

# Evaluation points for the survival function
q <- seq(0, 15, 0.01)
# Plot of the survival function
plot(q, cdf(claimant.iph, q, lower.tail = F),
  type = "l", col = "blue", lwd = 1.3,
  ylab = "1-F(y)", xlab = "y", main = "Survival function"
)

We refer readers to Section 6 and the documentation of matrixdist for further available methods for PH and IPH distributions.

The matrixdist package provides an efficient approach to working with DPH distributions. To exemplify its application, we will use the NMES1988 data set included in the R package AER. This data is the result of a US National Medical Expenditure Survey that took place from 1987 to 1988 and contains various count variables. For our illustration, we will focus on the annual number of visits to a physician’s office (visits variable) for US residents aged above 66. Note, however, that when working with DPH distributions, it is necessary to add 1 to data starting at 0, as DPH random variable’s support commences at 1.

# Load the data
data("NMES1988", package = "AER")
df <- NMES1988
# Yearly physician office visits
visit <- df$visits + 1
df$y <- visit

The process of fitting a DPH distribution to data closely parallels that of fitting PH distributions. It starts by defining an initial dph object, which must then be passed to the fit() function along with the data (Figure 5.5 depicts the fitted distribution).

# Create a DPH distribution
set.seed(23)
d <- dph(structure = "general", dimension = 3)
# Fit the distribution to the data we are interested in
d.fit <- fit(x = d, y = visit, stepsEM = 1500)

The data set also provides many other covariates. Here, we will concentrate on examining the influence of age, gender, and number of chronic conditions on the frequency of visits to a physician. For such a purpose, we will employ the MoE framework to integrate the covariate information into the DPH specification. This can be done by employing the MoE() method. The first step to performing this regression method is to specify a formula (f) that details how the covariates affect the number of physician visits.

# Convert the gender variable into binary format
df$gender <- ifelse(df$gender == "male", 1, 0)
# Regression formula
f <- y ~ age + gender + chronic

Next, such a formula has to be passed on to the MoE() method along with the data (including covariates) and an initial DPH distribution.

# MoE regression
d.reg <- MoE(x = d.fit, formula = f, data = df, stepsEM = 500)

The fitted model allows us to examine the impact of different covariate profiles on the number of physician visits. For example, let us compare the probability mass functions for two subjects. Both are 69 years old, female, but they differ in the number of chronic conditions they possess: one has five chronic conditions while the other does not have any.

# DPH distribution for a subject with 5 chronic conditions
dph.sub1 <- dph(alpha = d.reg$alpha[12, ], S = d.reg$S)
# DPH distribution for a subject with 0 chronic conditions
dph.sub2 <- dph(alpha = d.reg$alpha[38, ], S = d.reg$S)
# Plot of density function for both subjects
q <- seq(1, 20, 1)
plot(q, dens(dph.sub1, q),
  type = "p", col = "blue", cex = 0.8,
  ylab = "P(N=n)", xlab = "n", ylim = c(0, 0.3)
)
points(q, dens(dph.sub2, q), col = "red", cex = 0.8)
legend("topright",
  legend = c("5 chronic conditions", "0 chronic conditions"),
  col = c("blue", "red"), bty = "n", pch = 1, cex = 0.8
)

Furthermore, we can quantify the impact of chronic conditions on the expected number of visits to a physician’s office. Let us compute the expected number of visits for both individuals:

# Mean value for subject with 5 chronic conditions
mean(dph.sub1) - 1

#> [1] 8.098524

# Mean value for subject with 0 chronic conditions
mean(dph.sub2) - 1

#> [1] 3.656367

Our analysis unsurprisingly reveals that the subject with five chronic conditions is expected to visit more often than the one with none. This demonstrates the ability of our fitted model to quantify the impact of varying covariate profiles - a tool which is essential for healthcare professionals when predicting healthcare usage patterns.

To illustrate the capabilities of the implementations for MPH* distributions in matrixdist, we use the rdj data set of the copula package. More specifically, our focus is to use an MPH* for describing the joint behaviour of the absolute log returns of the three stock prices in this data set: Intel, Microsoft, and General Electric.

Firstly, we load the data, take the absolute value, and eliminate any data points with atoms at zero.

set.seed(24)
# Load the data
data("rdj", package = "copula")
# Absolute log return of stock prices
y <- as.matrix(abs(rdj[, -1]))
# Remove data points with atoms at zero
y <- y[y[, 1] != 0 & y[, 2] != 0 & y[, 3] != 0, ]

As in the univariate case, the fitting process requires the definition of an initial MPHstar object that is subsequently fed into the fit() function. In the present case, we use an MPHstar object with randomly generated parameters, a general form of the matrix parameters, and dimension $10$.

# Initial MPH* distribution
x <- MPHstar(structure = "general", dimension = 10, variables = 3)
# Fitting to data
x.fit <- fit(x = x, y = y, stepsEM = 200, uni_epsilon = 1e-6, r = 0.5)

In particular, the reward matrix of the fitted distribution can be accessed as follows:

# Estimated reward matrix
x.fit@pars$R

#>              [,1]       [,2]       [,3]
#>  [1,] 0.624656883 0.03162784 0.34371527
#>  [2,] 0.448807782 0.08398967 0.46720255
#>  [3,] 0.545561519 0.42181735 0.03262113
#>  [4,] 0.783031036 0.08518594 0.13178302
#>  [5,] 0.973586315 0.01054082 0.01587287
#>  [6,] 0.318181076 0.21878855 0.46303038
#>  [7,] 0.015947222 0.78138202 0.20267076
#>  [8,] 0.498338747 0.28584901 0.21581224
#>  [9,] 0.001025196 0.99897480 0.00000000
#> [10,] 0.471608399 0.48276616 0.04562544

Finally, we assess the quality of the fit by creating QQ plots of the fitted margins.

# Names for the QQ plot
name <- colnames(y)
# Quantiles’ levels
qq <- seq(0.01, 0.99, length.out = 100)
par(mfrow = c(1, 3))
for (j in 1:3) {
  # j-th marginal fitted distribution
  x <- marginal(x.fit, j)
  # Marginal observations
  ym <- y[, j]
  # quantiles
  y.q <- quantile(ym, probs = qq)
  x.q <- quan(x, qq)
  # QQ plot
  plot(y.q, x.q, main = name[j], xlab = "Sample", ylab = "Fitted", pch = 16, asp = 1)
  abline(0, 1, col = "blue")
}

Upon inspecting the QQ plots, we can observe that the fitted marginals provide a satisfactory approximation of the data.

Finally, let us illustrate the use of the mDPH class with a practical example. We will use the NMES1988 data set again, but this time, we will focus on jointly modelling two distinct variables: the yearly number of physician office visits (visits) and yearly non-physician office visits (nvisits). With the matrixdist package, the fitting of a mDPH distribution to this data can be done as follows.

# Load the data
data("NMES1988", package = "AER")
df <- NMES1988
# Observations in matrix form
mvisit <- cbind(df$visits, df$nvisits) + 1
# Initial mDPH distribution
set.seed(12345)
md <- mdph(structure = rep("general", 2), dimension = 4, variables = 2)
# Fitting procedure
md.fit <- fit(x = md, y = mvisit, stepsEM = 1500)

The expected values, variance-covariance matrix, and correlation matrix of the fitted distribution can be obtained using the subsequent code:

# Vector of marginal expectations
mean(md.fit) - 1

#> [1] 5.774399 1.618021

# Variance-covariance matrix
var(md.fit)

#>           [,1]      [,2]
#> [1,] 46.243751  5.726782
#> [2,]  5.726782 27.372432

# Correlation matrix
cor(md.fit)

#>           [,1]      [,2]
#> [1,] 1.0000000 0.1609635
#> [2,] 0.1609635 1.0000000

Note that the mean of the fitted mDPH distribution matches the sample’s mean, a property that also holds for the estimation methods for the PH, DPH, and mPH classes. Additionally, the fitted model captures the correlation between physician office visits and non-office visits. The mDPH class also supports the MoE framework, allowing us to incorporate covariate information. Next, we incorporate age, gender, and the number of chronic conditions into our analysis.

df$y <- mvisit
# Convert the gender variable into binary format
df$gender <- ifelse(df$gender == "male", 1, 0)
# Regression formula
f <- y ~ age + gender + chronic
# MoE regression
md.reg <- MoE(x = md.fit, y = mvisit, formula = f, data = df, stepsEM = 500)

The above output allows us to predict the joint distribution for new individuals with a set of specific characteristics. This is achieved by employing the predict() function in conjunction with the nnet object included in md.reg (accessible via md.reg$mm). For instance, for a 72-year-old male patient with one chronic condition, the mDPH distribution describing the joint behaviour of the two variables of interest can be specified as follows:

# Characteristics of a new patient
patient <- data.frame(age = 7.2, gender = 1, chronic = 1)
# Predicted starting probabilities
alpha <- predict(md.reg$mm, type = "probs", newdata = patient)
# New mDPH distribution
md.patient <- mdph(alpha = alpha, S = md.reg$S)

Figure 5.7 shows the joint probability mass function of yearly office and non-office visits for the new patient.