P-splines with an l1 penalty for repeated measures
Brian D. Segal Michael R. Elliott Thomas Braun Hui Jiang
stat.ME 62G08 (Primary), 62P10 (Secondary)



[
    Jorge Drumond Silva
    December 30, 2023
=======================




P-splines are penalized B-splines, in which finite order differences in coefficients are typically penalized with an ℓ_2 norm.  P-splines can be used for semiparametric regression and can include random effects to account for within-subject variability. In addition to ℓ_2 penalties, ℓ_1-type penalties have been used in nonparametric and semiparametric regression to achieve greater flexibility, such as in locally adaptive regression splines, ℓ_1 trend filtering, and the fused lasso additive model. However, there has been less focus on using ℓ_1 penalties in P-splines, particularly for estimating conditional means.

In this paper, we demonstrate the potential benefits of using an ℓ_1 penalty in P-splines with an emphasis on fitting non-smooth functions. We propose an estimation procedure using the alternating direction method of multipliers and cross validation, and provide degrees of freedom and approximate confidence bands based on a ridge approximation to the ℓ_1 penalized fit. We also demonstrate potential uses through simulations and an application to electrodermal activity data collected as part of a stress study.


















§ INTRODUCTION


Many nonparametric regression methods, including smoothing splines and regression splines, obtain point estimates by minimizing a penalized negative log-likelihood function of the form l_pen = -l(β) + λ P(β), where l is a log-likelihood, P is a penalty term, λ > 0 is a smoothing parameter, and β are the coefficients to be estimated. Typically, quadratic (ℓ_2 norm) penalties are used, which lead to straightforward computation and inference. In particular, ℓ_2 penalties typically lead to ridge estimators, which have both closed form solutions and are linear smoothers. The ℓ_2 penalty also has connections to mixed models, which allows the smoothing parameters to be estimated as variance components <cit.>.

However, nonparametric regression methods that use an ℓ_1-type penalty, such as ℓ_1 trend filtering 

<cit.> and locally adaptive regression splines <cit.>, are better able to adapt to local differences in smoothness and achieve the minimax rate of convergence for weakly differentiable functions of bounded variation <cit.>, whereas ℓ_2 penalized methods do not <cit.>. The trade-off is that ℓ_1 penalties generally lead to more difficult computation and inference because the objective function is convex but non-differentiable, and the fit is no longer a linear smoother.

In this article, we propose P-splines with an ℓ_1 penalty as a framework for generalizing ℓ_1 trend filtering within the context of repeated measures data and semiparametric (additive) models <cit.>. In Section <ref>, we discuss the connection between P-splines and ℓ_1 trend filtering which motivates the methodological development. In Section <ref>, we present our proposed model, and in Section <ref>, we discuss related work. In Section <ref>, we propose an estimation procedure using the alternating direction method of multipliers (ADMM) <cit.> and cross validation (CV). In Section <ref>, we derive the degrees of freedom and propose computationally stable and fast approximations, and in Section <ref>, we develop approximate confidence bands based on a ridge approximation to the ℓ_1 fit. In Section <ref>, we study our method through simulations and evaluate its performance in fitting non-smooth functions. In section <ref>, we demonstrate our method in an application to electrodermal activity data collected as part of a stress study. We close with a discussion in Section <ref>.



§ P-SPLINES AND ℓ_1 TREND FILTERING 


In this section, we give brief background on P-splines and ℓ_1 trend filtering, and show the relation between them when the data are independent and identically distributed (i.i.d.) normal.

P-splines <cit.> are penalized B-splines <cit.>. B-splines are flexible bases that are notable in part because they have compact support, which leads to banded design matrices and faster computation. This compact support can be seen in Figure <ref>, which shows eight evenly spaced first degree and third degree B-spline bases on [0,1]. We can define an order M (degree M-1) B-spline basis with j=1,…,p basis functions recursively as <cit.>

    ϕ_j^m(x)    = x - t_j/t_j+m-1 - t_jϕ_j^m-1(x) + t_j+m - x/t_j+m - t_j+1ϕ_j+1^m-1(x),        j = 1,…, 2M + c -m, 
              1 < m ≤ M 
    ϕ_j^1(x)    = 
    
    1     t_j ≤ x < t_j+1
    
    0    otherwise,        j = 1,…,2M +c -1

where t_j are the knots, division by zero is taken to be zero, and c is the number of internal knots. For order M B-splines defined on the interval [x_min, x_max], in order to obtain j=1,…,p basis functions, we set 2M boundary knots (M knots on each side) and c=p-M interior knots. In general, one can set t_1 ≤ t_2 ≤⋯≤ t_M = x_min < t_M+1 < ⋯ < t_M+c < x_max = t_M+c+1≤ t_M+c+2≤⋯≤ t_2M + c. In order to ensure continuity at the boundaries, we set t_1 < t_2 < ⋯ < t_M-1 < t_M = x_min and x_max = t_M+c+1 < t_M+c+2 < ⋯ < t_2M + c. We also use equally spaced interior knots, which is important for the P-spline penalty, and drop the superscript on ϕ designating order when the order does not matter or is stated in the text.



B-spline bases can be used to fit nonparametric models of the form y(x) = f(x) + ϵ(x), where y(x) is the outcome y at point x, f(x) is the mean response function at x, and ϵ(x) is the error at x. To that end, let y=(y_1,…,y_n)^T be an n × 1 vector of outcomes and x = (x_1 …, x_n)^T be a corresponding n × 1 vector of covariates. Also, let ϕ_1,…,ϕ_p be B-spline basis functions and let F be an n × p design matrix such that F_ij = ϕ_j(x_i), i.e., the j^th column of F is the j^th basis function evaluated at x_1,…,x_n. Equivalently, the i^th row of F is the i^th data point evaluated by ϕ_1,…,ϕ_p. For i.i.d. normal y, a simple linear P-spline model with the standard ℓ_2 penalty can be written as

    β̂_0, β̂ = _β_0 ∈ℝ, β∈ℝ^p1/2y - β_0 1 - Fβ_2^2 + λ/2D^(k+1)β_2^2,

where β_0 is the intercept, β is a p × 1 vector of parameter estimates, 1 is an n × 1 vector with each element equal to 1, λ>0 is a smoothing parameter, and D^(k + 1)∈ℝ^(p-k-1) × p is the k+1 order finite difference matrix. For example, for k=1

    D^(2)  =
    [  1 -2  1   ;     1 -2  1;        ⋱  ⋱  ⋱;           1 -2  1;              1 -2  1 ]∈ℝ^(p-2) × p

In general, as described by <cit.>, D^(k+1) = D^(1) D^(k) where D^(1) is the (p-k-1) × (p-k) upper left matrix of:

    D^(1) = 
    [ -1  1   ;    -1  1   ;        ⋱  ⋱;          -1  1 ]∈ℝ^(p-1) × p.

Our proposed model builds on one in which the ℓ_2 penalty in (<ref>) is replaced with an ℓ_1 penalty:

    β̂_0, β̂ = _β_0 ∈ℝ, β∈ℝ^p1/2y - β_0 1 - Fβ_2^2 + λD^(k+1)β_1.


Letting f(x) = ∑_j=1^p β_j ϕ_j^M(x), for order M=4 B-splines, <cit.> show that 

    ∫_x_min^x_max( d^2/dx^2 f(x) )^2 dx = c_1 D^(2)β_2^2 + c_2 ∑_j=4^p ∇^2 β_j ∇^2 β_j-1

where ∇^2 is the second-order backwards difference and c_1 and c_2 are constants. As shown in Appendix <ref>, a similar result holds for P-splines with an ℓ_1 penalty. In particular, for 0 ≤ k < M-1,

    ∫_x_min^x_max| d^k+1/dx^k+1 f(x) | dx ≤ C_M,k+1D^(k+1)β_1

where C_M, k+1 is a constant given in Appendix <ref> that depends on the order M of the B-splines and order k+1 of the finite difference. In other words, controlling the ℓ_1 norm of the (k+1)^th order finite differences in coefficients also controls the total variation of the k^th derivative of the function.

ℓ_1 trend filtering is similar to (<ref>). In the case where x_1 < x_2 < ⋯ <x_n are unique and equally spaced, ℓ_1 trend filtering solves the following problem (the intercept is handled implicitly):

    β̂ = _β∈ℝ^n1/2y - β_2^2 + λD^(k+1)β_1.

Problem (<ref>) differs from (<ref>) in that (<ref>) has one parameter per data point, and the design matrix is the identity matrix. D^(k+1) is also resized appropriately by replacing p with n in the dimensions of (<ref>) and (<ref>). However, under certain conditions noted in Observation <ref>, (<ref>) and (<ref>) are identical.

[Continuous representation]
For second order (first degree) B-splines with n basis functions, equally spaced data x_1 < x_2 < ⋯ < x_n with knots at t_1 < x_1, t_2 = x_1, t_3=x_2,…,t_n = x_n-1, t_n+1 = x_n, t_n+2 > x_n, and centered outcomes such that y(0) = 0, P-splines with an ℓ_1 penalty are a continuous analogue to ℓ_1 trend filtering.
  



Under these conditions, for i=1,…,n

    ϕ^2_j(x_i) = 
    
    1     i = j 
    
    0    otherwise.

To see this, note that

    ϕ_j^2(x_i)    = x_i - t_j/t_j+1 - t_jϕ_j^1(x_i) + t_j+2 - x_i/t_j+2 - t_j+1ϕ_j+1^1(x_i) 
       = t_i+1 - t_j/t_j+1 - t_jϕ_j^1(t_i+1) + t_j+2 - t_i+1/t_j+2 - t_j+1ϕ_j+1^1(t_i+1).

Now, 

    ϕ_j^1(t_i+1) = 
    
    1     t_j ≤ t_i+1 < t_j+1
    
    0    otherwise  and  ϕ_j+1^1(t_i+1) = 
    
    1     t_j+1≤ t_i+1 < t_j+2
    
    0    otherwise.

We have ϕ_j^1(t_i+1) = 1 for i=j-1 and 0 otherwise, but for i=j-1, we have t_i+1 - t_j = t_j - t_j = 0. We also have ϕ_j+1^1(t_i+1)=1 for i=j and 0 otherwise, and for i=j, we have t_j+2 - t_i+1 = t_j+2 - t_j+1 > 0. It follows that for i=1…,n, (<ref>) evaluates to 1 if i=j and 0 otherwise.

Let F be the design matrix in (<ref>), where F_ij = ϕ^2_j(x_i). Then from the previous result, we have F=I_n, where I_n is the n × n identity matrix. This, together with the assumption that β_0 = y(0) = 0, implies that the objective functions (<ref>) and (<ref>) are identical, which proves Observation <ref>.


We note that <cit.> shows that ℓ_1 trend filtering has a continuous representation when expressed in the standard lasso form, and Observation <ref> gives a continuous representation of ℓ_1 trend filtering when expressed in generalized lasso form.

ℓ_1 trend filtering can be applied to irregularly spaced data, such as with the algorithm developed by <cit.>. It might also be possible to extend ℓ_1 trend filtering to repeated measures data to account for within-subject correlations. However, due to Observation <ref>, we think it is beneficial to view ℓ_1 trend filtering as a special case of P-splines with an ℓ_1 penalty. We think this approach has the potential to be a general framework, because higher order B-splines could be used in combination with different order difference matrices just as can be done with P-splines that use the standard ℓ_2 penalty. Furthermore, expressing ℓ_1 trend filtering as P-splines with an ℓ_1 penalty may facilitate the development of confidence bands (see Section <ref>), which could help to fill a gap in the ℓ_1 penalized regression literature.

In addition, there are connections between P-splines with an ℓ_1 penalty and locally adaptive regression splines. In particular, as <cit.> shows, the continuous analogue of ℓ_1 trend filtering is identical to locally adaptive regression splines <cit.> for k=0,1, and asymptotically equivalent for k ≥ 2.



§ PROPOSED MODEL: ADDITIVE MIXED MODEL USING P-SPLINES WITH AN ℓ_1 PENALTY 


To introduce our model, let y_i = (y_i1, …, y_i n_i)^T be an n_i × 1 vector of responses for subject i=1,…, N, and let y = (y_1^T, …, y_N^T)^T be the stacked n × 1 vector of responses for all N subjects, where n= ∑_i=1^N n_i. Let x_i = (x_i1, …, x_i n_i)^T be a corresponding n_i × 1 vector of covariates for subject i, and x = (x_1^T,…,x_N^T)^T be the n × 1 stacked vector of all covariate values. In many contexts, x is time. To account for the within-subject correlations of y_i, we can incorporate random effects into the P-spline model. To that end, let Z_i be an n_i × q_i design matrix for the random effects for subject i (possibly including a B-spline basis), and let b_i = (b_i1, …, b_iq_i)^T be the corresponding q_i × 1 vector of random effect coefficients for subject i. Also, let

    Z = [ Z_1        ;       ⋱;         Z_N ]
 be the n × q block diagonal random effects design matrix for all subjects, where q=∑_i=1^N q_i, and let b = (b_1^T,…,b_N^T)^T be the q × 1 stacked vector of random effects for all subjects. We propose an additive mixed model with j=1,…,J smooths (tildes denote quantities that will be subject to additional constraints, as described below):

    β_0 ∈ℝ, b∈ℝ^q, β̃_j ∈ℝ^p_j, j = 1,…,Jminimize     1/2y - β_0 1 - ∑_j=1^J F̃_j β̃_j  - Z b_2^2 + ∑_j=1^J λ_j D̃_j^(k_j+1)β̃_j _1 
       + τ1/2b^T S b

where F̃_j is a n × p_j design matrix of B-spline bases for smooth j, D̃_j^(k_j+1) is the k_j + 1 finite difference matrix, and σ^2_b S is the covariance matrix of the random effects b. For example, if b are random intercepts, then S=I_N and Z would be an n × N matrix such that Z_il = 1 if observation i belonged to subject l and zero otherwise. Alternatively, to obtain random curves using smoothing splines and a B-spline basis, we could set

    S = [ S_1        ;       ⋱;         S_N ]

where S_j, il = ∫ϕ”_ji(t) ϕ”_jl(t) dt, and ϕ”_j1,…,ϕ”_jp_j are the second derivatives of the B-spline basis functions for the j^th smooth. We would then set Z to be the corresponding B-splines evaluated at the input points.

We note that (<ref>) includes varying-coefficient models <cit.>. For example, as pointed out by <cit.>, if F̃_1 are B-splines evaluated at x, we could have F̃_2 = diag(x') F̃_1, where x'x is another covariate vector and diag(x') is a diagonal matrix with x'_i at the i^th leading diagonal position.

As written, (<ref>) is not generally identifiable. To see this, suppose ŷ(x) = β̂_0 + f̂_1(x) + f̂_2(x), where neither f_1 nor f_2 are varying-coefficient terms. Then letting f̂'_1(x) = f̂_1(x) + δ and f̂'_2(x) = f̂_2(x) - δ for δ∈R, we also have ŷ(x) = β̂_0 + f̂'_1(x) + f̂'_2(x). To make (<ref>) identifiable, we follow <cit.> and introduce a centering constraint on each non-varying coefficient smooth, i.e. ∫f̂_j(x) dx = 0 for all smooths j=1,…, J such that F̃_j diag(x')F̃_l for some x' and l  j. To this end, let ℰ = {j ∈{1,…, J} : F̃_j diag(x')F̃_l  for some x', l  j} be the indices of the non-varying coefficient smooths, and let ℰ̅ = {j ∈{1,…, J} : j ∉ℰ} be its complement. We constrain 1^T F̃_j β̃_j = 0 for j ∈ℰ. We accomplish this by defining new p_j × (p_j - 1) orthonormal matrices Q_j, j=1, …, J, such that 1^T F̃_j Q_j = 0. If desired, one can also define a q × (q-1) matrix Q_J+1 such that 1^T Z Q_J+1 = 0. However, this last centering constraint is not necessary, because the penalty on the random effect terms pulls the coefficients themselves towards zero, as opposed to the finite order differences in coefficients.

As <cit.> shows, Q can be obtained by taking the QR decomposition of F̃^T_j 1 and retaining the last p_j -1 columns of the left orthonormal matrix.[The matrices 1^T F̃_j, j=1,…, J are of rank 1, so the remaining p_j-2 columns are arbitrary orthonormal vectors. In  <cit.>, when taking the QR decomposition of F̃^T 1, an appropriate matrix Q can be obtained as .] We can then re-parameterize the p_j constrained parameters β̃_j in terms of the p_j-1 unconstrained parameters β_j, such that β̃_j = Q_j β_j. For j ∈ℰ, let F_j = F̃_j Q_j and D_j = D̃_j^(k_j+1)Q_j. For j ∈ℰ̅, let F_j = F̃_j and D_j = D̃^(k_j+1)_j. If centering the random effects, then we redefine S := Q_J+1^T S Q_J+1 and Z := Z Q_J+1. Then we can re-write (<ref>) in the identifiable form

    β_0 ∈ℝ, b∈ℝ^q, β_j ∈ℝ^p'_j, j = 1,…,Jminimize     1/2y - β_0 1 - ∑_j=1^J F_j β_j  - Z b_2^2 + ∑_j=1^J λ_j D_j β_j _1 
       + τ1/2b^T S b

where p_j' = p_j - 1 for j ∈ℰ and p_j' = p_j for j ∈ℰ̅.

We note that the penalty matrix S given above for random subject-specific splines defines non-zero correlation between nearby within-subject random effect coefficients. This is in contrast to the approach of <cit.> for estimating subject-specific random curves, which focuses on the case in which nearby within-subject coefficients are not correlated. To see this, let d̂_i(x) = ∑_j=1^q_ib̂_ijϕ_ij(x) be the estimated difference between the i^th subject-specific curve and the marginal mean at point x. The smoothing spline approach above constrains ∫ (d̂”)^2(x)dx = b_i^T S_i b_i <C for some constant C>0, whereas the approach of <cit.> constrains b_i^T I_q_ib_i = ∑_j=1^q_ib̂_j^2 <C. Whereas the non-diagonal penalty matrix S implies correlations between nearby coefficients, the identity matrix in the approach of <cit.> implies zero correlation.

Similar to the equivalence between Bayesian models and ℓ_2 penalized smoothing splines <cit.>, there is an equivalence between Bayesian models and ℓ_1 penalized splines. In particular, (<ref>) is equivalent to the following distributional assumptions, which we can use to obtain Bayesian estimates:

    y|b   = β_0 1 + ∑_j=1^J F_j β_j + Z b + ϵ
    ϵ   ∼ N ( 0, σ^2_ϵ I_n ) 
    b   ∼ N(0, σ^2_b S^-1)  for σ^2_b = σ^2_ϵ / τ
    ϵ   ⊥b
    ( D_j β_j )_l    ∼Laplace(0, a_j)  for  a_j = σ^2_ϵ / (2 λ_j), l=1,…, p_j - k_j - 1, j=1,…,J

The last distributional assumption is an element-wise Laplace prior on the k_j+1 order differences in coefficients.

In some cases, the random effects penalty matrix S may be positive semidefinite but not invertible. For example, the smoothing spline random curves outlined above lead to a penalty matrix S that is not strictly positive definite, but that is still positive semidefinite. This does not cause problems for the ADMM algorithm, but some changes are required for other algorithms as well as for Bayesian estimation. Following <cit.>, let S=U Λ U^T be the eigendecomposition of a positive semidefinite matrix S, where U U^T = I_q and Λ is a diagonal matrix with eigenvalues in descending order in the diagonal positions. Let b̆ = U^T b and Z̆ = Z U, so that b^T S b = b̆^T Λb̆ and Z̆b̆ = Z b. Let q_r be the number of strictly positive eigenvalues of S, where 0 < q_r < q, and let Λ_r be the q_r × q_r upper left portion of Λ. We can partition b̆ as b̆ = (b̆_r^T, b̆_f^T)^T, where b̆_r^T is a q_r × 1 vector of penalized coefficients and b̆_f^T is a q_f × 1 vector of unpenalized coefficients, where q_r + q_f = q. Then b̆^T Λb̆ = b̆^T_r Λ_r b̆_r, and it follows that b̆_r ∼ N(0, σ^2_bΛ_r^-1) and b̆_f ∝1.

However, allowing for unconstrained random effect parameters leads to identifiability issues. Therefore, in practice if q_f > 0, we recommend using a normal or Cauchy prior on b̆_f. In particular, b̆_f,l∼ N(0, σ_f) or b̆_f,l∼Cauchy(0, σ_f), l=1,…, q_f with either a diffuse prior on σ_f and constraints to ensure σ_f > 0, or a diffuse prior on log(σ_f) without constraints. The Cauchy prior may be a preferable first choice, as it provides a weaker penalty and is similar to the recommendations of <cit.> for logistic regression. However, in some cases, such as in Section <ref>, it is necessary to use a normal prior.

To further improve the computational efficiency of Monte Carlo sampling methods, we can partition Z̆ into Z̆ = [Z̆_r, Z̆_f] where Z̆_r contains the first q_r columns of Z̆ and Z̆_f contains the remaining q_f columns. We then set b̌_r = Λ_r^-1/2b̆_r and Ž_r = Z̆_r Λ_r^1/2, so that Ž_r b̌_r = Z̆_r b̆_r and b̌_r ∼ N (0, σ^2_b I), which allows for more efficient sampling <cit.>.



§ RELATED WORK 


There are many nonparametric and semiparametric methods for analyzing repeated measures data. For an overview, please see <cit.>. However, most existing methods use an ℓ_2 penalty <cit.>. 

Focusing on the optimization problem, our method puts a generalized lasso penalty <cit.> on the fixed effects and a quadratic penalty on the random effects. Unlike the elastic net <cit.>, we do not mix the ℓ_1 and ℓ_2 penalties on the same parameters, though this could be done in the future.

The additive model with trend filtering developed by <cit.> is similar to our approach. <cit.> optimize 

    θ_1,…, θ_J ∈ℝ^nminimize  1/2y - y̅1 - ∑_j=1^J θ_j _2^2 + λ∑_j=1^J D^(k+1)θ_j _1 
       subject to  1^T θ_j = 0, j=1,…,J.

In contrast to (<ref>), (<ref>) has one smoothing parameter and constrains all smooths to be zero-centered. From Observation <ref>, we see that (<ref>) is equivalent to (<ref>) when there are is J=1 smooth and no random effects, in which case there would be only one smoothing parameter λ and no varying-coefficient smooths.

<cit.> develop the theoretical and computational aspects of additive models with trend filtering, including the extension of the falling factorial basis to additive models. Similar to the B-spline basis, the falling factorial basis allows for linear time multiplication and inversion, which leads to fast computation <cit.>.

When smooths j=1,…,J are expected to have similar degrees of freedom and n is not large enough to require dimension reduction, then (<ref>) with the addition of random effects and the relaxation of the zero-constraints for varying-coefficient smooths may be a viable alternative to (<ref>) that could potentially adapt better to local differences in smoothness because it would have one knot per data point.

While not developed for analyzing repeated measures, the fused lasso additive model (FLAM) <cit.> is also similar to (<ref>). FLAM optimizes the following problem:

    θ_0 ∈ℝ, θ_j ∈ℝ^n, 1 ≤ j ≤ Jminimize  1/2y - θ_0 1 - ∑_j=1^J θ_j _2^2 + αλ∑_j=1^J D^(1)θ_j _1 + (1-α) λ∑_j=1^J θ_j _2

where 0 ≤α≤ 1 specifies the balance between fitting piecewise constant functions (α = 1) and inducing sparsity on the selected smooths (α = 0). From Observation <ref>, we see that (<ref>) is equivalent to our model (<ref>) when: α = 1, there is J=1 smooth, our design matrix has p=n columns, our B-spline bases have appropriately chosen knots, and our model has no random effects. As <cit.> show, FLAM can be a very useful method for modeling additive phenomenon, and as with the fused lasso <cit.>, jumps in the piecewise linear fits have the advantage of being interpretable.

We also mention the sparse additive model (SpAM) <cit.> and sparse partially linear additive model (SPLAM) <cit.>. SpAM fits an additive model and uses a group lasso penalty <cit.> to induce sparsity on the number of active smooths. SPLAM fits a partially linear additive model and uses a hierarchical group lasso penalty <cit.> to induce sparsity in the selected predictors and to control the number of nonlinear features.

One notable difference between our model and that of 
<cit.>, as well as FLAM, SpAM, and SPLAM, is that we allow for multiple smoothing parameters. In our applied experience with additive models and standard ℓ_2 penalties, we have found that in practice it can be important to allow for multiple smoothing parameters, particularly when the quantities of interest are the individual smooths as opposed to the overall prediction. This is equivalent to allowing each smooth to have different variance. However, this flexibility comes at a cost: estimating multiple smoothing parameters is currently the greatest challenge in fitting our proposed model. Perhaps due in part to these computational difficulties, several other authors also assume a single smoothing parameter in high-dimensional additive models <cit.>.

There are fast and stable methods for fitting multiple smoothing parameters for ℓ_2 penalties paired with exponential family and quasilikelihood loss functions, notably the work of <cit.> using generalized cross validation (GCV) and <cit.> using restricted maximum likelihood. Furthermore, <cit.> extend these methods to larger datasets, and <cit.> extend these methods to likelihoods outside the exponential family and quasilikelihood form. However, similarly computationally efficient methods do not yet exist for fitting multiple smoothing parameters for ℓ_1 penalties.

In addition to allowing for multiple smoothing parameters, we also propose approximate inferential methods, which is not typically provided for ℓ_1 penalized models. <cit.>, <cit.>, <cit.>, and <cit.> focus on prediction and provide bounds on the prediction risk and related quantities. These are important results, and we think that distributional results for individual parameters and smooths will also be useful to practitioners.

We also note that <cit.> and <cit.> discuss a variant of P-splines for quantile regression, in which the ℓ_1 norm is used in both the loss and penalty function. However, we are not aware of existing P-spline methods that combine an ℓ_1 penalty with an ℓ_2 loss function.



§ POINT ESTIMATION 




 §.§ Regression parameters and random effects 


To fit (<ref>), we use the alternating direction method of multipliers (ADMM) <cit.>. ADMM has the advantage of being scalable to large datasets. To formulate (<ref>) for ADMM, we introduce constraint terms w_j and re-write the optimization problem as

    minimize      1/2y - β_0 1 - ∑_j=1^J F_j β_j -Zb_2^2 + ∑_j=1^J λ_j w_j _1 + τ/2b^T S b
    subject to       D_j β_j - w_j = 0,   j = 1,…, J


The augmented Lagrangian in scaled form (using u to denote the scaled dual variable) is

    L_ρ(β, b, w, u)    ∝1/2y - β_0 1 - ∑_j F_j β_j - Z b_2^2 + ∑_j λ_j w_j _1 
          + ρ/2∑_j  D_j β_j - w_j + u_j _2^2 + τ/2b^T S b

where ρ > 0 is the penalty parameter. The dimensions are 
y∈ℝ^n × 1,
β_0 ∈ℝ,
F_j ∈ℝ^n × p'_j,
β_j ∈ℝ^p_j' × 1,
Z ∈ℝ^n × q,
b∈ℝ^q × 1,
D_j ∈ℝ^(p_j-k_j-1) × p'_j,
w_j ∈ℝ^(p_j-k_j-1) × 1,
u_j ∈ℝ^(p_j-k_j-1) × 1, and 
S ∈ℝ^q × q,
where p'_j = p_j - 1 if j ∈ℰ (non-varying coefficient smooths) and p'_j = p_j if j ∈ℰ̅ (varying coefficient smooths).

ADMM is an iterative algorithm, and we re-estimate the parameters for updates m=1,2,… until convergence.[We use m to denote the iteration of the ADMM algorithm. This is unrelated to our use of m in Section <ref> to denote the order of the B-spline basis.] It is straightforward to derive the m+1 updates <cit.>:

    β_0^m+1   = 1/n1^T ( y - ∑_j F_j β_j^m - Z b^m ) 
    β^m+1_j    := _β_j L_ρ(β_0^m+1, β_j, β^m+1_l<j, β^m_l>j, b^m, w^m, u^m) 
       = (F_j^T F_j + ρ D_j^T D_j )^-1(F_j^T y^(j,m) + ρ D_j^T(w^m_j - u^m_j) ) 
    b^m+1   := _b L_ρ(β_j=1,…, J^m+1, b, w^m, u^m) 
       = (Z^T Z + τ S)^-1Z^T(y - β_0^m+11 - ∑_j F_j β^m+1_j) 
    w_j^m+1   := _w_j L_ρ(β^m+1_j=1,…,J, b^m+1, w_j, u^m) 
       = ψ_λ_j/ρ(D_j β_j^m+1 + u_j^m) 
    u_j^m+1   := u_j^m + D_j β_j^m+1 - w_j^m+1

where y^(j,m) = y - β_0^m+11 - ∑_l < j F_l β^m+1_l - ∑_l > j F_j β^m_l - Zb^m and ψ_λ/ρ is the element-wise soft thresholding operator, where for a single scalar element x

    ψ_λ/ρ(x) = 
    
    x - λ/ρ    x > λ/ρ
    
    0     |x| ≤λ / ρ
    
    x + λ / ρ    x < -λ / ρ

To initialize the algorithm, we set β_0 := y̅, b := 0, and β_j := 0, w_j := 0, and u_j := 0, for j = 1,…, J.

As an alternative to the closed-form update (<ref>) for the random effects, it is also possible to update the random effects via a linear mixed effects (LME) model that is embedded into the ADMM algorithm. In particular, an LME model is fit to the residuals y - β_0^m+11 - ∑_j F_j β^m+1_j, and b^m+1 are updated as the best linear unbiased predictors (BLUPs). This update occurs at each step of the ADMM algorithm and replaces the update given by (<ref>). The LME update has the additional benefit of simultaneously estimating the variance of the random effects σ^2_b. In simulations, we have found that using an LME update leads to more accurate estimates of σ^2_b, which is important for subsequent estimates of degrees of freedom and confidence intervals.

For stopping criteria, we use the primal and dual residuals (r^m and s^m, respectively):

    r^m    = [ D_1 β_1^m - w_1^m;                 ⋮; D_J β_J^m - w_J^m ]∈ℝ^(p-k - J ) × 1
    
    s^m    = -ρ[ D_1^T (w_1^m - w_1^m-1);                       ⋮; D_J^T (w_J^m - w_J^m-1) ]∈ℝ^p × 1

where k=∑_j=1^J k_j, p = ∑_j=1^J p_j - |ℰ|, and |ℰ| is the cardinality of ℰ.

Following the guidance of <cit.>, we stop when r^m_2 ≤ϵ^pri and s^m_2 ≤ϵ^dual, where

    ϵ^pri   = ϵ^abs√(p - k - J) + ϵ^relmax{
    D_1 β_1^m 
    ⋮
    
    D_J β_J^m
    _2,
    w^m_1 
    ⋮
    w^m_J
    _2
    }
    ϵ^dual   = ϵ^abs√(p) + ϵ^relρ
    D_1^T u_1^m 
    ⋮
    
    D_J^T u_J^m
    _2.

By default, we set ϵ^rel = ϵ^abs = 10^-4 and the maximum number of iterations at 1,000.



 §.§ Smoothing parameters 


To estimate λ_1,…,λ_J we compute cross validation (CV) error for a path of values one smoothing parameter at a time. In the CV, we split the sample at the subject level, as opposed to individual observations, and ensure that there are at least two subjects in each fold per unique combination of factor covariates. First, we estimate a path for τ with λ_1,…, λ_J set to 0. Then we fix τ at the value that minimizes CV error and compute a path for λ_1, setting it to the value that minimizes CV error, and so on.

We fit a path for each λ_j from λ_j^max to 10^-5λ_j^max evenly spaced on the log scale, where λ_j^max is the smallest value at which D_j β_j = 0. As shown in Appendix <ref>,

λ_j^max = (D_j D_j^T)^-1 D_j (F_j^T F_j)^-1 F_j^T r_j_∞, where r_j = y - β_0 1 - ∑_ℓ j F_ℓβ_ℓ  - Z b are the j^th partial residuals and for a vector a, a_∞ = max_j |a_j|.


We also use warm starts, passing starting values separately for each fold, though warm starts appear to be minimally beneficial with ADMM. We set ρ = min( max ( λ_1, …, λ_J), c) at each iteration for some constant c>0 (e.g. c=5). When the number of smooths J is small (e.g. J ≤ 2) a grid search is also feasible.

To estimate τ, we can either use CV and the close-form update given by (<ref>), or an LME update that is embedded in the ADMM algorithm, as described in Section <ref>. In simulations, we have found that the overall computation time to estimate the smoothing parameters is greater when using the LME update, and that the estimates of λ_1, …, λ_J do not appear sensitive to updates for b. However, the final estimates of σ^2_b, and consequently the width of confidence intervals can be improved by using the LME update. Consequently, we recommend using cross validation to estimate τ for the purposes of then estimating λ_1,…,λ_J, but using an LME update when estimating the final model.

With both the closed-form and LME update, we cannot use the training sample to estimate the random effect parameters b for the test sample, because these parameters are subject-specific and the test subjects are not included in the training sample. Instead, we use the training sample to obtain estimates for the fixed effect parameters β_0, β_j, j=1,…,J and then use the test sample to estimate the random effects.

To make our approach clear, we first fix notation. Let 𝒯^r ⊆{1,…,n} be the row indices for the observations in the test sample for both the fixed and random effect design matrices F_j, j=1,…,J, and Z. Also, let 𝒯^c ⊆{1,…,q} be the column indices of Z for observations in the test sample, and let 𝒯 = (𝒯^r, 𝒯^c) be the tuple of row and column indices designating the test sample. Let matrices F_j,𝒯 and F_j,-𝒯 be matrix F_j with only rows indexed by 𝒯^r retained and removed, respectively. Similarly, let matrices Z_𝒯 and Z_-𝒯 be matrix Z with only rows and columns indexed by 𝒯^r and 𝒯^c, respectively, retained and removed, respectively. Let matrices S_𝒯 and S_-𝒯 be matrix S with only rows and columns indexed by 𝒯^c retained and removed, respectively. Also, let y_𝒯 and y_-𝒯 be vector y with elements indexed by 𝒯^r retained and removed, respectively.

We obtain out-of-sample marginal estimates as μ̂_𝒯 = β̂_0 1 + ∑_j=1^J F_j,𝒯β̂_j, where β̂_0 and β̂_j, j=1,…,J are estimated with y_-𝒯, F_j, -𝒯, and Z_-𝒯. If using the closed-form update (<ref>), we estimate subject-specific random effects as b̂_𝒯 = ( Z_𝒯^T Z_𝒯^T + τ S_𝒯)^-1 Z_𝒯^T (y_𝒯 - μ̂_𝒯) and obtain the out-of-sample prediction residuals as r_𝒯 = y_𝒯 - μ̂_𝒯 - Z_𝒯b̂_𝒯. Letting 𝒯_k be the tuple of indices for test sample (fold) k=1,…,K, we obtain the CV error as ∑_k=1^K r_𝒯_k_2^2.



§ DEGREES OF FREEDOM 


In this section, we obtain the degrees of freedom, with the primary goal of estimating variance (see Section <ref>). However, we note that degrees of freedom does not always align with a model's complexity in terms of its tendency to overfit the data <cit.>.

In each of the approaches described in this section, the degrees of freedom (df) is a function of the smoothing parameters λ_1,…λ_J and τ. We always obtain the fixed effects smoothing parameters λ_1,…,λ_J from CV, but when using an LME update for the random effects b as described in Sections <ref> and <ref>, we do not directly obtain τ. Consequently, we cannot directly apply the results in this section to estimate df. However, from (<ref>), we have that τ = σ^2_b / σ^2_ϵ. Writing df = df(τ), and letting r = y - ∑_j=1^J F_j β̂_j - Z b̂ be an n × 1 vector of residuals and σ̂^2_ϵ = r_2^2 / (n - df(τ)) be an estimate of variance, we have that
 

    τ̂ = σ̂^2_b/σ̂^2_ϵ = σ̂^2_b/r_2^2(n - df(τ̂) ).

Therefore, letting

    ψ(τ) = τ - σ̂^2_b/r_2^2(n - df(τ) ),

we numerically solve for τ̂ such that ψ(τ̂) = 0 and set df = df(τ̂).



 §.§ Stein's method 

 
Let g(y) = ŷ, where g : ℝ^n →ℝ^n is the model fitting procedure. For y∼ N(μ, σ^2I), the degrees of freedom is defined as <cit.>

    df = 1/σ^2∑_i=1^n  (g_i(y), y_i).

As <cit.> notes, (<ref>) is motivated by the fact that the risk Risk(g) = 𝔼g(y) - μ_2^2 can be decomposed as

    Risk(g) = 𝔼 g(y) - y_2^2 - n σ^2 + 2 ∑_i=1^n (g_i(y), y_i).

Therefore, the degrees of freedom (<ref>) corresponds to the difference between risk and expected training error. Furthermore, if g is continuous and weakly differentiable, then df = 𝔼[∇· g(y)] <cit.> where ∇· g = ∑_i=1^n ∂ g_i / ∂ y_i is the divergence of g. Therefore, an unbiased estimate of df (also used in Stein's unbiased risk estimate <cit.>) is 

    d̂f̂ = ∑_i=1^n ∂ g_i / ∂ y_i.


To obtain an estimate of degrees of freedom, we transform the generalized lasso component of our model to standard form, similar to the approach of <cit.>. To do so, we use the following matrices described by <cit.>. Let

    D̃^*_j = [   D̃^(0)_j,1;            ⋮; D̃^(k_j)_j,1; D̃^(k_j+1)_j ]∈ℝ^p_j × p_j

be an augmented finite difference matrix, where D̃^(i)_j,1 is the first row of the finite difference matrix D̃^(i)_j, and D̃^(0)_j = I_p_j is the identity matrix.

As shown by <cit.>, the inverse of D̃^*_j is given by M_j = M_j^(0) M_j^(1)⋯ M_j^(k) where[We denote the inverse matrix as M_j. This is unrelated to our use of M in Section <ref> to denote the order of the B-spline basis.]

    M_j^(i) = 
    [                     I_i;                         L_(p_j - i) × (p_j - i) ]∈ℝ^p_j × p_j,

where L_(p_j - i) × (p_j - i) is the (p_j - i) × (p_j - i) lower diagonal matrix of 1s.

Assuming our outcome y is centered, so that β_0 = y(0) = 0, and letting V_j = F̃_j M_j, D_j^* = D̃_j^* Q_j for j ∈ℰ and D_j^* = D̃_j^* for j ∈ℰ̅, and α_j = D^*_j β_j, we can write the penalized log likelihood (<ref>) as 

    l_pen = 1/2y - ∑_j V_j α_j - Z b_2^2 + ∑_j=1^J λ_j ∑_l=k_j+2^p_j |α_jl | + 1/2τb^T S b.


To avoid difficulties later differentiating with respect to the ℓ_1 norm, we remove the non-active ℓ_1 penalized coefficients from (<ref>). We also form the concatenated design matrix V = [V_1,…,V_J] and will need to index the active set of V. To these ends, let 𝒜_j = {l ∈{k_j+2, …, p'_j} : α̂_j,l 0 } be the active set of the penalized coefficients for smooth j, and let 𝒜_j^* = {1,…,k_j + 1 }∪𝒜_j be the active set for smooth j augmented with the unpenalized coefficients. Also, for a set 𝒜_j and constant c ∈ℝ, let 𝒜_j + c = {i + c : i ∈𝒜_j } be the set of elements in 𝒜_j shifted by c. Now let 𝒜^* = ⋃_j=1^J (𝒜^*_j + ∑_l =0^j-1 p'_l) be the augmented active set of V, where p'_0 = 0 and p'_j, j=1,…,J are the number of columns in V_j (equivalently F_j). Finally, let V_𝒜^* be matrix V subset to retain only those columns indexed by 𝒜^*. Similarly, let α̂ = (α̂_1^T, …, α̂_J^T)^T be the concatenated vector of estimated coefficients, and let α̂_𝒜^* be vector α̂ subset to retain only elements indexed by 𝒜^*. Then we can write the estimated penalized loss (<ref>) as

    l̂_pen = 1/2y - [V_𝒜^*, Z] [ α̂_𝒜^*;     b̂ ]_2^2 + ∑_j=1^J λ_j ∑_l=k_j+2^p_j |α̂_jl | + 1/2τb̂^T S b̂


Taking the derivative of (<ref>) and keeping in mind that the first k_j + 1 elements of each α̂_j are unpenalized and |α̂_jl| > 0 for all l ∈𝒜_j, we have

    0_(|𝒜^*| + q) × 1 = ∂ l_pen/∂ (α̂_𝒜^*^T, b̂^T)^T = [ V_𝒜^*^T;     Z^T ]([V_𝒜^*, Z] [ α̂_𝒜^*;     b̂ ] - y)
    + 
    [      η; τ S b̂ ]

where 

    η = 
    [                           0_k_1+1;                 λ_1  sign(α̂_𝒜_1);                           0_k_2+1;           λ_2  sign(α̂_𝒜_2 + p_1);                                 ⋮;                           0_k_J+1; λ_J  sign(α̂_𝒜_J + ∑_j=1^J-1 p_j) ],

0_k_j + 1 is a (k_j + 1) × 1 vector of zeros, and the sign operator is taken element-wise. 

From <cit.>, we know that within a small neighborhood of y, the active set 𝒜 and the sign of the fitted terms α̂_𝒜 are constant with respect to y except for y in a set of measure zero. Therefore, ∂η / ∂y = 0_|𝒜^*| × n, where 0_|𝒜^*| × n is an |𝒜^*| × n matrix of zeros and |𝒜^*| is the cardinality of 𝒜^*. Then taking the derivative of (<ref>) with respect to y, we have 

    0_(|𝒜^*| + q) × n = ∂^2 l_pen/∂ (α̂_𝒜^*^T, b̂^T)^T ∂y   = 
    [ V_𝒜^*^T;     Z^T ] 
    [V_𝒜^*, Z] 
    [ ∂α̂_𝒜^* / ∂y;     ∂b̂ / ∂y ] - 
    [ V_𝒜^*^T;     Z^T ]
         + 
    [    0_|𝒜^*| × n; τ S (∂b̂ / ∂y) ].

Solving for the derivatives of the estimated coefficients, we have

    [ ∂α̂_𝒜^* / ∂y;     ∂b̂ / ∂y ]
    =
    (
    [ V_𝒜^*^T;     Z^T ] 
    [V_𝒜^*, Z] 
    +
    [ 0_|𝒜^*| × |𝒜^*|     0_|𝒜^*| × q;     0_q × |𝒜^*|             τ S ])^-1[ V_𝒜^*^T;     Z^T ].

Now let A = [V_𝒜^*, Z] and 

    Ω = [ 0_|𝒜^*| × |𝒜^*|     0_|𝒜^*| × q;     0_q × |𝒜^*|             τ S ].

Then since
ŷ = A (α̂_𝒜^*^T, b̂^T)^T
we have

    ∂ŷ/∂y   = 
    ∂ŷ/∂ (α̂_𝒜^*^T, b̂^T)^T∂ (α̂_𝒜^*^T, b̂^T)^T/∂y
       =
    A ( A^T A + Ω)^-1 A^T.


From <cit.>, we know that g(y) = ŷ is continuous and weakly differentiable. Also, ∇ g = (∂ŷ / ∂y). Therefore, we can use Stein's formula (<ref>) to estimate the degrees of freedom as

    d̂f̂ = 1 + (A ( A^T A + Ω)^-1 A^T) = 1 + ( (A^T A + Ω)^-1 A^T A ),

where we add 1 for the intercept. We note that this result is similar to the degrees of freedom for the elastic net <cit.> as well as for FLAM <cit.>.

To obtain degrees of freedom for individual smooths j=1,…, J, let E_j be an (|𝒜^*|+q) × (|𝒜^*|+q) matrix with 1s on the diagonal positions indexed by 𝒜^*_j + ∑_l=0^j-1 |𝒜^*_l| and zero elsewhere, where |𝒜^*_j| is the cardinality of 𝒜^*_j and 𝒜^*_0 = ∅. Also, let f̂_j = V_j α̂_j be the estimate of the j^th smooth. Then as <cit.> note, f̂_j = A E_j(A^TA + Ω)^-1 A^T y. Therefore, 

    d̂f̂_j = (A E_j(A^T A + Ω)^-1 A^T ) = (E_j(A^T A + Ω)^-1 A^T A ).

In other words, the degrees of freedom for smooth j is the sum of the diagonal elements of (A^T A + Ω)^-1 A^T A indexed by 𝒜_j^* + ∑_l=0^j-1 |𝒜^*_l|.

We note that when using the ADMM algorithm, or most likely any proximal algorithm, the fitted D_j β̂_j, or equivalently α̂_j, will typically have several very small non-zero values, but will not typically be sparse. However, the vector ŵ_j is sparse, where in the ADMM algorithm we constrain w_j = D_j β_j. Therefore, in practice we use w_j to obtain the active set 𝒜_j.



 §.§ Stable and fast approximations 


In some cases, such as the application in Section <ref>, the estimates based on Stein's method (<ref>) and (<ref>) cannot be computed due to numerical instability. In this section, we propose alternatives that are more numerically stable and which are also more computationally efficient.



  §.§.§ Based on restricted derivatives


In this approach, we take derivatives of the fitted values restricted to individual smooths. In particular, from Section <ref>, we see that

    ∂ŷ/∂α̂_𝒜^*_j∂α̂_𝒜^*_j/∂y   = V_𝒜_j^* (V_𝒜_j^*^T V_𝒜_j^*)^-1 V_𝒜_j^*^T 
    ∂ŷ/∂b̂∂b̂/∂y   = Z (Z^T Z + τ S)^-1 Z^T.

We can then approximate the degrees of freedom for each individual smooth and the random effects by

    d̃f̃_j = 
    ((V_𝒜_j^*^T V_𝒜_j^*)^-1 V_𝒜_j^*^T V_𝒜_j^*)     j=1,…,J 
    ((Z^T Z + τ S)^-1 Z^T Z )     j=J+1

We estimate the overall degrees of freedom as

    d̃f̃ = 1 + ∑_j=1^J+1d̃f̃_j

where we add 1 for the intercept.

This approach is similar to one described by <cit.>, though in a different context and for a different purpose. In particular, whereas we use this approach to approximate the degrees of freedom after fitting the model, <cit.> use it to set the degrees of freedom before fitting the model in the context of ℓ_2 penalized loss functions.



  §.§.§ Based on ADMM constraint parameters


In this approach, we propose estimates of degrees of freedom specific to the ADMM algorithm. As in the previous section, this approach is based on estimates for the individual smooths. Consider the model with J=1 smooth, no random effects, and centered y:

    y - F β_2^2 + λ D β_1.

Suppose we make the centering constraints described Section <ref>, i.e. we set F = F̃ Q and D = D̃^(k+1) Q for an n × p design matrix F̃, a k+1 order finite difference matrix D^(k+1), and an orthonormal p × (p-1) matrix Q. Let 𝒜 = {l ∈{1,…,p - k-1 }: (D β̂)_l  0} be the active set, and let |𝒜| be its cardinality. In our context, we expect the design matrices F to be full rank, in which case Theorem 3 of <cit.> (see the first Remark) states that the degrees of freedom is given by df = 𝔼[nullity(D_-𝒜)]. Here, nullity(D) is the dimension of the null space of matrix D, and D_-𝒜 is matrix D with rows indexed by 𝒜 removed. Now, D has dimensions (p-k-1) × (p-1), and we can see by inspection that for all k<p-1 the columns of D are linearly independent. Therefore, the rank of D_-𝒜 is equal to the number of rows p - k - 1 - |𝒜|, and the nullity is equal to the number of columns p-1 minus the number of rows. This gives d̂f̂ = nullity(D_-𝒜) = k + |𝒜| for centered smooths, i.e. the number of non-zero elements of D β̂ plus one less than the order of the difference penalty. This is similar to the result for ℓ_1 trend filtering, but we have lost one degree of freedom due to the constraint that 1^T F̃β̃ = 0. For uncentered smooths, D has dimensions (p-k-1) × p, which gives d̂f̂ = nullity(D_-𝒜)) = k + 1 + |𝒜|.

As before, we note that in the ADMM algorithm, D β̂ will not generally be sparse, as ADMM is a proximal algorithm. However, the corresponding w is sparse, where in the optimization problem we constrain D β = w. Now considering a model with smooths j=1,…,J, a numerically stable and fast alternative to (<ref>) is given by

    d̃f̃_j^ADMM   = 1[j ∈ℰ̅] + k_j + ∑_l=1^p - k - 11[w_jl 0].

where ℰ̅ indexes the un-centered smooths and 1 is an indicator variable. We then combine (<ref>) with the restricted derivative approximation for the degrees of freedom of the random effects given above to obtain the overall degrees of freedom

    d̃f̃^ADMM = 1 + ∑_j=1^J d̃f̃^ADMM_j + ((Z^T Z + τ S)^-1 Z^T Z ),

where we add 1 for the intercept.



 §.§ Ridge approximation


Let U = [F_1, …, F_J, Z] be the concatenated design matrix of fixed and random effects and

    Ω^ridge = [ λ_1 D_1^T D_1;                           ⋱;               λ_J D_J^T D_J;                                       τ S ]

be the penalty matrix. Then the hat matrix from the linear smoother approximation (see Section <ref>) is given by H = U (U^T U + Ω^ridge)^-1 U^T. Similar to before, we can get the overall degrees of freedom as

    d̂f̂^ridge = 1 + ((U^T U + Ω^ridge)^-1 U^T U ),

where we add 1 for the intercept. To obtain degrees of freedom for individual smooths j=1,…, J, let E_j be a (p+q) × (p+q) matrix with 1s on the diagonal positions indexed by the columns of F_j and zero elsewhere. Also, let f̂_j = F_j β̂_j be the estimate of the j^th smooth. Then the ridge approximation for smooth j is given by f̂_j ≈ U E_j(U^T U + Ω^ridge)^-1 U^T y. Therefore, 

    d̂f̂^ridge_j = ( E_j(U^T U + Ω^ridge)^-1 U^T U )


Similar to before, we also propose stable and fast approximations to the ridge estimate of degrees of freedom based on restricted derivatives. In particular, let

    d̃f̃_j^ridge = 
    ((F_j^T F_j + λ_j D_j^T D_j)^-1 F_j^T F )     j=1,…,J 
    ((Z^T Z + τ S)^-1 Z^T Z )     j=J+1

Then we can estimate the overall degrees of freedom as

    d̃f̃^ridge = 1 + ∑_j=1^J+1d̃f̃_j^ridge

where we add 1 for the intercept.

As noted above, this approach is similar to one described by <cit.>, though for a different purpose. Whereas we use this approach to obtain the degrees of freedom after fitting the model, <cit.> use it to set the degrees of freedom before fitting the model.



§ APPROXIMATE INFERENCE 


In this section, we discuss approximate inferential methods based on ridge approximations to the ℓ_1 penalized fit and conditional on the smoothing parameters λ_j, j=1,…,J and τ. We use the ADMM algorithm to analyze the approximation. In particular, we note that we can write the ADMM update for β_j as

    β_j^m+1 = (F_j^T F_j + ρ D_j^T D_j )^-1 F_j^T y^(j,m) + δ_j^m

where δ_j^m = ρ(F_j^T F_j + ρ D_j^T D_j)^-1F^T_j D^T_j(w^m_j - u^m_j) and y^(j,m) = y - β_0^m+1 - ∑_l<j F_l β_l^m+1 - ∑_l>j F_l β_l^m - Z b^m. As we note in Observation <ref>, δ_j loosely represents the difference in the estimate of β_j obtained with the ℓ_1 and ℓ_2 penalties.


With the ℓ_1 penalty, i.e. D_j β_j_1, in general δ_j^m0. However, with the ℓ_2 penalty, i.e. D_j β_j_2^2, and λ_j = ρ, we have δ_j^m = 0.




Similar to the ridge update for b, if we changed λ_j D_j β_j_1 to (λ_j / 2) D_j β_j_2^2 in (<ref>) we could remove the w_j term and the constraint that D_j β_j^m = w_j from (<ref>) to obtain the ridge update β_j^m+1 = (F_j^T F_j + λ_j D_j^T D_j )^-1 F_j^T y^(j,m). Then since we assumed λ_j = ρ, we have β_j^m+1 = (F_j^T F_j + ρ D_j^T D_j )^-1 F_j^T y^(j,m). By comparison with (<ref>), we see that δ_j^m = 0.


Observation <ref> motivates our approximate inferential strategy. Letting f̂_j be the j^th fitted smooth, and letting y^(j) = y - β̂_0 - ∑_l  j F_l β̂_l - Z b̂, we have

    f̂_j = F_j β̂_j    = F_j(F_j^T F_j + ρ D_j^T D_j)^-1 F_j^T y^(j) + F_j δ̂_j 
       ≈ F_j(F_j^T F_j + ρ D_j^T D_j)^-1 F_j^T y^(j)       (assuming  F_j δ̂_j ≈0) 
       ≈ F_j(F_j^T F_j + λ_j D_j^T D_j)^-1 F_j^T y^(j)       (assuming λ_j ≈ρ) 
       =H_j y^(j)

where H_j = F_j(F_j^T F_j + λ_j D_j^T D_j)^-1 F_j^T. We obtain confidence intervals for the linear smoother (<ref>) centered around the estimated fit (<ref>), ignore F_j δ_j when estimating variance, and assume λ_j ≈ρ. We also condition on the smoothing parameters λ_1,…,λ_J and τ.

Figure <ref> gives a visual demonstration of the approximation for the simulation presented in Section <ref> and the application shown in Section <ref>. As seen in Figure <ref>, in these examples the ℓ_1 fit and ridge approximation are very similar. If this holds in general, then this would suggest that 1) the approximate inferential procedures we propose might have reliable coverage probabilities, and 2) there may be minimal practical advantage to using an ℓ_1 penalty instead of the standard ℓ_2 penalty. However, as shown in Section <ref>, the ℓ_1 penalty appears to perform noticeably better in certain situations, including the detection of change points.



Before presenting the confidence bands in greater detail, we discuss our approach for estimating variance in Section <ref>, which we then use to form confidence bands in Section <ref>.



 §.§ Variance 


Let r = y - ∑_j=1^J F_j β̂_j - Z b̂ be an n × 1 vector of residuals. We estimate the overall variance as σ̂^2_ϵ = r^2_2 / d̂f̂_resid, where d̂f̂_resid is the residual degrees of freedom. When possible, we use the estimate based on Stein's method (<ref>) and set d̂f̂_resid = n - d̂f̂. If Stein's method is not numerically stable, then we use the restricted derivatives approximation (<ref>) and set d̂f̂_resid = n - d̃f̃. As another alternative, we could also use the ADMM approximation and set d̂f̂_resid = n - d̃f̃^ADMM.



 §.§ Confidence bands 


In this section, we obtain confidence bands for typical subjects, i.e. for subjects for whom b_i = 0. Since we assume a normal outcome, this is equivalent to the marginal population level response.



  §.§.§ Frequentist confidence bands


Ignoring the distribution on D_j β_j and treating β_l, l  j as fixed, y^(j) is normal with variance Var(y^(j)) = σ^2_ϵ I_n + σ^2_b Z S^+ Z^T, where S^+ is the Moore-Penrose generalized inverse of matrix S (as noted in Section <ref>, S may not be positive definite). Therefore, Var(f̂_j) ≈ H_j Var(y^(j)) H_j^T where Var(y^(j)) is an n × n estimate of Var(y^(j)) with σ̂^2_ϵ and σ̂^2_b plugged in for σ^2_ϵ and σ^2_b respectively, and f̂_j ·∼ N(f̂_j, H_j Var(y^(j)) H_j^T). The estimated variance of the fit at a single point x, which we denote as Var(f̂_j(x)), is the corresponding diagonal element of H_j Var(y^(j)) H_j^T. Therefore, asymptotic pointwise 1-α confidence bands take the form f̂_j(x) ± z_1-α/2√(Var(f̂_j(x))) where Φ(z_a) = a and Φ is the standard normal CDF, e.g. z_1-α/2 = 1.96 for α = 0.05.

For the purposes of interpretation, we include the intercept term in the confidence band for the j=1 smooth, but not for the remaining smooths.



  §.§.§ Bayesian credible bands 


Many authors, including <cit.>, recommend using Bayesian confidence bands for nonparametric and semiparametric models, because the point estimates are themselves biased. While Bayesian credible bands do not remedy the bias, they are self consistent.

To this end, we replace the element-wise Laplace prior with the (generally improper) joint normal prior that is equivalent to the standard ℓ_2 penalty: β_j ∼ N (0, ( λ_j D_j^T D_j)^-1). This leads to the posterior

    β_j | y·∼ N (β̂_j, (W_jF^T_j Var(y^(j))^-1 F_j + λ_j D^T_j D_j)^-1).

We can then form simultaneous Bayesian credible bands for f_j | y by simulating from the posterior (<ref>) and taking quantiles from  F_j β_j^b, b = 1,…,B. Alternatively, for a faster approximation we use frequentist confidence bands with F_j W_j^-1 F^T_j in place of H_j Var(y^(j)) H_j^T. In practice, we have found the simultaneous credible bands and the faster approximation to be nearly indistinguishable.[It appears that the latter (faster) method is the default in the  package <cit.>. As in , we only need to compute the diagonal elements of F_j W_j^-1 F^T_j as ((F_j W_j^-1) ∘ F_j), where ∘ is the Hadamard (element-wise) product.]

As before, for the purposes of interpretation, we include the intercept term in the credible band for the j=1 smooth, but not for the remaining smooths.



§ SIMULATION 


We simulated data from a piecewise linear mean curve as shown in Figure <ref>. Each subject had a random intercept and is observed over only a portion of the domain. There are 50 subjects, each with between 4 and 14 measurements (450 total observations). The random intercepts were normally distributed with variance σ_b^2 = 1, and the overall noise was normally distributed with variance σ_ϵ^2 = 0.01.



In all models, we used order 2 (degree 1) B-splines with p=21 basis functions.



 §.§ Frequentist estimation


We fit models with J=1 smooth term and random intercepts. To obtain estimates for the ℓ_1 penalized model, we used ADMM and 5-fold CV to minimize

    β_0 ∈ℝ, β∈ℝ^p-1, b∈ℝ^Nminimize1/2y - β_0 1 - F β - Z b_2^2 + λD^(2)β_1 + τb^T b.

where Z_il = 1 if observation i belongs to subject l and zero otherwise. As noted above, we used order 2 (degree 1) B-splines with p=21 basis functions, i.e. F ∈ℝ^n × (p-1) where n=450 and p=21. After estimating λ and τ via CV, we used LME updates to estimate σ^2_b and b in the final model. We also fit an equivalent model with an ℓ_2 penalty using the  package <cit.>, i.e. with (λ / 2) D^(2)β_2^2 in place of λD^(2)β_1 in (<ref>). Figure <ref> shows the marginal mean with 95% credible intervals, and Figure <ref> shows the subject-specific predicted curves.





As seen in Figures <ref> and <ref>, the results from the ℓ_1 and ℓ_2 penalized models are very similar. However, the ℓ_1 penalized model does slightly better at identifying the change points and the line segments. We explore this further in Section <ref>.

Table <ref> compares the degrees of freedom and variance estimates from the ℓ_1 penalized fit against those from the ℓ_2 penalized fit. From Table <ref>, we see that the ridge degrees of freedom d̂f̂^ridge appears reasonable, as it is near the estimate for the ℓ_2 penalized model. The true degrees of freedom d̂f̂ also seems reasonable. Ideally, the degrees of freedom for the ℓ_1 penalized fit should equal six, as there are four change points and we are using a second order difference penalty (see Section <ref>).



Table <ref> compares the different estimates of degrees of freedom. In this simulation, the degrees of freedom based on the ridge approximation is larger than that from Stein's formula, and the approximations based on restricted derivatives are equal or near the estimate with Stein's formula.





 §.§ Bayesian estimation


We modeled the data as y | b = β_0 1 + F β + b + ϵ where

    ϵ   ∼ N(0, σ^2_ϵ I) 
    b   ∼ N(0, σ^2_b I) 
    
    D^(2)β   ∼Laplace(0, σ^2_λ I) 
    
    p(σ_ϵ)    ∝ 1 
    
    p(σ_b)    ∝ 1 
    
    p(log(σ_λ))    ∝ 1.

We also fit models with normal and diffuse priors for D^(2)β.

We fit all models with  <cit.>, each with four chains of 2,000 iterations with the first 1,000 iterations of each chain used as warmup. The MCMC chains, not shown, appeared to be reasonably well mixing and stationary, and had R̂ values under 1.1 <cit.>.[As described by <cit.>, for each scalar parameter, R̂ is the square root of the ratio of the marginal posterior variance (a weighed average of between- and within-chain variances) to the mean within-chain variance. As the number of iterations in the MCMC chains goes to infinity, R̂ converges to 1 from above. Consequently, R̂ can be interpreted as a scale reduction factor, and <cit.> recommend ensuring that R̂<1 for all parameters.] Figure <ref> shows the marginal mean with 95% credible intervals, and Figure <ref> shows point estimates.





As seen in Figures <ref> and <ref>, all models performed well and gave similar fits as above. Similar to before, the Laplace prior appears to better enforce a piece-wise linear fit, particularly around x=0.2.



 §.§ Change point detection 


We simulated 1,000 datasets with the same generating mechanism used to produce the data shown in Figure <ref> and measured the performance of the ℓ_1 and ℓ_2 penalized models on two criteria: 1) the number of inflection points found, and 2) the distance between the estimated inflection points and the closest true inflection point. To that end, let 𝒯={τ_1,…,τ_4} be the set of true inflection points, and M = max_x ∈𝒳|f̂”(x)| be the maximum absolute second derivative of the estimated function, where 𝒳={x_1,x_2,…} is the ordered set of unique simulated x values. We approximate f̂” by

    f̂”(x_i) ≈(f̂(x_i+1) -f̂(x_i))/(x_i+1 - x_i) - (f̂(x_i) -f̂(x_i-1))/(x_i - x_i-1)/x_i+1 - x_i.

Then let ℐ = {x ∈𝒳 : |f̂”(x)| ≥ c M } be the set of estimated inflection points, where c ∈ (0,1) is a cutoff value defining how large the second derivative must be to be counted as an inflection point. Also, let n_ℐ = |ℐ| be the number of estimated inflection points, and d̅ = n^-1_ℐ∑_x ∈ℐmin_τ∈𝒯 |x - τ| be the mean absolute deviance of the estimated inflection points.

Figure <ref> shows the results from 1,000 simulated datasets. The ℓ_1 penalized model was better able to 1) find the correct number of inflection points, and 2) determine the location of the inflection points.





 §.§ Coverage probability


We simulated 1,000 datasets with the same generating mechanism used to produce the data shown in Figure <ref> and measured the coverage probability of the approximate Bayesian credible bands described in Section <ref> for the ℓ_1 penalized model, and simultaneous Bayesian credible bands for the ℓ_2 penalized model <cit.>. Figure <ref> shows the coverage probabilities for both approaches. As seen in Figure <ref>, the confidence bands perform similarly and are near the nominal rate over most of the x domain. Both approaches have difficulty maintaining nominal coverage at the edges of the x domain.





§ APPLICATION 




 §.§ Data description and preparation


In this section, we analyze electrodermal activity (EDA) data collected as part of a stress study. In brief, all subjects completed a written questionnaire prior to the study, which categorized the subjects as having either low vigilance or high vigilance personality types. During the study, all participants wore wristbands that measured EDA while undergoing stress-inducing activities, including giving a public speech and performing mental arithmetic in front of an audience. The scientific questions were: 1) Is EDA higher among high vigilance subjects, and 2) when did trends in stress levels change? In this section, we demonstrate how P-splines with an ℓ_1 penalty can address both questions.

The raw EDA data are shown in Figure <ref>. After excluding subjects who had EDA measurements of essentially zero throughout the entire study, we were left with ten high vigilance subjects and seven low vigilance subjects.



To remove the extreme second-by-second fluctuations in EDA, which we believe are artifacts of the measurement device as opposed to real biological signals, we smoothed each curve separately with a Nadaraya–Watson kernel estimator using the  function in . We then thinned the data to reduce computational burden, taking 100 evenly spaced measurements from each subject. Figure <ref> shows the results of this process for a single subject, and Figure <ref> shows the prepared data for all subjects. Because of the limited number of subjects, as well as issues of misalignment in the time series across individuals, the results presented here should be considered as illustrative rather than of full scientific validity.






 §.§ Models 


In all models, we fit the structure

    y_i(x) = β_0 + β_1(x) + 1_high[i] β_2(x) + b_i(x) + ϵ_i(x)

where x represents time in minutes, 1_high[i]=1 if subject i has high vigilance and 1_high[i]=0 if subject i has low vigilance, b_i(x) are random subject-specific curves, and ϵ_i(x) ∼ N(0, σ^2_ϵ). For β_1(x), β_2(x), and b_i(x), we used a fourth order B-spline basis with 31 basis functions each and a second order difference penalty (k=1).

Written in matrix notation, the ℓ_1 penalized model is

    min1/2y - β_0 1 - ∑_j=1^2 F_j β_j - Z b_2^2 + ∑_j=1^2 λ_j D^(2)β_j _1 + b^T S b

where y is a stacked vector for subjects i=1,…,17, F_1 is an n × p design matrix where n=1,700 and p=31, and F_2 = diag( 1_high[i]) F_1 where i is an n × 1 vector of subject IDs. In other words, F_2 is equal to F_1, but with rows corresponding to low vigilance subjects zeroed out. We set

    Z = [  Z_1;         ⋱;      Z_17 ]

where each Z_i is an n_i × 31 random effects design matrix of order 4 B-splines evaluated at the input points for subject i, and 

    S = [  S_1;         ⋱;      S_17 ]

where S_i,jl = ∫ϕ”_ij(t) ϕ”_il(t) dt are smoothing spline penalty matrices. We also mean-centered F_1 as described in Section <ref>, with the corresponding changes in dimensions.

To fit a comparable ℓ_2 penalized model, in which λ_j D^(2)β_j _1 in (<ref>) is replaced with (λ_j/2) D^(2)β_j _2^2, we rotated the random effect design and penalty matrices Z and S as described in Section <ref>. To facilitate the use of existing software, we used a normal prior for the “unpenalized" random effect coefficients, i.e. b̆_f ∼ N(0, σ^2_f I).

We also fit a Bayesian model using the same rotations and equivalent penalties as above. In particular, we modeled the data as y|b = β_0 1 + ∑_j=1^J F_j β_j + Ž_r b̌_r + Z̆_f b̆_f + ϵ where

    b̌_r    ∼ N(0, σ^2_r I) 
    b̆_f    ∼ N(0, σ^2_f I) 
    ( D_j β_j )_l    ∼Laplace(0, a_j)  for  a_j = σ^2_ϵ / (2 λ_j), l=1,…, p_j - k_j - 1, j=1,…,J 
    ϵ   ∼ N ( 0, σ^2_ϵ I_n ).




 §.§ Results 




  §.§.§ Frequentist estimation


We tried to use CV to estimate the smoothing parameters for the ℓ_1 penalized model. However, with only 17 subjects split between two groups, we only did 3-fold CV. CV did not find a visually reasonable fit so we set the tuning parameters by hand.

Figure <ref> shows the estimated marginal mean and 95% credible bands for the ℓ_1 penalized model, and Figure <ref> shows the subject-specific predicted curves for the ℓ_1 penalized model. As seen in Figure <ref>, our model identified a few inflection points, particularly near minutes 40, 50, and 60. From Figure <ref> it appears that the difference in EDA between the low and high vigilance subjects was not statistically significant. Also, as seen in Figure <ref>, the subject-specific predicted curves are shrunk towards the mean, which is expected, because the predicted curves are analogous to BLUPs, although they are not linear smoothers.





Figure <ref> shows the estimated marginal mean and 95% credible bands for the ℓ_2 penalized model, and Figure <ref> shows the subject-specific predicted curves for the ℓ_2 penalized model. The estimate shown in Figure <ref> is similar to that shown in Figure <ref>, though the inflection points are slightly less pronounced in Figure <ref>. The results in Figure <ref> are for the most part substantively the same as those in Figure <ref>; the ℓ_2 penalized model does not show a statistically significant difference between the low and high vigilance subjects, with the possible exception of minutes 45 to 66. As seen in Figure <ref>, the predicted subject-specific curves from the ℓ_2 penalized model are also shrunk towards the mean.





Table <ref> shows the estimated degrees of freedom for the ℓ_1 penalized model. Stein's method d̂f̂ ((<ref>) and (<ref>)) and the ridge approximation d̂f̂^ridge ((<ref>) and (<ref>)) were numerically instable (A^T A + Ω and U^T U + Ω^ridge were computationally singular). Therefore we used the restricted derivative approximation d̃f̃ to estimate the variance, as described in Section <ref>. In the ℓ_2 penalized model, smooth F_1 had 14.2 degrees of freedom, and smooth F_2 had 6.96 degrees of freedom.





  §.§.§ Bayesian estimation


We fit the model described in Section <ref> with an element-wise Laplace prior on D β given by (<ref>). To fit the model, we used  <cit.> with four chains of 5,000 iterations each, with the first 2,500 iterations of each chain used as warmup. The MCMC chains, not shown, appeared to be reasonably well mixing and stationary with R̂ values under 1.1 <cit.>. Figure <ref> shows the marginal means with 95% credible bands, and Figure <ref> shows the subject-specific curves. Similar to the ℓ_2 penalized model, the Bayesian model found a slightly statistically significant difference between low and high vigilance between minutes 42 and 65.







 §.§ Alternative correlation structure 


For comparison, we also fit ℓ_1 and ℓ_2 penalized models with alternative correlation structures similar to that recommended by <cit.>.

For the ℓ_1 penalized model, in place of the correlation structure implied by the penalty matrix S described above, we set the penalty matrix to S:= I_q. While this is a simplification of the correlation structure recommended by <cit.>, we think it offers a similar amount of flexibility.

Figure <ref> shows the estimated marginal mean and 95% credible bands, and Figure <ref> shows the subject-specific predicted curves. The point estimates shown in Figure <ref> are similar to that shown in Figure <ref>. However, the confidence intervals in Figure <ref> appear more reasonable. The subject-specific predicted curves shown in <ref> are not shrunk towards the mean as much as in Figure <ref>.





For the ℓ_2 Penalized model, in place of the correlation structure implied by the penalty matrix S described above, we augmented each Z_i matrix on the left with the columns [1, x_i], where x_i is an n_i × 1 vector of measurement times for subject i. We then replaced Z_i b_i with [1, x_i, Z_i] (u_i^T, b_i^T)^T, and assumed (u_i^T, b_i^T)^T ∼ N(0, Σ_i) where 

    Σ_i = [          Σ';             σ_b^2 I_q_i ]

and Σ' is a common 2 × 2 unstructured positive definite matrix. To model the within-subject correlations, we used a continuous autoregressive process of order 1. In particular, Cor(y_i(x_ij), y_i(x_ij')) = ζ^|x_ij - x_ij'| for a common parameter ζ>0.

Figure <ref> shows the estimated marginal mean and 95% credible bands, and Figure <ref> shows the subject-specific predicted curves. The estimates shown in Figure <ref> are similar to that shown in Figure <ref>. While estimates of the difference between low and high vigilance subjects differs between this model and the ℓ_2 penalized model in Section <ref>, the more notable difference is in the subject-specific predicted curves. As seen in Figure <ref>, the predicted subject-specific curves are not shrunk towards the mean as much as in Figure <ref>.





Table <ref> shows the mean squared error (MSE) and computing time for the ℓ_1 penalized and ℓ_2 penalized models. In Table <ref>, computing time for the ℓ_1 penalized model does not include cross-validation, because the parameters were hand-tuned (with only 17 subjects and a complex random effects structure, cross-validation did not find reasonable parameter values). As can be seen in Table <ref>, the alternative correlation structure led to smaller MSE for both the ℓ_1 and ℓ_2 penalized models, and less computing time for the ℓ_2 penalized model.





§ DISCUSSION AND POTENTIAL EXTENSIONS 


As demonstrated in this article, P-splines with an ℓ_1 penalty can be useful for analyzing repeated measures data. Compared to related work with ℓ_1 penalties, our model is ambitious in that we allow for multiple smoothing parameters and propose approximate inferential procedures that do not require Bayesian estimation. However, these are also the two aspects of our proposed approach that require additional future work. For P-splines with an ℓ_2 penalty, in most cases the knot placement is not critical so long as the number of knots is large enough <cit.>. We believe this also holds for P-splines with an ℓ_1 penalty, though further experimentation is needed to support this assumption. In practice, we recommend fitting models with a few different knot placements and widths to determine whether the model is sensitive to those choices for the data at hand.

Regarding estimation, our current approach of using ADMM and CV appears to work reasonably well for random intercepts, but is not yet reliable for random curves. In the future, we plan to develop more robust estimation techniques, particularly for smoothing parameters. As one possibility, we have done preliminary work to minimize quantities similar to GCV and AIC instead of the more computationally intensive CV, though these approaches do not seem as promising as their ℓ_2 counterparts. It may also be helpful to set the degrees of freedom prior to fitting the model. When possible, Bayesian estimation may be the most reliable way to currently fit these models. Bayesian estimation also opens the possibility of using other sparsity inducing priors, such as spike and slab models <cit.>.

Regarding inference, in future work it may be possible to use the δ quantity to bound difference between ℓ_1 and ℓ_2 penalized fits under certain assumptions on the data. It may also be helpful to investigate the use of post-selection inference methods to develop confidence bands for linear combinations of the active set, and to further investigate through simulations the performance of our proposed approximations of degrees of freedom. However, we note that our primary use of the degrees of freedom estimate d̂f̂ is to obtain the residual degrees of freedom d̂f̂_resid = n - d̂f̂, which we then use to estimate the variances σ̂^2_ϵ = r_2^2 / d̂f̂_resid. Therefore, when n ≫d̂f̂, σ̂^2_ϵ is not very sensitive to d̂f̂, in which case it is not critical for our purposes to obtain an exact estimate of degrees of freedom.

As for P-splines with an ℓ_2 penalty, users must select both the order M of the B-splines and the order k+1 of the finite differences. These choices will depend on the scientific problem and analytical goals. Using k=1 (2^nd order differences) is likely an appropriate starting point for most applications, and larger k could be used to increase the amount of smoothness.

For P-splines with an ℓ_2 penalty, in most cases the knot placement is not critical so long as the number of knots is large enough <cit.>. We believe this also holds for P-splines with an ℓ_1 penalty, though further experimentation is needed to support this assumption. In practice, we recommend fitting models with a few different knot placements and widths to determine whether the model is sensitive to those choices for the data at hand.

Regarding the rate of convergence, from Observation <ref> and the work of <cit.>, we know that for equally spaced data and F = I_n, P-splines with an ℓ_1 penalty achieve the minimax rate of convergence for the class of weakly differentiable functions of bounded variation. When there are less knots than data points, we do not think it is possible to achieve the minimax rate of convergence. However, if the knots are selected well, it may be possible to achieve the same performance in practice.

It could also be useful to extend these results to generalized additive models to allow for non-normal responses, and to extend the approach of <cit.> to include random effects and multiple smoothing parameters.



§ SUPPLEMENTARY MATERIAL


We have implemented our method in the R package  available at <https://github.com/bdsegal/psplinesl1>. All code for the simulations and analyses in this paper are available at <https://github.com/bdsegal/code-for-psplinesl1-paper>.



§ ACKNOWLEGEMENTS


We thank Margaret Hicken for sharing the data from the stress study.





§ SIMULATED DEMONSTRATION WITH TWO SMOOTHS


In this appendix, we simulated data similar to that in Section <ref>, but with an additional varying-coefficient smooth. In particular, we simulated data for two groups with 50 subjects in each group and between 4 and 14 measurements per subject (900 total observations). The data for subject i at time t was generated as y_it = β_0 + f_1(x_it) + f_2(x_it) 1[subject  i  in Group 2] + b_i + ϵ_it where b_i ∼ N(0, σ_b^2) and ϵ_it∼ N(0, σ^2_ϵ) for σ_b^2 = 1 and σ_ϵ^2 = 0.01. The true group means f_1(x) and f_2(x) are shown in Figure <ref> and the simulated data are shown in Figure <ref>.





We fit a varying-coefficient model with J=2 smooths to the data. In particular, we used ADMM and 5-fold CV to minimize

    β_0 ∈ℝ, β_1 ∈ℝ^p-1, β_2 ∈ℝ^p, b∈ℝ^Nminimize   1/2y - β_0 1 - F_1 β_1 - F_2 β_2 - Z b_2^2 
        + λ_1 D^(2)β_1 _1 + λ_2 D^(2)β_2 _1 + τb^T b.

where F_1, F_2 ∈ℝ^n × p were formed with second order (first degree) B-splines and p=21 basis functions, F_2 = diag(u) F_1 where u_i = 1[subject  i  in Group 2], and Z_il = 1 if observation i belongs to subject l and zero otherwise. We also fit an equivalent model with an ℓ_2 penalty using the  package <cit.>, i.e. with (λ_j / 2) D^(2)β_j _2^2 in place of λ_j D^(2)β_j _1 in (<ref>), j = 1, 2.

The estimated curves are shown in Figure <ref> for the ℓ_1 penalized model and in Figure <ref> for the ℓ_2 penalized model. We used 5-fold CV to estimate the smoothing parameters λ_1, λ_2 and τ in the ℓ_1 penalized model, and LME updates to estimate σ^2_b and b in the final model. As seen in Figures <ref> and <ref>, the fits are similar, but the results with the ℓ_1 penalized model are slightly closer to the truth.





Table <ref> shows the degrees of freedom and variance estimates with the ℓ_1 penalized and ℓ_2 penalized models. As seen in Table <ref>, variance estimates from both the ℓ_1 and ℓ_2 penalized models are very near the true values.





§ DETAILS FOR Λ_J^MAX 


Letting r_j = y - β_0 1 - ∑_ℓ j F_ℓβ_ℓ  - Z b be the j^th partial residuals, we can write the terms in (<ref>) that involve β_j as (1/2) r_j - F_j β_j_2^2 + λ_j D_j β_j _1. Then taking the sub-differential of (<ref>) with respect to β_j, we have

    0   = -F_j^T(r_j - F_j β̂_j) + D_j^T λ_j s_j

for some s_j = (s_j, 1,…,s_j, p_j - k_j - 1)^T where

    s_j,ℓ∈{1}    if  (D β̂_j)_ℓ > 1 
    {-1}    if  (D β̂_j)_ℓ < 1 
    
    [-1, 1]     if  (D β̂_j)_ℓ = 0.

Solving (<ref>) for β̂_j, we have β̂_j = (F_j^T F_j)^-1 F_j^T r_j - D_j^T λ_j s_j. Multiplying through by D_j and noting that D_j D_j^T is full rank and thus invertible, we have

    (D_j D_j^T)^-1 D_j β̂_j = (D_j D_j^T)^-1 D_j (F_j^T F_j)^-1 F_j^T r_j - λ_j s_j.

Setting D_j β̂_j = 0 in (<ref>), we get that (D_j D_j^T)^-1 D_j (F_j^T F_j)^-1 F_j^T r_j = λ_j s_j where s_j, ℓ∈ [-1, 1] for all ℓ. This can only hold if λ_j = (D_j D_j^T)^-1 D_j (F_j^T F_j)^-1 F_j^T r_j_∞, which gives us λ_j^max.



§ CONTROLLING TOTAL VARIATION WITH THE ℓ_1 PENALTY 


Let f(x) = ∑_j=1^p β_j ϕ_j^M(x). Suppose the knots are equally spaced, and let h_M-k-1 = (M-k-1) / (t_j+M-k-1 - t_j) for all j and 0 ≤ k < M - 1. Then on the interval [t_M = x_min, t_p+1 = x_max], from <cit.> we have

    d^k+1/dx^k+1f(x)    = h_M-1⋯ h_M-k-1∑_j=k+2^p∇^k+1β_j ϕ_j^M-k-1(x)

where ∇^k+1 is the (k+1)^th order backwards difference.

Let a^M_k+1 = max_j ∈{k + 2,… p }∫_x_min^x_maxϕ_j^M-k-1(x) dx. We note that a^M_k+1 is finite and positive for all 0 ≤ k < M-1. Then from (<ref>), we have

    1/h_M - 1⋯ h_M-k-1   ∫_x_min^x_max| d^k+1/dx^k+1 f(x) | dx 
       =  ∫_x_min^x_max| ∑_j=k+2^p ∇^k+1β_j ϕ_j^M-k-1(x) |  dx 
       = ∫_x_min^x_max| ∑_j=k+2^p (D^(k+1)β)_j-k-1ϕ_j^M-k-1(x) | dx 
       ≤∫_x_min^x_max∑_j=k+2^p | (D^(k+1)β)_j-k-1ϕ_j^M-k-1(x) | dx 
       = ∑_j=k+2^p ∫_x_min^x_max| (D^(k+1)β)_j-k-1ϕ_j^M-k-1(x) | dx 
       = ∑_j=k+2^p | (D^(k+1)β)_j-k-1| ∫_x_min^x_maxϕ_j^M-k-1(x) dx 
       ≤ a^M_k+1∑_j=k+2^p | (D^(k+1)β)_j-k-1| 
       = a^M_k+1 D^(k+1)β_1

where (<ref>) follows because ϕ_j^M-k-1(x) ≥ 0  ∀ x ∈ℝ.

Rewriting (<ref>), for 0 ≤ k < M-1 we have

    ∫_x_min^x_max| d^k+1/dx^k+1 f(x) | dx ≤ C_M,k+1D^(k+1)β_1

where C_M, k+1 = a^M_k+1 h_M-1⋯ h_M-k-1 is a constant. This shows that controlling the ℓ_1 norm of the (k+1)^th order finite differences in coefficients also controls the total variation of the k^th derivative of the function.


apalike


