Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,651 characters · 12 sections · 87 citation commands
Robust Ranking of Happiness Outcomes: A Median Regression Perspective
\setcounter{page}{1}
The study of human happiness has been cited as one of the fastest growing research fields in economics over the last two decades (kahneman_developments_2006, clark_relative_2008, stutzer_recent_2013). By looking at what socioeconomic and other factors predict (or cause) people to report higher or lower scores on a subjective well-being (SWB) scale, researchers have been able to add new insights to what have become standard views in economics. For example, studies of job and life satisfaction have shown that people tend to care far more about relative income rather than absolute income (clark_satisfaction_1996,ferrer-i-carbonell_income_2005), while unemployment is likely to hurt less when there is more of it around (clark_unemployment_2003,powdthavee_are_2007). The use of SWB data has therefore enabled economists to test many of the assumptions in conventional economic models that might have been untestable before the availability of proxy utility data (e.g., di_tella_preferences_2001, stevenson_subjective_2013, gruber_cigarette_2006, boyce_money_2013). It has also led many policy makers to start redefining what it means to be successful as a community and as a nation (kahneman_toward_2004, stiglitz_measurement_2009, de_neve_sdgs_2020).
One objective of happiness research is to understand the determinants of SWB and how they compare across groups. The predominant approach to the SWB data analysis is either through linear regression (OLS) or ordinal parametric methods (e.g. ordered probit or logit); see ferrericarbonell_how_2004 for a comprehensive discussion on the validity of both approaches. Conclusions are then drawn, as is common in applied fields of economics and other social sciences, based on conditional or unconditional mean comparisons using these estimates. Recent studies, however, have suggested that there might be a serious methodological problem associated with these standard estimation approaches.
The issue traces to the fact that SWB is an ordinal measure. For example, consider the 3-point ordinal happiness scale in the General Social Survey (GSS). In the GSS, respondents were asked whether they were \textquotedblleft 1. not too happy\textquotedblright , \textquotedblleft 2. pretty happy\textquotedblright , or \textquotedblleft 3. very happy\textquotedblright . We know that a score of 3 is higher than a score of 2 or 1. But we cannot interpret the extent of which each category is higher than another without further assumptions. Suppose we are interested in comparing a particular statistic between two ordinal variables such as ranking the mean SWB of two different groups. The rank order of any statistic is identified if that order relation remains unchanged when the ordinal variables undergo any increasing transformation. In particular, this implies, the mean ranking of ordinal variables is identified if and only if there holds between them a first order stochastic dominance (FOSD) relation. However, a FOSD relation between two variables does not always exist in which case the mean ranking between them is not identified. In contrast, the mean ranking between cardinal variables is always identified as long as they have finite first moment.
When an ordinal variable is used as the dependent variable in a linear regression, the sign of a slope parameter may be used to indicate the mean ranking across groups. However, the mean ranking is not identified if the sign of such parameter can change by relabelling the ordinal outcomes that preserve the order of the categories. See schroder_revisiting_2017 for a detailed analysis on this issue. A somewhat analogous problem exists in the context of an ordinal regression model. To see this, as we shall do throughout this paper, we interpret observed ordinal outcomes through a threshold crossing latent variable model. In this setting, Bond and Lang (2019, BL hereafter) point out in Section II.B of their paper that any pair of latent variables that follow a continuous two-parameter distribution and have different means, there is a FOSD relation between them if and only if their variances are identical. This result\footnote{ See Theorem 1 in bond_2014 for a general statement of this result. }, in particular, implies the mean rankings in ordered probit and logit models are identified if and only if the respective latent variables are homoskedastic. The following example illustrates its relevance in a simple context of male-female happiness comparison using a probit model.
Example 1: Suppose we observe $Y$ and $D$ that respectively denote reported happiness from a 3-point scale and a female gender dummy so that
If $0<\Pr \left[ D=1\right] <1$, then we can identify $\beta _{1}=E\left[ H|D=1\right] -E\left[ H|D=0\right] $. Let $\beta _{1}\neq 0$ and suppose we want to use the sign of $\beta _{1}$ to determine group comparison. Consider these two situations.
BL demonstrates the practical implication of non-identification by taking on nine of the most well-known findings from the happiness literature. They test and reject the homoskedasticity hypothesis for 8 out of 9 cases\footnote{ bond_sad_2019 did not have individual level data used in ludwig_neighborhood_2012 for comparing happiness between the Control and Experimental groups in the Moving to Opportunity program. They used shares of different happiness levels reported in the paper to back out the mean and variance of latent happiness under normality assumption, from which they can perform a transformation to reverse the result that suggests subjects in the experimental group were happier than those in the control group.}, and then they show conclusions previously drawn from mean comparisons can be reversed by applying some exponential transformation to the latent happiness or SWB variable in all cases.
The discussion above raises an important question about how we should interpret empirical results obtained from ordered probit and logit models. This concern is highly relevant because many applications use standardized ordered probit and logit models , which assume the conditional variance of latent happiness are respectively $1$ and $\frac{\pi ^{2}}{3}$, when the homoskedasticity assumption may not hold in the data. By contrast, the mean ranking in any generalized ordered probit or logit model that presumes heteroskedasticity cannot be identified. Nevertheless, irrespective of whether the mean ranking is identified, many existing empirical findings appear intuitive and have been widely accepted in the literature. This suggests we may be able to learn something about group ranking even if the mean ranking is fundamentally not identified.
Our main goal is to provide a pragmatic view on how to interpret results estimated from ordered logits and probits regardless whether or not the mean rank is identified, as well as to suggest an alternative method to compare ordinal variables generally. Our argument focuses on using the median as a statistic for comparisons instead of the mean. We have the following messages.
We focus on the median, an alternative well-known measure of central tendency, because it is \textquotedblleft equivariant\textquotedblright\ to all increasing transformations. I.e., letting $Med(Z)$ denote the median of a random variable $Z$ and $\tau $ be an increasing function, then $\tau \left( Med(Z)\right) =Med(\tau \left( Z\right) )$. Therefore once a median rank for any pair of latent variables has been established for a particular cardinalization, it cannot be reversed by any monotone transformation.
Example 1 (cont'): Let $Med\left( H|D\right) $ denote the conditional median of $H$ given $D$. By symmetry of the normal distribution, $Med\left( H|D\right) =E\left[ H|D\right] =\beta _{0}+\beta _{1}D$. Therefore $\beta _{1}$ identifies $Med\left( H|D=1\right) -Med\left( H|D=0\right)$. Furthermore, by equivariance, the sign of $\beta _{1}$ is the same as that of $Med\left( \tau \left( H\right) |D=1\right) -Med\left( \tau \left( H\right) |D=0\right) $ for every increasing function $\tau $. Thus the median rank between men and women is identified by the sign of $\beta _{1}$.
The continuation of Example 1 highlights perhaps the most empirically important point of our paper. That is: the median rank of an ordinal variable in an ordered probit or logit model can be identified even when the mean rank is not; furthermore, the median rank is identified by the sign of a model parameter that is routinely estimated in practice. This is due to the fact that the normal and logistic distributions are symmetric and the median identified in logit and probit models coincides with the mean. Specifically, researchers can practically perform group comparisons in the same fashion as previously\footnote{For examples, with Stata, oprobit and ologit can be used to estimate standardized probit and logit models respectively, and oglm can be used to estimate their generalized counterparts where users can specify the form of heteroskedasticity.}, the only change is to interpret rankings in terms of the median instead of the mean. To see why this simple change of stance can be very powerful, let us consider again the empirical illustrations in BL. There, BL show all of their ordered probit estimates that allow for heteroskedasticity deliver qualitatively the same conclusions as in previous studies under homoskedasticity before they apply exponential transformations to reverse the mean rankings. Therefore, under the median interpretation, their estimates would support existing results in the happiness literature rather than dispute them.
Our paper also explores robust estimation of the median. When the analysis is conducted with probit or logit models, it is implicitly assumed that latent happiness follows normal or logistic distribution. Furthermore, in the heteroskedastic case, when the conditional variance of happiness distribution is allowed to vary across respondents with different characteristics, the researcher typically chooses a functional form of the conditional variance. If these distributional or functional form assumptions do not hold, fully parametric models are misspecified and subsequent estimators would be inconsistent. From the econometrics literature, lee_median_1992 has shown it is possible to identify and consistently estimate the median semiparametrically without any parametric distributional assumption as well as functional form of heteroskedasticity. Estimating Lee's estimator in practice, however, can be a very challenging task. His estimator, in particular, is a generalization of the maximum score estimator (MSE), which manski_semiparametric_1985 introduced for estimating a binary choice model and is well-known for difficult implementation, for estimating a model with multiple categories.
Another contribution in this paper is that we provide a new and computationally efficient procedure to compute Lee's estimator. Our computational approach is based on the method of mixed integer optimization, which was used by florios_exact_2008 to implement Manski's MSE. We adapt their procedure designed for a binary choice to a multiple choice setting. Importantly, our procedure can also be used to estimate at any quantile in addition to the median. Since every quantile is equivariant to increasing transformations, quantile ranks are identified under weak conditions. We illustrate in the paper how results estimated at other quantiles can provide additional insights to supplement the conventional mean-median analysis. We apply our estimation procedure and show that a standard happiness equation structure continues to hold under our new approach (blanchflower_well-being_2004).
While there are compelling statistical reasons for pursuing median ranks over mean ranks, especially when the former is identified and the latter is not, a question of normative interest is whether, instead of the mean, a policy maker would prefer to use the median or other quantiles to represent aggregate well-being measures. A recent article by sechel_share_2021 makes an interesting argument that the emphasis on using averages, which implicitly indicates a preference for the conventional mean welfare maximization, may not be a sustainable social goal with finite resources. She provides anecdotal evidence in support of a sufficientarian welfarism view where a policy maker's aim is to have a sufficient level of welfare reached instead of improving welfare for everyone. More specifically, she suggests a planner may want to analyze a measure constructed from averaging \textquotedblleft headcounts\textquotedblright\ of reported well-being scores that is no less than a targeted threshold. I.e., the threshold corresponds to the $\alpha$-quantile of the reported well-being score where $\alpha $ is the average headcount in that sample. The ordinal regression framework we study in this paper can be particularly useful for a sufficientarian planner as it is able to perform group comparisons within a population at any targeted quantile level. In addition, from the decision theory literature, an economic agent whose goal is to maximize the quantile of their utility distribution has a well-established foundation. We refer readers to manski1988 and, for more recent developments, to rostek2010 and de2019 for decision theoretical justifications.
We emphasize that our work is not a critique of BL. Their theoretical point that identification of the mean ranking of ordinal data in familiar parametric models is possible only when homoskedasticity holds is insightful and cannot be disputed. We hope, however, to prevent the potentially extreme interpretation of their results that nothing about ordinal ranking can be learnt from popular parametric models and the value of prior works that used probit or logit models rests on the knife-edge condition whether the model is homoskedastic or not.
BL also point out another challenge for ordinal data analysis that is relevant for median ranking. They question the often assumed notion in empirical studies of a common reporting function that puts the latent variables from different groups on the same cardinal scale. In a threshold crossing model, this corresponds to the possibility that threshold values are heterogeneous across groups. We also maintain the common threshold assumption in the present paper. One may relax this assumption by using a parametric compound hierarchical ordered response model (see KinMurSal04). To the best of our knowledge, extensions to semiparametric ordinal median regression models with heterogeneous thresholds have not yet been developed and would be an important topic for further research.
The remainder of this paper proceeds as follows. Section 2 gives an account on how statistical analysis for discrete ordinal data have been developed in economics and other disciplines, and on some recent developments on estimating the median. Section 3 introduces an ordered response model and formalizes our argument to identify the median and other quantile ranks. Section 4 compares parametric and semiparametric estimation, and introduces the mixed integer optimization approach to the computation of the median and quantile estimators in a semiparametric ordered response model. Section 5 revisits estimation of a happiness equation using GSS data. Section 6 concludes. The Appendix provides the details of the mixed integer optimization based implementation algorithm in the context of semiparametric median and quantile estimation problems.
We consider discrete ordinal outcome that represents an individual's ordered categorical response in the data. The defining property of an ordinal variable is that there is a rank order over values it can take but the distances between these values are arbitrary and carry no information. A discrete ordinal variable can therefore, in constrast to cardinal variables, be put on a scale like $ \left\{ 1,2,..,J\right\} $ without any loss of generality. Such data measurements are common in social and biomedical sciences. Examples include individual happiness (unhappy, neither happy nor unhappy, happy), severity of injury in the accident (fatal injury, incapacitating injury, non-incapacitating, possible injury, and non-injury), and lethality of an insecticide (unaffected, slightly affected, morbid, dead insects) among many others.
Ordinal data analysis in a regression framework is widely acknowledged to have been co-founded by two independent sources. One originates from the contribution of mckelvey_statistical_1975, who developed the well-known ordered probit model to study Congressional voting on the 1965 Medicare Bill. The other is due to mccullagh_regression_1980, who focused on modelling proportional odds and proportional hazards that become prominent in the biomedical fields. Huge literature on ordinal data analysis has since grown from these influential works. We refer interested readers to greene_hensher_2010 and reference therein for the developments in social science, agresti_modelling_1999 for the medical science, and ananth_regression_1997 for epidemiology.
Researchers from different fields take different approaches to analyzing ordinal data. Applied researchers in biomedical fields pay a great deal of attention to choosing an appropriate model for their data (goodness of fit) but place less importance on the interpretation of individual parameters (e.g. coefficients in a generalized linear model). On the other hand, researchers in social sciences often focus on the model parameters. Economists, for example, typically work with linear regression or ordered probit and logit models that are very convenient for interpreting parameters, especially in a setting with many covariates. One research area in economics that has utilized ordered discrete response models the most is perhaps the economics of well-being. Given the ordinal nature of SWB data, many early and classic studies in this field have exclusively used ordered probit and logit models to analyze different predictors of human happiness, including unemployment (clark_unhappiness_1994), political institutions (frey_happiness_2000), income inequality (alesina_inequality_2004), and relative income (ferrer-i-carbonell_income_2005). More recently, many researchers have also been using linear regression models to conduct their analysis. For example, OLS model has been used to estimate the relationship between SWB and macroeconomic factors such as inflation and unemployment rate (di_tella_preferences_2001), comparison income (luttmer_neighbors_2005), and fertility (kohler_population_2005). One particular advantage for using linear models is the ease in incorporating and accommodating fixed effects when panel data are used.
The linear regression method treats ordinal variables as if they were cardinal. Several studies have shown linear regression and ordered probit/logit models could deliver similar qualitative results empirically. For example, in a study of vehicle driver injury severity, gebers_exploratory_1998 compared both OLS and ordered logit estimation results and found that estimated coefficients of both models were of the same sign and generally agreed in magnitude and statistical significance. In an empirical analysis of the effect of economics sanctions, major_timing_2012 compared the OLS and ordered probit estimates and found that the results were similar across these two estimating models. See also ferrericarbonell_how_2004, who found both the OLS and ordered probit models produced comparable results in their study on the sources of individual well-being.
One practical issue with the linear regression approach is the dependency on the reported scale. Particularly, it is possible for OLS estimates to change signs if the reported data are monotonically transformed. See, e.g. schroder_revisiting_2017 for a sufficient condition of this; also see kaiser_how_nodate for a sufficient condition for which the signs of OLS estimates cannot be reversed. Some recent works have suggested complementing OLS estimates obtained from the reported scale with those that undergo some transformations as a sensitivity analysis. E.g., bloem_how_nodate proposes a class of one-parameter monotone functions and suggesting to estimate a set of OLS estimates indexed by that parameter; bloem_analysis_2021 suggest dichotomizing the dependent variable, e.g. into high and low happiness groups using the median of reported happiness to form the split, to mitigate the dependency of the reported scale.
For latent variable threshold crossing ordinal regression models, the regression coefficients are by construction invariant to any order preserving relabeling of the ordinal outcome. These parameters naturally have an interpretation of the conditional mean difference of latent variables across different groups. Their signs can thus be used to identify the group mean ranking whenever the ranking order is invariant for all increasing transformations of the latent variable. This requirement amounts to the condition of there being a FOSD relation between latent variable distributions of the groups. In a fully parametric model such as ordered probit and logit, given a specific form of heteroskedasticity and assuming the means across groups differ, a hypothesis test of homoskedasticity would determine if there is a FOSD relation. It is also possible to test the null against a nonparametric alternative in this setting because the conditional variance function can be nonparametrically identified within a probit or logit framework (e.g. see oparina_analyzing_2021 ). In these cases, rejecting homoskedasticity means there is no FOSD relation. This knife-edge condition makes identification of the mean rank order in probit and logit models very fragile as BL have illustrated. We refer readers to their paper for further discussions. More generally one can test FOSD directly in a semiparametric or nonparametric setting that does not assume distributional assumption of the latent variable. See e.g. carneiro_estimating_2003, cunha_identification_2007, lewbel_constructing_1997, lewbel_semiparametric_2000, lewbel_simple_2007, honore_semiparametric_2002 and KaplanZhuo2021b.
In this paper we propose that a natural alternative for ranking ordinal outcomes is to focus on the median instead of the mean. For commonly used ordered response models such as the ordered probit and logit, the median and mean are identical. Maximum likelihood estimation of the median in these models can therefore be performed as readily as the mean using standard statistical softwares. The median can also be estimated semiparametrically without distributional assumption. The seminal work of manski_maximum_1975, manski_semiparametric_1985 develops the maximum score estimation in this setting for a binary choice model and lee_median_1992 extends this to multiple ordered choice data. While theoretically appealing, performing maximum score estimation for the discrete choice model is computationally difficult. More specifically, the ordered response median regression estimator proposed by lee_median_1992 is a solution to a non-smooth and non-covex least absolute deviation (LAD) optimization problem. The corresponding LAD objective function, which is akin to Manski's maximum score objective function, is non-convex -- it is piecewise constant with numerous local solutions. The computational challenges for solving the maximum score estimation problem is well noted in the econometrics literature (e.g. see manski_operational_1986, pinkse_computation_1993, skouras_algorithm_2003).
Recently florios_exact_2008 propose a mixed integer optimization (MIO) based approach for the computation of maximum score estimators. In particular, they show that Manski's binary choice maximum score estimation problem can be equivalently reformulated as a mixed integer linear programming problem (MILP). chen_best_2018 provide an alternative MILP formulation that complements the approach of florios_exact_2008 for solving the maximum score estimation problem. These reformulations enable exact computation of maximum score estimators through modern efficient MIO solvers. Well-known numerical solvers such as CPLEX and Gurobi can be used to effectively solve the MIO problems. We refer readers to bertsimas_optimization_2005 and conforti_integer_2014 for recent and comprehensive texts on the MIO methodology and applications.
The estimators of manski_semiparametric_1985 and lee_median_1992 have been shown to be consistent under very weak conditions. On the other hand, maximum score type estimators converge at a cube-root rate and have non-standard asymptotic distributions. See kim_cube_1990 and seo_local_2018. When there are continuous covariates, one can use the smoothed maximum score (SMS) estimator proposed by horowitz_smoothed_1992 that employs a smooth approximation of the original maximum score objective function. The SMS estimator is asymptotically normally distributed and can have a faster rate of convergence than the unsmoothed maximum score estimator. However, for implementation of the SMS estimator, users have to choose tuning parameters that might induce a smoothing bias which could be difficult to correct (kotlyarova_robust_2009). We refer the reader to horowitz_semiparametric_2009 for a detailed review on the theoretical aspects of the maximum score and SMS estimators.
We now present an empirical model in the context of happiness application. Suppose we are interested in using SWB data taken from two groups, say $A$ and $B$, to draw conclusions on whether people in group $A$ are happier than those in group $B$. Examples of group identities include gender, martial status, country, time etc. Previous analyses have been focusing on the mean as a statistic to compare happiness across groups. We will focus on the median and discuss analogous results for other quantiles.
Let observable data for an individual be $\left( Y,X,D\right) $ where $Y$ \ denotes a reported happiness scale taking values from $\mathcal{Y}=\left\{ 1,\ldots ,J\right\} $ for $J\geq 3$, $X$ is a vector of covariates other than the group identity of interest, and the latter is denoted by a dummy variable $D$ that takes value $1$ for group $A$ and $0$ for group $B$. We assume $Y$ is derived from a threshold crossing model based on latent happiness index, $H$, such that
where
for some strictly increasing real thresholds $\left\{ \gamma _{j}\right\} _{j=1}^{J-1}$ with $\gamma _{0}=-\infty ,\gamma _{J}=+\infty $, $\beta _{0}$ is a vector of parameters associated with $X$, $\beta _{1}$ is a scalar parameter associated with $D$, $U$\ is an unobserved scalar accounting for other factors, and the symbol $^{\top }$\ denotes the transpose.
The model described in this section is frequently used in practice with an additional assumption on the distribution of $U$. For instance, a standardized probit model assumes $U|X,D \sim N(0,1)$ and a generalized probit model assumes $U|X,D \sim N(0,\sigma^{2}(X,D))$ where $\sigma^{2}(X,D)$ is specified up to some unknown parameters. In what follows, we shall focus on the use of the sign of $\beta _{1}$\ for ranking the median across groups. Particularly, the definition of the median rank can be stated without specifying the distribution of $U|X,D$. However, as we shall see, some assumption on the distribution of $U|X,D$ is required to identify the median rank.
The parameter of interest for median comparison across groups with observed characteristics $X$\ is:
I.e. $\lambda \left( X\right) $\ is the difference between median levels of happiness for individuals from groups $A$ and $B$ with the same characteristics $X$. The median is equivariant to increasing transformation, so the signs for \
will be the same for all increasing function $\tau $. The sign of $\lambda \left( X\right) $\ can thus be used to identify the median rank between $ H|X,D=1$\ and $H|X,D=0$. We state this as a proposition.
Proposition 1.
The median rank of:
Proposition 1 shows the sign of $\lambda \left( X\right) $\ can be used to establish the median rank regardless whether $H$\ is homoskedastic or not, and its validity does not require any knowledge of the distributional assumption.
We can impose some assumptions on the distribution of $U|X,D$ for the identification of $\lambda \left( X\right) $. The simplest approach is to assume $U|X,D$ follows either a normal or logistic distribution. Example 2 below illustrates this using a more general probit model than Example 1.
Example 2: Consider model ((ref)) where $U|X,D\sim N\left( 0,\sigma ^{2}\left( X,D\right) \right) $ and $\gamma _{1}=0$ and $ \gamma _{2}=1$. Then $H|X,D\sim N\left( X^{\top }\beta _{0}+\beta _{1}D,\sigma ^{2}\left( X,D\right) \right) $ and $\lambda \left( X\right) =\beta _{1}$. The sign of $\beta _{1}$\ determines the median happiness rank between groups $A$ and $B$. This holds irrespective of the form of $\sigma ^{2}\left( X,D\right)$.
Estimation of $\lambda \left( X\right) $\ for probit and logit models can be performed using standard statistical softwares. We emphasize, however, such parametric assumptions are not necessary for identification and estimation. We will describe a semiparametric estimator in Section 4 where the only thing we assume about $U$ is a median restriction.
The previous discussion focuses on median rank. The median is a special case of an $\alpha$-quantile with $\alpha =0.5$. More specifically, let us denote the conditional $\alpha $-quantile of a continuous random variable $Z$\ given $W$ by $Q_{Z}\left( \alpha |W\right) $\ so that it satisfies,
In happiness applications it is commonly assumed that $U$, and hence $H$, is a continuous random variable. In this case $Q_{H}\left( \alpha |X,D\right) $\ is equal to the inverse of the CDF of $H$ conditional on $\left( X,D\right) $\ evaluated at the quantile level $ \alpha $.
The parameter for performing $\alpha$-quantile comparison across groups with observed characteristics $X$, which is the counterpart to ((ref)), is:
I.e. $\lambda \left( \alpha |X\right) $\ is the difference between $\alpha$- quantile levels of happiness for individuals from groups $A$ and $B$ with the same characteristics $X$. Like the median, all quantiles are equivariant to increasing transformations\footnote{ That is, for any increasing function $\tau $, $Q_{\tau \left( H\right) }\left( \alpha |X,D\right) =\tau \left( Q_{H}\left( \alpha |X,D\right) \right) $.} property, so we know the signs of
will be invariant for every increasing function $\tau $. Thus, the sign of $\lambda \left( \alpha |X;\tau \right) $\ can be used to identify the $\alpha$-quantile rank between $H|X,D=1$\ and $H|X,D=0$.
Proposition 2. For any $\alpha\in(0,1)$, the $\alpha$-quantile rank of:
Similar to Proposition 1, quantile ranks is determined by the sign of $ \lambda \left( \alpha |X\right) $\ without homoskedasticity or distributional assumptions.
In practice, we can use a linear quantile model for latent happiness $H $\ to identify $\lambda \left( \alpha |X\right) $. For example, suppose we assume
Then $\lambda \left( \alpha |X\right) =\beta _{1}\left( \alpha \right)$. Estimation of $\beta _{1}\left( \alpha \right)$ and $\beta _{0}\left( \alpha \right)$ is not straightforward even under parametric assumptions as these coefficients depend on the specification of heteroskedasticity of the distribution of $H$ conditional on $X$ and $D$. If we were able to observe $H$ and assume it has a continuous distribution, then we could estimate $\left( \beta _{0}\left( \alpha \right),\beta _{1}\left( \alpha \right) \right)$ using the quantile regression method of koenker_simple_1978. However, such approach is not applicable in our case because the outcome $Y$ is discrete. In the next section, building on the ordinal median regressor estimator of lee_median_1992, we will present an estimation approach that allows for estimating the quantile regression coefficients of ((ref)) in the context of discrete ordered response model ((ref)).
We begin this section by outlining the main challenges relating to the estimation of semiparametric median and quantile regressions for ordinal outcomes. We propose that a mixed integer optimization procedure can be useful in this setting and propose an estimation procedure for them. The details on mixed integer optimization can be found in the appendix. As in the previous section, we will focus on the median case and discuss the cases for other quantiles afterwards.
We begin with the following model of semiparametric median regression for ordinal outcomes:
where the covariate vector $X\in\mathbb{R}^{p+1}$ subsumes group dummy variables. Equation ((ref)) amounts to rewriting ((ref)) using indicator functions. Let $\gamma _{0}=-\infty ,\gamma _{J}=+\infty $\ and $ \gamma :=\{\gamma _{j}\}_{j=1}^{J-1}$ denotes the vector of other threshold parameters that are strictly increasing. The parameters of the model are $\left( \theta ,\gamma \right) $.
For presenting our estimation procedure, it will be convenient to write $ X=(X_{1},\widetilde{X})$, where $X_{1}\in \mathbb{R}$ is the first element of $X$\ and $\widetilde{X}\in \mathbb{R}^{p}$ is a vector containing the remaining components of $X$. Correspondingly, we let $\theta =(\theta _{1},\beta )$ so that $ X^{\top }\theta =\theta _{1}X_{1}+\widetilde{X}^{\top }\beta $. Since parameters of the ordinal choice model ((ref))-((ref)) are not identified without location and scale normalizations, we assume $X$\ does not contain a constant (location normalization) and $\left\vert \theta _{1}\right\vert =1$ (scale normalization). There are other ways to normalize. For example, another convenient choice is to set $\gamma _{1}=0$ and $\gamma _{2}=1$ as done in Example 1. We refer the reader to greene_hensher_2010 for further discussions on normalization in ordered choice models.
To estimate the median semiparametrically we assume that $Med(U|X)=0$ for all $X$ .\ This allows for nonparametric specification of the distribution of latent unobservable $U$, which admits general unknown form of heteroskedasticity conditional on the covariates. In particular, it can be shown the zero median condition implies that
which subsequently yields,
We want to highlight that the previous two equations do not depend on the variance or other distributional knowledge about $U|X$. This feature is attractive because we can then use ((ref)) to estimate the median without any risk of misspecifying aspects about the distribution of $U|X$. It is in this sense that the semiparametric estimator is robust.
Let $\Theta \subset \mathbb{R}^{p+J-1}$\ denote the parameter space containing $\left( \beta ,\gamma \right) $. We use $\left( b,c\right) $ to denote a generic point of $\Theta $. Then, given a random sample $\left( Y_{i},X_{i}\right) _{i=1}^{n}$, letting $ c_{0}=-\infty $ and $c_{J}=\infty $, consistent estimator $(\widehat{\theta } _{1},\widehat{\beta },\widehat{\gamma })$ of $(\theta _{1},\beta ,\gamma )$ can be obtained as a solution to the following least absolute deviation (LAD) estimation problem\footnote{ We refer the reader to lee_median_1992 for a comprehensive discussion of this LAD estimator.}:
The LAD estimator that solves the minimization problem in ((ref)) is a generalization of the maximum score estimation approach of manski_semiparametric_1985 for the binary response model to that for the discrete ordered response case. As in the problem of maximum score estimation, the objective function in ((ref)) is piecewise constant with numerous local minima. Thus, although this LAD estimator is theoretically well-defined, finding a global solution in practice is known to be a notoriously difficult task.
We propose an algorithm which enables exact computation of a global solution to the LAD problem ((ref)) using mixed integer linear programming. The idea is to turn the difficult optimization problem with a non-convex objective function over a convex domain ($\Theta $) to an optimization problem with a linear objective function over a domain that embodies integrality restrictions, which means the parameters are passed on as a set of linear inequality constraints that need to be satisfied. Specifically, we show that the integer values to be considered are finite, which can then be efficiently searched using modern mixed integer optimization (MIO) methods. This is the insight that florios_exact_2008 have exploited to compute Manski's maximum score estimator. The details of the mixed integer optimization procedures that can be used to estimate the semiparametric median as well as other quantiles are given in the Appendix.
The theoretical advantages of the semiparametric estimator is that the results will be robust to model misspecification, particularly on the form of heteroskedasticity and distribution of $U$ conditional on $X$. In particular, if the functional form of the variance of $U$ given $X$ is misspecified in probit or logit models, or if the conditional distribution of $U$ given $X$ is in fact not normal or logistic respectively, these parametric estimators are inconsistent.\ On the other hand, if the parametric assumptions are correct then the corresponding parametric estimator will be asymptotically more efficient than the semiparametric estimator in terms of having smaller asymptotic variance.
There are other considerations to bear in mind when using the semiparametric estimator. First, the researcher has to choose the variable $X_1$ that is associated with $\theta_1$, which we normalize to be -1 or 1, carefully. The chosen $X_1$ is assumed a priori to have non-zero effect on latent happiness though the sign of its effect, positive or negative, can be estimated. Ideally $X_1$ should also have a rich support as it facilitates parameter identification.\footnote{Standard sufficient conditions for identification in maximum score type estimation problems assume one of the covariates has full support on $\mathbb{R}$ conditional on other explanatory variables. The full support condition is, however, not necessary. It can be reduced to bounded or even finite support as explained in manski_identification_1988 and horowitz_semiparametric_2009.} On the practical front, estimating the semiparametric median regression is also computationally more involved than estimating ordered probit or logit models. In particular, despite using an efficient solver, it generally can take a long time to solve an MIO problem and the computational time increases with the sample size. This latter point can be seen by inspecting equation ((ref)) in the appendix where it is apparent that the number of the optimizing integer variables grows with the sample size. In contrast, it is much easier to compute the maximum likelihood estimators for the logit and probit models. Inference on parametric models is also relatively straightforward. For instance, the Stata commands mentioned earlier automatically provide standard errors and confidence intervals for the median estimators. Inference based on the semiparametric estimator here can be performed by resampling, which further exacerbates the computational cost.
By design, it is more difficult to identify and estimate the model without additional parametric structures. Yet the semiparametric approach can deliver more robust results against model misspecifications. Our view is that parametric and semiparametric estimators should be complementary rather than used as a substitute for each other in conducting empirical studies.
We now extend the ordinal median regression approach of Section 4.1 to the setting of estimation at other quantiles. Let us consider a linear quantile model for the latent variable $H$ in ((ref)). Specifically, using the notation introduced in Section 3.3, we assume that for $\alpha \in (0,1) $, $Q_{H}\left( \alpha |X\right)=X^{\top }\theta $. This assumption implies that there is a residual $U$ such that ((ref)) holds with $Q_{U}\left( \alpha |X\right)=0$. In this case it can be shown that,
which encompasses ((ref)) as a special case where $\alpha=0.5$. Using the same location and scale normalizations as in Section 4.1, we aim to estimate $\theta=(\theta_{1},\beta)$ where we can write $X^{\top }\theta=\theta_{1}X_{1}+\tilde{X}^{\top}\beta$ with $\vert\theta_{1}\vert=1$. In particular, analogously to ((ref)), we can estimate $(\theta_{1},\beta,\gamma)$ by solving the following minimization problem:
where $\rho_{\alpha}$ is the check function used in quantile regression (koenker_simple_1978) defined as $\rho_{\alpha}(u):=u\left(\alpha-1\left\{u\leq0\right\}\right)$ for any $u\in\mathbb{R}$. Indeed, the median estimator that minimize ((ref)) coincides with that of ((ref)) when $\alpha = 0.5$. It is not a simple task to solve ((ref)) for the same reason that makes solving ((ref)) difficult. Conveniently, the computational strategy we proposed based on mixed integer linear programming (MILP) to estimate the median can readily be adapted to estimate other quantiles. Thus, the semiparametric ordinal quantile regression estimator can be effectively computed through the method of mixed integer optimization for any quantile level.
To demonstrate the practical usefulness of our approaches, we apply them to estimate a standard happiness equation for studying correlation between individual well-being and demographic and socioeconomic characteristics (see e.g. blanchflower_well-being_2004). Previous studies have focused on understanding how certain socio-demographic characteristics, e.g. unemployment and martial status, are related to the average well-being. We are interested in learning how these factors affect the median and other quantiles of happiness distribution. More importantly, we are interested in whether the structure of the happiness equation remains qualitatively similar when the median instead of the mean is the outcome of interest. In particular, we consider ordered probit and logit models, where the underlying happiness distribution is assumed to be homoskedastic as well as heteroskedastic. We also estimate the semiparametric quartiles of the happiness distribution, viz. the 25th, 50th, and 75th percentiles, to investigate heterogeneous effects at different parts of conditional happiness distribution.
We use bi-annual data taken from the US General Social Survey (GSS) between the years 2000 and 2018 inclusive. It consists of 19,275 observations. The GSS happiness variable comes from respondents being asked whether they are “1. not too happy”, “2. pretty happy”, or “3. very happy”. The explanatory variables for the happiness equations are income, age, age-squared, sex, level of education, marital status, race, employment status, and time fixed effects. We use the logarithm of equivalence scale adjusted household income. A set of education dummies includes indicators for those who left high school, completed a bachelor and had a graduate degree, with high school/junior college being the omitted category. The marital status dummies include married, divorced, widowed, with the omitted category being never married. Race dummies include black and other race, white is used as the omitted category. Finally, employment status dummies include unemployed, non-labour status, and the omitted category is for those that work fulltime or parttime.
We estimate our parametric models under normal and logistic distributions using Stata. For the standardized probit and logit models, which assume homoskedasticity, we use the commands oprobit and ologit respectively. For the generalized ordered probit and logit models that allow for parametric forms of heteroskedasticity we use the oglm command.\footnote{Due to the many explanatory variables involved we specify the heteroskedastic function parametrically as the nonparametric approach in oparina_analyzing_2021 does not work well in this scenario. Specifically, the nonparametric approach requires nonparametric estimation of the conditional probabilities for different happiness levels. We would have very few observations in some cells leading to imprecise estimates because we have many discrete explanatory variables to condition on.} Specifically, the vector of covariates in the latter models does not contain a constant, as a location normalization, and the conditional variance function is an exponential function that contains the same vector of covariates to incorporate scale normalization.\footnote{Using the notation in ((ref)), we use oglm to estimate the generalized ordered probit under the assumption that $H=X^{\top }\theta +U$ where $U=\sigma ( X^{\top }\delta ) \varepsilon $ with $\sigma^{2}(t)=e^t$ and $ \varepsilon \sim N\left( 0,1\right) $ that is independent of $X$. We estimate the generalized ordered logit analogously where $\varepsilon $ would then follow a standard logistic distribution. Further details on oglm can be found in williams_fitting_2010.}
To estimate the semiparametric median and quartiles, we impose the normalizations described in Section 4.1. We do not have the constant term in the regression and the magnitude of the coefficient on income is set to one. We note that our income variable has a rich support\footnote{The GSS collected household income data in 23 to 26 bands (depending on the year). We computed their mid-points and debased them accordingly. Subsequently, the income variable in our sample can take 246 possible values. Since we used equivalence scale adjusted household income, it takes 1425 possible values altogether.}. Our semiparametric estimates are obtained using the MILP procedure described in the appendix. Specifically, we used the MATLAB implementation of the Gurobi Optimizer to solve the MILP problems.\footnote{The codes and data used for obtaining the results in our paper are publicly available from the following repository \url{https://data.mendeley.com/datasets/dj5mpwypsc}.}
Table 1 presents the ordered logit and probit estimation results. The top half of the table shows those estimates that are associated with the median, which is the same as the mean in these instances, and the bottom half are concerned with the variance estimates. We only report estimates for the socioeconomic factors. Time fixed effects, which are typically not the focus when analyzing the happiness equation, are omitted\footnote{For a particular model, the time effects are statistically significant at $1\%$ for some years but not significant at $10\%$ for some other years; there is no clear pattern for statistical significance across years. However, the signs of the effects and significance results for each year are the same in all models.} for the sake of space.
The median/mean results for the socioeconomic effects are in line with previous findings often reported in the well-being literature. For example, respondents with higher income report higher levels of well-being (e.g., easterlin_does_1974, ferrer-i-carbonell_income_2005); well-being follows a U-shape throughout one’s life (e.g., blanchflower_is_2008) where the convexity in age has been attributed to the midlife nadir (e.g. see cheng_longitudinal_2017); women are happier than men (stevenson_paradox_2009); university graduates are happier than those who do not have a degree (blanchflower_well-being_2004); being married is positively related with well-being, while being black or of other race, unemployed or not in labour force has negative correlation (e.g., easterlin_explaining_2003). We find almost no qualitative difference in the median estimates from the homoskedastic model and the heteroskedastic one despite there being evidence for heteroskedasticity. On the latter, we find income, low educational attainments, race, and employment status have significant effects on the variance of conditional happiness distribution of respective groups. Formal tests, based on the Wald statistic, also reject the null hypothesis of homoskedasticity\footnote{The null hypothesis corresponds to the setting where all the variance coefficients in the conditional variance function are zero.} with p-values below $0.01$. The probit and logit results described above hold almost identically except that the effect of the other race status is significant at $10\%$ for the heteroskedastic logit model but not significant at that level in all other cases. The fact that we find evidence of income affecting happiness, even in parametric models where the results may be biased due to misspecification, is a useful precursor for semiparametric estimation. This is because in the semiparametric model we assume the income coefficient is non-zero when we impose the scale normalization through this parameter.
Table 2 tabulates the median estimates of parameters and their confidence intervals from the generalized probit and logit models and the semiparametric counterpart. We note that the parametric estimates in Table 2 differ from those in Table 1 because we divide the parametric estimates by their respective estimated coefficient of income to facilitate comparisons with the semiparametric estimates. The reported semiparametric estimation results correspond to the set of estimates obtained when the sign of the income coefficient is positive, which is the sign that minimizes the LAD problem, so it is consistent with the positive income effects seen in the parametric model. The $95\%$ confidence intervals for the parametric and semiparametric estimators are computed using normal approximation and $m$ out of $n$ bootstrap respectively. We note that for Tables 2 and 3, we only signify statistical significance at $5\%$ level as the only confidence intervals constructed for our semiparametric estimators are the $95\%$ confidence intervals.
The most striking aspect from Table 2 is that all parametric and semiparametric estimates have the same sign for all socioeconomic variables. The insignificant/significant effects at the 5% level also coincide in most of these cases. The semiparametric model suggests more pronounced effects, in terms of statistical significance, for being widowed or other race than the parametric models, and finds the unemployed status to be negative but insignificant. These differences may be due to the symmetric distributional assumption imposed by the parametric models.
To better understand how well-being is related with individual characteristics across happiness distributions, Table 3 provides the estimates across the first (25th percentile), second (50th percentile) and third (75th percentile) quartiles. Most socioeconomic variables affect happiness across the three quartiles in similar ways if we look at the signs of the estimates. There are, however, some differences across quartiles in terms of statistical significance for some factors, but none in which the signs differ for different quartiles that are both statistically significant. Two interesting differences, using the median as the benchmark, are: the negative effect from being unemployed, which is not significant for the median, is significant for the lower and upper quartiles; and, the significant negative impact for being black or other race also stand in contrast with the non-significant effects for the lower quartile (negative) and upper quartile (positive). Other than these we find the factors that drive happiness at the 50th and 75th percentiles to be the same. In contrast, a particular pattern we observe is that, at the 25th percentile of happiness distribution, several drivers (female, bachelor, widowed, \textit{black}, \textit{other race}) of happiness that are significant at higher percentiles become insignificant.
Taken together, Tables 1 - 3 show our results largely align with previous findings in the well-being literature, which will be reassuring for many researchers working in this field. This is particularly clear for median regression analysis, and also true to a large extent at other quartiles of the happiness distribution. We summarize our results as follows.
A group ranking of ordinal outcomes is identified only when the ranking order is invariant across all increasing transformations on the ordinal variables. For the mean ranking, this invariance is equivalent to there being a first order stochastic dominance (FOSD) relation between the variables across groups. The usefulness of the probit and logit based mean ranking for discrete ordinal outcomes has in particular been put under question, as illustrated by bond_sad_2019, because FOSD in this setting requires the model to be homoskedastic, otherwise the mean rank is not identified.
In this paper we propose focusing on the median as a pragmatic alternative to the mean. Firstly, the median rank of ordinal outcomes can be identified even when the mean rank is not. Secondly, probit and logit based median ranks can be identified by the conditional means of latent variables of these parametric models and can hence be easily estimated using standard statistical softwares. Thirdly, we also propose a mixed integer optimization procedure to perform median regression in a semiparametric ordered response model that can accommodate unknown distribution of the latent unobservable. Furthermore, our mixed integer optimization procedure can be used to estimate other quantiles of the distribution in addition to the median. Quartile ranks are also identified under weak conditions similar to the median. Estimates of quantile regressions can be useful for researchers to understand effects of happiness driving factors across different parts of the happiness distribution.
In our empirical study, we revisit the happiness equation for the US using the GSS data. We find that median estimates of the ordered probit and logit results are qualitatively very similar to the semiparametric ones. These results, which are in line with the conventions in the happiness literature, suggest that structures, e.g. symmetry, imposed by familiar parametric models have limited role in determining how different socioeconomic factors affect median of the happiness distribution. While many qualitative results of the estimated semiparametric median are also found in the estimation of the lower and upper quartiles of the happiness distribution, there are some notable differences. Compared to the higher quartiles, the happiness distribution at the lower quartile seem to be affected by fewer factors. For example, being widowed, black or other race seem to matter less. But having dropped out of high school or being unemployed seem to have particularly high negative association with happiness in the lower quartile. Thus, a policy maker who wants to improve welfare of people at lower quantiles may want to focus on these factors for further investigations.
The model we consider in the paper is applicable in the context with single or repeated cross-sectional data. There are also many well-being applications that use panel data. The median in these models can be well estimated with a suitable procedure. For example, consider a panel probit or logit model of happiness with fixed effects that satisfies: $E \left[ H_{it}|X_{it},\mu _{i} \right] = X_{it}^{\top }\beta+\mu _{i},$ where the indices $i$ and $t$ denote individual and time respectively, and $\mu _{i}$\ denotes the unobserved individual fixed effect. By symmetry, the conditional median of $H_{it}$ given $\left( X_{it} , \mu_i \right)$ is also $X_{it}^{\top }\beta + \mu _{i}$. However, it is well-known that $\left( \beta,\mu _{i} \right)$ cannot be consistently estimated by maximum likelihood (ML) with fixed time periods $T$ due to the incidental parameter problem. Applications often focus on $ \beta$ and the finite sample bias of ML estimators can be substantial when $T$ is small. A popular way to deal with this in the econometrics literature is to conduct a bias-correction procedure. E.g., hahn_jackknife_2004 and bester_penalty_2012 provide methods to do this for general non-linear parametric panel data models that include ordered probit and logit. Another approach, which is based on dichotomization (chamberlain_analysis_1980), specific for an ordered logit model can yield a consistent estimator of $\beta$ with finite $T$. Dichotomization-based estimators have been quite popular in well-being applications. We refer the reader to baetschmann_consistent_2015 for an examples of such estimators that include a consistent and efficient version. All these estimation methods are viable options for estimating the conditional median of probit and logit models with fixed effects. Semiparametric estimation of median and quantile models with fixed effects for ordinal outcomes can also be performed using maximum score type estimation methods that attempt to estimate $\mu _{i}$ in a large $T$ framework. However, bias correction in these settings is more challenging, which is an interesting topic for further research.