Moment-based Estimation of Mixtures of Regression Models

TL;DR

This paper introduces a moment-based method for estimating parameters in finite mixture of regression models, requiring minimal assumptions and demonstrating high accuracy through simulations and real data application.

Contribution

It develops unbiased, assumption-light estimators for mixture regression models, expanding the toolkit for flexible modeling without distributional constraints.

Findings

Estimators are unbiased and consistent.

Method performs well in simulations.

Application to wine data demonstrates practical utility.

Abstract

Finite mixtures of regression models provide a flexible modeling framework for many phenomena. Using moment-based estimation of the regression parameters, we develop unbiased estimators with a minimum of assumptions on the mixture components. In particular, only the average regression model for one of the components in the mixture model is needed and no requirements on the distributions. The consistency and asymptotic distribution of the estimators is derived and the proposed method is validated through a series of simulation studies and is shown to be highly accurate. We illustrate the use of the moment-based mixture of regression models with an application to wine quality data.

Tables

Table 1. Table 1: Results from simulation study of the moment mixture model based on 1000 simulations. The model provides unbiased estimates and the estimted SE match the empirical SEs. 95% CP is the 95% coverage probability and the efficiency is the relative efficiency of the SE for the moment mixture model relative to the SE of the optimal model.

	N	True value	Estimate	Emp. SE	Est. SE	95% CP	Rel. efficiency
Gaussian mixture
$\beta$	300	1.0	1.004	0.331	0.323	0.939	3.96
p	300	0.5	0.481	0.197	0.174	0.850
$\beta$	2000	1.0	1.002	0.120	0.121	0.948	3.83
p	2000	0.5	0.497	0.066	0.065	0.948
Zero-inflated Gaussian
$\beta$	300	1.0	0.977	0.315	0.307	0.943	3.76
p	300	0.5	0.506	0.181	0.156	0.870
$\beta$	2000	1.0	0.995	0.117	0.118	0.956	3.73
p	2000	0.5	0.502	0.059	0.056	0.957
Exponential-Gaussian mixture
$\beta$	300	1.0	0.998	0.441	0.394	0.929	4.83
p	300	0.5	0.476	0.231	0.204	0.712
$\beta$	2000	1.0	1.006	0.153	0.154	0.953	4.87
p	2000	0.5	0.495	0.085	0.080	0.956
Zero-inflated exponential
$\beta$	300	1.0	0.996	0.404	0.375	0.942	4.59
p	300	0.5	0.490	0.212	0.185	0.773
$\beta$	2000	1.0	0.999	0.154	0.151	0.946	4.78
p	2000	0.5	0.500	0.080	0.073	0.947

Equations

$$Y=P{Y}_1+(1-P)Y_2,$$

$$E(Y|X)=p\cdot ({\mu }_1+{\beta }^{T}X)+(1-p)\cdot {\mu }_2,$$

$$E(Y^2|X)=p\cdot (({\mu }_1+{\beta }^{T}X{)}^2+{\sigma }_1^2)+(1-p)\cdot ({\mu }_2^2+{\sigma }_2^2).$$

$$E(Y_i^2|X_i)=\tilde{\alpha }+{\lambda }_2\cdot {\eta }_i+{\lambda }_3\cdot {\eta }_i^2,$$

$${\hat{\lambda }}_1=A(m)({\tilde{X}}^T\tilde{X}{)}^{-1}{\tilde{X}}^{T}Y,$$

$$\sqrt{n}({\hat{\lambda }}_1-{\lambda }_1)=A(m)(\frac{1}{n}{\tilde{X}}^T\tilde{X}{)}^{-1}\frac{1}{\sqrt{n}}\sum _{i=1}^n{\tilde{X}}_i\{Y_i-E(Y_i|{\tilde{X}}_i)\}=\frac{1}{\sqrt{n}}\sum _{i=1}^n{\epsilon }_i+o_P(1)$$

$${\epsilon }_i=A(m)E({\tilde{X}}_i{\tilde{X}}_i^T{)}^{-1}{\tilde{X}}_i\{Y_i-E(Y_i|{\tilde{X}}_i)\}.$$

$$E(Y_i^2|X_i)=\tilde{\alpha }+{\lambda }_2\cdot {\eta }_i+{\lambda }_3\cdot {\eta }_i^2$$

$${\mathbb{P}}_n^w(f)=\frac{{\sum }_{i=1}^n{f}_i{w}_i}{{\sum }_{i=1}^n{w}_i}.$$

$$a_1^n({\lambda }_1)={\mathbb{P}}_n^w({\eta }^4)-{\mathbb{P}}_n^w({\eta }^2{)}^2,$$

$$a_2^n({\lambda }_1)={\mathbb{P}}_n^w({\eta }^3)-{\mathbb{P}}_n^w(\eta ){\mathbb{P}}_n^w({\eta }^2),$$

$$a_3^n({\lambda }_1)={\mathbb{P}}_n^w({\eta }^2)-{\mathbb{P}}_n^w(\eta {)}^2,$$

$$b_1^n({\lambda }_1)={\mathbb{P}}_n^w(\eta Y^2)-{\mathbb{P}}_n^w(\eta ){\mathbb{P}}_n^w(Y^2),$$

$$b_2^n({\lambda }_1)={\mathbb{P}}_n^w({\eta }^2{Y}^2)-{\mathbb{P}}_n^w({\eta }^2){\mathbb{P}}_n^w(Y^2).$$

$${\lambda }_2^n({\lambda }_1)=\frac{a_1^n({\lambda }_1)\cdot b_1^n({\lambda }_1)-a_2^n({\lambda }_1)\cdot b_2^n({\lambda }_1)}{a_1^n({\lambda }_1)\cdot a_3^n({\lambda }_1)-a_2^n({\lambda }_1{)}^2},$$

$${\lambda }_3^n({\lambda }_1)=\frac{a_3^n({\lambda }_1)\cdot b_2^n({\lambda }_1)-a_2^n({\lambda }_1)\cdot b_1^n({\lambda }_1)}{a_1^n({\lambda }_1)\cdot a_3^n({\lambda }_1)-a_2^n({\lambda }_1{)}^2}.$$

$${\hat{\lambda }}_2={\lambda }_2^n({\hat{\lambda }}_1),$$

$${\hat{\lambda }}_3={\lambda }_3^n({\hat{\lambda }}_1).$$

$$g({\lambda }_1)=\frac{\mathbf{P}\{f\cdot w\}}{\mathbf{P}w}$$

$$D_{{\lambda }_1}g({\lambda }_1)=(\mathbf{P}w{)}^{-1}P\{D_{{\lambda }_1}(f\cdot w)\}-(\mathbf{P}w{)}^{-2}\mathbf{P}\{f\cdot w\}\mathbf{P}(D_{{\lambda }_1}w).$$

$$\sqrt{n}\{{\mathbb{P}}_n^{\hat{w}}(\hat{f})-{\mathbb{P}}_n^w(f)\}=\sqrt{n}\{{\hat{\lambda }}_1-{\lambda }_1{\}}^T{D}_{{\lambda }_1}g({\lambda }_1)+o_P(1).$$

$$P\{f\cdot w{\}}^2<\infty .$$

$$\sqrt{n}\big{\{}\mathbb{P}^{w}_{n}(f)-\frac{\mathbf{P}\{f\cdot w\}}{\mathbf{P}w}\big{\}}=\{\mathbf{P}w\}^{-1}\mathbb{G}_{n}(\{f-\frac{\mathbf{P}(f\cdot w)}{\mathbf{P}w}\}\cdot w)+o_{P}(1),$$

$$\sqrt{n}({\hat{\lambda }}_2-{\lambda }_2)=\sqrt{n}\{{\hat{\lambda }}_1-{\lambda }_1{\}}^T{D}_{{\lambda }_1}{\lambda }_2({\lambda }_1)+\{g_3({\lambda }_1)\mathbf{P}w{\}}^{-1}{\mathbf{G}}_{\mathbf{n}}(w\{a_1({\lambda }_1)f^1-a_2({\lambda }_1)f^2\})+o_P(1),$$

$$\sqrt{n}({\hat{\lambda }}_3-{\lambda }_3)=\sqrt{n}\{{\hat{\lambda }}_1-{\lambda }_1{\}}^T{D}_{{\lambda }_1}{\lambda }_3({\lambda }_1)+\{g_3({\lambda }_1)\mathbf{P}w{\}}^{-1}{\mathbf{G}}_{\mathbf{n}}(w\{a_3({\lambda }_1)f^2-a_2({\lambda }_1)f^1\})+o_P(1),$$

$$f_i^1=({\eta }_i-\frac{\mathbf{P}(\eta \cdot w)}{\mathbf{P}w})\{Y_i^2-E(Y_i^2|X_i)\},$$

$$f_i^2=({\eta }_i^2-\frac{\mathbf{P}({\eta }^2\cdot w)}{\mathbf{P}w})\{Y_i^2-E(Y_i^2|X_i)\}.$$

$$\sqrt{n}(\hat{\beta }-\beta )=\frac{1}{\sqrt{n}}\sum _i^n{\xi }_i+o_P(1),$$

$${\xi }_i=[{\lambda }_3\cdot I_m+{\lambda }_1\{D_{{\lambda }_1}{\lambda }_3({\lambda }_1){\}}^T]{\epsilon }_i+{\lambda }_1\cdot \{g_3({\lambda }_1\mathbf{P}w{\}}^{-1}\cdot w_i\cdot \{a_1({\lambda }_1)f_i^1-a_2({\lambda }_1)f_i^2\}.$$

$$\frac{1}{n}\sum _{i=1}^n{\xi }_i{\xi }_i^T\stackrel{P}{\to }\text{var}({\xi }_i)$$

$$v_i=v({\eta }_i)=\frac{1}{(1+{\eta }_i^2)}$$

$$w_i=w({\eta }_i)=\frac{1}{(1+{\eta }_i^4)}$$

$$\text{Var}(\hat{\beta })=\frac{{\sigma }_1^2}{(p\cdot N-1)\cdot \text{Var}(X_i)}$$

$$g_1^n({\lambda }_1)=a_1^n({\lambda }_1)\cdot b_1^n({\lambda }_1)-a_2^n({\lambda }_1)\cdot b_2^n({\lambda }_1),$$

Full text

1 Introduction

The class of finite mixtures of regression models provides a flexible approach to model a wide range of phenomena and to handle non-standard data analysis problems such as excess zeros or heterogeneity (McLachlan and Peel, 2000). Mixture models can also be used to accommodate data that is ”contaminated” due to poor quality data, laboratory errors, or situations where data originates from multiple sources. In particular, zero-inflated regression models, and hurdle models can be considered special cases of the class of finite mixture of regression models with two components and both of these types of models see frequent use (Lambert, 1992; Kwagyan and Apprey, 2016).

Finite mixture regression models are relevant for analyzing many problems such as sudden-infant-death-syndrome (Dalrymple and others, 2003), HIV-risk reduction trials (Hu and others, 2011), medical care (Deb and Trivedi, 1997). However, problems relevant for finite mixtures are especially common in genomics where for example population admixture (certain subgroups of the population do not segregate a phenotype-influencing mutation), or in the presence of gene-environment interactions, where the effect of some genes are never triggered because the individual is living under specific environmental conditions.

Other recent approaches of using finite mixture models in genetics include the paper by Xu and others (2015) where they compare different model types (standard parametric, non-parametric, zero-inflated and hurdle models) to microbiome data in order to assess and infer the best model for these types of data. They find that mixture models in general fit the data best but that the choice of parametric model (Poisson or negative-binomial) can have substantial impact on the results and on the convergence on the estimation algorithm.

Typically, parameters in mixture regression models are estimated by specifying a fully parametric model for each of the components. This leads to efficient estimates provided that the parametric models are specified correctly. In this article we will consider the class of finite mixture regressions where the predictors are influencing exactly one of the regression model components. Within this framework we wish to estimate the association between the predictors and the outcome. We place no restriction on the distribution of the components but only require that the relationship between the predictors and the outcome can be modeled as a linear model. Moment-based estimators allow us to obtain unbiased estimates of the regression parameters with a minimum of assumptions. In particular, we do not require a specification of the distribution, and we only need to specify the average regression model for one of the components in the mixture model.

The recent paper by Kong and others (2017) uses an approach for identifying interaction that most closely resembles the approach in this manuscript. They also model the first and second moments but they are focused on interaction in particular and do not use a moment-based method for estimation.

The manuscript is organized as follows: The next section explains the problem, derives the moment-based estimators, and proves the large-sample properties relevant for inference. In section 3, we present properties of the proposed method for moment-based mixture of regression models and use simulations to compare the proposed moment-based mixture regression model to the analogous Gaussian mixture model. Finally, we apply the moment-based mixture of regression models to a dataset of wine quality before we discuss the findings along with possible extensions. The approaches presented in this manuscript are available in the R package mommix which can be found on github at www.github.com/ekstroem/mommix.

2 Methods

Let $Y$ be given as a mixture of two distributions, that is;

[TABLE]

where $\text{Pr}(P=1)=1-\text{Pr}(P=0)=p$ for some $p\in(0,1)$ and $Y_{1}$ and $Y_{2}$ are independent stochastic variables from the two underlying distributions. We are interested in modeling the relationship between a set of $m$ predictors represented by the design matrix $X$ and $Y_{1}$ and can think of $Y_{2}$ as a contamination of the response.

For the distribution of interest we assume that $E(Y_{1}|X)=\mu_{1}+\beta^{T}X$, and $V(Y_{1}|X)=\sigma_{1}^{2}$, but for the other component we only require that $Y_{2}\perp X$ and that it has well-defined mean and variance $E(Y_{2})=\mu_{2}$, and $V(Y_{2})=\sigma_{2}^{2}$.

The mean of $Y$ is

[TABLE]

where $\beta\in\mathbb{R}^{m}$ is the parameter vector of interest. Simple algebra yields that the second order moment is

[TABLE]

Note that the formula for the second order moment has a form corresponding to a multiple linear regression model with a mean that can be reparameterized as $2\cdot\mu_{1}\cdot p\cdot(\beta^{T}X)+p\cdot(\beta^{T}X)^{2}+\tilde{\alpha}$. Thus, if $X$ is univariate we can obtain estimates of $2\cdot\mu_{1}\cdot p\cdot\beta$ and $p\cdot\beta^{2}$ directly from least-squares regression of the second order moment and estimate $p\cdot\beta$ from least-squares regression of the first moment. These estimates can be combined to extract individual estimates of $\mu_{1}$, $\beta$, and $p$. In the following we formalize and extend this idea.

Let $\lambda_{1}=p\cdot\beta$ and set $\eta_{i}=\lambda_{1}^{T}X_{i}$ (the linear predictor), $\lambda_{2}=2\cdot\mu_{1}$, and $\lambda_{3}=p^{-1}$. Then equation (2) may be rewritten as

[TABLE]

where the intercept $\tilde{\alpha}$ contains contributions from both $Y_{1}$ and $Y_{2}$ but holds no direct information about $X$.

From this parameterization we have that $\beta=\lambda_{1}\cdot\lambda_{3}$, $\mu_{1}=\lambda_{2}/2$, and $p=\frac{1}{\lambda_{3}}$. Thus, if we can estimate $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ then we can also estimate the mixture proportion and all mean parameters related to the distribution of interest.

2.1 Inference and large sample

properties

Note that $\lambda_{1}=p\cdot\beta$ corresponds to the regression parameter vector when we regress $Y$ on $X$. Consequently by least squares regression we can obtain a closed form estimator $\hat{\lambda}_{1}$ of $\lambda_{1}$. Specifically with $X_{i}\in\mathbb{R}^{m}$ and $Y_{i}\in\mathbb{R}$ denoting the data from the $i$th individual and with $\tilde{X}=(\tilde{X}_{1},\ldots,\tilde{X}_{n})^{T}$, $\tilde{X}_{i}=(1,X_{i}^{T})^{T}$, and $Y=(Y_{1},\ldots,Y_{n})$ we have:

[TABLE]

where $A(m)$ is a $m\times(m+1)$ matrix projecting an $m+1$ vector to its last $m$ components.

The residual $Y_{i}-E(Y_{i}|\tilde{X}_{i})$ is not ensured to follow a Gaussian distribution — even if the original distributions of $Y_{1}$ and $Y_{2}$ were Gaussian — but we may still resort to a characterization of $\sqrt{n}(\hat{\lambda}_{1}-\lambda_{1})$ as a sum of i.i.d. zero mean terms ensuring consistency and asymptotic normality by means of the central limit theorem. Specifically the characterization is as follows:

[TABLE]

with

[TABLE]

Recall that equation (2) may be rewritten as

[TABLE]

with $\eta_{i}=\lambda_{1}^{T}X_{i}$ (the linear predictor), $\lambda_{2}=2\cdot\mu_{1}$, and $\lambda_{3}=p^{-1}$. From this observation and the fact that we have an estimator of $\lambda_{1}$ it seems only natural to estimate $\lambda_{2}$ and $\lambda_{3}$ by a multiple linear regression of $Y_{i}^{2}$ on $\hat{\eta}_{i}=\hat{\lambda}_{1}^{T}X_{i}$. However, since the variance of $Y_{i}^{2}$ given $X_{i}$ depends on $\eta_{i}$ and $\eta_{i}^{2}$ it is more appropriate to estimate $\lambda_{2}$ and $\lambda_{3}$ using a weighted regression, where the weight is a function of $\eta_{i}$.

Specifically let $w$ be a weight function, put $w_{i}=w(\eta_{i})$, $\hat{w}_{i}=w(\hat{\eta}_{i})$, and let $\mathbb{P}_{n}^{w}$ denote the weighted empirical mean, that is for $f_{i}=f(Y_{i},X_{i})$

[TABLE]

Also define

[TABLE]

Finally define

[TABLE]

Then our estimators are given as

[TABLE]

We now proceed to derive large sample properties of the estimators. In what follows we adopt the notation $\mathbf{P}f=Ef(X_{i},Y_{i})$ for any function of $f$ of our data samples $(X_{i},Y_{i})$

To derive large sample results we need to assume that the weight function $w$ is twice continuously differentiable and that $\lambda_{1}\rightarrow w(\lambda_{1}^{T}X)$ and the derivatives $\lambda_{1}\rightarrow D_{\lambda_{1}}w(\lambda_{1}^{T}X)$ and $\lambda_{1}\rightarrow D_{\lambda_{1}}^{2}w(\lambda_{1}^{T}X)$ are bounded by functions with finite means in some open neighbourhood of the true value of $\lambda_{1}$.

We also need the following lemmas

Lemma 2.1**.**

Let $f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ be a two times continuously differentiable function in the first coordinate. Assume that $\lambda_{1}\rightarrow f(\lambda_{1}^{T}X,Y)\cdot w(\lambda_{1}^{T}X)$ and the derivatives $\lambda_{1}\rightarrow D_{\lambda_{1}}\{f(\lambda_{1}^{T}X,Y)\cdot w(\lambda_{1}^{T}X)\}$ and $\lambda_{1}\rightarrow D^{2}_{\lambda_{1}}\{f(\lambda_{1}^{T}X,Y)\cdot w(\lambda_{1}^{T}X)\}$ are bounded by functions with finite mean in an open neighbourhood of the true value of $\lambda_{1}$. Also put $f_{i}=f(\eta_{i},Y_{i})$ and $\hat{f}_{i}=f(\hat{\eta}_{i},Y_{i})$. Then with

[TABLE]

it holds that $g$ is differentiable with derivative

[TABLE]

Moreover it holds that

[TABLE]

Lemma 2.2**.**

Let $f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ and put $f_{i}=f(\eta_{i},Y_{i})$. Assume that

[TABLE]

Then

[TABLE]

where $\mathbb{G}_{n}(f)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(f_{i}-\mathbf{P}f)$.

Theorem 2.1**.**

The estimators $\hat{\lambda}_{2}$ and $\hat{\lambda}_{3}$ are consistent, asymptotically normal, and have the following asymptotic characterizations:

[TABLE]

[TABLE]

where $a_{1}(\lambda_{1})$, $a_{2}(\lambda_{1})$, $a_{3}(\lambda_{1})$, $g_{3}(\lambda_{1})$, $\lambda_{2}(\lambda_{1})$, and $\lambda_{3}(\lambda_{1})$ are the limits in probability of $a^{n}_{1}(\lambda_{1})$, $a^{n}_{2}(\lambda_{1})$, $a^{n}_{3}(\lambda_{1})$, $g^{n}_{3}(\lambda_{1})$ $\lambda_{2}^{n}(\lambda_{1})$, and $\lambda_{3}^{n}(\lambda_{1})$, respectively, and

[TABLE]

As an immediate consequence of (4) and Theorem 1 we have the that the obvious estimator $\hat{\beta}=\hat{\lambda}_{3}\hat{\lambda}_{1}$ of $\beta$ is both consistent and asymptotically normal. In particular, we have the following iid decomposition of $\sqrt{n}(\hat{\beta}-\beta)$:

[TABLE]

where

[TABLE]

Notice that the iid decomposition above enables consistent estimation of standard errors based on the law of large numbers. Specifically

[TABLE]

and accordingly a consistent estimator of $\text{var}(\xi_{i})$ can be obtained by inserting the empirical counterpart of $\xi_{i}$ in the left hand side of (8).

2.2 Choice of weights

The optimal weights for weighted least squares regression are $w_{i}\propto\frac{1}{V(Y|X)}$. For the present we suggest using

[TABLE]

as weights for the first order moment regression, (1), and

[TABLE]

as weights for the second order moment regression (2). Both choices mimic the order of conditional variances as functions of $\eta_{i}$ and fulfills the requirements that the weight functions are twice continuously differentiable and locally bounded by functions with finite means. In the rest of the following we will use these weights.

3 Simulation study

We illustrate the empirical properties of the proposed estimation procedure through a series of simulations. Consider the following four scenarios:

Simple Gaussian mixture: $Y_{1}\sim N(1+X,1^{2}),Y_{2}\sim N(0,1^{2})$ and $X\sim N(0,1^{2})$. 2. 2.

Gaussian distribution with zero-inflation: $Y_{1}\sim N(1+X,1^{2}),Y_{2}=0$ and $X\sim N(0,1^{2})$. 3. 3.

Exponential-Gaussian mixture: $Y_{1}\sim\text{exp}(1)+X$ (shifted exponential distribution with rate $1$), $Y_{2}\sim N(0,0.5^{2})$ and $X\sim N(0,1^{2})$. 4. 4.

Exponential distribution with zero-inflation: $Y_{1}\sim\text{exp}(1)+X$, $Y_{2}=0$ and $X\sim N(0,1^{2})$.

For each scenario we simulate 100 datasets each with $N\in\{300,500,800,100,1500,2000,3000\}$ observations and with a constant mixture proportion of $p=0.7$ (i.e., 30% contamination). The estimates $\hat{\beta}$ and $\hat{p}$ are computed using both the moment mixture estimation procedure as well as a standard two-component Gaussian regression mixture model where $X$ can only influence one of the components.

Figure 1 shows the empirical estimates of $\hat{\beta}$ (and the corresponding 95% pointwise standard errors) for the four different scenarios while Figure 1 shows the corresponding plot for the mixing proportion $\hat{p}$. Not surprisingly, the Gaussian regression mixture model provides an unbiased estimate of both $\hat{\beta}$ and $\hat{p}$ with the smallest variance when the data follows a mixture of Gaussian distributions. However, it is also clear that the moment mixture model — which places fewer assumptions on the distributions involved — is also unbiased and has a variance that is only slightly larger than the Gaussian mixture model. Substantial differences are seen when the distribution of $Y_{1}$ is not Gaussian where the moment mixture model quickly converges to the true value while the Gaussian mixture model results in biased estimates of both $\beta$ and $p$. The estimates from the moment mixture model are unbiased even when the distributions are non-normal, but Figure 2 suggests that large samples may be needed in order to have sufficient power to use the moment mixture model to infer whether there are two or only one component.

The variance of the mixing proportion estimate, $\text{var}(\hat{p})$, is larger for the moment mixture model than for the Gaussian regression mixture model which essentially is the “price that is paid” by having fewer assumptions about the distribution and estimating the parameters from the moments of the distribution.

Table 1 shows results from the same four scenarios (although with mixing proportion $p=0.5$) to illustrate the precision of the estimates and to compare the estimated standard errors from (8) to the empirical standard errors. The estimates for both $\beta$ and $p$ are close to the true values even for smaller sample sizes and when the error distribution is non-normal. Table 1 also shows that the estimated standard errors are close to the empirical standard errors which suggests that the asymptotic estimate in (8) provides a useful measure of the standard errors of the regression and mixing proportion parameters even for smaller sample sizes. The 95% coverage probabilities for $\beta$ are close to the true value even for $N=300$ regardless of the underlying distributions while the coverage probabilities for $p$ are quite unstable for smaller sample sizes while they achieve right right level for $N=2000$. Note that the mixing proportion is 50% so there is a substantial amount of noise in the data.

To investigate the efficiency of the proposed estimates we consider the case where we know the individual allocations, $P_{i}$. If we make no assumptions on the distribution of $Y_{2}$ but assume that $Y_{1}$ follows a normal distribution with mean $\beta_{0}+\beta_{1}X$ and variance $\sigma_{1}^{2}$ conditional on $X\in\mathbb{R}$, then the most efficient way to estimate $\beta_{1}$ would be by simple linear regression based only on the pairs $(X_{i},Y_{i})$ where $P_{i}=1$. For this estimator the asymptotic variance is equal to

[TABLE]

Accordingly, the maximum likelihood estimator $\hat{\beta}_{1}$ of $\beta_{1}$ from any parametric mixture model of $Y$ with the above specification of $Y_{1}$ will have an asymptotic variance that is at least (9). We compare the estimated standard error from the moment mixture model to this lower bound to see the efficiency loss incurred by removing the assumptions about the parametric distributions. The result is seen in the right-most column of Table 1, and we see that the efficiency loss from using the moment mixture model is around a factor 4-5.

4 Application: The effect of pH on wine volatile

acidity

Volatile acidity (VA) refers to the process when lactic acid bacteria and acetic acids turns wine into vinegar, and the process takes place mainly due to growth of bacteria, the oxidation of ethanol, or the metabolism of acids/sugars. Wines with a high level of pH are supposedly more susceptible to oxidation and the antibacterial effects of sulfur dioxide and of fumaric acid are reduced rapidly as the pH level increases. Consequently, wines are thought to lose their quality as they become less acidic (increased pH) since the volatile acidity increases.

The paper by Cortez and others (2009) considers 11 physicochemical properties of a selected sample of Portuguese vinho verde wines. Samples from 1599 red wines and 4898 white wines are available and the relationship between volatile acidity and pH are shown separately for red and white wines in Figure 3. While it is apparent that the available red wines generally have slightly higher levels of pH, it also appears as if the impact of pH on volatile acidity is largest for the red wines in the sample: Individual regression lines for the two types of wine show an almost horizontal line for white wine while there is an effect of pH on VA for the red wines (slopes $-0.022$ and $0.272$, respectively, and the slopes are significantly different, $p<0.0001$).

If we did not know that the full 6497 samples were comprised of two different types of wine we might pursue regressing volatile acidity on pH for the full data. This gives the dashed regression shown in Figure 3 (slope $\hat{\beta}$ =0.27, 95% CI 0.24-0.29), which suggests an overall effect of pH on volatile acidity.

When we fit the moment-based mixture model then we see that the data are likely to consist of a mixture, $\hat{p}=0.28$ (95% CI: 0.24-0.33), which suggests that only a proportion of the wines are influenced by pH. Also, the effect driving the association between the volatile acidity and pH appears to be much stronger than what was observed from analyzing the full data with a simple linear regression model, $\hat{\beta}=$ 0.94 (95% CI 0.51-1.37). Consequently, the moment mixture model is able to identify that the wine data is likely to consist of two types of wines that respond differently to changes in levels of pH with only placing very few restrictions on the underlying distributions.

We estimate the mixing proportion to be 28% whereas the dataset contains 24.6% red wines. While the wine colour classification need not correspond to the separation we estimate from the moment mixture model it may indeed be the case that the different wine types contain a set of features that influence the impact of pH which is what we observe for these data.

5 Discussion

In this article we address the problem of estimating regression parameters for a two-component mixture regression model where one component is influenced by a set of predictors. Our proposed method imposes no assumptions on the distributions of the components and only a minimum of restrictions on the regression effects. The moment-based mixture of regression models estimators can be used to detect or account for unspecified mixtures in regression problems, and since estimation is fast it could be used for large-scale studies such as in genome-wise association studies.

As shown by the simulation study, the price for the flexibility and lack of assumptions comes with a larger variance of the estimators but the variance is not prohibitively larger and the proposed methods provides consistent estimates even when the mixture distributions are highly skewed and other models such as the Gaussian mixture model fails.

In conclusion, we have introduced an estimation technique for mixtures of regression models that can be applied to a large number of situations. The moment-based estimator is very versatile: it can be used not only to estimate the regression parameters, but for larger datasets it provides a foundation for detecting the number of mixture components. If $\hat{p}$ is different from both 0 and 1 then this suggests that there indeed are two components where only one of them is influenced by the predictors of interest.

6 Software

Software in the form of the mommix R package and code, together with a sample input data set and complete documentation is available on GitHub at www.github.com/ekstroem/mommix.

Conflict of Interest: None declared.

7 Appendix

7.1 Proof of Theorem 2.1

Denote

[TABLE]

with corresponding limits in probability

[TABLE]

where

[TABLE]

By repeatedly using Lemma 2.1, the consistency of $\hat{\lambda}_{1}$, and standard arguments we have

[TABLE]

Then$\sqrt{n}(\hat{\lambda}_{2}-\lambda_{2})$ may be rewritten as:

[TABLE]

For the first term of the right hand side of equation (10)

[TABLE]

by the consistency of $\hat{\lambda}_{1}$ and the law of large numbers. Combining the above results we obtain

[TABLE]

For the second term on the right hand side of equation (10) first note that by repeatedly using Lemma 2.2 and the law of large numbers

[TABLE]

Next, by the law of large numbers

[TABLE]

Now note that in the true value of $\lambda$ we have

[TABLE]

and thus

[TABLE]

Accordingly

[TABLE]

Finally, combining all of the above and using the law of large numbers:

[TABLE]

where

[TABLE]

Adding (12) and (14) we obtain

[TABLE]

Similarly one may show

[TABLE]