Mixture of Censored Partially Linear Experts Models for Symmetric or Asymmetric Data with or without Outlying Observations: A novel approach with two-piece scale-mixture of normal distributions

Introduction

The concept of finite mixture of regressions (FMR), first put forward by Quandt (1972), offers a way to model data that arise from multiple hidden subgroups. In this approach, each latent group is assumed to follow its own regression relationship, allowing the overall model to capture unobserved structural differences in the population. The literature on FMR has grown considerably, with substantial methodological developments reported by Jiang and Tanner [20]; McLachlan and Peel [25]; Neykov et al. (2007); Garcia-Escudero et al. [12]; Bai et al. [3]; Bashir and Carter [4]; Hunter and Young [16]; Yao et al. (2014); Song et al. (2014); Ma and Genton (2021); Zarei et al. (2022); Oh and Seo (2024).

A well-recognized drawback of traditional FMR models is their reliance on the assumption that covariates play no role in determining cluster membership ([15]). When this independence assumption is not appropriate, the resulting regression–based clustering may be unreliable. To overcome this issue, Jacobs et al. [19] introduced the mixture of linear experts (MoE) framework, in which the mixing proportions are allowed to vary as functions of the covariates. Since then, the MoE model has undergone numerous extensions. Nguyen and McLachlan (2016) explored the use of Laplace errors, while Chamroukhi [8, 9] examined Student-t and skew-t error structures. Murphy and Murphy (2020) proposed a more parsimonious MoE architecture aimed at improving estimation efficiency. Mirfarah et al. [27] incorporated scale mixtures of normal to model the error distribution. More recently, Oh and Seo (2024) introduced a tailored MoE formulation, adding to this expanding line of research.

For data arising from a single homogeneous population, Engle et al. [11] proposed the partially linear model, which blends a parametric linear component with a nonparametric function. In this framework, the response variable $Y$ is modeled as:

where $\mathbf{x} \in\mathbb{R}^{p}$ is a $p$-dimensional covariate vector, $t \in\mathbb{R}$ is an additional predictor entering through an unknown smooth function $g(\cdot)$, and $\varepsilon$ is an error term. We rewrite the nonparametric elements in model (1) as

where $\mathbf{n}$ is a $q \times1$ incidence vector with the $j$-th element equal to the indicator function $I(t=t_m^0)$ for $m=1,\ldots,q$, and $\mathbf{g}=(g(t_1^0),\ldots,g(t_q^0))$ with $t_1^0,\ldots,t_q^0$ being the distinct and ordered values of $t$ (see e.g., Bazrafkan et al. [5]). Building on this idea, Wu and Liu (2017); Skhosana et al. (2023); Hwang et al. [17] extended the finite mixture of regressions framework to incorporate partially linear structures within heterogeneous populations.

In many practical settings, researchers must contend with censored and often missing observations. Such censoring commonly arises due to detection limits of measuring devices or restrictions built into study designs. For example, clinical trials, economic data collections, and environmental monitoring programs frequently exhibit left- or right-censoring when the actual values lie beyond the range that can be reliably recorded. In mixture models, where data come from several hidden subpopulations censored data (partially observed values) and missing data (unobserved values) make it harder to identify each component. Methods like the EM algorithm handle this by treating the unobserved parts as latent variables and adjusting the likelihood accordingly. In this work, we introduce a more adaptable and resilient framework: a mixture of censored partially linear expert models built upon the two-piece scale mixture of normal (TP–SMN) distributions (Maleki and Mahmoudi [24]). The TP–SMN family extends traditional SMN models by allowing asymmetry, making it a powerful tool for capturing data features such as skewness and heavy-tailed behavior. Compared with symmetric scale mixture of normal (SMN; Andrews and Mallows [2]) or asymmetric scale mixture of skew-normal (SMSN; Branco and Dey [6]) formulations, TP–SMN distributions offer enhanced flexibility, which is especially valuable when modelling censored outcomes that may include outliers or asymmetric patterns. One of the main strengths of the proposed finite mixture of censored partially linear expert models built on TP–SMN distributions is that all conditional expectations needed in the E–step of the ECME algorithm can be obtained in closed form. This removes the need for MCMC–based approximations, which are often computationally intensive and prone to numerical inaccuracy, and consequently yields more stable and precise maximum likelihood estimates compared with other mixture models employing skewed distributions. Overall, the proposed approach is well suited for applications where censoring and departures from normality occur simultaneously. Extensive simulation studies and analyses of real datasets illustrate the effectiveness and practicality of the proposed semiparametric mixture censored partially linear expert framework. The remainder of this paper is organized as follows: § 2 reviews relevant properties of the TP–SMN family and the associated components needed for an EM–type algorithm to obtain maximum penalized likelihood (MPL) estimates. § 3 introduces the semiparametric mixture of censored partially linear expert (MoCPLE) models based on the TP–SMN family and outlines the EM–type algorithm for parameter estimation along with the observed information matrix. § 4 presents simulation results and real data analyses to evaluate the performance of the proposed method. Finally, § 5 concludes with some remarks and future research directions.

The TP-SMN class of distributions (An overview)

In this section, we review some properties of the TP–SMN family of distributions.The well known SMN family introduced by Andrews and Mallows [2] has pdf

where $\phi(\cdot;\mu,\sigma^{2})$ represents the density of the normal distribution $\mathcal{N}(\mu,\sigma^{2})$, and $H(\cdot;\nu)$ plays the role of an arbitrary cumulative distribution function (cdf) containing $\nu$, which includes several members and is denoted by $X \sim\mathcal{SMN}(\mu,\sigma,\nu)$.

The TP-SMN is a rich class of distributions constructed based on the SMN family that covers symmetric/asymmetric and light/heavy-tailed distributions in terms of their pdf and cdf. For $x \in\mathbb{R}$, this family is represented as

where $0<\gamma<1$ is the slant parameter, $f_{\mathrm{SMN}}(\cdot;\cdot)$ and $F_{\mathrm{SMN}}(\cdot;\cdot)$ denote the pdf and cdf of the SMN distribution, respectively. A random variable with pdf and cdf given in (4) and (5) is denoted by $X\sim\mathcal{TP}$-$\mathrm{SMN}(\mu,\sigma,\gamma,\nu)$, which has the following stochastic representation:

where $Z\sim N(0,\sigma^2)$, the scalar random variable $U\sim\mathcal{H}(\cdot;\nu)$ and the shape random variable $F$, with the following probability mass function (pmf):

The advantage of the above representation of the binary random variable $F$ becomes evident in the construction of the penalized log-likelihood function of the model.

To illustrate the richness of the TP-SMN class of distributions, note that various members of this family can be obtained by choosing different distributions for the scale-mixing random variable $U$ in (6). Specifically, we have:

• Two-piece Normal (TP–N): $U=1$ with probability one.

• Two-piece t (TP–T): $U\sim\operatorname{Gamma}(\nu/2,\nu/2)$, $\nu>0$.

• Two-piece Slash (TP–SL): $U\sim\operatorname{Beta}(\nu,1)$, $\nu>0$.

• Two-piece Contaminated-Normal (TP–CN): $P(U=u)=\nu I_{(u=\tau)}+(1-\nu)I_{(u=1)}$; $0<\nu<1$, $0<\delta\le1$.

Finally, the proposed TP-SMN random variable $X$ has mean and variance given by:

where $b=\sqrt{2/\pi}k_1$, $\Delta=\sigma(1-2\gamma)$, $c_r=\gamma^{r+1}+(-1)^r(1-\gamma)^{r+1}$ and $k_r=E(U^{-r/2})$, $r\in\mathbb{Z}$.

The essential conditional expectations for the non-truncated and interval-truncated TP-SMN family used for parameter estimation are provided in the Appendix, with further details available in Maleki and Mahmoudi [24].

The TP–SMN Mixture of Censored Partially Linear Expert Model

Our focus is on semiparametric mixtures of partially linear models, where the responses may be observed or censored and the errors follow independent and identically distributed (i.i.d.) TP-SMN distributions. So we introduce a semiparametric mixture of censored partially linear experts model which hereafter so called TP-SMN-MoCPLE model.

Model specification

Mixtures of partially linear experts (MoPLE) models are constructed by treating each expert network of the MoE model, for $k = 1,\ldots,K$ sub-populations, as a partial linear model (2). The expert components are defined using TP-SMN distributions and specifying the component means via polynomial regression and the mixing proportions via a multinomial logistic function, the resulting pdf can be expressed as:

where $\mathbf{Y}=(Y_1,\ldots,Y_n)^{\mathsf{T}}$ is the vector of response variables, $\mathbf{x}_i\in\mathbb{R}^p$, and $\mathbf{r}_i=(1,r_1,\ldots,r_{d-1})^{\mathsf{T}}\in\mathbb{R}^d$ are the vector of explanatory, additional predictor and covariate variables corresponding to $Y_i$ for $i=1,\ldots,n$, $f_k(\cdot\mid\cdot)$ is the TP-SMN pdf (4), and for the gating parameters $\boldsymbol{\tau}=(\boldsymbol{\tau}_1^{\mathsf{T}},\ldots,\boldsymbol{\tau}_{K-1}^{\mathsf{T}})^{\mathsf{T}}$ with $\boldsymbol{\tau}_k=(\tau_{k0},\ldots,\tau_{k(d-1)})^{\mathsf{T}}$, and

for which the model’s parameter set is $\boldsymbol{\Theta}=(\boldsymbol{\Theta}_1^{\mathsf{T}},\ldots,\boldsymbol{\Theta}_K^{\mathsf{T}},\boldsymbol{\tau})^{\mathsf{T}}$ such that $\boldsymbol{\Theta}_k=(\boldsymbol{\beta}_k^{\mathsf{T}},\mathbf{g}_k^{\mathsf{T}},\sigma_k,\gamma_k,\nu_k)$, $k=1,\ldots,K$.

In Mixture of Experts (MoE) models, the gating network utilizes a Softmax function to assign weights to experts. Since the Softmax output is invariant to the addition of a constant vector $\mathbf{c}$ to all expert parameters (i.e., $\boldsymbol{\tau}_k\to\boldsymbol{\tau}_k+\mathbf{c}$ yields identical probabilities), the model parameters become non-identifiable, leading to infinite equivalent solutions. To resolve this redundancy and ensure a unique parameter estimation during optimization, we constrain the parameter vector of the first expert to zero, setting $\boldsymbol{\tau}_1=0$ as a reference baseline. Furthermore, along with this identifiability constraint on the gating parameters, assuming that the smooth functions within each expert are continuous, the covariate $\mathbf{x}_i$ is non-constant and not a deterministic function of $t$, and the parameter sets across experts are distinct, ensures the overall identifiability of the proposed model.

In the proposed TP-SMN-MoPLE models, we assume also that the response variables are censored observations and called it by TP-SMN mixture of censored partially linear experts (MoCPLE) models. In other words, some of the responses are subject to censoring, which may take the form of interval-, left-, or right-censoring. Thus, let the available response variable $Y_i$ be presented as the joint variables $(\mathbf{W}_i,\rho_i)$, where $\mathbf{W}_i$ represents the uncensored reading $(\mathbf{W}_i=Y_{i0})$ or interval censoring $(\mathbf{W}_i=(C_{i1},C_{i2}))$, and $\rho_i$ is the censoring indicator as:

Note that in this setting, if $C_{i1}=-\infty$ (or $C_{i2}=+\infty$), the left-censoring (or right-censoring) is occurred, if $C_{i1}=-\infty$ and $C_{i2}=+\infty$, the missing is occurred, and in the case $-\infty<C_{i1}<C_{i2}<+\infty$ the interval-censoring realization is observed. We establish our methodology based on the interval-censoring scheme; however, the left/right-censoring and missing schemes are also investigated in the simulation and real-data analyses. Under this framework, $\mathbf{Y}$ is partitioned into observed values and censored observations. Therefore, $\mathbf{Y}$ is treated as a latent variable due to its partial observability. Given these assumptions and using Remark 1, the penalized log-likelihood function of the MoCPLE model can be written as

where $\boldsymbol{\mu}_{ik}=\mathbf{x}_i^{\mathsf{T}}\boldsymbol{\beta}_k+\mathbf{n}_i^{\mathsf{T}}\mathbf{g}_k$, $\boldsymbol{\lambda}=(\lambda_1,\ldots,\lambda_K)^{\mathsf{T}}$, $\mathbf{w}=(w_1,\ldots,w_n)^{\mathsf{T}}$ and $\boldsymbol{\rho}=(\rho_1,\ldots,\rho_n)^{\mathsf{T}}$ denote the realizations of $\mathbf{W}=(W_1,\ldots,W_n)^{\mathsf{T}}$ and $\boldsymbol{\rho}=(\rho_1,\ldots,\rho_n)^{\mathsf{T}}$, respectively. The term $J(\mathbf{g}_k)$ is the penalty function over $\mathbf{g}_k(\cdot)$ and $\alpha_k$ is a smoothing parameter that controls the trade-off between goodness-of-fit, measured by large values of $\ell_p(\boldsymbol{\Theta},\lambda\mid\boldsymbol{w},\boldsymbol{\rho})$, and smoothness, measured by small values of $J(g_k)$. Following Ibacache-Pulgar and Paula [18], we use a natural cubic spline to estimate $g_k(\cdot)$. The general form of the penalty function is

where $g_k''(t)$ denotes the second derivative of $g_k(t)$ on $[a,b]$ which contains the values $t_1^0,\ldots,t_q^0$. According to Theorem 2.1 of Green and Silverman [14], $\int_a^b [g_k''(t)]^2\,dt=\boldsymbol{g}_k^{\mathsf{T}}\boldsymbol{L}\boldsymbol{g}_k$, for which $\boldsymbol{L}\in\mathbb{R}^{q\times q}$ is a non-negative definite matrix depending only on the knot differences.

The maximum penalized likelihood (MPL) estimates are obtained by maximizing the penalized log-likelihood. Because the log-likelihood function in (8) is analytically intractable, closed-form MPL estimates are unavailable, and numerical optimization is required. By exploiting the stochastic representation in (6) together with the hierarchical representation of the TP-SMN distribution, we develop an expectation conditional maximization either (ECME; Liu and Rubin [22]) algorithm to estimate parameters from the TP-SMN-MoCPLE models. The ECME algorithm is an extension of the expectation conditional maximization (ECM; Meng and Rubin [26]) that not only inherits its stable properties (e.g., monotone convergence and implementation simplicity) but also can be faster than ECM. The iterative ECME algorithm replaces some CM-steps of the ECM with CML-steps that maximize the corresponding constrained log-likelihood function instead.

MPL estimation of the MoCPLE models parameters via an EM-Algorithm

Beginning with pdf (7), we define the component label vector

such that $Z_{ik}=1$ if and only if $Z_i^*=k$. Under this formulation, we have

and under the above approach, the latent random vector $\boldsymbol{Z}_i$, $i=1,\ldots,n$, has the following multinomial pmf:

Let $\boldsymbol{y}_{\mathrm{com}}=(\boldsymbol{w}^{\mathsf{T}},\boldsymbol{\rho}^{\mathsf{T}},\boldsymbol{\gamma}^{\mathsf{T}},\boldsymbol{r}^{\mathsf{T}},\boldsymbol{u}^{\mathsf{T}},\boldsymbol{f}^{\mathsf{T}},\boldsymbol{z}^{\mathsf{T}})^{\mathsf{T}}$, then $\ell_{\mathrm{cp}}(\boldsymbol{\Theta})=\sum_{i=1}^{n}\ell_{icp}(\boldsymbol{\Theta})$, where

Let the superscript $(m)$ indicates the estimate of the related parameter at $m$ stage of the algorithm. In the E-step of the algorithm, we must obtain the so-called Q-function $Q(\boldsymbol{\Theta}|\boldsymbol{\Theta}^{(m)}) = \mathbb{E}\left[\ell_{cp}(\boldsymbol{\Theta}) \mid\boldsymbol{w},\boldsymbol{\rho},\boldsymbol{\Theta}^{(m)}\right]$ is given by:

where

In what follows, we derive expressions for conditional expectations under both uncensored and censored scenarios,

Uncensored case ($\rho_i=0$):

Censored case ($\rho_i=1$):

where $\widehat{\mu}_{ik}^{(m)}=\boldsymbol{x}_i^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_k^{(m)}+\boldsymbol{n}_i^{\mathsf{T}}\widehat{\boldsymbol{g}}_k^{(m)}$.

Computing $\widehat{\xi}_{j,ik}^{(m)}$, $j=0,1,2$, and $\widehat{s}_{ik}^{(m)}$ for uncensored and censored responses are obtained via Lemmas 1, 2 and 3 in the Appendix § 7.

• E-step: Knowing $\boldsymbol{\Theta}=\widehat{\boldsymbol{\Theta}}^{(m)}$ and using the above calculations, we compute the corresponding conditional expectations. Then we maximize $Q(\boldsymbol{\Theta}|\widehat{\boldsymbol{\Theta}}^{(m)})$ in the EM approach by solving the first-order condition

where $\nabla$ denotes the gradient operator, yielding the following updates.

• CM-step: Update $\boldsymbol{\Theta}^{(m)}$ (denoted by $\boldsymbol{\Theta}^{(m+1)}$) by maximizing $Q(\boldsymbol{\Theta}|\boldsymbol{\Theta}^{(m)})$ over $\boldsymbol{\Theta}$. For $k=1,\ldots,K$ we obtain:

• CM-step 1: According to the proposition 2 of Nguyen and McLachlan (2016), the update of $\boldsymbol{\tau}_k$ can be update as

• CM-step 2:

where $\widehat{\mu}_{ik}^{(m+1)}=\mathbf{x}_i^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_k^{(m+1)}+\mathbf{n}_i^{\mathsf{T}}\widehat{\boldsymbol{\gamma}}_k^{(m+1)}$.

• CML-step: Update $\boldsymbol{\nu}=(\nu_1^{\mathsf{T}},\ldots,\nu_K^{\mathsf{T}})^{\mathsf{T}}$ by maximizing the actual penalized log-likelihood:

where $\ell_{\rho}(\cdot\mid\mathbf{w},\rho)$ is defined in (8) and $\boldsymbol{\Theta}_{(-\boldsymbol{\nu})}^{(m+1)}$ denotes the updated parameter vector without $\boldsymbol{\nu}$ at iteration $m+1$.

To start the algorithm, initial values are required, so First, the response data are partitioned into $K$ groups using the $k$-means clustering algorithm. The proportion of observations belonging to each group is considered as the initial mixing weights (component probabilities). Then, the initial values of the regression coefficients for each component are obtained by applying the least squares method within each corresponding partitioned group. Finally, based on the computed residuals, the initial values of the parameters of the TP-SMN family are obtained using the procedure proposed by Maleki and Mahmoudi [24]. The algorithm iterates between the E-step and the CML-step until convergence is achieved. Convergence is assessed by monitoring the difference between two successive evaluations of the log-likelihood function. Specifically, the algorithm stops when $|\ell^{(m+1)}-\ell_{\infty}^{(m+1)}|\leq\epsilon$, where $\epsilon$ is a predetermined tolerance level. Here, $\ell^{(m+1)}$ denotes the value of the log-likelihood evaluated at $\boldsymbol{\Theta}^{(m+1)}$, and $\ell_{\infty}^{(m+1)}$ is the asymptotic estimate of the log-likelihood at the $(m+1)$th iteration, computed using Aitken’s acceleration as

where the Aitken acceleration at iteration $m$ is defined as $c^{(m)}=\frac{\ell^{(m+1)}-\ell^{(m)}}{\ell^{(m)}-\ell^{(m-1)}}$.

To address the label switching problem and ensure parameter identifiability, we impose an ordering constraint on the gating parameters. Since the parameter vectors $\boldsymbol{\tau}_k$ are multidimensional, we define a fixed reference vector $\mathbf{a}$ to project them onto scalars. A common choice is the vector of ones, $\mathbf{a}=(1,\ldots,1)^{\mathsf{T}}$, in which case $\boldsymbol{\tau}_k^{\mathsf{T}}\mathbf{a}$ represents the sum of the elements of $\boldsymbol{\tau}_k$. We then constrain the parameters such that $\boldsymbol{\tau}_1=\mathbf{0}$ and $\boldsymbol{\tau}_2^{\mathsf{T}}\mathbf{a}<\cdots<\boldsymbol{\tau}_K^{\mathsf{T}}\mathbf{a}$. This constraint eliminates the symmetry of the likelihood function and guarantees a unique solution during the estimation process.

Effective degrees of freedom and model selection

In additive models, the Akaike information criterion (AIC; Akaike [1]) can be applied to select an appropriate $\alpha_k$, $k=1,\ldots,K$. The AIC for MoCPLE models are given by

Let $\mathbf{Q}=\mathbf{N}^{\mathsf{T}}\mathbf{N}$ and $\mathbf{E}_k=\hat{\sigma}_k^{-2}\mathbf{Q}^{-1/2}\mathbf{L}\mathbf{Q}^{-1/2}$, then, from Eilers and Marx [10], we may express

where $e_{kj}$ is the $j$th eigenvalue of $\mathbf{E}_k$ and $\mathbf{I}_q$ is the identity matrix of order $q$. Thus, it follows that $\operatorname{df}(\lambda_k)$ is inversely proportional to $\lambda_k$. Therefore, one has a total of $p+\psi+\operatorname{df}(\lambda_k)$ parameters to be estimated, for which $p$ is the dimension of the regression parameters, and $\psi$ is the number of parameters of the proposed TP-SMN distributions in the $k$-th component. The Akaike information criterion is useful for selecting an appropriate model for the given data with adequate sample size, and applied in minimizing the function given by

where $\ell_p(\hat{\boldsymbol{\Theta}},\lambda_k\mid\mathbf{w},\boldsymbol{\rho})$ is the component-wise of the penalized log-likelihood function (8) available at $\hat{\boldsymbol{\Theta}}$ for a fixed $\lambda_k$. We employ a grid search strategy over the $\lambda$-space to determine the optimal value of the penalty parameter $\boldsymbol{\lambda}=(\lambda_1,\ldots,\lambda_K)^{\mathsf{T}}$.

Standard error approximation

In this section, we obtain the standard errors of the MPL estimates of the parameters of the proposed TP-SMN-CPLE models. Under the usual regularity conditions, we have that the asymptotic covariance matrix of $\boldsymbol{\Theta}$ can be approximated by the inverse of the empirical information matrix defined as

where $\mathbf{s}(y_i\mid\boldsymbol{\Theta})=\frac{\partial Q_i(\boldsymbol{\Theta}\mid\boldsymbol{\Theta}^{(r)})}{\partial\boldsymbol{\Theta}}$, and $Q_i(\boldsymbol{\Theta}\mid\boldsymbol{\Theta}^{(r)})$ is the summand in (12). Substituting $\boldsymbol{\Theta}$ by the MPL estimates $\hat{\boldsymbol{\Theta}}$ in $\mathbf{I}_0$, conclude its approximation in the form of $\hat{\mathbf{I}}_0=\sum_{i=1}^{n}\hat{\mathbf{s}}_i\hat{\mathbf{s}}_i^{\mathsf{T}}$, for which $\hat{\mathbf{s}}_i=\mathbf{s}(y_i\mid\hat{\boldsymbol{\Theta}})$ is an individual score vector given by

with explicit expressions for the elements of $\hat{\mathbf{s}}_i$ given by

Note the block part of $\hat{s}_{(i,\nu_k)}$ has some integrals which calculated numerically. Standard errors of $\hat{\boldsymbol{\Theta}}$ are extracted from the square root of the diagonal elements of the inverse of $\hat{\mathbf{I}}_0$.

Simulation studies based on Monte-Carlo

In this section some Monte-Carlo simulation studies are considered to verify the performances of the MPL estimates. The proposed ECME algorithm has been implemented in the R software (R Core Team (2025)) version 4.5.2 with a core i7 760 processor 2.8 GHz.

In particular, we conduct three simulation experiments using the proposed censored model (3), generating censored samples with censoring levels of $p=10\%$, $15\%$, and $30\%$. It is noteworthy that the results across these censoring levels were comparable; therefore, for brevity, we report the results for a representative censoring level in each study.

To explain the mechanism of fulfills the non-informative condition of censored sampling, we consider the following algorithm recommended by Gomez et al. [13].

Recovery and asymptotic properties of the MPL estimates

To examine the recovery of model’s parameters and asymptotic properties of the MPL parameter estimates of the TP-SMN-MoCPLE model obtained via the ECME algorithm, a simulation study with $K=2$ and following presumed parameters is performed,

$\nu_1=3$ and $\nu_2=4$ for the TP-T-MoCPLE and TP-SL-MoCPLE models, and $(\nu_1,\nu_2)^{\top}=(0.7,0.7)^{\top}$ and $(\nu_1,\nu_2)^{\top}=(0.3,0.3)^{\top}$ for the TP-CN-MoCPLE model. The non-parametric $g_1(t)=\cos(0.3\pi t)$ and $g_2(t)=\sin(0.16\pi t)$ for $t\in[0,12]$, for the first and second components, receptively. We set up $\mathbf{x}_i=(1,x_{i1})^{\top}$ and $\mathbf{r}_i=(1,r_{i1},r_{i2})^{\top}$, such that $x_{i1}\sim U(2,4)$, $r_{i1}\sim U(-2,1)$ and $r_{i2}\sim U(-1,1)$, for sample sizes $n=100,250,400,700,1000$, and then using the Remark 3 with $c=1$ and level of censoring $p=15\%$ on the simulated datasets, the ECME algorithm described in § 3.2 is preformed to carry out the MPL parameter estimates. To investigate the parameters recovery with 1000 replications, we consider the bias (BIAS) and the mean squared error (MSE):

where $\widehat{\theta}_j$ denotes the estimate of a specific model’s parameter at the $j$th replication.

Figure 1 displays the BIAS and MSE plots of the MPL estimates of the model’s parameters $\beta_k$, $\sigma_k$ and $\gamma_k$ for $k=1,2$, of the TP-N-MoCPLE, TP-T-MoCPLE, TP-SL-MoCPLE and TP-CN-MoCPLE models for the censored datasets. It can be observed that the MPL estimates of the parameters have very small (around zero) BIAS for all sample sizes, and, as $n$ increases the MSE of MPL estimates tend to zero. In the case of TP-T-MoCPLE model, the plots in Figure 2 also shows the descending trends and the variability of the nonparametric function estimates increases, and diminishes with increasing sample size. Similar results were obtained for the non-parametric parts of the other TP–SMN–CPLMs.These behaviours indicate the empirically consistency of the our MPL estimates of the model parameters via the ECME algorithm.

Model selection performance in the presence of outliers and asymmetry

Optimal number of experts $K$ in the MoE models is one the main challenges. In this simulation study we generate 1000 sample of size $n=300$ from the MoCPLE model that error terms on each expert follow the skew-t (SMSN; Branco and Dey [6]), which is a skew and heavy-tailed distribution with similar properties of the TP-T distribution in the class of the TP-SMN family. We call the proposed model as ST-MoCPLE model, which same as some members of our models, cover the censorious, asymmetry and outliers. But note that ST-MoCPLE model has not conditional expectations and so the MPL estimates are not exist, but we expect our MoCPLE models could detect the true number of experts.

We consider ST-MoCPLE model with three number of experts $K=3$ and following presumed parameters:

with exponential, harmonic and inverse power non-parametric parts of experts for the first, second and third experts given by:

for $t\in[0,10]$, $\boldsymbol{x}_i=(1,x_{i1})^{\mathsf{T}}$ and $\boldsymbol{r}_i=(1,r_{i1})^{\mathsf{T}}$, such that $x_{i1}\sim U(2,4)$, $r_{i1}\sim U(1,2)$, with $\nu_1=\nu_2=\nu_3=4$ and other remaining parameters are chose such that have three following scenarios:

• Strongly separated experts (S.S.E.):

• Moderately separated experts (M.S.E.):

• Weakly separated experts (W.S.E.):

where $\alpha_k$, $k = 1, 2, 3$ are the skewness parameters of skew-$t$ distribution on each expert (Example of generated samples of the above three scenarios are shown in Figure 3), and then using the Remark 3 with $c = 1$ and level of censoring $p = 10\%$ on the simulated datasets.

In the following analyses, the number of mixture components (experts), $K$, is assumed to be unknown. Accordingly, we fitted the TP-N-MoCPLE, TP-T-MoCPLE, TP-SL-MoCPLE, and TP-CN-MoCPLE models to the simulated data. This fitting was performed for each model assuming $K = 1, 2, 3, 4$, and $5$ experts in each replication. The results obtained from the AIC criterion, which indicate the percentage of correct specification for the number of experts among the different replications, are presented in Table 1. As shown in Table 1, the TP-T-MoCPLE, TP-SL-MoCPLE, and TP-CN-MoCPLE models outperform the TP-N-MoCPLE model in determining the number of components, likely due to the heavy-tailed nature of the generated data. Notably, the TP-T-MoCPLE model achieves a better fit than the other models, particularly when the optimal number of experts is three ($K = 3$).

Separating Level	Model	$K = 1$	$K = 2$	$K = 3$	$K = 4$	$K = 5$
S.S.E.	TP-N-MoCPLE	2.3%	4.1%	2.7%	1.2%	1.8%
	TP-T-MoCPLE	67.2%	60.4%	67.3%	58.4%	47.9%
	TP-SL-MoCPLE	20.3%	20.5%	21.8%	18.9%	27.6%
	TP-CN-MoCPLE	10.2%	15.0%	8.2%	21.5%	22.7%
M.S.E.	TP-N-MoCPLE	4.3%	3.8%	1.9%	1.1%	0.9%
	TP-T-MoCPLE	58.7%	59.2%	62.8%	39.1%	40.6%
	TP-SL-MoCPLE	24.1%	20.5%	12.4%	22.6%	20.8%
	TP-CN-MoCPLE	12.9%	16.5%	22.9%	37.2%	37.7%
W.S.E.	TP-N-MoCPLE	1.8%	2.3%	1.5%	2.1%	1.9%
	TP-T-MoCPLE	34.1%	30.8%	69.3%	34.1%	39.2%
	TP-SL-MoCPLE	28.4%	30.1%	30.6%	33.3%	30.1%
	TP-CN-MoCPLE	35.7%	37.1%	5.2%	30.5%	28.8%

Real datasets analysis

In this section to asses the performance the proposed model and estimates we consider two real applications. The first dataset is New York Air Quality Measurements, well known as airquality which is meteorological data and available in an R package datasets, and second is Labor Force Participation Data, well-known as PSID1976 dataset which is economical data, and available in an R package AER.

New York Air Quality Measurements Data

For the first empirical analysis, we utilized the well-known airquality dataset, which contains daily air quality measurements recorded in New York City from May 1st to September 30th, 1973. This dataset, originally provided by the New York State Department of Conservation (for ozone data) and the National Weather Service, consists of a data frame with 153 observations (Chambers et al. [7]). After removing missing values of auxiliaries and gating variables, our analysis is conducted on 111 complete observations. We consider the following variables:

• Response variable (**Ozone; $y$): Mean ozone concentration in parts per billion (ppb), measured between 1300 and 1500 hours at Roosevelt Island.

• First auxiliary variable (**Solar.R; $x_1$): Solar radiation in Langleys, measured in the frequency band of 4000–7700 Angstroms between 0800 and 1200 hours at Central Park.

• Second auxiliary variable (**Month; $x_2$): The month of observation (May to September), included as a categorical covariate in the linear component.

• First gating variable (**Wind; $r_1$): Average wind speed in miles per hour (mph), recorded at 0700 and 1000 hours at LaGuardia Airport.

• Functional variable (**Temp; $t$): Maximum daily temperature in degrees Fahrenheit, recorded at LaGuardia Airport, used to construct the nonparametric basis matrix.

Parameter $\downarrow$	TP-N-MoCPLE Estimates	TP-N-MoCPLE SE	TP-T-MoCPLE Estimates	TP-T-MoCPLE SE	TP-SL-MoCPLE Estimates	TP-SL-MoCPLE SE	TP-CN-MoCPLE Estimates	TP-CN-MoCPLE SE
$\beta_{10}$	14.8654	0.7343	12.2424	0.6104	12.8632	0.7014	13.0142	0.9132
$\beta_{11}$	0.1201	0.0432	0.0480	0.0065	0.0491	0.0036	0.0501	0.0044
$\beta_{12}$	0.1067	0.0945	0.5709	0.0213	0.5811	0.0196	0.6002	0.0201
$\beta_{20}$	-1.0533	0.1456	-0.6205	0.0229	-0.6511	0.0184	-0.6801	0.0252
$\beta_{21}$	0.4711	0.0565	0.2201	0.0321	0.2196	0.0252	0.2297	0.0364
$\beta_{22}$	7.1028	0.9231	5.6859	0.8524	5.7532	0.9432	5.7604	0.8603
$\sigma_1$	20.2011	1.9309	16.2606	1.0763	16.7109	1.1201	16.7301	1.1172
$\sigma_2$	39.7212	3.3028	34.8594	2.2010	35.0172	2.2654	35.2191	2.4021
$\nu_1$	–	–	14.7762	0.7036	14.5673	0.8001	0.8389	0.0021
$\nu_2$	–	–	4.0576	0.1453	4.2011	0.1351	0.9001	0.0073
$\delta_1$	–	–	–	–	–	–	0.8243	0.0312
$\delta_2$	–	–	–	–	–	–	0.8637	0.0211
$\gamma_1$	0.4548	0.0076	0.4985	0.0024	0.4891	0.0032	0.4893	0.0031
$\gamma_2$	0.3683	0.0089	0.4759	0.0029	0.4702	0.0037	0.4731	0.0040
$\tau_{20}$	6.8746	2.2191	5.2849	0.6222	5.1873	0.4563	5.5631	0.6101
$\tau_{21}$	-1.263	0.0431	-0.7127	0.0283	-0.7990	0.0131	-0.7488	0.0236
AIC	1119.203		1055.603		1053.803		1061.203
BIC	1314.289		1258.818		1257.018		1267.127
RI	0.5123		0.5214		0.5323		0.5118
JCI	0.3676		0.3847		0.4029		0.3713

The trends observed in the pairwise plots in Figure 4, visualized via nonparametric regression curves (LOESS smoothers), indicate that a purely linear assumption may be inadequate to capture the complex relationships, particularly the nonlinear behavior of temperature and wind speed. To address this, we consider our semi-parametric framework. Specifically, we define the linear component using Solar.R and Month, while the nonlinear effects are modeled through the functional variable (Temp) and the gating variable (Wind). In the next step, we apply the TP-SMN-MoCPLE models with several experts to estimate the optimal model structure under an interval-censoring scheme.

In our modelling approach, we employed auxiliary regression variables comprising Solar.R and Month to define the linear component of the model. Furthermore, the gating expert variables selected for the model were Wind and the functional basis derived from Temp. These specific variables were chosen due to their established influence on air quality dynamics and their frequent use in related atmospheric studies. This dataset, along with the selected variables, provides a comprehensive basis for investigating the relationships between meteorological factors and air pollution levels.

Using the procedure described in Remark 3 with c=1 and an interval censoring level of 10% (including 20% of the censored observations assigned to wider intervals to simulate potential outliers), we fitted the TP-SMN-MoCPLE models to the data, considering components ranging from $K=1$ to $K=4$. Based on the AIC criterion, the two-expert ($K=2$) model was selected as the optimal fit. To the best of our knowledge, Mixture of Experts models based on semi-parametric frameworks and non-normal distributions have not yet been introduced in the literature, precluding a direct comparative analysis with existing studies. However, the incorporation of both skew and symmetric distributions within our proposed model provides a comprehensive representation of various data behaviors, thereby ensuring the model’s flexibility and robustness.

Table 2 presents the MPL estimates for the four considered models. The penalty coefficients $(\hat{\lambda}_1,\hat{\lambda}_2)$ for both experts were set to 0.01 and 0.02, using the $\lambda$-grid approach. The substantial disparity between the intercepts points to different underlying distributions for the two experts. Both solar radiation (Solar.R) and month of observation (Month) have positive influence on the mean ozone in both experts; however, the influence of solar radiation is more strong than that of month of observation, and influence of auxiliaries are more effective in the second expert. Regarding the shape parameters, the gamma estimates reveal slight skewness (near symmetry) for the both experts. Additionally, the small values of degrees of freedom on each experts (especially on the second) characterized by heavy tails and the presence of outliers. The coefficients in the gating function also demonstrate that average wind speed (Wind) has positive influence on expert selection. Finally, based on the AIC and BIC (Schwarz et al. (1978)) values recorded in the Table 3, the TP-T-MoCPLE and TP-SL-MoCPLE models demonstrate a superior fit compared to the other models. Figure 5 (left) illustrates the non-parametric functions $g_k$, $k=1,2$ for both experts, plotted against the response observations, and Figure 5 (right) displays the observed response values plotted against the estimated fitted values. The close alignment of the points along a straight line on Figure 5 (right), indicates that the model provides a good fit to the data.

Parameter $\downarrow$	TP-N-MoCPLE Estimates	TP-N-MoCPLE SE	TP-T-MoCPLE Estimates	TP-T-MoCPLE SE	TP-SL-MoCPLE Estimates	TP-SL-MoCPLE SE	TP-CN-MoCPLE Estimates	TP-CN-MoCPLE SE
$\beta_{10}$	15.3431	0.9873	13.8067	0.8112	13.4463	0.9038	13.9983	0.9432
$\beta_{11}$	-0.0412	0.0023	-0.0012	0.0003	-0.0014	0.0005	-0.0013	0.0003
$\beta_{12}$	0.9082	0.0887	0.5812	0.0232	0.5814	0.0306	0.5902	0.0311
$\beta_{13}$	0.1067	0.0945	-0.0602	0.0054	-0.0644	0.0065	-0.0612	0.0069
$\beta_{20}$	9.9109	0.8575	8.8633	0.2429	8.9018	0.2871	-0.6801	0.0252
$\beta_{21}$	0.4711	0.0565	-0.0019	0.0004	-0.0021	0.0004	-0.0019	0.0005
$\beta_{22}$	7.1028	0.9231	0.2907	0.0231	0.3010	0.0306	0.2958	0.0149
$\beta_{23}$	7.1028	0.9231	0.0028	0.0006	0.0030	0.0007	0.0029	0.0007
$\beta_{30}$	7.3918	0.9001	5.1371	0.6739	5.2101	0.6881	5.3010	0.7781
$\beta_{31}$	-0.0722	0.0146	-0.0011	0.0003	-0.0012	0.0005	-0.0014	0.0006
$\beta_{32}$	0.5093	0.0882	0.1690	0.0332	0.1705	0.0367	0.1686	0.0300
$\beta_{33}$	-0.0542	0.0655	-0.0080	0.0005	-0.0082	0.0004	-0.0089	0.0004
$\sigma_1$	18.3029	1.0817	5.4143	0.5387	5.4396	0.5611	5.4146	0.6048
$\sigma_2$	6.6355	0.9983	2.3106	0.1872	2.4018	0.1998	2.2827	0.2018
$\sigma_3$	7.0948	1.0008	2.7178	0.1099	2.8018	0.1222	2.6902	0.1100
$\nu_1$	–	–	2.0500	0.0653	2.2201	0.0809	0.0162	0.0043
$\nu_2$	–	–	6.9094	0.1908	7.2011	0.2019	0.0210	0.0088
$\nu_3$	–	–	12.9973	0.2201	13.2031	0.2353	0.3001	0.0104
$\delta_1$	–	–	–	–	–	–	0.0621	0.0041
$\delta_2$	–	–	–	–	–	–	0.1031	0.0113
$\delta_3$	–	–	–	–	–	–	0.2109	0.0344
$\gamma_1$	0.3309	0.0276	0.2738	0.0108	0.2409	0.0099	0.2693	0.0110
$\gamma_2$	0.3209	0.0201	0.4133	0.0138	0.4042	0.0221	0.3939	0.0202
$\gamma_3$	0.7033	0.0332	0.6267	0.0200	0.6301	0.0248	0.6132	0.0313
$\tau_{20}$	6.0983	1.4533	7.1261	0.5984	7.8736	0.7001	7.5847	0.5666
$\tau_{21}$	-0.8981	0.0431	-0.0001	0.0000	-0.0002	0.0000	-0.0001	0.0000
$\tau_{22}$	-1.263	0.0431	-0.1217	0.0109	-0.1345	0.0145	-0.1298	0.0121
$\tau_{30}$	20.1981	1.8716	14.7412	0.9817	15.0911	1.0081	15.0918	0.9987
$\tau_{31}$	-0.8711	0.0601	-0.0005	0.0000	-0.0008	0.0000	-0.0010	0.0000
$\tau_{32}$	-1.5409	0.1675	-0.3375	0.0312	-0.4012	0.0406	-0.3987	0.0298
AIC	3756.218		3538.333		3600.380		3658.380
BIC	4422.084		4218.071		4280.118		4351.990
RI	0.5123		0.5214		0.5323		0.5118
JCI	0.3676		0.3847		0.4029		0.3713

Labor Force Participation Data

To further assess the applicability of our proposed methods, we employ the PSID1976 dataset, a widely used cross-section dataset originating from the 1976 Panel Study of Income Dynamics (PSID). This dataset focuses on labor force participation and related socio-economic factors for individuals in 1975. It comprises 753 observations offering a rich context for analyzing factors influencing economic decisions, that for our analysis we choose following variables:

• Response variable (hwage; $y$): Husband’s average hourly wage in 1975 dollars. This is often used as the dependent variable in seminal studies using this dataset.

• First auxiliary variable (hhours; $x_1$): Husband’s hours of work in 1975.

• Second auxiliary variable (heducation; $x_2$): Husband’s education in years.

• Second auxiliary variable (hage; $x_3$): Husband’s age in years.

• First gating variable (fincome; $r_1$): Family income in 1975 dollars, used for constructing property income.

• Second gating variable (unemp; $r_2$): Unemployment rate in county of residence, in percentage points. (This is taken from bracketed ranges.)

• Functional variable (experience; $t$): Actual years of wife’s previous labor market experience.

This dataset, also known as the ”wage rates” of 753 married women dataset Mroz (1987); Karlsson and Laitila [21], is particularly valuable for analyzing labor supply decisions. The trends observed in the pairwise plots in Figure 6, visualized via nonparametric regression curves (LOESS smoothers), indicate that a purely linear assumption may be inadequate to capture the complex relationships between the observed variables. We assume the husband’s average hourly wage depends on their hours of work, education, age and experience, and using the Remark 3 with $c=1$ and an interval censoring level of 10% (including 20% of the censored observations assigned to wider intervals

to simulate potential outliers). Therefore, we fit our proposed TP-SMN-MoCPLE models for several possible number of experts. Optimal number of experts $K = 3$ were obtained by AIC criteria between possible values $K = 1, 2, 3, 4$. Also using the $\lambda$-grid approach, the penalty coefficients $(\hat{\lambda}_1,\hat{\lambda}_2,\hat{\lambda}_3)$ were set to $(0.05,0.05,0.05)$, and MPL estimates of proposed models are provided in the Table 3.

Based on the AIC and BIC criteria, the TP-SMN-MoCPLE model demonstrates a superior fit to the data compared to competing models. The results indicate that education has a significant positive influence on the wages of husbands. This effect is most pronounced in the first expert, followed by the second, and finally the third, in a decreasing order of magnitude. Furthermore, the first expert exhibits pronounced asymmetry ($\hat{\gamma}_1 = 0.27$) and heavy-tailed behavior ($\hat{\nu}_1 = 2.05$). In contrast, the second and third experts display weak asymmetry ($\hat{\gamma}_2 = 0.41$ and $\hat{\gamma}_3 = 0.62$) but remain relatively heavy-tailed ($\hat{\nu}_2 = 7$ and $\hat{\nu}_3 = 13$) (see also right part of Figure 7).

The left part of Figure 7 illustrates the non-parametric density estimates of the fitted values for each of the three experts. The distinct separation and minimal overlap between these density curves provide clear visual evidence of a mixture structure comprising three distinct experts. Moreover, the shapes of these densities effectively capture the underlying probabilistic characteristics of each component, visually confirming the presence of skewness and heavy-tailed behaviour within the distributions of the individual experts.

Figure 8 displays the scatter plots of the response variable against the covariates—hours, education, and age—where the observations are categorized into three distinct clusters based on the fitted model. To visualize the underlying relationship within each group, linear regression lines are superimposed on the plots. These regression lines demonstrate appropriate trends within each cluster, showing consistency with the estimated experts. The clear separation of the clusters and the distinct trajectories of the fitted lines across the different covariates provide compelling visual evidence that the proposed three-component mixture model effectively captures the heterogeneity of the data. This structure confirms that the clustering derived from the model is well-defined and meaningful.

Conclusions and discussions

This paper introduces a novel robust class of Mixture of Partially Linear Experts (MoCPLE) models tailored for censored data, grounded in the family of two-piece scale-mixture of normal distributions. The proposed TP-SMN-MoCPLE framework extends the classical Mixture of Experts (MoE) paradigm, effectively addressing complex data challenges such as non-linearity, asymmetry, heavy-tailed disturbances, outliers, and censoring mechanisms. This model is highly general and encompasses various existing frameworks as special cases, including symmetric and asymmetric classical MoE models, finite mixtures of Partially Linear Models (PLM), Mixture Regression Models (MRM), and models designed for censored data. Furthermore, the incorporation of covariate-dependent gating functions facilitates the automatic selection of experts, thereby enhancing classification performance. By leveraging the TP–SMN family, the proposed model accommodates both light- and heavy-tailed distributions as well as skewness; notably, the zero-skewness constraint naturally recovers the well-known standard SMN family. A key advantage of this class is the availability of closed-form expressions for conditional expectations, which significantly streamlines the definition and estimation of the model. Notably, prior to this work, censored models based on skewed distributions possessing closed-form conditional expectations and allowing for exact Maximum Likelihood Estimators (MLE) without resorting to Monte Carlo approximations were non-existent. This contrasts with the present study, which offers this distinct computational advantage. The proposed methodology has been implemented in R software, and the corresponding computational code is available from the authors upon request.

Conditional expectations of the TP-SMN distributions

Uncensored cases

For the uncensored data $y$, we have $\rho= 0$. Therefore, we have the following results:

Censored cases

In the censored cases, we have $\rho=1$. First, we define the expectation functions

which are needed at follow, and for each special cases of the TP-SMN class of distributions are given by:

• TP-N:

• TP-T:

where $F_{PVII}(\cdot;\nu,\delta)$ denotes the cdf of Pearson type VII distribution,

• TP-SL:

where $\Gamma(a,b)=\int_0^\infty e^{-t}t^{a-1}dt$ is the incomplete gamma function and $F_{SL}(\cdot)$ is the cdf of the slash distribution,

• TP-CN:

then for both cases, we have:

Also for case 1, we have $\hat{s}=E[S\mid Y\in\mathcal{A}]=1$ and for case 2, $\hat{s}=E[S\mid Y\in\mathcal{A}]=-1$.