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.
99,883 characters · 10 sections · 29 citation commands
Symmetric generalized Heckman models
{ { Keywords.} {Generalized Heckman models $\cdot$ Symmetric distributions $\cdot$ Variable dispersion $\cdot$ Variable correlation.}}
It is common in the areas of economics, statistics, sociology, among others, that in the sampling process there is a relationship between a variable of interest and a latent variable, in which the former is observable only in a subset of the population under study. This problem is called sample selection bias and was studied by heckman76. The author proposed a sample selection model by joint modeling the variable of interest and the latent variable. The classical Heckman sample selection (classical Heckman-normal model) model received several criticisms, due to the need to assume bivariate normality and the difficulty in estimating the parameters using the maximum likelihood (ML) method, which led to the introduction of an alternative estimation method known as the two-step method; see heckman79. Some studies on Heckman models have been done by Nelson1984, Paarsch1984, Manning1987, Stolzenberg1990 and Yu1996. These works suggested that the Heckman sample selection model can reduce or eliminate selection bias when the assumptions hold, but deviation from normality assumption may distort the results.
The normality assumption of the classical Heckman-normal model heckman76 has been relaxed by more flexible models such as the Student-$t$ distribution Genton2012,Ding2014,lachosetal21, the skew-normal distribution Ogundimu_2016 and the Birnbaum-Saunders distribution bastosbarretosouza:21. Moreover, the classical Heckman-normal model assumes that the dispersion and correlation (sample selection bias parameter) are constant, which may not be adequate. In this context, the present work aims to propose generalized Heckman sample selection models based on symmetric distributions fkn:90. In the proposed model, covariates are added to the dispersion and correlation parameters, then we have covariates explaining possible heteroscedasticity and sample selection bias, respectively. Our proposed methodology can be seen as a generalization of the generalized Heckman-normal model with varying sample selection bias and dispersion parameters proposed by Bastos2021 and the Heckman-Student-$t$ model proposed by Genton2012. We demonstrate that the proposed symmetric generalized Heckman model outperforms the generalized Heckman-normal and Heckman-Student-$t$ models in terms of model fitting, making it a practical and useful model for modelling data with sample selection bias.
The rest of this work proceeds as follows. In Section (ref), we briefly describe the bivariate symmetric distributions. We then introduce the symmetric generalized Heckman models. In this section, we also describe the maximum likelihood (ML) estimation of the model parameters. In Section (ref), we derive the generalized Heckman-Student-$t$ model, which is a special case of the symmetric generalized Heckman models. In Section (ref), we carry out a Monte Carlo simulation study for evaluating the performance of the estimators. In Section (ref), we apply the generalized Heckman-Student-$t$ to two real data sets to demonstrate the usefulness of the proposed model, and finally in Section (ref), we provide some concluding remarks.
Let $\boldsymbol{Y}=(Y_1,Y_2)^{\top}$ be a random vector following a bivariate symmetric (BSY) distribution fkn:90 with location (mean) vector $\boldsymbol{\mu}=(\mu_1,\mu_2)^{\top}$, covariance matrix
and density generator $g_c$, with $\mu_i\in\mathbb{R}$, $\sigma_i>0$, for $i=1,2$. We use the notation $\boldsymbol{Y}\sim {\rm BSY}(\boldsymbol{\mu},\boldsymbol{\Sigma},g_c)$. Then, the probability density function (PDF) of $\boldsymbol{Y}\sim {\rm BSY}(\boldsymbol{\mu},\boldsymbol{\Sigma},g_c)$ is given by
{ where $|\boldsymbol{\Sigma}|=\sigma_1^2\sigma_2^2(1-\rho^2)$ and $Z_{g_c}$ is a normalization constant so that $f_{\boldsymbol{Y}}$ is a PDF, that is,
} The density generator $g_c$ in (ref) leads to different bivariate symmetric distributions, which may contain an extra parameter (or extra parameter vector).
We propose a generalization of the classical Heckman-normal model heckman76 by considering independent errors terms following a BSY distribution with regression structures for the sample selection bias ($0<\rho<1$) and dispersion ($\sigma>0$) parameters:
In the above equation $\mu_{1i}, \mu_{2i}, \sigma_i$ and $\rho_i$ are are the mean, dispersion and correlation parameters, respectively, with the following regression structure $g_1(\mu_{1i}) = \boldsymbol{x}_i^\top \boldsymbol{\beta} $, $g_2(\mu_{2i}) = \boldsymbol{w}_i^\top \boldsymbol{\gamma} $, $h_1(\sigma_{i}) = \boldsymbol{z}_i^\top \boldsymbol{\lambda}$ and $h_2(\rho_{i}) = \boldsymbol{v}_i^\top \boldsymbol{\kappa}$, where $\boldsymbol{\beta} = (\beta_1, \ldots, \beta_k)^\top \in \mathbb{R}^{k}$, $\boldsymbol{\gamma} = (\gamma_1, \ldots, \gamma_l)^\top \in \mathbb{R}^{l}$, $\boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_p)^\top \in \mathbb{R}^{p}$ and $\boldsymbol{\kappa} = (\kappa_1, \ldots, \kappa_q)^\top \in \mathbb{R}^{q}$ are vectors of regression coefficients, ${\boldsymbol{x_i}}=(x_{i1},\ldots,x_{ik})^{\top} $, ${\boldsymbol{w_i}}=(w_{i1},\ldots,w_{il})^{\top}$, ${\boldsymbol{z_i}}=(z_{i1},\ldots,z_{ip})^{\top}$ and ${\boldsymbol{v_i}}=(v_{i1},\ldots,v_{iq})^{\top}$ are the values of $k$, $l$, $p$ and $q$ covariates, and $k+l+p+q < n$. The links $g_1(\cdot), g_2(\cdot), h_1(\cdot)$ and $h_2(\cdot)$ are strictly monotone and twice differentiable. The link functions $g_1: \mathbb{R}\rightarrow \mathbb{R}$, $g_2: \mathbb{R} \rightarrow \mathbb{R}$, $h_1: \mathbb{R}^+ \rightarrow \mathbb{R}$ and $h_2: [-1,1] \rightarrow \mathbb{R}$ must be strictly monotone, and at least twice differentiable, with $g_1^{-1}(\cdot)$, $g_2^{-1}(\cdot)$, $h_1^{-1}(\cdot)$, and $h_2^{-1}(\cdot)$ being the inverse functions of $g_1(\cdot)$, $g_2(\cdot)$, $h_1(\cdot)$, and $h_2(\cdot)$, respectively. For $g_1(\cdot)$ and $g_2(\cdot)$ the most common choice is the identity link, whereas for $h_1(\cdot)$ and $h_2(\cdot)$ the most common choices are logarithm and arctanh (inverse hyperbolic tangent) links, respectively.
We can agglutinate the information from $U_{i}^{*}$ in the following indicator function $U_i=\mathds{1}_{\{U_{i}^{*}>0\}}$.
Let $Y_{i}=Y_{i}^{*}U_i$ be the observed outcome, for $i=1,\ldots,n$. Only $n_1$ out of $n$ observations $Y_{i}^{*}$ for which $U_{i}^{*} > 0$ are observed. This model is known as “Type 2 tobit model” in the econometrics literature. Notice that $U_i\sim {\rm Bernoulli}(\mathbb{P}(U_{i}^*> 0))$. By using law of total probability, for $\boldsymbol{\theta}=(\boldsymbol{\beta}^{\top},\boldsymbol{\gamma}^{\top},\sigma,\rho)^{\top}$, the random variable $Y_i$ has distribution function
The function $F_{Y_i}$ has only one jump, at $y_i=0$, and $\mathbb{P}(Y_{i}=0)=\mathbb{P}(U_{i}^*\leq 0)$. Therefore, $Y_i$ is a random variable that is neither discrete nor absolutely continuous, but a mixture of the two types. In other words,
where $F_{\rm d}(y_i)=\mathds{1}_{[0,+\infty)}(y_i)$ and $F_{\rm ac}(y_i)=\mathbb{P}(Y_{i}^*\leq y_i|U_{i}^*>0)$. Hence, the PDF of $Y_i$ is given by
wherein $\mathbb{P}(U_i=0)=1-\mathbb{P}(U_i=1)=\mathbb{P}(U_{i}^*\leq 0)$ for $i=1,\ldots,n$, and $\delta_0$ is the Dirac delta function. That is, the density of $Y_i$ is composed of a discrete component described by the probit model $\mathbb{P}(U_i = u_i) = (\mathbb{P}(U_{i}^*\leq 0))^{1-u_i} (\mathbb{P}(U_{i}^*> 0))^{u_i} $, for $u_i=0,1$, and a continuous part given by the conditional PDF $f_{Y_{i}^*|U_{i}^*>0 }(y_i;\boldsymbol{\theta})$.
In the context of sample selection models, the interest lies in finding the PDF of ${Y}^*_i|{U}^*_i > 0$ given that $({Y}^*_i, {U}^*_i)^\top$ follows the BSY distribution (ref); see Theorem (ref).
Before stating and proving the main result (Theorem (ref)) of this section, throughout the paper we will adopt the following notations:
and
where
have the joint PDF $f_{Z_{1i},Z_{2i}}$ and $f_{X}$ denotes the PDF corresponding to a random variable $X$. Here, the random variables $V_{1i}$, $V_{2i}$, $R$, and $D$ are mutually independent and $\mathbb{P}(V_{ki} = -1) = \mathbb{P}(V_{ki} = 1) = 1/2$, $k=1,2$. The random variable $D$ in (ref) is positive and has PDF
The random variable $R$ in (ref) is positive and is called the generator of the random vector $(Y_i^*,U_i^*)^{\top}$. Moreover, $R$ has PDF given by
where $g_c$ is the density generator in (ref).
We now proceed to establish the main result of this section.
{
}
Given the density generator $g_c$ we can directly determine the PDF of ${Y}^*_i|{U}^*_i > 0$ as follows.
By combining Equation (ref) with Theorem (ref), the following formula for the PDF of $Y_i$ is valid:
where $\alpha_i=\rho_i/\oldsqrt[\ ]{1-\rho_i^2}$, $\tau_i=\mu_{2i}/\oldsqrt[\ ]{1-\rho_i^2}$, $u_i = 1$ if $u_i^{*}>0$ and $u_i = 0$ otherwise, $ g_1(\mu_{1i}) = \boldsymbol{x}_i^\top \boldsymbol{\beta} $, $g_2(\mu_{2i}) = \boldsymbol{w}_i^\top \boldsymbol{\gamma} $, $h_1(\sigma_{i}) = \boldsymbol{z}_i^\top \boldsymbol{\lambda}$ and $h_2(\rho_{i}) = \boldsymbol{v}_i^\top \boldsymbol{\kappa}$.
The log-likelihood of the symmetric generalized Heckman model for $\boldsymbol{\theta} = (\boldsymbol{\beta}^\top, \boldsymbol{\gamma}^\top, \boldsymbol{\lambda}^\top, \boldsymbol{\kappa}^\top)^{\top}$ is given by
To obtain the ML estimate of $\boldsymbol{\theta}$, we maximize the log-likelihood function (ref) by equating the score vector $\dot{\ell}(\boldsymbol{\theta})$ to zero, providing the likelihood equations. They are solved by means of an iterative procedure for non-linear optimization, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method.
The likelihood equations are given by {
} where
The generalized Heckman-normal model proposed by Bastos2021 is a special case of (ref) when the underlying distribution is bivariate normal. In this work, we focus on the generalized Heckman-$t$ model, which is based on the bivariate Student-$t$ (B$t$) distribution. This distribution is a good alternative in the symmetric family of distributions because it possesses has heavier tails than the bivariate normal distribution. From (ref), if $\boldsymbol{Y}=(Y_1,Y_2)^{\top}$ follows a B$t$ distribution, then the associated PDF is given by
where $\nu$ is the number of degrees of freedom. Here, the density generator of the B$t$ distribution is given by $g_c(x)=(1+x/\nu)^{-(\nu+2)/2}$ {, $|\boldsymbol{\Sigma}|=\sigma_i^2(1-\rho_i^2)$ and $Z_{g_c}= [\nu\pi\Gamma(\nu/2)]/\Gamma((\nu+2)/2)$ is a normalization constant.} Therefore, if $({Y}^*_i, {U}^*_i)$ follow a B$t$ distribution, then, by Corollary (ref), the PDF of ${Y}^*_i|{U}^*_i>0$ is written as
where $f_{\nu}$ and $F_\nu$ are the PDF and CDF, respectively, of a univariate Student-$t$ distribution with $\nu$ degrees of freedom, $\alpha_i=\rho_i/\oldsqrt[\ ]{1-\rho_i^2}$ and $\tau_i=\mu_{2i}/\oldsqrt[\ ]{1-\rho_i^2}$. The log-likelihood for $\boldsymbol{\theta} = (\boldsymbol{\beta}^\top, \boldsymbol{\gamma}^\top, \boldsymbol{\lambda}^\top, \boldsymbol{\kappa}^\top,\nu)^{\top}$ is given by
where $u_i = 1$ if $u_i^{*}>0$ and $u_i = 0$ otherwise, $\mu_{1i}$, $\mu_{2i}$, $\sigma_{i}$ and $\rho_{i}$ are as in (ref). The ML estimate of $\boldsymbol{\theta}$ is obtained by maximizing the log-likelihood function (ref), that is, by equating the score vector $\dot{\ell}(\boldsymbol{\theta})$ (given in Subsection (ref)) to zero, providing the likelihood equations. They are solved using an iterative procedure for non-linear optimization, such as the BFGS quasi-Newton method.
In this section, we carry out Monte Carlo simulation studies to evaluate the performance of the ML estimators under the symmetric generalized Heckman model. We focus on the generalized Heckman-$t$ model and consider three different set of true parameter value, which leads to scenarios covering moderate to high censoring percentages. The studies consider simulated data generated from each scenario according to
for $i = 1, \ldots, n$, $x_{1i}$, $x_{2i}$ and $x_{3i}$ are covariates obtained from a normal distribution in the interval (0,1). Moreover, the simulation scenarios consider sample size $n \in \{ 500, 1000, 2000\}$ and $\nu=4$, with $\text{NREP}=1000$ Monte Carlo replicates for each sample size. In the structure presented in (ref) - (ref), $\mu_{1i}$ is the primary interest equation, while $\mu_{2i}$ represents the selection equation. The R software has been used to do all numerical calculations; see rmanual.
The performance of the ML estimators are evaluated through the bias and mean squared error (MSE), computed from the Monte Carlo replicas as
where $\theta$ and $\widehat{\theta}^{(i)}$ are the true parameter value and its respective $i$-th ML estimate, and $\text{NREP}$ is the number of Monte Carlo replicas.
We consider the following sets of true parameter values for the regression structure in (ref)-(ref):
To keep the censoring proportion around 50%, in Scenario 1 a threshold greater than zero was used, so $U_{i}^{*} > a$. According to Bastos2021, in general the value of $a$ is zero, as any other value would be absorbed by the intercept, so considering another value does not cause problems for the model. In Scenario 2, the dispersion and correlation parameters were changed and the censoring proportion was maintained around 30%. In Scenario 3, the censoring rate around 50% was obtained by changing the parameters of the selection equation $\mu_{2i}$.
The ML estimation results for the Scenarios 1), 2) and 3) are presented in Tables (ref)-(ref), respectively, wherein the bias and MSE are all reported. As the ML estimators are consistent and asymptotically normally distributed, we expect the bias and MSE to approach zero as $n$ grows. Moreover, we expect that the performances of the estimates deteriorate as the censoring proportion (%) grows. A look at the results in Tables (ref)-(ref) allows us to conclude that, as the sample size increases, the bias and MSE both decrease, as expected. In addition, the performances of the estimates decrease when the censoring proportion increases.
In this section, two real data sets, corresponding to outpatient expense and investments in education, are analyzed. The outpatient expense data data set has already been analyzed in the literature by Heckman models Genton2012, whereas the education investment data is new and is analyzed for the first time here.
In this subsection, a real data set corresponding to outpatient expense from the 2001 Medical Expenditure Panel Survey (MEPS) is used to illustrate the proposed methodology. This data set has information about the cost and provision of outpatient services, and is the most complete coverage about health insurance in the United States, according to the Agency for Healthcare Research and Quality (AHRQ).
The MEPS data set contains information collected from 3328 individuals between 21 and 64 years. The variable of interest is the expenditure on medical services in the logarithm scale ($Y_{i}^{*}=lnambx$), while the latent variable ($U_{i}^{*} = dambexp$) is the willingness of the individual to spend; $U_i=\mathds{1}_{\{U_{i}^{*}>0\}}$ corresponds the decision of the individual to spend. It was verified that 526 (15.8%) of the outpatient costs are identified as zero (censored). The covariates considered in the are: $age$ is the age measured in tens of years; $fem$ is a dummy variable that assumed value 1 for women and 0 for men; $educ$ is the years of education; $blhisp$ is a dummy variable for ethnicity (1 for black or Hispanic and 0 if non-black and non-Hispanic); $totcr$ is the total number of chronic diseases; $ins$ is the insurance status; and $revenue$ denotes the individual income.
Table (ref) reports descriptive statistics of the observed medical expenditures, including the minimum, mean, median, maximum, standard deviation (SD), coefficient of variation (CV), coefficient of skewness (CS) and coefficient of (excess) kurtosis (CK) values. From this table, we note the following: the mean is almost equal to the median; a very small negative skewness value; and a very small kurtosis value. The symmetric nature of the data is confirmed by the histogram shown in Figure Figure (ref)(a). The boxplot shown in Figure (ref)(b) indicates some potential outliers. Therefore, we observe that a symmetric distribution is a reasonable assumption, more specifically a Student-$t$ model, since since we have to accommodate outliers.
We then analyze the medical expenditure data using the generalized Heckman-$t$ model, expressed as
We initially compare the adjustments of the generalized Heckman-$t$ (GH$t$) model, in terms of Akaike (AIC) and Bayesian information (BIC), with the adjustments of the classical Heckman-normal (CHN) heckman76 and generalized Heckman-normal (GHN) Bastos2021 models; see Table (ref). The AIC and BIC values reveal that the GH$t$ model provides the best adjustment, followed by the GHN model.
Table (ref) presents the estimation results of the GHN and GH$t$ models. From this table, we observe the following results: the explanatory variables $totchr$ and $ins$ that model the dispersion are significant, in both models, indicating the presence of heteroscedasticity in the data. The explanatory variable $age$ is only significant in the GHN model. For the correlation term, the covariates $fem$ and $totchr$ are significant for both models, which indicates the presence of selection bias in the data.
In the outcome equation, when we look at both models, $age$, $fem$, $blhisp$ and $totchr$ are significant at the 5% level, and $educ$ is not significant. The explanatory variable $ins$ is significant at the 5% and 10% levels in the GHN and GH$t$ models, respectively; see Table (ref). We can interpret the estimated coefficients in terms of the effect on the expenditure on medical services; see weisberg:14. For example, a 1-year increase in $age$ rises by $(\exp(0.1838)-1)\times100 = 20.18\%$ and $(\exp(0.1895)-1)\times100 = 20.86\%$ the expected value of the expenditure on medical services according to the GHN and GH$t$, respectively. Moreover, a 1-unit increase in $totchr$ rises by $(\exp(0.4306)-1)\times100 = 53.82\%$ and $(\exp(0.4464)-1)\times100 = 56.27\%$ the expected value of the response according to the GHN and GH$t$ models, respectively.
In the case of the selection equation, we observe that, for both models, the explanatory variables $age$, $fem$, $educ$, $blhisp$, $totchr$ and $ins$ are significant at the 5% level and $revenue$ is significant at the 10% level; see Table (ref). The interpretation is made in terms of odds ratio, which is obtained by exponentiating the estimated explanatory variable coefficient. For example, for the GH$t$ case, the odds ratio for $age$ is $\exp(0.0930) = 1.0975$, suggesting that each additional year of age raises the likelihood of an individual having expenditures on medical services by $ ((1.0975 - 1)\times 100)=9.75\%$.
Figure (ref) displays the quantile versus quantile (QQ) plots of the martingale-type (MT) residuals for the GHN and GH$t$ models. This residual is given by
where $r^{_{\textrm{\tiny M}_{i}}} = u_i + \log(\widehat S(t_i))$, $\widehat S(t_i)$ is the fitted survival function, and $u_i=0$ or $1$ indicating that case $i$ is censored or not, respectively; see tgf:90. The MT residual is asymptotically standard normal, if the model is correctly specified whatever the specification of the model is. From Figure (ref), we see clearly that the GH$t$ model provides better fit than GHN model.
In this subsection, data on education investments are used to illustrate the proposed methodology. We consider the investments made by municipalities of two Brazilian states: Sao Paulo (SP) and Minas Gerais (MG). This data set was obtained from the Brazilian National Fund for Educational Development (FNDE) website\footnotemark[1], which is a federal agency under the Ministry of Education. The origin of these investments comes from the Fund for the Maintenance and Development of Basic Education and the Valorization of Education Professionals (FUNDEB). The resources are distributed to 27 federative units (26 states plus the Federal District), according to the number of students enrolled in their basic education network. This rule is established for the previous year's school census data, e.g., 2018 resources were based on 2017 student numbers. This method helps to better distribute resources across the country, as it takes into account the size of education networks.
The variable of interest is education investments with 1503 observations, of which 102 (7%) correspond to unobserved investment values identified as zero investment. The explanatory variables considered in the study were: $revenue$\footnotemark[3] represents per capita revenue collected by the municipality; $gnp$ \footnotemark[3] is the Gross National Product of the municipality; $distribute$\footnotemark[2] is a dummy variable indicating if the municipality receives resources from the Financial Compensation for Exploration of Mineral Resources (CFEM); this resource must be destined to investments in the areas of health, education and infrastructure for the community; $sp$ is an indicator variable for state ($sp$ receives value 1); $enrollment$\footnotemark[4] is the school census enrollment numbers.
As in the previous study, the response variable, investment in education, is in the logarithm scale $Y_{i}^{*}=lninvest$. The latent variable ($U_{i}^{*} = dinvest$) denotes the willingness of the $i$th municipality to invest education; $U_i=\mathds{1}_{\{U_{i}^{*}>0\}}$ corresponds to the decision or not of the $i$th municipality to invest in education.
\footnotetext[1]{Filtered data are available at \hbox{https://repositorio.shinyapps.io/plataforma_de_dados_municipais}.}
\footnotetext[2]{\hbox{https://dados.gov.br/dataset/sistema-arrecadacao}.}
\footnotetext[3]{\hbox{http://www.ipeadata.gov.br/Default.aspx}.}
\footnotetext[3]{\hbox{https://www.ibge.gov.br/estatisticas/sociais/populacao/9103-estimativas-de-populacao.html?=&t=resultados}.}
\footnotetext[4]{\hbox{https://www.gov.br/inep/pt-br/areas-de-atuacao/pesquisas-estatisticas-e-indicadores/censo-escolar}.}
The descriptive statistics for the investments in education are reported in Table (ref). From this table, we note the following: the mean is almost equal to the median; a very small negative skewness value, and a high kurtosis value. The symmetric nature of the data is confirmed by the histogram shown in Figure Figure (ref)(a). The boxplot shown in Figure (ref)(b) indicates potential outliers. Therefore, we observe that a Student-$t$ model is a reasonable assumption.
We then analyze the education investment data using the GH$t$ model, expressed as
Table (ref) reports the AIC and BIC of the CHN, GHN and GH$t$ models. From Table (ref), we observe that the GH$t$ model provides better adjustment than other models based on the values of AIC and BIC.
Table (ref) presents the estimation results of the GHN and GH$t$ models. From this table, we note that the explanatory variables $revenue$ and $sp$ associated with the dispersion are significant, in both models, suggesting the presence of heteroscedasticity. For the explanatory variables related to the correlation parameter, $revenue$ is significant at the 5% and 10% levels in the GHN and GH$t$ models, respectively, while $distribute$ is significant at the 10% and 5% levels in the GHN and GH$t$ models, respectively. Therefore, there is evidence of the presence of sample selection bias in the data.
From Table (ref), we observe that in the outcome equation $revenue$, $sp$, $enrollment$ and $gnp$ are significant, according to the GHN and GH$t$ models. Note that increasing $revenue$ by one unit is associated with $(\exp(-0.2683)-1)\times100 = -23.53\%$ and $(\exp(-0.2167)-1)\times100 = -19.48\%$ decrease in the average investment in education according to the GHN and GH$t$ models, respectively. Note also that when the municipality is located in the state of Sao Paulo ($sp$), the average investment in education increases by $(\exp(1.0714)-1)\times100 = 191.95\%$ and $(\exp(1.0063)-1)\times100 = 173.55\%$, according to the GHN and GH$t$ models, respectively. These results show that the municipalities of the state of Sao Paulo invest much more in education than the municipalities of the state of Minas Gerais.
In the selection equation, we observe that the explanatory variables $revenue$, $sp$ and $enrollment$ are significant in the GH$t$ model. On the other hand, $revenue$ and $sp$ are not significant in the GHN model; see Table (ref). We observe that, for the GH$t$ case, the odds ratio for $sp$ is $\exp(0.3674) = 1.4439$, that is, the likelihood of a municipality to spend on medical services increases by $((1.4439 -1) \times 100)=44.39\%$ when the municipality is located in the state of Sao Paulo
Figure (ref) displays the QQ plots of the MT residuals. This figure indicates that the MT residuals in the GH$t$ model shows better agreement with the reference distribution.
In this paper, a class of Heckman sample selection models were proposed based symmetric distributions. In such models, covariates were added to the dispersion and correlation parameters, allowing the accommodation of heteroscedasticity and a varying sample selection bias, respectively. A Monte Carlo simulation study has showed good results of the parameter estimation method. We have considered high/low censoring rates and the presence of strong/weak correlation. We have applied the proposed model along with some other two existing models to two data sets corresponding to outpatient expense and investments in education. The applications favored the use of the proposed generalized Heckman-$t$ model over the classical Heckman-normal and generalized Heckman-normal models. As part of future research, it will be of interest to propose sample selection models based on skew-symmetric distributions. Furthermore, the behavior of the Wald, score, likelihood ratio and gradient tests can be investigated. Work on these problems is currently in progress and we hope to report these findings in future.