Reducing bias and mitigating the influence of excess of zeros in regression covariates with multi-outcome adaptive LAD-lasso

Abstract

Zero-inflated explanatory variables, as opposed to outcome variables, are common, for example, in environmental sciences. In this article, we address the problem of having excess of zero values in some continuous explanatory variables, which are subject to multi-outcome lasso-regularized variable selection. In short, the problem results from the failure of the lasso-type of shrinkage methods to recognize any difference between zero value occurring either in the regression coefficient or in the corresponding value of the explanatory variable. This kind of confounding will obviously increase the number of false positives – all non-zero regression coefficients do not necessarily represent true outcome effects. We present here the adaptive LAD-lasso for multiple outcomes, which extends the earlier work of multi-outcome LAD-lasso with adaptive penalization. In addition to well-known property of having less biased regression coefficients, we show that the adaptivity also improves method’s ability to recover from influences of excess of zero values measured in continuous covariates.

Keywords:

• Multivariate analysis
• $p\gg n$ regression
• penalized regression
• robust procedures
• variable selection
• zero-inflated continuous data

1 Introduction

In high-dimensional regression problems, the number of parameters is often larger than the number of individuals in the sample ($p\gg n$). Ordinary least squares and other conventional estimation methods do not work in such situations. Simultaneous estimation and variable selection using lasso (least absolute shrinkage and selection operator) (Tibshirani 1996; Li and Sillanpää 2012) is a popular shrinkage-estimation approach to obtain sparse estimates for regression coefficients in high-dimensional regression problems with desirable statistical properties (Fan and Li 2001; Fan and Peng 2004). However, the generally known drawback of shrinkage-inducing methods, including lasso, is that they improve estimation accuracy but also introduce downward bias to the estimates (due to bias-variance tradeoff). To alleviate this problem, a re-estimation procedure has been proposed, where effects of selected predictors are re-estimated using no penalty (Efron et al. 2004; Meinshausen 2007). Adaptive lasso has been also presented (Zou 2006), so that the selected predictors are subject of reduced penalty and non-selected positions obtain heavier penalty. This is performed either by using two-stage or iterative estimation strategy. See also related work on thresholded lasso (Zhou 2010; Van De Geer, Bühlmann, and Zhou 2011). One would expect to see larger variance as a consequence of reduced bias in the adaptive lasso. However, it has been observed that the opposite could happen in the adaptive lasso context meaning that we may see improved accuracy together with reduced bias, especially in the $p\gg n$ situation (Huang, Ma, and Zhang 2008).

When the data set contains some outlying observations, it is possible to apply robust least absolute deviation (LAD) regression instead of ordinary regression. In lasso context, LAD-lasso has been proposed for univariate (Wang, Li, and Jiang 2007) and multi-outcome cases (Möttönen and Sillanpää 2015; Li, Möttönen, and Sillanpää 2015). The multi-outcome LAD-lasso is related to the group lasso (Yuan and Lin 2006) because the object function of multi-outcome LAD-lasso also contains group lasso penalty. As common shrinkage-inducing methods, LAD-lasso methods are also suffering from downward bias of the estimates. To alleviate this, adaptive LAD-lasso has been proposed for univariate LAD-lasso (Arslan 2012). In this article, we study the multi-outcome adaptive LAD-lasso and make comparisons to the non-adaptive or standard version, which we refer to as ordinary LAD-lasso. Note that adaptive group lasso has been presented by Wang and Leng (2008), and it has already been used in empirical application and in modeling old-age retirement (Lähderanta et al. 2022).

In our earlier work with empirical data, we located a potential problem in lasso methods in general, when there are excess of real zero values in some covariates (Lähderanta et al. 2022). We describe this problem here in more detail in LAD-lasso and provide some potential solutions using controlled simulation and real data example. In short, the problem occurs as the lasso-type of shrinkage methods do not recognize any difference between zero value occurring either at the regression coefficient or in the corresponding value of the covariate. As a result, there is a danger of misinterpretation because all non zero regression coefficients do not necessarily represent true effects on the outcomes. As a potential solution, we show empirically how the use of adaptive LAD-lasso will minimize this kind of confounding to occur in practice.

In this article, we show the potential benefits of the multi-outcome adaptive LAD-lasso methods with three examples. In the first example, we utilize a public genotype data set (Crooks et al. 2009) to simulate new phenotypes with tri-variate traits. In the second example, we simulate several general data sets with zero-inflated (excess of zeros) values in continuous explanatory variables, which is typical in many fields of science. Examples of zero-inflated data can be found in environmental sciences (Maldonado, Aguilera, and Salmerón 2016; Kamarianakis et al. 2008). We demonstrate that this problem can be mitigated with the adaptive version of the lasso. In the third example, we use a real-life data set to study the main socioeconomic factors that affect pensions and length of working lives of disability pension retirees, both having highly skewed outcome distributions.

2 Methodology

In this section, we present the multi-outcome ordinary LAD-lasso (Möttönen and Sillanpää 2015; Li, Möttönen, and Sillanpää 2015) and the novel adaptive version of the method.

2.1 Multi-outcome ordinary LAD-lasso

Consider multiple regression model with multiple outcomes$$\mathbf{Y}=\mathbf{XB}+\mathbf{E},$$where $\mathbf{Y}={({\mathbf{y}}_1,\dots ,{\mathbf{y}}_n)}^{\prime }$ is an $n\times q$ matrix of $n$ observed values of $q$ outcome variables, $\mathbf{X}={({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)}^{\prime }$ is an $n\times p$ matrix of $n$ observed values of $p$ explanatory variables, $\mathbf{B}={({\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_p)}^{\prime }$ is a $p\times q$ matrix of regression coefficients, and $\mathbf{E}={({\mathit{\epsilon }}_1,\dots ,{\mathit{\epsilon }}_n)}^{\prime }$ is an $n\times q$ matrix of residuals. We further assume that ${\mathit{\epsilon }}_1,\dots ,{\mathit{\epsilon }}_n$ is a random sample of size $n$ from a $q$-variate distribution centered at the origin.

The multi-outcome lasso estimation method is based on the penalized objective function(1) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i|{|}^2+\lambda \sum _{j=2}^p|\left|{\mathit{\beta }}_j\right||$$(1)

(e.g., Turlach et al. 2005; Yuan and Lin 2006; Yuan et al. 2007). The minimizer of the objective function (1)(1) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i|{|}^2+\lambda \sum _{j=2}^p|\left|{\mathit{\beta }}_j\right||$$(1) gives now multi-outcome lasso estimate for the regression coefficient matrix $\mathbf{B}$. Note that Equation (1)(1) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i|{|}^2+\lambda \sum _{j=2}^p|\left|{\mathit{\beta }}_j\right||$$(1) is also an objective function of group lasso where each of $p-1$ explanatory variables (the intercept terms are omitted) have $q$ regression coefficients (one for each outcome variable) which form a group.

The multi-outcome lasso method gives sparse solutions but it is obviously not very robust. A more robust version of multi-outcome lasso, the multi-outcome LAD-lasso (Möttönen and Sillanpää 2015; Li, Möttönen, and Sillanpää 2015), can be obtained by replacing the squared norm with an L1 norm in the objective function (1)(1) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i|{|}^2+\lambda \sum _{j=2}^p|\left|{\mathit{\beta }}_j\right||$$(1) :(2) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i\left|\right|+\lambda \sum _{j=2}^p\left|\right|{\mathit{\beta }}_j\left|\right|.$$(2)

It is easily seen that if we define$$\left(\begin{array}{c}{\mathbf{y}}_i^{*}\\ {\mathbf{x}}_i^{*}\end{array}\right)=\{\begin{array}{cc}\left(\begin{array}{c}{\mathbf{y}}_i\\ {\mathbf{x}}_i\end{array}\right),\hfill & \text{if}i=1,\dots ,n,\hfill \\ \hfill & \hfill \\ \left(\begin{array}{c}\mathbf{0}\\ n\lambda {\mathbf{e}}_{i-n+1}\end{array}\right),\hfill & \text{if}i=n+1,\dots ,n+p-1,\hfill \end{array}$$where ${\mathbf{e}}_i$ is a $p\times 1$ unit vector with $i$ th element as one and all the other elements are zero and $\mathit{0}$ is a $q\times 1$ vector of zeros, the objective function (2)(2) $$\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i\left|\right|+\lambda \sum _{j=2}^p\left|\right|{\mathit{\beta }}_j\left|\right|.$$(2) reduces to the LAD estimation objective function (Oja 2010) of the form(3) $$\frac{1}{n+p-1}\sum _{i=1}^{n+p-1}\left|\right|{\mathbf{y}}_i^{*}-{\mathbf{B}}^{\prime }{\mathbf{x}}_i^{*}\left|\right|,$$(3) which shows that we can use any multi-outcome LAD regression estimation routine to find the multi-outcome LAD-lasso estimates. One can, for example, use the function mv.l1lm of the R-package MNM (Nordhausen, Möttönen, and Oja 2016; Nordhausen and Oja 2011).

2.2 Multi-outcome adaptive LAD-lasso

The multi-outcome LAD-lasso estimate for a fixed $\lambda$ can be defined as$$\widehat{\mathbf{B}}(\lambda )=\underset{{}_{\mathbf{B}}}{\text{argmin}}\left[\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i\left|\right|+\lambda \sum _{j=2}^p\left|\right|{\mathit{\beta }}_j\left|\right|\right],$$where argmin stands for” the value minimizing”. This brings up the question of how to choose the tuning parameter $\lambda$? If we are mainly concerned about recovering the right model, then we can use, for example, Akaike’s information criterion (AIC) or Bayesian information criterion (BIC). On the other hand, if we are mainly concerned about prediction accuracy, then cross-validation technique (CV) is often a good choice. For a general survey on CV methods, see Arlot and Celisse (2010). We continue this study with a CV criterion as robust AIC or BIC are not readily available in the multi-outcome LAD-lasso context. We have used fivefold cross-validation where the original sample of $n$ observations is randomly partitioned into five approximately equal sized subsamples. Let $C_j$ denote the set of indices of the observations that belong to the $j$ th subsample. The CV criterion is then defined as the mean absolute error over the five subsamples$$\text{CV}(\widehat{\mathbf{B}}(\lambda ))=\frac{1}{5}\sum _{j=1}^5\left[\frac{1}{n_j}\sum _{i\in C_j}\left|\right|{\mathbf{y}}_i-{\widehat{\mathbf{B}}}_{(-j)}^{\prime }{\mathbf{x}}_i\left|\right|\right],$$ where $n_j$ is the number of observations in the $j$ th subsample and ${\widehat{\mathbf{B}}}_{(-j)}$ is the multi-outcome LAD-lasso estimate using all but the observations in the $j$ th subsample. For a BIC-like criterion in multi-outcome LAD-lasso context, see Möttönen and Sillanpää (2015).

Let$$\widehat{\lambda }=\underset{{}_{\lambda }}{\text{argmin}}\text{\hspace{0.17em}CV}(\widehat{\mathbf{B}}(\lambda ))$$be the value of the tuning parameter $\lambda$ which minimizes the cross-validation criterion function for the multi-outcome LAD-lasso. The multi-outcome LAD-lasso estimate can then be defined as$$\widehat{\mathbf{B}}=\widehat{\mathbf{B}}(\widehat{\lambda }).$$

It has been shown that lasso estimation tends to underestimate the regression coefficients and the same is true for the multi-outcome LAD-lasso estimate $\widehat{\mathbf{B}}$. Zou (2006) proposed an adaptive lasso method, which gives an estimate whose bias is smaller than that of the ordinary lasso. Similarly, the adaptive method can also be used in the LAD-lasso (Arslan 2012). We extend this method here for multi-outcome LAD-lasso case. Since the performance of the adaptive LAD-lasso might be sensitive to the initial weights, we propose an iterative $k$-step adaptive estimation method. The $s$ th ($s=1,\dots ,k$) step multi-outcome adaptive LAD-lasso estimate for fixed $\lambda$ is defined as$${\widehat{\mathbf{B}}}^{(s)}(\lambda )=\underset{{}_{\mathbf{B}}}{\text{argmin}}\left[\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i\left|\right|+\lambda \sum _{j=2}^p{w}_j^{(s-1)}\left|\right|{\mathit{\beta }}_j\left|\right|\right],$$where$$w_j^{(s-1)}=\frac{1}{\left|\right|{\widehat{\mathit{\beta }}}_j^{(s-1)}\left|\right|+1/n},j=2,\dots ,p,$$

${\widehat{\mathit{\beta }}}_j^{(s-1)}$ is the $j$ th row of the $(s-1)$ th step multi-outcome LAD-lasso estimate ${\widehat{\mathbf{B}}}^{(s-1)}={\widehat{\mathbf{B}}}^{(s-1)}({\widehat{\lambda }}^{(s-1)})$ and ${\widehat{\lambda }}^{(s-1)}={\text{argmin}}_{\lambda }\text{\hspace{0.17em}CV}({\widehat{\mathbf{B}}}^{(s-1)}(\lambda ))$. The multi-outcome LAD-lasso estimate $\widehat{\mathbf{B}}$ is used as the initial value ${\widehat{\mathbf{B}}}^{(0)}$. It can be easily shown that the $s$ th step of the multi-outcome adaptive LAD-lasso estimate for fixed $\lambda$ can be written in the form$${\widehat{\mathbf{B}}}^{(s)}(\lambda )=\underset{{}_{\mathbf{B}}}{\text{argmin}}\left[\frac{1}{n}\sum _{i=1}^n\left|\right|{\mathbf{y}}_i-{\mathbf{B}}^{\prime }{\mathbf{x}}_i\left|\right|+\frac{1}{n}\sum _{j=2}^p\left|\right|\mathbf{0}-{\mathbf{B}}^{\prime }n\lambda w_j^{(s-1)}{\mathbf{e}}_j\left|\right|\right],$$which further implies that it can be written in the multi-outcome LAD regression form(4) $${\widehat{\mathbf{B}}}^{(s)}(\lambda )=\underset{{}_{\mathbf{B}}}{\text{argmin}}\left[\frac{1}{n+p-1}\sum _{i=1}^{n+p-1}\left|\right|{\mathbf{y}}_i^{*}-{\mathbf{B}}^{\prime }{\mathbf{x}}_i^{*}\left|\right|\right],$$(4) where$$\left(\begin{array}{c}{\mathbf{y}}_i^{*}\\ {\mathbf{x}}_i^{*}\end{array}\right)=\{\begin{array}{cc}\left(\begin{array}{c}{\mathbf{y}}_i\\ {\mathbf{x}}_i\end{array}\right),\hfill & \text{if}i=1,\dots ,n,\hfill \\ \hfill & \hfill \\ \left(\begin{array}{c}\mathbf{0}\\ n\lambda w_j^{(s-1)}{\mathbf{e}}_{i-n+1}\end{array}\right),\hfill & \text{if}i=n+1,\dots ,n+p-1.\hfill \end{array}$$

The estimate ${\widehat{\mathbf{B}}}^{(s)}(\lambda )$ can then be found using the function mv.l1lm of the R-package MNM (Nordhausen, Möttönen, and Oja 2016; Nordhausen and Oja 2011).

3 Simulation studies

In this section, we present two simulation experiments with simulated data sets to provide insights to the multi-outcome adaptive LAD-lasso in certain scenarios. R code for the simulations can be found on GitHub (Lähderanta and Möttönen 2021).

3.1 Bias reduction

The first experiment stems from medical sciences, where single and multi-trait quantitative trait locus (QTL) mapping is commonly applied. In QTL mapping, for each individual in the population sample, we have quantitative trait phenotypes and genotype patterns measured in a high number of marker loci distributed over the genome. Based on this data, the target is to find small subsets of trait-associated loci out of many candidates. The trait-associated findings, as well as, the true causative loci, are both interchangeably called QTLs. More details of QTL mapping can be found in Balding (2006) and Li and Sillanpää (2012).

As a first step, we simulated new phenotypes with tri-variate traits using the public genotype data set from the 12th QTL-MAS workshop in Uppsala, Sweden, in 2008 (Crooks et al. 2009). The original genotype data set contains of 5865 individuals and 6000 markers. We took a random sample of 300 individuals and chose every 30th marker with a resulting total number of 200 markers. We then simulated tri-variate traits by using a multi-outcome multiple regression model$$\mathbf{Y}=\mathbf{XB}+\mathbf{E},$$where $\mathbf{Y}$ is a $300\times 3$ matrix of tri-variate traits, $\mathbf{X}$ is a $300\times 200$ matrix with $ij$ th element$$x_{ij}=\{\begin{array}{cc}-1,& \text{if}\mathrm{ }\text{individual}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{homozygote}\mathrm{ }\mathrm{(}11\mathrm{)}\mathrm{ }\text{at}\mathrm{ }\text{marker}\mathrm{ }j,\\ 0,& \text{if}\mathrm{ }\text{individual}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{heterozygote}\mathrm{ }\text{at}\mathrm{ }\text{marker}\mathrm{ }j,\\ 1,& \text{if}\mathrm{ }\text{individual}\mathrm{ }i\mathrm{ }\text{is}\mathrm{ }\text{homozygote}\mathrm{ }\mathrm{(}22\mathrm{)}\mathrm{ }\text{at}\mathrm{ }\text{marker}\mathrm{ }j,\end{array}$$

$\mathbf{B}$ is a $200\times 3$ matrix with four QTLs indicated as non-zero rows$$\begin{array}{c}{\mathit{\beta }}_{50}^{\prime }=(1.5,1.5,1.5),{\mathit{\beta }}_{75}^{\prime }=(0,0.5,1.5),\\ {\mathit{\beta }}_{100}^{\prime }=(0.5,0.5,0.5)\text{and}{\mathit{\beta }}_{150}^{\prime }=(1,0.5,0.5),\end{array}$$and $\mathbf{E}$ is a $300\times 3$ matrix with i.i.d. rows normally distributed as$$N\left(\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right),\left(\begin{array}{ccc}1.0& 0.5& 0.3\\ 0.5& 1.0& 0.2\\ 0.3& 0.2& 1.0\end{array}\right)\right).$$

We estimated the tuning parameter $\lambda$ by using fivefold cross-validation. In Figure 1a, we show the marker effects ${\beta }_{ij}$ corresponding to ordinary LAD-lasso. We see that the LAD-lasso method correctly finds all four QTLs.

Fig. 1 Marker effects ${\widehat{\beta }}_{ij}$ with absolute values greater than ${10}^{-6}$. (a) Multi-outcome ordinary LAD-lasso estimates and (b) Multi-outcome adaptive LAD-lasso estimates (five steps).

In the next step, we studied the bias of the LAD-lasso estimates of the non-zero coefficient vectors ${\mathit{\beta }}_{50}$, ${\mathit{\beta }}_{75}$, ${\mathit{\beta }}_{100}$ and ${\mathit{\beta }}_{150}$:$$\left(\begin{array}{c}{\widehat{\mathit{\beta }}}_{50}^{\prime }-{\mathit{\beta }}_{50}^{\prime }\\ {\widehat{\mathit{\beta }}}_{75}^{\prime }-{\mathit{\beta }}_{75}^{\prime }\\ {\widehat{\mathit{\beta }}}_{100}^{\prime }-{\mathit{\beta }}_{100}^{\prime }\\ {\widehat{\mathit{\beta }}}_{150}^{\prime }-{\mathit{\beta }}_{150}^{\prime }\end{array}\right)\mathrm{ }=\left(\begin{array}{ccc}-0.13& -0.22& -0.20\\ 0.01& -0.26& -0.60\\ -0.22& -0.18& -0.08\\ -0.16& -0.18& -0.10\end{array}\right).$$

It is clear that the ordinary LAD-lasso estimates are biased in most cases. In the final step, we used the multi-outcome adaptive LAD-lasso estimation method. The Figure 2a shows that the biases decrease and the regression coefficients stabilize after some adaptive steps. We can see from Figure 2b that there are still small biases after five adaptive steps but the biases are smaller than the biases corresponding to the ordinary LAD-lasso. Note that with adaptive LAD-lasso the accuracy of the estimates has improved on average. The biases of the non-zero coefficient vectors corresponding to the adaptive LAD-lasso with five steps were$$\left(\begin{array}{c}{\widehat{\mathit{\beta }}}_{50}^{\prime }-{\mathit{\beta }}_{50}^{\prime }\\ {\widehat{\mathit{\beta }}}_{75}^{\prime }-{\mathit{\beta }}_{75}^{\prime }\\ {\widehat{\mathit{\beta }}}_{100}^{\prime }-{\mathit{\beta }}_{100}^{\prime }\\ {\widehat{\mathit{\beta }}}_{150}^{\prime }-{\mathit{\beta }}_{150}^{\prime }\end{array}\right)\mathrm{ }=\left(\begin{array}{ccc}-0.01& -0.06& -0.01\\ 0.10& -0.15& -0.10\\ -0.11& -0.06& 0.07\\ -0.02& -0.13& -0.03\end{array}\right).$$

Fig. 2 Biases of the marker effects $\text{bias}({\widehat{\beta }}_{ij})={\widehat{\beta }}_{ij}-{\beta }_{ij}$, $i=50,75,100,150$, $j=1,2,3$. (a) Biases of the marker effects with different numbers of adaptive steps. Zero adaptive steps corresponds to the ordinary LAD-lasso estimate. (b) Box-plots of the biases of the marker effects for the ordinary LAD-lasso and the adaptive LAD-lasso with five steps.

In Figure 1b, we show the marker effects ${\beta }_{ij}$ corresponding to the adaptive LAD-lasso with five steps. We see that the multi-outcome adaptive LAD-lasso method correctly finds the true QTLs. Note that the true zeros indicated by the adaptive technique are not shown in Figure 1b.

3.2 Excess of zeros in covariates

In the next experiment, we simulate an other kind of data sets to demonstrate the excess of zeros in covariates scenario with multi-outcome adaptive LAD-lasso and multi-outcome ordinary LAD-lasso.

Results from multiple data sets are shown to illustrate the performance of the algorithms in wide variety of situations. In practice, we generate two types of data sets, one with high number of observations compared to covariates ($n=100,p=10$), and contrarily one with high number of covariates when compared to observations ($n=25,p=50$). Moreover, we alternate the proportion of zeros (${\pi }_{\mathit{zeros}}$) in the covariates from 0.1 to 0.4 and apply two different error E distributions: uniform and asymmetric Laplace. When all the combinations of these properties are considered, 16 different types of scenarios are examined in total. In each scenario, we simulate 100 data sets to further evaluate the performance of the methods. The data sets and the alternating proportion of zeros are shown in Figure 3.

Fig. 3 Simulated data replicates in the excess of zeros study with different proportions of zeros (${\pi }_{\mathit{zeros}}$).

In practice, the data sets are generated with a multi-outcome regression model$$\mathbf{Y}=\mathbf{XB}+\mathbf{E},$$with matrix $\mathbf{X}$ of a size $n\times p$, response matrix $\mathbf{Y}$ of a size $n\times 2$. The observed values $x_{ij}$ are simulated from a normal distribution such as$$x_{ij}\sim \{\begin{array}{cc}N(3,1),& \text{if}\mathrm{ }{u}_{ij}>{\pi }_{\mathit{zeros}},\\ 0,& \text{else},\end{array}$$ where $u_{ij}\sim \mathit{Unif}(0,1)$ and ${\pi }_{\mathit{zeros}}$ is the proportion of zeros.

$\mathbf{B}$ is a $p\times 2$ matrix where$$b_{jk}\sim \{\begin{array}{cc}N(0,1),& \text{when}\mathrm{ }j=1,2,3,\\ 0,& \text{else}.\end{array}$$

$\mathbf{E}$ is a $n\times 2$ matrix where$$e_{jk}\sim \mathit{Unif}(0,1)$$or$$e_{jk}\sim \mathit{ALaplace}(\mu =[3,6],\Sigma =\left(\begin{array}{cc}0.5& 0\\ 0& 0.5\end{array}\right)),$$ depending of the choice for error distribution. Above, Unif refers to uniform distribution and ALaplace stands for Asymmetric Laplace distribution.

From the simulations, we can see that the adaptive LAD-lasso is superior to ordinary LAD-lasso in every scenario, when we compare the correctly found zero coefficients (Figure 4). In scenario $p\gg n$, the difference between the methods is somewhat smaller.

Fig. 4 Percentage of correctly found zero coefficients with adaptive LAD-lasso and ordinary LAD-lasso in multi-outcome context.

4 Empirical data example: retirement analysis

In this section, we apply the presented adaptive method to model retirement behavior. The data consist of disability pension retirees in Finnish statutory earnings-related pension system for the public sector employees. A random sample of 3598 municipal employees is drawn from the comprehensive administrative registers. The study-design includes two outcome variables: length of working life at retirement in 2017 (years) and amount of final pension (Eur/month), and 40 socioeconomic explaining variables preceding retirement in 2005/2010–2016 (see Table 1).

Table 1 Yearly explaining variables used in the study with the retirement data.

Variable name	Years	Description
Unemp.	2010–2016	Unemployment spells
Occ. reh.	2010–2016	Occupational rehabilitation spells
Sickness	2010–2016	Long-term sick leaves
Accident	2010–2016	Occupational accident recovery spells
Wage	2005–2016	Salary wage in Euros

The distributions of both outcomes have long tails, both left and right, and no zeros. Explaining variables include measures (spells) of social security benefits available for all employees – unemployment, occupational rehabilitation, long-term sickness and occupational accident compensation (days in year), and wages (Eur/year). The distributions of explaining variables have long right tails and excess number of zeros, as some of the aforementioned benefits are rarely used. A note is in order on the fact that in practice zero yearly wages follow from absence from work for different reasons, such as, working in private sector or as self-employed, or drawing social security benefits. Wages and pensions are transformed using asinh transformation (Bellemare and Wichman 2020).

The results indicate that the adaptive technique yields more non-selected variables for several social security benefits compared to ordinary LAD-lasso (see Figures 5 and 6 and Table A1). From practical point of view and variable selection perspective, it is important to separate which covariates become truly zero and which are nearly so. From Figure 6, we can also see that benefits drawn near retirement (2017) tend to be mild (in terms of values of the estimates), yet meaningful estimates on both outcomes. Occupational rehabilitation and unemployment benefits have positive effect on work ability and working life, whereas long-term sickness indicates weakening work ability and risk on working life. On the other hand, drawing social security benefits in general indicates negative effect on final pension. Wages have positive effect on pension, which is as expected in earnings-related pension system. The effects of yearly wages on the length of working life are somewhat mixed.

Fig. 5 Multi-outcome ordinary LAD-lasso coefficient estimates for earnings-related disability pension (left) and length of working life (right). Filled circles refer to the value of the coefficient while triangles refer to zero-valued coefficients.

Fig. 6 Multi-outcome adaptive LAD-lasso coefficient estimates for earnings-related disability pension (left) and length of working life (right). Filled circles refer to the value of the coefficient while triangles refer to zero-valued coefficients.

The spike in the estimate for wages in 2005 on pension is a result of high correlation between consecutive wages (see Figure A1). This multicollinearity in regression analysis can generally lead to instability in the estimates, and to a problem where the explained variance is allocated arbitrarily among the correlated variables (Farrar and Glauber 1967).

5 Concluding remarks

We have shown that the multi-outcome adaptive LAD-lasso is a versatile and robust tool, capable of reducing bias from effect estimates as well as alleviating the influence from the problem of an excess of zero values in continuous covariates. The competitive performance of the adaptive version compared to the ordinary LAD-lasso has been illustrated with simulated and real-life examples.

The empirical data example illustrates the practical importance of multi-outcome modeling strategy with highly skewed outcome distributions (including outliers) and high share of zeros in covariates, which is common situation, for example, in social sciences analyzing complex socioeconomic phenomena. Researchers in other fields of science, with high-dimensional study-designs ($n\gg p$, $p\gg n$ or $q>1$) could certainly also benefit from robust estimates provided by multi-outcome adaptive LAD-lasso regression technique.

As future steps, we have considered multiple ways to improve the LAD-lasso modeling approach. First, it would be beneficial to develop a robust BIC criterion for a multi-outcome LAD-lasso context in order to take into account linear dependencies between outcomes. Second, as we demonstrated with the empirical data, the problem of multicollinearity in the explaining variables can lead to some unintuitive estimates. To tackle this problem, some interesting approaches have been developed, such as, fused lasso (Tibshirani et al. 2005) and clustering algorithm in regression via data-driven segmentation (CARDS) (Ke, Fan, and Wu 2015). Implementing these concepts into the LAD-lasso framework might be beneficial to applications with correlated explaining variables. Third, the inclusion of categorized variables in the adaptive LAD-lasso modeling approach would be intriguing in many applications, especially in the life-course studies, where several categorical factors exist in the data sets. This requires a careful implementation of coding scheme such that the re-coded variables are managed correctly in the adaptive algorithm (O’Grady and Medoff 1988). Last, the size of the data sets in many fields bring some challenges to the estimation. Naturally, a more scalable algorithm is often needed to perform the estimation in a reasonable time.

Disclosure statement

No potential competing interest was reported by the authors.

Additional information

Funding

This research is supported by the Infotech Oulu project funding.

Appendix A

Retirement analysis

Fig. A1 Correlation matrix of yearly wages in Section 4.

Table A1 Model coefficient estimates (cf. Section 4). Dropped variables are marked in italics.

	Ordinary LAD-lasso		Adaptive LAD-lasso
	asinh(Pension)	Working life	asinh(Pension)	Working life
Intercept	5.1958	6.6389	4.8417	6.6413
Unemp. 10	–0.0050	0.0005	–0.0045	0.0006
Unemp. 11	–0.0045	0.0010	–0.0048	0.0010
Unemp. 12	–0.0015	0.0002	–0.0000	0.0000
Unemp. 13	–0.0029	0.0019	–0.0033	0.0019
Unemp. 14	–0.0024	0.0011	–0.0007	0.0007
Unemp. 15	–0.0041	0.0000	–0.0047	0.0003
Unemp. 16	–0.0062	0.0015	–0.0060	0.0015
Occ. reh. 10	–0.0040	–0.0008	–0.0013	–0.0004
Occ. reh. 11	–0.0052	0.0002	–0.0047	0.0001
Occ. reh. 12	–0.0082	–0.0003	–0.0062	0.0004
Occ. reh. 13	0.0030	0.0008	0.0000	0.0001
Occ. reh. 14	0.0002	0.0005	0.0000	0.0000
Occ. reh. 15	–0.0069	0.0012	–0.0065	0.0014
Occ. reh. 16	–0.0073	0.0012	–0.0071	0.0012
Sickness 10	0.0000	–0.0000	0.0000	–0.0000
Sickness 11	–0.0000	–0.0000	–0.0000	–0.0000
Sickness 12	–0.0031	0.0045	–0.0000	0.0000
Sickness 13	–0.0000	–0.0001	0.0000	–0.0000
Sickness 14	0.0079	–0.0040	0.0012	–0.0009
Sickness 15	–0.0124	–0.0082	–0.0096	–0.0079
Sickness 16	0.0001	0.0001	–0.0000	0.0000
Accident 10	–0.0025	–0.0003	–0.0000	–0.0000
Accident 11	0.0000	0.0000	0.0000	0.0000
Accident 12	0.0086	0.0002	0.0000	–0.0000
Accident 13	0.0100	0.0051	0.0000	0.0000
Accident 14	–0.0147	0.0033	–0.0116	0.0037
Accident 15	0.0063	–0.0029	0.0000	–0.0000
Accident 16	0.0000	–0.0000	0.0000	–0.0000
Wage 05	1.3254	0.0024	1.3464	0.0031
Wage 06	0.1100	0.0120	0.0913	0.0147
Wage 07	0.1698	0.0057	0.1837	0.0058
Wage 08	0.0907	0.0025	0.0822	–0.0028
Wage 09	0.1613	–0.0002	0.1704	0.0044
Wage 10	0.1155	–0.0085	0.1224	–0.0127
Wage 11	0.0593	–0.0002	0.0692	0.0062
Wage 12	0.0278	0.0086	0.0016	0.0010
Wage 13	0.0516	0.0066	0.0445	0.0085
Wage 14	0.1583	0.0302	0.1950	0.0339
Wage 15	0.0028	0.0065	–0.0000	0.0000
Wage 16	0.0417	0.0187	0.0459	0.0232