Tensor product P-splines using a sparse mixed model formulation

Abstract

A new approach to represent P-splines as a mixed model is presented. The corresponding matrices are sparse allowing the new approach can find the optimal values of the penalty parameters in a computationally efficient manner. Whereas the new mixed model P-splines formulation is similar to the original P-splines, a key difference is that the fixed effects are modelled explicitly, and extra constraints are added to the random part of the model. An important feature ensuring that the entire computation is fast is a sparse implementation of the Automated Differentiation of the Cholesky algorithm. It is shown by means of two examples that the new approach is fast compared to existing methods. The methodology has been implemented in the R-package LMMsolver available on CRAN (https://CRAN.R-project.org/package=LMMsolver).

1 Introduction

Penalized regression using B-splines (known as P-splines; Eilers and Marx, 1996, 2021) is computationally efficient, because of the local character of B-splines. The corresponding linear equations are sparse and can be solved quickly. However, the main problem is to find the optimal value for the tuning or penalty parameter. A good way to approach this problem is to use mixed models and restricted maximum likelihood (REML; Patterson and Thompson, 1971). Several methods have been proposed to transform the original P-spline model to a mixed model (Eilers, 1999; Currie and Durban, 2002; Lee and Durbán, 2011). A disadvantage of existing transformations to mixed models is that the local character of the B-splines is lost, which reduces computational efficiency.

The proposed transformations have in common that they start from the P-splines formulation and decompose into a fixed and a random term to obtain a mixed model. The same idea can be used for tensor product P-splines (Currie et al., 2006; Lee and Durbán, 2011; Lee et al., 2013; Rodríguez-Ávarez et al., 2015; 2018; Wood et al., 2013; Wood, 2017). There are many different ways to represent tensor product P-splines as a mixed model; Piepho et al. (2022) give an overview of the different two-dimensional P-splines models.

Here, a new and efficient approach is outlined. In Boer et al., (2020) it was shown that the linear variance (LV) mixed model (Williams, 1986) is equivalent to P-splines with first-degree B-splines and first-order differences. This idea will be generalized here, to show that P-splines with $k$ th-order difference penalties are equivalent to a class of mixed models. The extension to tensor product P-splines is relatively straightforward.

Another property of the LV model is that it is computationally attractive, as the precision matrix is sparse. It will be shown that the P-splines model is equivalent to a class of mixed models with a sparse precision matrix. In this article we will show that the sparse mixed model P-splines are computationally efficient.

The layout of the article is as follows. In Section 2 we will discuss in detail the equivalence between the linear variance mixed model and a special case of P-splines. Section 3 outlines a brief overview of B-splines, and gives some results that will be used to derive the connection between mixed models and P-splines. In Section 4 it will be shown how the idea presented in Section 2 can be generalized to P-splines. Section 5 shows how this can be further extended to mixed model tensor product P-splines. In Section 6 an implementation of the sparse mixed model in the LMMsolver R-package is presented, and its computational efficiency is demonstrated by means of two examples. The article concludes with a brief discussion in Section 7. The R-scripts to reproduce the figures are available in the Supplementary Material.

2 Linear variance model and Whittaker smoother

The Whittaker smoother (Whittaker, 1923) and the Linear Variance (LV) mixed model introduced by Williams (1986) are equivalent, as shown by Boer et al. (2020). Here we will use a new approach to show that the two models are equivalent with the aim to develop a computational efficient method to estimate the optimal penalty parameters using mixed models that can be further extended to P-splines. The Whittaker smoother is a special case of P-splines (Frasso and Eilers, 2015), and originally third-order differences were used by Whittaker (1923). Here we will use a simpler variant of Whittakers proposal, using first-order differences.

Suppose there are $n$ observations at equidistant points $i=\mathrm{1,2,}\dots \mathrm{,}n$ . Let $y_i$ be the observation at point $i$ . The LV mixed model is given by

$$y_i=\mu +u_i+e_i\mathrm{,}i=1,\dots \mathrm{,}n\mathrm{,}$$

(2.1)

where $\mu$ is the intercept, the parameters $e_i$ represent the residual error and are assumed to be independent and identically distributed, $e_i\sim N(\mathrm{0,}{\sigma }_e^2)$ , where ${\sigma }_e^2$ is the residual variance. The covariances between $u_i$ and $u_j$ are given by

$$\text{cov}(u_i\mathrm{,}{u}_j)=\alpha (n-1-|i-j|){\sigma }_u^2\mathrm{,}$$

(2.2)

where $\alpha$ is a scaling factor, and ${\sigma }_u^2$ is the variance corresponding to the LV model. Whereas Williams (1986) used $\alpha =1$ , here we will use $\alpha =\frac{1}{2}$ , to simplify the notation to show the equivalence with the Whittaker smoother. As can be seen from equation (2.2) the covariance decreases linearly with distance between two observation points $i$ and $j$ , and the covariance between first and last observations is equal to zero.

In this article we will formulate mixed models in terms of penalty and precision matrix (the inverse of the covariance matrix), which is a more modern representation (Lindgren et al., 2022). The parameter $\lambda ={\sigma }_e^2/{\sigma }_u^2$ is called the tuning or penalty parameter (Eilers and Marx, 1996; Hastie et al., 2009). From equations (7) and (A2) in Boer et al. (2020) it follows that the objective function to be minized for equations (2.1) with respect to parameters $\mu$ and $u_i$ is given by

$$S_1(\mu \mathrm{,}{u}_1\mathrm{,}\dots \mathrm{,}{u}_n)={\displaystyle \sum _{i=1}^n}{\left(y_i-(\mu +u_i)\right)}^2+\text{λ}{\displaystyle \sum _{i=2}^n}{\left(\Delta u_i\right)}^2+\text{λ}{\left({\displaystyle \sum _{i=1}^n}{c}_i{u}_i\right)}^2\mathrm{,}$$

(2.3)

where $\Delta u_i=u_i-u_{i-1}$ is the first-order difference, $c_1=c_n=1/\sqrt{n-1}$ and $c_i=0$ for $i=\mathrm{2,}\dots \mathrm{,}n-1$ .

To make the connection between the LV model and the Whittaker smoother, we define the following invertible linear transformation:

$$\begin{array}{ll}{a}_i\hfill & =\mu +u_i\mathrm{,}i=\mathrm{1,2,}\dots \mathrm{,}n\mathrm{,}\hfill \\ \eta \hfill & ={\displaystyle \sum _{i=1}^n}{c}_i{u}_i\mathrm{,}{\displaystyle \sum _{i=1}^n}{c}_i\ne 0.\hfill \end{array}$$

(2.4)

Using equation (2.4) and

$$\Delta a_i=\Delta (\mu +u_i)=\Delta u_i$$

(2.5)

we obtain

$$S_2(a_1\mathrm{,}\dots \mathrm{,}{a}_n\mathrm{,}\eta )={\displaystyle \sum _{i=1}^n}{\left(y_i-a_i\right)}^2+\text{λ}{\displaystyle \sum _{i=2}^n}{\left(\Delta a_i\right)}^2+\text{λ}{\eta }^2\mathrm{.}$$

(2.6)

Equation (2.6) can be decomposed as $S_2(a_1\mathrm{,}\dots \mathrm{,}{a}_n\mathrm{,}\eta )=S_3(a_1\mathrm{,}\dots \mathrm{,}{a}_n)+S_4(\eta )$ where

$$S_3(a_1\mathrm{,}\dots \mathrm{,}{a}_n)={\displaystyle \sum _{i=1}^n}{\left(y_i-a_i\right)}^2+\text{λ}{\displaystyle \sum _{i=2}^n}{\left(\Delta a_i\right)}^2$$

(2.7)

and $S_4(\eta )=\text{λ}{\eta }^2$ . To find the minimum of $S_2$ the functions $S_3$ and $S_4$ can be minimized separately. The objective function $S_4$ has trivial solution $\widehat{\eta }=0$ . Equation (2.7) is the objective function to be minimized for the Whittaker smoother with first-order differences. This shows that equations (2.3) and (2.7) are equivalent. The LV model is not the only mixed model equivalent to the Whittaker smoother: other values of $c_i$ can also be chosen, provided that ${\displaystyle {\sum }_i}{c}_i\ne 0$ . An example is the Random Walk (RW) mixed model, with $c_1=1$ and $c_i=0$ for $i=\mathrm{2,}\dots \mathrm{,}n$ , cf. equation (A3) in Boer et al. (2020).

The LV and the RW mixed models have a sparse precision matrix and they can be solved in a computational efficient way (Boer et al., 2020). This approach can be further extended to show that there are mixed models, with sparse precision matrices, that are equivalent to P-splines. The method can be generalized to tensor product P-splines, as shown in Section 5.

3 B-splines

In this section we will give a brief overview of B-splines; for more details see for example De Boor (1978) and Lyche et al. (2018). Let the domain $[x_{\text{min}}\mathrm{,}{x}_{\text{max}}]$ , with $x_{\text{min}}<x_{\text{max}}$ and $x_{\text{min}}\mathrm{,}{x}_{\text{max}}\in ℝ$ , be divided in $n_{\text{seg}}$ intervals, each of length $h=(x_{\text{max}}-x_{\text{min}})/n_{\text{seg}}$ . The knots for the B-spline basis are given by

$${\kappa }_j=x_{\text{min}}-(p+1)h+jh\mathrm{,}j=\mathrm{1,}\dots \mathrm{,}{n}_{\text{seg}}+2p+\mathrm{1,}$$

(3.1)

where $p$ is the degree of the B-splines functions, denoted by $B_{j\mathrm{,}p}(x)$ . The B-spline functions have local support, that is, $B_{j\mathrm{,}p}(x)=0$ if $x\notin [{\kappa }_j\mathrm{,}{\kappa }_{j+p+1})$ . The number of B-splines is equal to $q=n_{\text{seg}}+p$ . Let $g(x)$ be a linear combination of the B-splines

$$g(x)={\displaystyle \sum _{j=1}^q}{u}_j{B}_{j\mathrm{,}p}(x),x\in [x_{\text{min}}\mathrm{,}{x}_{\text{max}}].$$

(3.2)

The $k$ th-order derivative of $g(x)$ can be expressed as follows (De Boor, 1978; Eilers and Marx, 1996)

$$\frac{\text{d}g(x)}{\text{d}{x}^k}=\frac{1}{h^k}{\displaystyle \sum _{j=k+1}^q}{\Delta }^k{u}_j{B}_{j\mathrm{,}p-k}(x),k\le p\mathrm{,}x\in [x_{\text{min}}\mathrm{,}{x}_{\text{max}}],$$

(3.3)

where ${\Delta }^k{u}_j$ is the $k$ -th order difference, recursively defined as ${\Delta }^k{u}_j={\Delta }^{k-1}(\Delta u_j)$ .

Let $f(x)$ be a $(k-1)$ -degree polynomial with coefficients ${\beta }_0\mathrm{,}\dots \mathrm{,}{\beta }_{k-1}$ . The representation of $f(x)$ in terms of B-splines is given by

$$f(x)={\displaystyle \sum _{i=0}^{k-1}}{\beta }_i{\displaystyle \sum _{j=1}^q}{\xi }_{i\mathrm{,}j\mathrm{,}p}{B}_{j\mathrm{,}p}(x),k-1\le p\mathrm{,}x\in [x_{\text{min}}\mathrm{,}{x}_{\text{max}}].$$

(3.4)

where ${\xi }_{i\mathrm{,}j\mathrm{,}p}$ are constants. Analytic expressions for ${\xi }_{i\mathrm{,}j\mathrm{,}p}$ can be found in Lyche et al. (2018). Taking the $k$ -th order derivative of equation (3.4) and using equation (3.3) gives

$${\Delta }^k{\xi }_{i\mathrm{,}j\mathrm{,}p}=0\text{\hspace{1em}for\hspace{1em}}i<k\mathrm{,}$$

(3.5)

a result that will be used to establish, in Sections 4 and 5, the equivalence of a class of mixed models with P-splines.

The idea of representing polynomials in terms of B-splines can be extended a bit further, and this idea will be used to derive our main result in Section 4. Suppose we are interested in the interpolating $(k-1)$ -degree polynomial (Press et al., 2002, p. 123) defined by $k$ pairs $(x_{*,i}\mathrm{,}{y}_{*,i})$ . For $k=1$ we will use $x_{*\mathrm{,1}}=(x_{\text{min}}+x_{\text{max}})/2$ . For $k>1$ the points $x_{*,i}$ are taken equidistant on the interval $[x_{\text{min}}\mathrm{,}{x}_{\text{max}}]$ :

$$x_{*,i}=x_{\text{min}}+(i-1)\frac{x_{\text{max}}-x_{\text{min}}}{k-1}\mathrm{,}i=\mathrm{1,}\dots \mathrm{,}k\mathrm{.}$$

(3.6)

Let $B_{*}$ be a $k\times q$ matrix, with $(i\mathrm{,}j)$ th entry $B_{j\mathrm{,}p}(x_{*,i})$ , and let $G$ be a $q\times k$ matrix with $(j\mathrm{,}i)$ th entry ${\xi }_{i\mathrm{,}j\mathrm{,}p}$ . The $k\times k$ matrix $M=B_{*}G$ is a Vandermonde matrix

$$M=B_{*}G=\left(\begin{array}{ccccc}1& x_{*\mathrm{,1}}& x_{*\mathrm{,1}}^2& \cdots & x_{*\mathrm{,1}}^{k-1}\\ 1& x_{*\mathrm{,2}}& x_{*\mathrm{,2}}^2& \cdots & x_{*\mathrm{,2}}^{k-1}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& x_{*,k}& x_{*,k}^2& \cdots & x_{*,k}^{k-1}\end{array}\right)$$

with determinant (Harville, 1998, section 13.6)

$$\left|M\right|=\left|B_{*}G\right|={\displaystyle \prod _{1\le i<j\le k}}(x_{*,j}-x_{*,i})\ne 0.$$

(3.7)

The determinant is unequal to zero, since the points $x_{*,i}$ are distinct.

4 A new formulation for P-splines

Using the B-splines properties described in the previous section, we can generalize the results from Section 2 to P-splines. First we introduce some notation. Let $y=(y_1\mathrm{,}{y}_2\mathrm{,}\dots \mathrm{,}{y}_n{)}^{\prime }$ be the response variable, depending on the variable $x=(x_1\mathrm{,}{x}_2\mathrm{,}\dots \mathrm{,}{x}_n{)}^{\prime }$ with $x_i\in [x_{\text{min}}\mathrm{,}{x}_{\text{max}}]$ for $i=\mathrm{1,}\dots \mathrm{,}n$ . Let $B$ be a $n\times q$ matrix, with $(i\mathrm{,}j)$ th entry $B_{j\mathrm{,}p}(x_i)$ . The matrix $X$ is defined by $X=BG=[1|x^1|\dots |x^{k-1}]$ .

In the next step we will define a model equivalent with P-splines, following the same steps as in Section 2. To make this connection more clear, we will use the same subscripts for the objective functions to be minimized, cf. equations (2.3), (2.6), and (2.7). The objective function $S_1$ to be minimized is given by

$$S_1(\beta \mathrm{,}u)=\left|\right|y-B\left(G\beta +u\right){\left|\right|}^2+\text{λ}u\text{'}D\text{'}Du+\text{λ}u\text{'}{B}_{*}\text{'}{B}_{*}u\mathrm{,}$$

(4.1)

where $\beta =({\beta }_0\mathrm{,}\dots \mathrm{,}{\beta }_{k-1}{)}^{\prime }$ , $u=(u_1\mathrm{,}\dots \mathrm{,}{u}_q{)}^{\prime }$ , and $D$ is a $(q-k)\times q$ matrix, with $k$ th-order differences. Next, we define the linear transformation

$$\left(\begin{array}{c}a\\ \eta \end{array}\right)=\left(\begin{array}{cc}{I}_q& G\\ B_{*}& 0\end{array}\right)\left(\begin{array}{c}u\\ \beta \end{array}\right)=\left(\begin{array}{c}u+G\beta \\ B_{*}u\end{array}\right)\mathrm{,}$$

(4.2)

with $a=(a_1\mathrm{,}\dots \mathrm{,}{a}_q{)}^{\prime }$ and $\eta =({\eta }_1\mathrm{,}\dots \mathrm{,}{\eta }_k{)}^{\prime }$ . This transformation is invertible since

$$\left|\left(\begin{array}{cc}{I}_q& G\\ B_{*}& 0\end{array}\right)\right|=\left|0-B_{*}{I}_{q}G\right|=(-{1)}^k\left|B_{*}G\right|=(-{1)}^k\left|M\right|\ne 0,$$

where we use the properties of block matrices in the first step (Harville, 1998, section 13.3), and equation (3.7) in the second and third steps. For the differences we have

$$Da=D(G\beta +u)=Du\mathrm{,}$$

(4.3)

using $DG=0$ , which is equation (3.5) in matrix notation.

By substitution of equation (4.2) into equation (4.1) and using equation (4.3) we obtain

$$S_2(a\mathrm{,}\eta )=\left|\right|y-Ba{\left|\right|}^2+\text{λ}a\text{'}{D}^{\prime }Da+\text{λ}\eta \text{'}\eta \mathrm{.}$$

(4.4)

Equation (4.4) can be decomposed as $S_2(a\mathrm{,}\eta )=S_3(a)+S_4(\eta )$ , where

$$S_3(a)=\left|\right|y-Ba{\left|\right|}^2+\text{λ}a\text{'}{D}^{\prime }Da$$

(4.5)

is the objective function to be minimized for P-splines (Eilers and Marx, 1996), and $S_4(\eta )=\text{λ}\eta \text{'}\eta$ achieves its minimum at $\widehat{\eta }=0$ . It follows that equations (4.1) and (4.5) are equivalent. The advantage of the new formulation is that we now have a computationally efficient formulation as a mixed model. This will be explained in more detail in Section 6.

5 Extension to tensor product P-splines

In this section we will extend the new formulation of P-splines to higher dimensions. We generalise the notation that was used for the one-dimensional case. Let $M$ be the number of dimensions, and $x_m$ the covariate for dimension $m$ . The matrix $B_m$ is the $n\times q_m$ matrix corresponding to dimension $m$ , $D_m$ is a $(q_m-k_m)\times q_m$ matrix, $G_m$ is a $q_m\times k_m$ matrix, and $B_{*,m}$ is a $k_m\times q_m$ matrix. Let $q={\displaystyle {\prod }_{m=1}^M}{q}_m$ be the total number of parameters, and $k={\displaystyle {\prod }_{m=1}^M}{k}_m$ . The Kronecker product is denoted by $\otimes$ and the row-wise Kronecker product (Eilers and Marx, 2003) is denoted by $\square$ . The matrices $B$ , $B_{*}$ and $G$ are defined by

$$B=B_1\square \cdots \square B_M\mathrm{,}{B}_{*}=\underset{m=1}{\overset{M}{{\displaystyle \otimes }}}{B}_{*,m}\mathrm{,}G=\underset{m=1}{\overset{M}{{\displaystyle \otimes }}}{G}_m\mathrm{.}$$

The fixed part of the model is given by $X=BG$ . For the difference penalty corresponding to dimension $m$ we have:

$${\mathcal{D}}_m=\underset{j=1}{\overset{M}{{\displaystyle \otimes }}}{\mathcal{H}}_j\mathrm{,}{\mathcal{H}}_j=\left\{\begin{array}{cc}{D}_j& j=m\\ I_{q_j}& j\ne m\end{array}\right.\mathrm{,}$$

For the product between $G$ and ${\mathcal{D}}_m$ we obtain:

$${\mathcal{D}}_{m}G=\left(\underset{j=1}{\overset{M}{{\displaystyle \otimes }}}{\mathcal{H}}_j\right)\left(\underset{j=1}{\overset{M}{{\displaystyle \otimes }}}{G}_j\right)=\underset{j=1}{\overset{M}{{\displaystyle \otimes }}}{\mathcal{H}}_j{G}_j=0\mathrm{,}m=\mathrm{1,}\dots \mathrm{,}M\mathrm{,}$$

using that ${\mathcal{H}}_m{G}_m=D_m{G}_m=0$ .

Define the tensor product objective function $S_1$ to be minimized:

$$S_1(\beta \mathrm{,}u)=\left|\right|y-BG\beta -Bu{\left|\right|}^2+{\displaystyle \sum _{m=1}^M}{\text{λ}}_{m}u\text{'}{\mathcal{D}}_m^{\text{'}}{\mathcal{D}}_{m}u+{\displaystyle \sum _{m=1}^M}{\text{λ}}_{m}u\text{'}{B}_{*}\text{'}{B}_{*}u$$

(5.1)

Using the linear transformation defined in equation (4.2) and

$${\mathcal{D}}_{m}a={\mathcal{D}}_m\left(G\beta +u\right)={\mathcal{D}}_{m}u\mathrm{,}m=\mathrm{1,}\dots \mathrm{,}M$$

we obtain

$$S_2(a\mathrm{,}\eta )=\left|\right|y-Ba{\left|\right|}^2+{\displaystyle \sum _{m=1}^M}{\text{λ}}_{m}a\text{'}{\mathcal{D}}_m^{\text{'}}{\mathcal{D}}_{m}a+\left({\displaystyle \sum _{m=1}^M}{\text{λ}}_m\right)\eta \text{'}\eta \mathrm{.}$$

(5.2)

Since $S_2$ can be decomposed as $S_2(a\mathrm{,}\eta )=S_3(a)+S_4(\eta )$ and $S_4$ attains its minimum at $\widehat{\eta }=0$ , we obtain the tensor product P-splines model introduced by Eilers and Marx (2003)

$$S_3(a)=\left|\right|y-Ba{\left|\right|}^2+{\displaystyle \sum _{m=1}^M}{\text{λ}}_{m}a\text{'}{\mathcal{D}}_m^{\text{'}}{\mathcal{D}}_{m}a\mathrm{.}$$

(5.3)

This shows that equations (5.1) and (5.3) are equivalent. In Section 6 we will present an efficient mixed model formulation of equation (5.1). Rodríguez-Ávarez et al. (2015) use another transformation of tensor P-splines to a mixed model, which is based on the spectral decomposition of the difference penalty matrices $D_m^{\text{'}}{D}_m$ . For a detailed explanation of this mixed model representation for data on a multidimensional grid, see Section 6 in Currie et al., (2006). As will be shown in Section 6, the new approach presented in this article is considerably faster then the algorithm in Rodríguez-Ávarez et al. (2015).

6 The LMMsolver R-package

The LMMsolver R-package on CRAN is a general Linear Mixed Model (LMM) solver. The aim of the LMMsolver package is to provide an efficient and flexible system to estimate variance parameters using restricted maximum likelihood or REML (Patterson and Thompson, 1971). The package was developed specifically for models where the mixed model equations are sparse, including the mixed model P-splines introduced in the previous sections. An example of the use of LMMsolver for a standard mixed model and allowing for heterogeneous residual errors is in the R-package statgenMPP (Li et al., 2022). In this article we will only give the mixed model for the tensor P-splines, although extra fixed or random terms can be added to the mixed model in LMMsolver.

6.1 Solving the linear mixed equations

The penalized regression model defined by equation (5.1) can be formulated as a mixed model:

$$y=X\beta +Zu+e\mathrm{,}u\sim N(0\mathrm{,}\Sigma ),e\sim N(0\mathrm{,}R),$$

(6.1)

with $X=BG$ , $Z=B$ and

$$R^{-1}={\theta }_0{I}_n\mathrm{,}{\Sigma }^{-1}={\displaystyle \sum _{m=1}^M}{\theta }_m\left({\mathcal{D}}_m^{\text{'}}{\mathcal{D}}_m+B_{*}\text{'}{B}_{*}\right)\mathrm{,}$$

(6.2)

where ${\theta }_0=1/{\sigma }_e^2$ and ${\theta }_m={\lambda }_m/{\sigma }_e^2(m=\mathrm{1,}\dots \mathrm{,}M)$ are precision parameters to be estimated. For $M>1$ equation (6.1) is not a standard mixed model, it has overlapping penalties (Currie et al., 2006; Rodríguez-Ávarez et al., 2015; Rodríguez-Ávarez et al., 2019). The mixed model equations corresponding to equation (6.1) are defined by

$$\left(\begin{array}{cc}X\text{'}{R}^{-1}X& X\text{'}{R}^{-1}Z\\ Z\text{'}{R}^{-1}X& Z\text{'}{R}^{-1}Z+{\Sigma }^{-1}\end{array}\right)\left(\begin{array}{c}\widehat{\beta }\\ \widehat{u}\end{array}\right)=\left(\begin{array}{c}X\text{'}{R}^{-1}y\\ Z\text{'}{R}^{-1}y\end{array}\right)\mathrm{.}$$

(6.3)

The matrix on the left hand side of equation (6.3) is called the mixed model coefficient matrix and will be denoted by $C$ . This matrix is sparse and therefore equation (6.3) can be solved in a computational efficient way using the Cholesky decomposition of $C$ (Furrer and Sain, 2010).

The followings equations are used to update the precision parameters ${\theta }_i$ (see Appendix A for the derivation):

$$\begin{array}{ll}{\theta }_0\hfill & =\frac{n-{\rho }_0}{\widehat{e}\text{'}\widehat{e}}\hfill \\ {\theta }_i\hfill & =\frac{{\tau }_i-{\rho }_i}{\widehat{u}\text{'}({\mathcal{D}}_i^{\text{'}}{\mathcal{D}}_i+B_{*}\text{'}{B}_{*})\widehat{u}}i=\mathrm{1,}\dots \mathrm{,}M\mathrm{,}\hfill \end{array}\mathrm{,}$$

(6.4)

where $\widehat{e}=y-X\widehat{\beta }-Z\widehat{u}$ and

$$\begin{array}{ll}{\rho }_i\hfill & ={\theta }_i\frac{\partial \log \left|C\right|}{\partial {\theta }_i}={\theta }_i\text{trace}\left(C^{-1}\frac{\partial C}{\partial {\theta }_i}\right)\mathrm{,}i=0,\dots \mathrm{,}M\hfill \\ {\tau }_i\hfill & ={\theta }_i\frac{\partial \log \left|{\Sigma }^{-1}\right|}{\partial {\theta }_i}={\theta }_i\text{trace}\left(\Sigma \frac{\partial {\Sigma }^{-1}}{\partial {\theta }_i}\right)\mathrm{,}i=\mathrm{1,}\dots \mathrm{,}M\mathrm{.}\hfill \end{array}$$

(6.5)

The iterative algorithm defined by equations (6.3) and (6.5) is an extension of the modified Henderson algorithm described in Harville (1977).

An important element needed is to calculate ${\rho }_i$ and ${\tau }_i$ defined by equation (6.5) in an efficient way, avoiding the calculation of the matrices $C^{-1}$ and $\Sigma$ , which are not sparse. One way to do this is to calculate the so-called sparse inverse (Takahashi, 1973; Misztal and Perez-Enciso, 1993; Johnson and Thompson, 1995; Lindgren et al., 2022). In LMMsolver another, but closely related, method was implemented, using Automated Differentiation of the Cholesky algorithm (Smith, 1995). Backward differentiation was implemented, which calculates the partial derivatives of the likelihood efficiently (Meyer and Smith, 1996; Smith, 1995). The chol() function in the spam package was used to permute the row and columns of $C$ and ${\Sigma }^{-1}$ in an optimal way, see Furrer and Sain (2010) for details. The automated differentation was implemented using supernodal Cholesky factorization (Ng and Peyton, 1993, Furrer and Sain, 2010). The implementation was written in C++ using the Rcpp package (Eddelbuettel and Balamuta, 2018).

6.2 Example 1: United States precipitation data set

In this section we will compare the performance of LMMsolver with the R-packages mgcv (Wood, 2017) and SOP (Rodríguez-Ávarez et al., 2019). The method in the SOP package uses the same tensor P-splines model as in LMMsolver. For a detailed comparison between tensor P-splines in mgcv and SOP see Rodríguez-Ávarez et al. (2015). We will use two-dimensional P-splines for the USA precipitation data set (Rodríguez-Ávarez et al., 2015). There are $n=5,906$ observations, with longitude-latitude positions of monitoring stations, and total precipitation in millimeters and a standardization of this raw observation for April 1948. This data set can be found in the R-package spam (Furrer and Sain, 2010), under the name USprecip.

All computations were performed in R4.2.3 (R Core Team 2023) and a 2.90GHz Intel Core i5-9400 CPU with 24GB of RAM and Windows10 operating system. Version 1.8.42 of mgcv, version 1.0 of SOP, and version 1.0.5 of LMMsolver were used. For mgcv we used the bam() function, with method=‘fREML’. Cubic B-splines with second-order differences were used for both latitude and longitude. The differences between the estimated effective dimensions (Rodríguez-Ávarez et al., 2015) for the SOP and LMMsolver packages were less than $5\cdot {10}^{-3}$ . We compared the computation times for different number of segments. The same number was used in both dimensions.

Figure 1A shows the sparse structure for the mixed model coefficient matrix corresponding to equation (5.1) for $40$ segments.

Figure 1

Panel A: The mixed model coefficient matrix $C$ for the new approach for the US precipitation example with $40$ segments for both latitude and longitude. The non-zeros are indicated in red, showing that $C$ is sparse. Panel B: Computation time on a logarithmic scale for the three methods as function of the number of segments in both dimensions: LMMsolver is much faster, especially if the number of segments is high.

In Figure 1B the computation time is compared for different number of segments, showing that LMMsolver is much faster than mgcv and SOP. For example, if the number of segments is 40, the computation time for mgcv is 600 seconds (= 10 minutes), for SOP 140 seconds, and for LMMsolver only one second. The difference in computation time between mgcv and SOP is a factor $4$ , consistent with the results in Rodríguez-Ávarez et al., (2015).

There are two main reasons why LMMsolver is fast for this two-dimensional mixed model P-splines example. First, it retains the sparse structure of the original two-dimensional P-splines as in Eilers and Marx (2003). The second reason is an efficient implementation of the Automated Differentiation of the Cholesky Algorithm in LMMsolver, which avoids the calculations of the inverse of matrices for the derivatives of the REML log-likelihood (Smith, 1995; Meyer and Smith, 1996).

6.3 Example 2: Three-dimensional simulated data

In this section we will use a simulated three-dimensional example to illustrate that LMMsolver can analyze large datasets in an efficient way. One hundred thousand ( $n=\mathrm{100,000}$ ) values of covariates $z_1$ , $z_2$ , $z_3$ were simulated independently from a uniform distribution on the interval $[\mathrm{0,1}]$ , and the response was generated from the equation (Wood, 2006; Rodríguez-Ávarez et al., 2015)

$$y=1.5\text{\hspace{0.33em}exp}\left(-\frac{{(z_1-0.2)}^2}{5}-\frac{{(z_2-0.5)}^2}{3}-\frac{{(z_3-0.9)}^2}{4}\right)$$

$$+0.5\text{\hspace{0.33em}exp}\left(-\frac{{(z_1-0.3)}^2}{4}-\frac{{(z_2-0.7)}^2}{2}-\frac{{(z_3-0.4)}^2}{6}\right)$$

$$+\text{exp}\left(-\frac{{(z_1-0.1)}^2}{5}-\frac{{(z_2-0.3)}^2}{5}-\frac{{(z_3-0.7)}^2}{4}\right)+ϵ\mathrm{,}$$

with $ϵ\sim N(\mathrm{0,1})$ . Figure 2A shows the structure of mixed model coefficient matrix $C$ , which is less sparse than the two-dimensional case shown in Figure 1A. Therefore the differences in computation time between LMMsolver and the packages mgcv and SOP is less than in the two-dimensional case, but still substantial: if ten segments are used in each direction, the computation time for mgcv is 38 minutes, for SOP 10 minutes, and for LMMsolver half a minute.

7 Discussion

In this article we have shown that there is a mixed model formulation for P-splines that adheres closely to the original proposal by Eilers and Marx (1996). The main difference is that in the new approach there is an explicit term for the fixed part, with an additional term for the penalty. An important feature is that the precision matrices and the mixed model equations are sparse. Therefore the sparse matrix algebra developed for mixed models can be used.

An additional advantage of the new approach is that the mixed model P-splines formulation can be further extended to Generalized Additive Models (GAM; Hastie et al., 2009). In standard P-splines a small ridge penalty is added to make the system identifiable (Eilers and Marx, 2021). In the approach described here, this problem does not arise: the column corresponding to ${\beta }_0$ can be left out of the model, and one intercept term is modelled. This idea is similar to leaving out one column for a factor in the fixed part of the model. GAMs are implemented in LMMsolver.

The mixed model P-splines formulation extends the P-splines framework developed by Paul Eilers and Brian Marx (1996, 2021). An interesting aspect of the mixed model P-splines formulation is that it takes advantage of the properties of the two key components of P-splines: B-splines and discrete difference penalties.

Figure 2

Panel A: The mixed model coefficient matrix $C$ for the new approach for the three-dimensional simulated data with $5$ segments for each dimension. The non-zeros are indicated in red. Panel B: Computation time on a logarithmic scale for LMMsolver, mgcv, and SOP as function of the number of segments in each dimension.

Appendix A

In this appendix we will give a derivation of equation (6.4). The restricted log-likelihood of equation (6.1) is given by (Meyer and Smith, 1996):

$$\log \mathcal{L}=-\frac{1}{2}\left(\text{const}+\log \left|R\right|+\log \left|\Sigma \right|+\log \left|C\right|+y\text{'}Py\right)$$

(A.1)

where $P=V^{-1}-V^{-1}X{(X\text{'}{V}^{-1}X)}^{-1}X\text{'}{V}^{-1}$ , and $V=Z\Sigma Z\text{'}+R$ . The matrix $P$ has the following properties (Johnson and Thompson, 1995): $X\text{'}P=0$ , $Py=R^{-1}\widehat{e}$ , and $Z\text{'}Py={\Sigma }^{-1}\widehat{u}$ . Using these properties of $P$ it follows that

$$y\text{'}Py=(X\widehat{b}+Z\widehat{u}+\widehat{e}{)}^{\prime }Py=\widehat{u}\text{'}{\Sigma }^{-1}\widehat{u}+\widehat{e}\text{'}{R}^{-1}\widehat{e},*6pt$$

(A.2)

and therefore equation (A.1) can be rewritten as:

$$\log \mathcal{L}=\frac{1}{2}\left(\text{const}+\log \left|R^{-1}\right|+\log \left|{\Sigma }^{-1}\right|-\log \left|C\right|-\widehat{e}\text{'}{R}^{-1}\widehat{e}-\widehat{u}\text{'}{\Sigma }^{-1}\widehat{u}\right)*6pt$$

(A.3)

with matrices $R^{-1}$ , ${\Sigma }^{-1}$ , and $C$ being sparse for the mixed model P-splines model, and therefore $\log \mathcal{L}$ and the partial derivatives with respect to precision parameters ${\theta }_m$ can be calculated in a computationally efficient way (Meyer and Smith, 1996).

The partial derivatives of equation (A.3) are given by

$$\begin{array}{ll}\frac{\partial \log \mathcal{L}}{\partial {\theta }_0}\hfill & =\frac{1}{2{\theta }_0}\left(n-{\rho }_0\right)-\widehat{e}\text{'}\widehat{e}\hfill \\ \frac{\partial \log \mathcal{L}}{\partial {\theta }_i}\hfill & =\frac{1}{2{\theta }_i}\left({\tau }_i-{\rho }_i\right)-\widehat{u}\text{'}({\mathcal{D}}_i^{\text{'}}{\mathcal{D}}_i+B_{*}\text{'}{B}_{*})\widehat{u}i=\mathrm{1,}\dots \mathrm{,}M\mathrm{,}\hfill \end{array}*6pt$$

(A.4)

where the parameters ${\rho }_i$ and ${\tau }_i$ are defined by equation (6.5).

Setting the partial derivatives in equation (A.4) equal to zero gives the expressions for ${\theta }_i$ ( $i=\mathrm{0,1,}\dots \mathrm{,}M$ ) in equation (6.4).

Supplementary materials

Supplementary materials for this article are available online.

Supplemental Material for Tensor product P-splines using a sparse mixed model formulation by Martin P. Boer, in Statistical Modelling