Next Article in Journal
The Quantum Geometric Tensor in a Parameter-Dependent Curved Space
Next Article in Special Issue
Targeted L1-Regularization and Joint Modeling of Neural Networks for Causal Inference
Previous Article in Journal
A Conditional Mutual Information Estimator for Mixed Data and an Associated Conditional Independence Test
Previous Article in Special Issue
High Resolution Treatment Effects Estimation: Uncovering Effect Heterogeneities with the Modified Causal Forest
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simultaneous Maximum Likelihood Estimation for Piecewise Linear Instrumental Variable Models

1
Department of Statistics, The Pennsylvania State University, University Park, PA 16801, USA
2
Department of Statistics and Actuarial Science, University of Waterloo, Waterloo, ON N2L 3G1, Canada
*
Author to whom correspondence should be addressed.
Entropy 2022, 24(9), 1235; https://doi.org/10.3390/e24091235
Submission received: 26 July 2022 / Revised: 28 August 2022 / Accepted: 31 August 2022 / Published: 2 September 2022
(This article belongs to the Special Issue Causal Inference for Heterogeneous Data and Information Theory)

Abstract

:
Analysis of instrumental variables is an effective approach to dealing with endogenous variables and unmeasured confounding issue in causal inference. We propose using the piecewise linear model to fit the relationship between the continuous instrumental variable and the continuous explanatory variable, as well as the relationship between the continuous explanatory variable and the outcome variable, which generalizes the traditional linear instrumental variable models. The two-stage least square and limited information maximum likelihood methods are used for the simultaneous estimation of the regression coefficients and the threshold parameters. Furthermore, we study the limiting distribution of the estimators in the correctly specified and misspecified models and provide a robust estimation of the variance-covariance matrix. We illustrate the finite sample properties of the estimation in terms of the Monte Carlo biases, standard errors, and coverage probabilities via the simulated data. Our proposed model is applied to an education-salary data, which investigates the causal effect of children’s years of schooling on estimated hourly wage with father’s years of schooling as the instrumental variable.

1. Introduction

In observational studies, the measured confounders can be controlled by a variety of methods such as propensity score based matching and regression adjustment. However, when the confounding variable is unmeasured, the traditional causal inference methods usually lead to biased estimators since changes in the unmeasured confounder will lead to changes in the explanatory variable, both of which will result in changes in the response variable. Failing to adjust such a confounder will lead to spurious association between the explanatory variable and the outcome. Analysis of instrumental variables (IV) has gained popularity in causal inference, such as investigating causal graphical structures [1,2] and controlling for unmeasured confounding [3,4]. An instrument is a variable that is correlated with the explanatory variable but not associated with any unmeasured confounders. In addition, the instrumental variable is supposed to have influence on the response variable only through the explanatory variable, i.e., there is no direct effect of this variable on the response. Instrumental variable analysis can be applied to many areas and disciplines, such as economics and epidemiology. For example, causality between the years of schooling and earnings in economics has been studied in the literature [5]. This example exploits the college proximity as the instrumental variable because it is revealed that those living near college or university usually have significantly higher level of education than others. On the other hand, it is believed that college proximity may improve earnings only by increasing the subject’s years of schooling. Both indicate that college proximity is a useful instrumental variable. In biomedical and epidemiological research, the main interest is to investigate the causal effect of an exposure variable on a certain disease outcome. A gene can be assumed as a good instrument if it is closely linked to the exposure but has no direct effect on the disease [6]. The study of genetic variants as instrumental variables is known as Mendelian randomization, which is discussed extensively in the literature (e.g., [7,8]). For instance, a set of 32 recently identified genetic variants are used as instrumental variables to study whether child fat mass causally affects academic achievement and blood pressure [9].

1.1. Related Work

Since the development of instrumental variables, a plenty of instrumental variable estimation methods have been proposed for the causal effect estimation. Two-stage least squares (2SLS) [10] is one of the most commonly used methods for the instrumental variable estimation. Theoretical analyses such as consistency and asymptotic normality also exist in the literature. When the response variable is binary, the second stage can be modified with logistic regression in mendelian randomization studies [11]. Another method is the likelihood-based method, particularly the limited information maximum likelihood (LIML) [12]. It is proved that the LIML method is more effective in dealing with the weak instruments [13]. The phenomenon of weak instruments occurs when the correlation between the instrument(s) and the explanatory variable is close to zero. When there are weak instruments, 2SLS is generally unstable and the causal effect estimators are badly biased. The typical rule of thumb to detect weak instruments is the F-statistic, which states that an instrument may be weak if the first-stage F-statistic is less than 10 [14].
Most of the IV approaches impose linear assumptions among the instrument, explanatory and response variables. However, this is not always the case. For example, a subject’s years of schooling may only have a positive effect on subsequent earnings if the subject obtained at least a high-school degree. There would be no difference in the earnings if the subject obtained either an elementary or middle school degree. In this hypothetical scenario, a linear regression model between the explanatory and response variables is clearly misspecified. When the null hypothesis of linearity relationship is rejected, one strategy could be to develop piecewise linear models, which is more interpretable compared to the completely nonlinear models.
In this paper, we propose a piecewise linear instrumental variable (PLIV) model for estimating the causal effect via a continuous threshold function. The continuous threshold function assumes that both the explanatory variable and the instrumental variable are continuous. Instrumental variable models with continuous variables have been studied extensively in the literature. For example, continuous instruments have been used in the classical IV models, developed in a structural equation modeling framework [15]. A recent paper proposes semiparametric doubly robust estimators of causal effects with the continuous instruments [16]. Moreover, some discussions about continuous exposure and a continuous response for Mendelian randomization can be found in a review paper [8].
A threshold in a variable occurs when there is a sudden change in the values of this variable. We call the point where the change happens as a cut-off point or a threshold. The subset causal effect exists when there is a threshold in the explanatory variable. The proposed PLIV model is useful because it can study the subset causal effect when the true model is not linear and it can also degenerate to a linear instrumental variable model when the relationship among the variables is indeed linear. In other words, by using piecewise linear functions, we can quantitatively find the subset effects of the explanatory and the instrumental variables.
We use the Rectified Linear Unit (ReLU) function, mathematically defined in Equation (1), to incorporate the piecewise relationships. Utilization of ReLU function for defining the subset effects have been studied in the literature, such as a regression kink model that tests the presence of the threshold [17] and the segmented and hinge models to study the subset effects in logistic regression [18]. Besides, the continuous threshold models via the ReLU function with two-way interactions is considered in the Cox’s proportional hazards model, where the asymptotic normality under mild conditions is established [19]. In this paper, we use a continuous threshold function with multiple thresholds to formulate the piecewise linear instrumental variable models. A similar study of the piecewise linear instrumental variable model through the random slope approach is studied in the literature [20]. It divides the data into a few segments and analyzes the data in each segment individually. However, this method suffers from huge efficiency and accuracy loss.

1.2. Contribution of This Article

In this paper, we consider a piecewise linear model when the linearity assumption of the data is inappropriate and provide a rigorous treatment of the statistical properties of the model. Our contributions can be summarized as follows.
  • We simultaneously estimate the coefficients and thresholds of the piecewise linear instrumental variable model by the limited information maximum likelihood (LIML) method, assuming the number of thresholds is known.
  • The proposed piecewise linear instrumental variable model will degenerate to the linear instrumental variable model if there are no thresholds. Therefore, it provides a generalization to the linear instrumental variable model. To our best knowledge, this is the first work on the piecewise linear extension to the traditional linear instrumental variable models.
  • We also study the theoretical properties of the PLIV model, including the consistency and asymptotic normality of the estimators.

2. Piecewise Linear Instrumental Variable Model

Notations: we denote scalars by unbolded lowercase letters (e.g., sample size n and the i-th observation of outcome variable y i ), random variable by unbolded capital letter (e.g., X), random vectors by boldface lowercase letters (e.g., x i and β ), and matrices with boldface capital letters (e.g., X ).
In the ordinary linear regression model y i = x i β + ϵ i , there is an assumption that the explanatory variables are uncorrelated with the error term, i.e., cov( x i , ϵ i ) = 0. However, there are some situations where the covariance between the explanatory variables and error term exists. This leads to inconsistent estimation of ordinary least squares due to the phenomenon of endogeneity in x . One way to deal with this issue is to introduce an instrument variable, whose changes are related to changes in the explanatory variable but do not lead to the change in the response variable directly.
Let ( x i , y i , z i ) , i = 1 , , n , denotes the observed data for ( X , Y , Z ) , where X is the explanatory variable, Y is the response variable, and Z is the instrumental variable. To estimate the subset causal effect and establish the piecewise linear relationship, for any threshold parameter t R , we use a continuous threshold function which is defined as:
φ ( x i , t ) = ( x i t ) I ( x i > t ) = ( x i t ) + ,
where I ( · ) is an indicator function. ReLU function, commonly used as an activation function in deep learning, is a special case with t = 0 such that φ ( x i , 0 ) = ( x i 0 ) I ( x i > 0 ) = ( x i 0 ) + .
The proposed model provides sparsity and computational efficiency compared to the smoothing or approximation approach in the literature. The estimation stage involves indicator functions but it does not require an approximation of the indicator function. Let K and J denote the number of thresholds in Z and X, respectively. Denote c = ( c 1 , , c K ) T as the vector of thresholds in Z and denote t = ( t 1 , , t J ) T as the vector of thresholds in X. We propose the following piecewise linear instrumental variable model:
x i = α 0 + α 1 φ ( z i , c 1 ) + + α K φ ( z i , c K ) + α K + 1 z i + v i
y i = β 0 + β 1 φ ( x i , t 1 ) + + β J φ ( x i , t J ) + β J + 1 x i + u i ,
where β = ( β 0 , , β J + 1 ) T is the vector of coefficients representing the causal effect of X on Y; α = ( α 0 , , α K + 1 ) T is the vector of coefficients representing the instrumental effect of Z on X; u i and v i are the error terms for the ith observation. In the context of causal inference, we interpret β as the causal effect of x on y. More specifically, for t j < x t j + 1 , 1 j J with t J + 1 denoting the maximum value of x, one unit increase in x leads to β J + 1 + j = 1 j β j units change in y. Besides, β J + 1 represents the change in y that is caused by one unit increase in x for t 0 < x t 1 where t 0 is the minimum value of x. To better understand this, in Figure 1, we plot the function y = φ ( x , 2 ) + 3 × φ ( x , 3 ) + 2 x where β 1 = 1 , β 2 = 3 , β 3 = 2 as an example. When 2 < x 3 , the slope is β 1 + β 3 = 3 . When 3 < x 4 , the slope is β 1 + β 2 + β 3 = 6 .
Here, we assume K and J are prespecified according to some prior knowledge or theoretical justifications. Practically, we may use the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) [21] to select them. A more elegant examination of the condition for the number of thresholds can be found in Newey [22]. In particular, when α 1 = = α K = 0 and β 1 = = β J = 0 , our proposed model degenerates to the traditional linear instrumental variable model.
For instrumental variable analysis, an instrumental variable is correlated with the explanatory variable but not correlated with the error term. In our model, ( Z c ) + = { ( Z c 1 ) + , , ( Z c K ) + } is the vector of instrumental variables with the following properties:
  • Instrument relevance: cov { ( Z c ) + , X } 0 : ( Z c ) + is correlated with the explanatory variable X.
  • Instrument exogeneity: cov { ( Z c ) + , U } = 0 : ( Z c ) + is uncorrelated with the error term U.
We assume K J for identifiability, i.e., the number of instruments should be larger than or equal to the number of endogenous variables.
Remark 1.
Note that intensive research about nonlinear instrumental variable models has been conducted in the literature, such as the nonparametric instrumental regression [23,24,25]. We point out that the target of our method is to quantitatively find the thresholds and estimate the subset causal effects. We aim to generalize the traditional linear IV model and fit an interpretable model rather than approximate the data by a nonlinear function.
To estimate the unknown parameters in (2) and (3), we utilize the two-stage least square (2SLS) method and the limited information maximum likelihood (LIML) method. Details about the proposed estimation methods are discussed below.

3. Simultaneous Maximum Likelihood Estimation

We first introduce how the LIML method is used in our model and initialize the naive estimators by the 2SLS method.

3.1. Limited Information Maximum Likelihood

As discussed in the introduction about the advantages, limited information maximum likelihood is another popular approach for estimation in the instrumental variable models. Here, we assume the error terms ( U , V ) are jointly normally distributed and correlated to some extent due to the unmeasured confounding effect. Let 0 be the zero-mean vector and ρ be the correlation of ( U , V ) . Denote σ u 2 and σ v 2 as the variance of the error terms U and V, respectively. Then the probability density function of the bivariate normal (U,V) can be written as:
f ( U , V ) = 1 2 π σ u σ v 1 ρ 2 exp 1 2 ( 1 ρ 2 ) Q ( U , V ) ,
where the quadratic form Q ( U , V ) = U T U σ u 2 2 ρ U T V σ v σ u + V T V σ v 2 . For a single observation, the log-likelihood is
( u i , v i ; θ ) log σ u σ v 1 2 log ( 1 ρ 2 ) 1 2 ( 1 ρ 2 ) u i 2 σ u 2 2 ρ u i v i σ u σ v + v i 2 σ v 2 ,
where θ = ( α T , β T , c T , t T , ρ , σ u , σ v ) T denote all the model parameters and
v i = x i α 0 α 1 φ ( z i , c 1 ) α K φ ( z i , c K ) α K + 1 z i
u i = y i β 0 β 1 φ ( x i , t 1 ) β J φ ( x i , t J ) β J + 1 x i .
To simplify notations, we let ( θ ) = ( u i , v i ; θ ) denote the log-likelihood. The maximum likelihood estimates for θ is obtained by maximizing the log-likelihood within the compact set Θ R D ( θ ) such that θ ^ n = arg max θ Θ n ( θ ) , where n ( θ ) = 1 / n i = 1 n ( θ ) . However, there is no closed-form solution for θ , so we take the gradient-based algorithm for estimation. This yields approximate M-estimators. To speed up estimation, we use the two-stage least square method to initialize the estimators.

3.2. Initialization: Two-Stage Least Square

The traditional two-stage least squares method regresses the explanatory variable on the instrumental variable and computes the predictions x ^ in the first stage. In the second stage, it regresses the response variable on the predictions x ^ . The causal effect of interest is estimated from the second stage. In our method, we employ 2SLS to obtain the initial values of the parameters of the piecewise linear instrumental variable model. Below we describe the 2SLS procedures for initializations:
Stage 1: First, we regress x on { ( z c ) + , z } and then obtain the fitted values x ^ , where ( z c ) + = { ( z c 1 ) + , , ( z c K ) + } .
Stage 2: We regress y on { ( x ^ t ) + , x ^ } , where ( x ^ t ) + = { ( x ^ t 1 ) + , , ( x ^ t J ) + } . Thus, in the second stage, we fit the following regression model:
y i = β 0 + β 1 φ ( x ^ i , t 1 ) + + β J φ ( x ^ i , t J ) + β J + 1 x ^ i + u i .
For each combination of the number of thresholds in X and Z, we could pick c , t and the regression coefficients simultaneously through grid search when the sum of squared errors (SSE) of Y is minimized. However, for J 2 or K 2 , it is slightly computationally expensive to conduct grid search. Since we only need 2SLS to provide the initialization of the parameters in our method, we choose c to be a vector of the points that are evenly spaced between the 5% to 95% quantiles of Z. Similarly, we choose t to be a vector of the points that are evenly spaced between the 5% to 95% quantiles of X. We ignore points below and above the 5% to 95% quantiles in order to avoid boundary effects. The regression coefficients are obtained accordingly.

3.3. Theoretical Analysis

Under mild conditions, we study the statistical properties of the proposed model and establish the robust variance-covariance estimators for the estimated parameters under the correctly specified and misspecified models, separately. To investigate the theoretical properties, we consider the following regularity conditions:
  • C1. Observations ( X i , Y i , Z i ) , i = 1 , , n are independently and identically distributed on a compact set X Y Z R 1 R 1 R 1 . Furthermore, E ( X 2 ) < , E ( Y 2 ) < , and E ( Z 2 ) < .
  • C2. The explanatory variable X and the instrumental variable Z are continuous in the parameter space, i.e., they have continuous probability density functions f X ( · ) and f Z ( · ) . The density functions are uniformly bounded, that is, there exist constants c ̲ 1 , c ̲ 2 , c ¯ 1 , and c ¯ 2 such that
    c ̲ 1 inf Z Z f Z ( · ) sup Z Z f Z ( · ) c ¯ 1 and c ̲ 2 inf X X f X ( · ) sup X X f X ( · ) c ¯ 2 .
    Furthermore, the true value of the coefficients for the threshold effects satisfy α 0 0 and β 0 0 , where α 0 = ( α 20 , , α ( K 1 ) 0 ) and β 0 = ( β 20 , , β ( J 1 ) 0 ) .
  • C3. ( θ ) is upper-semicontinuous for almost all ( X , Y , Z ) , that is, for every θ ,
    lim sup θ n θ ( X , Y , Z ; θ n ) ( X , Y , Z ; θ ) , a . s .
Remark 2.
Condition C1 is commonly used in regression models. Condition C2 is used for estimating the unknown thresholds and ensures the model is identifiable. The continuity requirements of X and Z are used to estimate the thresholds. Condition C3 is used to establish the consistency and the asymptotic normality of the maximum likelihood estimator.
In terms of estimation, we take the gradient-based method which depends on the first order derivative ˙ ( θ ) = ( θ ) / θ (details can be found in Appendix A) with the initialized estimators by 2SLS. In this paper, we do not approximate the indicator function by the logistic function as some researchers do (e.g., [18,26,27]). The gradient-based algorithm for the ReLU function has shown success in the context of deep learning and machine learning. Compared to the approximation techniques as discussed in Section 1, model estimation with the ReLU function is computationally cheaper since no approximation of the indicator function is required. In fact, as long as Condition C2 is satisfied which requires variables X and Z to be continuous, the gradients composed of the indicator functions converge to a continuous function of the threshold parameters as n , for example,
1 n i = 1 n I ( z i > c k ) P E I ( z i > c k ) = P z i > c k ,
for k = 1 , , K by the law of large numbers. Therefore, the second order derivative of the ReLU function with respect to the thresholds can be derived based on the resulting continuous probability function. More specifically, the second order derivative with respect to c k is simply f Z ( c k ) .
To prove the asymptotic normality, we first need to show the consistency of the proposed estimators.
Theorem 1.
Under conditions C1–C4, assume that Θ is compact and the true parameter vector θ 0 = arg max θ Θ E ( θ ) is unique. Furthermore, for every sufficiently small ball B Θ , sup θ B ( θ ) is measurable with E sup θ B ( θ ) < , then θ ^ n p θ 0 .
Proof. 
The proof follows the Theorem 5.7 of van der Vaart [28]. For completeness, we include it as Theorem A1 in Appendix B. To utilize Theorem 5.7, we need to check the condition that ( θ ^ n ) ( θ 0 ) o P ( 1 ) for some θ 0 Θ 0 . This is true since n ( θ ) is continuous in θ , n ( θ ) converges to ( θ ) uniformly, and θ ^ n (approximately) maximizes n ( θ ) . Thus, all the conditions are satisfied and the result follows. □
Theorem 2.
Under conditions C1–C4, let θ 0 be the true value of θ . Let ˙ ( θ ) be a measurable function with E ˙ ( θ ) ˙ ( θ ) T ( i , j ) < for i , j = 1 , , | θ | , where | θ | denotes the number of elements in θ , then
n θ ^ n θ 0 d N 0 , V θ 0 1 M θ 0 V θ 0 1 ,
where M θ 0 = E ˙ ( θ 0 ) ˙ ( θ 0 ) T and ˙ ( θ 0 ) is the first order derivative of ( θ ) with respect to θ evaluated at θ 0 and V θ 0 is the second order derivative of E { ( θ ) } with respect to θ evaluated at θ 0 (derivations in Appendix A). V θ has the form
V θ = V θ ( 1 ) + V θ ( 2 ) = V θ ( 1 ) + 0 0 V α c ( 2 ) 0 0 0 0 0 0 V β t ( 2 ) 0 0 0 V c c ( 2 ) 0 0 0 0 V t t ( 2 ) 0 0 0 0 0 0 0 0 sym . 0 ,
where 0 denotes a zero vector or a zero matrix and 0 denotes a scalar. Details of V θ ( 1 ) and V θ ( 2 ) are given in the Appendix A.
Proof. 
First, note that ( θ ) is Lipschitz continuous in θ . Moreover, the fact that V θ is continuous in θ admits the Taylor expansion of E X Y Z ( θ ) :
E ( X , Y , Z ) ( θ ) = E ( X , Y , Z ) ( θ 0 ) + 1 2 θ θ 0 V θ 0 θ θ 0 T + o p θ θ 0 2 .
Since θ ^ is the maximum likelihood estimate of θ , 1 n i = 1 n ( θ ^ ) sup θ 1 n i = 1 n ( θ ) o P ( 1 n ) . Plus the result from Theorem 1 that θ ^ n p θ 0 , we conclude from Theorem 5.14 of van der Vaart [28] that:
n θ ^ n θ 0 = V θ 0 1 1 n i = 1 n ˙ i ( θ 0 ) + o P ( 1 ) ,
which implies an asymptotic normal distribution with mean 0 and variance-covariance matrix V θ 0 1 M θ 0 V θ 0 1 . □
For completeness, we include Theorem 5.14 of van der Vaart [28] (2000) as Theorem A2 in Appendix B. When the model is correctly specified, V θ 0 = M θ 0 , the asymptotic variance is the inverse of Fisher information. Matrices V θ 0 and M θ 0 are estimated through the replacement of θ 0 by the MLE θ ^ n . Thus, for the correctly specified model, the variance-covariance matrix is estimated by the inverse of M θ ^ n . For the misspecified model, the variance-covariance matrix is estimated by V θ ^ n 1 M θ ^ n V θ ^ n 1 . Let us define V n as the second derivative of n ( θ ) with respect to θ , then we can decompose V n the same way as V θ into two matrices V n ( 1 ) and V n ( 2 ) . Note that V n is the empirical process of V θ and V n p V θ by the law of large numbers, so we use the estimated probability densities f ^ Z ( c ^ k ) and f ^ X ( t ^ j ) for f Z ( c k ) and f X ( t j ) for k = 1 , , K and j = 1 , , J , respectively.

4. Simulation Studies

In this section, we evaluate the performance of the proposed model using simulated datasets. We consider two scenarios with the same sample size n = 500 . We let error terms U and V be jointly normally distributed with mean 0 and correlation ρ { 0.2 , 0.5 , 0.8 } . Here, we consider a common standard deviation σ u = σ v = 0.3 . Besides, we simulate the instrumental variable Z N ( 0 , 1 ) . The first scenario has one threshold in X and one threshold in z, and it takes the following form:
x i = 1 + 0.5 × ( z i 0.5 ) + + z i + v i y i = 0.2 + ( x i 0 ) + + 0.5 × x i + u i .
The true values of the parameters in PLIV models are α = ( 1 , 0.5 , 1 ) , β = ( 0.2 , 1 , 0.5 ) , c = 0.5 , and t = 0 . The second scenario has two thresholds in x and two thresholds in z, and it takes the following form:
x i = 1 + 0.5 × ( z i + 1 ) + + ( z i 1 ) + + z i + v i y i = 1 + 1.2 × ( x i + 1 ) + + ( x i 2 ) + + 0.5 × x i + u i .
The true parameters are α = ( 1 , 0.5 , 1 , 1 ) , β = ( 1 , 1.2 , 1 , 0.5 ) , c = ( 1 , 1 ) , and t = ( 1 , 2 ) . We show the simulated piecewise linear instrumental variable models for scenario 1 and scenario 2 in Figure 2. We replicate the simulation 1000 times to evaluate the finite sample properties of the proposed model by the PLIV method.
Table 1 summarizes the biases, standard errors of θ ^ and coverage probabilities of θ by the proposed PLIV method for scenario 1, where tse is the theoretical standard error and ese is the empirical standard error. As we can see in the table, all the biases of θ ^ are close to zero. We also find that the theoretical standard error and the empirical standard error are close enough, which confirms the validity of our theoretical results in Section 3. The results show that our model estimation is quite accurate and therefore provides unbiased and consistent estimators. Besides, we notice that the coverage probabilities are around 95% under different values of ρ . Moreover, biases and the standard errors decrease as we increase ρ because the instrumental variables becomes stronger.
Table 2 summarizes the biases, standard errors of θ ^ and 95% coverage probabilities of θ by the PLIV method for scenario 2, where tse is the theoretical standard error and ese is the empirical standard error. We find the similar patterns as in Table 1 from scenario 1. For instance, all the biases are small. Theoretical standard errors and the empirical standard errors are close to each other. Most coverage probabilities are around 95% when ρ = 0.2 and ρ = 0.5 . We also observe that the coverage probabilities of the thresholds are slightly low when ρ = 0.8 . The reason might be due to the high correlation between errors. With multiple thresholds and high correlation, it poses challenges to estimate the exact locations.
We include results with a sample size of 1000 in Appendix C, while fixing ρ = 0.5 . Overall, as n increases, we observe that both biases and standard errors drop.

5. Application

In this section, we revisit the Card’s education data [5]. We apply the proposed model to study the causal effect of years of schooling on hourly wage in cents with father’s years of schooling as the instrumental variable. The interest here is to find a threshold and study the threshold effect of the years of schooling. It is generally believed that a child’s years of schooling has a direct effect on the child’s wage and parents’ education only affects the child’s income by affecting the child’s education level. In other words, parents’ education level has no direct effect on child’s wage. Therefore, the father’s years of schooling can be treated as a valid instrumental variable.
In Card’s data, we remove the missing values and include a total of n = 2657 observations. The explanatory variable X (child’s years of education) is between 1 and 18 with median 13, and the instrumental variable Z (father’s years of education) has minimum 0, maximum 18, and median 12. Figure 3 indicates that variables X and Y are skewed and have heavy tails so transformations are needed before the analysis. A log transformation is applied to both.
Table 3 shows the point estimate, standard error, and associated 95% confidence interval of θ by the proposed model with K = 1 and J = 0 , which are selected by BIC. In the table, α 1 and c are the coefficient and threshold for the transformed father’s years of schooling, respectively. β 1 is the causal effect of years of schooling on earnings. The estimated causal effect of interest β ^ 1 is 0.87, which results in a difference of exp ( 0.87 × a ) units increase in wage if there are a units increase in the log of years of schooling. In economics, β ^ 1 is interpreted as “elasticity". That is, if years of education increases by 1%, the person’s income will increase by 0.87% by our estimation. In terms of the instrumental variable, we notice that the threshold c is estimated to be 7.86. The corresponding p-value is not calculated since testing c = 0 is meaningless in this context. It shows that there exists a threshold at around 8 in the father’s years of schooling. That is, the father’s years of schooling only has a positive effect on the child’s years of schooling if father receives at least 8 years of education. This information can not be observed if the traditional 2SLS method or nonparametric approaches are applied to analyze the data. The threshold effect as well as the thresholds are all statistically significant since their corresponding p-values are far less than 0.05.

6. Discussion, Limitations, and Future Research

In this paper, we propose a simultaneous maximum likelihood estimation for a piecewise linear instrumental variable model. We use the two-stage least square estimators as the initial values and the limited information maximum likelihood methods to estimate the regression coefficients and the threshold parameters simultaneously. We also provide a robust inference of the proposed model. The proposed model with the piecewise linear functions allows us to find the thresholds for both the explanatory and the instrumental variables, which generalizes the traditional linear instrumental variable models. In the simulation study, we evaluate the performance of the proposed model and find that it behaves well in terms of the biases, standard errors, and coverage probabilities in different settings.
In our model, we include a single continuous explanatory variable and a single continuous instrumental variable. We assume the explanatory variable and the instrumental variable are continuous. More complicated cases can be considered. For example, developing a piecewise linear model with count data might be interesting. However, finding the optimal number of thresholds as well as the locations is challenging from the theoretical side. Furthermore, we assume the number of thresholds K and J are prespecified. Treating the numbers of thresholds as random variables, finding the optimal values, and investigating the theoretical properties can be future research.

Author Contributions

Conceptualization, S.S.L. and Y.Z.; methodology, S.S.L. and Y.Z.; experiments and analysis, S.S.L.; original draft writing, S.S.L.; writing review and editing, S.S.L. and Y.Z. All authors have read and agreed to the published version of the manuscript.

Funding

Zhu’s research is supported by the National Sciences and Engineering Research Council of Canada (Grant No. RGPIN-2017-04064).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data used in the application section come from the ivmodel package of CRAN, which can be downloaded from https://github.com/hyunseungkang/ivmodel/tree/master/data (accessed on 31 August 2022). Codes to simulate data, generate tables and plots in Section 4 can be found at https://github.com/shuoshuoliu/PLIV (accessed on 31 August 2022).

Acknowledgments

We thank the editor, the associate editor, and the three reviewers for careful reviews and insightful comments, which have improved this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of the Information and Hessian Matrices

The likelihood to be minimized is
θ = 1 n i = 1 n log σ u σ v 1 2 log ( 1 ρ 2 ) 1 2 ( 1 ρ 2 ) u i 2 σ u 2 2 ρ u i v i σ u σ v + v i 2 σ v 2 .
When the model is specified,
E X Y Z θ = log σ u σ v 1 2 log ( 1 ρ 2 ) 1 2 ( 1 ρ 2 ) E X Y Z U T U σ u 2 2 ρ U T V σ v σ u + V T V σ v 2 .
To write out the first order derivative ˙ ( θ ) of θ with respect to θ , we define the following notations. θ / α c is the row concatenation of the first order derivative of θ with respect to α and c . θ / β t is the row concatenation of the first order derivative of θ with respect to β and t . For notation simplicity, we drop the subscription i. Let α I ( z > c ) = { α 1 I ( z > c 1 ) , , α k I ( z > c K ) } and β I ( x > t ) = { β 1 I ( x > t 1 ) , , β j I ( x > t J ) } . Then we can divide the first order derivative ˙ ( θ ) as following
θ α c = 1 n i = 1 n 1 , ( z c ) + , z , α I ( z > c ) T 1 ( 1 ρ 2 ) ( v σ v 2 ρ u σ u σ v ) θ β t = 1 n i = 1 n 1 , ( x t ) + , x , β I ( x > t ) T 1 ( 1 ρ 2 ) ( u σ u 2 ρ v σ u σ v ) θ ρ = 1 n i = 1 n ρ 1 ρ 2 ρ ( 1 ρ 2 ) 2 u 2 σ u 2 2 ρ u v σ u σ v + v 2 σ v 2 + u v σ v σ u ( 1 ρ 2 ) θ σ u = u 2 ( 1 ρ 2 ) σ u 3 ρ u v ( 1 ρ 2 ) σ v σ u 2 1 σ u θ σ v = v 2 ( 1 ρ 2 ) σ v 3 ρ u v ( 1 ρ 2 ) σ u σ v 2 1 σ v .
The interchangeability of expectation and differentiation is satisfied here and it implies E X Y Z ( θ ) / θ = E X Y Z ˙ ( θ ) . It is easy to check E X Y Z θ / θ = 0 at θ 0 as it should be. We next derive the second order derivative V θ of E X Y Z θ when the model is specified. We partition the symmetric matrix V θ as two symmetric matrices V 1 , θ and V 2 , θ such that
V θ = V θ ( 1 ) + V θ ( 2 ) = V θ ( 1 ) + 0 0 V α c ( 2 ) 0 0 0 0 0 0 V β t ( 2 ) 0 0 0 V c c ( 2 ) 0 0 0 0 V t t ( 2 ) 0 0 0 0 0 0 0 0 sym . 0 .
For the derivation of V θ ( 1 ) , let z c = 1 , ( z c ) + , z and x t = 1 , ( x t ) + , x . Since the matrix V θ ( 1 ) is symmetric, we only need to derive the upper diagonal elements. The first row of V θ ( 1 ) is is the row concatenation of 2 E X Y Z ( θ ) / α 2 , 2 E X Y Z ( θ ) / α β , 2 E X Y Z ( θ ) / α c , 2 E X Y Z ( θ ) / α t , 2 E X Y Z ( θ ) / α ρ , 2 E X Y Z ( θ ) / α σ u , and 2 E X Y Z ( θ ) / α σ v , such that
V 1 , θ ( 1 ) = 1 ( 1 ρ 2 ) E X Y Z ( z c ) T z c σ v 2 , ρ x t σ v σ u , α I ( z > c ) σ v 2 , ρ β I ( x > t ) σ v σ u , 2 ρ v σ v 2 ( 1 ρ 2 ) u ( 1 + ρ 2 ) ( 1 ρ 2 ) σ v σ u , ρ u σ v σ u 2 , ρ u σ v 2 σ u 2 v σ v 3 .
The second row of V θ ( 1 ) is the row concatenation of 2 E X Y Z ( θ ) / β 2 , 2 E X Y Z ( θ ) / β c , 2 E X Y Z ( θ ) / β t , 2 E X Y Z ( θ ) / β ρ , 2 E X Y Z ( θ ) / β σ u , and 2 E X Y Z ( θ ) / β σ v such that
V 2 , θ ( 1 ) = 1 ( 1 ρ 2 ) E X Y Z ( x t ) T x t σ u 2 , ρ α I ( z > c ) σ u σ v , β I ( x > t ) σ v 2 , 2 ρ u σ u 2 ( 1 ρ 2 ) v ( 1 + ρ 2 ) ( 1 ρ 2 ) σ v σ u , ρ v σ v 3 2 u σ u 3 , ρ v σ u σ v 2 .
The third row of V θ ( 1 ) is the row concatenation of 2 E X Y Z ( θ ) / c 2 , 2 E X Y Z ( θ ) / c t , 2 E X Y Z ( θ ) / c ρ , 2 E X Y Z ( θ ) / c σ u , and 2 E X Y Z ( θ ) / c σ v such that
V 3 , θ ( 1 ) = 1 ( 1 ρ 2 ) E X Y Z α I ( z > c ) T α I ( z > c ) σ v σ u , β I ( x > t ) σ u 2 , v ( ρ 2 + 1 ) σ u σ v ( 1 ρ 2 ) 2 ρ v σ v 2 ( 1 ρ 2 ) , ρ u σ v σ u 2 , 2 v σ v 3 ρ u σ u σ v 2 .
The fourth row of V θ ( 1 ) is the row concatenation of 2 E X Y Z ( θ ) / t 2 , 2 E X Y Z ( θ ) / t ρ , 2 E X Y Z ( θ ) / t σ u , and 2 E X Y Z ( θ ) / t σ v such that
V 4 , θ ( 1 ) = 1 ( 1 ρ 2 ) E X Y Z β I ( x > t ) T β I ( x > t ) σ v 2 , v ( 1 + ρ 2 ) σ v σ u ( 1 ρ 2 ) 2 ρ u ( 1 ρ 2 ) σ u 2 , 2 u σ v 3 ρ v σ v σ u 2 , ρ v σ u σ v 2 .
The remaining terms in V θ ( 1 ) is given by
2 E X Y Z ( θ ) / ρ 2 = 1 + ρ 2 ( 1 ρ 2 ) 2 4 u v ρ ( ρ 2 + 1 ) σ u σ v ( ρ 2 1 ) 3 + 2 ρ u v σ u σ v ( ρ 2 1 ) 2 ,
2 E X Y Z ( θ ) / ρ σ u = 2 ρ u 2 ( 1 ρ 2 ) 2 σ u 3 2 u v ρ 2 σ u 2 σ v ( ρ 2 1 ) 2 u v σ u 2 σ v ( 1 ρ 2 ) ,
2 E X Y Z ( θ ) / ρ σ v = 2 ρ v 2 ( 1 ρ 2 ) 2 σ v 3 2 u v ρ 2 σ u σ v 2 ( ρ 2 1 ) 2 u v σ u σ v 2 ( 1 ρ 2 ) ,
2 E X Y Z ( θ ) / σ u 2 = 2 ρ u v ( 1 ρ 2 ) σ v σ u 3 3 u 2 σ u 4 ( 1 ρ 2 ) + 1 σ u 2 ,
2 E X Y Z ( θ ) / σ u σ v = ρ u v ( 1 ρ 2 ) σ v 2 σ u 2 ,
2 E X Y Z ( θ ) / σ v 2 = 2 ρ u v ( 1 ρ 2 ) σ u σ v 3 3 v 2 σ v 4 ( 1 ρ 2 ) + 1 σ v 2 .
In terms of the matrix V θ ( 2 ) , we decompose the following elements
V α c ( 2 ) = E X Y Z v σ v 2 ( 1 ρ 2 ) ρ u σ u σ v ( 1 ρ 2 ) × 0 0 I ( z > c 1 ) 0 0 0 I ( z > c K ) 0 0 ,
V β t ( 2 ) = E X Y Z u σ u 2 ( 1 ρ 2 ) ρ v σ u σ v ( 1 ρ 2 ) × 0 0 I ( x > t 1 ) 0 0 0 I ( x > t J ) 0 0 ,
V c c ( 2 ) = E X Y Z 1 σ u σ v ( 1 ρ 2 ) × α 1 f Z ( c 1 ) 0 0 0 α K f Z ( c K ) ,
V t t ( 2 ) = E X Y Z 1 σ v 2 ( ρ 2 1 ) × β 1 f X ( t 1 ) 0 0 0 β J f X ( t J ) .
It is easy to check that when the model is correctly specified, V θ ( 2 ) = 0 and V θ = E X Y Z ˙ ( θ ) ˙ ( θ ) T .

Appendix B. Theorems

Define P f as the expectation E f ( X ) = f d P and abbreviate the average n 1 i = 1 n f ( X i ) to P n f , an empirical distribution. Furthermore, we define
M n ( θ ) = 1 / n i = 1 n m θ ( X i ) = P n m θ and Ψ n ( θ ) = 1 / n i = 1 n ψ θ ( X i ) = P n ψ θ .
Theorem A1
(Theorem 5.7 of van der Vaart [28]). Let M n be random functions and let M be a fixed function of θ such that for every ϵ > 0
sup θ Θ M n ( θ ) M ( θ ) P 0 , sup θ : d θ , θ 0 ε M ( θ ) < M θ 0 .
Then every sequence of estimators θ ^ n with M n ( θ ^ n ) M n ( θ 0 ) o P ( 1 ) converges in probability to θ 0 .
Theorem A2
(Theorem 5.14 of van der Vaart [28]). For each θ in an open subset of Euclidean space, let θ ψ θ ( x ) be twice continuously differentiable for every x. Suppose that P ψ θ 0 = 0 , that P ψ θ 0 2 < and that the matrix P ψ ˙ θ 0 exists and is nonsingular. Assume that the second-order partial derivatives are dominated by a fixed integrable function ψ ¨ ( x ) for every θ in a neighborhood of θ 0 . Then every consistent estimator sequence θ ^ n such that Ψ n ( θ ^ n ) = 0 for every n satisfies
n θ ^ n θ 0 = P ψ ˙ θ 0 1 1 n i = 1 n ψ θ 0 X i + o P ( 1 ) .
In particular, the sequence n θ ^ n θ 0 is asymptotically normal with mean zero and covariance matrix P ψ ˙ θ 0 1 P ψ θ 0 ψ θ 0 T P ψ ˙ θ 0 1 .

Appendix C. Additional Simulation Results

Table A1. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 1 with sample size 1000.
Table A1. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 1 with sample size 1000.
ρ = 0.2 ρ = 0.5 ρ = 0.8
bias tse ese cp bias tse ese cp bias tse ese cp
α 0 −8.3027.0429.78928−6.4825.4127.35933−3.2822.0022.78942
α 1 3.0868.2970.769502.9664.9967.549322.4653.8855.17949
α 2 −7.9030.9232.50936−6.0528.9030.11938−2.7923.0523.46955
β 0 −3.6138.7639.70949−2.8036.6637.62938−1.7730.7431.19945
β 1 −0.4655.4454.45956−0.1852.0051.809480.6541.7942.11939
β 2 −1.2124.1824.78938−0.8822.7623.35928−0.4318.3818.42949
c−41.08123.92167.07873−31.07111.23148.14873−12.7079.4798.09886
t−7.6368.1876.36919−4.9061.1366.47920−1.5043.6946.31935
ρ 1.0834.2334.639481.0826.5326.779480.7212.4112.49946
σ 2 −0.869.829.68949−0.6410.9610.75949−0.2612.6812.42946
Note: all numbers are multiplied by 1000. These results are based on 1000 replications.
Table A2. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 2 with sample size 1000.
Table A2. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 2 with sample size 1000.
ρ = 0.2 ρ = 0.5 ρ = 0.8
bias tse ese cp bias tse ese cp bias tse ese cp
α 0 −25.84176.86168.03943−15.74155.65161.09929−7.23104.42114.68927
α 1 8.53115.79106.559567.00103.98101.289474.0173.8275.63944
α 2 8.49112.53105.089645.55108.98107.919583.3198.3893.69957
α 3 −11.29108.1499.84951−5.5195.9895.87934−1.8467.3771.38935
β 0 −2.8783.3184.73942−2.0377.0478.41945−1.0959.8462.78929
β 1 3.8649.6950.239452.7246.3246.969422.3436.2737.33941
β 2 5.7772.6467.829603.8867.5663.169632.3253.6552.82939
β 3 −0.6939.9240.14940-0.4837.1437.55943−0.1129.1431.00944
c 1 −16.09171.89185.999230.96152.10212.499072.67103.26158.39891
c 2 −2.6981.7695.519124.3775.10125.849037.8656.60131.88894
t 1 2.1853.2857.749332.0847.8252.289212.3334.1738.55921
t 2 20.13111.57136.4692513.6199.30108.4593013.3070.2285.85927
ρ 1.2132.9633.189531.5225.6325.739500.9312.1412.31942
σ 2 −1.419.819.64948−1.1710.8810.59951−0.5812.5712.27948
Note: all numbers are multiplied by 1000. These results are based on 1000 replications.

References

  1. Sokolovska, N.; Wuillemin, P.H. The Role of Instrumental Variables in Causal Inference Based on Independence of Cause and Mechanism. Entropy 2021, 23, 928. [Google Scholar] [CrossRef]
  2. Zander, B.; Liśkiewicz, M. On searching for generalized instrumental variables. In Proceedings of the Artificial Intelligence and Statistics (PMLR), Cadiz, Spain, 9–11 May 2016; pp. 1214–1222. [Google Scholar]
  3. Angrist, J.D.; Imbens, G.W.; Rubin, D.B. Identification of causal effects using instrumental variables. J. Am. Stat. Assoc. 1996, 91, 444–455. [Google Scholar] [CrossRef]
  4. Greenland, S. An introduction to instrumental variables for epidemiologists. Int. J. Epidemiol. 2000, 29, 722–729. [Google Scholar] [CrossRef] [PubMed]
  5. Card, D. Using Geographic Variation in College Proximity to Estimate the Return to Schooling; Technical Report; National Bureau of Economic Research: Cambridge, MA, USA, 1993. [Google Scholar]
  6. Didelez, V.; Sheehan, N. Mendelian randomization as an instrumental variable approach to causal inference. Stat. Methods Med. Res. 2007, 16, 309–330. [Google Scholar] [CrossRef] [PubMed]
  7. Lawlor, D.A.; Harbord, R.M.; Sterne, J.A.; Timpson, N.; Davey Smith, G. Mendelian randomization: Using genes as instruments for making causal inferences in epidemiology. Stat. Med. 2008, 27, 1133–1163. [Google Scholar] [CrossRef] [PubMed]
  8. Burgess, S.; Small, D.S.; Thompson, S.G. A review of instrumental variable estimators for Mendelian randomization. Stat. Methods Med. Res. 2017, 26, 2333–2355. [Google Scholar] [CrossRef]
  9. von Hinke, S.; Smith, G.D.; Lawlor, D.A.; Propper, C.; Windmeijer, F. Genetic markers as instrumental variables. J. Health Econ. 2016, 45, 131–148. [Google Scholar] [CrossRef] [PubMed]
  10. Theil, H. Economic Forecasts and Policy, 2nd ed.; Palgrave Macmillan: Amsterdam, The Netherlands, 1961. [Google Scholar]
  11. Palmer, T.M.; Holmes, M.V.; Keating, B.J.; Sheehan, N.A. Correcting the Standard Errors of 2-Stage Residual Inclusion Estimators for Mendelian Randomization Studies. Am. J. Epidemiol. 2017, 186, 1104–1114. [Google Scholar] [CrossRef]
  12. Davidson, R. Estimation and Inference in Econometrics; Oxford University Press: New York, NY, USA, 1993. [Google Scholar]
  13. Angrist, J.; Pischke, J. Instrumental Variables in Action: Sometimes You get What You Need. Most. Harmless Econom. Empiricist’s Companion 2009, 113–220. [Google Scholar]
  14. Stock, J.; Wright, J.H.; Yogo, M. A Survey of Weak Instruments and Weak Identification in Generalized Method Of Moments. J. Bus. Econ. Stat. 2002, 20, 518–529. [Google Scholar] [CrossRef]
  15. Wooldridge, J.M. Econometric Analysis of Cross Section and Panel Data; MIT Press: Cambridge, MA, USA, 2010. [Google Scholar]
  16. Kennedy, E.H.; Lorch, S.; Small, D.S. Robust causal inference with continuous instruments using the local instrumental variable curve. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 2019, 81, 121–143. [Google Scholar] [CrossRef] [Green Version]
  17. Hansen, B.E. Regression kink with an unknown threshold. J. Bus. Econ. Stat. 2017, 35, 228–240. [Google Scholar] [CrossRef]
  18. Fong, Y.; Di, C.; Huang, Y.; Gilbert, P.B. Model-robust inference for continuous threshold regression models. Biometrics 2017, 73, 452–462. [Google Scholar] [CrossRef] [PubMed]
  19. Liu, S.S.; Chen, B.E. Continuous threshold models with two-way interactions in survival analysis. Can. J. Stat. 2020, 48, 751–772. [Google Scholar] [CrossRef]
  20. Scheines, R.; Cooper, G.; Yoo, C.; Chu, T. Piecewise Linear Instrumental Variable Estimation of Causal Influence. PMLR 2001, 265–271. [Google Scholar]
  21. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978, 6, 461–464. [Google Scholar] [CrossRef]
  22. Newey, W.K. Efficient instrumental variables estimation of nonlinear models. Econom. J. Econom. Soc. 1990, 48, 809–837. [Google Scholar] [CrossRef]
  23. Darolles, S.; Fan, Y.; Florens, J.P.; Renault, E. Nonparametric instrumental regression. Econometrica 2011, 79, 1541–1565. [Google Scholar] [CrossRef]
  24. Florens, J.P.; Johannes, J.; Van Bellegem, S. Identification and estimation by penalization in nonparametric instrumental regression. Econom. Theory 2011, 27, 472–496. [Google Scholar] [CrossRef]
  25. Carroll, R.J.; Ruppert, D.; Crainiceanu, C.M.; Tosteson, T.D.; Karagas, M.R. Nonlinear and nonparametric regression and instrumental variables. J. Am. Stat. Assoc. 2004, 99, 736–750. [Google Scholar] [CrossRef]
  26. Seo, M.H.; Linton, O. A smoothed least squares estimator for threshold regression models. J. Econom. 2007, 141, 704–735. [Google Scholar] [CrossRef]
  27. Lin, H.; Zhou, L.; Peng, H.; Zhou, X.H. Selection and combination of biomarkers using ROC method for disease classification and prediction. Can. J. Stat. 2011, 39, 324–343. [Google Scholar] [CrossRef]
  28. Van der Vaart, A.W. Asymptotic Statistics; Cambridge University Press: Cambridge, UK, 2000; Volume 3. [Google Scholar]
Figure 1. Plot of the function y = φ ( x , 2 ) + 3 × φ ( x , 3 ) + 2 x .
Figure 1. Plot of the function y = φ ( x , 2 ) + 3 × φ ( x , 3 ) + 2 x .
Entropy 24 01235 g001
Figure 2. Piecewise linear instrumental variable models with simulated data for scenario 1 and scenario 2. The upper panel plots the simulated X versus Z, Y versus X for scenario 1, respectively. The lower panel plots the simulated X versus Z, Y versus X for scenario 2, respectively.
Figure 2. Piecewise linear instrumental variable models with simulated data for scenario 1 and scenario 2. The upper panel plots the simulated X versus Z, Y versus X for scenario 1, respectively. The lower panel plots the simulated X versus Z, Y versus X for scenario 2, respectively.
Entropy 24 01235 g002
Figure 3. Histogram plots of the raw data X, Y, and Z.
Figure 3. Histogram plots of the raw data X, Y, and Z.
Entropy 24 01235 g003
Table 1. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 1.
Table 1. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 1.
ρ = 0.2 ρ = 0.5 ρ = 0.8
bias tse ese cp bias tse ese cp bias tse ese cp
α 0 −19.2541.2545.80937−16.4338.2641.56939−9.1032.0833.78940
α 1 7.6598.27102.669276.3693.1397.029244.1077.3281.80919
α 2 −16.9546.2047.71931−14.7942.8243.64933−8.2833.5234.34943
β 0 −7.8655.4154.87950−6.8852.3752.74944−4.2843.9244.80945
β 1 0.4880.5877.07955−0.3575.4874.69942−0.5860.3762.50940
β 2 −4.3534.5734.06947−3.8432.4932.60945−2.3826.2126.57933
c−95.15178.21247.82839−82.89159.34224.83846−46.25113.96165.49864
t−14.8897.77108.77922−12.7187.80101.10908−6.7662.6971.68908
ρ 2.8248.9947.549512.6737.9136.819471.6217.7017.22941
σ 2 −2.3214.0013.72954−1.8515.6515.40953−1.1018.1217.82956
Note: all numbers are multiplied by 1000. These results are based on 1000 replications.
Table 2. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 2.
Table 2. Empirical biases, theoretical standard errors (tse), and empirical standard errors (ese) of θ ^ , as well as 95% coverage probabilities (cp) on θ for scenario 2.
ρ = 0.2 ρ = 0.5 ρ = 0.8
bias tse ese cp bias tse ese cp bias tse ese cp
α 0 −51.88268.22247.08946−38.92232.37226.53939−20.83158.06169.46921
α 1 29.20176.58157.4696624.67157.87143.2696513.44110.56107.65949
α 2 15.11172.47166.4094311.80178.03163.6394911.40146.19143.76955
α 3 −26.32164.95147.35945−19.39144.98135.53931−9.21101.13101.32934
β 0 −8.36120.42116.63944−8.23111.05108.00950−0.8485.3182.56958
β 1 6.6171.8271.499476.5766.8466.579483.3952.0752.12950
β 2 6.44115.1399.079665.38106.2990.789693.3083.0575.06962
β 3 −4.1457.8956.20947−4.3353.6952.40950−1.1041.8040.31955
c 1 −3.01253.38246.839309.41221.21257.369246.90152.06218.68898
c 2 2.15120.17138.809135.07139.96140.179019.1084.42134.44880
t 1 0.7976.2579.609441.0468.3172.989394.5748.7049.52935
t 2 18.65168.54189.8192617.60149.74174.5491116.26104.90158.56922
ρ 2.8747.4445.589503.4036.8135.359532.1417.3716.77948
σ 2 −3.6414.0013.64939−2.9915.5515.21946−1.8417.9917.63955
Note: all numbers are multiplied by 1000. These results are based on 1000 replications.
Table 3. Summary table of θ by the SML-PLIV model.
Table 3. Summary table of θ by the SML-PLIV model.
ParameterEstimateStd. Errorz Value95% C.I.p-Value
α 0 : intercept2.250.013168.8(2.222, 2.274)≈0
α 1 : ( Z c ) + −0.020.003−4.8(−0.023, −0.009)≈0
α 2 : Z0.040.00314.3(0.033, 0.043)≈0
β 0 : intercept4.040.21718.6(3.613, 4.464)≈0
β 1 : log X 0.870.08410.4(0.705, 1.033)≈0
c7.860.9398.4(6.016, 9.696)-
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, S.S.; Zhu, Y. Simultaneous Maximum Likelihood Estimation for Piecewise Linear Instrumental Variable Models. Entropy 2022, 24, 1235. https://doi.org/10.3390/e24091235

AMA Style

Liu SS, Zhu Y. Simultaneous Maximum Likelihood Estimation for Piecewise Linear Instrumental Variable Models. Entropy. 2022; 24(9):1235. https://doi.org/10.3390/e24091235

Chicago/Turabian Style

Liu, Shuo Shuo, and Yeying Zhu. 2022. "Simultaneous Maximum Likelihood Estimation for Piecewise Linear Instrumental Variable Models" Entropy 24, no. 9: 1235. https://doi.org/10.3390/e24091235

APA Style

Liu, S. S., & Zhu, Y. (2022). Simultaneous Maximum Likelihood Estimation for Piecewise Linear Instrumental Variable Models. Entropy, 24(9), 1235. https://doi.org/10.3390/e24091235

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop