{'fraction_non_alphanumeric': 0.05129685556937289, 'fraction_numerical': 0.01558891454965358, 'mean_word_length': 4.358638743455497, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of Relaxation Methods for Feature Selection in Mixed Effects Models

Linear Mixed-Effects (LME) models are a fundamental tool for modeling clustered data, including cohort studies, longitudinal data analysis, and meta-analysis. The design and analysis of variable selection methods for LMEs is considerably more difficult than for linear regression because LME models are nonlinear. The approach considered here is motivated by a recent method for sparse relaxed regularized regression (SR3) for variable selection in the context of linear regression. The theoretical underpinnings for the proposed extension to LMEs are developed, including consistency results, variational properties, implementability of optimization methods, and convergence results. In particular we provide convergence analyses for a basic implementation of SR3 for LME (called MSR3) and an accelerated hybrid algorithm (called MSR3-fast). Numerical results show the utility and speed of these algorithms on realistic simulated datasets. The numerical implementations are available in an open source python package pysr3.


Introduction
Linear mixed-effects (LME) models are used to analyze nested or combined data across a range of groups or clusters. These models use covariates to separate the total population variability (the fixed effects) from the group variability (the random effects). LMEs borrow strength across groups to estimate key statistics in cases where the data within groups may be sparse or highly variable, and play a fundamental role in population health sciences [25,21], meta-analysis [9,33], life sciences, and in many others domains [34].
This paper develops the theoretical bases for the algorithmic approach to variable selection within LME context presented in [27]. Although there are many successful algorithms and software for variable selection for linear regression, e.g. the Lasso method ( [28,12]) and related extensions, approaches to variable selection for LMEs are far less settled with few open source software tools available despite this being an active research area for over 20 years [7].
Variable selection for LMEs is significantly complicated by the underlying nonlinear structure associated with estimating variance parameters induced by the group structure. In the context of this paper, the nonlinearities come from the logs of the determinants of the within group variances as well as the regularizers used for variable selection. Current approaches make key design decisions including the choice of likelihood (e.g. marginal/restricted/h-likelihood), regularizer (e.g. 1 [5] or SCAD [15]), and information criteria ( [29,16]). The wide variety of these decision choices has likely contributed to small number of standardized software tools that allow for a comparison of different regularizers. The goal in [27] is to fill this gap by proposing a unified methodological framework that accommodates a variety of variable selection strategies based on a range of easily implementable regularizers. Here we provide a theoretical justification for the algorithmic approach to the solution of the marginalized maximum likelihood estimation problems presented in [27].
The approach is motivated by the sparse relaxed regularization regression (SR3) strategy developed in [32]. Both the approach of [32] and the MSR3/MSR3-fast extensions for LMEs described in [27] use auxiliary variables to decouple the smooth terms from the variable selection regularizer. The original and auxiliary variables are tethered by adding to the objective their norm squared difference. While the analysis of [32] relies on the least squares data-fitting term, here we develop the algorithmic design and analysis required for the nonlinear and nonconvex LME extension.
The paper proceeds as follows. The mathematical description of the LME model is given in Section 2 along with two results on the existence of solutions and a brief description of a naive PGD algorithm for their solution. The proposed relaxation strategy described in [27] is given in Section (3) where the the relaxation depends on a decoupling parameter η and log-barrier smoothing parameter µ. In particular, we introduce the optimal value function uη,µ used to exploit the decoupling of the likelihood function from the sparsity regularizer. Here uη,µ is obtained by partially minimizing the decoupled variables in the likelihood while keeping those in the regularizer fixed. Section 4 introduces the MSR3 algorithm as the PGD algorithm applied to the sum of uη,µ and the sparsity regularizer. A brief discussion the basic assumptions typically required for establishing the viability of the PGD algorithm for this formulation is given. Section 5 is the theoretical core of the paper. In this section w show that the optimal value function function uη,µ satisfies the properties necessary for the application of the PGD algorithm. In particular, we establish the Lipschitz continuity of ∇uη,µ (Lemma 14). The convergence results for MSR3 are presented in Section 6 for fixed values of η and µ. In Section 7 we address the key issues surrounding the initialization of the coupling and smoothing parameters η and µ when only approximate values for uη,µ and ∇uη,µ are known. Here we appeal to both variable metric ideas as well as properties of the interior point algorithm. In Section 8 we give a briefly synopsis of some of these results obtained in [27]. These results indicate that the use of the optimal value function can dramatically improve both the efficiency and the performance of the numerical solution procedure in terms of computational speed and accuracy in variable selection.

Models and Notation
Consider m groups of observations indexed by i, with sizes ni, so that the total number of observations is n := n1 + n2 + · · · + nm. Each group is paired with a design matrix of fixed features Xi ∈ R n i ×p and a matrix of random features Zi ∈ R n i ×q along with vectors of outcomes Yi ∈ R n i . Set X : [23,24], we define a Linear Mixed-Effects (LME) model as where β ∈ R p is a vector of fixed (mean) covariates, ui ∈ R q are unobservable random effects assumed to be distributed normally with zero mean and the unknown covariance matrix Γ, with S k + and S k ++ denoting the sets of real symmetric k × k positive semidefinite and positive definite matrices, respectively. We assume that the observation error covariance matrices Λi are known and that the random effects covariance matrix is an unknown diagonal matrix, Γ(γ) := Diag (γ) , γ ∈ R s + , where, for any vector γ ∈ R q , Diag (γ) is the diagonal matrix with diagonal γ.
Define ωi := Ziui + εi to be the unknown cluster-specific error vectors. Then model (1) can also be viewed as a correlated noise model with Yi = Xiβ + ωi, ωi ∼ N (0, Ωi(Γ)), Ωi(Γ) = ZiΓZ T i + Λi. This yields the marginalized negative log-likelihood function of a linear mixed-effects model ( [23]): Maximum likelihood estimates for β and Γ solve the problem and when Γ = Diag (γ) the problem becomes min β∈R p ,γ∈R q + L(β, γ) := LML(β, Diag (γ)). (4) In this setting, an entry γj takes the value 0 when the corresponding coordinates of all random effects uij are identically 0 for all i, or equivalently, the randomness in Yi is completely explained by i. The existence of solutions to (4) and, more generally, (3) follows from the techniques developed in [33].
Theorem 1 (Existence of a Minimizer). Let the assumptions in the statement of problem (3) hold. Then optimal solutions to (3) exist.
A standalone proof of Theorem 1 that allows us to give a simple extension to the variable selection case is given in Appendix A.
We approach feature selection for the model (1)-(4) by adding a regularizer to the objective yielding an optimization problem of the form min x L(β, γ) + R(β, γ) + δC(β, γ), where C := R p × R q + , R : R P × R q + → R+ := R+ ∪{+∞} is a lower semi-continuous (lsc) regularization term, and δC is the convex indicator function In practice, it is often advisable to include a constraint of the form γ ≤ γmax for γ ∈ R q ++ chosen sufficiently large since an excessively large variance usually indicates that the model is poorly posed and needs review. Such a constraint is also numerically expedient since it prevents γ from diverging. We return to this issue when our algorithm is specified. We have the following extension to Theorem 2 which tells us that solutions to (5) exist whenever R is level compact. The proof appears in Appendix A.
Theorem 2. Let the assumptions in the statement of problem (4) hold. Suppose R : . Then L + R is level compact and solutions to the following optimization problem exist: Algorithm 1: Proximal Gradient Descent (PGD) for LMEs Since L is smooth on its domain, a standard approach to solving (5) is the proximal gradient descent (PGD) algorithm described in Algorithm 1, where, for x = (β, γ) and w = (β,γ), Here the parameter L is intended to be a global Lipschitz constant for ∇ L over its domain. Unfortunately, ∇ L is not globally Lipschitz on its domain. Nonetheless, it is possible to obtain convergence results with the inclusion of a line search or trust region strategy [6].
In this paper a different approach is explored that uses global variational information on L rather than the local linearizations for L which form the basis of the PGD algorithm and its variants.

Relaxation of the mixed-effects variable selection model
Our strategy for obtaining approximate solutions to the mixed-effects variable selection problem (5) is motivated by the sparse relaxed regularization regression (SR3) strategy developed in [32]. That is, we introduce auxiliary variables to decouple two competing goals -variable selection and data fitting. In addition, we add a barrier term to relax the constraint 0 ≤ γ. The decoupling uses the coupling function κη : R p × R q → R given by (7) κη(y, z) : where η ∈ R+, while the constraint 0 ≤ γ is relaxed using the perspective of the negative log, i.e ϕ : R q × R → R ∪{∞} given by The mapping ϕ is known to be a closed proper convex function and, for µ > 0, it is essentially equivalent to the well-known log-barrier function. For more information on the perspective mapping, its calculus, and perspective duality, we refer the reader to [2,3]. We call η the coupling parameter and µ the log-barrier parameter and write φµ(·) := ϕ(·, µ). The relaxed problem employs auxiliary variables (β,γ) and relaxation parameters 0 ≤ η and 0 ≤ µ to obtain the problem We rewrite (8) so as to separate the smooth and nonsmooth components to obtain min (β,γ),(β,γ) where Lη,µ((β, γ), (β,γ)) := L(β, γ) + φµ(γ) + κη(β −β, γ −γ).
This is the formulation of the mixed-effects variable selection problem we study. Our focus is on the optimal value function uη,µ which captures global variational information about the function L over its domain. We show that uη,µ has a locally Lipschitz continuous gradient and that the evaluation of uη,µ and ∇uη,µ is accomplished by optimizing a well conditioned strongly convex function. This allows us to apply the PGD algorithm to the function uη,µ rather than the function L. Our numerical studies show that the global information captured by uη,µ significantly improves both the accuracy of the solution obtained and the overall numerical efficiency of the algorithm.
We follow the analysis of the PGD algorithm given in [4,Chapter 10] as it applies to the objective (13) Φη,µ(β,γ) := uη,µ(β,γ) Since uη,µ is nonconvex, one typically applies a line search method to select stepsize. However, this is often not required in practice. For this reason we state the algorithm with and without a line search.
Algorithm 3: Proximal Gradient Descent fo Φη,µ with Backtracking In Algorithm 4, the parameter L is assumed to be a global Lipschitz constant for ∇uη,µ. In Section 6, we show that the existence of L is not needed. In both algorithms we introduce the requirement that γ k ≤ γmax. While it is possible to include an explicit constraint of this form in the optimal variable selection problem (5), we do not do so since we assume that γmax is chosen so large that, from a practical perspective, the violation of this constraint indicates that the model is poorly posed and the algorithm needs to be terminated. We base our analysis of the convergence properties of Algorithms 1 and 2 on [4, Theorem 10.15] which makes use of the following three basic assumptions: Basic Assumptions for the PGD Algorithm (A) R : R p × R q → R is a closed proper convex function.
(C) Problem (11) has an optimal solution with optimal value ΦOPT.
We assume that (A) holds. This is not an overly restrictive assumption since it is satisfied by most of the standard variable selection regularizers. We show that (C) holds when R satisfies an additional coercivity hypothesis (Theorem 5). On the other hand, establishing that (B) holds in a concrete setting such as ours can be quite difficult. In particular, just as with L, uη,µ may fail to be globally Lipschitz. Validating Assumption (B) as well as developing a technique for circumventing the need for a global Lipschitz constant for ∇uη,µ consumes the majority of the theoretical development.

The Smoothness of uη,µ
We investigate the relationship between the problems (5) and (9), the existence of solutions to (9), and the properties of the function uη,µ and its derivative.
Proof. The formula for ∇ 2 L is given in Appendix B (see (45)). By [1, This establishes that the matrix in (15) is positive semidefinite. Since φµ is convex, the mapping is strongly convex for any choice ofη > νm, where For the remainder of the paper, we assume that η > νm =:η (17) so that the mapping (β, γ) → Lη,µ((β, γ), (β,γ)) is strongly convex with positive definite Hessian regardless of the choice ofγ ∈ R q . With this in mind, the function uη,µ defined by (12) resembles a Moreau envelope. However, this is misleading since, in particular, we are not even assured of the existence of solutions to the optimization problem defining uη,µ.

Existence and consistency.
To establish the existence of solutions to the relaxed optimization problems (9) and the problems defining the parametrized family uη,µ in (12), we assume that R is 1-coercive.
Proof. Let v * be the optimal value in (9) and let {( (41) and (42), it must be the case that If v * = −∞, then (18) tells us that This in turn implies that µ > 0, φµ(γ k ) → −∞ and γ k → ∞. Since R is 1-coercive and γ k → ∞, we can assume with no loss in generality that there is anᾱ > 0 such that which is a contradiction. Hence we may assume with no loss in generality that γ k → +∞. If µ = 0, then, by (18), ρ > n+1 2 ln(α) + R(β k ,γ k ) ↑ +∞, a contradiction, and so we can assume that µ > 0 and R is 1-coercive. Using (18) we may proceed as in (20) to find that again a contradiction, so the sequence {γ k } is bounded. Therefore, the first inequality in (18) tells us that the entire sequence {((β k , γ k ), (β k ,γ k ))} is necessarily bounded. Consequently, a limit point of the sequence {((β k , γ k ), (β k ,γ k ))} exists and, since R is lsc, any such limit point is a solution to (9).
Next we fix µ ≥ 0 and show that as η ↑ ∞ the solutions to (9) converge to solutions of In particular, for µ = 0, they converge to solutions of (5).
Proof. With no loss in generality η k >η for all k. Set where x = (β, γ) and w = (β,γ) with κη defined in (7). Set x k = (β k , γ k ),x = (β,γ), w k = (β k ,γ k ) andŵ = (β,γ). By Lemma 4 and Theorem 2 with R = φµ +R, there is an optimal solution xµ to (22) yielding an optimal value of vµ for which a k (x k , w k ) ≤ a k (xµ, xµ) = vµ for all k ∈ N. Hence, the sequence {a k (x k , w k )} is upper bounded by vµ. Since By adding these inequalities together we find that We now pair Theorem 6 with a consistency result for the barrier parameter µ.
Theorem 8 (Continuity of uη,µ and SSη,µ). Let the assumptions of Theorem 5 hold. For every (µ, η) ∈ R+ × R++, the function uη,µ defined in (12) is well-defined and continuous on R p × R q . In addition, the solution mapping SSη,µ is well-defined, single-valued and continuous on R p × R q .
Proof. Since η >η = νm, Lemma 3 tells us that the objective in (12) is strongly convex, and so (12) has a unique solution. Consequently, SSη,µ is well-defined and single-valued on R p × R q . This implies that uη,µ is also well defined on R p × R q since The result follows once it is shown that SSη,µ is continuous.
Plugging this into (24) and simplifying gives Therefore the sequence {(β k ,γ k )} must be bounded.
Our approach to establishing the Lipschitz continuity of SSη,µ is to first show that SSη,µ is differentiable and then obtain a bound on its Jacobian. As usual, diffentiability follows by applying the implicit function theorem to Gη,µ.
The formula for the inverse follows from (31) which tells us that (32) H Alternatively, one can apply the formulas in [20].
where Φ OPT η,µ := inf Φη,µ. Proof. As observed prior to the statement of the theorem, both Algorithm 4 and 3 behave as if they were applied to the the functionsûη,µ and R defined in (38). It was also shown that the functionsûη,µ and R satisfy the Assumptions (A)-(C) required by [4,Theorem 10.15]. Hence, the consequences of [4,Theorem 10.15] hold. By translating the notion of [4,Theorem 10.15] to that of this paper, we obtain the result.

A Hybrid Algorithms for Feature Selection in Mixed Effects Models
In the previous section we established the convergence properties of the PGD algorithm applied to the function Φη,µ for fixed values of η and µ. In subsection 5.2, two consistency results are established for the relaxed problem (11). Theorem 6 shows that, for fixed µ ≥ 0, every limit point of solutions to (11) as η ↑ ∞ is a solution to (22), while Theorem 7 tells us that every limit point of solutions (22) as µ ↓ 0 is a solution to the variable selection problem (5). These results suggest a range of numerical approachs to obtaining approximate solutions to the target problem (5). The issue of foremost concern is the method for approximating solutions to (12) since the accuracy in this approximation determines the accuracy in both uη,µ and ∇uη,µ. To address this concern, we view the algorithm from an interior point perspective where every point on the central path is a solution to the optimization problem (12) defining uη,µ for the associated value of the homotopy parameter µ. An approximate solution is then considered acceptable if it is sufficiently close to the central path where proximity to the central path is measured in terms of the notion of the neighborhood of the central path , e.g. see [30]. Due to the convexity of the optimization problems (12), this is an efficient algorithm for approximating uη,µ to high accuracy.
The next issue we addressed is the method for initializing and adjusting the parameter η. This is particularly significant since the initial value of η must be chosen to assure the convexity of the problems in (12). Lemma 3 gives us guidance in this regard, but the necessary computations to obtain a lower bound on η can be arduous, and, in general, produce a wildly pessimistic lower bound. For this reason, we take a somewhat different approach by proposing a variable metric strategy for solving the optimization problems in (12). In this approach, we replace the Hessian matrix ∇ 2 Lη,µ(β, γ) in the Newton equation Gη,µ((β, γ, v), (β,γ)) + ∇Gη,µ((β, γ, v), (β,γ))[dv, dβ, dγ] = 0 by the positive semi-definite approximation which is motivated by the expression for ∇ 2 Lη,µ(β, γ) given in (45). That is, we simply With this modification, the subproblems we solve are strongly convex for all η > 0. Consequently, the problem of initializing η is less problematic. Our numerical experiments indicate that the performance of the algorithm is robust with respect to η. For this reason, we choose an initial value for η and then leave it fixed over all iterations. Our method for choosing η is described in [27, Section 4, Figure 5]. Briefly, we maximize the Baysian Information Criterion (BIC) over a grid of values for η. The resulting BIC response curve shows that the method is robust with respect to the choice of η and choosing η ∈ [1, 10] yields accurate solutions for our selected test problems. Once η is fixed the PGD algorithm can be applied to solve the problem (11) for decreasing values of µ.

Numerical results
A detailed numerical study and comparison Algorithms 1, 4, and 5 as well as other algorithms for variable selection in LME models is given in [27]. Here we give one illustration from [27].
Experimental Setup. In this experiment we take the number of fixed effects p and random effects q to be 20. We set β = γ = [ 1 2 , 2 2 , 3 2 , . . . , 10 2 , 0, 0, 0, . . . , 0], i.e. the first 10 covariates are increasingly important and the last 10 covariates are not. The data is generated as yi = Xiβ + Ziui + εi, εi ∼ N (0, 0.3 2 I) Xi ∼ N (0, I) p , Zi = Xi ui ∼ N (0, Diag (γ)) We generated 9 groups with the sizes of [10,15,4,8,3,5,18,9,6] to capture a variety of group sizes. To estimate the uncertainty bounds, each experiment is repeated 100 times. Table 1 compares the performance of algorithms 1, 4, and 5 for four different feature selection regularizers: L0 (the 0-norm), L1 (the 1-norm), ALASSO (adaptive LASSO [5,17,31,18,10,22]), and SCAD (smoothed clipped absolute deviation [11,8,13]). Figure 1 gives a more detailed picture of the statistical performance of the algorithms over the set of 100 test problems. The L1 and ALASSO regularizers are convex while the L0 and SCAD are not. Despite the non-convexity of the L0 and SCAD regularizers, they exhibit superior accuracy in identifying the correct features. There are closed form expressions for the prox operator for all of these regularizers [27]. The hybrid MSR3-fast Algorithm 5 is the clear winner in terms of efficiency in that it produces highly accurate solutions in a tiny fraction of the time it takes Algorithms 1 and 4. As expected, the vanilla PGD Algorithm 1 is the least accurate in identifying the correct features since it is only a first-order method while Algorithms 4 and 5 both use higher-order information as well as incorporating global variational information on the relaxed objectives Lη,µ. The whisker plots in Figure 1 show that although Algorithm 4 has a slight edge in accuracy, Algorithm 5 strongly dominates both Algorithms 1 and 4 in speed.  Proof for Theorem 1. The bound (42) tell us that D := (r, Ω(Γ)) r ∈ R n , Γ ∈ S q + and f (r, Ω(Γ)) ≤ ρ ⊂ Dα,ρ .