# On the Use of Minimum Penalties in Statistical Learning

Ben Sherwood
School of Business, University of Kansas
and
Bradley S. Price
John Chambers College of Business and Economics, West Virginia University  The
authors gratefully acknowledge the support of the West Virginia Clinical and
Translational Science Institute (NIH/NIGMS Award Number U54GM104942) and note
that this work is partially supported by by the National Science Foundation
(NSF) Major Research Instrumentation Program (MRI Award 1726534) and the John
Chambers College of Business and Economics Summer Research Program.

###### Abstract

Modern multivariate machine learning and statistical methodologies estimate
parameters of interest while leveraging prior knowledge of the association
between outcome variables. The methods that do allow for estimation of
relationships do so typically through an error covariance matrix in
multivariate regression which does not scale to other types of models. In this
article we proposed the MinPEN framework to simultaneously estimate regression
coefficients associated with the multivariate regression model and the
relationships between outcome variables using mild assumptions. The MinPen
framework utilizes a novel penalty based on the minimum function to exploit
detected relationships between responses. An iterative algorithm that
generalizes current state of the art methods is proposed as a solution to the
non-convex optimization that is required to obtain estimates. Theoretical
results such as high dimensional convergence rates, model selection
consistency, and a framework for post selection inference are provided. We
extend the proposed MinPen framework to other exponential family loss
functions, with a specific focus on multiple binomial responses. Tuning
parameter selection is also addressed. Finally, simulations and two data
examples are presented to show the finite sample properties of this framework.

Keywords: Non-Convex Optimization, Post-Selection Inference, High Dimensional
Convergence, Graph Constrained Models, Selection Consistency

##  1 Introduction

Multivariate, also known as multiple response or multi-task, regression
jointly models multiple responses (outcomes) given a common set of predictors
(inputs). Joint modeling of responses, as opposed to modeling each response
separately, is only interesting if the responses are related. Leveraging these
relationships should improve estimation of the model, but in practice these
relationships may not always be known. In this work we propose a framework to
simultaneously estimate the regression coefficients and detect relationships
between the response variables. The framework also allows for simultaneous
variable selection.

Consider the sample of independent data, $\\{{\bf x}_{i},{\bf
y}_{i}\\}_{i=1}^{n}$, where ${\bf x}_{i}\in\mathbb{R}^{p}$ and ${\bf
y}_{i}\in\mathbb{R}^{r}$ with a model of

${\bf y}_{i}={\bf B}^{*T}{\bf x}_{i}+\mbox{\boldmath$\epsilon$}_{i},$ (1)

where ${\bf
B}^{*}=(\mbox{\boldmath$\beta$}_{1}^{*},\ldots,\mbox{\boldmath$\beta$}_{r}^{*})\in\mathbb{R}^{p\times
r}$ and $\mbox{\boldmath$\epsilon$}_{i}\in\mathbb{R}^{r}$ are i.i.d. random
vectors with mean zero and covariance matrix $\Sigma$. A popular estimator for
${\bf B}^{*}$ is

$\widehat{{\bf
B}}=\mathop{arg\,min}\displaylimits_{\mbox{\boldmath$\beta$}_{1},\ldots,\mbox{\boldmath$\beta$}_{r}\in\mathbb{R}^{p}}\frac{1}{2n}\sum_{k=1}^{r}\sum_{i=1}^{n}(y_{ik}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})^{2},$ (2)

which has desirable qualities such as consistency and asymptotic normality
under mild conditions, and is the maximum likelihood estimator if the errors
are from a multivariate normal distribution. However, there is something
unsatisfying about not using correlation of the responses in the estimation.
To improve efficiency in estimation of coefficients, while accounting for
correlation across the errors Rothman et al. (2010) (MRCE), Witten &
Tibshirani (2009) (SCOUT) and Lee & Liu (2012) have all proposed methods which
simultaneously estimates the regression coefficients and the inverse
covariance matrix of the errors. These methods all aim to improve estimation
of regression coefficients by exploiting correlation across the responses that
is not explained by the predictors.

We propose an alternative approach that focuses on relationships across the
responses that can be explained, due to similarity in the regression
coefficients. Specifically, we consider the kth and lth response related if
$\mbox{\boldmath$\beta$}_{k}^{*}$ is similar to
$\mbox{\boldmath$\beta$}_{l}^{*}$ or $-\mbox{\boldmath$\beta$}_{l}^{*}$. If
relationships were known _a priori_ , then for each response $k$ one could
define three disjoint sets: (1) $P_{k}$, a set of responses positively related
to response $k$; (2) $N_{k}$, a set of responses negatively related to
response $k$; and (3) $Z_{k}$ a set of responses that are not related to $k$.
These relationships could be used to improve estimation using a penalized
likelihood such as

$\mathop{arg\,min}\displaylimits_{{\bf B}\in\mathbb{R}^{p\times
r}}\frac{1}{2n}\sum_{k=1}^{r}({\bf y}_{k}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})^{T}({\bf y}_{k}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})+\frac{\gamma}{2}\sum_{k=1}^{r}\left(\sum_{l\in
P_{k}}\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{m\in
N_{k}}\|\mbox{\boldmath$\beta$}_{m}+\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{s\in
Z_{k}}\|\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}\right),$ (3)

where, for a vector $a\in\mathbb{R}^{d}$ define
$\|a\|_{q}=(\sum_{j=1}^{d}a_{j}^{q})^{\frac{1}{q}}$ as the $L_{q}$ norm and
$\gamma>0$ is a tuning parameter to promote similarities between groups of
related coefficients.

The method we propose is a generalization of (3) when the sets $P_{k}$,
$N_{k}$ and $Z_{k}$ are not known _a priori_ and uses an $L_{1}$ penalty for
feature selection. The proposed method simultaneously estimates these unknown
sets, while using the framework of (3) to include these detected relationships
in the estimation. To simultaneously estimate the group structure we propose a
novel penalty based on the minimum function. The penalized objective function
is non-convex but can be re-stated as solving a fixed number of convex
problems. However, this quickly becomes computationally intractable for even
small values of $r$. Unlike the earlier cited work which focuses on structure
due to relationships in the errors, this approach can be easily generalizable
to simultaneous modeling of binary response variables or other non-continuous
values by replacing the least squares loss function with the appropriate
generalized linear model loss function.

In the least squares case we provide rates of convergence, model selection
consistency and post selection inference results which hold for high-
dimensional predictors, $p>>n$. The results hold for any grouping of the
response variables and thus do not depend on the relationships between
responses to be estimated correctly. To derive rates of convergence under
standard conditions we generalize the results of Negahban et al. (2012), which
do not directly apply because of the ridge fusion penalty used in this
framework.

Our proposal builds on methods which simultaneously estimate clusters and
univariate regression models (Witten et al., 2014), multivariate regression
models (Price & Sherwood, 2018), and precision matrices (Price, Molstad &
Sherwood, 2021). All of which propose iterative algorithms that alternate
between estimating the clustering structure, using k-means, and the model of
interest. These approaches only accommodate positive relationships and rely on
k-means, which provides distinct sets of clusters and may be unreliable for
high-dimensional problems. The proposed method allows for more complex
structures in the response variables and provides a global minimum at each
iteration.

Extensive work has been done on exploiting relationships in multivariate
regression to improve efficiency. For a recent survey of the field see Price,
Allenbrand & Sherwood (2021) and citations within. Reduced rank regression
reduces the dimension of the problem by constraining the rank of the
coefficient matrix ${\bf B}$ while minimizing the multivariate least squares
objective function to find a set of latent variables that increase prediction
accuracy (Anderson et al., 1951; Velu & Reinsel, 2013; Chen & Huang, 2016).
Various group lasso penalties have been used to promote structure across
response or a combination of response and predictors. (Kim et al., 2012; Li et
al., 2015) Instead of the proposed approach of simultaneously estimating
relationships and coefficients, many have used a two-stage clustering
approach. In the first stage the clusters are estimated and in the second
stage the structure is incorporated into the estimation of the coefficients.
Examples include hierarchical clustering for tree guided lasso (Kim et al.,
2012), convex clustering in conjoint analysis (Chen et al., 2016), separating
global and response specific features (Xu et al., 2015) and replacing high
dimensional responses or predictors with clusters (Bühlmann et al., 2013; Zhou
et al., 2017).

Other methods in univariate regression, such as the GRACE method, utilize
known relationships, or relationships that are estimated _a priori_ , to
assist in increasing accuracy and select relevant variables (Li & Li, 2008,
2010b). Others such as Zhao & Shojaie (2016) have investigated statistical
inference frameworks for univariate graph constrained models, that also
accommodate error in the estimation of the graph Laplacian.

The key difference between the proposed method and state of the art methods in
the literature, is that our method simultaneously detects relationships
between the responses, regardless of sign, and estimates regression
coefficients without needing to estimate a covariance matrix of the errors. In
Section 2 we introduce the MinPen framework with the least squares loss
function along with theoretical results and the proposed algorithm. In Section
3 we extend the framework to multiple binomial responses with discussion of
how it may generalized to other exponential family based loss functions.
Sections 4 and 5 present simulations to investigate finite sample properties
of the proposed methods. Finally we present examples of our methodology in
applications in genomics and substance abuse overdoses in Section 6.

##  2 Least Squares Model

###  2.1 Method

We consider estimating ${\bf B}^{*}$ in (1) when there may be similarities in
the mean functions of the different responses. Define,

$P(\mbox{\boldmath$\beta$}_{l},\mbox{\boldmath$\beta$}_{k})=\min(\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{k}\|_{2}^{2},\|\mbox{\boldmath$\beta$}_{l}+\mbox{\boldmath$\beta$}_{k}\|_{2}^{2},\|\mbox{\boldmath$\beta$}_{l}\|_{2}^{2}).$
(4)

If
$P(\mbox{\boldmath$\beta$}^{*}_{l},\mbox{\boldmath$\beta$}^{*}_{k})=\|\mbox{\boldmath$\beta$}_{l}^{*}-\mbox{\boldmath$\beta$}^{*}_{k}\|_{2}^{2}$
this implies that the mean functions for response $k$ and $l$ are positively
related. If
$P(\mbox{\boldmath$\beta$}^{*}_{l},\mbox{\boldmath$\beta$}^{*}_{k})=\|\mbox{\boldmath$\beta$}_{l}^{*}+\mbox{\boldmath$\beta$}_{k}^{*}\|_{2}^{2}$
this implies negative relationship between the two mean functions, while
$P(\mbox{\boldmath$\beta$}^{*}_{l},\mbox{\boldmath$\beta$}^{*}_{k})=\|\mbox{\boldmath$\beta$}^{*}_{l}\|_{2}^{2}$
implies little or no relationship between the mean functions. Therefore a
penalty based on the minimum function can be used to simultaneously identify
and leverage relationships between mean functions. Motivated by this, we
propose the minimum penalty elastic net multivariate regression (MinPen)
estimator as the solution to

$\mathop{arg\,min}\displaylimits_{{\bf B}\in\mathbb{R}^{p\times
r}}\frac{1}{2n}\sum_{i=1}^{n}\sum_{k=1}^{r}(y_{ik}-{\bf
x}^{T}_{i}\mbox{\boldmath$\beta$}_{k})^{2}+\delta\sum_{k=1}^{r}\|\mbox{\boldmath$\beta$}_{j}\|_{1}+\frac{\gamma}{2}\sum_{l=1}^{r}\sum_{k=1,k\neq
l}^{r}P(\mbox{\boldmath$\beta$}_{l},\mbox{\boldmath$\beta$}_{k}),$ (5)

where $\delta$ and $\gamma$ are non-negative tuning parameters specified by
the user. We refer to this as an elastic net estimator as it uses the
combination of the lasso and ridge penalties as first introduced by Zou &
Hastie (2005). The lasso penalty, associated with tuning parameter $\delta$,
is used to simultaneously perform variable selection and estimate regression
coefficients. While, the proposed ridge type penalty, associated with
$\gamma$, is a non-convex penalty that simultaneously identifies and exploits
relationships between coefficient vectors of different responses.

Both our theory and algorithms will deal with a vectorized version of the
solution to (5). For a matrix ${\bf A}$, define $\mbox{vec}({\bf A})$ as the
vector formed by stacking the columns of ${\bf A}$ on top of one another.
Define ${\bf Y}=({\bf y}_{1},\ldots,{\bf y}_{n})^{T}\in\mathbb{R}^{n\times
r}$, ${\bf X}=({\bf x}_{1},\ldots,{\bf x}_{n})^{T}\in\mathbb{R}^{n\times r}$,
${\bf E}=(\epsilon_{1},\ldots,\epsilon_{n})^{T}\in\mathbb{R}^{n\times r}$,
$\tilde{{\bf y}}=\mbox{vec}({\bf Y})\in\mathbb{R}^{nr}$, $\tilde{{\bf X}}={\bf
I}_{r}\otimes{\bf X}\in\mathbb{R}^{nr\times pr}$,
$\mbox{\boldmath$\beta$}^{*}=\mbox{vec}({\bf B}^{*})\in\mathbb{R}^{pr}$,
$\tilde{\mbox{\boldmath$\epsilon$}}=\mbox{vec}({\bf E})\in\mathbb{R}^{nr}$ and
$\tilde{\mbox{\boldmath$\Sigma$}}=\mbox{\boldmath$\Sigma$}\otimes{\bf
I}_{n}\in\mathbb{R}^{nr\times nr}$. Then the vectorized version of (1) is

$\tilde{{\bf y}}=\tilde{{\bf
X}}\mbox{\boldmath$\beta$}^{*}+\tilde{\mbox{\boldmath$\epsilon$}}.$ (6)

Define the set

$\operatorname{\mathcal{A}}=\left\\{{\bf A}\in\\{-1,0,1\\}^{r(r-1)p\times
rp}\middle|\,||{\bf
A}\mbox{\boldmath$\beta$}||_{2}^{2}=\sum_{l=1}^{r}\sum_{m=1,l\neq
m}^{r}||\mbox{\boldmath$\beta$}_{l}-d_{lm}\mbox{\boldmath$\beta$}_{m}||_{2}^{2}\mbox{
where }d_{lm}\in\\{-1,0,1\\}\right\\}.$

The estimator in (5) is equivalent to

$(\hat{\mbox{\boldmath$\beta$}},\hat{{\bf
A}})=\underset{\mbox{\boldmath$\beta$}\in\mathbb{R}^{pr}}{\mbox{argmin}}\,\underset{{\bf
A}\in\operatorname{\mathcal{A}}}{\min}\left(\frac{1}{2n}||\tilde{{\bf
Y}}-\tilde{{\bf
X}}\mbox{\boldmath$\beta$}||_{2}^{2}+\delta||\mbox{\boldmath$\beta$}||_{1}+\frac{\gamma}{2}||{\bf
A}\mbox{\boldmath$\beta$}||_{2}^{2}\right).$ (7)

Due to the minimum function the objective function is non-convex, which
typically makes theoretical study of the global minimizer challenging.
However, the global solution to (7) is the pair
$(\hat{\mbox{\boldmath$\beta$}},\hat{{\bf A}})$ that minimizes across the
$|\operatorname{\mathcal{A}}|$ potential convex objective functions. So while
the objective function is nonconvex, the solution
$\hat{\mbox{\boldmath$\beta$}}$ minimizes a convex function. If the matrix
${\bf A}$ is known _a priori_ the estimator of $\mbox{\boldmath$\beta$}^{*}$
is

$\hat{\mbox{\boldmath$\beta$}}({\bf
A})=\mathop{arg\,min}\displaylimits_{\mbox{\boldmath$\beta$}\in\mathbb{R}^{pr}}\frac{1}{2n}\sum_{i=1}^{n}\sum_{k=1}^{r}(y_{i}-{\bf
x}^{T}_{i}\mbox{\boldmath$\beta$}_{k})^{2}+\delta\sum_{k=1}^{r}\|\mbox{\boldmath$\beta$}_{k}\|_{1}+\frac{\gamma}{2n}\mbox{\boldmath$\beta$}^{T}{\bf
A}^{T}{\bf A}\mbox{\boldmath$\beta$}.$ (8)

For our theoretical results we first analyze the estimator
$\hat{\mbox{\boldmath$\beta$}}({\bf A})$. Then, using that the solution is one
of the potential $|\operatorname{\mathcal{A}}|$ solutions, we provide rates of
convergence for the global solution, $\hat{\mbox{\boldmath$\beta$}}$. It may
be desirable to replace $\operatorname{\mathcal{A}}$ with some subset, say
$\operatorname{\mathcal{A}}_{1}\subset\operatorname{\mathcal{A}}$. If it is
known that two responses are positively related then it does not makes sense
to consider the negative or no relationship for those two responses and thus a
suitable $\operatorname{\mathcal{A}}_{1}$ could be used instead of
$\operatorname{\mathcal{A}}$. For this paper we consider the general case, but
the results and algorithm presented here will hold if a smaller set is
considered.

###  2.2 Theoretical Results

Mild or standard conditions were used to derive the rates of convergence.
Before presenting the conditions, we define some notation that is used in our
theorems and conditions. The subspace for the active predictors is defined as
$\mathcal{M}\equiv\\{{\bf a}\in\mathbb{R}^{pr}|a_{j}=0\mbox{ if
}\mbox{\boldmath$\beta$}^{*}_{j}=0\\},$ with cardinality $s=|\mathcal{M}|$.
The parameter space is separated using projections of vectors into orthogonal
complements. Define ${\bf u}_{\mathcal{M}}$ as a projection of a vector ${\bf
u}$ into space $\mathcal{M}$, $\mathcal{M}^{\perp}$ as the orthogonal
complement of $\mathcal{M}$ and $\mathcal{C}\equiv\\{{\bf
a}\in\mathbb{R}^{pr}|\,||{\bf a}_{\mathcal{M}^{\perp}}||_{1}\leq 3||{\bf
a}_{\mathcal{M}}||_{1}\\}$.

######  Condition 1.

Define ${\bf x}_{j}\in\mathbb{R}^{p}$ to be the $j$th column vector of ${\bf
X}$. For all $j\in\\{1,\ldots,p\\}$ $\frac{||{\bf x}_{j}||_{2}^{2}}{n}\leq 1$.

######  Condition 2.

There exists positive constants $\kappa_{l}$ and $\kappa_{u}$ such that for
all $\mbox{\boldmath$\theta$}\in\mathcal{C}$ that
$\kappa_{l}\leq\mbox{\boldmath$\theta$}^{T}\frac{1}{n}\tilde{{\bf
X}}^{T}\tilde{{\bf X}}\mbox{\boldmath$\theta$}\leq\kappa_{u}.$

######  Condition 3.

The error vector $\mbox{\boldmath$\epsilon$}_{j}$ has a mean of zero and sub-
Gaussian tails for all $j\in\\{1,\ldots,r\\}$. That is, for all
$j\in\\{1,\ldots,r\\}$ there exists a constant $\sigma_{j}$ such that for any
${\bf a}\in\mathbb{R}^{n}$, with $||{\bf a}||_{2}=1$,

$P(|\langle\mbox{\boldmath$\epsilon$}_{j},{\bf a}\rangle|>t)\leq
2\mbox{exp}\left(-\frac{t^{2}}{2\sigma_{j}^{2}}\right).$

In addition, there exists a positive constant $\sigma$ such that
$\sigma\geq\underset{j}{\max}\,\sigma_{j}$.

######  Condition 4.

There is at least one non-zero entry in $\mbox{\boldmath$\beta$}^{*}$.

Conditions 1-3 are commonly made for penalized estimators, see Negahban et al.
(2012) and the citations within. The framework in Negahban et al. (2012) does
not directly apply to (7), because they assume the penalty is a norm which is
not the case for our penalty, even for ${\bf A}$ fixed, because of the
inclusion of the Ridge penalty. However, under these conditions we can extend
the approach of Negahban et al. (2012) and prove rates of convergence under
milder conditions then some similar methods (Price & Sherwood, 2018; Li & Li,
2010b). In our results the upper bound of the tuning parameter $\gamma$
depends on the inverse of $||\mbox{\boldmath$\beta$}^{*}||_{\infty}$, which is
well defined under Condition 4.

######  Theorem 1.

Assume Conditions 1-4 hold and
$\gamma<\delta/(16r||\mbox{\boldmath$\beta$}^{*}||_{\infty})$ then for any
${\bf A}\in\operatorname{\mathcal{A}}$

$P\left[||\hat{\mbox{\boldmath$\beta$}}({\bf
A})-\mbox{\boldmath$\beta$}^{*}||_{2}^{2}\leq\frac{9\delta^{2}}{4\kappa_{l}^{2}}s\right]\geq
1-2\mbox{exp}\left[\frac{-n\delta^{2}}{32\sigma^{2}}+\log(rp)\right].$

Proof of Theorem 1, and all other theoretical results, are provided in the
supplemental material. Theorem 1 can be used to derive convergence rates for
$\hat{\mbox{\boldmath$\beta$}}({\bf A})$ for a set sequence of $\delta$.

######  Corollary 1.

Assume Conditions 1-3 hold,
$\gamma<\delta/(16r||\mbox{\boldmath$\beta$}^{*}||_{\infty})$ and
$\delta=\sqrt{\frac{64\sigma^{2}\log(rp)}{n}}$ then for any ${\bf
A}\in\operatorname{\mathcal{A}}$

$P\left[||\hat{\mbox{\boldmath$\beta$}}({\bf
A})-\mbox{\boldmath$\beta$}^{*}||_{2}^{2}\leq
144\frac{\sigma^{2}\log(rp)}{n\kappa_{l}^{2}}s\right]\geq
1-2\mbox{exp}[-\log(rp)].$

Corollary 1 provides that for fixed ${\bf A}$ and $rp\rightarrow\infty$ that
the estimator $\hat{\mbox{\boldmath$\beta$}}({\bf A})$ achieves the same rate
of converge as the lasso estimator (Negahban et al., 2012) and elastic net
(Hebiri & van de Geer, 2011). Using Corollary 1 and noting that
$\hat{\mbox{\boldmath$\beta$}}$ is equivalent to a
$\hat{\mbox{\boldmath$\beta$}}(\hat{{\bf A}})$ we can derive a rate of
convergence for $\hat{\mbox{\boldmath$\beta$}}$, the global minimum of a
nonconvex objective function.

######  Corollary 2.

Assume Conditions 1-3 hold,
$\gamma<\delta/(16r||\mbox{\boldmath$\beta$}^{*}||_{\infty})$ and
$\delta=\sqrt{\frac{64\sigma^{2}\log(rp)}{n}}$ then

$P\left[||\hat{\mbox{\boldmath$\beta$}}-\mbox{\boldmath$\beta$}^{*}||_{2}^{2}\leq
144\frac{\sigma^{2}\log(rp)}{n\kappa_{l}^{2}}s\right]\geq
1-2\mbox{exp}[-\log(rp)+\log(|\operatorname{\mathcal{A}}|)].$

###  2.3 Model Selection Consistency

Let $\dot{{\bf X}}_{(1)}\in\mathbb{R}^{nr\times s}$ be the submatrix of
$\tilde{{\bf X}}$ consisting of the active predictors of the regression
defined in (6), and $\dot{{\bf X}}_{(2)}\in\mathbb{R}^{nr\times pr-s}$ be the
submatrix of the remaining predictors. Define $\dot{{\bf X}}=(\dot{{\bf
X}}_{(1)},\dot{{\bf X}}_{(2)})$ and let $\dot{\mbox{\boldmath$\beta$}}^{*}$
and $\dot{\mbox{\boldmath$\beta$}}({\bf A})$ be the similarly re-arranged
versions of $\mbox{\boldmath$\beta$}^{*}$ and
$\hat{\mbox{\boldmath$\beta$}}({\bf A})$, respectively. In addition, let
$\dot{\mbox{\boldmath$\beta$}}^{*}_{(1)}$ be the first $s$ entries of
$\dot{\mbox{\boldmath$\beta$}}^{*}$. Note, that we still have

$\tilde{y}=\dot{{\bf
X}}\dot{\mbox{\boldmath$\beta$}}^{*}+\tilde{\mbox{\boldmath$\epsilon$}}=\dot{{\bf
X}}_{(1)}\dot{\mbox{\boldmath$\beta$}}^{*}_{(1)}+\tilde{\mbox{\boldmath$\epsilon$}}.$

Thus model selection consistency results found for
$\dot{\mbox{\boldmath$\beta$}}({\bf A})$, also hold for
$\hat{\mbox{\boldmath$\beta$}}({\bf A})$.

Define $\dot{{\bf A}}$ such that it relates to $\dot{{\bf X}}$ the same way
${\bf A}$ corresponds to ${\bf X}$. The signed graph Laplacian of $\dot{{\bf
A}}\in\operatorname{\mathcal{A}}$ is

${\bf L}=\dot{{\bf A}}^{T}\dot{{\bf A}}=\left(\begin{array}[]{cc}{\bf
L}_{11}&{\bf L}_{12}\\\ {\bf L}_{21}&{\bf L}_{22}\end{array}\right),$

where ${\bf L}_{11}$ and ${\bf L}_{22}$ are the signed graph Laplacians of the
active and inactive sets respectfully. Finally, define ${\bf
C}_{11}=\frac{1}{n}\dot{{\bf X}}_{(1)}^{T}\dot{{\bf X}}_{(1)}$, ${\bf
C}_{22}=\frac{1}{n}\dot{{\bf X}}_{(2)}^{T}\dot{{\bf X}}_{(2)}$, and ${\bf
C}_{21}={\bf C}^{T}_{12}=\frac{1}{n}\dot{{\bf X}}_{(2)}^{T}\dot{{\bf
X}}_{(1)}$.

Zhao & Yu (2006) first presented conditions for model selection consistency
for lasso. While, Li & Li (2010b) present model selection results for a lasso-
type estimator that includes a penalty for a graph structure in the
predictors. The estimator in that paper is similar to the one in this paper,
but with ${\bf A}$ fixed and a univariate response. The following conditions
are generalizations of conditions presented in those papers. For vectors ${\bf
a},{\bf b}\in\mathbb{R}^{d}$, let $\mbox{sign}({\bf a})$ be a vector of the
signs of the entries of ${\bf a}$, $|{\bf a}|$ be the vector formed by taking
the absolute values of each entry of ${\bf a}$, $\min({\bf a})$ is the
smallest entry of ${\bf a}$ and ${\bf a}<{\bf b}$ is true if and only if
$a_{i}<b_{i}$ for all $i\in\\{1,\ldots,d\\}$.

######  Condition 5.

Assume $\mbox{\boldmath$\epsilon$}_{i}\sim N_{r}(\mathbf{0}_{r},\sigma^{2}{\bf
D}_{0})$, are $i.i.d$ and ${\bf D}_{0}\in\mathbb{R}^{r\times r}$ is a diagonal
matrix where diagonal values are on the support $(0,1]$ and thus the maximum
variance of the errors is $\sigma^{2}$.

######  Condition 6.

The minimum eigenvalue of ${\bf C}_{11}$, $k_{1}$, is positive for all $n$.

######  Condition 7.

There exists a fixed $\eta$, that does not change with $n$, such that

$\left|\left({\bf C}_{21}+\frac{\gamma}{n}{\bf L}_{21}\right)\left({\bf
C}_{11}+\frac{\gamma}{n}{\bf
L}_{11}\right)^{-1}\left\\{\operatorname{sign}[\dot{\mbox{\boldmath$\beta$}}^{*}_{(1)}]+\frac{\gamma}{\delta}{\bf
L}_{11}\dot{\mbox{\boldmath$\beta$}}^{*}_{(1)}\right\\}-\frac{\gamma}{\delta}{\bf
L}_{21}\dot{\mbox{\boldmath$\beta$}}^{*}_{(1)}\right|\leq\vec{1}-\eta.$

Condition 7 is multivariate response version of the graph constrained
irrepresentable condition from Li & Li (2010b), which is an extension of the
original irrepresentable condition in Zhao & Yu (2006). Using these
conditions, Lemma 1 from Li & Li (2010a) and additional conditions on the
tuning parameters, we can prove the estimator is model selection consistent.

######  Theorem 2.

Assume condition 1, 5, 6, and 7 hold. Define
${\color[rgb]{1,0,0}\bf\rho}=min\left[|({\bf C}_{11}+\gamma/n{\bf
L}_{11})^{-1}({\bf C}_{11}\mbox{\boldmath$\beta$}_{(1)})|\right]$.If

$\frac{\delta^{2}}{\log(rp-s)(n+\gamma^{2}r/k_{1})}\rightarrow\infty,$

and

$\frac{1}{\rho}\left\\{\sqrt{\frac{\log(s)}{nk_{1}}}+\delta\left\|({\bf
C}_{11}+\gamma/n{\bf
L}_{11})^{-1}\operatorname{sign}[\dot{\mbox{\boldmath$\beta$}}_{(1)}^{*}]\right\|_{\infty}\right\\}\rightarrow
0,$

then for any ${\bf A}$,

$\lim_{n\rightarrow\infty}P\left\\{\operatorname{sign}[\dot{\mbox{\boldmath$\beta$}}({\bf
A})]=\operatorname{sign}(\dot{\mbox{\boldmath$\beta$}}^{*})\right\\}=1.$

Thus the estimator is is strongly sign consistent and therefore model
selection consistency holds, see Zhao & Yu (2006) for a more detailed
discussion of this relationship.

###  2.4 Post Selection Inference

This subsection provides details how to construct confidence intervals for
non-zero coefficients of $\hat{\mbox{\boldmath$\beta$}}(\hat{{\bf A}})$ using
the framework provided in Lee et al. (2016) for lasso estimators.
Interpretation of models depends on the variables in the model and we restrict
our analysis to inference conditional on the selected model. A consequence of
this is $\mbox{\boldmath$\beta$}^{*}$ will change depending on the model. Let
$Q\subseteq\\{1,\ldots,pr\\}$ be a subset of selected predictor-response
combinations. For simplicity of notation we assume that each response has at
least one predictor selected. Define $\tilde{{\bf
X}}_{Q}\in\mathbb{R}^{nr\times|Q|}$ to be the matrix made from the columns of
$\tilde{{\bf X}}$ that belong to $Q$. Define
$\mbox{\boldmath$\mu$}=E(\tilde{{\bf y}})\in\mathbb{R}^{nr}$, $\tilde{{\bf
X}}_{Q}^{+}=\left(\tilde{{\bf X}}^{T}_{Q}\tilde{{\bf
X}}_{Q}\right)^{-1}\tilde{{\bf X}}_{Q}^{T}\in\mathbb{R}^{|Q|\times nr}$ and

$\mbox{\boldmath$\beta$}^{*}_{Q}=\underset{\mbox{\boldmath$\beta$}\in\mathbb{R}^{|Q|}}{\mbox{arg
min}}E\left|\left|\tilde{{\bf y}}-\tilde{{\bf
X}}_{Q}\mbox{\boldmath$\beta$}\right|\right|^{2}_{2}=\tilde{{\bf
X}}_{Q}^{+}\mbox{\boldmath$\mu$}.$ (9)

Where
$\mbox{\boldmath$\beta$}^{*}_{Q}=(\beta^{*}_{Q,1},\ldots,\beta^{*}_{Q,|Q|})^{T}\in\mathbb{R}^{|Q|}$,
is the true vector of the coefficients conditional on using only the
covariates in $Q$. For purposes of post selection inference, we assume that
$\tilde{{\bf y}}$ follows a normal distribution of $\tilde{{\bf y}}\sim
N(\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\epsilon$}}\otimes{\bf
I}_{n})$.

Define the vector $\hat{{\bf
s}}=(\hat{s}_{1},\ldots,\hat{s}_{pr})^{T}\in\mathbb{R}^{pr}$, where
$\hat{s}_{j}=\mbox{sign}(\hat{\beta}_{j})$ if $\hat{\beta}_{j}\neq 0$ and
$\hat{s}_{j}\in[-1,1]$ if $\hat{\beta}_{j}=0$. As defined in Tibshirani (2013)
for the lasso estimator, for technical convenience we consider the
equicorrelation set

$\hat{Q}\equiv\\{j\in\\{1,\ldots,pr\\}\mid|\hat{s}_{j}|=1\\}.$ (10)

The index of nonzero coefficients of $\hat{\mbox{\boldmath$\beta$}}$ is a
subset of $\hat{Q}$ and typically these sets are the same.

Let $\hat{Q}$ represent the selected model and let $\hat{{\bf
A}}\in\mathbb{R}^{pr(r-1)/2\times pr}$ be the graph selected from (7).
Consider a specific index $j\in Q$, the goal of this work is to define an
interval $C_{q}^{{\bf A},Q}$ for a given level $\alpha$, such that

$P(\beta_{Q,j}^{*}\in C_{j}^{{\bf A},Q}|\hat{Q}=Q,\hat{{\bf A}}={\bf A})\geq
1-\alpha.$ (11)

See Lee et al. (2016) for a justification of coverage of post-selection
intervals being conditional on the selected model and a review of other
approaches to post-selection inference. To derive this interval we consider an
arbitrary vector of $\mbox{\boldmath$\eta$}_{Q}\in\mathbb{R}^{nr}$ examine the
conditional distribution of

$\mbox{\boldmath$\eta$}^{T}_{Q}\tilde{{\bf y}}|(\hat{Q}=Q,\hat{{\bf A}}={\bf
A}).$ (12)

Note that this definition does include all covariates that have non-zero
coefficients, but could potentially include a non-zero coefficient as
$\hat{s}_{i}\in[-1,1]$. For $\hat{\mbox{\boldmath$\beta$}}$ and $\hat{s}$ to
solve (8), conditional on $\hat{{\bf A}}={\bf A}$, the following sufficient
and necessary KKT conditions must be satisfied,

$\left(\frac{1}{n}\tilde{{\bf X}}^{T}\tilde{{\bf X}}+\gamma{\bf A}^{T}{\bf
A}\right)\hat{\mbox{\boldmath$\beta$}}({\bf A})-\frac{1}{n}\tilde{{\bf
X}}^{T}\tilde{{\bf y}}+\delta\hat{s}=\mathbf{0}_{pr},$ (13)
$\hat{s}_{j}=\mbox{sign}(\hat{\beta}_{j})\mbox{ if }\hat{\beta}_{j}\neq 0,$
(14)

and

$\hat{s}_{j}\in[-1,1]\mbox{ if }\hat{\beta}_{j}=0.$ (15)

For a given ${\bf A}$, $\hat{\mbox{\boldmath$\beta$}}({\bf A})$ is a unique
minimizer of (8) and confidence intervals are constructed using these KKT
conditions.

For the following results we need the OLS solution to be well defined and thus
need $\tilde{{\bf X}}_{Q}^{T}\tilde{{\bf X}}_{Q}$ to be positive definite.

######  Condition 8.

The matrix $\tilde{{\bf X}}_{Q}^{T}\tilde{{\bf X}}_{Q}$ is positive definite.

Define ${\bf e}_{j}$ as the vector with zeros everywhere, but with a one in
the $j$th position. The following theorem presents how to construct confidence
intervals conditional on $\hat{Q}$ and $\hat{{\bf A}}$.

######  Theorem 3.

Let $F^{[a,b]}_{\mu,\sigma^{2}}(x)$ be the CDF of a random normal variable
with mean $\mu$, variance $\sigma^{2}$, truncated to the interval $[a,b]$. Let
$\mbox{\boldmath$\eta$}=\tilde{{\bf X}}(\tilde{{\bf X}}^{T}\tilde{{\bf
X}})^{-1}{\bf e}_{j}$ and define $L$ and $U$ such that

$F_{L,\mbox{\boldmath$\eta$}^{T}\tilde{\mbox{\boldmath$\Sigma$}}\mbox{\boldmath$\eta$}}^{[V_{{\bf
s}_{Q},{\bf A}}^{-},V_{{\bf s}_{Q},{\bf
A}}^{+}]}(\mbox{\boldmath$\eta$}^{T}\tilde{{\bf y}})=1-\frac{\alpha}{2}\mbox{
and
}F_{U,\mbox{\boldmath$\eta$}^{T}\tilde{\mbox{\boldmath$\Sigma$}}\mbox{\boldmath$\eta$}}^{[V_{{\bf
s}_{Q},{\bf A}}^{-},V_{{\bf s}_{Q},{\bf
A}}^{+}]}(\mbox{\boldmath$\eta$}^{T}\tilde{{\bf y}})=\frac{\alpha}{2}.$

Where $V_{{\bf s}_{Q},{\bf A}}^{-}$ and $V_{{\bf s}_{Q},{\bf A}}^{+}$ are
constants that depend on ${\bf s}_{Q}\in\mathbb{R}^{|Q|}$, ${\bf A}$, $Q$,
$\tilde{\mbox{\boldmath$\Sigma$}}$, $\tilde{{\bf X}}_{Q}$ and ${\bf y}$. If
Condition 8 holds then

$P\left(\beta_{Q,j}^{*}\in[L,U]|\hat{Q}=Q,\hat{{\bf s}}_{\hat{Q}}={\bf
s}_{Q},\hat{{\bf A}}={\bf A}\right)=1-\alpha.$

Theorem 3 provides a way to construct confidence intervals for a least squares
model that accounts for the fact that variable selection was first done by
satisfying the KKT conditions of the estimator
$\hat{\mbox{\boldmath$\beta$}}(\hat{{\bf A}})$. The result uses the framework
proposed in Lee et al. (2016), which provides a more detailed discussion of
the issue of post-selection inference and the intuition behind the interval
construction. Theorem 3 assumes that $\tilde{\mbox{\boldmath$\Sigma$}}$ is
known, but usually in practice it will need to be estimated. If $p<n$, then
the covariance matrix of the residuals from a saturated multivariate
regression model could be used. An alternative approach would be to make some
simplifying assumption to make it easier to estimate $\Sigma$, for instance
assume it is diagonal.

###  2.5 Algorithm

To obtain the estimates for the estimator defined in (5) we propose using the
formulation of the estimator that is similar to the optimization defined in
(3). Define the following sets

$\begin{split}P_{l}=&\left\\{m:\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{m}\|_{2}^{2}\leq\|\mbox{\boldmath$\beta$}_{l}+\mbox{\boldmath$\beta$}_{m}\|_{2}^{2},\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{m}\|_{2}^{2}<\|\mbox{\boldmath$\beta$}_{l}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\}\\\
N_{l}=&\left\\{m:\|\mbox{\boldmath$\beta$}_{l}+\mbox{\boldmath$\beta$}_{m}\|_{2}^{2}<\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{m}\|_{2}^{2},\|\mbox{\boldmath$\beta$}_{l}+\mbox{\boldmath$\beta$}_{m}\|_{2}^{2}<\|\mbox{\boldmath$\beta$}_{l}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\}\\\
Z_{l}=&\left\\{m:\|\mbox{\boldmath$\beta$}_{l}\|_{2}^{2}\leq\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{m}\|_{2}^{2},\|\mbox{\boldmath$\beta$}_{l}\|_{2}^{2}\leq\|\mbox{\boldmath$\beta$}_{l}+\mbox{\boldmath$\beta$}_{m}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\}\\\
\end{split}$

such that $P_{l}\cup N_{l}\cup Z_{l}=\\{1,\ldots,r\\}\setminus\\{l\\}$,
$P_{l}\cap N_{l}=\emptyset$, $P_{l}\cap Z_{l}=\emptyset$, and $Z_{l}\cap
N_{l}=\emptyset$ $\forall l\in\\{1,\ldots,r\\}$. Solving (5) is equivalent to
solving

$\begin{split}\mathop{arg\,min}\displaylimits_{{\bf B}\in\mathbb{R}^{p\times
r},P_{1},N_{1},Z_{1},\ldots,P_{r},N_{r},Z_{r}}&\frac{1}{2n}\sum_{k=1}^{r}\sum_{i=1}^{n}({\bf
y}_{i}-{\bf x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})^{T}({\bf y}_{i}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})+\delta\sum_{k=1}^{r}\|\mbox{\boldmath$\beta$}_{j}\|_{1}\\\
+&\frac{\gamma}{2}\sum_{k=1}^{r}\left(\sum_{l\in
P_{k}}\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{m\in
N_{k}}\|\mbox{\boldmath$\beta$}_{m}+\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{s\in
Z_{k}}\|\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}\right).\end{split}$ (16)

By using this formulation we are able to propose an iterative algorithm that
iterates between estimating the sets $P_{l},N_{l},$ and $Z_{l}$, for each
$l\in\\{1,\ldots,r\\}$, with the regression coefficients fixed, and estimating
the regression coefficients with the sets fixed. To initialize the algorithm
estimates of ${\bf B}$ or $P_{l},N_{l},$ and $Z_{l}$, $l=1,\ldots,r$ are
needed. We define initial values for the coefficients of the $c$th response
are set by

$\hat{\mbox{\boldmath$\beta$}}^{(0)}_{c}=\mathop{arg\,min}\displaylimits_{\mbox{\boldmath$\beta$}_{c}\in\mathbb{R}^{p}}\frac{1}{2n}\sum_{i=1}^{n}(y_{ic}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{c})^{2}+\delta\|\mbox{\boldmath$\beta$}_{c}\|_{1}+\gamma\|\mbox{\boldmath$\beta$}_{c}\|_{2}^{2},$
(17)

and define $\hat{{\bf
B}}^{(w)}=\left(\hat{\mbox{\boldmath$\beta$}}_{1}^{(w)},\ldots,\hat{\mbox{\boldmath$\beta$}}^{(w)}_{r}\right)$
to be the $w$th iterative estimate of ${\bf B}^{*}$. Note we will also use
$\hat{\mbox{\boldmath$\beta$}}^{(w)}$ to represent the vectorized version of
$\hat{{\bf B}}^{(w)}$, consistent with previous notation. Given a fixed
$(\delta,\gamma)$ we propose the following algorithm as a solution to (16).

  1. 1.

Initialize
$\hat{\mbox{\boldmath$\beta$}}_{1}^{(0)},\ldots,\hat{\mbox{\boldmath$\beta$}}_{r}^{(0)}$
as defined in (17).

  2. 2.

For the $w$th iteration, where $w>0$, repeat the following steps until the
estimated sets $P_{l},N_{l}$ and $Z_{l}$ do not change from iterate $w-1$ to
iterate $w$ for all $l\in\\{1,\ldots,r\\}$.

    1. (a)

Holding $\hat{{\bf B}}^{(w-1)}$ fixed obtain estimates for
$\hat{P}^{(w)}_{k},\hat{N}^{(w)}_{k},\hat{Z}^{(w)}_{k}$ for all
$k\in\\{1,\ldots,r\\}$ by solving the optimization

$\begin{split}\underset{P_{k},N_{k},Z_{k}\forall
k\in\\{1,\ldots,r\\}}{\mbox{minimize}}&\sum_{k=1}^{r}\left(\sum_{l\in
P_{k}}\|\hat{\mbox{\boldmath$\beta$}}^{(w-1)}_{l}-\hat{\mbox{\boldmath$\beta$}}^{(w-1)}_{k}\|_{2}^{2}+\sum_{m\in
N_{k}}\|\hat{\mbox{\boldmath$\beta$}}^{(w-1)}_{m}+\hat{\mbox{\boldmath$\beta$}}^{(w-1)}_{k}\|_{2}^{2}+\sum_{s\in
Z_{k}}\|\hat{\mbox{\boldmath$\beta$}}^{(w-1)}_{k}\|_{2}^{2}\right).\end{split}$
(18)

    2. (b)

Holding $\hat{P}^{(w)}_{k},\hat{N}^{(w)}_{k},\hat{Z}^{(w)}_{k}$ for all
$k\in\\{1,\ldots,r\\}$ fixed obtain the estimate $\hat{{\bf B}}^{(w)}$ by
solving the optimization

$\begin{split}\mathop{arg\,min}\displaylimits_{{\bf B}\in\mathbb{R}^{p\times
r}}&\frac{1}{2n}\sum_{k=1}^{r}\sum_{i=1}^{n}({\bf y}_{i}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})^{T}({\bf y}_{i}-{\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}_{k})+\delta\sum_{k=1}^{r}\|\mbox{\boldmath$\beta$}_{j}\|_{1}\\\
+&\frac{\gamma}{2}\sum_{k=1}^{r}\left(\sum_{l\in\hat{P}^{(w)}_{k}}\|\mbox{\boldmath$\beta$}_{l}-\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{m\in\hat{N}^{(w)}_{k}}\|\mbox{\boldmath$\beta$}_{m}+\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}+\sum_{s\in\hat{Z}^{(w)}_{k}}\|\mbox{\boldmath$\beta$}_{k}\|_{2}^{2}\right).\end{split}$
(19)

The update steps, shown in (18) and (19) respectively, break the non-convex
problem into two problems that can be solved directly with well studied
solutions. The optimization in (18) can actually be solved directly in a
single pass, where the set assignments to the $w+1$ iteration for
$l=1,\ldots,r$ are defined as

$\begin{split}\hat{P}^{(w+1)}_{l}=&\left\\{m:\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}-\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2}\leq\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}+\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2},\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}-\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2}<\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\},\\\
\hat{N}^{(w+1)}_{l}=&\left\\{m:\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}+\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2}<\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}-\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2},\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}+\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2}<\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\},\\\
\hat{Z}^{(w+1)}_{l}=&\left\\{m:\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}\|_{2}^{2}\leq\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}-\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2},\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}\|_{2}^{2}\leq\|\hat{\mbox{\boldmath$\beta$}}^{(w)}_{l}+\hat{\mbox{\boldmath$\beta$}}^{(w)}_{m}\|_{2}^{2},m\in\\{1,\ldots,r\\}\setminus\\{l\\}\right\\}.\\\
\end{split}$

Note that this update requires calculating $r^{2}$ terms and performing
$r(r-1)$ comparisons, but the sets are deterministic and no iterates are
needed.

The optimization in (19) can be solved using a gradient descent algorithm
where each $\beta_{mk}$ is solved for iteratively with all other regression
coefficients held fixed. Let $S(x,a)=\mbox{sign}(x)(|x|-a)_{+}$ be the soft-
thresholding operator. To solve (19) we use a coordinate descent algorithm
where the update for the $m$th predictor for the $k$th response,
$\bar{\beta}_{mk}$, is

$\bar{\beta}_{mk}=\frac{S\left(n^{-1}\\{\sum_{i=1}^{n}x_{im}y_{ik}-\sum_{j=1,j\neq
m}^{p}(\sum_{i=1}^{n}x_{im}x_{ij})\bar{\beta}_{jk}+\gamma
H_{mk}\\},\delta/2\right)}{n^{-1}\left(\sum_{i=1}^{n}x_{im}^{2}+\gamma\left\\{(r-1)+\sum_{g=1}^{r}I\left[k\in\hat{P}_{g}^{(w)}\right]+\sum_{v=1}^{r}I\left[k\in\hat{N}_{v}^{(w)}\right]\right\\}\right)},$

where

$H_{mk}=\sum_{l\in\hat{P}_{k}^{(w)}}\bar{\beta}_{ml}-\sum_{h\in\hat{N}_{k}^{(w)}}\bar{\beta}_{mh}+\sum_{g=1}^{r}\bar{\beta}_{mg}I\left[k\in\hat{P}_{g}^{(w)}\right]-\sum_{v=1}^{r}\bar{\beta}_{mv}I\left[k\in\hat{N}_{v}^{(w)}\right],$

and $\bar{\beta}_{jv}$ are the current iterates of the coordinate descent
algorithm. The algorithm iterates through all $m\in\\{1,\ldots,p\\}$ and
$k\in\\{1,\ldots,p\\}$ in a similar fashion to other coordinate descent
algorithms and converges with the change in the estimates across iterations is
small. We propose selecting $\delta$ and $\gamma$ using $k$-fold cross
validation minimizing the validation residual sum of squares. The proposed
coordinate descent update is similar to GRACE, but differs due to structured
use of the minimum penalty across multivariate responses.

While this algorithm is related to the iterative algorithms proposed by Witten
et al. (2014), Price & Sherwood (2018), and Price, Molstad & Sherwood (2021),
there are some differences. First, the penalties investigated by those authors
are considered _cluster fusion penalties_. Second, the algorithms proposed in
those papers are two stage procedure that require a solution to the well
studied $k$-means problem, which can be unreliable and computationally
burdensome in high-dimensional settings. The proposed algorithm provides an
efficient algorithm that is guaranteed to find the optimal sets of the
responses given the coefficients. An alternative approach would be an
exhaustive search of all possible ${\bf A}\in\operatorname{\mathcal{A}}$ which
would provide a global minimizer. However, this becomes computational
intractable for medium or large $r$ as there are $3^{r(r-1)}$ combinations
that must be searched before tuning parameter selection is considered. The
iterative method we propose is not constrained in this way but is not
guaranteed to achieve a global solution to the non-convex optimization
problem. The proposed algorithm can be easily adjusted to accommodate if other
penalty functions, besides the proposed ridge fusion, are used within the
minimum function.

##  3 Extension to Multiple Binomial Responses

###  3.1 Method

This section details how to extend the proposed method to the setting where
conditional on the predictors each response follows a binomial distribution.
The focus will be on the binomial setting, but the method presented in this
section can be generalized to multiple responses from the same exponential
family with different mean functions. Define
$\mbox{\boldmath$\theta$}_{k}=(\alpha_{k},\mbox{\boldmath$\omega$}_{k}^{T})$,
where $\alpha_{k}$ is a response specific intercept, and ${\bf u}_{i}=(1,{\bf
x}_{i})$ be the $p+1$ dimensional vector of covariates for the $i$th
observation. Let $y_{ik}$ be a realization of the random variable

$Y_{ik}\sim\operatorname{Bin}(n_{ik},\pi_{ik}),\mbox{ where
}\pi_{ik}=\frac{\exp({\bf u}_{i}^{T}\mbox{\boldmath$\theta$}_{k})}{1+\exp({\bf
u}_{i}^{T}\mbox{\boldmath$\theta$}_{k})}.$

We define the penalized likelihood for the minimum penalized model as

$\begin{split}\sum_{k=1}^{r}\sum_{i=1}^{n}&y_{ik}{\bf
u}_{i}^{T}\mbox{\boldmath$\theta$}_{k}-n_{ik}\log\left[1+\exp({\bf
u}_{i}^{T}\mbox{\boldmath$\theta$}_{k})\right]+\delta\sum_{k=1}^{r}\|\mbox{\boldmath$\omega$}_{k}\|_{1}\\\
+&\frac{\gamma}{2}\sum_{l=1}^{r}\sum_{m=1,l\neq
m}^{r}\min\left(\|\mbox{\boldmath$\theta$}_{l}-\mbox{\boldmath$\theta$}_{m}\|_{2}^{2},\|\mbox{\boldmath$\theta$}_{l}+\mbox{\boldmath$\theta$}_{m}\|_{2}^{2},\|\mbox{\boldmath$\theta$}_{l}\|_{2}^{2}\right).\end{split}$
(20)

The lasso penalty does not include the intercept, because we assume that an
intercept is part of the model for each response. However, the intercept is
part of the penalty in the minimum function because we wish to group responses
based on the relationships between fitted values.

###  3.2 Algorithm

We propose solving (20) by approximating it with a penalized quadratic
function similar to the glmnet algorithm (Friedman et al., 2010). Define,

$g(\pi_{ik})=\log\left(\frac{\pi_{ik}}{1-\pi_{ik}}\right)={\bf
u}_{i}^{T}\mbox{\boldmath$\theta$}_{k}.$ (21)

To implement our penalized quadratic approximation we define the following,

$\displaystyle z_{ik}$ $\displaystyle=$ $\displaystyle
g(y_{ik})=g(\pi_{ik})+\frac{y_{ik}-\pi_{ik}}{\pi_{ik}(1-\pi_{ik})},$ (22)
$\displaystyle w_{ik}$ $\displaystyle=$ $\displaystyle\pi_{ik}(1-\pi_{ik}),$
(23) $\displaystyle-l_{k}(\mbox{\boldmath$\theta$}_{k})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}w_{ik}(z_{ik}-{\bf
u}_{i}^{T}\mbox{\boldmath$\theta$}_{k})^{2}.$ (24)

Note that $z_{ik}$ is the first order Taylor approximation of $g(y_{ik})$, and
that $w_{ik}$ is the conditional variance of $z_{ik}$ given ${\bf u}_{i}$.
Define ${\bf Z}_{k}=(z_{1k},\ldots,z_{nk})^{T}\in\mathbb{R}^{n}$ and ${\bf
W}=(w_{1k},\ldots,w_{nk})^{T}\in\mathbb{R}^{n}$.

Let $P_{l},N_{l}$ and $Z_{l}$ be defined as in Section 2.1, but now with
respect to $\mbox{\boldmath$\theta$}_{l}$ instead of
$\mbox{\boldmath$\beta$}_{l}$. We propose an iterative algorithm similar to
that proposed in Section 2.5 to minimize (20).

Initial estimates of $\mbox{\boldmath$\theta$}_{l}$ or $P_{l},N_{l},$ and
$Z_{l}$ are needed for all $l\in\\{1,\ldots,r\\}$. Similar to the method
proposed in Section 2.5 we propose initializing the regression coefficients
for each of the $r$ responses separately using the elastic net estimator,
which we define as $\hat{\mbox{\boldmath$\theta$}}_{k}^{(0)}$. Given a fixed
$(\delta,\gamma)$ we propose the following algorithm as a solution to (24).

  1. 1.

Initialize
$\hat{\mbox{\boldmath$\theta$}}_{1}^{(0)},\ldots,\hat{\mbox{\boldmath$\theta$}}_{r}^{(0)}$.

  2. 2.

For the $w$th iteration where $w>0$, repeat steps until the estimated sets
$P_{l}$, $N_{l}$ and $Z_{l}$ do not change from iterate $w-1$ to iterate $w$:

    1. (a)

Holding $\hat{\mbox{\boldmath$\theta$}}_{k}^{(w-1)}$ fixed for all
$k\in\\{1,\ldots,r\\}$, obtain $\hat{P}^{(w)}_{l}$, $\hat{N}^{(w)}_{l}$ and
$\hat{Z}^{(w)}_{l}$ for all $l\in\\{1,\ldots,r\\}$ by solving the optimization

$\begin{split}\underset{P_{1},N_{1},Z_{1},...,P_{r},N_{r},Z_{r}}{\mbox{minimize}}&\sum_{k=1}^{r}\left(\sum_{l\in
P_{k}}\|\hat{\mbox{\boldmath$\theta$}}^{(w-1)}_{l}-\hat{\mbox{\boldmath$\theta$}}^{(w-1)}_{k}\|_{2}^{2}+\sum_{m\in
N_{k}}\|\hat{\mbox{\boldmath$\theta$}}^{(w-1)}_{m}+\hat{\mbox{\boldmath$\theta$}}^{(w-1)}_{k}\|_{2}^{2}+\sum_{s\in
Z_{k}}\|\hat{\mbox{\boldmath$\theta$}}^{(w-1)}_{k}\|_{2}^{2}\right).\end{split}$
(25)

    2. (b)

Holding $\hat{P}^{(w)}_{k},\hat{N}^{(w)}_{k},\hat{Z}^{(w)}_{k}$ fixed for all
$k\in\\{1,\ldots,r\\}$ obtain the estimates of
$\hat{\mbox{\boldmath$\theta$}}_{k}^{(w)}$ for all $k\in\\{1,\ldots,r\\}$ by
solving the optimization

$\begin{split}&\mathop{arg\,min}\displaylimits_{\mbox{\boldmath$\theta$}_{k}\in\mathbb{R}^{r},\forall
k\in\\{1,\ldots,r\\}}\sum_{k=1}^{r}-l_{k}(\mbox{\boldmath$\theta$}_{k})+\delta\sum_{k=1}^{r}\|\omega_{k}\|_{1}\\\
&+\frac{\gamma}{2}\sum_{l=1}^{r}\left(\sum_{m\in\hat{P}^{(w)}_{l}}\|\mbox{\boldmath$\theta$}_{m}-\mbox{\boldmath$\theta$}_{l}\|_{2}^{2},\sum_{h\in\hat{N}^{(w)}_{l}}\|\mbox{\boldmath$\theta$}_{h}+\mbox{\boldmath$\theta$}_{l}\|_{2}^{2},\sum_{v\in
Z^{(w)}_{l}}\|\mbox{\boldmath$\theta$}_{l}\|_{2}^{2}\right).\end{split}$ (26)

Similar to the algorithm in Section 2.5, this could be replaced with an
exhaustive search algorithm for the global minimum. However, even for small
$r$ this will be computationally burdensome. The step for estimating the sets
is identical to the previous algorithm and thus once again if the coefficients
are considered given there is an easily derived global solution for the sets,
unlike the K-means algorithm which has been used in similar problems.

For the update in (26) we propose using a proximal gradient descent method
similar to the glmnet algorithm proposed by Friedman et al. (2008). Let
$\hat{w}_{ik}$ and $\hat{z}_{ik}$ be the most recent estimates of $w_{ik}$ and
$z_{ik}$. The update of the $j$th variable of the $m$th response of the
coefficients for the proposed proximal gradient descent method is

$\bar{\theta}_{jm}=\frac{S\left(\sum_{i=1}^{n}\hat{w}_{ik}\\{x_{ij}\hat{z}_{ik}-\sum_{l=1,l\neq
j}^{p}u_{il}u_{ij}\bar{\theta}_{lm}\\}+\gamma M_{jm},I(j\neq
1)\delta/2\right)}{\left(\sum_{i=1}^{n}u_{ij}^{2}+\gamma\\{(r-1)+\sum_{g=1}^{r}I(m\in\hat{P}_{g}^{(w)})+\sum_{v=1}^{r}I(m\in\hat{N}_{v}^{(w)})\\}\right)},$

where

$M_{jm}=\sum_{l\in\hat{P}_{m}^{(w)}}\bar{\theta}_{jl}-\sum_{h\in\hat{N}_{m}^{(w)}}\bar{\theta}_{jh}+\sum_{g=1}^{r}\bar{\theta}_{jg}I(m\in\hat{P}_{g}^{(w)})-\sum_{v=1}^{r}\bar{\theta}_{jv}I(m\in\hat{N}_{v}^{(w)}).$

The algorithm iterates through all $j\in\\{1,\ldots,p\\}$ and
$m\in\\{1,\ldots,r\\}$ until convergence and then updates the quadratic
approximation and continues to solve the IRWLS optimization until convergence.
This entire process is the solution to (26). To select $\delta$ and $\gamma$
we propose using $k$-fold cross validation minimizing validation likelihood
loss.

##  4 OLS Simulations

###  4.1 Data Generating Model and Evaluation Metrics

In this section we investigate the performance of the MinPen estimator defined
by (5), the MinPen estimator with the correct indices matrix ${\bf A}$ known
_a priori_ (T-MinPen), the multivariate cluster elastic net estimator (MCEN)
Price & Sherwood (2018), the separate elastic net estimator (SEN) which fits
the elastic net on each response using common tuning parameters, and joint
elastic net estimator (JEN) which utilizes a group penalty on the same
predictor variable across all responses (Friedman et al., 2008).

The SEN and JEN estimators are fit using the glmnet package in R. Tuning
parameters are selected using a test train procedures, where each model is
trained on $n=100$ observations and then evaluated on a different set of $100$
test observations. The selected tuning parameters minimize the predicted sum
of squares error on the test observations. We establish the T-MinPen estimator
as a baseline for the MinPen estimator when the relationships would be known
by a practitioner prior to fitting the model.

Let $\tilde{\mbox{\boldmath$\Sigma$}}_{X}\in\mathbb{R}^{4\times 4}$ with
entries $\tilde{\sigma}_{ii}=1$ and $\tilde{\sigma}_{ij}=\rho$, for $i\neq j$.
The covariates are generated by ${\bf x}_{i}\sim
N(\mathbf{0}_{p},\mbox{\boldmath$\Sigma$}_{x})$, where
$\mbox{\boldmath$\Sigma$}_{x}$ is a block diagonal matrix with $p/4$ blocks of
$\tilde{\mbox{\boldmath$\Sigma$}}_{X}$ with $\rho=0.7$ on the diagonal blocks
and all other entries are set to 0.

Define that data generating model to the be same as proposed in (1) such that
$\mbox{\boldmath$\epsilon$}_{i}\sim N_{15}(\mathbf{0}_{15},{\bf I}_{15})$. In
all simulations we perform 100 replications where each replication consists of
using training data set with $n=100$ and tuning parameters are selected using
predicted squared error loss on 100 independently sampled observations. To
evaluate the model we generate an independent validation set consisting of
1000 observations and calculate the average squared prediction error (SPE)

$\frac{1}{15000}\sum_{k=1}^{15}\sum_{i=1}^{1000}(y_{ik}-\hat{y}_{ik})^{2},$

where $y_{ik}$ is the $k$th element of ${\bf y}_{i}$ and $\hat{y}_{ik}$ is a
prediction for the $i$th observation and $k$th response. We also use the mean
square error (MSE) of the estimators of the regression coefficient matrix
${\bf B}$, which we defined as

$\frac{1}{15p}\sum_{k=1}^{15}\|\hat{\mbox{\boldmath$\beta$}}_{k}-\mbox{\boldmath$\beta$}_{k}^{*}\|_{2}^{2}.$

We also report on the number of true variables and false variables selected
through true positive (TP) and false positive (FP) rates respectively over the
100 replications.

###  4.2 Block Structure in Regression Coefficients

In this simulation we will investigate how the methods perform when
encountering a structured set of regression coefficients where

${\bf
B}^{*}=\left(\begin{array}[]{lll}\mbox{\boldmath$\Delta$}_{10}(\eta,\lambda)&\vec{0}_{10}&\vec{0}_{10}\\\
\vec{0}_{10}&\mbox{\boldmath$\Delta$}_{10}(\eta,\lambda)&\vec{0}_{10}\\\
\vec{0}_{10}&\vec{0}_{10}&\mbox{\boldmath$\Delta$}_{10}(\eta,\lambda)\\\
\vec{0}_{p-30}&\vec{0}_{p-30}&\vec{0}_{p-30}\end{array}\right).$

We define the regression coefficients for subsets of responses as
$\mbox{\boldmath$\Delta$}_{q}(\eta,\lambda)=(-\mbox{\boldmath$\eta$}_{q}-\lambda,\mbox{\boldmath$\eta$}_{q},\mbox{\boldmath$\eta$}_{q}+\lambda,-\mbox{\boldmath$\eta$}_{q}-2\lambda,\mbox{\boldmath$\eta$}_{q}+3\lambda)\in\mathbb{R}^{q\times
5}$ where $\mbox{\boldmath$\eta$}_{q}$ is a $q$-dimensional vector with each
element set equal to $\eta$ and $\lambda$ is a constant.

The proposed method should identify both the positive and negative
relationships across the responses and leverage that information to improve
prediction accuracy. Of the methods we compare to, only MCEN can identify
relationships and it is limited to only identifying positive relationships. We
investigate the settings where
$(\eta,\lambda)\in\\{0.5,1.0\\}\times\\{0.02,0.05,0.10\\}$ for $p=40,100$, and
$300$. Figure 1 presents the results of the simulation for $p=300$, while the
results for $p=40$ and $100$ are shown in the supplemental material. The
results show that MCEN and MinPen perform competitively in MSE and PSE for all
values of $\eta$, while JEN is superior when $\gamma=1$. With respect to true
and false positive rates, MinPen out performs or performs as well as the other
approaches regardless of the values of $\eta$. It may be somewhat surprising
that MCEN would perform well in this setting, but after investigation MCEN
selects two clusters, a positive and negative cluster. As the effect size is
consistent across each of these models, MCEN is competitive with MinPen though
does tend to have a higher false positive rate.

Figure 1: Results of the simulation described in Section 4.2 for the case of
$p=300$.

The results presented show with respect to SPE all methods except for SEN
perform similarly for all values of $p$, $\eta$, and $\lambda$. A similar
results is found for both MSE and TP, with MinPen having a slight advantage
for smaller $p$. The real differentiation between methods is shown when
comparing the FP rate where we see MCEN and MinPen perform as well as T-MinPen
and out perform all competitors SEN and JEN in the case of $\eta=0.5$. In the
case of $\eta=1.0$ MinPen is able to perform closely to T-MinPen. For larger
values of $\lambda$ we see a much more consistent result in FP in MinPen than
MCEN.

###  4.3  Overlapping Variables in Regression Coefficients

We next investigate a simulation where there is overlap between the non-zero
variables in each coefficient vector. We define ${\bf B}^{*}$ to be be
structured such that the regression coefficients for the $k$th response, with
$r=15$, are defined as

$\mbox{\boldmath$\beta$}^{*}_{k}=\left(\begin{array}[]{l}\vec{0}_{v(k-1)}\\\
(-1)^{k}\vec{0.5}_{10}\\\ \vec{0}_{p-v(k-1)-10}\end{array}\right),$

We investigate the settings $(v,p)\in\\{0,2,4\\}\times\\{100,300\\}$ over 100
replications. Figure 2 presents the case of $p=300$, while the results for
$p=100$ are available in the supplemental material. When $v=0$ each predictor
is either active or inactive for all 15 responses. In this setting JEN out
performs competitors in all metrics for every $p$ studied, because it uses a
group lasso penalty for each predictor across the fifteen different responses.
As the amount of overlap increases MCEN and MinPen methods perform the best
with regard to MSE and PSE. Similar to the simulations presented in Section
4.2 T-MinPen and MinPen out preform competitors finding fewer false positive
variables. Thus using MinPen is superior as the method is able to detect the
underlying structure in the regression coefficients more often than
competitors.

Figure 2: Results of the simulation described in Section 4.3 for the case of
$p=300$

##  5 Binomial Simulations

Next we study the impact of the minimum penalty methodology in binomial
logistic regression setting and compare it to SEN and MCEN which also have
binomial logistic regression cases for the respective methods. We will also
compare against T-MinPen which is a version of MinPen where we fix the true
${\bf A}$ to be known by a practitioner _a priori_. We generate $\\{{\bf
x}_{i}\\}_{i=1}^{n}$ in the exact same way as proposed in Section 4.1. The
responses are generated by

$y_{ik}\sim Bin(1,\pi^{*}_{ik})\mbox{ where }\pi^{*}_{ik}=\frac{\exp({\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}^{*}_{k})}{1+\exp({\bf
x}_{i}^{T}\mbox{\boldmath$\beta$}^{*}_{k})}.$

In this simulation study we define $\mbox{\boldmath$\beta$}^{*}_{k}$
$k\in\\{1,\ldots,r\\}$ to be equivalent to the $k$th column of the ${\bf
B}^{*}$ studied in section 4.2 for $p=40,100,$ and $300$ respectively. Again
we study $(\eta,\lambda)\in\\{0.5,1.0\\}\times\\{0.02,0.05,0.10\\}$. For each
of the 100 iterations tuning parameters are selected using a train test
procedure with 100 observations in the training set, and 100 observations in
the test set. A validation of 1000 observations is to used evaluate methods
using Kullback-Leibler divergence (KL), that is

$\sum_{k=1}^{15}\sum_{i=1}^{1000}\left\\{\log\left(\frac{\hat{\pi}_{ik}}{\pi^{*}_{ik}}\right)\hat{\pi}_{ik}+\log\left(\frac{1-\hat{\pi}_{ik}}{1-\pi^{*}_{ik}}\right)(1-\hat{\pi}_{ik})\right\\},$

where $\hat{\pi}_{ik}$ is the resulting estimated probability. We also compare
methods using MSE, TP, and FP.

Figure 3 presents the results for the case of $p=300$. The results for $p=40$
and $100$ are available in the supplementary material. The results show that
MinPen performs as well as TMinPen and out performs the other competitors with
respect to MSE and KL. We see that with regard to variable selection, MinPen
is comparable, if not better, than the other methods with respect to TP and
outperforms all approaches, even TMinPen, with respect to FP. TMinPen performs
better the MinPen on TP but worse on FP indicating it is selecting more
variables than necessary.

Figure 3: Results of the simulation described in Section 5 for $p=300$.

##  6 Applied Examples

###  6.1 Genomics Data

In this section we compare the post model selection inference approach
outlined in this paper, after using the proposed least squares method, with
the post model selection inference approach proposed by Lee et al. (2016),
after using lasso for model selection. The models are fit to data analyzed by
Votavova et al. (2011), who collected demographic, birth and gene expression
data from 72 postpartum women and their newborns, but our analysis is limited
to 64 of the women, 65 had complete data and one was dropped due to outliers
in the demographic data. Four response variables are modeled: placental
weight, newborn weight and two measures of cotinine level, one from the
mother’s peripheral blood and a second from the umbilical cord. The predictors
are Smoking status, mother’s age, mother’s BMI, parity, gestational age and 33
gene expression probes. The 33 probes were selected by taking the absolute
value of the correlation for each response and the 24,526 probes measured in
the study and using a union of the top ten for each response.

The lasso models are fit separately for each response. For both methods tuning
parameters are selected to minimize the mean squared prediction error from
five folds cross validation. For post selection inference,
$\tilde{\mbox{\boldmath$\Sigma$}}$ is estimated using the covariance matrix of
residuals from the full multivariate linear regression model. The four
responses are related. First, the two cotinine measurements are measurements
of the same variable at nearly the same time, but from different samples. In
addition, cotinine levels are high in people who smoke and it has been shown
that lower cotinine levels are associated with larger birth weights (Li et
al., 1993). In addition, Wang et al. (2014) found both smoking has a negative
relationship with both birth and placenta weights. Finally, newborn and
placental weights are positively correlated (Molteni et al., 1978). The
analysis of this data set serves two purposes. First, we want to compare post
selection inference results with an existing method. Second, we want to verify
that minimizing (7) with real data can provide a sensible graph of the
responses. For the latter issue the answer is yes. Table 1 provides the
entries of $\hat{{\bf A}}$ and aligning with the scientific literature the
optimal penalty is to penalize a difference in the weight measurements and a
difference in the cotinine measurements, while all other comparisons the
penalty enforces a negative correlation between the measurements.

| Placental Weight | Newborn Weight | Blood Cotinine | Cord Cotinine
---|---|---|---|---
Placental Weight | 0 | -1 | 1 | 1
Newborn Weight | -1 | 0 | 1 | 1
Blood Cotinine | 1 | 1 | 0 | -1
Cord Cotinine | 1 | 1 | -1 | 0

Table 1:  The ${\bf A}$ matrix for the genomics data

Figure 4 provides the post selection confidence intervals for the lasso and
MinPen methods. The selected models are very similar with the biggest model
difference being in the model of the weight for placenta, where MinPen selects
four more variables include the mother’s age which is the only demographic
variable difference between the models. In addition, for this model the post
selection results for MinPen find a negative relationship between smoking and
placenta weight, but the lasso results are inconclusive. The previously cited
literature appears indicates that for this particular relationship the lasso
approach suffers from a Type II error, while MinPen correctly identifies this
relationship. However, for modeling cotinine level in the umbilical cord, the
one setting where the variables selected are the same, the lasso confidence
intervals are noticeably smaller.

Figure 4: Post Selection Confidence Intervals: Clockwise starting in the top
left corner; (1) Newborn weight; (2) Placenta Weight; (3) Umbilical cord
cotinine level; (4) Peripheral blood cotinine level

###  6.2 Connecticut Multiple Substance Use Data

We investigate the use of the minimum penalty to identify co-occurring drug
use based on accidental drug related deaths in Connecticut from 2012-2018. The
data set consists of 5097 accidental overdoses in the state of Connecticut and
corresponding toxicology, death certificate, and scene investigation
information from the Chief Medical Examiner Office, publicly available at
https://data.ct.gov and included in the supplemental material. We identify 17
response variables which are 16 individual drugs that were investigated plus a
response that identifies if any opiate is involved. The goal of this analysis
is to identify how different drugs appear (or do not appear) in cases of
overdose, with a hope of better understanding and being able to predict co-
occurrence of drugs. By being able to predict co-occurrence of drugs, both
clinical practitioners and law enforcement are better able to understand
emerging trends in drug usage allowing them to quickly identify, intervene,
and mitigate issues.

The covariates in this data are indicators of the city and state of the
residence and death, along with using a bag of words approach to find common
language used in the cause of death and description of injury provided by
police and the medical examiner. As the goal is to better understand the
relationships between drugs being used, information in the cause of death may
not be relevant to all drugs found in the system. In total 111 covariates were
used in this analysis. A validation set was created from 500 randomly selected
observations from 2018 tuning parameters were selected using a 90/10
train/test approach using total classification error rate. Once the tuning
parameters were obtained the model was fit using the chosen tuning parameters
on all data less the validation set. We fit MinPen, MCEN, and elastic net
models using a common sparsity parameter for all 17 responses, and compare
methods using ROC curves on the validation set, see Figure 5. The results show
the MinPen is competitive or out performs MCEN and the elastic net approach.
Furthermore, MinPen provides relationships between responses. Only positive or
no relationships were found between responses, that is no negative
relationships were found. The supplementary material includes a graph of these
relationships. Most notably heroin is not predictive or related to fentanyl,
fentanyl analogue, benzodiazapine, or inform the presence of other opioids.
The opioid variable was mostly related to prescription based drugs such as
oxymorphone, hydrocodone, methadone, among others. Methadone was positively
related to every other variable which is notable because it is used as a
medical assisted treatment for opiate abuse. Further, the model finds
benzodiazapine is not related to cocaine, which is of interest as
benzodiasipine is used to treat cocaine toxicity. The results show MinPen
provides better prediction accuracy than both MCEN and elastic net, providing
evidence that use of the minimum penalty can be beneficial in the presence of
complex and potentially unknown relationships. In the supplemental material we
also present a heatmap of the coefficients produced by MinPen for each
response.

Figure 5: Results of the validation set of 500 observations from 2018 with
regard to four substances: Cocaine, Fentanyl Analogue, Methadone, and the
ability to detect any opioid. The results show that MinPen performs as well if
not better with regard to AUC in classification.

## References

  * (1)
  * Anderson et al. (1951) Anderson, T. W. et al. (1951), ‘Estimating linear restrictions on regression coefficients for multivariate normal distributions’, The Annals of Mathematical Statistics 22(3), 327–351.
  * Bühlmann et al. (2013) Bühlmann, P., Rütimann, P., van de Greer, S. & Zhang, C.-H. (2013), ‘Correlated variables in regression: Clustering and sparse estimation’, Journal of Statistical Planning and Inference 143(11), 1835–1858.
  * Chen & Huang (2016) Chen, L. & Huang, J. Z. (2016), ‘Sparse reduced-rank regression with covariance estimation’, Statistics and Computing 26(1-2), 461–470.
  * Chen et al. (2016) Chen, Y., Iyengar, R. & Iyengar, G. (2016), ‘Modeling multimodal continuous heterogeneity in conjoint analysis—a sparse learning approach’, Marketing Science 36(1), 140–156.
  * Friedman et al. (2008) Friedman, J., Hastie, T. & Tibshirani, R. (2008), ‘Regularized paths for genearlized linear models via coordinate descent’, Journal of Statistial Softwawre 33(1).
  * Friedman et al. (2010) Friedman, J., Hastie, T. & Tibshirani, R. (2010), ‘Regularization paths for generalized linear models via coordinate descent’, Journal of Statistical Software 33(1), 1–22.
  * Hebiri & van de Geer (2011) Hebiri, M. & van de Geer, S. (2011), ‘The smooth-lasso and other $\ell_{1}$+$\ell_{2}$-penalized methods’, Electronic Journal of Statistics 5, 1184 – 1226.
  * Kim et al. (2012) Kim, S., Xing, E. P. et al. (2012), ‘Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eQTL mapping’, The Annals of Applied Statistics 6(3), 1095–1117.
  * Lee et al. (2016) Lee, J. D., Sun, D. L., Sun, Y. & Taylor, J. E. (2016), ‘Exact post-selection inference, with application to the lasso’, The Annals of Statistics 44(3), 907–927.
  * Lee & Liu (2012) Lee, W. & Liu, Y. (2012), ‘Simultaneous multiple response regression and inverse covariance matrix estimation via penalized gaussian maximum likelihood’, Journal of Multivariate Analysis 111, 241–255.
  * Li & Li (2008) Li, C. & Li, H. (2008), ‘Network-constrained regularization and variable selection of genomic data’, Bioinformatics 24(9), 1175–1182.
  * Li & Li (2010a) Li, C. & Li, H. (2010a), ‘Supplamental material for “Variable selection and regression analysis for graph-structure covariates with an application to geneomics”’, The Annals of Applied Statistics 4(3), 1–5.
  * Li & Li (2010b) Li, C. & Li, H. (2010b), ‘Variable selection and regression analysis for graph-structure covariates with an application to geneomics’, The Annals of Applied Statistics 4(3), 1498–1516.
  * Li et al. (1993) Li, C. Q., Windsor, R. A., Perkins, L., Goldenberg, R. L. & Lowe, J. B. (1993), ‘The Impact on Infant Birth Weight and Gestational Age of Cotinine-Validated Smoking Reduction During Pregnancy’, JAMA 269(12), 1519–1524.
  * Li et al. (2015) Li, Y., Nan, B. & Zhu, J. (2015), ‘Multivariate sparse group lasso for multivariate multiple linear regression with arbitrary group sparsity’, Biometrics 71, 354–363.
  * Molteni et al. (1978) Molteni, R. A., Stys, S. J. & Battaglia, F. C. (1978), ‘Relationship of fetal and placental weight in human beings: fetal/placental weight ratios at various gestational ages and birth weight distributions’, The Journal of Reproductive Medicine 21, 327–334.
  * Negahban et al. (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J. & Yu, B. (2012), ‘A unified framework for high-dimensional analysis fo $m$-estimators with decomposable regualrizers’, Statistical Science 27(4), 538–557.
  * Price, Allenbrand & Sherwood (2021) Price, B. S., Allenbrand, C. & Sherwood, B. (2021), ‘Dececting clusters in multivariate regression’, WIREs Computational Statistics To Appear.
  * Price, Molstad & Sherwood (2021) Price, B. S., Molstad, A. J. & Sherwood, B. (2021), ‘Estimating multiple precision matrices using cluster fusion regularization’, Journal of Computational and Graphical Statistics, To Appear .
  * Price & Sherwood (2018) Price, B. S. & Sherwood, B. (2018), ‘A cluster elastic net for multivariate regression’, Journal of Machine Learning Research 19, 1–37.
  * Rothman et al. (2010) Rothman, A. J., Levina, E. & Zhu, J. (2010), ‘Sparse multivariate regression with covariance estimation’, Journal of Computational and Graphical Statistics 19(4), 947–962.
  * Tibshirani (2013) Tibshirani, R. J. (2013), ‘The lasso problem and uniqueness’, Electronic Journal of Statistics 7, 1456–1490.
  * Velu & Reinsel (2013) Velu, R. & Reinsel, G. C. (2013), Multivariate reduced-rank regression: theory and applications, Vol. 136, Springer Science & Business Media.
  * Votavova et al. (2011) Votavova, H., Merkerova, M. D., Fejglova, K., Vasikova, A., Krejcik, Z., Pastorkova, A., Tabashidze, N., Topinka, J., Jr., M. V., Sram, R. & Brdicka, R. (2011), ‘Transcriptome alterations in maternal and fetal cells induced by tobacco smoke’, Placenta 32, 763–770.
  * Wang et al. (2014) Wang, N., Tikellis, G., Sun, C., Pezic, A., Wang, L., Wells, J., Cochrane, J., Ponsonby, A.-L. & Dwyer, T. (2014), ‘The effect of maternal prenatal smoking and alcohol consumption on the placenta-to-birth weight ratio’, Placenta 35(7), 437–441.
  * Witten et al. (2014) Witten, D. M., Shojaie, A. & Zhang, F. (2014), ‘The cluster elastic net for high-dimensional regression with unknown variable grouping’, Technometrics 56(1), 112–122.
  * Witten & Tibshirani (2009) Witten, D. & Tibshirani, R. (2009), ‘Covariance regularized regression and classification for high-dimensional problems’, Journal of Royal Statistical Society, Series B 71(3), 615–636.
  * Xu et al. (2015) Xu, L., Huang, A., Chen, J. & Chen, E. (2015), Exploiting task-feature co-clusters in multi-task learning, in ‘Twenty-Ninth AAAI Conference on Artificial Intelligence’.
  * Zhao & Yu (2006) Zhao, P. & Yu, B. (2006), ‘On model selection consistency of lasso’, Journal of Machine Learning Research 7, 2541–2563.
  * Zhao & Shojaie (2016) Zhao, S. & Shojaie, A. (2016), ‘A significance test for graph constrained estimation’, Biometrics 72(2), 484–493.
  * Zhou et al. (2017) Zhou, W., Sherwood, B., Ji, Z., Xue, Y., Du, F., Bai, J., Ying, M. & Ji, H. (2017), ‘Genome-wide prediction of DNase I hypersensitivity using gene expression’, Nature Communications 8, 1–17.
  * Zou & Hastie (2005) Zou, H. & Hastie, T. (2005), ‘Regularization and variable selection via the elastic net’, J. R. Statist. Soc. B 67, 301–320.