Semiparametric correction for endogenous truncation bias with Vox Populi based participation decision

SEMIPARAMETRIC CORRECTION FOR ENDOGENOUS TRUNCATION BIAS WITH VOX POPULI BASED PARTICIPATION DECISION

\history{This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2018.2888575, IEEE Access} \doi{10.1109/ACCESS.2018.2888575}

\address[1]{University of Haifa, Haifa, Israel (e-mail: [email removed])} \address[2]{University of Haifa, Haifa, Israel (e-mail: [email removed])}

\markboth {This article has been accepted for publication in a future issue of this journal, but has not been fully edited.} {This article has been accepted for publication in a future issue of this journal, but has not been fully edited.}

\corresp{Corresponding author: Moshe kim (e-mail: [email removed]).}

\titlepgskip=-15pt

\IEEEpeerreviewmaketitle

Introduction

An important fact, but one that is largely overlooked or taken-for-granted, is that researchers hardly ever have access to the entire data distribution pertaining to their specific research and rely, instead, on a truncated form of such data. The truncated data employed probably has different characteristics than the latent non-truncated full distribution and may result in biased parameter estimates generated by the specific investigated models. The problem is further aggravated when truncation is endogenously propagated by various decision units, or observations. Examples of a straight-forward endogenous truncation emerge from aspects of some type of discouragement. For instance, in labor markets, long-term unemployed persons are often discouraged workers who are afraid that they will not find employment and therefore do not seek employment; hence, they will be absent from reported unemployment rates. Another, gender-related, labor market example is discouraged women: it has been shown that most women do not apply for jobs that require a high degree of aggressiveness; hence we have discouraged women who fail to participate in specific sectors of the economy, which affects female labor supply. In financial markets, discouraged borrowers, such as some small and medium size enterprises, do not apply for loans and result in biased modeling of default probabilities, which hampers optimal credit allocation. Endogenous truncation also is involved in the measurement of social problems, such as the crime or divorce rates, because measurement may represent a latent rate of reporting, rather than the variables of interest. Similarly, endogenous truncation severely impacts the measurement of important economic indices, such as economic growth, productivity, income distribution, and welfare.

The notion of truncation is different than the concept of the known selection bias. In the known selection bias, information on data (observations) has been censored but still observable or, alternatively, information regarding the counterfactual (e.g., the rejected rather than the discouraged borrower) has been censored but still observable. Selection bias under censoring has already been remedied by Heckman's seminal contribution heckman1979sample. Under Heckman's model, the selection process is entirely observed and selectivity bias can be alleviated. Under endogenous truncation, however, the selection rule is completely unobserved and no information is available concerning the truncated observations. Thus, statistical biases are myriad and interwoven to the extent that researchers may not even be able to assess their magnitude and direction, and the problem becomes extremely challenging. Given the potential severity of the aforementioned problem, it is surprising that the endogenous truncation problem has attracted hardly any scientific investigation, assessment, or suggestions for proper remedies.

In the few existing interests in the literature, the identification of the semiparametric truncated sample selection model is achieved by observing the selection variable (which is modeled as continuous), while imposing different restrictions on the disturbances powell1994estimation,honore1997estimation, or by utilizing information regarding some of the non-participants' characteristics khan2007weighted. Both of these studies rely on available data regarding either the covariates' joint distribution function or the selection variable, implying that the variable is not treated as a latent binary response variable (unlike the approach taken in the present paper). Further, these studies model the selection rule of each datum as a function of its observed characteristics. Yet, the selection rule might be affected by unobserved (truncated) characteristics, as well. Ignoring these characteristics may lead to misspecification of the selection equation, potentially biasing the estimates. Additionally, the estimation and identification of semiparametric truncated sample selection models with a latent binary selection variable are known to be difficult, due to the absence of observed variation in exactly this selection variable. The various estimation procedures that utilize a continuous selection variable to alleviate this difficulty use different kernel estimators lee1993quadratic. The closest approach to the proposed methodology, dealing with a latent binary selection variable, is ichimura1993semiparametric, which also employs a kernel to estimate the bias term in the substantive equation ichimura1991semiparametric,ichimura1993semiparametric. However, the resulting estimates can still be biased, as the kernel estimator's accuracy depends on selecting the optimal bandwidth, which is hard to find in the semiparametric context lewbel2007simple.

An additional, important weakness of the existing literature dealing with endogenous truncation problems is the assumption of similar behavior on the part of the truncated and non-truncated distributions ichimura1993semiparametric, an assumption which is referred to as a population regression, in the econometric literature heckman1979sample, and a covariate shift, in the computer science literature gretton2009covariate. The various truncated sample selection models treat the data as if they all consist of a single, homogeneous, monolithic cohort sharing identical actions, such that the selection rule of each datum is not affected by the participation decisions of other members. This restrictive assumption, however, can introduce selection bias by itself. In fact, as manski2010consensus describes it: “If agents knew the state of nature, they would make the same decision. However, they may have different beliefs or may use different decision criteria to cope with their incomplete knowledge. Hence, they may use different actions even though they share the same objective” (p.187).

Taking into consideration that we are unable to observe the selection variable, we propose an estimation procedure. In order to rectify the aforementioned potential bias and to improve upon the covariate shift assumption, which is frequently used in machine learning, the data in our model are treated as a mixture of sub-populations, each characterized by its own action regarding the participation decision. Thus, we build on the vox populi concept galton1907vox or, in its modern term, “The Wisdom of Crowds” surowiecki2005wisdom, as the basis by which data points “sort themselves” in the truncation process. As such, each data point's “decision” to allow itself to be truncated from the original distribution is an important building block that generates our offered algorithm.

The vox populi concept relies on the idea that aggregates of opinions measuring the central tendency will be more accurate than individual opinions davis2014crowd. budescu2005confidence suggests that an aggregate of multiple sources maximizes the amount of information available and reduces the potential impact of unreliable information sources. The implication is that the combination of the various sources leads to error cancellation.

Further, we refine the concept of “the wisdom of crowds” to be a non monolithic concept and apply it to truncation. Each observation “decides” whether to allow itself to be truncated depending on its reference group's (rather than on the entire crowd's) decision opinion space average forecast. This, in turn, is inspired by the similarity-based classification in machine learning and Cybernetics chen2009similarity,hummel1996statistical and management science models of decision making budescu2014identifying. This enables the various opinion spaces, generated by the various reference groups, to provide expert opinion with rather superior average forecast, by eliminating poorly-performing individuals from the crowd prelec2017solution,budescu2014identifying. Such treatment is also inspired by economic theories of ethnic capital Borjas1992 and informational cascades bikhchandani1992theory, highlighting the fact that individual characteristics depend on the average characteristics of the group to which they belong. Recently, we have witnessed an upsurge of interest in the relationship between culture and genetic diversity desmet2017culture, through the process of endogenous group selection ashraf2013Genetic.

Building upon this insight, we model the number and type of reference groups to be endogenously determined, rather than arbitrarily imposed. A Latent Class Analysis (LCA) clogg1984latent is used to estimate the latent characteristics (type) of the various reference groups. This is implemented by integrating Machine Learning concepts and providing a Fourier-based Sieve semiparametric estimator, which is distribution-free. Our estimator uses a penalized non-linear regression schuurmans2002metric, an important characteristic emphasizing the generality and applicability of the offered methodology. The Fourier series is a functional of the Orthonormal polynomials sequence family, which allows for efficient estimation of functions with non-smoothness, discontinuities in derivatives, sharp spikes and discontinuities in the function itself. Thus, it is useful in nonparametric regression for approximating a much broader class of functions ogden2012essential than the kernel approach.

The most attractive feature of our proposed estimator is that it intrinsically prevents potential multicollinearity problems. Even though the multicollinearity might arise in certain circumstances, we can prevent it. For example, multicollinearity might arise if we extend the model by incorporating an endogenous covariate in the substantive equation and estimate sequentially a system of partially linear equations. The first equation is the endogenous covariate regression, which linearly depends on the selection bias term, while the second equation is the substantive equation (of interest), which linearly depends both on the endogenous covariate, as well as on a similar selection bias term. These two selection bias terms depend on the same covariate vector and thus they might be correlated. However, this problem is alleviated by the fact that each selection bias term is approximated by a different orthonormal polynomial sequence (a different number of mutually orthogonal basis functions), which implies, by definition of orthonormality, that these two approximated functions cannot be perfectly multicollinear. This result is required for identification.\footnote{The identification can be achieved due to the fact that some of the mutually orthogonal basis functions (covariates) are not common to both series expansions. These non-common covariates play the role of an exclusion restriction which is commonly used to assure identification.} We note that the classical kernel estimator does not possess this advantageous orthonormality feature and consequently may produce biased estimates due to cross- equation correlation.

Another aspect that our estimator must consider is the optimal number of groups. In order to find the optimal number of groups that best fits the data generation process, we perform variable selection (also referred to as "sparse regression" friedman2012fast,yang2016sparse) by employing the smoothly clipped absolute deviation (SCAD) penalty function.\footnote{The SCAD penalty function is superior to the often employed least absolute shrinkage and selection operator (LASSO), because it is general and nests the LASSO as a special case.} We develop a generic non-linear penalized regression estimation method, in the sense that it can easily be extended to enable a wide collection of penalty functions to be estimated. The novelty of our modeling lies in the integration and synthesis of knowledge present in various scientific disciplines, such as: (i) computer science (pattern recognition, unsupervised machine learning,\footnote{For a constructive overview of the field of unsupervised learning, see ghahramani2004unsupervised.} artificial intelligence and self-organizing maps in neural networks); (ii) electrical engineering (signal extraction); (iii) economics and; (iv) management for the creation of new algorithms correcting for truncation bias, due to the endogenous self-selection of observations into a sample. This integration enriches the algorithms' accuracy, efficiency and applicability and hopefully can be of use in economics and many other disciplines.

We offer a three-stage procedure to correct for the endogenous truncation bias: in the first stage, latent classes analysis is employed based on results from an auxiliary survey data, consisting of experts' (binary) opinions, as well as of their observed group characteristics, to recover the unobserved latent reference groups.\footnote{For example, the Small Business Credit Survey administered by FederalNY is an annual survey of firms with fewer than 500 employees reporting on financing needs and choices and borrowing experiences. Based on the small business credit survey 2016, a total of $17\%$ of the non-applicants are discouraged borrowers.} A given expert's opinion captures his belief regarding the expected participation decision in his reference group.\footnote{Since the opinion is binary, each expert is asked what is the most likely decision for a random member belonging to her reference group being a participant or a non-participant.} In the second stage, each participant share is obtained by averaging the members' opinions belonging to the specific reference group.\footnote{The average of opinions belonging to a particular reference group reflects a refined version of the (monolithic) wisdom of crowds.} In the third stage, a semiparametric truncated sample selection model is estimated, consisting of a selection equation and a substantive equation. The estimated participants' share in the group, conditional on the reference group's observed and unobserved characteristics, is included in the selection equation as an additional covariate. We run Monte Carlo simulations in order to examine our estimator's performance in the presence of a truncated sample selection model. Further, for sake of generality of the offered estimator, we subject it to various distributions in which the disturbances are neither jointly nor marginally normally distributed. These disturbances are constructed as realizations of non-symmetric and non-unimodal distribution functions.\footnote{Unlike the practice in some other studies applying only normally distributed disturbances.}

The rest of the paper is organized as follows: Section (ref) introduces the model consisting of substantive and selection equations; Section (ref) deals with model estimation; Section (ref) recovers the number of reference groups; Section (ref) examines our truncated selection model's performance, employing Monte Carlo simulations; and Section (ref) concludes by summarizing the main findings, as well as our estimator's performance.

Next we present our suggested methodology for estimation of a truncated endogenous sample selection model in the presence of a reflection problem, when the entire data consist of participants only.

The model

Binary choice selection equation

We describe the participation choice of each individual observation $i\in\left \{{1,...,N} \right\}$, as a function of its reference group's participation decision which is captured by the participants' share in the particular reference group. A model in which an individual's decision is affected by the average decision made by all its group members is referred to as the “reflection problem” manski1993dynamic,\footnote{The reflection problem arises manski1993dynamic “when a researcher observing the distribution of behavior in a population tries to infer whether the average behavior in some group influences the behavior of the individuals that comprise the group. The term reflection is appropriate, because the problem is similar to that of interpreting the almost simultaneous movements of a person and his reflection in a mirror. Does the mirror image cause the person's movements or reflect them?” (p. 532)} or Manski's notion of role models/emulation manski1993identification. These concepts may touch on an earlier idea of “ethnic capital” Borjas1992, showing individual characteristics to be dependent on the average characteristics of the group they belong to, and a tendency to follow the decision of others bikhchandani1992theory.

Let the number of reference groups (unknown to the researcher) be denoted by $G$. The individual choices given a membership in reference group $g\in\left \{{1,...,G} \right\}$ are coded by $\omega_{i,g}\in\left \{{0,1} \right\}$ and are defined as:

These individual choices are determined by two sets of factors. The first set consists of the observed group-level characteristics $\boldsymbol{{x}}_{g}\in\mathbb{R}^{L_{\boldsymbol{{x}}}}$\footnote{The notation $\mathbb{R}^{L_{\boldsymbol{{x}}}}$ stands for a vector of size $L_{\boldsymbol{{x}}}\times 1$.} and the unobserved group-level characteristics, captured by a latent categorical variable $\psi_{i,g}$ of $\tilde{G}$ different outcomes, where $\tilde{G}$ is not arbitrarily imposed (as will be depicted in section (ref) to follow).\footnote{We allow for (but do not require) a dependence between the observed and unobserved group's characteristics, determined by some unknown joint distribution function (as depicted in section (ref) to follow).} The second set consists of the observed individual-level characteristics $\boldsymbol{{z_i}}\in\mathbb{R}^{L_{\boldsymbol{{z}}}}$ and an individual random disturbance $\xi_{2i}$.\footnote{Each reference group $g$ is a unique combination of observed and unobserved characteristics ($\boldsymbol{{x}}_g$ and $\psi_{i,g}$, respectively) which are common to all of the $g$'th reference group's members. However, the presence of unobserved characteristics $\psi_{i,g}$, implies that in order to assign observations into reference groups, $\psi_{i,g}$ is required to be estimated (as will be discussed in section (ref) to follow).}

These factors are assumed to produce payoffs for the possible participation choices, $u_{i,g}(1)$ and $u_{i,g}(0)$, the utility of participation and non-participation, respectively. The difference between these payoffs is additive in the various factors. A participation choice is made when the following difference is positive brock2007identification:\footnote{$T$ is defined everywhere in the manuscript as the transpose operator.}

where $\boldsymbol{{x}}_{g}^c$ is a subset of $\boldsymbol{{x}}_{g}$ consisting of contextual factors,\footnote{ This decomposition is intended to satisfy the exclusion restriction in (ref) for the sake of identification of the $\beta$ and $\delta$ parameters.}$^,$\footnote{A contextual effect exists whenever the propensity of a person to behave in some way varies with the characteristics of the reference group members.} $m^{\mathlcal{e}}(\boldsymbol{{x}}_g,\psi_{i,g})$ is the expectation (forecast) of individual $i$ with reference group's characteristics $\boldsymbol{{x}}_{g}$ and $\psi_{i,g}$ regarding the participants' share in his group, and the super-script $\mathlcal{e}$ represents expectation (forecast).\footnote{The subjective belief (forecast) is a mapping from group's (observed and unobserved) characteristics to a scalar representing a participation probability (participants' share).}

It is worth noting that the difference $u_i(1)-u_i(0)$ in (ref) is positive iff the following inequality holds:\footnote{Instead of employing merely the average participation decision $m^{\mathlcal{e}}(\boldsymbol{{x}}_g,\psi_{i,g})$, an interesting extension of this model would be to allow for each datum to be affected by a vector of moments (various dispersion measures) obtained from the survey.}

which implies that the conditional participation probability given the reference group and individual level characteristics $\left \{{\boldsymbol{{z_i}},\boldsymbol{{x}}_g,\psi_{i,g}} \right\}$ with $(\boldsymbol{{x}}_{g}^c\subset \boldsymbol{{x}}_g)$ is:

where $F_{\mathrm{\xi_2}}$ stands for the distribution function of the random disturbance $\xi_{2i}$, which is unknown to the researcher, $\mathlarger{\mathlarger{\omega_{i,g}}}$ is a random variable that is conditionally Bernoulli-distributed, given the individual-level and group-level covariates, while $\omega_{i,g}$ depicted in (ref) stands for its realization.

Each individual is small relative to the population blume2010identification. Using (ref), the following condition is obtained:

where $m(\boldsymbol{{x}}_g,\psi_{i,g})$ is the actual participants' share, given a membership in a reference group characterized by observed and unobserved characteristics $\boldsymbol{{x}}_g$ and $\psi_{i,g}$, respectively; $F_{\boldsymbol{{\mathrm{z}}}|\boldsymbol{{x_{g}}}}$ is the conditional distribution function of $\boldsymbol{{\mathrm{z}}}$ (given $\boldsymbol{{x}}_{g}$) which is unknown to the researcher.

We next present the theoretical model equations.

The sample selection model

The underlying model consists of two equations in which the latent (population) dependent variables $y_{1i,g}^*$ and $y_{i2,g}^*$ are defined as follows:

$\hspace{3em}$ and

where $\boldsymbol{{w_i}}\in\mathbb{R}^{L_w}$ and $\boldsymbol{{\theta}}$ denote the substantive equation's covariate vector and a $L_w\times 1$ parameter vector, respectively. The substantive equation's random disturbance is $\xi_{1i}$, and the selection equation's disturbance satisfies $\tilde{\xi}_{2i}\equiv -\xi_{2i}$.\footnote{Using the definition in (ref), $y_{2i,g}^*$ is the difference between the participation and non-participation utilities, which includes a random disturbance $\xi_{2i}$ followed by a minus sign.} $\text{ }$The random disturbances $\xi_{1i}$ and $\xi_{2i}$, with their respective marginal distribution functions $F_{\xi_1}$ and $F_{\xi_2}$, are jointly distributed. Their joint distribution function is $F_{\xi_1, \xi_2}$. The model is semiparametric as neither the marginals nor the joint distribution function are required to be specified by the researcher. $y_{ji,g}^{*}$ denotes a realization of the latent random variable $\mathrm{y}_j^*$ for $j=1,2$.\footnote{Asterisk implies a latent (population) variable.} The group-level characteristics $\boldsymbol{{x}}_g$ and $\psi_{i,g}\in\left \{{1,...,\tilde{G}} \right\}$ constitute the $i$'th observation's specific reference group; $m^{\mathlcal{e}}(\boldsymbol{{x}}_g,\psi_{i,g})$ is the latent participants' share given a membership in latent reference group $(\boldsymbol{{x}}_g,\psi_{i,g})$. $\beta$ captures the endogenous effect\footnote{The presence of an endogenous effect implies that the propensity of a person to behave in some way varies with the behavior of the reference group manski2000economic. } and $\delta$ the contextual effect.

In the truncated sample the $i$'th observation in group $g$ is denoted by the sequence $\left \{{y_{1i,g},\boldsymbol{{x}}_{g}^T,\boldsymbol{{w}}_i^T,\boldsymbol{{z}}_i^T} \right\}$, where $y_{1i,g}$ is defined as:

A binary random variable indicating participation is denoted by $\mathrm{S}^*$ defined as:

However, $\mathrm{S}^*$ in (ref) is unobserved, and only $\mathrm{S}$ is observed:

Let $n<N$ denote the number of observations in the truncated data set. The participants' share $m^{\mathlcal{e}}(\boldsymbol{{x}}_g,\psi_{i,g})$ is a forecast of the actual participants' share $m(\boldsymbol{{x}}_g,\psi_{i,g})$, and they are interrelated through (ref).

Next we discuss the model estimation.

Model Estimation

In this section, we propose an estimation procedure for a truncated selection model, consisting of a substantive equation and a selection equation.

The estimation is a three steps sequential procedure: (i) A Latent Classes Analysis to estimate the reference groups' unobserved characteristics, as will be discussed in section (ref); (ii) Evaluation of the participation probability in each reference group, controlling for its unobserved characteristics, by utilizing experts' opinions; and (iii) Estimating a partially linear index model using Sieve (series) estimator for the non linear component, which is referred to as the “bias term”.

Next we discuss the main idea behind the assignment of each observation into latent classes, utilizing a survey data set consisting of experts' opinions and group-level covariates (a combination of continuous and categorical variables).

Identification of latent reference groups

Latent Class Analysis (LCA) is a statistical method for matching a set of manifest (observed) variables to a set of latent variables referred to as classes goodman1974exploratory,clogg1984latent,greene2003latent. A specific realization of the manifest variables is referred to as a “response pattern”. Let $\mathfrak{Y}$ denote the set of response patterns consisting of all possible realizations of a $J\times 1$ categorical variable vector, defined as:

where the number of outcomes in the $j$'th categorical variable is $K_j$.

The role of the realizations of manifest variables in (ref) is for identification purposes, by means of classifying observations into their most likely latent class utilizing recruitment probabilities. A recruitment probability is the probability that a specific response pattern $\scaleobj{0.8}{\boldsymbol{{\mathlcal{Y}}}}\in\mathfrak{Y}$ will be observed for a randomly selected member of a given latent class.\footnote{The response pattern of the $i$'th observation is its set of responses to all the manifest variables. These responses are conditionally independent of each other in a given class.} The a posteriori probability of being a member in a given class is obtained by using Bayes' theorem as a function of the estimated recruitment probabilities and the estimated prevalence of each latent class (the class membership prior probability). Each observation is assigned to the latent class that has the highest a posteriori probability.

The number of latent classes (labels), $\tilde{G}$, is recovered by the model rather than arbitrarily imposed. The latent classes analysis is employed repeatedly for a given specific number of latent classes $\mathlcal{k}\in\left \{{2,...,\tilde{G}_{\max}} \right\}$. $\tilde{G}_{\max}$ is the largest possible number of classes and is in the spirit of the Bayesian Information Criterion (BIC) schwarz1978estimating, reported to perform well by finding the correct number of components in the mixture roeder1997practical. Other authors suggest using Bayesian-based graphical techniques to aid in deciding on the number of classes garrett2000latent. We depart from the aforementioned literature in that we apply the BIC criterion directly to the substantive equation, in order to find the best model specification under truncation by using SCAD (section (ref)). To achieve this goal, we employ a penalized non-linear regression model, using the SCAD penalty function, to select the best solution obtained from the latent classes analysis.

We distinguish between two cases: (i) the class membership prior probabilities varies among observations, as these probabilities are determined by a covariates set; (ii) the class membership prior probabilities are constant across observations and there is no dependence on covariates. In the former, a parametric multinomial response model, such as the multinomial logistic regression, is employed to estimate the prior class membership given the covariates.\footnote{The justification for a parametric model is to reduce the complexity of calculations.} For the latter, one only needs to estimate $\mathlcal{k}-1$ class membership proportions (given $\mathlcal{k}$ classes) that characterize the entire data.\footnote{Without loss of generality, these unknown proportions can be estimated, nonparametrically, by a logistic multinomial response model characterized by a unique intercept per class. This is a nonparametric estimation procedure, due to the absence of covariates.} As we focus on the endogenous determination of class membership, covariates are involved in the estimation of the prior class membership probabilities. The model parameters that are required for the estimation of the labels' a posteriori distribution (in section (ref) to follow) are: (i) the parameters which affect the class membership prior probability (in section (ref) to follow); and (ii) the parameters which affect the response pattern given the class membership (the conditional response probabilities in section (ref) to follow).

The latent classes estimation

In this section, we introduce an estimation procedure to recover the outcomes of the sequence $\left \{{\psi_{i,g}} \right\}_{g=1}^G$, which are the reference groups' latent characteristics. Each outcome has its own label, and the labels are estimated by employing latent classes analysis clogg1984latent,goodman1974exploratory, a procedure to estimate their labels' a posteriori probability density function. Once this posterior function is estimated, the sequence of fitted labels are the arguments that maximize ($\arg \max$) the estimated posterior distribution function.

Suppose that the population consists of $\mathlcal{k}$ latent classes, such that the class membership of each observation $i=1,...N$ is denoted by an unobserved categorical variable $\psi_{i,g}$ with $\mathlcal{k}$ possible outcomes. We treat each observation as a random realization of the conditional labels' distribution function, given its groups' observed covariates $\boldsymbol{{x}}_{g}$. This methodology is based on the non-random assignment into classes (heterogeneous class membership prior probabilities) introduced by clogg1984latent.

Next we present the prior class membership probabilities under non-random assignment.

The prior distribution function of the labels

Let $\Psi_{\mathlcal{k}}$ be a categorical random variable of $\mathlcal{k}$ potential outcomes. We denote the prior probability of belonging to label $t$, given the group's observed characteristics $\boldsymbol{{x}}_g$ by $\lambda_{t|\boldsymbol{{x}}_g}$, satisfying $\sum_{t=1}^{\mathlcal{k}}\lambda_{t|\boldsymbol{{x}}_g}=1$ $\forall g\in\left \{{1,...,G} \right\}$ and defined as:

where $\boldsymbol{{x}}\in\mathbb{R}^{L_{\boldsymbol{{\mathrm{x}}}}}$ is a group-level covariates vector, $\mathlarger{\mathlarger{\varsigma}}_{_{0,1}},...,\mathlarger{\mathlarger{\varsigma}}_{_{0,\mathlcal{k}}}$ are intercepts and $\boldsymbol{{\mu}}_t\in\mathbb{R}^{L_{\boldsymbol{{\mathrm{x}}}}}$ for $t=1,...,\mathlcal{k}-1$ are parameter vectors.\footnote{Although this function can be formulated nonparametrically, we have opted for the present multinomial logistic formulation for computational simplification. Latent classes analysis involves an iterative estimation procedure, and thus each iteration requires a different optimal bandwidth. Since we estimate 10,000 different data sets, the number of bandwidths to be computed would requires 10,000 times the number of iterations. Computationally, this is extremely cumbersome.}

Suppose, also, that conditional on $\boldsymbol{{x}}_g$, $\Psi_{\mathlcal{k}}$ is jointly distributed with a vector of $J$ categorical variables (the group's manifest variables) $\boldsymbol{{\mathcal{Y}_i}}=\left[ {\mathcal{Y}_{1i},...,\mathcal{Y}_{Ji}} \right]^T$, which is referred to as the vector of responses and its realization is denoted by $\scaleobj{0.8}{\boldsymbol{{\mathlcal{Y}}}}_i=\left[ {\mathlcal{y}_{1i},...,\mathlcal{y}_{Ji}} \right]^T\in\mathfrak{Y}$. The $j$'th observed categorical variable (for each observation) $\mathcal{Y}_{ji}$ contains $K_j$ possible outcomes.\footnote{These categorical variables may have different numbers of outcomes, hence the indexing by $j$.} \color{black}

Let $\mathrm{D_{ijk}}$ be an indicator variable equal to unity, if respondent $i$ gives the $k$'th response to the $j'$th variable, and equals zero otherwise:

Next, we construct the recruitment (response) probabilities; each denotes the probability of observing a specific response pattern, given the class membership.

The recruitment (response) probabilities

Let $\pi_{jtk}$ be the probability that an observation in class $t$ produces the $k$'th outcome on the $j$'th variable. The recruitment probabilities are class-dependent, but are assumed to be homogeneous within classes, which implies that the following must hold:

Under conditional independence, which is a necessary condition for the class membership identification, the manifest variables are independent of each other, given the class membership and the group's observables characteristics. The probability that observation $i$ in class $t$ produces a particular set of $J$ outcomes on the observed categorical variables is the product:

For any given class $t$ and observed categorical variable $j$, the following requirement must be satisfied $\sum_{k=1}^{K_j}\pi_{jtk}=1$.

The probability density function across all classes is the total probability over the conditional probability in (ref):

where the parameters to be estimated by the latent class model are $\lambda_{t|\boldsymbol{{x}}_g}$ and $\pi_{jtk}$.

A posteriori distribution function of The latent labels

Given the prior and the recruitment probabilities' estimates for $\widehat{\lambda}_t$ and $\widehat{\pi}_{jtk}$, respectively, the posterior probability that a given individual belongs to a given class, conditional on the observed response pattern $[\mathlcal{y}_{1i},...,\mathlcal{y}_{Ji}]$ is:

where $t\in\left \{{1,...,\mathlcal{k}} \right\}$.

The log-likelihood function to be maximized, with respect to the parameters values in the prior and the recruitment probabilities $\lambda_{t|\boldsymbol{{x}}_g}$ and $\pi_{jtk}$ using (ref) is:

where the estimation procedure is expectation-maximization (EM) algorithm dempster1977maximum.\footnote{The EM algorithm enables us to maximize the log-likelihood function in (ref), iteratively, to simplify the estimation process. Moreover, in the absence of slope covariates in both the class-prior and class-conditional probability functions, these probabilities are estimated nonparametrically. However, in the present case, the class-prior probability functions are estimated parametrically, due to the non-random assignment embedded in the presence of covariates. This is important for satisfying the non-covariate shift notion, as has been discussed earlier. In an important paper by greene2003latent, a similar likelihood function is maximized, using a parametric technique.}

This log-likelihood function is identical in form to the standard finite mixture model log-likelihood. As with any finite mixture model, the EM algorithm is applicable, because each individual's class membership is unknown and may be treated as missing data mclachlan2000mixtures,mclachlan2007algorithm.

The EM algorithm is an iterative procedure involving two sequential steps: an expectation and maximization. First, initial parameter values $\widehat{\lambda}_t^{\mathrm{old}}$ are arbitrarily chosen and $\widehat{\pi}_{jtk}^{\mathrm{old}}$ for each $t\in\left \{{1,..,\mathlcal{k}} \right\}$ and $j\in\left \{{1,..,J} \right\}$ for all $k\in\left \{{1,..,K_j} \right\}$. In the expectation step, calculate the "missing" class membership probabilities using (ref):

In the maximization step, we update the parameter estimates by maximizing the log-likelihood function in (ref), given the estimated posterior in (ref). The new-prior probabilities are:

and the new class conditional probabilities are:

We replace the old estimates $\widehat{\lambda}_t^{\mathrm{old}}$ and $\widehat{\pi}_{jtk}^{\mathrm{old}}$ with the new estimates $\widehat{\lambda}_t^{\mathrm{new}}$ and $\widehat{\pi}_{jtk}^{\mathrm{new}}$, respectively, and repeat the expectation and maximization steps in (ref)-(ref), until a convergence criterion is satisfied for these new parameter values.

Using the estimated posterior function in (ref), the sequence of fitted labels $\widehat{\psi}_{i,g}$ are the arguments that maximize ($\arg \max$) the estimated posterior distribution function. Thus, given a response pattern $\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}}_i=[\mathlcal{y}_{1i},...,\mathlcal{y}_{Ji}]$ and group's observed characteristics $\boldsymbol{{x}}_g$, the fitted label for the latent $i$'th datum is:

Next, we utilize the experts' opinions in each reference group to evaluate the participants' share. The reference groups are identified by using both the group's observed characteristics and the fitted labels in (ref), capturing its unobserved characteristics.

The opinion space

We introduce an opinion space composed of a set of experts defined as:

in which an expert $\varphi\in\Lambda$, a set of observed characteristics $(\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E,\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}}^E)$ and unobserved characteristics ($\psi^E$), has a discretized opinion $\chi_{\varphi}\in\left \{{0,1} \right\}$ regarding the expected participation decision of a member belonging to own reference group.\footnote{An expert opinion reflects the decision that a member of his group is more likely to make. That is, being a participant or a non-participant.}

Let $\varphi_e=(\varphi_e^*,\psi_e^E)\in\Lambda$, where $\varphi_e^*=(\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_e,\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}}^E_e)$. Denote a random sample $\Lambda^S=\left \{{\varphi_e^*,\chi_{\varphi_e}} \right\}_{e=1}^{N^E}$ consisting of $N^E$ experts $\left \{{\varphi_e^*} \right\}_{e=1}^{N^E}$ and their respective opinions $\left \{{\chi_{\varphi_e}} \right\}_{e=1}^{N^E}$, where $\chi_{\varphi_e}\in\left \{{0,1} \right\}$. Each of the opinions in $\left \{{\chi_{\varphi_e}} \right\}_{e=1}^{N^E}$ is an independent realization of a Bernoulli random variable $\mathlarger{\mathlarger{\mathlarger{\omega}}}$, with probability of success defined by the function $m^{\mathlcal{e}}(\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_e, \psi_e^E)$.\footnote{Not to be confused with the actual participants' share $m(\boldsymbol{{x}}_g^E,\psi_{i,g}^E)$ depicted in (ref).}

It follows that $\Lambda^S$ consists entirely of the experts' observed characteristics and their opinions. The unobserved characteristics are essential for being able to assign the experts into their respective reference groups. However, the unobserved and observed characteristics are interrelated, through the a posteriori probability density function depicted in (ref). The former are substituted with their fitted values, which are the arguments maximizing the posterior probability density function, given the observed characteristics. Using the aforementioned interrelationship and given the sample $\Lambda^S$, the set of expert opinions that are assigned to latent class $t$ is denoted by:

The entire experts' opinions data set is denoted by the sequence $\left \{{\mathcal{O}_{p_{(t)}}} \right\}_{t=1}^{\mathlcal{k}}$.

The participants' shares given $\mathlcal{k}$ latent classes are obtained by Bayes' rule:\footnote{The expression $\mathrm{Pr}\left( {\Psi_{\mathlcal{k}}=t} \right)$ is canceled out and thus, is not presented in either the numerator or the denominator in (ref). }

where $m_{\mathlcal{k}}^e(\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E,t)\equiv \mathrm{Pr}\left( {\mathlarger{\mathlarger{\omega}}=1|\scaleobj{1.2}{\boldsymbol{{\mathrm{x}}}}=\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E,\Psi_{\mathlcal{k}}=t} \right)$.

Neither of the density functions $f_{\scaleobj{1.2}{\boldsymbol{{\mathrm{x}}}}|\mathlarger{\mathlarger{\omega}}=0,\Psi_{\mathlcal{k}}=t}$ nor $f_{\scaleobj{1.2}{\boldsymbol{{\mathrm{x}}}}|\mathlarger{\mathlarger{\omega}}=1,\Psi_{\mathlcal{k}}=t}$ is known or specified by the researcher, and they are substituted with their respective estimates: $\widehat{f}_{\boldsymbol{{\mathrm{x}}}|\mathlarger{\mathlarger{\omega}}=0}^{t}$ and $\widehat{f}_{\boldsymbol{{\mathrm{x}}}|\mathlarger{\mathlarger{\omega}}=1}^{t}$, as described in (ref), to follow. Similarly, the probabilities $\mathrm{Pr}\left( {\mathlarger{\mathlarger{\omega}}=1|\Psi_{\mathlcal{k}}=t} \right)$ and $\mathrm{Pr}\left( {\mathlarger{\mathlarger{\omega}}=0|\Psi_{\mathlcal{k}}=t} \right)$ are replaced by their estimates $\mathfrak{p}_t$ and $1-\mathfrak{p}_t$, respectively. Thus,

where $N_{t}^E$ is the cardinality (number of elements) of the set $\mathcal{O}_{p_{(t)}}$.

Using the Parzen-Rosenblatt rosenblatt1956remarks,parzen1962estimation window method for a nonparametric density estimation given a $L_{\boldsymbol{{x}}}\times L_{\boldsymbol{{x}}}$ bandwidth matrix $\mathcal{H}$,\footnote{The multivariate gaussian kernel density estimator is employed due to its applicability to multivariate data. Unlike in the case of semiparametric estimation, in the case of nonparametric estimation there is a “protocol” for finding the optimal bandwidth for instance, scott1991feasibility's rule. } we denote a conditional density estimator of the random variable vector $\boldsymbol{{\mathrm{x}}}\in \mathbb{R}^{L_{\boldsymbol{{x}}}}$ given the opinion $\omega\in\left \{{0,1} \right\}$ and an estimated membership in latent class $t$:

where

is a subset of $\left \{{\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_e} \right\}_{e=1}^{N^E}$, consisting only of the observed characteristics related to experts assigned to latent class $t$ with the opinion $\omega$ and $N_{t,\omega}^E$ is the cardinality (number of elements) of the set $\mathcal{O}_{\boldsymbol{{x}}_{(t)}}^{\omega}$. The determinant of $\mathcal{H}$ is $\left| {\mathcal{H}} \right|$. The main idea behind the mapping from an expert set to a sequence of opinions is to take advantage of auxiliary data (e.g., survey data, training data and the like), in which each data point depicts an opinion of a specific expert. Averaging the opinions in each reference group obtained from (ref) generates the share of participants belonging to that reference group. Thus, the best forecast, resulting from the various reference groups is a refinement of the wisdom of crowd (galton1907vox,surowiecki2005wisdom,prelec2017solution,hummel1996statistical,budescu2014identifying). The type and number of reference groups are unobserved and are estimated by the posterior class membership probability density function.

The proposed implementation relies on the utilization of two data sets: (i) a survey data set consisting of experts' opinions $\left \{{\mathcal{O}_{p_{(t)}}} \right\}_{t=1}^{\mathlcal{k}}$, with the observed group-level covariates; and (ii) a truncated data set consisting of both individual-level covariates as well as group-level covariates. The proposed procedure is closely related to the similarity-based classification, which is referred to as “nearest neighbor algorithm”, in the field of machine learning (e.g., chen2009similarity). Nearest neighbor algorithm assigns labels in the truncated (test) data set based on the similarities between this data set and the non-truncated labeled (training) data set. However, in the present case, the labels are unobserved not only in the truncated data set, but in both data sets.\footnote{This phenomenon is termed “unlabeled data” in the field of machine learning.} Thus, the purpose of the survey data set is to estimate the posterior distribution function in order to fit the labels in the truncated data set.

Estimating the substantive equation

We formulate the estimation procedure in terms of a non-linear least squares (NLS) minimization. Although the substantive equation is a linear function of its covariates, it can be reformulated as a partially linear single-index model in order to correct for the endogenous selection bias. The single-index modeling draws on the johnson1984extensions Lemma, alleviating the complexity present in high-dimension covariates space. In our model, the single index function is referred to as the bias term heckman1979sample and is constructed using (ref), by taking its conditional expectation, given the covariates and being a participant:

where $\boldsymbol{{w}}_i$ is a vector of the substantive equation's covariates; while $\boldsymbol{{x}}_i$ and $\boldsymbol{{z}}_i$ are specific group-level and individual-level characteristics, respectively.

The residual $\epsilon_i$ between $y_{1i,g}$ and its conditional expectation, given participation (ref) in the truncated data is constructed as:

Using (ref) and denoting $\mathbb{E}[\xi_{1i}|\mathrm{S}=1,\boldsymbol{{w}}_i,\boldsymbol{{z}}_i,\boldsymbol{{x}}_g,\psi_{i,g}]\equiv\mathcal{M}(\beta m^{\mathlcal{e}}(\boldsymbol{{x}}_g,\psi_{i,g})+\boldsymbol{{z}}_i^T\boldsymbol{{\eta}}+\left( {\boldsymbol{{x}}_{g}^{c}} \right)^T\boldsymbol{{\delta}})$ we arrive at the partially linear single-index model:

where $\boldsymbol{{x}}_{g}^{c}$ are the contextual covariates.

However, neither the function $m^{\mathlcal{e}}(\cdot)$ nor its $\psi_{i,g}$ argument is observed. Thus, they are substituted with their respective estimates $\widehat{m}_{\mathlcal{k}}^e(\cdot)$ and $\widehat{\psi}_{i,g}$ given $\mathlcal{k}$ possible latent classes (labels), obtained from the survey data or any other auxiliary data. The former is constructed using (ref), which is a refinement of the vox populi (average forecast) mechanism (see section (ref)):

where $N_{\widehat{\psi}_{i,g}}^E$ is the cardinality (number of elements) of the set $\mathcal{O}_{p_{(\widehat{\psi}_{i,g})}}$ and $\widehat{\psi}_{i,g}=\underset{t}{\arg\max}\hspace{0.5em}\mathfrak{P}_t(\boldsymbol{{x}}_g,\scaleobj{0.8}{\mathlcal{Y}}_i)$.

Our objective is to estimate the substantive equation (ref), which includes the function $\mathcal{M}(.)$ as an additional covariate controlling for the endogenous sample selection. However, the function $\mathcal{M}(.)$ in (ref) is unknown and has to be approximated. In the next section we attend to this issue.

Transformation of both Cosine and Fourier series for unknown \\ functions estimation

The function $\mathcal{M}(.)$ in (ref) is approximated using its conditional moment expansion by employing either The Cosine or The Fourier sequence. The Cosine sequence requires that the support of the index variable in (ref) will be on the $[0,1]$ domain, while the Fourier sequence requires that the support will be on the $[-1,1]$ domain.\footnote{Fourier series decomposes a “periodic” signal into a sum of an infinite number of harmonics (sine and cosine functions) of different frequencies and amplitudes, while Fourier transform decomposes a “non-periodic” signal into an infinite number of harmonics having different frequencies and amplitudes. } This assumption does not entail loss of generality, because it is satisfied by utilizing a different monotone transformation function on the index variable horowitz2014adaptive in each one of the Cosine and Fourier series. The series generated by the transformation is referred to as a transformed Cosine (or Fourier) series.

In the case of the (transformed) Cosine sequence the conditional moment expansion of $\mathcal{M}(.)$ is denoted by $\widehat{\mathcal{M}}^c(\mathlcal{b};\boldsymbol{{\vartheta_{c}}})$ and is defined $\forall\mathlcal{b}\in\mathbb{R}$ as:

where $\varphi(.)$ is some known, arbitrarily chosen, strictly monotonic twice differentiable mapping $\mathbb{R}\mapsto(0,1)$, $\boldsymbol{{\vartheta_{c}}}\equiv\left \{{\alpha_{c},\boldsymbol{{\tau^c}}} \right\}$ and $\boldsymbol{{\tau^c}}\equiv[\tau_{1}^c,...,\tau_{\mathcal{K}}^c]$, with $\mathcal{K}$ being the number of elements in the expansion.

Similarly, in the case of the (transformed) Fourier sequence the conditional moment expansion of $\mathcal{M}(.)$ is denoted by $\widehat{\mathcal{M}}^f(\mathlcal{b};\boldsymbol{{\vartheta_{c}}})$ and is defined $\forall\mathlcal{b}\in\mathbb{R}$ as:

where $\zeta(.)$ is some known, arbitrarily chosen, strictly monotonic twice differentiable mapping $\mathbb{R}\mapsto(-1,1)$,\footnote{The main drawback of Fourier series, however, is the requirement of the approximated function to be periodic on a bounded interval. This is problematic, as we are interested in approximating a non-periodic function defined on an unbounded interval. To alleviate this problem, we use monotonic mapping of the function's argument from the real line to the [-1,1] domain to make it periodic only at infinity and bounded on this domain. The aforementioned transformation results in enhanced accuracy of the estimates. due to the flexibility of Fourier series estimator, without being restricted to the family of periodic functions.} $\boldsymbol{{\vartheta_{f}}}\equiv\left \{{\alpha_{f},\boldsymbol{{\tau_{1}^f}},\boldsymbol{{\tau_{2}^f}}} \right\}$ and $\boldsymbol{{\tau_{m}^f}}\equiv[\tau_{m_1}^f,...,\tau_{m_{\mathcal{K}}}^f]$, $m=1,2$ representing Sine or Cosine, respectively.

For brevity, we denote the parameter vector $\boldsymbol{{\theta}}^*\equiv\left[ {\boldsymbol{{\theta}}^T,\boldsymbol{{\eta}}^T,\boldsymbol{{\delta}}^T,\beta} \right]^T$. Following racine2014oxford, given the non-linear function $\widehat{\mathcal{M}}^{\mathcal{G}}(\mathlcal{b};\boldsymbol{{\vartheta_{\mathcal{G}}}})$ with $\mathcal{G}\in\left \{{c,f} \right\}$ an index model can be estimated as follows:

where $y_{1i,g}$ is the substantive equation's dependent variable; $\mathcal{K}$ is the number of elements in the expansion; $\boldsymbol{{w_i}}$ and $\boldsymbol{{\theta}}$ stand for the covariates set and the parameter set, respectively in the linear part of the substantive equation. Note that the combination in (ref) of the linear component $\boldsymbol{{w}}_i^T\boldsymbol{{\theta}}$ and the non-linear component $\mathcal{M}(\cdot)$ implies partial linearity of the model.\footnote{This is where we depart from racine2014oxford, who introduce only the non-linear component, as they did not deal with truncation.}

We require that the expectation of the objective function in (ref) is finite for all values of the parameters $\boldsymbol{{(\theta^*,\vartheta_{\mathcal{G}})}}$\footnote{This assumption can be relaxed using a positive weight function $\mathcal{K}(x)$ on $(0,\infty)$ in the nonlinear minimization (see, racine2014oxford).} that is,

Next we have to modify (ref) and accommodate it for the presence of latent reference groups. This is done by introducing a penalization into the model.

The optimal number of latent reference groups

In practice, the number of latent classes (labels) in the truncated data is unknown. Arbitrarily choosing the number of latent classes may amount to misspecification.\footnote{A non-feasible solution is to assume that any individual observation is its own advisor (reference group) based on his past experience. This is problematic (unless an auxiliary data set with historical individual level participation probabilities is accessible), as the individual data consists of participants only, and consequently one cannot estimate the probability to participate for a specific data-point using only one observation, which is the participant herself.} To alleviate probable misspecification, we propose an estimation procedure generating the “optimal” number of latent classes to fit the correct model without arbitrarily assuming the number of reference groups. This procedure specifies the participation decision that best fits the data generation process in the truncated data, which is related to some criterion function (to be defined in (ref) to follow). The aforementioned participation decision is chosen from a menu consisting of selection equations differentiated by $\mathlcal{k}$ number of available (latent) reference groups. This is achieved by minimizing an additive penalized objective function $\Upsilon(\boldsymbol{{\varphi}})$ racine2014oxford:

where $\boldsymbol{{\varphi}}$ is a vector of estimated parameters and $\lambda_n$ is a tuning parameter.\footnote{When $\lambda_n$ approaches zero the penalty function is not effective, leading to the parameter estimates that would have been obtained without penalization.}

We note that increasing the number of reference groups decreases the model bias, due to enhanced information (explanatory ability), however, it is at a cost of higher variance in the model (low accuracy). To overcome this bias-variance trade-off, a penalization procedure is applied, as is depicted by the penalty function in (ref). The penalized regression is also termed “sparse regression”, where “sparsity” implies that only a small fraction of the predictor variables has an influence on the dependent variable friedman2012fast. These regression methods are intended to find the subset of the most influent predictors by shrinking down the parameter estimates toward zero and reducing the number of non-zero parameter estimates.

The most popular choice of loss-functions are Mean Squared Error (MSE), negative log-likelihood and profiled least squares. In our case, we employ the Mean Squared Error (MSE) loss function to be consistent with the nonlinear least squares problem depicted in (ref). In order to select the optimal number of latent reference groups, we use the SCAD penalty function, as it nests the LASSO as a special case, defined as:\footnote{ The SCAD performs well in partially linear index models racine2014oxford,liang2010estimation.}

where $a>2$ is a constant. For practical use we set $a=3.7$ racine2014oxford.\footnote{This has been shown to facilitate computation time.}

Next, we present an algorithm for determining the optimal number of latent reference groups, employing the SCAD penalty function.

Estimating the number of latent reference groups

We construct a sequence of functions $\left \{{\widehat{m}_{\mathlcal{k}}^e(\boldsymbol{{x}}_g,\psi_{i,g})} \right\}_{\mathlcal{k}=1}^{\tilde{G}_{\max}}$, where $\tilde{G}_{\max}\in\mathbb{N}$ is the largest latent labels (outcomes) number, $\widehat{m}_{\mathlcal{k}}^e(\boldsymbol{{x}}_g,\psi_{i,g})$ is the conditional participation probability, given $k\le \tilde{G}_{\max}$ latent labels, the observed group's characteristics $\boldsymbol{{x}}_g$ and belonging to label (being a member in class) $\psi_{i,g}\in\left \{{1,...,k} \right\}$.\footnote{In a recent contribution diebold2017beating also utilize a penalty function to combine forecasts. However, they utilize the LASSO penalty function which is restrictive in that it forces most of the covariates to have zero coefficients, instead of allowing for a combination of the covariates to be utilized like the SCAD penalty function employed here. } For brevity, we define $\widehat{\rho}_{\mathlcal{k}_{i,g}}\equiv\widehat{m}_{\mathlcal{k}}^e(\boldsymbol{{x}}_g,\widehat{\psi}_{i,g})$, where $\widehat{\psi}_{i,g}$ is the estimated label membership for the $i$'th observation, given $\mathlcal{k}$ latent classes. The partially linear single index regression is represented as:

where $\boldsymbol{{\varphi}}=\left \{{\boldsymbol{{\theta,\beta,\eta,\delta}},\boldsymbol{{\vartheta_{\mathcal{G}}}}} \right\}$, such that $\boldsymbol{{\beta}}\equiv\left \{{\beta_1,...,\beta_{\tilde{G}}} \right\}$.

In the first step, a solution path $\boldsymbol{{\varphi}}_{\lambda_n}=\scaleobj{1.2}{\left\{\right.}{\boldsymbol{{\theta_{\lambda_n},\beta_{\lambda_n},\eta_{\lambda_n}}}}$, ${\beta_{\delta_{\lambda_n}},\boldsymbol{{\vartheta_{\mathcal{G}}}}_{\lambda_n}}\scaleobj{1.2}{\left.\right\}}$ indexed by a tuning parameter, $\lambda_n$, is estimated as a penalized partially linear single index model liang2010estimation:\footnote{Unlike the (penalized) partially linear single index model (PLSIM) estimation procedure introduced by liang2010estimation which utilizes kernel estimator (suffering from bandwidth selection consideration) to approximate the unknown function $\mathcal{M}$, we use a Sieve estimator.}

where $\boldsymbol{{\beta_{\lambda_n}}}\equiv\left[ {\beta_{1_{\lambda_n}},...,\beta_{J_{\lambda_n}}} \right]^T$, $\boldsymbol{{y_{_1}}}=[\boldsymbol{{y_{_{1,1}}^T}},...,\boldsymbol{{y_{_{1,G}}^T}}]^T$ such that $\boldsymbol{{y_{_{1,g}}}}=[y_{_{11,g}},...,y_{_{1n_{g},g}}]^T$ and the $\left\lVert\boldsymbol{{\cdotp}}\right\rVert_{2}$ is the usual $\ell_2$ (Euclidean) norm.\footnote{The $\ell_p$ norm definition is:

}

In the second step, a criterion $\mathcal{C}_p$ is computed for the solution path $\boldsymbol{{\widehat{\varphi}}}_{\lambda_n}$. The conventionally chosen criterion is BIC (Bayesian information criterion greene2003latent) computed as:

where $\text{MSE}(\lambda_n)=n^{-1}\sum_{i=1}^{n} \left( {y_{1i,g}-\mathcal{F}_{i,g}(\boldsymbol{{\widehat{\varphi}}}_{\lambda_n})} \right)^2$ and $df_{\lambda_n}$ is the number of non-zero coefficients in $\boldsymbol{{\widehat{\varphi}}}_{\lambda_n}$.

The algorithm for finding the correct model requires estimating (ref) repeatedly, each time given a different tuning parameter value $\lambda_n$, and computing $\text{MSE}(\lambda_n)$ in order to find $\lambda_n$ which minimizes (ref).

Estimating (ref) involves the utilization of a non-convex penalty function optimization, which enhances computational complexity. To alleviate this complexity and without loss of accuracy, we transform the optimization problem into a constrained one with a convex penalty function figueiredo2007gradient.

Thus, we introduce the parameter vectors $\boldsymbol{{\beta_{\lambda_n}^{+}}}\equiv\left[ {\beta_{1_{\lambda_n}}^{+},...,\beta_{J_{\lambda_n}}^{+}} \right]^T$ and $\boldsymbol{{\beta_{\lambda_n}^{-}}}\equiv\left[ {\beta_{1_{\lambda_n}}^{-},...,\beta_{J_{\lambda_n}}^{-}} \right]^T$ where $\beta_{k_{\lambda_n}}^{+}=\max\left \{{0,\beta_{k_{\lambda_n}}} \right\}$ and $\beta_{k_{\lambda_n}}^{-}=\max\left \{{0,-\beta_{k_{\lambda_n}}} \right\}$ $\forall k$ and make the following substitution:

Using $\boldsymbol{{\beta_{\lambda_n}^{+}}}$ and $\boldsymbol{{\beta_{\lambda_n}^{-}}}$ the optimization becomes:

where the modified solution path is $\boldsymbol{{\varphi}}_{\lambda_n}^*=\scaleobj{1.2}{\left\{\right.}{\boldsymbol{{\theta_{\lambda_n},\boldsymbol{{\beta_{\lambda_n}^{+}}}}}a_n}^{+}}}}}$, ${\boldsymbol{{\beta_{\lambda_n}^{-}}},\boldsymbol{{\eta_{_{\lambda_n}}}},\boldsymbol{{\delta_{\lambda_n}}},\boldsymbol{{\alpha_{_{\lambda_n}}}}}\scaleobj{1.2}{\left.\right\}}$.

Note that the the sum of squares term in (ref) is unaffected, if we set $\boldsymbol{{\beta_{\lambda_n}^{+}}}\longleftarrow\boldsymbol{{\beta_{\lambda_n}^{+}}}+\boldsymbol{{s}}$ and $\boldsymbol{{\beta_{\lambda_n}^{-}}}\longleftarrow\boldsymbol{{\beta_{\lambda_n}^{-}}}+\boldsymbol{{s}}$ $\forall\boldsymbol{{s}}\ge \boldsymbol{{0}}$, because $\boldsymbol{{s}}$ is canceled out in (ref).\footnote{$\boldsymbol{{s}}$ cannot contain negative elements, because the set $(\beta_{k_{\lambda_n}}^{+},\beta_{k_{\lambda_n}}^{-})=(\beta_{k_{\lambda_n}},0)$ implies $\beta_{k_{\lambda_n}}>0$, while the set $(\beta_{k_{\lambda_n}}^{+},\beta_{k_{\lambda_n}}^{-})=(0,-\beta_{k_{\lambda_n}})$ implies $\beta_{k_{\lambda_n}}<0$. The intuition being that if $\boldsymbol{{s}}<0$, the requirements $\boldsymbol{{\beta_{\lambda_n}^{+}}}+\boldsymbol{{s}}\ge 0$ and $\boldsymbol{{\beta_{\lambda_n}^{-}}}+\boldsymbol{{s}}\ge 0$ are not satisfied. If $\boldsymbol{{s}}>0$ the penalty function is not minimized. } However, the argument in the penalty function term increases by $2\boldsymbol{{s}}$. As a result, $\boldsymbol{{s}}=0$ minimizes the penalty function, implying that the solution of problem (ref) for a given $k$ is either $\beta_{k_{\lambda_n}}^{+}=0$ or $\beta_{k_{\lambda_n}}^{-}=0$. Problem (ref) is equivalent to the original problem (ref), where $\left| {\beta_{k_{\lambda_n}}} \right|=\beta_{k_{\lambda_n}}^{+}+\beta_{k_{\lambda_n}}^{-}$ and $\beta_{k_{\lambda_n}}=\beta_{k_{\lambda_n}}^{+}-\beta_{k_{\lambda_n}}^{-}$ $\forall k$. The aforementioned argument points to the possibility of using simple constrained convex penalty function algorithms for the estimation of the optimal number of reference groups; this is embedded in the selection equation, which is affected by the number of reference groups. Technical details appear in Appendix (ref).

Simulations

We examine our truncated selection model's performance in the presence of various latent classes, capturing the unobserved characteristics of each datum. A sequence $\left \{{(\Lambda_{k}^S,\Lambda_{k}^T)} \right\}_{k=1}^{10,000}$ consisting of $10,000$ elements is generated. The $k$'th element is composed of a survey data set and a truncated data set denoted by $\Lambda_k^S$ and $\Lambda_k^T$, respectively. For simplicity, the data sets are generated using three latent classes.

Next we discuss the data generation process (DGP) used to construct these distribution functions.

Data Generation Process (DGP)

Let $\mathrm{z}$ be a continuous random variable, such that a given realization of this random variable represents specific individual level characteristics. Our objective is to characterize a sequence of distribution functions $\left \{{\mathcal{D}_{\mathrm{z}|\boldsymbol{{x}}_g}} \right\}_{g=1}^{G}$. Each is a conditional distribution function of $\mathrm{z}$, given a specific realization of the group-level observed characteristics, $\boldsymbol{{x}}_g$. These distribution functions are not restricted to being unimodal or symmetric (e.g., the normal distribution function).\footnote{Unlike the Monte Carlo simulations in breunig2017nonparametric for censored sample selection models implemented by using normally distributed disturbances, we consider a truncated sample selection model characterized by non-normally distributed disturbances.} By employing such an algorithm, each datum in the truncated data set to be generated is a random draw from its group-specific distribution function. We arbitrarily set $G=2,000,000$ indicating the number of distribution functions in the sequence. The conditional density function of $\mathrm{z}$ given $\boldsymbol{{x}}_g$ is denoted by $\mathlcal{d}_{z|\boldsymbol{{x}}_g}(z|\boldsymbol{{x}}_g)$ and satisfies $\forall t=1,...,\tilde{G}$:

where $m(\boldsymbol{{x}}_g,t)$ is the actual participants' share given $\boldsymbol{{x}}_g$ and being a member in class $t$, and $\boldsymbol{{x}}_g^c$ is a subset of $\boldsymbol{{x}}_g$, consisting of the contextual covariates only.

Finding a density function $\mathlcal{d}_{z|\boldsymbol{{x}}_g}(z|\boldsymbol{{x}}_g)$ that satisfies (ref) is computationally cumbersome, due to the presence of the integral. In order to facilitate the computation process, this density is expressed as a finite mixture of arbitrarily chosen continuous density functions, such that only the weights (mixture coefficients) are required to uncover.

Thus, let $\left \{{\phi_l(.): l=1,...,L} \right\}$ be an arbitrary set of distinct continuous probability density functions on the real line\footnote{We utilize a mixture consisting of various density functions, including normal, gamma and log-normal. Each density function is characterized by a unique set of parameters.} and $\left \{{\Phi_l(.): l=1,...,L} \right\}$ be the corresponding distribution functions. The mixture of the probability density functions using the weights $w_l$, satisfying $\sum_{l=1}^{L}w_l=1$, is defined as follows:

and the mixture of the distribution functions is

Denote $\mathfrak{q}(p,z)\equiv 1-F_{\xi_{2}}(\alpha+\beta p +\left( {\boldsymbol{{x}}_g^c} \right)^T\delta + \eta z)$. The sequence of optimal weights, $\left \{{w_l} \right\}_{l=1}^L$, consists of $L$ elements such that $\sum_{l=1}^{L}w_l=1$ and $0\le w_l \le 1$. The following condition must hold for all $t=1,...,\tilde{G}$:

A matrix $\mathbf{M}_{\scaleto{\mathcal{P}}{4pt}}$ of size $\tilde{G}\times L$ and a vector $\mathbf{w}$ of size $L\times 1$ are defined as:

These optimal weights, which solve (ref), can be obtained as a solution to the following minimization problem:

where $\mathcal{P}=[m(\boldsymbol{{x}}_g,1),...,m(\boldsymbol{{x}}_g,\tilde{G})]^T$ is a $\tilde{G}\times 1$ vector, $\mathbf{M_{\scaleto{\mathcal{P}'}{4pt}}}$ and $\mathbf{M_{\scaleto{\mathcal{P}''}{4pt}}}$ are matrices constructed in a similar fashion to $\mathbf{M_{\scaleto{\mathcal{P}}{4pt}}}$ (described in (ref)) to ensure that the expected participants' shares vector, $\mathcal{P}$, is the unique solution for (ref) such that any other expected participant share depicted in either vector $\mathcal{P}'$ of size $J'\times 1$ or vector $\mathcal{P}''$ of size $J''\times 1$ will not constitute a solution.

The main idea is that for any vector $\mathcal{P}$ consisting of participants' shares (calculated for a given observed group's characteristics), we match a distribution function of individual characteristics, such that (ref) is satisfied.\footnote{By construction, the generated distribution function leads to a unique solution characterized by vector $\mathcal{P}$.} Formally, the $k$'th data consisting of $N$ observations is generated by using a sequence $\left \{{\boldsymbol{{x_i}}} \right\}_{i=1}^N$ of randomly drawn vectors from a joint distribution function $F_{\mathrm{\boldsymbol{{x}}}}$ (as will be described in (ref) to follow) and a sequence of distribution functions $\left \{{\mathcal{D}_{z|\boldsymbol{{x_i}}}(z|\boldsymbol{{x_i}})} \right\}_{i=1}^N$. Each data point $\boldsymbol{{z_i}}$ of individual level characteristics is a realization of a unique random variable $\boldsymbol{{\mathrm{z|x}=x_i}}$, drawn from $\mathcal{D}_{z|\boldsymbol{{x_i}}}(z|\boldsymbol{{x_i}})$, where $\boldsymbol{{x}}_i=[x_{i}, x_{i}^c]$ and represents the observed group's characteristics of the $i$'th observation.\footnote{$z_i$ is a single covariate, so we use a univariate distribution function. However, in cases where it is a covariate vector, the index function $\boldsymbol{{z}}_i^T\boldsymbol{{\eta}}$ can be treated as a random realization from a univariate distribution function $\mathcal{D}_{\boldsymbol{{z}}_i^T\boldsymbol{{\eta}}|\boldsymbol{{x_i}}}(\boldsymbol{{z}}_i^T\boldsymbol{{\eta}}|\boldsymbol{{x_i}})$.}

Based on (ref), each of the distribution functions to be found is required to satisfy a restriction concerning a specific participant shares vector. Therefore, these shares must be known in the data generation process.\footnote{Thus, we define the $t$'th element of $\mathcal{P}(\boldsymbol{{x}}_g)$ as:

where $\left( {\mathlcal{p}_{lt}, \mathlcal{p}_{ht}} \right)$ is the range of the conditional participation probability given a membership in class $t$. The ranges are arbitrarily determined to be: $\left( {\mathlcal{p}_{l1}, \mathlcal{p}_{h1}} \right)=\left( {0.05, 0.4} \right)$, $\left( {\mathlcal{p}_{l2}, \mathlcal{p}_{h2}} \right)=\left( {0.4, 0.75} \right)$ and $\left( {\mathlcal{p}_{l3}, \mathlcal{p}_{h3}} \right)=\left( {0.65, 0.95} \right)$. These numbers enable us to verify the model performance in cases where the participants' share is a non-smooth function of the group's observed characteristics ($\boldsymbol{{x}}_g$). This non-smoothness stems from the presence of latent classes, which are determined as a function of $\boldsymbol{{x}}_g$ (due to non-random assignment). The cumulative standard normal distribution function is denoted by $\Phi_{_{\mathcal{N}}}$. }

Generation of survey and truncated data sets utilizing latent classes

The group's characteristics covariate vector $\boldsymbol{{x}}_g=[x_{g}, x_{g}^c]$ is a realization of the random variables vector $\boldsymbol{{\mathrm{x}}}=[\mathrm{x,x^c}]$, which is jointly distributed $F_{x}$. For simplicity we characterize $F_x$ as follows:

where $\mathcal{N}_2$ denotes the bivariate normal distribution function.

The classes and a frequency table consisting of manifest variables (depicted in Table (ref) to follow) are generated in $R$, using the latent classes packages $`poLCA'$ and $`SimCorMultRes'$. For each $k=1,...,10,000$, two data sets are generated: (i) a survey data set $\Lambda_k^S=\left \{{\varphi_i^*,\chi_{\varphi_i}} \right\}_{i=1}^{N^E}$ with $\varphi_i=(\varphi_i^*,\psi_i^E)\in\Lambda$ (constructed as depicted in section (ref)) consisting of a sequence $\left \{{\varphi_i^*} \right\}_{i=1}^{N^E}$, in which $\varphi_i^*=(\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_i,\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}}^E_i)$ are the experts' observed characteristics and $N^E=10,000$; (ii) a truncated data set $\Lambda_k^T=\left \{{(\boldsymbol{{x}},\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}},y_{1}^*,z)|(\boldsymbol{{x}},\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}},y_{1}^*,y_{2}^*,z,\psi)\in \Lambda_k^C,y_{2}^*\ge0} \right\}$, where $\Lambda_k^C=\left \{{\boldsymbol{{x}}_i,\scaleobj{0.8}{\boldsymbol{{\mathlcal{{Y}}}}}_{i},y_{1i}^*,y_{2i}^*,z_i,\psi_i} \right\}_{i=1}^{N}$ denotes the complete (non-truncated) data set. The number of observations in the truncated data set is denoted by $n_{k}$, which is the cardinality of the set $\Lambda_k^T$.\footnote{The truncated data set is produced by generating a complete (non-truncated) data set and keeping only the observations that satisfies the selection equation.} The set $\Lambda_k^T$ consists of classes, manifest variables and observed group characteristics $\boldsymbol{{x}}$ randomly and independently drawn from (ref). These characterize the group's observed characteristics and individual-level covariates, as denoted by the sequence $\left \{{z_i} \right\}_{i=1}^n$. Using the estimated posteriori classes density function and the survey data set, the predicted class, $\widehat{\psi}_{i,g}$, for each observation $i$ is calculated in the truncated data set (given the manifest variables and the group's observed characteristics\footnote{The observed group's characteristics affect the prior class membership assignment probabilities.}). These predicted classes are intended to be used later, in the estimation stage, and not in the data generation process. The selection model's equation will be generated by using the true classes, as will be described in section (ref).

The characteristics $\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_i=[x_i^E,(x_i^c)^E]$ are randomly and independently drawn from (ref). The $i$'th observation's latent class is generated as a random realization from (ref), which is a function of $\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_i$.\footnote{Due to the non-random assignment, the $i$'th observation's latent class membership depends on $\scaleobj{0.8}{\boldsymbol{{\mathlcal{{X}}}}}^E_i$.} The manifest variables are determined using the frequencies in Table (ref) and are generated for each $i\in{1,...,N^E}$.

Table (ref) exhibits the construction of the manifest categorical variables in the survey data set, which are denoted by $\mathlcal{y}_{_1}$,...,$\mathlcal{y}_{_7}$. These variables enable the recovery of the latent classes. Recovery is feasible due to the conditional independence (given the class) property. It follows that the probability of observing a specific response (the alternative being chosen) in $\mathlcal{y}_{_k}$ is independent of the response to $\mathlcal{y}_{_l}$ for all $k\ne l$, given the class membership.

Next we construct the selection model's joint disturbances distribution function, in order to examine our model's performance in cases of a non-standard distribution function of the disturbances (such as the normal distribution function).

The disturbances' joint distribution function

Each pair of disturbances $\left \{{\xi_{1i},\xi_{2i}} \right\}$ is randomly and independently drawn from $F_{\xi_1,\xi_2}$, which is the joint distribution function of the substantive and participation equations' disturbances. The aforementioned joint density function consists of two components: a Copula function\footnote{According to Sklar's Theorem sklar1959fonctions, any continuous joint distribution function can be characterized by a set of marginal distribution functions and a joint distribution function determining the dependence structure which is referred to as a Copula function.} characterizing the disturbances' dependence structure and two marginal distribution functions $F_{\xi_1}$ and $F_{\xi_2}$ for the substantive equation and selection equation, respectively. In order to verify our model's performance in the presence of random disturbances' distribution functions which are not restricted to the family of symmetric and unimodal distribution functions, each one of these two disturbances is marginally distributed according to a mixture of three different distribution functions: (i) a normal distribution function with expectation and standard deviation parameters $(\mu, \sigma_a)$, denoted by $\mathcal{N}(\mu,\sigma_{a}^2)$; (ii) a normal distribution function with expectation and standard deviation parameters $(-\mu, \sigma_b)$, denoted by $\mathcal{N}(-\mu,\sigma_{b}^2)$; (iii) a gamma distribution function with scale and shape parameters $(\mu\varphi,\varphi)$, denoted by $\Gamma_{\scaleto{\mathcal{\text{Gamma}}}{4pt}}\left( {\mu\varphi, \varphi} \right)$\footnote{The scale and shape parameters implies that the expectation and standard deviation parameters are $(\mu ,\sqrt{\mu/\varphi})$.}. This mixture distribution function is defined as:

where $\mathbb{E}\left[ {\xi_{j}} \right]=0$, $j=1,2$.

The parameters set $(\mu,\sigma_a,\sigma_b,\varphi)=(4,3.5,2.5,2)$ is arbitrarily chosen. Due to its simplicity, the Clayton Copula, with a degree of dependence parameter $4$ (to assure that the disturbances are highly correlated) is used for controlling the dependence structure.

Constructing the selection model's equations

In the last step, we construct the latent selection equation's dependent variable $y_{2i,g}^*$ for $i=1,...,N$:

where the $i$'th observation's true class membership is denoted by a categorical variable\\ $\psi_{i,g}\in\left \{{1,2,3} \right\}$ and $\boldsymbol{{x}}_g=[x_g,x_g^c]$.

In a similar fashion, we construct the substantive equation's dependent variable:

where any covariates pair $(w_i, D_i)$ is an independent realization of a random variable vector $(\mathrm{w,D})|z_i$, conditionally distributed given $z_i$, such that $\mathrm{w}|z_i$ and $\mathrm{D}|z_i$ are independent.\footnote{This assumption is termed a conditional independence given $z_i$ and is imposed to reduce the complexity of the data generation process.} These conditional random variables are distributed according to a normal and a Bernoulli probability distribution function, respectively. We arbitrarily set $(\alpha,\beta,\delta)=(-30, 25, 1.5)$ and $(\theta_1, \theta_2)=(2,4)$.

Denote the selection variable by an indicator function $y_{2i,g}=I(y_{2i,g}^*\ge 0)$ or alternatively, calculate the $i$'th observation's probability of being selected $p_i=\Pr(\mathlarger{\mathlarger{\omega}}_{i,g}=1|x_{g},x_{g}^c,\psi_{i,g},z_i)$, which is a function of $(x_{g},x_{g}^c,\psi_{i,g},z_i)$ and defined for $i=1,...,N$ as:

Let $\left \{{u_1,...,u_N} \right\}$ be a sequence of continuous and independent uniform random variables on the support $[0,1]$. The indicator variable for the $i$'th observation is:

that is, $y_{2i,g}$ is the realization of a Bernoulli distributed random variable with a probability of success, $p_{i,g}$, with $p_{i,g}$ the probability of observation $i$, in reference group $g$, to be observed in the truncated data.

The final truncated data set consists of two sequences of the self-selected observation (satisfying $y_{2i}=1$). The truncated data set consists of 50$\%$ of the observations.\footnote{Following arabmazar1982investigation, employing a parametric (censored or truncated) sample selection model and misspecifying the random disturbances to be joint normal distributed might lead to bias in the estimates, such that its magnitude depends on the degree of censoring or truncation of the sample. They find evidence that the bias is substantial, especially for truncated samples that are 50 percent complete. Because our model is distribution-free, it is important to verify its performance given those conditions in which the parametric models underperformed. Thus, we use a truncated data set which is 50$\%$ complete.} The first includes the individual level covariates $\left \{{y_{1i}^*,z_i,w_i,D_i|y_{2i}=1} \right\}_{i=1}^N$, and the second includes the group's observed characteristics and the manifest categorical variables: $\left \{{x_{1i},x_{2i},\mathlcal{y}_{1i},..,\mathlcal{y}_{7i}|y_{2i}=1} \right\}_{i=1}^N$. Similarly, the survey data set consists of the sequence $\left \{{\bar{\mathcal{P}},x_{1i}^{\text{survey}},x_{2i}^{\text{survey}},\mathlcal{y}_{1i}^{\text{survey}},..,\mathlcal{y}_{7i}^{\text{survey}}} \right\}_{i=1}^N$.\footnote{Of course, the true class sequence in both the truncated data set $\left \{{c_i} \right\}_{i=1}^N$ as well as the survey data set $\left \{{c_i^\text{survey}} \right\}_{i=1}^N$ is unobserved. Thus, they are excluded from these data sets.} The aggregate participants' share in the entire population is denoted by $\bar{\mathcal{P}}$.

The main results regarding the estimates obtained using $10,000$ Monte Carlo simulations for sample sizes $N\in\left \{{2000,5000,8000,10000} \right\}$, as depicted in sections (ref)-(ref), are summarized in Tables (ref) and (ref) to follow.

The entries in Table (ref) indicate that for a sample size of $2000$ observations, the mean estimate of $\theta_1$ in the full sample is $2.000$, while in the truncated sample, without correction for bias is $2.5873$. Similarly, for the same sample size, the mean estimate of $\theta_2$ in the full sample is $4.000$, while in the truncated sample without correction for bias is $4.2916$. The estimates for $\theta_1$ and $\theta_2$ are remained biased (upward) as the sample size increases.

\newlength\q \setlength\q{\dimexpr .1\textwidth -2\tabcolsep} \newlength\qq \setlength\qq{\dimexpr .12\textwidth -2\tabcolsep}

Entries presented in Table (ref) indicate that the substantive equation's estimated parameters, $\widehat{\theta}_1=2.0083$ and $\widehat{\theta}_2=3.9909$, in the refined Vox Populi model (given the correct number of latent classes) almost mimic the parameters estimates that would have been generated from a non-truncated sample, which are $\theta_1=2.000$ and $\theta_2=4.000$ respectively, for a sample size of 2,000 observations.\footnote{The parameter estimates start deteriorating below $2000$ observations in both the penalized and the non-penalized refined Vox Populi models, and thus are not presented in the table.} However, treating the truncated data as if it consist of a single class (a monolithic Vox Populi) in the presence of multiple latent classes, the substantive equation's parameters' estimates are $\widehat{\theta}_1=2.0509$ and $\widehat{\theta}_2=3.9083$, given a sample of 2,000 observations. The standard deviation for each estimated parameter is calculated over the estimates obtained from all the Monte Carlo simulations,\footnote{The standard deviations are calculated using the same methodology for each of the models, to make it easier to compare results from different regression models.} for each one of the refined and monolithic Vox Populi models. The means of the parameter estimates obtained for $\theta_1$ and $\theta_2$ in the refined Vox Populi model are $2.007$ and $3.9978$, respectively, for a sample size of 10,000 observations. While in the monolithic Vox Populi model, the means of the parameter estimates obtained for $\theta_1$ and $\theta_2$ in the selection equation, are $2.0388$ and $3.9157$, respectively, using the same sample size.\footnote{The nuisance parameters' estimates, consisting of the selection equation's estimated coefficients, will be supplied upon request.}

Entries presented in Table (ref) indicate that the substantive equation's estimated parameters obtained in the refined Vox Populi model, using a SCAD penalty function (in the absence of a prior knowledge regarding the number of latent classes). These estimated parameters, $\widehat{\theta}_1=2.0049$ and $\widehat{\theta}_2=4.0037$, almost mimic the parameters estimates that would have been generated from a non-truncated sample, which are $\widehat{\theta}_1=2.000$ and $\widehat{\theta}_2=4.000$, respectively, for a sample size of 5,000 observations. It is worth noting that the parameters estimates' accuracy is improved by employing the SCAD penalty function (in terms of proximity to the true parameters values) relative to the parameters' estimates obtained in the absence of a penalty function, which are $\widehat{\theta}_1=2.0333$ and $\widehat{\theta}_2=3.9518$, using the same sample size. For a given sample size, the estimated standard deviations are slightly smaller in the model without a penalty function relative to the model with the SCAD penalty function. This implies that the penalty function reduces the bias in the estimates at the cost of a minor increase in dispersion.

As reflected by entries in the above tables, correcting for endogenous truncation bias is accurately achieved by applying our semiparametric Sieve estimator, which embeds the notion of refined Vox-Populi decision making.

Conclusion

\raggedbottom The primary purpose of this paper is to correct for selectivity bias generated by endogenous truncation. Incorporating behavioral aspects from economics, psychology and management science to introduce cognition into the participation decision-making allows for endogeneity to take place, to be modeled and to be controlled. We treat each data point's truncation decision based on the decision made by its reference group's opinion space. To accomplish this, we refine the monolithic notion of vox populi (wisdom of the crowd) by treating the data as a mixture of reference groups. We offer a three-stage procedure to correct for this endogenous selectivity bias. In the first stage, latent classes analysis is employed to estimate the various reference groups' memberships based on results from an auxiliary survey data. In the second stage, estimates for the groups' participation decisions are obtained by averaging their respective group members opinions. In the third stage, a semiparametric truncated sample selection model is estimated, consisting of a substantive equation and a selection equation, in which the estimated group's participation decision is an additional covariate.

The number of reference groups is not arbitrarily imposed but rather estimated using the smoothly clipped absolute deviation (SCAD) penalization mechanism. Monte Carlo simulations involve 2,000,000 different distribution functions, which are not restricted to the unimodal symmetric family of distribution function. This practically generates 100 million realizations which are not i.i.d. They attest to a very high accuracy of the model, as depicted by the parameter estimates, which quite accurately mimic the true parameters. \appendices

Latent classes model's assumptions

In this section we impose the latent classes model assumptions.

Homogeneity (H) The core assumption in latent class analysis is that the population consists of a set of mutually exclusive and homogeneous subgroups called classes. The individuals within a sub-group are homogeneous in the sense that the probability for a particular response on a particular item depends only on the latent class to which the individual belongs.

Local Independence (LI) Local independence assumes that the observed manifest variables, $\mathlcal{y}_1,...,\mathlcal{y}_j$ are related only due to the latent class $\Psi_{\mathlcal{k}}$. Under this assumption, the joint probability of $\mathlcal{y}_1,...,\mathlcal{y}_j$ given $\Psi_{\mathlcal{k}}$ can be written as the product of probabilities of $Y_t$ given the latent class $t$.

Unidimensionality (U) The assumption of unidimensionality posits that the observed categorical variables Y are assumed to measure only one ability, attitude, trait, or attribute.

Monotonicity (M) To obtain stochastic ordering among the latent classes within an item, Croon (1991) proposed an ordinal latent class model by imposing inequality restrictions:

for all $j$ and $k$, and for all $t_1$ and $t_2$ such that $t_1<t_2$.

The penalty function

The penalty function $p_{\lambda_n}(\cdotp)$ in (ref) is decomposed as $p_{\lambda_n}(\cdotp)=\mathlcal{h}_1(\cdotp)-\mathlcal{h}_2(\cdotp)$ which is a difference of two convex functions:

Let $b^{(t)}$ be a parameter value obtained at iteration $t$. The best affine approximation of $\mathlcal{h}_2$ at $b^{(t)}$ is:

The penalty function is approximated using (ref) as:

It is worth noting the following equivalence which must be satisfied:

The algorithm is to solve iteratively the problem using the affine approximation in (ref):

where $t$ is the iteration number and $\nabla \mathlcal{h}_2$ is the gradient of $\mathlcal{h}_2$ evaluated at $\boldsymbol{{\beta_{\lambda_n}^{+,(t)}}}+\boldsymbol{{\beta_{\lambda_n}^{-,(t)}}}$, such that $\left( {\boldsymbol{{\beta_{\lambda_n}^{+,(t)}}},\boldsymbol{{\beta_{\lambda_n}^{-,(t)}}}} \right)\in \boldsymbol{{\varphi}}_{\lambda_n^*}^{(t)}$.

After decomposing the coefficient vector $\boldsymbol{{\beta}}$ to enable difference of convex functions (DC) programming formulation, problem (ref) can be formulated as a weighted LASSO problem:

where $\boldsymbol{{\tilde\lambda}}=n(\lambda_{n}-\nabla \mathlcal{h}_2)$ and its $j$th element is $\tilde\lambda_{j}$.

Let $\mathlarger{\mathlarger{\mathfrak{f}}}(\boldsymbol{{x}})\equiv \frac{1}{2}\Big\lVert{\boldsymbol{{y_{_1}}}-\mathcal{F}(\boldsymbol{{x}})}\Big\rVert_{2}^{2}$, the update rule to minimize (ref) is computed using a second-order approximation of $\mathlarger{\mathlarger{\mathfrak{f}}}(.)$ at $\boldsymbol{{\varphi}}_{\lambda_n}^{(k)}$ yang2013review:

To simplify the minimization problem in (ref) we let $\alpha I\cong\nabla^2\mathfrak{f}(\boldsymbol{{\varphi}}_{\lambda_n}^{({k})})$, the following approximation is used:

where $I$ is the identity matrix.

to get:

After canceling out the constant terms in (ref) (which are not function of $\boldsymbol{{\varphi}}_{\lambda_n}$), the minimization problem is reformulated as follows:

where $\boldsymbol{{u^{(t)}}}=\boldsymbol{{\varphi}}_{\lambda_n}^{(t)}-\frac{1}{\alpha_t}\nabla \mathcal{F}(\boldsymbol{{\varphi}}_{\lambda_n}^{(t)})$ and $u_j^{(t)}$ is its $j$'th element.

The algorithm for solving (ref) is to update each $j$ component in the parameter vector $\boldsymbol{{\varphi}}_{\lambda_n}$ using the well-known soft-threshold algorithm (yang2013review and yang2016sparse for linear and non-linear treatments respectively):\footnote{Since the $\ell_1$ norm (the weighted LASSO penalty function) is separable, the computation of $\boldsymbol{{\varphi}}_{\lambda_n}^{(k+1)}$ reduced to solve a one dimensional minimization problem for each of its components.}

where $\text{soft}(u,a)\equiv \mathrm{sign}(u)\max(\left| {u} \right|-a,0)$.

For any parameter $j$ which is not part of the penalization $\tilde\lambda_j=0$.

The modified Soft thresholding algorithm

Let define $p_{\lambda_n}^{\prime}(\cdot)$ as:

where $\mathrm{sign}(v)=1\left \{{v>0} \right\}-1\left \{{v<0} \right\}$ such that $1\left \{{\cdot} \right\}$ is an indicator function.

We characterize the second order approximation for the non-linear function in (.):

Taking the derivative with respect to $\boldsymbol{{\beta}}$ to obtain:

where $\nabla\left( {\sum_{k=1}^{\tilde{G}_{\max}}p_{\lambda_n}\left( {\left| {\beta_{k_{\lambda_n}}^{(t+1)}} \right|} \right)} \right) \right)}$ is a vector of size $\tilde{G}_{\max}\times 1$ and its $k$'th element is $p_{\lambda_n}^{\prime}(\beta_{k_{\lambda_n}}^{(t+1)})$.

Let $\mathcal{A}=\scaleobj{1.8}{\left\{\right.}{\left[ {\mathfrak{s}_1,...,\mathfrak{s}_{\tilde{G}_{\max}}} \right]^T| \hspace{1em}\mathfrak{s}_{k}\in\left \{{-1,0,1} \right\}},\\ {\hspace{1em} k=1,...,\tilde{G}_{\max}}\scaleobj{1.8}{\left.\right\}}$ be the set consisting of all possible signs for a real number vector of size $\tilde{G}_{\max}\times 1$. Using this set, we denote a sign operator $\mathfrak{S}:\mathbb{R}^{\tilde{G}_{\max}}\mapsto \mathcal{A}$. It follows that the $k$'th element of $\mathfrak{S}(\boldsymbol{{\beta}}^{(t+1)})$ is $\mathrm{sign}(\beta_k^{(t+1)})$ and the matrix representation of $\nabla\left( {\sum_{k=1}^{\tilde{G}_{\max}}p_{\lambda_n}\left( {\left| {\beta_{k_{\lambda_n}}^{(t+1)}} \right|} \right)} \right) \right)}$ is:\footnote{The notation $\circ$ represents the Hadamard product.}

where $\boldsymbol{{\mathcal{R}}}$ and $\boldsymbol{{\mathlcal{r}}}$ are a matrix of size $\tilde{G}_{\max}\times \tilde{G}_{\max}$ and a $\tilde{G}_{\max}\times 1$ vector, defined in (ref) and (ref), respectively.

The matrix representation depicted in (ref) implies that $\boldsymbol{{\beta}}$ can be isolated from (ref) to get:

However, the expression in (ref) depends on the sign operator which is a function of $\boldsymbol{{\beta}}_{\lambda_n}^{(t+1)}$ to alleviate the recursive nature of the formula for $\boldsymbol{{\beta}}_{\lambda_n}^{(t+1)}$ we substitute $\mathfrak{S}(\boldsymbol{{\beta}}^{(t+1)})$ with a vector of signs $\mathfrak{g}\in\mathcal{A}$. It worth noting that $\mathfrak{g}=\mathfrak{S}(\boldsymbol{{\beta}}_{\lambda_n}^{(t+1)})$ if and only if:

In cases where a solution $\mathfrak{g}$ to (ref) is found, the updated $\boldsymbol{{\beta}}$ is:

The justification for the proposed non-linear penalized regression estimation algorithm is based on a unified algorithm introduced by fan2001variable which optimizes various linear penalized regression problems via local quadratic approximations.

In cases where there is no solution, we find the largest subset of signs for which there is a solution to (ref), and set to zero all the rest of the parameter values as if we employed the original soft-thresholding algorithm.

Binary Response model assumptions

Our objective here is present the necessary conditions for identification of the binary response model unknown parameters. These necessary conditions are depicted in the following assumptions:\footnote{The assumptions are taken from brock2007identification.}

Based on Proposition 1 in brock2007identification, under assumptions (ref)-(ref), the parameters of the binary choice model are identified up to scale.

The contextual effect identification:

The assumption (ref) implies that the objective function can change across reference groups, but not within reference group. It also implies Homogeneity within reference group.

There are two core assumptions in latent classes analysis: Assumption 1: (Homogeneity) The population consists of a set of mutually exclusive homogeneous subgroups.

Assumption 2: (Local Independence) The vector of observed characteristics, $[\mathlcal{y}_1,...,\mathlcal{y}_J]^T$ are related only due to the latent classes.\footnote{Latent classes can be modeled nonparametrically. However, for computational convenience, we model the outcome variables, determined by the latent classes (labels), parametrically using an ordinal logistic regression for each outcome variable.}

Acknowledgment

We thank Larry Manevitz for very constructive comments and Omiros Papaspiliopoulos for very constructive conversations.

\EOD