Model specification

Our modeling strategy is intended to model the level y of expenditures of a household for a given good or service during a given period of observation. To this purpose, we use up to three zero expenditure generating mechanisms, called hurdles, and a demand function.

Each hurdle is represented by a qualitative binary response model resting on one of the following three latent dependent variables relations:

where y₁^*, y₂^*, y₃^* stand for continuous latent dependent variables x₁, x₂, x₃ for column-vectors of explanatory variables (called covariates in the followings), β₁, β₂, β₃ for column-vectors of the impact coefficients of the explanatory variables on the latent dependent variables, z₁, z₂, z₃ for standard normal random disturbances. Since variables y₁^* and y₃^* are never observed, contrary to y₂^*, the units of measurement of y₁^* and y₃^* are not identified. Hence, the disturbances in the corresponding binary response models are normalized by setting their variances equal to 1.

Relying on a Johnson (1949)’s proposal, T(⋅) denotes a differentiable one-to-one transformation intended to rescale a potentially not normally distributed random variable y₂^* into a normally distributed ones, with σ standard deviation. As the domain of definition of the inverse transformation T⁻¹(⋅) my be restrained to a subset of real numbers, the standard normal random disturbance z₂ is potentially truncated at the bounds of an interval [B₁, B₂] with B₁ and/or B₂ finite. To account for heteroscedasticity of z₂, in mhurdle the standard deviation σ can be specified as a positive monotonic transformation of a linear function of a vector of covariates x₄, namely β₄^⊤x₄. However, for notational simplicity, we will adopt in what follows the notation of an homoscedastic standard deviation σ rather than the more general one of an heteroscedastic standard deviation σ(β₄^⊤x₄). ³

These equations model the hurdles, namely the sequence of binary decisions a consumer faces when buying a good. Each hurdle decision is coded by a binary variable I_j taking value 1 if y_j^* > 0 and 0 otherwise. Therefore, the outcome of each of such decisions is modeled as a probit model with outcome probabilities: P(I_j = 1) = P(y_j^* > 0) and P(I_j = 0) = 1 − P(y_j^* > 0).

By assuming that consumption and purchases are uniformly distributed over time, but according to different timetables entailing a frequency of consumption higher than that of purchasing, we can also interpret the probability P(I₃ = 1) = P(y₃^* > 0) as measuring the share of purchasing frequency to that of consumption during the observation period. This allows to relate the observed level of expenditures y to the unobserved level of consumption y₂^* during the observation period, using the following identity:

This modeling strategy provides an explanatory framework where censoring of the observed dependent variable results from the effect of up to three different censoring mechanisms, namely: I₁ = 0 and/or I₂ = 0 and/or I₃ = 0.

A priori information (theoretical or real-world knowledge) may suggest that one or more of these censoring mechanisms are not in effect. For instance, we know in advance that all households purchase food regularly, implying that the first two censoring mechanisms are inoperative for food. In this case, the relevant model is defined by only two relations: one defining the desired consumption level of food and the other the decision to purchase food during the observation period.

When the first and/or the third hurdles are supposed not to be in effect, the previous binary response models for I₁ and/or I₃ must be replaced by singular probability response models where P(I₁) = 1 and/or P(I₃) = 1. When the second hurdle is supposed to be inoperative, we must not only replace the binary response model explaining I₂ by a singular probability response model where P(I₂ = 1) = 1 but also select a monotonic transformation T(y₂^*) enforcing non-negative values to y₂^*.

Cragg (1971) suggested two types of functional forms of the distribution of y₂^* enforcing such constraint, namely the (inverse) log-normal distribution:

and the normal distribution left-truncated at 0:

where max {B₁; z₂} is a standard normal random variable left-truncated at B₁ = −β₂^⊤x₂/σ.

To account for these two functional forms within our general formal framework @ref(eq:threeeqs) and illustrate in section 5 their empirical use with mhurdle, in this paper we consider two transformations T(y₂^*).⁴

The first transformation is the identity transformation which generates an un-truncated normal distribution of y₂^* when hurdle 2 is in effect, and the left-truncated normal distribution @ref(eq:trunc) when hurdle 2 is not in effect.

The second transformation is a slight generalization of the logarithmic transformation leading to the (inverse) log-normal distribution @ref(eq:lognorm), namely:

with α a location parameter which shifts the left-truncation bound of y₂^* towards negative or positive values according to the sign of α. As the inverse of this transformation is written as:

it turns out that this transformation allows to model skewed distributions of y₂^* according to a translated (inverse) log-normal functional form where the support of z₂ is un-truncated. By the same token, testing the statistical significance of parameter α, against the alternative α > 0, amounts to testing the assumption that hurdle 2 is not in effect. Conversely, when α < 0 the transformation @ref(eq:lognorm2) holds for y₂^* > −α meaning that −α stands for a “committed” consumption of a basic necessity, while α > 0 typify a luxury good whose consumption occurs only above a given income threshold. ⁵

To derive the form of the probability distribution of the observable dependent variable y, we must still specify the joint distribution of the random disturbances entering the hurdle relations assumed to be in effect and that of the demand function, which in turn depends on the domain of definition of transformation T⁻¹(⋅).

For trivariate hurdle models, where hurdles 1 and 3 are in effect (but not necessarily hurdle 2), the joint density function of z₁, z₂ and z₃ is a possibly truncated standard trivariate normal density function written as:

where Π = Φ(B₂) − Φ(B₁) with Φ(z) the distribution function of a standard univariate normal distribution, ϕ(z₁, z₂, z₃; ρ) the density function of a standard trivariate normal distribution with ρ the vector of the three symmetric correlation coefficients ρ₁₂, ρ₁₃, ρ₂₃ between the couples of normal standard random variables z₁ and z₂, z₁ and z₃, z₂ and z₃, respectively.

For bivariate hurdle models, where either hurdle 1 or hurdle 3 is not in effect, the joint density function of z_i, with i = 1 or 3, and z₂ is a possibly truncated standard bivariate normal density function written as:

where ϕ(z_i, z₂; ρ_i2) denotes the density function of a standard bivariate normal distribution.

Finally, for univariate hurdle models, where both hurdles 1 and 3 are not in effect, the density function of z₂ is a possibly truncated standard univariate normal density function:

where ϕ(z₂) denotes the density function of a standard univariate normal distribution.

While the assumption of correlated disturbances is intended to account for the interdependence between latent variables y₁^*, y₂^* and y₃^* unexplained by covariates x₁, x₂ and x₃, a priori information (theoretical or real-world knowledge) may also suggest to set to zero some or all correlations between the random disturbances entering these models, entailing a partial or total independence between model relations. The use of this a priori information generates, for each trivariate or bivariate hurdle model, a subset of special models all nested within the general model from which they are derived. For a trivariate hurdle model the number of special models so derived is equal to seven, but for a bivariate hurdle model only one special model is generated, namely the model obtained by assuming the independence between the two random disturbances of the model.

Probability distribution of censored dependent variable

As for the standard Tobit model, the probability distribution of the observed censored variable y of our hurdle models is a discrete-continuous mixture, which assigns a probability mass P(y = 0) to y = 0 and a density function f₊(y) to any y > 0, with:

To compute the probability mass P(y = 0) = 1 − P(y > 0) and the density function f₊(y) we must first establish the joint density function of the latent variables entering the hurdle model.

For trivariate hurdle models, the joint density function of latent variables y₁^*, y₂^* and y₃^* can be derived from the density function @ref(eq:trivdens) of variables z₁, z₂ and z₃ by means of the change of variables:

which leads to:

where T^′(y₂^*) stands for the derivative of T(y₂^*).

To compute P(y > 0), we integrate this density function over the three-dimensional positive octant using the change of variables @ref(eq:threez) and the symmetry property of density function @ref(eq:trivdens). Performing this integration for a strictly increasing transformation T(⋅)⁶, which is the case for the identity and logarithmic transformations we focus in this paper, and z₂ truncated over the interval [B₁, B₂] leads to:

Therefore, for z₂ left-truncated at B₁ = −β₂^⊤x₂/σ, which is the case for the truncated normal distribution @ref(eq:trunc), this closed form for P(y > 0) simplifies to:

while for the un-truncated normal and the translated (inverse) log-normal distribution @ref(eq:lognorm2i) it is written as:

To compute the density function f₊(y) we first derive the joint density function of variables y₁^*, y and y₃^* by performing the change of variable y₂^* = Φ₃y with Φ₃ = Φ(β₃^⊤x₃) = P(I₃ = 1) on the joint density function @ref(eq:threedens), which leads to :

We must then integrate this transformed density function over the positive quadrant of the bi-dimensional support of variables y₁^* and y₃^*. To this purpose we rewrite this trivariate normal density function as the product of the marginal density function of y, written as:

and of the conditional joint density function f(y₁^*, y₃^*|y), which is bivariate normal with expectations, standard-deviations and correlation coefficient written as:

Integrating this factorization of density function @ref(eq:denstriv) with respect to positive values of y₁^* and y₃^* leads to the following closed form for f₊(y):

The probability distribution of the observed censored variable y for bivariate and univariate hurdle models can be derived from that of the trivariate model by eliminating hurdles 1, 3 and 1 and 3, respectively. By eliminating hurdle 1 we obtain:

and

while the elimination of hurdle 3 leads to:

and

Finally, removing both hurdles 1 and 3 leads to the univariate censured regression model with:

and

The econometric framework described in the previous section provides a theoretical background for tackling the problems of model estimation, evaluation, selection and prediction within the statistical theory of classical inference. In mhurdle these statistical issues are tackled by assuming that data at hand are those observed on a sample of individuals selected in a large population using a sampling design generating a data base having the features of a random sample. ⁷

To appraise the results of a model estimation two fundamental principles should be used, namely its economic relevance and its statistical and predictive adequacy. The first principle deals with the issues of accordance of model estimate with the economic rationale underlying the model specification and of its relevance for answering the questions for which the model has been built. These issues are essentially context specific and, therefore, cannot be dealt with by means of generic criteria. The second principle refers to the issues of empirical soundness of an estimated model and of its ability to predict sample or out-of-sample observations. These issues can be tackled by means of measures of goodness of fit, by formal tests of significance based on the distribution of model parameter estimators, and by the quality of predictions provided by an estimated model, respectively.

Model evaluation and selection using goodness of fit measures

To assess the goodness of fit of a mhurdle model estimate we can not rely on a single measure, as a sample of censored observations of dependent variable y consists of two sub-samples of incomparable observations, namely the sub-sample of zero observations and that of positive observations. The former is a sample of observations of a qualitative binary variable coding the outcome of the censoring mechanisms at work, while the latter is a sample of observations of a quantitative variable measuring the outcome of the demand function, when this demand is uncensured. As a consequence, a measure of goodness of fit for each of these two samples of observations must be used as the goodness of fit is not measured the same way in qualitative binary response models as in quantitative response models.

To this purpose, we reformulate a mhurdle model as a two-part model. The first-part specifies a model for the censoring mechanism explaining the outcome of the qualitative binary variable:

according to a probit model assigning the probability P(y = 0) to the censored observation (d = 0) and the probability P(y > 0) = 1 − P(y = 0) to the uncensored observation (d = 1). Therefore, the log-likelihood for this first-part model is written as:

The second-part specifies, conditionally to an uncensored response being observed, the demand model explaining the outcome of quantitative variable y according to the following conditional density function:

For this second-part model, the log-likelihood is written as:

Adding-up the log-likelihoods of these two-part model leads to:

meaning that the two-part model is entirely consistent with an integrated specification of mhurdle model.

To develop a measure of fit for binary response models which may be intuitively interpreted in a similar way to the coefficient of determination R² in the linear regression model, Estrella (1998) suggests to rely on the relationship between R² and the classical likelihood-based test statistics, namely the likelihood ratio (LR), Wald (W) and Lagrange multiplier/score (LM) statistics. In the normal linear regression model, where such relationships holds exactly for the likelihood ratio statistics¹⁰, it is written as:

with A = −(2/n)ln (L_c/L_u) the average value of the likelihood ratio statistic, L_u the maximum value of the unconstrained likelihood, and L_c the maximum value of the likelihood constrained by the hypothesis of a model without covariates (intercept-only model).

From this relationships, it follows that R² = 0 for A = 0, which corresponds to the case where the fit of the unconstrained model is identical to that of the constrained model (worst fit), and R² = 1 when A = ∞, the upper bound of the likelihood ratio statistic achieved when the fit of the unconstrained model is perfect (best fit). Between these two limits R² increases with A according to a rescaling of A determined by a differential equation written as:

These properties of relationships @ref(eq:eq33) were adopted by Estrella (1998) as a set of requirements for developing a R² analog of the measure of fit for a binary response model. Taking into account that, contrary to the case of a normal linear regression model, for a binary response model the average likelihood ratio statistic A varies between 0 (when L_u = L_c) and an upper finite bound B = −(2/n)ln L_c (when L_u = 1), a R² analog for a binary response model is obtained by solving the differential equation:

with the initial condition R²(0) = 0. This solution:

provides a rescaling of the likelihood ratio statistic for the binary response model that satisfies the condition R²(B) = 1 and converge towards the rescaling @ref(eq:eq33) when B → ∞.

To conclude, for assessing the goodness of fit of the two parts of a mhurdle model we shall use the following R²-rescaling of likelihoods L₀ and L₊:

with $n_+=\sum_{i=1}^n d_i$ the sample size of uncensored observations.

Note that the estimation of an intercept-only mhurdle model is not generally possible in its original parametrization because in the absence of covariates many parameters are under-identified, and thus cannot be consistently estimated. To compute the maximum values of the constrained likelihoods L_0c and L_+c, we must respecify the probability distribution of dependent variable y by means of identifiable parameters. To this end we observe from formula @ref(eq:pzero) that P(y > 0) is a constant P₊, and from formula @ref(eq:denspos) that f₊(y) is proportional to the density function of T(y) which is normal with a constant expectation μ and a constant variance s². Hence, the log-likelihood of the constrained first-part model is written as:

and provides the following maximum likelihood estimator for P₊: P̂₊ = n₊/n. Similarly, the log-likelihood of the constrained second-part model is written as:

leading to the following estimators for parameters μ and s²: $\hat{\mu} = \sum_{i=1}^n d_i T(y_i)/n_+$ and $\hat{s}^2=\sum_{i=1}^n d_i (T(y_i)-\hat{\mu})^2/n_+$.

The maximized values of constrained likelihoods L_0c and L_+c are obtained by inserting these estimators for parameters P₊, μ and s² into formulas @ref(eq:eq29a) and @ref(eq:eq31a). Note that for the logarithmic transformation @ref(eq:lognorm2) with a non-zero location parameter α, the constrained log-likelihood @ref(eq:eq31a) must also be maximized with respect to this parameter. This maximization may be performed effectively by considering μ̂ and ŝ² as functions of α, and maximizing the concentrated log-likelihood obtained by inserting these conditional estimators into constrained log-likelihood @ref(eq:eq31a), which becomes a function of the single parameter α written as:

Goodness of fit measures @ref(eq:eq37) can also be used for model selection, a decision problem dealing with the identification, among several model specifications used to explain the same dependent variable, of the one that is best suited to explain the sample of available observations. From this point of view, goodness of fit measures @ref(eq:eq37) allow to identify the model specification achieving the best in-sample fit. To be fair, the goodness of fit measure used for model selection must account for differences in the degree of model parametrization. Indeed, the value of the above goodness of fit measures can be improved by increasing model parametrization, because model estimation is performed by maximizing the model likelihood function from which functionally depend the goodness of fit measures we consider. Consequently, a penalty that increases with the number of model parameters should be added to formulas @ref(eq:eq37). To this purpose, we shall rely on THEIL(1971)’s correction of the coefficient of determination in the linear regression model. This leads to the following adjusted R²-rescaling of likelihoods L₀ and L₊:

with K and K₀ the number of parameters in the unconstrained model and the number of parameters in the constrained second-part model, respectively. Such adjustments may not be relevant in large samples but may be relevant in small samples.

Model selection using Vuong tests

Model evaluation can also be tackled from the point of view of the model specification that is favored in a formal test comparing two model alternatives.

This second model selection criterion relies on the use of a test proposed by Vuong (1989). According to the rationale of this test, the “best” parametric model specification among a collection of competing specifications is the one that minimizes the Kullback-Leibler Information Criterion (KLIC), namely a measure of the distance between the conditional probability mass/density function¹¹ f(y|x; θ) of a possibly misspecified parametric model and that of the true unknown model denoted by h(y|x), defined by the following formula:

where H(y, x) denotes the distribution function of the true joint distribution of (y, x) and θ_* the probability limit, with respect to H(y, x), of the so called quasi-maximum likelihood estimator θ̂ obtained by applying the maximum likelihood method when f(y|x; θ) is misspecified. Note that KLIC criterion takes a minimum value of 0 when there is a value θ₀ of parameter vector θ such that f(y|x; θ₀) = h(y|x) almost everywhere¹², meaning that conditional probability function f(y|x; θ) is correctly specified. At the opposite, KLIC will take large values when f(y|x; θ) very poorly specifies the true conditional probability function h(y|x).

As the application of this model selection criterion equals selecting the model specification for which the quantity:

is the largest, given two competing models with conditional probability functions f(y|x; θ) and g(y|x; π) and parameter vectors θ and π of size K and L, respectively, Vuong suggests to discriminate between these models by testing the null hypothesis:

$$ H_0 : E[\ln f(y|x;\theta_\ast)]=E[\ln g(y|x;\pi_\ast)] \Longleftrightarrow E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]=0, $$

meaning that the two models are equivalent, against either:

$$ H_f : E[\ln f(y|x;\theta_\ast)]>E[\ln g(y|x;\pi_\ast)] \Longleftrightarrow E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]>0, $$

meaning that specification f(y|x; θ) is better than g(y|x; π), or:

$$ H_g : E[\ln f(y|x;\theta_\ast)]<E[\ln g(y|x;\pi_\ast)] \Longleftrightarrow E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]<0, $$

meaning that specification g(y|x; π) is better than f(y|x; θ).

The quantity E[ln f(y|x; θ_*)] is unknown but it can be consistently estimated, under some regularity conditions, by 1/n times the log-likelihood evaluated at the quasi-maximum likelihood estimator for θ. Hence, 1/n times the log-likelihood ratio statistic:

is a consistent estimator for E[ln (f(y|x; θ_*)/g(y|x; π_*))].

Therefore, an obvious test of H₀ consists in verifying whether a LR-based statistic differs from zero. The distribution of such a statistic can be worked out even when the true model is unknown, as the quasi-maximum likelihood estimators θ̂ and π̂ converge in probability to the pseudo-true values θ_* and π_*, respectively, and have asymptotic normal distributions centered on these pseudo-true values.

In his seminal paper Vuong (1989) suggests to use different LR-based statistic according to the relation linking the two competing models, which may be nested, strictly non-nested or overlapping models ¹³.

In a recent paper Shi (2015) has criticized Vuong testing procedures for discriminating between competing models by showing that they can generate severe size distortions in finite samples. She develops a new test statistic that corrects these size distortions and which universally applies to nested, strictly non-nested and overlapping models. The rationale of Shi test statistic consists in correcting two finite sample size distortions present in the test statistic used by Vuong to discriminate two strictly non-nested models, namely:

with ω̂² a sample analogue consistent estimator of the variance of the asymptotic distribution of the numerator of statistic T_LR, which is N(0, ω²) under the null hypothesis H₀.

The first size distortion is due to a bias of order O(1/n) present in the expectation of the asymptotic distribution of the numerator of T_LR under H₀. This expectation is written as −tr(V)/(2n) with tr(V) = tr(A_F⁻¹B_F) − tr(A_G⁻¹B_G), where A_F and B_F are the limits of the Hessian and of the gradient version of the Fisher information matrix of model F_θ, respectively, and A_G and B_G those of model G_π. Therefore, for a model for which the information identity A + B = O holds δ = −tr(A⁻¹B) = tr(I) is equal to the number of model parameters. This suggests to consider δ as a measure of the degree of model parametrization which increases with the number of model parameters and with the discrepancy of the information identity resulting from the nonlinearity of the model with respect to its parameters. As a consequence, the finite sample bias of statistic LR(θ̂, π̂)/n is given by (−δ_F + δ_G)/(2n). It will be negligible for models of similar degrees of parametrization but negative and of possibly significant magnitude when model F_θ is more heavily and non-linearly parameterized than model G_π. Hence, in order to avoid an over-rejection of H₀ in favor of the more parameterized model, Shi suggests to correct the Vuong statistic T_LR by subtracting to LR(θ̂, π̂) a consistent estimator of its asymptotic bias under H₀, namely −tr(V̂)/2 with V̂ a consistent estimator of matrix V.

The second adjustment of Vuong statistic T_LR suggested by Shi is intended to prevent the denominator ω̂ of this statistic from taking values close to zero with non negligible probability, which contributes to create a fat tail for the distribution of T_LR. This adjustment consists in adding a positive constant term to ω̂², written as ctr(V̂²/n) with c a positive constant chosen according to a data-dependent method ¹⁴, in order to avoid large values of T_LR when its denominator is close to zero while its numerator is not.

Finally the modified Vuong test by Shi, called by her a , uses the following statistic:

To completely remove the asymptotic size distortion, Shi proposes to compute a simulated critical value cv(a), which has a well controlled asymptotic size, in the sense that, asymptotically and for any c ≥ 0, the probability of rejecting H₀ when it is true is upper bounded by a ¹⁵. By systematically changing the value of a in order to equalize the simulated critical value cv(a) to the observed value of statistic T_LR^mod, it is possible to derive a simulated p-value allowing to perform a non-degenerate Vuong test at any controlled asymptotic size level.

For non-nested models the test is two sided. It is carried out by rejecting H₀ if the simulated p-value is less than the selected size level a, and in favor of either H_f if T_LR^mod > 0, or of H_g if T_LR^mod < 0.

For nested models the test is one sided because the nested model cannot be closer to the true unknown model than the nesting model according to the KLIC criterion. Supposing without loss of generality that model G_π is nested into model F_θ, the test is carried out by rejecting H₀ in favor of H_f if the simulated p-value is less than the selected size level a. Note that this test is a robust alternative to the usual LR-test for testing the statistical significance of a priori restrictions on the parameters of model F_θ, in the sense that it allows the unrestricted model to be misspecified.

A last comment deserves to be devoted to the choice of the size α of the test, which measures the probability of rejecting H₀ when it is true, and thus the risk of favoring, on grounds of statistical evidence, the use of one of two competing models while they are equivalent. This choice is routinely made by assuming a low value of α (5% or 1%) without referring to the operational meaning of the test result. But in a decision-making framework of selection between two models, the size of the test can be chosen as large as desired because if the two models are equivalent it is irrelevant to use one or the other of these two models. On the other end, the risk of the second kind, namely the probability of rejecting either H_f or H_g when one of these two models is true, must be minimized, because concluding that the two models are equivalent does not allow deciding which one of them should be used. Not rejecting H₀ simply point out that further sample information is needed to determine, with the desired degree of certainty, which of the two model alternatives is the more suitable for the intended uses.

But how to proceed if further data-gathering is impossible, which is the typical situation faced in econometric studies. In this situation the discrimination test between two models must be reformulated as a Simon (1943) symmetric test, namely a binary choice between models F_θ and G_π where none of these models is privileged as a null hypothesis of a classical asymmetric Neyman-Pearson test. This reformulation removes the irrelevant hypothesis of doubt H₀ from the decision-making framework by forcing the test to always conclude in favor of one of the two confronted models, which seems appropriate when there is no a priori evidence to consider one of the two models more plausible than the other¹⁶. The test is carried out by first computing the value of statistic T_LR or T_LR^mod, then by rejecting model G_π in favor of model F_θ when this value is positive, and conversely by rejecting F_θ in favor of G_π when the value is negative ¹⁷. The probability of committing an error with this symmetric test (accepting model F_θ when it is worse than G_π) ¹⁸ is a function of D = E[ln (f(y|x; θ_*)/g(y|x; π_*))], since the tested hypothesis is a composite one leaving unspecified the absolute value of the distance D between the distributions of models F_θ and G_π. This risk is maximum for D = 0 and decreases with |D|. The exact value of this risk can be found by integrating the density function of statistics T_LR or T_LR^mod over its negative values¹⁹.

Prediction and marginal effects

Prediction with a causal model like mhurdle of an observable dependent random variable z refers to defining a function of model covariates x providing, for any admissible value of x, a point estimate ẑ of z whose prediction error e = z − ẑ is minimized according to some criterion. In the frame of classical or frequentist statistical inference paradigm, an optimal conditional predictor of z, for a given value of covariates x, is defined as the one that minimizes the expectation of the squared conditional prediction error, namely E[e²|x] called the mean-square error of prediction. Using this criterion leads to select as optimal predictor for z its conditional expectation E[z|x], called the best mean-square error predictor for z.

As defined by identity @ref(eq:expcons), the dependent observable variable y of mhurdle is a mixture of a fqualitative and a quantitative variables whose expectation has no actual meaning, the zero-value of y denoting the numerical value taken by the binary dummy variable d when the observation of latent variable y₂^* is censored. Therefore, for a predictive purpose, one must disentangle the qualitative component of y from its purely quantitative component. This can be done by solving the prediction problem of y in the framework of the two-part respecification of mhurdle models presented in section 3.2.

The first-part of this model explains the outcome of the qualitative component of y using the binary variable d, whose best mean-square error predictor is given by:

with P(y > 0) specified by formulas @ref(eq:pzero), @ref(eq:pzero2), @ref(eq:pzero3), @ref(eq:pzerosm1), @ref(eq:pzerom3) and @ref(eq:pzerom13).

The second-part of the model explains the outcome of the quantitative component of y taking positive real values according to the density function @ref(eq:eq30). Accordingly, the best mean-square error predictor of y > 0 is given by:

with f₊(y) specified by formulas @ref(eq:denspos), @ref(eq:densposm1), @ref(eq:densposm3) and @ref(eq:densposm13).

For the special cases of transformations T(⋅) identity and logarithmic focused in this paper, closed analytical forms for E[y|x] can be established²⁰.

For the identity transformation T(y₂^*) = y₂^*, this analytical closed form for trivariate hurdle models with density function f₊(y) specified by formula @ref(eq:denspos) is written as:

where

$$ \begin{array}{l} {\displaystyle \Phi_{123}=\Phi\left(\beta_1^\top x_1, \frac{\beta_2^\top x_2}{\sigma}, \beta_3^\top x_3;\rho\right), \quad \phi_i=\phi(\beta_i^\top x_i), \quad i=1,3 , \quad \phi_2=\phi\left(\frac{\beta_2^\top x_2}{\sigma}\right), } \\[6pt] {\displaystyle \Phi_{ij\mid k} = \Phi\left(\frac{\mu_{i|k}}{\sigma_{i|k}},\frac{\mu_{j|k}}{\sigma_{j|k}};\rho_{ij|k}\right), \quad i\neq j\neq k,\quad \mu_{i|k}=\beta_i^\top x_i-\rho_{ik}\beta_k^\top x_k, \quad k=1,3 ,} \\[6pt] {\displaystyle \quad \mu_{i|2}=\beta_i^\top x_i-\rho_{i2}\frac{\beta_2^\top x_2}{\sigma}, \quad \sigma_{i|k}=\sqrt{1-\rho_{ik}^2}, \quad \rho_{ij|k}=\frac{\rho_{ij}-\rho_{ik}\rho_{jk}}{\sigma_{i|k}\sigma_{j|k}}.} \end{array} $$

The closed forms for bivariate and univariate hurdle models with density function f₊(y) specified by formulas @ref(eq:densposm1), @ref(eq:densposm3) and @ref(eq:densposm13), are written as:

respectively.

Similarly, closed analytical forms for E[y|x] can be established in the case of the logarithmic transformation T(y₂^*) = ln (y₂^* + α). For trivariate, bivariate and univariate hurdle models with density function f₊(y) specified by formulas @ref(eq:denspos), @ref(eq:densposm1), @ref(eq:densposm3) and @ref(eq:densposm13) these closed forms are written, respectively:

with Φ₁₂₃ = Φ(β₁^⊤x₁, (β₂^⊤x₂ − ln α)/σ, β₃^⊤x₃; ρ),

The marginal effects we consider refer to the change in the value of predictors @ref(eq:eq47) and @ref(eq:eq48) generated by a ceteris paribus increase of a model covariate, the value of the other model covariates being kept unchanged, as in a laboratory experiment.

The definition of the marginal effect of a mhurdle covariate x₀ on a dependent variable of the model, depends on the scale of measure of x₀. When x₀ is a continuous quantitative variable measured either on an interval or on a ratio scale,²¹ the marginal effect of x₀ on predictor @ref(eq:eq47) or @ref(eq:eq48) is the first derivative of the predictor with respect to x₀, namely the rate of change in the value of the predictor generated by a ceteris paribus infinitesimal increase of x₀. Similarly, if x₀ is a count variable measuring the size of a population of individuals or objects which can only vary by integers, the marginal effect will be defined as the finite change in the predictor induced by a ceteris paribus unit increase of x₀.

When x₀ is a qualitative explanatory variable that can be in one and only one of K ≥ 2 different mutually exclusive and exhaustive situations called attributes, it is usual to code it as a column vector of K dummy variables x₀ = [d_k]_{k = 1, ..., K}, taking values d_h = 1 and d_k = 0, ∀k ≠ h when the observed attribute of x₀ is the one indexed by h. In this case, considering x₀ as the first (vector) element of one of mhurdle covariate vectors x_i, i = 1, ..., 4, the first term of the corresponding linear combination β_i^⊤x_i will take the value β_ih.

Yet, this dummy coding of a nominal explanatory variable ²² has an annoying drawback: because the K dummy variables are linearly related by the constraint $\sum_{k=1}^K d_k = 1$, the parameters β_ik, k = 1, ..., K are not identified and therefore not estimable when vector β_i does contain an intercept term. To avoid this inconvenience, it is suggested to suppress from a K-attribute dummy coding one of the K dummies, the discarded attribute playing the role of a “reference attribute” with respect to which the impact on β_i^⊤x_i of one of the K − 1 remaining attributes is measured. Indeed, substituting in $\sum_{k=1}^K \beta_{ik} d_k$ dummy variable d_K of reference attribute K by $d_K=1-\sum_{k=1}^K d_k$, we obtain $\sum_{k=1}^{K-1} (\beta_{ik}-\beta_{iK})d_k+\beta_{iK}$, showing that the impact on β_i^⊤x_i of attribute k is measured by parameter b_k = β_ik − β_iK called contrast with respect to attribute K, while the impact coefficient β_iK is included in vector β_i as an intercept, or incorporated to an already existing intercept of this vector. Consequently, the marginal effect of a change of x₀, from attribute h to attribute k, will be defined as the finite change in the predictor induced by a ceteris paribus change of the first term of β_i^⊤x_i from contrast b_ih to contrast b_ik, with obviously b_K = 0.

When x₀ is a nominal qualitative variable whose K attributes may be partially or totally ordered according to some meaningful criterion (for example the educational level measured and ranked according to the highest level of diploma awarded), it is possible to use a dummy coding allowing to assess the degree of monotonicity between a ceteris paribus marginal increment of x₀ from one attribute level to the next higher attribute level, and its impact on the dependent variable predictor. To this purpose, this impact generated by the marginal increment of x₀ from its attribute of rank k − 1 to its attribute of immediately higher rank k must be measured by the following marginal contrasts: γ_k = β_ik − β_{i, k − 1}, k = 2, ..., K. Hence, substituting in $\sum_{k=1}^K \beta_{ik} d_k$ parameters β_ik by the following solution of the former recurrence equations with respect to γ_k, k = 2, ..., K, namely $\beta_{i1}+\sum_{h=2}^k \gamma_h$ leads to the following recoding of the K-attribute dummy coding in terms of marginal contrasts: $\beta_{i1}+\sum_{k=2}^K \gamma_k \delta_k$ with $\delta_k=\sum_{h=k}^K d_h$. Consequently, the marginal effect of a change of x₀, from pre-ordered attributes k − 1 to k, will be defined as the finite change in the predictor induced by a ceteris paribus change of the first term of β_i^⊤x_i from marginal contrast γ_k − 1 to marginal contrast γ_k, with γ₁ = 0.

So far, punctual prediction and marginal effects for a mhurdle independent variable have been tackled by assuming known the value of the model parameters. However, to be operational these formulas must be quantified by means of an estimate of these theoretical parameters. Using the maximum likelihood estimate θ̂ of model parameters θ presented in section 3.1, provide an implied maximum likelihood estimator for a column-vector h(θ) of best mean-square error predictors @ref(eq:eq47), @ref(eq:eq48), and/or marginal effects of these predictors, namely h(θ̂). This estimator is, as θ̂, a consistent and asymptotically normal estimator for h(θ) whose asymptotic distribution can be derived from a linearization of h(θ̂) around θ called the “delta method”, which leads to the following asymptotic variance-covariance matrix of h(θ̂): (1/n)(∂h/∂θ^⊤)I_A(θ)⁻¹(∂h/∂θ^⊤)^⊤.

A last comment deserves to be devoted to how to present the estimates of marginal effects which are representative for the available sample. To this purpose it is customary to provide either the marginal effects computed at the sample average value of model covariates, or the average value of the sample of marginal effects computed at the model covariates values of the individuals of the sample. The first solution must be avoided when same model covariates are qualitative, because the sample average value of a dummy coded qualitative variable is no more a qualitative variable but a vector of proportions without a meaningful interpretation for a representative individual. Regarding the second solution, it provides an informative measure of location of the marginal distribution of a marginal effect, but that is not enough to characterize, by itself, the a priori unknown profileq of this distribution. Accordingly, it should be supplemented by other indicators allowing to characterize the typical shape of the distribution, in particular the bounds of the distribution range, the standard deviation and the three quartiles.

Estimation

The estimation is performed using the mhurdle function, which has the following arguments:

-formula: a formula describing the model to estimate. It should have between two and four parts on the right-hand side specifying, in the first part, the good selection equation covariates, in the second part, the desired consumption equation covariates, in the third part, the purchasing equation covariates and in the fourth part, the covariates of the variance equation.

-data: a data frame containing the observations of the variables present in the formula.

-subset, weights, na.action: these are arguments passed on to the model.frame function in order to extract the data suitable for the model. These arguments are present in the lm function and in most of the estimation functions.

-start: the starting values, that can be provided by the user, which may be necessary for complicated and highly parametric models, -dist: this argument indicates the functional form of the desired consumption equation, which may be either log-normal "ln" (the default), normal "n", Box-Cox "bc", or inverse Hyperbolic Sine "ihs",

-corr: this boolean argument indicates whether the disturbance of the different equations are correlated, the default value is FALSE,

-h2: this boolean argument indicates whether the second hurdle is effective or not,

-robust: if TRUE, transformations of some parameters are used, so that they lie in the required range (positive values for the standard deviation and for the position parameter, between -1 and +1 for the coefficients of correlation),

-...: further arguments that are passed to the optimisation function maxLik.

For sake of clarity, we’ll denote the estimated models using à 5 digits name. The first one is a capital letter indicating the distribution used (either N for a normal distribution and L for a log-normal distribution), the second, third and fourth ones are 1 if the first, second and third hurdles are present, 0 if they are not and the fifth digit is either D or I if a dependent or an independent model is estimated.

To illustrate the use of mhurdle, we start with single equation tobit models, using a normal and a log-normal specification.

N010I <- mhurdle(shows ~ 0 | linc + smsa + age + 
                     educ + size, data = Interview, 
                 h2 = TRUE, dist = "n", method = "bhhh")
L010I <- update(N010I, dist = "ln")

We then consider the two other potential hurdles, namely the selection process and the infrequency of purchase. We only keep linc as a covariate for the desired consumption and the other covariates (smsa, age, educ and size) are used either in the first (selection model) or the third part (infrequency of purchase model) of the formula. Note that the dist argument is not provided (so that the default log-normal distribution is chosen) and that the lack of ressource hurdle is either maintained or removed by setting the h2 argument respectively to TRUE or FALSE

L100D <- mhurdle(shows ~ smsa + age + educ + size | 
                     linc, data = Interview, 
                 h2 = FALSE, dist = "ln", corr = TRUE, method = "bhhh", 
                 finalHessian = TRUE)
L100D2 <- update(L100D, start = coef(L100D), robust = FALSE)
L110D <- update(L100D, h2 = TRUE)
L110D2 <- update(L110D, start = coef(L110D), robust = FALSE)
L001D <- mhurdle(shows ~ 0 | linc | smsa + age +  
                       educ + size, data = Interview, 
                   h2 = FALSE, corr = TRUE, method = "bhhh", 
                 finalHessian = TRUE)
L001D2 <- update(L001D, start = coef(L001D), robust = FALSE)
L011D <- update(L001D, h2 = TRUE)
L011D2 <- update(L011D, start = coef(L011D), robust = FALSE)

Finally, we estimate three equations models. The 4 socio-economic covariates that were previously either included in the first or the third parts of the formula are now split on these two equations. More precisely, we assume that educ and size are the determinants of the selection process and that smsa and size explain the frequency of purchase. The first model is a double-hurdle model ; selection and purchasing mechanism are effective and note that in this case, the h2 argument is set to FALSE). The last two ones triple-hurdle models, with correlated and uncorrelated errors:

L101D <- mhurdle(shows ~ educ + size | linc | 
                        smsa + age, data = Interview, 
                    h2 = FALSE, method = "bhhh", corr = TRUE, 
                 finalHessian = TRUE)
L101D2 <- update(L101D, start = coef(L101D), robust = FALSE)
L111D <- update(L101D, h2 = TRUE)
L111D2 <- update(L111D, start = coef(L111D), robust = FALSE)
L111I <- update(L111D, corr = FALSE)
L111I2 <- update(L111I, start = coef(L111I), robust = FALSE)

The results are presented in table @ref(tab:estimations).

Tests

The position parameter (μ in table @ref(tab:estimations) enables to test directly the hypothesis that the second hurdle is operative, using a Wald-like test. For all the models, the hypothesis that the coefficient is equal to 0 is highly rejected. Note also that the only model which impose that the second hurdle is not operative (denoted L101D in the table) has the worst value of log-likelihood.

The coefficient of correlation is highly significant for the double-hurdle selection model, but not for the double-hurdle selection model. For the most general three equations model, only the coefficient of correlation for the first two equations is significant. Formal tests for this two nested models can be performed using a traditional log-likelihood ratio test:

library("lmtest")
lrtest(L111D, L111I)

or a Vuong test, which don’t make the hypothesis that one of the competing model is the “true” model:

vuongtest(L111D, L111I, type = "nested")

Both tests don’t reject the independence hypothesis at the 5% level, especially the Vuong test, the p-values being respectively equal to 0.076 and 0.2154219.

For non-nested models, one can use Vuong test. For example, for the double-hurdles selection and ptobit models, the test is performed using:

vuongtest(L110D, L011D)

The statistic is -2.173. Its negative sign indicate that the p-tobit model better fits the data compared to the selection model. The difference is significant at the 5% level, as the (two-ways) p-values equal 0.015, which leads to the conclusion that the ptobit model is significantly better than the selection model.

The parcimony principal leads to the selection of the ptobit double-hurdle model, as it is not significantly worst than the more general three-hurdle model.