Doubly structured sparsity for grouped multivariate responses with application to functional outcome score modeling

1]Jared D. Huling*

2,3]Jennifer P. Lundine

4,5]Julie C. Leonard

Huling, J.D. et al

[1]Division of Biostatistics, University of Minnesota, Minnesota, U.S.A.

[2]Department of Speech and Hearing Science, The Ohio State University, Ohio, U.S.A.

[3]Division of Clinical Therapies and Inpatient Rehabilitation Program, Nationwide Children’s Hospital, Ohio, U.S.A.

[4]Division of Emergency Medicine, Department of Pediatrics, The Ohio State University College of Medicine, Ohio, U.S.A.

[5]Abigail Wexner Research Institute, Nationwide Children’s Hospital, Ohio, U.S.A.

*Jared D. Huling<EMAIL_ADDRESS>
[Summary]This work is motivated by the need to accurately model a vector of responses related to pediatric functional status using administrative health data from inpatient rehabilitation visits. The components of the responses have known and structured interrelationships. To make use of these relationships in modeling, we develop a two-pronged regularization approach to borrow information across the responses. The first component of our approach encourages joint selection of the effects of each variable across possibly overlapping groups related responses and the second component encourages shrinkage of effects towards each other for related responses. As the responses in our motivating study are not normally-distributed, our approach does not rely on an assumption of multivariate normality of the responses. We show that with an adaptive version of our penalty, our approach results in the same asymptotic distribution of estimates as if we had known in advance which variables were non-zero and which variables have the same effects across some outcomes. We demonstrate the performance of our method in extensive numerical studies and in an application in the prediction of functional status of pediatric patients using administrative health data in a population of children with neurological injury or illness at a large children's hospital.

§ INTRODUCTION

Each year in the United States more than 20,000 children are diagnosed with an acute neurologic injury or illness that result in debilitating physical and cognitive complications, and reduced quality-of-life [1, 2, 3, 4]. These events are costly.  As an example, it is estimated that one billion dollars are spent annually on management of pediatric traumatic brain injury (TBI)-associated hospitalizations [5]. Though neurologic illnesses and injuries as a category are a leading cause of morbidity in children, each individual diagnosis is uncommon. These low disease-specific incidences make the study of rehabilitation interventions for children with neurologic injuries or illnesses challenging. Further, since many of the material effects manifest in the long term and can change with later child development, it is important to be able to track outcomes over time.

The WeeFIM® is a validated scoring system for appraising functional ability in children and provides valuable information about a multitude of components of health in children with neurologic injury or illness. WeeFIM® is administered and scored by trained assessors and has been demonstrated to have high interrater reliability. WeeFIM® has been shown to be predictive of longitudinal functional recovery of children with neurological disorders and can be used for discharge planning, prediction of functional outcomes, and documentation of functional performance of children over time to assess recovery or decline.
Yet, because the WeeFIM® system requires trained assessors, scores for all children are not always available.

WeeFIM® training and subscription is expensive and time-intensive and thus there is significant interest in being able to assess broadly the functional ability of children across a health system without the need for explicit scoring by trained assessors.
One avenue for doing so is to relate administrative data sources to WeeFIM® scores to facilitate care management and identify individuals who may require additional care or rehabilitation. Our interest is thus in building predictive models to understand functional ability of children using administrative health data in scenarios where WeeFIM® is not used. WeeFIM® scores have particular characteristics and organization. In this work we seek to leverage these characteristics to facilitate the building of accurate and interpretable WeeFIM® risk models. We now describe these characteristics and how we aim to utilize them in this paper.

WeeFIM® is comprised of 18 component scores, each of which is measured on a 7 point Likert scale with 7 indicating complete independence and 1 indicating complete dependence on others to perform various tasks. The WeeFIM® component scores are categorized into three main domains: mobility, cognition, and self-care, and represent three distinct clinical concepts related to functional outcomes. We describe the components of each domain in Section <ref>.
In our dataset, the correlations of the WeeFIM® component scores align closely with these domains. We computed the distance correlations [6] of all pairs of the 18 WeeFIM® components to form an $18\times 18$ distance matrix and ran a hierarchical clustering algorithm on the resulting distance matrix. As can be seen in Figure <ref>, the scores in the cognition domain are immediately completely separated from the scores in mobility and self-care;  self-care and mobility separate from each other later in the hierarchical clustering.

Heatmap of the distance correlations of the WeeFIM® scores. The dendrogram from hierarchical clustering of the scores shows that they group naturally by the three major WeeFIM® domains: self care, mobility, and cognition; thus the pre-defined domains align well with the natural variation in the components in our motivating study.

In our motivating study, we aim to use information commonly present in administrative health data including International Classification of Diseases (ICD) -9 and -10 codes, Current Procedural Terminology (CPT) codes, pharmaceutical codes, data indicating durable medical equipment use, and basic demographic information to predict and evaluate functional outcomes in children admitted to an inpatient rehabilitation unit with a diagnosis of neurologic injury or illness at a large children's hospital. Our code information is high-dimensional, with well over 500 billing codes as potential predictors of a child's functional status.

To deal with high-dimensionality, one could fit lasso-penalized regression models [7] for each of the 18 scores separately, however this would ignore the inherent relationships and similarities between the components. First, as some WeeFIM® components are highly related (e.g. dressing lower and dressing upper, bowel management and bladder management, bathing and tub shower), it is plausible that some variables (e.g. billing codes) may have the same or similar effects across some of the WeeFIM® components. Similarly, some variables may not have the same effect size across some WeeFIM® components, but are likely to be jointly predictive for multiple components, especially for related components, or even for all components. A complication is that while linearity of effects of covariates on the components may be reasonable, an assumption of normality of the responses or error terms is not, as the WeeFIM® component scores are not normally-distributed.

There exists some work for building models for multivariate responses using high-dimensional predictors.
Rothman, et al. [8] introduced a lasso-based estimation procedure that incorporates variable selection and covariance estimation for multivariate normal responses and developed an iterative estimation procedure that makes use of the covariance structure among the responses by penalizing the inverse covariance matrix using the methods of Yuan and Lin [9].  However, this approach does not make use of similar predictor-response relationships across related responses and only borrows information across responses through their covariance structure. Sofer, et al. [10] extended these ideas to incorporate a wider variety of penalties such as the SCAD penalty [11] and the MC+ penalty[12] and developed new algorithms for estimation. Li, et al. [13] incorporated group structure both in responses and among predictors using an overlapping group lasso penalty [14], which allows for borrowing of strength across responses through joint selection of variables across related responses.

In our work, we borrow strength across responses in two ways through two different structured sparsity inducing penalties. In particular, we first aim to incorporate joint selection and removal of variables across all responses and joint selection and removal of variables across related responses (i.e. responses that are in the same functional domain). In our application, the group structure is the natural groups formed by the WeeFIM® components displayed in Figure <ref>.  Thus, we require use of a group lasso [15] with overlapping groups [14] similar to Li, et al. [13]. Second, we use a fused lasso [16] for shrinkage of the effects of a variable across related outcomes to borrow strength more explicitly in estimation by partially collapsing models across different responses into a single model for individual variables separately. The fused lasso penalty allows a variable's effects across related responses to be estimated to be exactly the same.
We prove that with adaptively chosen weights [17, 18] but for the overlapping group lasso as in Huling, et al. [19] and fused lasso as in Viallon, et al. [20], our doubly-structured sparsity inducing estimator has an oracle property. The theoretical results are general in that they allow for any arbitrary overlapping group structure for the group lasso penalty and any arbitrary fused lasso penalty. The oracle property we show demonstrates that estimation for the non-zero coefficients has the same asymptotic distribution as if we had known in advance both i) which coefficients are non-zero and ii) which coefficients are equal to each other and fit models with only the non-zero coefficients and forcing the truly equal coefficients to be equal to each other. We also show selection consistency for the non-zero coefficients and selection consistency for which coefficients are equal to each other. Neither our methodology nor our theory requires an assumption of normality or independence of the responses, as we work under a semiparametric multivariate linear model assumption, where we make no parametric distributional assumptions about the error terms. Further, our theory, methodology, and computational approach allows for general use of the overlapping group lasso in combination with the fused lasso beyond our application to multivariate response regression.

The remainder of our paper is organized as follows. In Section <ref> we introduce the key definition of our methodology. In Section <ref> we develop theory for our methodology using adaptive regularization under semiparametric multivariate linear models and in Section <ref> we investigate the operating characteristics of our method in small samples using simulation experiments. In Section <ref> we demonstrate the use of our methodology on our motivating application involving pediatric functional outcomes score modeling. Finally, we conclude with discussion.

§ DOUBLY STRUCTURED SPARSITY FOR MULTIVARIATE RESPONSES

§.§ The Overlapping Group + Fused Method

We assume that the observed multivariate response $\bfy_i = (y_{i1}, \dots, y_{iK})^\top$ follows the semiparametric linear model

\begin{align}
\bfy_i = {\betatrue}^\top\bfx_i + \boldsymbol\epsilon_i \text{ for } i=1,\dots,N, \label{lin_model}
\end{align}

where $\bfx_i = (x_{i1}, \dots, x_{ip})^\top$ is a length $p$ vector of predictors, $\betatrue$ is a $p\times K$ matrix of regression coefficients with $k$th column equal to $\betatrue_{\cdot,k} = (\beta^0_{1,k}, \dots, \beta^0_{p,k})^\top$, and $\bepsilon_i = (\epsilon_{i1}, \dots, \epsilon_{iK})^\top$ is a vector of errors with mean zero and finite, positive definite variance-covariance matrix $\bSigma$; no requirement that the response vector be continuous is required, as we do not assume a particular distribution for $\boldsymbol\epsilon_i$. Denote the $j$th variable's coefficients across the $K$ response as $\betatrue_{j,\cdot} = (\beta^0_{j,1}, \dots, \beta^0_{j,K})^\top$. Let $\bfX$ denote the $N\times p$ matrix of predictors with $i$th row as $\bfx_i^\top$ and $j$th column as $\bfx_{\cdot,j}$ and let $\bfY$ denote the $N\times K$ response matrix with $i$th row as $\bfy_i^\top$ and denote the $k$th column of $\bfY$ as $\bfY_{\cdot,k}$.

The typical least squares estimator is the solution of

\begin{align*}
\argmin_{\bbeta} (2N)^{-1}\tr\left[(\bfY - \bfX\bbeta)(\bfY - \bfX\bbeta)^\top\right],
\end{align*}

which is simply $(\bfX^\top\bfX)^{-1}\bfX^\top\bfY$, where for a matrix $\bfA$, $\tr(\bfA)$ indicates its trace.

In our setting, the $K$ outcomes are grouped into $M$ different, possibly overlapping groupings that describe inherent relationships between the outcomes. In our motivating data, the WeeFIM outcomes are grouped into three pre-defined domains that represent different categories of functional status. In other settings, these domains may be further subcategorized. In this manner, the groupings  generally correspond to hierarchically-organized categorizations of the outcomes. In this work, we make use of the natural grouping of the outcomes by adding regularizers that induce the coefficients of each different variable to be selected or removed jointly across an entire group of the responses.  The $j$th grouping $\calG_j$ consists of a set of $g_j$ groups $G_{j,1}, \dots, G_{j,g_j}$ such that the union of the groups includes all of the $K$ outcomes.
Generally, since the groupings will be hierarchically defined,  $g_j$ will be larger than $g_{j+1}$ and groups in $\calG_{j}$ can be expressed as unions of groups in $\calG_{j+1}$.

More formally, the groupings are written as $\calG_1 = \{G_{1,1}, \dots, G_{1,g_1}\}, \dots, \calG_M = \{G_{M,1}, \dots, G_{M,g_M}\}$, where $G_{m,g} \in \{1,\dots,K\} \equiv \calK$ and $\bigcup_{G\in\calG_m}G = \calK$ for $m=1,\dots, M$.

To encourage a variable to be selected jointly across all responses, we define the trivial group $\calG_0 = \{\calK\} $, thus $G_{0,1} = \calK  = \{1, \dots, K \}$. To allow variables to be selected or removed individually, we further include the group $\calG_{M+1} = \{\{1\}, \dots, \{K\}\}$, thus $G_{M+1,k} = \{k\}$. The group structure for a toy example with $K=8$ responses and $M=1$ groupings is displayed in Figure <ref>.
We note that in this toy example, the groups are overlapping (e.g. $G_{1,1} \subset G_{0,1}$) because the groups are hierarchically structured. In our motivating application, the group structure we use involves such overlapping groups.

In our setting the outcome groupings are pre-defined, but in general when such a pre-defined set of groups are not available a priori,
the $M$ groupings can be formed by iteratively refined, hierarchical groupings of the responses, as illustrated in Figure <ref>; for example, the $M$ groupings could be formed by clustering the responses at $M$ different levels in a hierarchical clustering of the responses.

For a particular group $G\in\calK$ with size $|G|$, define the length $|G|$ subvector of $\bbeta_{j,\cdot}$ limited to the outcomes in group $G$ as $\bbeta_{j,G}$.
We propose to add to the above objective function two regularizers that induce sparsity in a structured manner. Let $\betahat$ be the solution of the following problem

\begin{align}
\argmin_{\bbeta} (2N)^{-1}\tr\left[(\bfY - \bfX\bbeta)(\bfY - \bfX\bbeta)^\top\right] + \lambda_1 P_1(\bbeta) + \lambda_2 P_2(\bbeta) \label{overlapping_linear_model_objective}
\end{align}

$$P_1(\bbeta) = \sum_{j=1}^p\sum_{m=0}^{M+1}\sum_{G \in \mathcal{G}_m} \lambda_{1,j,G} ||\bbeta_{j,G}||_2$$

$$P_2(\bbeta) = \sum_{j=1}^p\sum_{G\in \calG_M}\sum_{l, o \in G: l \neq o } \lambda_{2,j,l,o}|{\beta}_{j,l} - {\beta}_{j,o}|,$$

where the penalty $\lambda_{1,j,G} ||\bbeta_{j,G}||_2$ is an overlapping group lasso penalty encourages joint selection or removal of the effects of variable $j$ across the responses in response group $G$. All possible zero patterns of coefficients can be represented as unions of groups due to the results of Jenatton, et al [14]; hence nonzero patterns of coefficients can be thought of as complements of unions of groups. Thus, the grouping $\calG_0$ encourages joint selection and removal of effects for a variable across all outcomes simultaneously and grouping $\calG_{M+1}$ allows for individual selection or removal for each outcome separately, while the groupings $\calG_{1}, \dots, \calG_M$ allow for selection or removal of effects across related groups.
The second penalty $\lambda_{2,j,l,o}|{\beta}_{j,l} - {\beta}_{j,o}|$ is a fused lasso penalty that encourages the effect of variable $j$ on response $l$ and on response $o$ to be more similar. Thus, $P_1(\bbeta)$ and $P_2(\bbeta)$, respectively, help borrow strength across the $K$ responses by 1) leveraging the natural groupings of the responses and taking advantage of joint significance of predictors across related responses and 2) encouraging the effect estimates for a particular variable to be similar across the most related responses according to the response groupings.  These two regularizers incorporate structural knowledge about outcomes in two manners by utilizing two different types of structured sparsity inducing penalties. The overall structure of the regularizers are depicted in Figure <ref> in a toy example with $K=8$ responses. In our formulation, we only add a fused lasso penalty within the final grouping $\calG_M$, reflective of the notion that with an iteratively refined grouping, the final groups are likely to have more similar outcomes within each group. However, the fused lasso penalty can be added within any grouping and can even be set so that fused lasso terms are included for every possible pair of responses; our asymptotic results cover any arbitrary set of fused pairs. In our application, the multivariate outcomes are measured on the same scale, but our approach can be adapted to outcomes that have different scales by standardizing the outcomes prior to analysis. Further, as we list below, our adaptive penalization allows for the fused lasso penalty to adapt to observed differences in effects across outcomes.

To perform adaptive penalization, we take $\lambda_{1,j,G} = ||\betahat^{OLS}_{j,G}||_2^{-\gamma_1}$ and $\lambda_{2,j,l,m} = |\hat{\beta}_{j,l}^{OLS} - \hat{\beta}_{j,m}^{OLS}|^{-\gamma_2}$, where $\gamma_1, \gamma_2 > 0$ and $\betahat^{OLS} = (\bfX^\top\bfX)^{-1}\bfX^\top\bfY$ is the ordinary least squares estimate of $\bbeta$. We study this adaptive version of our estimator in Section <ref>. For non-adaptive penalization, we take $\lambda_{1,j,G} = |G|^{1/2}$, as is used in Jenatton, et al. [14] and, Huling et al. [19] and $\lambda_{2,j,l,m}=1$, however for this choice selection consistency is not guaranteed without an irrepresentable condition [14].
In order to get reliable estimates of out-of-sample performance using cross validation when the group structure is derived in a data-driven manner, the procedure used to estimate the group structure should be applied in each cross validation fold instead of being fixed.

background rectangle=
[draw=blue!8,fill=blue!8,rounded corners=1.5ex]
[y=1cm, x=0.8cm, thick, font=,
vc/.style =
circle, draw, thick, fill=#1,
minimum width=12mm, opacity=0.7,
vg/.style args = #1/#2
minimum height=12mm,
minimum width=#1+/pgf/minimum height,
draw, rounded corners=/pgf/minimum height/2,

frames/.style args = #1/#2minimum height=#1,
minimum width=#2+/pgf/minimum height,
draw, thick, rounded corners=3mm, opacity=0.6, color=red, dashed,
framesbig/.style args = #1/#2minimum height=#1,
minimum width=#2+/pgf/minimum height,
draw, thick, rounded corners=4mm, opacity=0.8, color=blue, dashed,

number line/.style=

pen colour=blue,
line width=1.0pt,
decoration=calligraphic brace,amplitude=8pt,
pen colour=red,
line width=1.0pt,
decoration=calligraphic brace,amplitude=8pt,
pen colour=red,
decoration=calligraphic brace,amplitude=8pt, mirror,

ıin 1,...,8

[draw,thick,color=Green,circle, inner sep=0.1, minimum size=1.5em,scale=0.85] (coefı) at (ı-1,0)$\textcolor{black}{\beta_{j,\i}}$;

(start_week) at (-0.4,0.95) ;
(end_week) at (7.4,0.95) ;
[brace_top] (start_week.north) – node [above=4pt, pos=0.5] $G_{\textcolor{blue}{0},1}$ (end_week.north);

[label=[black, align=left,scale=1.15]right: $\textcolor{blue}{\calG_0}$] (fl) at (7.5,1.5);
[label=[black, align=left,scale=1.15]right: $\textcolor{red}{\calG_M = \calG_1}$] (fl) at (7.5,0.925);
[label=[black, align=left,scale=1.15]right: $\textcolor{Green}{\calG_{M+1} = \calG_2}$] (fl) at (7.5,0);

(start_week) at (-0.4,0.35) ;
(end_week) at (2.4,0.35) ;
[brace_top_red] (start_week.north) – node [above=4pt, pos=0.5] $G_{\textcolor{red}{1},1}$ (end_week.north);

(start_week) at (2.6,0.35) ;
(end_week) at (4.4,0.35) ;
[brace_top_red] (start_week.north) – node [above=4pt, pos=0.5] $G_{\textcolor{red}{1},2}$ (end_week.north);

(start_week) at (4.6,0.35) ;
(end_week) at (7.4,0.35) ;
[brace_top_red] (start_week.north) – node [above=4pt, pos=0.5] $G_{\textcolor{red}{1},3}$ (end_week.north);

[label=[black, align=right,rotate=90,scale=0.85]left:Fused
Lasso] (fl) at (-0.9,-0.5);

[label=[black, align=right,rotate=90,scale=0.85]left:Lasso] (fl) at (-0.9,0.35);

[label=[black, align=right,rotate=90,scale=0.85]left:Group
Lasso] (fl) at (-0.9,1.25);

[->,thin,Green] (-0.7,0)  to [bend right=-55] (-0.275,0.2);

[->,thin,red] (-0.6,0.85)  to [bend right=-45] (1.05,0.375);
[->,thin,blue] (-0.6,0.75)  to [bend right=-45] (0.25,0.45);

[->,thin] (coef1.south)++(0,0cm)  to [bend right] (0.4,-0.65) node[below right, draw=none] ;
[->,thin] (coef2.south)++(-0.05,0cm)  to [bend left] (0.6,-0.65) node[below left, draw=none] ;

[scale=0.425] (fuse12) at (0.5,-0.8) $|\beta_{j,1}-\beta_{j,2}|$;

[->,thin] (coef2.south)++(0.05,0cm)  to [bend right] (1.4,-0.65) node[below right, draw=none] ;
[->,thin] (coef3.south)++(0,0cm)  to [bend left] (1.6,-0.65) node[below left, draw=none] ;

[scale=0.425] (fuse23) at (1.5,-0.8) $|\beta_{j,2}-\beta_{j,3}|$;

[->,thin] (coef1.south)++(-0.075,0cm)  to [bend right=45] (0.5,-1.15) node[below right, draw=none] ;
[->,thin] (coef3.south)++(0.075,0cm)  to [bend left=45] (1.5,-1.15) node[below left, draw=none] ;

[scale=0.425] (fuse13) at (1,-1.15) $|\beta_{j,1}-\beta_{j,3}|$;

[->,thin] (coef4.south)++(0,0cm)  to [bend right] (3.4,-0.65) node[below right, draw=none] ;
[->,thin] (coef5.south)++(0,0cm)  to [bend left] (3.6,-0.65) node[below left, draw=none] ;

[scale=0.425] (fuse45) at (3.5,-0.8) $|\beta_{j,4}-\beta_{j,5}|$;

[->,thin] (coef6.south)++(0,0cm)  to [bend right] (5.4,-0.65) node[below right, draw=none] ;
[->,thin] (coef7.south)++(-0.05,0cm)  to [bend left] (5.6,-0.65) node[below left, draw=none] ;

[scale=0.425] (fuse67) at (5.5,-0.8) $|\beta_{j,6}-\beta_{j,7}|$;

[->,thin] (coef7.south)++(0.05,0cm)  to [bend right] (6.4,-0.65) node[below right, draw=none] ;
[->,thin] (coef8.south)++(0,0cm)  to [bend left] (6.6,-0.65) node[below left, draw=none] ;

[scale=0.425] (fuse78) at (6.5,-0.8) $|\beta_{j,7}-\beta_{j,8}|$;

[->,thin] (coef6.south)++(-0.075,0cm)  to [bend right=45] (5.5,-1.15) node[below right, draw=none] ;
[->,thin] (coef8.south)++(0.075,0cm)  to [bend left=45] (6.5,-1.15) node[below left, draw=none] ;

[scale=0.425] (fuse68) at (6.0,-1.15) $|\beta_{j,6}-\beta_{j,8}|$;

let 2 = ($(coef1.center)-(coef8.center)$),
2 = veclen(2,2) in
(coef1) – node[framesbig=8.5mm/2]  (coef8);
let 2 = ($(coef6.center)-(coef8.center)$),
2 = veclen(2,2) in
(coef6) – node[frames=7mm/2]  (coef8);
let 2 = ($(coef1.center)-(coef3.center)$),
2 = veclen(2,2) in
(coef1) – node[frames=7mm/2]  (coef3);
let 2 = ($(coef4.center)-(coef5.center)$),
2 = veclen(2,2) in
(coef4) – node[frames=7mm/2]  (coef5);

An illustration of the structure of the two components of the penalization methods used to induce group-wise selection of covariate effects and fuse them to be more similar. Blue indicates joint selection across all responses via the group $\calG_0$, red indicates joint selection by response group $\calG_M = \calG_1$, green indicates  individual coefficient selection achieved via the group $\calG_{M+1}$,  and black indicates shrinkage towards a common effect within response group $G_{1,k}$.

§.§ Implementation Details

In our implementation, we use a generic ADMM algorithm described in the next section for computation and feed into it the response vector $\wtbfY$ and design matrix $\wtbfX$, where the latter is treated as a sparse matrix object that only stores the values and locations of the non-zero entries of $\wtbfX$, which both dramatically saves space, memory used, and reduces computation time. The sparse matrix object is of the type provided in the  linear algebra library [21] with interface to an  package through  [22], allowing for highly efficient sparse-matrix manipulations. Our  package implementation of our method is  available in the open source repository <https://github.com/jaredhuling/groupFusedMulti>, which uses an interface similar to the interface used in .

For ease of tuning the tuning parameters, we utlize the following re-parameterization of $\lambda_1$ and $\lambda_2$. Instead of using the penalty $\lambda_1 P_1(\bbeta) + \lambda_2 P_2(\bbeta)$, we use $\lambda (1-\alpha) P_1(\bbeta) + \lambda\alpha P_2(\bbeta)$, where $\lambda \geq 0$ and $\alpha\in [0,1]$ so that $\alpha$ controls the proportion of the total penalty that the fused lasso comprises.

§.§ Computation via a Multi-block ADMM Algorithm

We utilize an alternating direction method of multipliers (ADMM) [23, 24, 25] algorithm for optimization. The ADMM algorithm works by decomposing an objective function and solving the decomposed subproblems iteratively, where each subproblem has a simple and computationally tractable form. ADMM solves problems of the form

\begin{align*}
& \mbox{minimize } f(\bbeta) + P(\bgamma) \mbox{ subject to } \bfA\bbeta + \bfB\bgamma = \bfc
\end{align*}

where $\bfA$ and $\bfB$ are constraint matrices, both $f$ and $P$ are convex functions, $\bbeta$ is the parameter vector of interest, and typically $f$ represents some loss function and $P$ represents a penalty. In the simplest case the constraint is of the form $\bbeta = \bgamma$ and the purpose of the constraint is to find a new variable for which the penalty is equivalent, but separable across loss and penalty in the terms of the new variable. To optimize our objective, we require the following multi-block version of ADMM

\begin{align*}
& \mbox{minimize } f(\bbeta) + P_1(\bgamma)  + P_2({\boldsymbol \eta}) \mbox{ subject to } \bfA\bbeta + \bfB\bgamma + \bfC{\boldsymbol \eta} = \bfc,
\end{align*}

where $\bfC$ is an additional constraint matrix that allows the second penalty to have a decomposed form with loss and the two penalties all separable from each other.
To solve the above problem, the augmented Lagrangian is formed as:

\begin{align*}
L_\rho(\bbeta, \bgamma, {\boldsymbol \eta}, \bnu) = {} &  f(\bbeta) + P_1(\bgamma)  + P_2({\boldsymbol \eta}) + \bnu^\top (\bfA\bbeta + \bfB\bgamma +  \bfC{\boldsymbol \eta} - \bfc) \\
& + (\rho/2)||\bfA\bbeta + \bfB\bgamma +  \bfC{\boldsymbol \eta}- \bfc||^2_2,
\end{align*}

where $\rho$ is any strictly positive number.
The multi-block ADMM algorithm iterates by alternatingly minimizing with respect to $\bbeta$, $\bgamma$, and ${\boldsymbol \eta}$ and following these minimizations, updating the Lagrangian parameter $\bnu$:

\begin{align}
\bbeta^{(t + 1)} = {} &  \argmin_{\bbeta}L_\rho(\bbeta, \bgamma^{(t)}, {\boldsymbol \eta}^{(t)}, \bnu^{(t)}) \label{eqn:betamin}\\
\bgamma^{(t + 1)} = {} & \argmin_{\bgamma}L_\rho(\bbeta^{(t + 1)}, \bgamma, {\boldsymbol \eta}^{(t)}, \bnu^{(t)})\label{eqn:gammamin} \\
{\boldsymbol \eta}^{(t+1)} = {} & \argmin_{\boldsymbol \eta}L_\rho(\bbeta^{(t + 1)}, \bgamma^{(t + 1)}, {\boldsymbol \eta}, \bnu^{(t)})\label{eqn:etamin} \\
\bnu^{(t + 1)} = {} &  \bnu^{(t)} + \rho (\bfA\bbeta^{(t + 1)}  + \bfB\bgamma^{(t + 1)}  + \bfC {\boldsymbol \eta}^{(t+1)} - \bfc) \nonumber
\end{align}

where $t$ indexes the iteration number. The standard ADMM has been shown to converge for any $\rho > 0$ and the multi-block ADMM algorithm above has been shown to converge under certain conditions on $\bfA$, $\bfB$, and $\bfC$. These conditions, from Chen, et al [26], are met if either $\bfA^\top \bfB = \bfzero$,  $\bfB^\top \bfC = \bfzero$, or $\bfA^\top \bfC = \bfzero$.

The following describes the multi-block ADMM algorithm applied to the overlapping group lasso with fused lasso problem. Let $m = \sum_{G \in \mathcal{G}}|G|$, $g = |\mathcal{G}|$, and suppose $\mathcal{G}=\{G_1,\cdots,G_g\}$, let $\bfF = (\bfF_1, \dots, \bfF_g)$ be a matrix of dimension $m \times Kp$ where $\bfF_l$ is a $|G_l|\times Kp$ matrix with $(i,j)$th entry equal to 1 if $j$ is the $i^{th}$ element of group $G_l$, and 0 otherwise, $\forall j=1,\dots,g$. Then $\bfF\bbeta$ is a vector of length $m$ comprised components of $\bbeta$ where each element of $\bbeta$ appears in $\bfF\bbeta$ the total number of times it appears in any group. For example, if p=1, K=3, $ \bbeta = (\beta_1, \beta_2, \beta_3)^\top $ and $\mathcal{G} = \{ \{ 1, 2\}, \{2, 3 \} \}$, then

\bfF =
\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{pmatrix}
\mbox{ and }
\bfF\bbeta =
\begin{pmatrix}
\beta_1 \\
\beta_2 \\
\beta_2 \\
\beta_3
\end{pmatrix}.

In this example, the penalty $P_1(\bgamma)$ is
$ which is a standard group lasso penalty on the variable $\bgamma$.
In general, the overlapping group lasso penalty $P_1(\bgamma)$ can be written as
where $\bgamma=(\bgamma_{1},\dots,\bgamma_{g})$ and $\bgamma_{l}$ is a $|G_l|$-dimensional vector.
Note that for the nonoverlapping group lasso, $\bfA = \bfI_{Kp}$, where $\bfI_{Kp}$ is the identity matrix of dimension $Kp\times Kp$.
To accommodate the fused lasso penalty, we allow for a more general penalty, the generalized fused lasso penalty. Denote the matrix $\bfD$ with $e$ rows and $Kp$ columns where each row represents a pair of effects which have a fused lasso penalty applied. Thus $e$ is the total number of fused lasso terms/pairs. For example, if the pair of effects $\beta_{j,\ell}$ and $\beta_{j,m}$ are have a fused lasso applied, then the row of $\bfD$ corresponding to this pair will have $\ell$th position taking the value 1 and the $m$th position taking the value $-1$. In the toy example in Figure <ref>, if $\beta_1$ has a fused lasso penalty with $\beta_2$ and $\beta_2$ has a fused lasso penalty with $\beta_3$, then

\bfD =
\begin{pmatrix}
1 & -1 & \hphantom{-}0 \\
0 & \hphantom{-}1 & -1
\end{pmatrix},
\mbox{ and }
\bfD\bbeta =
\begin{pmatrix}
\beta_1 - \beta_2\\
\beta_2 - \beta_3
\end{pmatrix}.

In this example, the $P_2(\boldsymbol\eta)=\lambda_2(\lambda_{2,1}|\eta_1| + \lambda_{2,2}|\eta_2|)$, which is a standard lasso penalty on the variable $\boldsymbol \eta$.
In general, the fused lasso penalty $P_2(\boldsymbol\eta)$ can be written as
where $\boldsymbol\eta$ is an $e$-dimensional vector.
The ADMM algorithm for the overlapping group lasso plus the fused lasso is constructed by taking
$$\bfA = \begin{pmatrix}
\bfF \\
\bfD
\end{pmatrix}, \bfB = \begin{pmatrix}
-\bfI_m \\
\boldsymbol 0
\end{pmatrix}, \bfC = \begin{pmatrix}
\boldsymbol 0 \\
\end{pmatrix},
\text{ and } \bfc = \bfzero,$$
which meets the condition that $\bfB^\top\bfC=\bfzero$, meaning a multi-block ADMM algorithm using this setup is valid and convergent. When $f(\bbeta) = \frac{1}{2}|| \wtbfY - \wtbfX\bbeta||^2_2$, step (<ref>) for the overlapping group lasso is simply the solution of $(\wtbfX^\top\wtbfX + \rho (\bfF^\top \bfF+ \bfD^\top \bfD) )\bbeta = \wtbfX^\top\wtbfY + (\bfF^\top, \bfD^\top) \bnu^{(t)} + \rho (\bfF^\top \bgamma^{(t)} + \bfD^\top{\boldsymbol \eta}^{(t)})$.  When $f(\bbeta)$ is the negative log-likelihood, step (<ref>) can be carried out by Newton-Raphson or other standard optimization techniques. As step (<ref>) is group-separable, it can be minimized by minimizing with respect to each group $\bgamma_{l}$ independently. This is achieved by the block soft-thresholding operator $S_{\lambda_1 \lambda_{G_l} / \rho}((\bfF\bbeta^{(t+1)})_{l} - \bnu_{l}^{(t)}/ \rho)$, where $S_\lambda(\bfu) = \bfu\left(  1 - \lambda / ||\bfu||_2  \right)_+$ and $(\bfF\bbeta^{(t+1)})_{l}$ and $\bnu_{l}^{(t)}$ are defined in the same way as $\bgamma_l$. Since (<ref>) is separable and is equivalent to a type of lasso penalization problem, it can be similarly minimized simply via soft thresholding: $S_{\lambda_2 \lambda_{j} / \rho}((\bfD\bbeta^{(t+1)})_{j} - \bnu_{j}^{(t)}/ \rho)$, where with some abuse of notation $(\bfD\bbeta^{(t+1)})_{j}$ is the $j$th element of $(\bfD\bbeta^{(t+1)})$,  $\bnu_{j}^{(t)}$ is the $j$th element of $\bnu^{(t)}$, and $S_\lambda(u) = \text{sign}(u)\left(  |u| - \lambda  \right)_+$, which is the equivalent of the block soft thresholding operator applied to a scalar.
Our convergence criterion is the same as suggested in Section 3.3.1 of Boyd, et al [25] with $\epsilon^{\mbox{abs}} = \epsilon^{\mbox{rel}} = 10^{-5}$. The parameter $\rho$ for the ADMM algorithm used is the adaptive scheme described in Section 3.4.1 of Boyd, et al [25].

§.§ Connections with existing literature

The methods of Rothman, et al [8] and Li, et al [13] both also address scenarios where high dimensional covariates are used to predict a set of multivariate outcomes in a penalized regression framework. These two methods in addition to ours borrow strength in information across outcomes through penalization techniques, albeit in different ways.

If the multivariate outcomes are only correlated due to correlations between the outcomes themselves, then the multivariate regression with covariance estimation (MRCE) approach  [8] may be most appropriate as it directly incorporates the covariance of the responses via a penalized log-likelihood; here, the lasso penalty induces the estimated coefficients to be a function of the covariance of the responses, unlike the unpenalized maximum likelihood estimator. However, when many covariates  are expected to have effects of similar magnitude across the multivariate outcomes, our additional fused lasso penalty is designed to make use of this by borrowing information across outcomes for such covariates at a more granular level. Further, when there is strong group structure of whether or not a given covariate has an effect across groups of outcomes, but the magnitudes of the effects are unrelated, then both our approach and the multinomial logistic regression with sparse group lasso (MSGL) approach  [13] are appropriate options and their performances should be expected to be similar, despite our additional fused lasso penalty, whereas MRCE is not designed explicitly to make use of such information. However, if indeed there are many covariates with similar effects across outcomes and a group structure, our approach may be expected to perform well. If there is no relationship whatsoever across the multivariate outcomes, an approach of simply fitting a lasso-penalized model separately for each outcome may have reasonable performance.

§.§ Asymptotic Properties

We now present new asymptotic results for our penalization methods pertaining to semiparametric linear models. We show results for a more general, encompassing form for the overlapping group lasso penalty $P_1$ and the fused lasso penalty $P_2$. We allow $P_1(\bbeta) = \sum_{j=1}^p\sum_{g\in\calG}\lambda_{1,j,G}||\bbeta_{j,G}||_2$ for an arbitrary, potentially overlapping group structure $\calG$ and we allow $P_2(\bbeta) = \sum_{j=1}^p\sum_{(l,o) \in \calF} \lambda_{2,j,l,o}|{\beta}_{j,l} - {\beta}_{j,o}|$ for an arbitrary fusing set $\calF$ that contains pairs of indices of coefficients to be fused, where the elements of $\calF$ are of the form $(l,o)$ with $l,o\in\calK$. The penalty form introduced in Section <ref> is a special case of the form we study in this section, as our Theorem <ref> applies to any group structure, overlapping or not.

We also denote the set of all pairs of coefficients which are equal and non-zero for the $j$th variable by $\calE_{j,\cdot} = \{ (l,m) \in \mathcal{F}: \beta_{j,l}^0 = \beta_{j,m}^0 \neq 0  \}$. Similarly, define $\hat{\calE}_{j,\cdot} = \{ (l,m) : l,m \in \mathcal{F} \mbox{ and } \hat{\beta}_{j,l} = \hat{\beta}_{j,m} \neq 0  \}$. Denote the union of these sets over all variables as $\calE = \bigcup_{j=1}^p\calE_{j,\cdot}$ and $\hat{\calE} = \bigcup_{j=1}^p\hat{\calE}_{j,\cdot}$.

To facilitate our explanation of the asymptotic results, for any vector ${\bf a}_{\cdot,k}\in \mathbbm{R}^p$ associated with response $k$ and a index set ${\calI}_{\cdot,k}\subset\{1,\dots,p\}$ of size $|{\calI}_{\cdot,k}|$ associated with response $k$, ${\bf a}_{{\calI}_{\cdot,k}}$ represents a $|{\calI}_{\cdot,k}|$ dimensional sub-vector with elements in ${\bf a}_{\cdot,k}$ indexed by ${\calI}_{\cdot,k}$.
Now let $\wtbfX = \bfI_K\ \otimes \bfX  = \mbox{diag}(\bfX, \dots, \bfX)$ be the block diagonal design matrix with $K$ blocks, one for each outcome, with the $k$th block as the design matrix $\bfX$, where $\otimes$ is the Kronecker product; similarly let $\widetilde{\bbeta}^0 = \vec(\bbeta^0) = ({\betatrue_{\cdot,1}}^\top, \dots, {\betatrue_{\cdot,K}}^\top)^\top$ be the vectorization of the true coefficients for all $K$ responses, let $\wtbfY^\top = (\bfY_{\cdot,1}^\top, \dots, \bfY_{\cdot,K}^\top)$. Using this notation, we can re-express the model (<ref>) as $\wtbfY = \wtbfX \widetilde{\bbeta}^0 + \wtbepsilon$, where $\wtbepsilon$ is the stacked error vector defined similarly as $\wtbfY$.

Now let $\calJ = (\calJ_{\cdot,1}, \dots, \calJ_{\cdot,K}) \subseteq \{1, \dots, Kp\}$, where $\calJ_{\cdot,k} = \{j+(k-1)p: j\in\{1,\dots,p\} \text{ and } \beta^0_{j,k} \neq 0 \}$, to be the set of indices of all non-zero effects in $\widetilde{\bbeta}^0$ and let $\calH = \text{Hull}(\calJ) = \left\{  \cup_{G\in \mathcal{G}, G\cap \calJ = \varnothing}G  \right\}^c$ be the hull of the nonzero pattern $\calJ$, where for a set $\calS$, $\calS^c$ denotes its complement. The hull of the non-zero pattern is essentially the smallest set of groups in $\calG$ that contains all elements in $\calJ$.  Similarly, denote $\hat{\calJ} = (\hat{\calJ}_{\cdot,1}, \dots, \hat{\calJ}_{\cdot,K})$ to be the set of indices of all nonzero variables in $\betahat$. Note that by construction of $P_1(\cdot)$ due to each effect of each variable across the outcomes having its own group in $\calG_{M+1}$,  $\text{Hull}(\hat{\calJ}) = \hat{\calJ}$, however we present our results in terms of the hull so as to be fully general with respect to the group structure. We denote $\wtbfX_H$ as the columns in $\wtbfX$ corresponding to variables in $\calH$ and similarly denote $\wtbbeta^0_H$ as the values in $\wtbbeta^0$ corresponding to elements in $\calH$.
Then denote $\wtbfX_\calH^*$  constructed by dropping and collapsing columns of $\wtbfX_\calH$ corresponding to the distinct, non-zero values of $\wtbbeta^0_\calH$. Specifically, for each $j$, let $\{c_{j,1}, \dots, c_{j,L}\}$ for $L\leq K$ denote the unique non-zero values of $\{\widehat{\beta}_{j,1}, \dots, \widehat{\beta}_{j,k}\}$ and denote $\calC_{j,\ell} = \{k\in\calK:  \widehat{\beta}_{j,k} = c_{j,\ell}  \}$. Then for each $c_{j,\ell}\in\{c_{j,1}, \dots, c_{j,L}\}$
all columns $j+(k-1)p$ of $\wtbfX$ such that $k\in \calC_{j,\ell}$, $j\in\{1,\dots,p\}$ are collapsed and added together into a single column, i.e. $\sum_{k\in \calC_{j,\ell}} \wtbfX_{\cdot,j+(k-1)p}$ . As an illustration, denote $\wtbfX_{\cdot,j+(k-1)p}$ as the $j$th column of the design matrix for response $k$. Then if the coefficient for variable $j$ is in $\calH$ for responses $k, \ell$, and $K$ and $\beta^0_{j,k} = \beta^0_{j,\ell} = \beta^0_{j,K} \neq 0$, this results in the following to the corresponding columns in $\wtbfX$ in the formulation of $\bfX^*_\calH$:

\stackrel{\mbox{$\vphantom{\wtbfX_{\cdot,j+(k-1)p}}$}}{%
\begin{matrix}
\vphantom{\boldsymbol 0} \\[15pt]
\vphantom{\vdots} \\
\mbox{Block } k\rightarrow \\
\vphantom{\vdots} \\
\mbox{Block } \ell\rightarrow \\
\vphantom{\vdots} \\
\mbox{Block } K\rightarrow \\  \\
\end{matrix}
\stackrel{\mbox{$\wtbfX_{\cdot,j+(k-1)p}$}}{%
\begin{pmatrix}
\boldsymbol 0 \\
\vdots \\
\bfx_{\cdot,j} \\
\vdots \\
\boldsymbol 0 \\
\vdots \\
\boldsymbol 0 \\
\end{pmatrix}
\stackrel{\mbox{$\wtbfX_{\cdot,j+(\ell-1)p}$}}{%
\begin{pmatrix}
\boldsymbol 0 \\
\vdots \\
\boldsymbol 0 \\
\vdots \\
\bfx_{\cdot, j} \\
\vdots \\
\boldsymbol 0 \\
\end{pmatrix}
\stackrel{\mbox{$\wtbfX_{\cdot,j+(K-1)p}$}}{%
\begin{pmatrix}
\boldsymbol 0 \\
\vdots \\
\boldsymbol 0 \\
\vdots \\
\boldsymbol 0 \\
\vdots \\
\bfx_{\cdot, j} \\
\end{pmatrix}
\rightarrow
\begin{pmatrix}
\boldsymbol 0 \\
\vdots \\
\bfx_{\cdot, j} \\
\vdots \\
\bfx_{\cdot, j} \\
\vdots \\
\bfx_{\cdot, j} \\
\end{pmatrix}

In other words, the columns in $\wtbfX$ corresponding to these three variables are collapsed and added together.  Note that there exist matrices $\bfH$ and $\bfE$ such that $\wtbfX_\calH^* = \wtbfX\bfH\bfE$, where $\bfH$ is formed by removing all columns of positions in the indices of $\calH$ from the identity matrix of dimension $pK\times pK$ and $\bfE$ is formed by taking the identity matrix of dimension $|\calH|\times |\calH|$ and collapsing and summing columns corresponding to non-zero coefficients that are equal to each other and are in $\calF$, the fusing set used in $P_2$.

We assume the following standard regularity conditions

(D.1) $\lim_{N \to \infty} \frac{1}{N}\bfX^\top \bfX \to \bQ$ where $\bQ$ is positive definite

(D.2) The errors $\boldsymbol\epsilon_i$ for all samples $i=1\dots,n$ are $i.i.d.$ random vectors with mean zero and finite, positive definite variance-covariance matrix $\bSigma$.

The following result pertains to cases where the group structure has been misspecified. As such, we define ${\bbeta^*}^0$ to be the distinct values of ${\bbeta}^0_\calH$ (i.e. the union of the distinct values in $\{\beta^0_{j,k},\beta^0_{j,\ell}: (k,\ell) \in \calE \text{ or } (\ell,k) \in \calE \text{ and }  (k,\ell) \in \calF \text{ or } (\ell,k) \in \calF\}$ and all remaining values $\beta^0_{j,k} \neq 0$)  appended with zeros corresponding to the elements of ${\bbeta}^0_{\calH^c}$. Similarly, we define $\betahat^*$ to be the distinct values of ${\betahat}_{\hat{\calJ}}$ appended with zeros corresponding to the elements of ${\betahat}_{\hat{\calJ}^c}$.

Assume the data are generated under the model described in (<ref>) and that our estimator is given by (<ref>). Furthermore,
assume conditions (D.1) and (D.2) and let $\lambda_{1,j,G} = ||\betahat^{OLS}_{j,G}||_2^{-\gamma_1}$ and $\lambda_{2,j,l,m} = |\hat{\beta}_{j,l}^{OLS} - \hat{\beta}_{j,m}^{OLS}|^{-\gamma_2}$, where $\gamma_1, \gamma_2 > 0$ such that $N^{(\gamma_\ell + 1)/2}\lambda \to \infty$ for $\ell=1,2$. If $\sqrt{N}\lambda \to 0$, then we have the following:

\begin{align}
& P(\hat{\calJ}_{\cdot,j} = {\calH}_{\cdot,j}) \to 1 \mbox{ as }  N \to \infty,  \label{linear_model_selection_consistency_mis} \\
& P(\hat{\calE}_{\cdot,j} = {\calE}_{\cdot,j}) \to 1  \mbox{ as }  N \to \infty, \label{linear_model_fused_selection_consistency_mis}
\end{align}

for each $j=1,\dots,p$	and

\begin{align}
\sqrt{N}(\betahat^* - {\bbeta^*}^0) \xrightarrow{d} \bfZ^*, \label{linear_model_asympt_distr_mis}
\end{align}

where $\bfZ^* = ({{}\bfZ_{\calH}^*}^\top, {\boldsymbol 0}^\top)^\top$  and $\bfZ_{\calH}^* \sim N_{|\calH^*|}(0,  {\bfQ^*}_{\calH}^{-1}{\bfV^*}_{\calH}{\bfQ^*}_{\calH}^{-1} )$, where  $\bfV^*_\calH = \bfE^\top\bfH^\top \left(\bSigma \otimes \bfQ\right)\bfH\bfE$, $\bfQ^*_\calH = \bfE^\top\bfH^\top \left(\bfI_K \otimes \bfQ\right)\bfH\bfE$, and $|\calH^*|$ denotes the number of distinct non-zero elements in the true coefficients ${\bbeta^*}^0$.

The proof of Theorem <ref> is provided in the Supplementary Material Appendix A. We reiterate that Theorem <ref> applies to any group structure, including those with overlapping groups.
We note that the variance-covariance matrix in the distribution of $\bfZ^*$ does not simplify due to the allowance of selection of variables individually by outcome in combination with the Kronecker-product variance structure of the response vector $\wtbfY$; without either it would simplify considerably. The terms $\bfH$ and $\bfE$ in the asymptotic variance of the estimates are due to selection of variables and collapsing of effects of a single variable across outcomes, respectively; without any selection these matrices are removed and the variance simplifies to the usual asymptotic variance of multivariate response regression via least squares. However, the asymptotic variance matches that of the estimator where all and only all truly non-zero variable effects are included in a model and all variables with equal effects across a subset of outcomes are accordingly collapsed. The result of Theorem <ref> does not require that $\calE \subseteq \calF$, however, if indeed $\calE \subseteq \calF$, i.e.
if $\calF$ contains indices corresponding to all pairs of coefficients in $\bbeta^0_{\calH}$ that are equal to each other, then due to (<ref>), when using adaptive terms for the tuning parameters, our approach yields selection consistency for the hull of the non-zero terms; when the group structure has individual groups for each of the individual coefficients in the model, the convex hull of the non-zero coefficients is simply the set of non-zero coefficients. Thus, when such individual groups are included, we have selection consistency for the non-zero coefficients. Our theory also shows that with probability tending to one, our approach will fuse together the coefficients that are truly equal to each other. Finally, our theory shows that our estimate will converge to the same asymptotic distribution as if we had known which coefficients were non-zero and which coefficients were equal to each other.

§ NUMERICAL EXPERIMENTS

§.§ Estimators used

In this section, we conduct simulation experiments to assess the small sample operating characteristics of our proposed doubly-structured sparsity inducing estimation approach in comparison with several other state-of-the-art approaches for multivariate regression in high-dimensions. In our simulations, we vary both the dimensionality of the problem and the sample size. We also consider data-generating mechanisms with varying degrees of sparsity of variable effects that align with grouping of the outcomes and further consider scenarios where the effects of each variable have either minimal similarity across outcomes or a varying degree of similarity across outcomes. This allows our studies to explore under what data-generating mechanisms our two penalties indeed help in estimation. We use adaptive and non adaptive versions of our method, denoted OGFM(adapt) and OGFM, respectively, where OGFM indicates Overlapping Group $+$ Fused Multivariate regression. For the adaptive version, we set $\gamma_1=\gamma_2=0.5$ and for high dimensional settings with $p\geq n$, we used marginal regression estimates for the adaptive weights as in Huang, et al. [27].
We compare our approach with the approach of Li, et al. [13] (denoted as MSGL), which allows for overlapping group lasso penalties. For MSGL, we use the   package version 2.1 We further compare with a simple approach of fitting a separate lasso-penalized linear regression model for each outcome, with the tuning parameter chosen separately for each outcome (denoted as Sep-Lasso) and do so using the   package version 4.1-3. We also compare with the approach of Rothman, et al. [8] (denoted MRCE) implemented in the    package version 2.1. For all methods, the tuning parameters are chosen by 10-fold cross validation; MSGL has two tuning parameters (one for a lasso penalty, one for a group lasso penalty), MRCE has two tuning parameters (one for a lasso penalty for the coefficients for the variable effects on responses, another for a lasso penalty on the elements of the inverse covariance matrix of the residuals), Sep-Lasso, for separate lasso has a tuning parameter for each outcome for a lasso penalty on coefficients, and the OGFM approaches have two tuning parameters (one for the overlapping group lasso penalty and another for the fused lasso penalty).  For both the OGFM approaches and MSGL, the group lasso penalty applied corresponds to the true underlying group structure of the data-generating mechanism.

§.§ Data generation

For each replication of the simulation, we generate data under model (<ref>), where $\bfx_i$ are generated as i.i.d. multivariate normal random variables with mean vector $\bzero$ and covariance matrix $\bSigma_{\bfX} = [\sigma_{Xjk}]_{j,k=1}^p$, where $\sigma_{Xjk} = 0.5^{|j-k|}$, as was used in the simulations of Yuan and Lin [9] and Rothman, et al. [8]. The responses are generated according to model (<ref>), where the error vectors $\bepsilon_i$ are generated as i.i.d. multivariate normal random variables with mean vector $\bzero$ and covariance matrix $\bSigma_{\bepsilon} = [\sigma_{\epsilon jk}]_{j,k=1}^K$, where $\sigma_{\epsilon jk} = 4\left(0.5^{|j-k|}\right)$. In the Supplementary Material Appendix B.2, we explore a simulation setting with binary covariates and a Likert scale outcome similar to our motivating data.
The dimensionality of the outcome/error vector is $K=8$ and the outcomes form 3 groups, with the first three outcomes forming group 1, the fourth and fifth forming group 2, and the last three outcomes forming group 3. The variable effects are generated as

\begin{equation*}
\bbeta^0 =
\begin{pmatrix}
\bovermat{$G_{1,1}$}{ \xi_{1,1} \xi^G_{1,1} \eta_{1,1} & \xi_{1,2} \xi^G_{1,1} \eta_{1,2} & \xi_{1,3} \xi^G_{1,1} \eta_{1,3} } & \bovermat{$G_{1,2}$}{\xi_{1,4} \xi^G_{1,2} \eta_{1,4} & \xi_{1,5} \xi^G_{1,2} \eta_{1,5}} & \bovermat{$G_{1,3}$}{\xi_{1,6} \xi^G_{1,3} \eta_{1,6} & \xi_{1,7} \xi^G_{1,3} \eta_{1,7} & \xi_{1,8} \xi^G_{1,3} \eta_{1,8}} \\
\vdots & & & & & & & \vdots \\
\xi_{z,1} \xi^G_{z,1} \eta_{z,1} & \xi_{z,2} \xi^G_{z,1} \eta_{z,2} & \xi_{z,3} \xi^G_{z,1} \eta_{z,3} & \xi_{z,4} \xi^G_{z,2} \eta_{z,4} & \xi_{z,5} \xi^G_{z,2} \eta_{z,5} & \xi_{z,6} \xi^G_{z,3} \eta_{z,6} & \xi_{z,7} \xi^G_{z,3} \eta_{z,7} & \xi_{z,8} \xi^G_{z,3} \eta_{z,8} \\
0 & \cdots & & & & &  \cdots & 0 \\
\vdots & \ddots & & & & & & \vdots \\
0 & \cdots & & & & &  \cdots & 0
\end{pmatrix},
\end{equation*}

where the last $p-z$ rows of $\bbeta^0$ have all elements as 0, the terms $\xi_{j,k}\sim$ Bernoulli$(0.9)$ induce sparsity at the individual effect level, the terms $\xi^G_{j,1}, \xi^G_{j,2}, \xi^G_{j,3} \sim$ Bernoulli$(1-p_{\text{HS}})$ induce sparsity at the group level, the variable effect size terms $\eta_{j,k}$ are distributed i.i.d. uniformly from $\{-1, -0.5, -0.25, -0.125, 0.125, 0.25, 0.5, 1\}$ to create both small and large effects, and for each variable $j$ separately the terms $\eta_{j,k}$ for $k$ in the same group are set to be all equal to each other with probability $p_{\text{GE}}/2$ and independently, the terms $\eta_{j,k}$ for all outcomes $k=1,\dots,8$ are set to be all equal to $\eta_{j,1}$ with probability $p_{\text{GE}}/2$. The latter process induces effects for some variables within a group to be equal to each other and induces effects for some variables to be equal for all outcomes.

We explore dimensions of $p=50, 100,$ and $200$ and set $z=25, 50, 50$ for each dimension setting, respectively. We explore hierarchical sparsity parameters from $p_{\text{HS}} = 0, 0.25$ and $0.5$ and fusing probabilities $p_{\text{GE}} = 0, 0.5, 0.75$ and $0.95$; when $p_{\text{HS}} = 0$, there is no group-wise sparsity and thus any group lasso penalty applied is superfluous. For each replication of the simulation, we additionally generate an independent test set of size 10000 for use in evaluation of predictive performance.

§.§ Performance evaluation

We evaluate methods in terms of three metrics: root mean squared error of predictions (RMSE) on a large, independent test set of size 10000 generated anew for each simulation replication, model error defined below, and balanced accuracy which is the average of the true positive rate (TPR) and true negative rate (TNR), described in detail in the Supplementary Material Appendix B.1 and not to be confused with standard classification accuracy. As our primary interest is performance of predictions on validation data,  it is our primary focus of evaluation. For the RMSE, we compute the RMSE for each outcome separately and then average the RMSEs across the outcomes. Model error for a given estimate is defined as ME$(\widehat{\bbeta}, \bbeta^0) = \text{tr}\left[  (\widehat{\bbeta} - \bbeta^0)^\top \bSigma_{\bfX} (\widehat{\bbeta} - \bbeta^0) \right]$.

§.§ Results

The RMSE results fixing the parameter at $p_{\text{HS}}=0$ are displayed in Figure <ref> and the RMSE results for $p_{\text{HS}}=0.25$ are displayed in Figure <ref>. When the sample size is small ($n=100$) and the hierarchical sparsity probability is 0 (Figure <ref>), our proposed approach OGFM results in the smallest validation RMSE on average across the replications across all settings including moderate dimensions ($p=50,100$) and high dimensions ($p=200$), with MSGL as second best across the majority of settings and OGFM(adapt) a close third. The separate lasso tends to perform worse than all competing approaches, except for in some settings (e.g. when $n=400$) where it performs marginally better than MRCE. We note that the code for MRCE throws an error when $p\geq n$, although in principle MRCE can work for high dimensional settings. In the same settings but with $n=200$, OGFM still tends to work the best across all settings, however when the probability of coefficients being equal to each other is zero ($p_{\text{GE}}=0$), and the dimensionality is moderate, the MRCE approach performs nearly as well. In the high dimensional setting ($p=200$), MSGL performs on par with OGFM and sometimes better, however this benefit slightly attenuates for the larger sample size setting ($n=400$). In general the adaptive version of our approach, OGFM(adapt) performs worse than OGFM in small sample size settings and performs better with larger sample sizes. The results with the hierarchical sparsity probability set to 0.25 (Figure <ref>) roughly mirror the results with no hierarchical sparsity, however the performance for all methods is slightly better in terms of RMSE as there are fewer overall non-zero coefficients in the data-generating process.

The model error results largely track with the validation RMSE results, albeit on a different scale, and are thus shown in the Supplementary Material Appendix B.1.

Validation RMSEs for all methods across 100 replications of the simulation experiment holding the parameter $p_{\text{HS}}=0$ so that there is no group-level sparsity.

Validation RMSEs for all methods across 100 replications of the simulation experiment holding the parameter $p_{\text{HS}}=0.25$ so that there is moderate group-level sparsity.

Figure <ref> shows the same results as Figures <ref> and <ref>, but is re-organized to focus on the effect of the fused coefficients probability ($p_{\text{GE}}$) and the hierarchical sparsity probability ($p_{\text{HS}}$). In this figure, we fix the sample size to be 200. In general, we can see that as the probability of coefficients being equal/fused within groups increases, the performance of OGFM and OGFM(adapt) relative to the Sep-Lasso, MSGL, and MRCE tends to improve, with the trend being pronounced in moderate dimensional settings. Only the OGFM methods improve as the probability of coefficients being equal increases, as they are the only approaches which explictly allow for shrinkage of effect sizes to each other across outcomes. The relative performance of all methods tends to stay relatively consistent as $p_{\text{HS}}$ is varied. However, we note that MRCE tends to perform best in the setting with the fewest true non-zero coefficients ($p_{\text{HS}}=0.25$) and smallest amount of truly equal coefficients ($p_{\text{GE}}=0$).

In the Supplementary Material Appendix B.1, we show simulation results in terms of the average of the TPR and TNR. From these results, it can be seen that while OGFM performs best in terms of prediction performance, it tends to over-select terms, resulting in a high TPR but low TNR and thus lower average TPR and TNR value. On the other hand, for larger sample sizes OGFM(adapt) performs very well in terms of the average TPR and TNR as it selects fewer variables with a higher proportion of selected variables being ones with truly non-zero coefficients. Overall, MSGL performs best in terms of the average TPR and TNR and MRCE worst, the latter trend having been concurrently observed in the original work of Rothman, et al. [8], who noted that MRCE tends to result in better model error but with less benefit in terms of TPR and TNR. In the Supplementary Material Appendix B.1 we also present computation times for all methods. The OGFM approaches are highly competitive computationally as the sample size increases, however their performance deteriorates as the dimensionality increases.

Validation RMSEs for all methods across 100 replications of the simulation experiment holding fixing the sample size $n=200$ and varying all other simulation parameters. The points are the average RMSEs across the 100 replications and error bars are plus and minus 1 standard deviation of this average.

§ ANALYSIS OF PEDIATRIC FUNCTIONAL STATUS DATA

Our motivating study consists of 1897 children, adolescents, and young adults who were admitted to an inpatient rehabilitation unit with a diagnosis of neurologic injury or illness at a major Midwestern children's hospital for neurological injury or illness between the years 2000 and August of 2020. The sample consisted of individuals between the ages of six months to 32 years with a median age of 12 years and 1st and 3rd quartiles of 6 and 16 years, respectively. The covariates used as predictors are (numbering 608 in total) included basic sociodemographic information including age and sex, and billing codes including CPT codes, ICD-9 and -10 codes, pharmaceutical codes, and data indicating durable medical equipment (DME) use. The billing codes represent information from a single episode of care involving an admission to an inpatient rehabilitation unit. The billing code information pertain only to information collected prior to the WeeFIM assessment. CPT and all ICD-9 and -10 codes were linked together to unique concept unique identifiers (CUIs) using the UMLS metathesaurus mapping system [28]. Concept Unique Identifiers (CUIs) are distinct medical concepts (codes, diseases, etc) identified by the UMLS metathesaurus. By linking all billing codes to CUIs, we are able to in many cases handle the switch-over from ICD-9 to ICD-10 codes. An added benefit of CUIs is that they can be linked to other coding systems.
We utilize all 18 WeeFIM® components as responses; these 18 components are divided into three domains. The self-care domain describes how well a child is able to feed themselves, groom, bath, dress, and complete toileting tasks including the management of their bowel and bladder. The mobility domain describes how well a child is able to transfer on and off a toilet, in and out of a bathtub or shower, and in and out of a chair or wheelchair. The mobility domain also describes a child's ability to walk, crawl, or use a wheelchair, and to move up or down stairs. The cognition domain describes how well a child can express themselves, understand information, interact with peers, solve daily problems, and recall information. Together, they describe the ability of children to function in routine and important aspects of daily life.

As our motivating use case of developing a model for the WeeFIM® components is to apply it to future data to assess functional ability across a health system population, we validate all developed models by splitting data into training and validation sub-datasets, where the training dataset is from years prior to data from the validation dataset. We use data prior to 2018 for training and data from 2018 to 2020 as validation data, leaving 1592 observations for training and 305 observations for validation.  For validation, we compare methods in terms of the response-specific mean squared error for each response and the average validated mean squared error across all 18 responses. We also compute the corresponding validation R-squared values, i.e. the proportion of the validation responses explained by the out-of-sample predictions.
We exclude any covariates that have no variation in either the training or validation datasets. After this screening, the combined dimensionality of all the predictors was 608, a high-dimensional scenario given the sample size.

We apply all methods described in Section <ref> and use the approaches described therein for selection of tuning parameters. For OGFM(adapt), we use $\gamma_1=\gamma_2=0.5$, as use of $\gamma_1=\gamma_2=1$ resulted in some extreme weights given the high dimensionality of the data. For both the OGFM approaches and MSGL, the group lasso penalty applied corresponds to the three pre-defined domains of the WeeFIM® components (self-care, mobility, and cognition).

The validation MSEs and R-squared values are displayed in Tables <ref> and <ref>, respectively. Our OGFM approach has the best performance on the validation data on average across the 18 responses and also most often performs best for the individual responses, with the adaptive version of OGFM with the second lowest validation MSE on average. For responses where OGFM does not perform best, most often its performance does not differ much from the best MSE across the remaining methods, excluding OGFM(adapt). When OGFM performs best for particular responses, its improvement in terms of RMSE and R-squared is often large, as is depicted in Figure <ref>, which shows the difference between OGFM and competing approaches in terms of MSE on the validation data.

For both OGFM approaches, the fused lasso mixing parameter $\alpha$ was chosen by cross-validation to be $1\times10^{-5}$, indicating only a small amount of fused lasso was required by the data. However, this small amount of fused lasso penalty had the effect of shrinking a large number of coefficients either close to each other or exactly to each other. In Figure <ref>, we display the coefficient paths versus $\lambda$ (described in Section <ref>) for 6 variables, where each plot is the coefficient path for a particular variable across the 18 outcomes. For each of the coefficient path plots, we fix the mixing tuning parameter $\alpha$ at the value that minimizes the cross validation error. The patterns of these coefficient paths demonstrate that our approach fuses coefficients for a variable to be the same across multiple responses when appropriate, allows effects to differentiate when warranted by the data, an allows for joint selection of effects across an entire domain of related responses. For the ICD code “shock, unspecified”, the fused lasso grouped the effects for all cognition outcomes together along the entire path, all mobility coefficients together, and grouped the self-care coefficients into two groups. As can be seen for other variables, different amounts fusing occurred for different groups of outcomes, and for some outcomes effects were shrunk towards each other but without exact fusing.

Domain      8cSelf-Care WeeFIM Responses

Method     Average     Ea     G     B     D-U     D-L     T-g     Bl-M     Bo-M

OGFM     2.153     2.074     1.816     2.034     1.673     2.099     2.052     3.533     3.492

OGFM(adapt)     2.165     2.070     1.824     2.032     1.663     2.112     2.165     3.573     3.521

MRCE     2.185     2.155     1.864     2.093     1.674     2.117     2.093     3.735     3.638

MSGL     2.222     2.145     1.860     2.149     1.712     2.218     2.190     3.729     3.617

Sep-Lasso     2.213     2.142     1.860     2.189     1.672     2.225     2.148     3.788     3.641

Domain     5|c|Mobility WeeFIM Responses     5cCognition WeeFIM Responses

C-W      T     T-S     W-W     1r|S     C     Ex     S-I     P-S     M

OGFM         1.483      1.770     1.811     1.344     1r|2.169     2.140     2.165     2.159     2.537     2.395

OGFM(adapt)        1.492     1.792     1.808     1.336     1r|2.187     2.155     2.159     2.180     2.495     2.407

MRCE        1.571      1.811     1.889     1.331     1r|2.164     2.126     2.081     2.140     2.446     2.399

MSGL        1.494     1.857     1.790     1.305     1r|2.182     2.156     2.185     2.264     2.685     2.453

Sep-Lasso        1.539     1.809     1.915     1.286     1r|2.186     2.178     2.152     2.174     2.491     2.432

Mean squared errors (MSEs) on the 18 WeeFIM components in the validation data and average MSE across the 18 components. Bold indicates the smallest MSE across the different methods for a particular component. Ea: Eating, G: Grooming, B: Bathing, D-U: Dressing Upper, D-L: Dressing Lower, T-g: Toileting, Bl-M: Bladder Management, Bo-M: Bowel Management, C-W: Chair Wheelchair, T: Toilet, T-S: Tub Shower, W-W: Walk Wheelchair, S: Stairs, C: Comprehension, Ex: Expression, S-I: Social Interaction, P-S: Problem Solving, M: Memory.

Comparisons of the validation MSE for each of the 18 WeeFIM components and on average (left) across the 18 components

Domain      8cSelf-Care WeeFIM Responses

Method     Average     Ea     G     B     D-U     D-L     T-g     Bl-M     Bo-M

OGFM     0.432     0.526     0.536     0.476     0.446     0.447     0.498     0.449     0.429

OGFM(adapt)     0.429     0.527     0.534     0.477     0.449     0.444     0.470     0.443     0.425

MRCE     0.425     0.507     0.524     0.461     0.446     0.442     0.488     0.418     0.405

MSGL     0.416     0.509     0.525     0.447     0.433     0.416     0.464     0.419     0.409

Sep-Lasso     0.419     0.510     0.525     0.437     0.446     0.414     0.474     0.410     0.405

Domain     5|c|Mobility WeeFIM Responses     5cCognition WeeFIM Responses

C-W     T     T-S     W-W     1r|S     C     Ex     S-I     P-S     M

OGFM         0.433      0.522     0.447     0.287     1r|0.298     0.410     0.405     0.400     0.375     0.384

OGFM(adapt)         0.430     0.516     0.448     0.290     1r|0.292     0.406     0.406     0.394     0.385     0.381

MRCE        0.400     0.510     0.423     0.293     1r|0.299     0.414     0.428     0.405     0.397     0.383

MSGL        0.429     0.498     0.454     0.307     1r|0.294     0.406     0.399     0.371     0.338     0.369

Sep-Lasso        0.412     0.511     0.415     0.317     1r|0.292     0.400     0.408     0.396     0.386     0.374

R-squared values on the 18 WeeFIM components in the validation data and average R-squared value across the 18 components. Bold indicates the larger R-squared across the different methods for a particular components.

From the coefficient paths in Figure <ref>, our approach yields benefit in terms of predictive performance and also results in clinically sensible estimated coefficients. We chose these 6 variables to demonstrate effects of several classes of billing codes in the data: ICD codes, DME codes, pharmaceutical codes, and CPT codes. The three figures on the left of Figure <ref> represent examples of diagnostic codes from billing data. The top panel represents the functional domains’ coefficients path plot for children coded with autistic disorder.  Children with autism spectrum disorders are diagnosed based on deficits in areas of social communication and behavior. The degree of impairment in social skills as well as comorbid neurodevelopmental disorders often leads to functional deficits in self-care.  Mobility is typically unaffected in children with autism. Another example is the diagnosis code for shock. Shock is a life-threatening condition in which the body is not getting enough blood flow, resulting in end-organ damage, including damage to the brain.   Consistent with this, in the figure we see from the coefficient path plots that the diagnosis code of shock is most related  with the functional domain of  cognition/communication as compared to mobility and self-care skills. Hypothyroidism is a disorder caused when the thyroid gland does not make enough hormone. In children, it has numerous etiologies including autoimmune reactions, brain injury, or radiation treatment. Hypothyroidism causes slow growth, decreased strength, and impaired cognition. The associated plot shows some relationship to the self-care and mobility domains, but a stronger relationship to the cognition/communication domain. The three figures in the right of the panel represent billing codes for a medication, a durable good, and a procedure. Methylphenidate Hydrochloride, a stimulant medication used to treat attention deficit hyperactivity disorder, showed minimal relationships with mobility and self-care.  However a pronounced relationship with cognition/communication is observed, as would be expected in a child with being treated with this medication. An example of durable medical equipment are enteral feeding supplies. These are ordered for patients who require feeding through a tube inserted into the stomach. This plot shows a clear relationship between feeding supplies and eating outcomes, but markedly less strong relationship with other functional outcome domains. Mechanical chest wall oscillation is used in patients who have an impaired lung function such as in children with cystic fibrosis or other diseases which lead to general physicial debilitation or loss of strength. The plot demonstrates that this procedure code holds associations with  mobility and self-care limitations as would be expected, but does not relate to cognition/communication.

Coefficient path plots across all 18 WeeFIM components for six variables from the pediatric functional status data. The dashed vertical line indicates the value of $\lambda$ that minimizes the cross-validation error.

§ DISCUSSION

In this paper we introduced a doubly structured sparsity regression approach for modeling multivariate outcomes with known groupings to borrow information across related outcomes. Our approach allows for overlapping groupings of the responses and borrows strength by joint selection of the effects of variables across groups of related outcomes and encourages shrinkage of the effects of variables to be more similar across related outcomes. We prove an oracle property for an adaptive version of our penalty, showing that our approach results in estimates with the same asymptotic distribution as if the true non-zero coefficients and which variable effects are truly equal across related responses were both known in advance. The results pertaining to our motivating study indicate our approach yields benefit in terms of predictive performance and yields estimated coefficients that are clinically sensible, as was demonstrated in our analysis of the pediatric functional status data.
Our work is motivated by modeling functional status scores in a pediatric population however the methods are applicable more generally. Since our approach does not rely on normality of the multivariate outcomes, it can be straightforwardly extended to more general models, such as generalized linear models both in theory and computationally. Incorporation of covariance information in addition to borrowing strength through joint selection and effect shrinkage may be a promising direction for additionally borrowing strength of information across outcomes. However, in this work we have not focused on this as it would introduce another tuning parameter and another layer of complexity to computation, which may limit the applicability of such an approach to large scale and high dimensional settings.

§ ACKNOWLEDGMENTS

Research reported in this publication was supported by the Eunice Kennedy Shriver National Institute of Child Health & Human Development of the National Institutes of Health under Award Number R03HD101083. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

[1]
Lo W, Stephens J, Fernandez S. Pediatric stroke in the United States and the
impact of risk factors. Journal of Child Neurology 2009;
24(2): 194–203.

[2]
Patel S, Bhatnagar A, Wear C, et al. Are pediatric brain tumors on the rise in
the USA? Significant incidence and survival findings from the SEER
database analysis. Child's Nervous System 2014; 30(1):

[3]
Taylor CA, Bell JM, Breiding MJ, Xu L. Traumatic brain injury–related
emergency department visits, hospitalizations, and deaths—United States,
2007 and 2013. MMWR Surveillance Summaries 2017; 66(9):

[4]
Dhillon JK, Shi J, Janezic A, Wheeler KK, Xiang H, Leonard JC. US estimates
of pediatric spinal cord injury: implications for clinical care and research
planning. Journal of Neurotrauma 2017; 34(12):

[5]
Schneier AJ, Shields BJ, Hostetler SG, Xiang H, Smith GA. Incidence of
pediatric traumatic brain injury and associated hospital resource utilization
in the United States. Pediatrics 2006; 118(2): 483–492.

[6]
Székely GJ, Rizzo ML, Bakirov NK. Measuring and testing dependence by
correlation of distances. The Annals of Statistics 2007;
35(6): 2769–2794.

[7]
Tibshirani R. Regression shrinkage and selection via the lasso. Journal of
the Royal Statistical Society. Series B (Methodological) 1996;
58: 267–288.

[8]
Rothman AJ, Levina E, Zhu J. Sparse multivariate regression with covariance
estimation. Journal of Computational and Graphical Statistics
2010; 19(4): 947–962.

[9]
Yuan M, Lin Y. Model selection and estimation in the Gaussian graphical model.
Biometrika 2007; 94(1): 19–35.

[10]
Sofer T, Dicker L, Lin X. Variable selection for high dimensional multivariate
outcomes. Statistica Sinica 2014; 24(4): 1633.

[11]
Fan J, Li R. Variable selection via nonconcave penalized likelihood and its
oracle properties. Journal of the American Statistical Association
2001; 96(456): 1348–1360.

[12]
Zhang CH. Nearly unbiased variable selection under minimax concave penalty.
The Annals of Statistics 2010; 38(2): 894–942.

[13]
Li Y, Nan B, Zhu J. Multivariate sparse group lasso for the multivariate
multiple linear regression with an arbitrary group structure. Biometrics 2015; 71(2): 354–363.

[14]
Jenatton R, Audibert JY, Bach F. Structured variable selection with
sparsity-inducing norms. The Journal of Machine Learning Research
2011; 12: 2777–2824.

[15]
Yuan M, Lin Y. Model selection and estimation in regression with grouped
variables. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 2006; 68(1): 49–67.

[16]
Tibshirani R, Saunders M, Rosset S, Zhu J, Knight K. Sparsity and smoothness
via the fused lasso. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 2005; 67(1): 91–108.

[17]
Zou H. The Adaptive Lasso and Its Oracle Properties. Journal of the
American Statistical Association 2006; 101(476): 1418–1429.

[18]
Wang H, Leng C. A note on adaptive group lasso. Computational Statistics
& Data Analysis 2008; 52(12): 5277–5286.

[19]
Huling JD, Yu M, Liang M, Smith M. Risk prediction for heterogeneous
populations with application to hospital admission prediction. Biometrics 2018; 74(2): 557–565.

[20]
Viallon V, Lambert-Lacroix S, Höfling H, Picard F. Adaptive generalized
fused-lasso: Asymptotic properties and applications. Technical Report

[21]
Guennebaud G, Jacob B, others . Eigen v3. http://eigen.tuxfamily.org;  2010.

[22]
Bates D, Eddelbuettel D. Fast and Elegant Numerical Linear Algebra Using the
RcppEigen Package. Journal of Statistical Software 2013;
52(5): 1–24.

[23]
Glowinski R, Marroco A. Sur l'approximation, par éléments finis d'ordre
un, et la résolution, par pénalisation-dualité d'une classe de
problèmes de Dirichlet non linéaires. Revue française
d'automatique, informatique, recherche opérationnelle. Analyse
numérique 1975; 9(2): 41–76.

[24]
Gabay D, Mercier B. A dual algorithm for the solution of nonlinear variational
problems via finite element approximation. Computers & Mathematics with
Applications 1976; 2(1): 17–40.

[25]
Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and
statistical learning via the alternating direction method of multipliers.
Foundations and Trends® in Machine Learning
2011; 3(1): 1–122.

[26]
Chen C, He B, Ye Y, Yuan X. The direct extension of ADMM for multi-block
convex minimization problems is not necessarily convergent. Mathematical
Programming 2016; 155(1-2): 57–79.

[27]
Huang J, Ma S, Zhang CH. Adaptive Lasso for sparse high-dimensional regression
models. Statistica Sinica 2008: 1603–1618.

[28]
Bodenreider O. The unified medical language system (UMLS): integrating
biomedical terminology. Nucleic Acids Research 2004;
32(suppl_1): D267–D270.

§ SUPPORTING INFORMATION

Web Appendix A, referenced in Section <ref> and Web Appendix B, referenced in Section <ref>, and code for our simulation studies are available online.