As a subclass of non-negative data, semicontinuous data comprises a continuous distribution of positive values and a probability mass function of a nonnegligible number of zero values. We encounter semicontinuous data widely in various applications. Health care expenditures and medical costs are well-known applications (Duan et al., 1983; Liu, 2009; Mullahy, 1998; Smith et al., 2014; Tu and Zhou, 1999), where a considerable portion of healthy people incur no health care expenditures, while clinic or hospital visitors spend a wide range of medical costs. Other applications include substance abuse (Brown et al., 2005), alcohol consumption (Olsen and Schafer, 2001), and smoking behavior (Zhao et al., 2016).

The presence of two heterogeneous distributions of semicontinuous data makes ordinary least squares (OLS) estimation biased and inefficient. The two-part model (TPM) has gained popularity as an alternative to OLS. The TPM separately models the binary response (either zero or positive) and the mean response given that it is positive (Duan et al., 1983; Cragg, 1971). Recent TPM studies include deep learning-based feature importance (Zou et al., 2023), TPM quantile regression (Merlo et al., 2022), and multivariate TPM model (Frees et al., 2013). Similar to the TPM, Tweedie models can model semicontinuous data using a single unified distribution, called a Tweedie distribution (Dunn and Smyth, 2005; Kokonendji et al., 2021; Tweedie, 1984). Despite the widespread use of the TPM for semicontinuous data, researchers have paid relatively little attention to variable selection and prediction performance of the TPM. Variable selection of identifying a best subset of predictors is a backbone of predictive regression modeling, and it is increasingly important to deal with high dimensionality in the era of big data. Our study aims to explore the variable selection and prediction performance of the two-part model using selected regularized regression methods in high-dimensional data.

Penalized (or regularized) regressions are useful in high-dimensional data in which a few number of predictors may contribute to modeling the response. Penalized regressions simultaneously perform variable selection and parameter estimation. Popular penalized regression methods based on soft thresholdinging include the least angle shrinkage and selection operator (LASSO) (Tibshirani, 1996), elastic net (ENET) (Zou and Hastie, 2005), adaptive LASSO (Zou, 2006), group LASSO (Yuan and Lin, 2006), Dantzig selector (Candes and Tao, 2007), smoothly clipped absolute deviation (SCAD) (Fan and Li, 2001), and minimax concave penalty (MCP) (Zhang, 2010). Since recent variations of penalized regressions (Hao et al., 2018; Fan and Lv, 2008) mainly incorporate LASSO, ENET, SCAD, and MCP, we focus on these methods for the variable selection and prediction of TPM. We will compare these methods to traditional variable selection methods such as Akaike information criterion (AIC) and Bayesian information criterion (BIC) in simulation and empirical studies. The selected penalized regressions are estimated using the coordinate descent algorithm (CDA) was introduced (Wu and Lange, 2008; Breheny and Huang, 2011). The CDA optimizes a single parameter with others held fixed, cycling until the solution path stabilizes.

The TPM offers several intriguing characteristics. First, the two models for binary outcome and continuous positive response utilize the same predictor variables. However, the significance of individual predictors may differ between the two response models. Second, modeling the binary response uses a larger number of observations than the continuous positive response. In other words, not only is the entire sample size important, but the proportion of the positive values is important to modeling TPM. The conditional linear model for the continuous positive response may suffer from a smaller sample size. Last, the marginal mean of the response is a function of the probability that the responses are positive and the marginal mean of the positive responses. These characteristics deserve thorough investigation while performing variable selection and prediction in the TPM.

The simulation study probes the effects of various statistical scenarios in controlled settings. Our empirical study considers the prediction of crime incidents. In recent years, there have been significant advances in predictive modeling of crime incidents because of its societal implications and importance. Crime prediction has employed a broad spectrum of machine learning algorithms and data sources. These prediction algorithms encompass nonparameteric regression, support vector machine, and deep neural network, and the crime data include spatiotemporal data, social media data, and community-based data (Kang and Kang, 2017 and references therein). Our predictive modeling incorporates penalized regression using community-based data. The regression approach has advantages of identification and interpretation of the predictors associated with crime incidents.

The remaining sections are as follows. Section 2 briefly introduces the two-part model and its mean squared error. Section 3 describes regularized regression-based variable method and estimation. Section 4 describes the design and result of simulation study using the methods. Section 5 implements empirical data analysis using the community-based crime data. Last, Section 6 discusses the findings and related issues from both simulation study and empirics, concludes the study, and suggests further research areas.

Let y = (y₁, …, y_n)^T be an n × 1 response vector, where y_i, i = 1, …, n, are independently observed semicontinuous responses such that, for 0 < π_i < 1,

Let X = (x₁, …, x_p) be the n × p matrix of predictors, where x_j is a n × 1 vector of the j^th predictor for j = 1, …, p. The TPM requires two sets of parameter spaces, α = (α₁, …, α_p)^T and β = (β₁, …, β_p)^T. Let with cardinality , and with cardinality , where j is an element of the index set . Then and are subvectors of α and β, respectively. Also, and be submatrices of X. We also define two indicator functions associated with the response as follows:

A conventional two-part model of y_i is defined as

where h(y_i|y_i > 0, x_i) is any continuous density function for y_i > 0. In the TPM in (2.1), π_i = P(Y_i > 0) can be modeled by a parametric binary probability model such as logit or probit using predictors, that is, π_i = P(Y_i > 0|x_i). In this article, the logit model is used as follows:

The logistic regression in (2.2) is the first part of the TPM. Consequently, the positive values of y_i are modeled by the parameterization of h(y_i|y_i > 0, x_i) as follows:

where the error terms ɛ_i are independent and identically normally distributed with mean 0 and variance σ². The positive values of the response are often skewed right, so they are modeled through a log-transformation of Y. The log-normal assumption is widely applied, but it is often restrictive, which can be relaxed by the log-skew-normal assumption (Smith et al., 2014). This linear model is the second part of the TPM. The TPM of (2.2) and (2.3) is sometimes called the Bernoulli-log-normal regression model (Neelon et al., 2016). The two regression models of the TPM are independently fitted. We can see the independence of fitting the two models via their likelihood functions. The likelihood function of the two-part model for the random sample of n independent observations is

The likelihood function of the two-part model in (2.4) is decomposed into two multiplicative terms (Duan et al., 1983; Min and Agresti, 2002). The first term solely depends on the parameters, α, in (2.2), and the second term solely depends on the parameters, β, as in (2.3). This multiplicity allows to separately achieve the maximum likelihood estimate (MLE) of α and the MLE of β and σ for the

TPM. The log-likelihood function of TPM, denoted l(·) = log L(·), is given as

The expected value of the response variable, Y, is defined as

Plugging (2.2) and the expected value of (2.3) into (2.6) yields

As mentioned in the introduction, the expected value in (2.7) is a function of the probability that Y > 0 and the expected value of Y for y > 0. Furthermore, the independence of the likelihood functions of two part model leads to the predicted value of Y for the TPM as follows:

where α̂ and β̂ can be estimated by various estimation methods, which are described in the following section. The mean squared error (MSE) is defined as

which will be used to evaluate prediction performance in the following sections.

In this section, we will discuss how we select and estimate and , each of which are subvectors of α and β, respectively. Variable selection is a key analytical component in predictive regression analysis. Numerous variable selection techniques have been proposed, and this trend will remain steady in the data-centered era. Variable selection techniques include stepwise methods (forward, backward, and stepwise selection methods), pretesting method (univariate t-test), prediction-oriented criteria (minR², maxR², Mallow’s C_p, and adjusted R²), information criterion methods (AIC, BIC, and AIC for a small-sample correction (AIC_c)), and penalized regression methods (LASSO, Ridge, ENET, SCAD, and MCP). Limited space restricts us to cover selected methods that have been the most used in practice, which include four penalized regression methods (LASSO, ENET, SCAD, and MCP).

For the purpose of variable selection and prediction via penalized regression, we consider a penalized log-likelihood of TPM, denoted l_p(·), using (2.5) as follows:

where the nonnegative penalty functions are defined as

in which the tuning parameter λ > 0 controls a trade-off between the likelihood (or loss function) and the penalty. The tuning parameter can be determined by many different methods such as BIC and cross-validation. The shrinkage parameter γ determines the degree of shrinkage of parameters in the concave penalty methods, which is either fixed or estimated. A similar setting of the penalized likelihood in (3.1) was introduced for zero-inflated Poisson model (Tang et al., 2014). The structure of the two penalty functions, and , are the same except the parameters α and β. Hence, we discuss only in this article. A type of the penalty function is the bridge penalty function which has the following functional form

where the value of γ is related to famous regression methods. When γ = 2 in penalty (3.2), the penalized likelihood in (3.1) becomes ridge regression which was proposed to solve multicollinearity among predictors in the 1970s. Unfortunately, the ridge regression does not achieve variable selection because the parameter estimate could be close to zero but never achieves the zero value. When γ = 1, (3.1) is called LASSO regression, which simultaneously achieves variable selection and parameter estimation. ENET regression proposed by Zou and Hastie (2006) incorporated a convex combination of ridge and LASSO penalties to take advantage of both methods.

where δ = λ₁/(λ₁ + λ₂).

Although the LASSO regression has substantial advantages, this method suffers from biased estimate under a certain condition (Zou, 2006). The solution of the bridge penalty is only continuous if γ ≥ 1, and sparse only if γ = 1. When γ = 1, the LASSO solution is shifted by a constant value of λ. In order to simultaneously satisfy mathematical conditions for unbiasedness, sparsity, and continuity, non-convex penalty functions were introduced. We consider two important non-convex methods, SCAD (Fan and Li, 2001) and MCP (Zhang, 2010). The SCAD penalty function is given

for λ ≥ 0 and γ > 2. The MCP penalty is given:

for λ ≥ 0 and γ > 1.

The original LASSO solution employed quadratic programming in Tibshirani (1996). Efron et al. (2004) demonstrated that the LASSO solution is a variant of the Least Angle Regression (LARS) algorithm, which is a computationally efficient stagewise procedure based on piecewise linearity. The coordinate descent algorithm proposed by Wu and Lange (2008) is significantly faster than LARS. This algorithm is incorporated in glmnet and ncvreg packages in R, where the former is for LASSO-type estimators, and the latter is for non-convex penalty estimators. Nonconvexity bares a burden of numerical optimization, but the coordinate descent algorithm provides faster, stable numerical solutions.

The regularized regression method simultaneously achieves variable selection and parameter estimation, which are encompassed in the oracle property that an oracle estimator must hold asymptotic consistency in both variable selection and parameter estimation. Fan and Li (2001) and Zou (2006) posited that a good selection and estimation procedure should hold these two oracle properties. If λ_n/n → 0 and λn/n→∞, the LASSO estimator achieves a correct selection of the set of true non-zero parameters. However, Zou (2006) showed that when λn=O(n), . To overcome this issue, Zou (2006) introduced the adaptive LASSO (ALASSO) under the assumption that λn/n→0, which satisfies the oracle property. The SCAD and MCP penalty methods also achieve the oracle property.

Numerical optimization is a serious challenge in the regularized regression. Tibshirani (1996) used quadratic programming (QP) with a convex constraint as a special case of convex optimization to find the LASSO solution. The LARS algorithm is another technique to find the piecewise linear path, where the LARS algorithm is a homotopy method in the sense that the piecewise linear path is sequentially constructed. Meanwhile, SCAD and MCP encountered more serious challenge due to the nonconvexcity of penalties. The local linear approximation (LLA) algorithm using the LARS algorithm was proposed to find the solution for SCAD and MCP (Zou and Li, 2008). Regardless of the convexity or nonconvexity of the penalties, the CDA method was proposed as a fast and computationally efficient optimization method. The CDA method achieves the optimization of an objective function for a single parameter with fixing all other parameters and iteratively cycling through all parameters until convergence is achieved. The computational efficiency of CDA is O(np), which is even better than the linear regression with QR decomposition (O(np²)). Detailed rationale behind the CDA approach to MCP and SCAD is available in Breheny and Huang (2011). In this study, we incorporated the CDA method to compare the regularized regression methods.

The goal of parameter estimation is to achieve an (asymptotically) unbiased and consistent estimator. However, the goal of variable selection is more complex. The possible goals of variable selection include sensitivity, specificity, predictability, and selection consistency (Dziak et al., 2012). Sensitivity measures how a method accurately includes the set of predictors with non-zero coefficients. Specificity measures how a method accurately excludes the set of predictors with zero coefficients. Predictability measures a certain degree of prediction error. Selection consistency indicates how a set of selected predictors approaches to the true set of predictors with non-zero coefficients as n → ∞. As pointed out in many studies (Ng, 2013 and references therein), these four goals may not be achieved simultaneously in the sense that a high sensitivity (or specificity) does not guarantee a high predictability, and vice versa. Therefore, we compare sensitivity, specificity, and MSPE in the subsequent simulation study.

Our simulation study aimed to investigate the performance of selected methods described in the previous section under various situations. The selected methods included LASSO, ENET, SCAD, and MCP. If p = 100, 2^p ≈ 1.27 × 10³⁰ in which subset selection based on an information criterion is almost prohibited due to computational complexity. LASSO and ENET were implemented in the R glmnet package (Friedman et al., 2009) and SCAD and MCP were implemented in the R ncvreg package (Breheny, 2013).

The design of the simulation study for the two-part model considered the following issues. First, we considered two different sample sizes, n = 500, and 1000, for both the training and test data, which means that the ratio of both the data is one to one. The different sample sizes allowed to check the behavior of variable selection and prediction performance as sample size varied. Second, the length of the parameter space, or the number of predictors, are p = 20, 100, 1, 000 without the intercept. We consider two parameter spaces as follows:

The first parameter space comprises few strong coefficients, which was used in the seminal LASSO paper (Tibshirani, 1996) and many other studies (Fan and Li, 2005; Zou and Hastie, 2001). The second parameter space was introduced to investigate how the selected penalized regression methods perform for diverse coefficient values.

Third, the covariance structure among predictors is one of the most important factors affecting the variable selection performance. We considered the independent and AR(1) covariance structures as follows:

• Independence: ∑_x = σ²I_p.

• AR(1) correlation: σ²σ_{i j} = ρ^|i−j| for i = 1, …, p and j = 1, …, p.

Our simulation study used that σ = 3 and ρ = 0.5. Fourth, all predictors were standardized after being generated from a multivariate normal distribution, which was implemented in the R mvrnorm package. Fifth, in TPM, LM and GLM share the same data set where LM uses a less portion of the data. Hence, the proportion of zero values could affect the performance. In our case, we set 50–60 % data positive by controlling the intercept value in the simulation because a small portion of zeros or positive values can cause biased estimation in logistic regression. The shrinkage parameter (γ) in SCAD and MCP used the default value in the package, and the tuning parameter was estimated using a 10-fold cross-validation in the training data. All other options used the default values in the packages.

We focus on the two types of performance evaluation: Variable selection and prediction performances. These performance measures are closely related to the oracle properties. The variable selection performance is measured by β-sensitivity and β-specificity, which are defined as:

where is either or described in section 2, and for the estimated parameter space θ̂ using any selected method, and R indicates the number of iterations, which is 200 in our study. The β-sensitivity measures how accurately a selection method includes the predictors with non-zero coefficients in the model, and the β-specificity measures how accurately a selection method exclude the predictors with zero coefficients from the model.

The prediction performance for the GLM and LM is measured separately using accuracy, sensitivity, and specificity as metrics for the GLM, and the mean squared prediction error (MSPE) for the LM. The GLM metrics are used for the classification between zero and positive values. The mean absolute deviation serves as an alternative to the MSPE. These metrics are all evaluated in the testing data set. It is considered that the lower the MSPE, the better the prediction performance of the LM, while for the GLM, values of accuracy, sensitivity, and specificity closer to one indicates better classification performance.

Tables 1–4 reported the average and standard error of the performance measures obtained from 200 simulations to evaluate variable selection and prediction performance in the test data set. Table 1 presented the average values of β-sensitivity and β-specificity for the logistic regression (GLM) and linear regression (LM) for the parameter space P1. Table 2 presented the two MSPE values for the positive values and all values of y and the classification measures such as accuracy, sensitivity, and sprecificity for the binary responses in the GLM. Similarly, we presented the variable selection and prediction performance for the paramter space P2 in Tables 3 and 4, respectively. The performance results for AIC and BIC were only reported for the small number of predictors, p = 20.

As evident in Table 1, all the methods achieved better β-sensitivity average values although they showed various results in β-specificity for the parameter space P1, which contains strong coefficient values. The AR(1) covariance strucucture showed slightly lower performance than the independent covariance structure. In Table 2, the prediction results showed slightly different performance among the selected methods. For the average test MSPE, the LASSO and ENET methods outperformed the SCAD and MCP methods. However, the methods demonstrated no clear differences in classification between the zero and positive values. All the standard error values had no significant differences across the selected methods.

For the parameter space P2, which contains coefficients generated from the uniform distribution between 0.5 and 2, the selected methods overall presented underperformance in variable selection and prediction assessment, as can be seen Tables 3 and 4. Unlike the results of P1, P2 demonstrated that the LASSO and ENET methods outperform the SCAD and MCP methods in terms of β-sensitivity.

Our simulation results showed that penalized TPMs achieve better performance in both variable selection and prediction in high-dimensional data, regardless of covariance structures and sample sizes. The selected nonconvex methods such as SCAD and MCP exhibited better performance in β-specificity than the LASSO-type methods. The LASSO-type methods exhibited better performance in prediction than the nonconvex methods, regardless of the parameter spaces. Meanwhile, traditional variable selection methods based on AIC and BIC achieved better performance in β-specificity and MSPE than penalized methods for low-dimensional data.

Our empirical study for TPM was conducted using community-based crime data. Crime data can be collected in many ways. One common way is to collect crime data for each community such as city or town. Each community possesses its idiosyncratic characteristics with respect to demographics and socioeconomic status. Redmond and Baveja (2002) generated a comprehensive community-based crime data which is available at UC Irvine Machine Learning Repository(https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime). The data set contains 2,215 observations (communities), 124 predictors, and 18 response variables (9 different types of crimes with the original frequency and the frequency per 100,000 inhabitants) from multiple original data sources such as socio-economic data from the 1990 Census, law enforcement data from the 1990 Law Enforcement Management and Administrative Statistics (LEMAS), and crime data from the 1995 FBI Uniform Crime Report (UCR) Statistics.

The Census data include age-, race-, income-, family-, and house-related variables. The UCR data contained the original counts and the counts per 100,000 inhabitants for murder, rape, robbery, assault, burglary, larceny, auto theft, and arson. The LEMAS data collected policing-related data from state and local law enforcement agencies, including all those that employ 100 or more sworn police officers and a nationally representative sample of smaller agencies. Many communities with a small number of sworn officers had missing values for the predictors from the LEMAS data.

Among the 9 types of crime, murder, arson, and rape showed semicontinuous distributions, where these types of crime contained a considerable number of zero values over 10% of the total number of observations and demonstrated a skewed right distribution as the mean is much greater than the median. This empirical study focused on the response variable of the murder incidents per 100,000 habitants of which semicontinuous property is demonstrated in Figure 1. In Figure 1, the right-side figure illustrated the zero values and a bell-shaped curve of log-transformed positive values, which was modeled via the TPM. After removing the 23 LEMAS variables and any communities with significant amount of missing values, our final analytical data set consisted of one dependent variable and 101 predictors in 1,901 communities.

In order to check multicollinearity among 101 predictors, first we identified the perfect linearity (and hence, multicollinearity) between OwnOccQrange and OwnOccLowQuart/OwnOccHiQuart as well as RentQrange and RentLowQ/RentHighQ using the alias function for the lm function in R. Therefore, we removed two variables, OwnOccQrange and RentQrange. Using the variance influence factor (VIF), we examined multicollinearity among 99 predictors. The variance inflation factor (VIF) analysis shows that the 83 out of 99 predictors had the squared VIF value greater than 2, which indicates that the community-based crime data presented severe multicollinearity among predictors. In summary, the community-based crime data was characterized by skewed-right responses with semi-continuity and a fairly large number of predictor with multicollinearity. We identified a set of predictors via various variable selection methods and evaluated their prediction performance described in the previous sections. We split the whole data into the training and test data sets with 1 : 1 ratio, which resulted in a slightly different sample sizes for the linear regression in the training and test data sets.

In Table 5, we demonstrated the MSPE and the number of predictors selected by each methods. This result should be cautiously interpreted because a different sampling of training and test data may lead to a different result. Overall, the folded concave penalty methods of SCAD and MCP outperformed other methods with respect to the MSPE. The MSPE of MCP is 12% lower than the ENET one. AIC and ENET methods tend to select a higher number of predictors, and BIC and MCP tend to select a smaller number of predictors, which is consistent to the simulation results reported in the previous section.

In Figure 2, we presented two Venn diagrams of predictors selected by the four penalized methods from both GLM and LM in Table 5 to closely look at the patterns among the selected predictors. The LASSO predictors were mostly selected by the ENET method, and the MCP predictors were mostly selected by SCAD. The LASSO and ENET selected more predictors than SCAD and MCP in both models. The four penalized methods in both GLM and LM commonly selected seven predictors, although their compositions were different. The commonly selected variables from the logistic regression among four methods included population, blackPerCap, PersPerOwnOccHous, PersPerRentOccHous, MedNumBR, LandArea, and racePctWhite, and the commonly selected variables from the linear regression among four methods included PersPerFam, NumKidsBornNeverMar, NumImmig, PersPerOwnOccHous, PopDens, racePctWhite, and PctWorkMomYoungKids where the two predictors, PersPerOwnOccHous(mean persons per owner occupied household) and racePctWhite (percentage of population that is Caucasian) were common for both models. For further information, refer to the dictionary of these predictors at the UCI Machine Learning Repository mentioned above.

In this study, we investigated penalized regression-based variable selection methods for two-part models. We conducted simulation studies under diverse statistical assumptions and an empirical study using community-based crime data. Our analytical results demonstrated that penalized TPMs achieve better performance in both variable selection and prediction in high-dimensional data, regardless of covariance structures and sample sizes. Moreover, the LASSO-type methods such as LASSO and ENET outperformed the nonconvex methods such as SCAD and MCP in mean squared error. In simulation studies, for a small number of predictors, for example, p = 20, traditional variable selection methods based on information criteria achieved better performance in β-specificity and MSE than penalized methods.

Variable selection in the TPM is affected by several unique features in addition to conventional matters, such as high dimensionality, multicollinearity, covariance structure, and sparsity. These TPM-specific features include the same pool of predictors for both models in (2.2) and (2.3), a smaller sample size for modeling positive values, and the mean squared error as a function of the probability of a positive value and the marginal mean of the positive values. Our simulation was designed such that positive values were 50–60% of the total sample size. Simulation results showed that the sensitivity of LM was as good as those of GLM while the specificity of LM was similar to or worse than those of GLM. As can be seen in (2.8), the probability of a positive value results in a smaller MSE for TPM compared to that of the linear regression for the positive values.

Additionally, the simulation study also showed that the best variable selection may not be associated with the least prediction error. The trade-off between consistent variable selection and efficient prediction is well addressed in Ng (2013) and references therein. The convex penalty methods and the folded concave penalty methods show different behaviors in sensitivity and specificity. When the coefficients are quite different from zero as in Table 1, both penalty methods performed well for sensitivity. On the other hand, when coefficients were generated from Uniform[0.25, 2] as in Tables 3 and 4, the two methods exhibited similar sensitivities, or the convex methods marginally performed better. The folded concave methods outperformed the convex penalty methods in the specificity performance regardless of parameter spaces and covariance structures. This result is partly because the convex penalty methods are inclined to include more predictors with non-zero coefficient estimates, as explained in Breheny and Huang (2011).

Our current study can be extended in several directions. First, as we only considered the selected number of methods, it is worthwhile to consider recently developed variable selection methods especially for high dimensional data such as interaction selection (Hao et al., 2018) and screening methods (Fan and Lv, 2008; Pan et al., 2019; Tibshirani et al., 2012). In particular, screening out insignificant predictors is expected to help improve the specificity issue of penalized regression methods. Second, the positive values of a semicontinuous variable might follow the log-skewed normal distribution, necessiating a generalized regression beyond a linear regression under the normality assumption. It is desirable to consider regularized GLM with gamma or skewed lognormal distribution for various degrees of positive skewness.