Diagnostics for 2SLS Regression

Introduction

The ivreg package extends a variety of standard numeric and graphical regression diagnostics to linear models fit by two-stage least-squares (2SLS) regression, a commonly employed method of instrumental-variables estimation for potentially overidentified structural equations in which there are endogenous regressors (see e.g., Greene 2003). The ivreg.fit() function in the package computes the 2SLS estimator employing a low-level interface not generally intended for direct use, and returns a list containing quantities that facilitate the computation of various diagnostics. The ivreg() function provides a user-friendly formula-based interface to ivreg.fit().

ivreg() is derived from and supersedes the ivreg() function in the AER package (Kleiber and Zeileis 2008), making additional provision for regression diagnostics. The principal subject of this vignette is the rationale for the extension of various standard regression diagnostics to 2SLS and the use of functions in the ivreg package to compute them, along with functions in other packages, specifically the base-R stats package (R Core Team 2020) and the car and effects packages (Fox and Weisberg 2019), that work with the "ivreg" objects produced by ivreg().

Review of 2SLS Estimation

We’ll need some basic results for 2SLS regression to develop diagnostics and so we review the method briefly here. 2SLS regression was invented independently in the 1950s by Basmann (1957) and Theil (as cited in Theil 1971), who took slightly different but equivalent approaches, both described below, to derive the 2SLS estimator.

We want to estimate the linear model \(y = X \beta + \varepsilon\), where \(y\) is an \(n \times 1\) vector of observations on a response variable, \(X\) is an \(n \times p\) matrix of regressors, typically with an initial columns of \(1\)s for the regression constant, \(\beta\) is a \(p \times 1\) vector of regression coefficients to be estimated from the data, and \(\varepsilon\) is an \(n \times 1\) vector of errors assumed to be distributed \(N_n(0, \sigma^2 I_n)\) where \(N_n\) is the multivariate-normal distribution, \(0\) is an \(n \times 1\) vector of zeroes, and \(I_n\) is the order-\(n\) identity matrix. Suppose that some (perhaps all) of the regressors in \(X\) are endogenous, in the sense that they are thought not to be independent of \(\varepsilon\). As a consequence, the ordinary least-squares (OLS) estimator \(b_{\mathrm{OLS}} = (X^\top X)^{-1} X^\top y\) of \(\beta\) is generally biased and inconsistent.

Now suppose that we have another set of \(q\) instrumental variables (IVs) \(Z\) that are independent of \(\epsilon\), where \(q \ge p\). If \(q = p\) we can apply the IVs directly to estimate \(\beta\), but if \(q > p\) we have more IVs than we need. Simply discarding IVs would be inefficient, and 2SLS regression is a procedure for reducing the number of IVs to \(p\) by combining them in a sensible way.

The first stage of 2SLS regresses all of regressors in the model matrix \(X\) on the IVs \(Z\) by multivariate ordinary least squares, obtaining the \(q \times p\) matrix of regression coefficients \(B = (Z^\top Z)^{-1} Z^\top X\), and the fitted values \(\widehat{X} = Z B\). The columns of \(B\) are equivalent to the coefficients produced by separate least-squares regressions of each of the columns of \(X\) on \(Z\). If some of the columns of \(X\) are exogenous, then these columns also appear in \(Z\), and consequently the columns of \(\widehat{X}\) pertaining to exogenous regressors simply reproduce the corresponding columns of \(X\).

Because the columns of \(\widehat{X}\) are linear combinations of the columns of \(Z\), they are (asymptotically) uncorrelated with \(\varepsilon\), making them suitable IVs for estimating the regression equation. This IV step is the second stage of 2SLS in Theil’s approach.

As an alternative, we can obtain exactly the same estimates \(b_{\mathrm{2SLS}}\) of \(\beta\) by performing an OLS regression of \(y\) on \(\widehat{X}\), producing \(b_{\mathrm{2SLS}} = (\widehat{X}^\top \widehat{X})^{-1} \widehat{X}^\top y\). This is Basmann’s approach and it motivates the name “2SLS.”

Whether we think of the second stage as IV estimation or OLS regression, we can combine the two stages into a single formula: \[ b_{\mathrm{2SLS}} = [X^\top Z(Z^\top Z)^{-1} Z^\top X]^{-1} X^\top Z (Z^\top Z)^{-1} Z^\top y \] This is what the tsls() function in the sem package (Fox, Nie, and Byrnes 2020) does, but from the point of view of developing regression diagnostics, it’s advantageous to compute the 2SLS estimator by two distinct OLS regressions, which is the approach taken by ivreg().

Unusual-Data Diagnostics for 2SLS Regression

As far as we can tell, diagnostics for regression models fit by 2SLS are a relatively neglected topic, but were addressed briefly by Belsley, Kuh, and Welsch (1980, 266–68). Deletion diagnostics directly assess the influence of each case on a fitted regression model by removing the case, refitting the model, and noting how the regression coefficients or other regression outputs, such as the residual standard deviation, change.

Case-deletion diagnostics for influential data can always be obtained by brute-force computation, literally refitting the model with each case removed in turn, but this approach is inefficient and consequently unattractive in large samples. For some classes of statistical models, such as generalized linear models (e.g., Pregibon 1981), computationally less demanding approximations to case-deletion diagnostics are available, and for linear models efficient “updating” formulas are available (as described, e.g., by Belsley, Kuh, and Welsch 1980) that permit the exact computation of case-deletion diagnostics.

As it turns out, and as Belsley, Kuh, and Welsch note, exact updating formulas for 2SLS regression permitting the efficient computation of case-deletion statistics were given by Phillips (1977, Equations 15 and 16). Phillips’s formulas, reproduced here in our notation (and fixing a couple of small typos in the original), are used in the case-deletion statistics computed in the ivreg package: \[ b_{\mathrm{2SLS}-i} = b_{\mathrm{2SLS}} +A^{-1}g_i \] where \(b_{\mathrm{2SLS}-i}\) is the 2SLS vector of regression coefficients with the \(i\)-th case removed, and

\[\begin{align} A &= X^\top Z(Z^\top Z)^{-1}Z^\top X \\ g_i &= u_i[(y_i - z_i^\top a) - (y_i - r_i^\top b_{\mathrm{2SLS}})] + (u_i + j_i)(y_i - x_i^\top b_{\mathrm{2SLS}}) \\ u_i &= \frac{1 - x_i^\top A^{-1} x_i}{[1 - c_i + (x_i - r_i)^\top A^{-1} (x_i - r_i)]D_i}(x_i - r_i) + \frac{(x_i - r_i)^\top A^{-1} x_i}{[1 - c_i + (x_i - r_i)^\top A^{-1} (x_i - r_i)]D_i}x_i \\ j_i &= \frac{(x_i - r_i)^\top A^{-1} x_i}{[1 - c_i + (x_i - r_i)^\top A^{-1} (x_i - r_i)]D_i}(x_i - r_i) - \frac{1}{D_i}x_i \\ r_i &= X^\top Z (Z^\top Z)^{-1} z_i \\ c_i &= z_i^\top (Z^\top Z)^{-1} z_i \\ D_i &= 1 - x_i^\top A^{-1} x_i + \frac{[(x_i - r_i)^\top A^{-1} x_i]^2}{1 - c_i + (x_i - r_i)^\top A^{-1}(x_i - r_i)} \\ a &= (Z^\top Z)^{-1} Z^\top y \end{align}\]

Here, \(y_i\) is the value of the response for the \(i\)-th case, \(x_i^\top\) is the \(i\)-th row of the model matrix \(X\), and \(z_i^\top\) is the \(i\)-th row of the instrumental-variables model matrix \(Z\).

Belsley, Kuh, and Welsch specifically examine (in our notation) the values of \(\mathrm{dfbeta}_i = b_{\mathrm{2SLS}} - b_{\mathrm{2SLS}-i}\). They discuss as well the deleted values of the residual standard deviation \(s_{-i}\). (Belsley, Kuh, and Welsch define the residual variances \(s^2\) and \(s_{-i}^2\) respectively as the full-sample and deleted residual sums of squares divided by \(n\); in the ivreg packages, we divide by the residual degrees of freedom, \(n - p\) for the full-sample value of \(s^2\) and \(n - p - 1\) for the case-deleted values.)

Belsley, Kuh, and Welsch then compute their summary measure of influence on the fitted values (and regression coefficients) \(\mathrm{dffits}\) as \[ \mathrm{dffits}_i = \frac{x_i^\top \mathrm{dfbeta_{i}}}{s_{-i} \sqrt{x_i^\top (\widehat{X}^\top \widehat{X})^{-1} x_i}} \] where (as before) \(x_i^\top\) is the \(i\)-th row of the model matrix \(X\) and \(\widehat{X}\) is the model matrix of second-stage regressors.

Let \[ H^* = X[X^\top Z(Z^\top Z)^{-1} Z^\top X]^{-1} X^\top Z (Z^\top Z)^{-1} Z^\top \] represent the \(n \times n\) matrix that transforms \(y\) into the fitted values, \(\widehat{y} = H^* y\). In OLS regression, the analogous quantity is the hat-matrix \(H = X(X^\top X)^{-1}X^\top\). Belsley, Kuh, and Welsch note that \(H^*\), unlike \(H\), is not an orthogonal-projection matrix, projecting \(y\) orthogonally onto the subspace spanned by the columns of \(X\). (They say that \(H^*\) isn’t a projection matrix, but that isn’t true: It represents an oblique projection of \(y\) onto the subspace spanned by the columns of \(X\).) In particular, although \(H^*\), like \(H\), is idempotent (\(H^* = H^* H^*\)) and \(\mathrm{trace}(H^*) = p\), \(H^*\), unlike \(H\), is asymmetric, and thus its diagonal elements can’t be treated as summary measures of leverage, that is, as hatvalues.

Belsley, Kuh, and Welsch recommend simply using the havalues from the second-stage regression. These are the diagonal entries \(h_i = h_{ii}\) of \(H_2 = \widehat{X}(\widehat{X}^\top \widehat{X})^{-1} \widehat{X}^\top\). We discuss some alternatives below.

In addition to hatvalues, \(\mathrm{dfbeta}\), \(s_{-i}\), and \(\mathrm{dffits}\), the ivreg packages calculates Cook’s distances \(D_i\), which are essentially a slightly differently scaled version of \(\mathrm{dffits}\) that uses the overall residual standard deviation \(s\) in place of the deleted standard deviations \(s_{-i}\): \[ D_i = \frac{s_{-i}^2}{s^2} \times \frac{\mathrm{dffits}_i^2}{p} \]

Because they have equal variances and are approximately \(t\)-distributed under the normal linear model, studentized residuals are useful for detecting outliers and for addressing the assumption of normally distributed errors. The ivreg package defines studentized residuals in analogy to OLS regression as \[ \mathrm{rstudent}_i = \frac{e_i}{s_{-i} \sqrt{1 - h_i}} \] where \(e_i = y_i - x_i^\top b_{2SLS}\) is the response residual for the \(i\)-th case.

As mentioned, Belsley, Kuh, and Welsch (1980) recommend using hatvalues from the second-stage regression. That’s a reasonable choice and the default in the ivreg package, but it risks missing cases that have high leverage in the first-stage but not the second-stage regression. Let \(h_i^{(1)}\) represent the hatvalues from the first stage and \(h_i^{(2)}\) those from the second stage. If the model includes an intercept, both sets of hatvalues are bounded by \(1/n\) and \(1\), but the average hatvalue in the first stage is \(q/n\) while the average in the second stage is \(p/n\). To make the hatvalues from the two stages comparable, we divide each by its average, \(h^{(1*)}_i = \frac{h_i^{(1)}}{q/n}\) and \(h^{(2*)}_i = \frac{h_i^{(2)}}{p/n}\). Then we can define the two-stages hatvalue either as the (rescaled) larger of the two for each case, \(h_i = (p/n) \times \max \left( h^{(1*)}_i, h^{(2*)}_i \right)\), or as their (rescaled) geometric mean, \(h_i = (p/n) \times \sqrt{h^{(1*)}_i \times h^{(2*)}_i}\). The ivreg package provides both of these options.

library("ivreg")
data("Kmenta", package = "ivreg")
Kmenta

##            Q       P     D     F  A
## 1922  98.485 100.323  87.4  98.0  1
## 1923  99.187 104.264  97.6  99.1  2
## 1924 102.163 103.435  96.7  99.1  3
## 1925 101.504 104.506  98.2  98.1  4
## 1926 104.240  98.001  99.8 110.8  5
## 1927 103.243  99.456 100.5 108.2  6
## 1928 103.993 101.066 103.2 105.6  7
## 1929  99.900 104.763 107.8 109.8  8
## 1930 100.350  96.446  96.6 108.7  9
## 1931 102.820  91.228  88.9 100.6 10
## 1932  95.435  93.085  75.1  81.0 11
## 1933  92.424  98.801  76.9  68.6 12
## 1934  94.535 102.908  84.6  70.9 13
## 1935  98.757  98.756  90.6  81.4 14
## 1936 105.797  95.119 103.1 102.3 15
## 1937 100.225  98.451 105.1 105.0 16
## 1938 103.522  86.498  96.4 110.5 17
## 1939  99.929 104.016 104.4  92.5 18
## 1940 105.223 105.769 110.7  89.3 19
## 1941 106.232 113.490 127.1  93.0 20

deq <- ivreg(Q ~ P + D | D + F + A, data=Kmenta)     # demand equation
summary(deq)

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4305 -1.2432 -0.1895  1.5762  2.4920 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 94.63330    7.92084  11.947 1.08e-09 ***
## P           -0.24356    0.09648  -2.524   0.0218 *  
## D            0.31399    0.04694   6.689 3.81e-06 ***
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16    88.025 2.32e-09 ***
## Wu-Hausman         1  16    11.422  0.00382 ** 
## Sargan             1  NA     2.983  0.08414 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.966 on 17 degrees of freedom
## Multiple R-Squared: 0.7548,  Adjusted R-squared: 0.726 
## Wald test: 23.81 on 2 and 17 DF,  p-value: 1.178e-05

seq <- ivreg(Q ~ P + F + A | D + F + A, data=Kmenta) # supply equation
summary(seq)

## 
## Call:
## ivreg(formula = Q ~ P + F + A | D + F + A, data = Kmenta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8724 -1.2593  0.6415  1.4745  3.4865 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.53244   12.01053   4.124 0.000795 ***
## P            0.24008    0.09993   2.402 0.028785 *  
## F            0.25561    0.04725   5.410 5.79e-05 ***
## A            0.25292    0.09966   2.538 0.021929 *  
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   1  16    256.34 2.86e-11 ***
## Wu-Hausman         1  15     36.14 2.38e-05 ***
## Sargan             0  NA        NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.458 on 16 degrees of freedom
## Multiple R-Squared: 0.6396,  Adjusted R-squared: 0.572 
## Wald test:  10.7 on 3 and 16 DF,  p-value: 0.0004196

par(mfrow=c(2, 2))
plot(deq)

library("car") # for diagnostic generic functions

## Loading required package: carData

qqPlot(deq)

## 1937 1929 
##   16    8

influencePlot(deq)

##      StudRes        Hat      CookD
## 1 -1.7359357 0.09079703 0.06956671
## 2 -1.3686682 0.26453459 0.21973049
## 3 -2.0995532 0.13849570 0.17147564
## 4 -0.2010944 0.39711512 0.01508349
## 5 -0.4505155 0.46498004 0.05257374

Kmenta1 <- Kmenta
Kmenta1[20, "Q"] <- 95

deq1 <- update(deq, data=Kmenta1)
compareCoefs(deq, deq1)

## Calls:
## 1: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta)
## 2: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1)
## 
##             Model 1 Model 2
## (Intercept)   94.63  117.96
## SE             7.92   11.64
##                            
## P           -0.2436 -0.4054
## SE           0.0965  0.1417
##                            
## D            0.3140  0.2351
## SE           0.0469  0.0690
## 

qqPlot(deq1)

## 1941 1940 
##   20   19

outlierTest(deq1)

##       rstudent unadjusted p-value Bonferroni p
## 1941 -4.599583         0.00029602    0.0059204

influencePlot(deq1)

##      StudRes       Hat     CookD
## 1 -1.4737565 0.2645346 0.2447875
## 2 -0.9139638 0.3971151 0.2269833
## 3  1.6021281 0.1280020 0.1155278
## 4 -4.5995825 0.4649800 2.8361307

avPlots(deq1)

deq1.20 <- update(deq1, subset = -20)
compareCoefs(deq, deq1, deq1.20)

## Calls:
## 1: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta)
## 2: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1)
## 3: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1, subset = -20)
## 
##             Model 1 Model 2 Model 3
## (Intercept)   94.63  117.96   92.42
## SE             7.92   11.64    9.67
##                                    
## P           -0.2436 -0.4054 -0.2300
## SE           0.0965  0.1417  0.1047
##                                    
## D            0.3140  0.2351  0.3233
## SE           0.0469  0.0690  0.0527
## 

H <- cbind(hatvalues(deq1), hatvalues(deq1, type="both"), 
           hatvalues(deq1, type="maximum"))
colnames(H) <- c("stage2", "geom.mean", "maximum")
head(H)

##          stage2  geom.mean    maximum
## 1922 0.10349313 0.12269459 0.14545857
## 1923 0.11215042 0.12972476 0.15005306
## 1924 0.08882553 0.10233878 0.11790784
## 1925 0.09515432 0.10207539 0.10949987
## 1926 0.06166289 0.07959715 0.10274748
## 1927 0.05684346 0.06794727 0.08122009

scatterplotMatrix(H, smooth=FALSE)

cbind(dfbeta(deq1)[20, ], coef(deq1) - coef(deq1.20))

##                    [,1]        [,2]
## (Intercept) 25.53936742 25.53936742
## P           -0.17547231 -0.17547231
## D           -0.08827334 -0.08827334

c(influence(deq1)$sigma[20], sigma(deq1.20))

##     1941          
## 2.028434 2.028434

Nonlinearity Diagnostics

The theoretical properties of component-plus-residual plots as nonlinearity diagnostics were systematically explored by Cook (1993) and Cook and Croos-Dabrera (1998). Following these authors and focusing on the explanatory variable \(x_1\), let’s assume that the partial relationship of the response \(y\) to \(x_1\) is potentially nonlinear, as represented by the partial regression function \(f(x_1)\), and that the partial relationships of \(y\) to the other \(x\)s are linear, so that an accurate model for the data is: \[ E(y) = \alpha + f(x_1) + \beta_2 x_2 + \cdots + \beta_k x_k \]

We don’t know \(f()\) and so instead fit the working model \[ E(y) = \alpha^\prime + \beta_1^\prime x_1 + \beta_2^\prime x_2 + \cdots + \beta_k^\prime x_k \] in our case by 2SLS regression, obtaining estimated regression coefficients \(a^\prime, b_1^\prime, b_2^\prime, \ldots, b_k^\prime\). Cook and Croos-Dabrera’s work shows that as long as the regression estimator is consistent and the \(x\)s are linearly related, the partial residuals \(b^\prime_1 x_1 + e\) can be plotted and smoothed against \(x_1\) to visualize an estimate of \(f()\), where \(e = y - (a^\prime + b_1^\prime x_1 + b_2^\prime x_2 + \cdots b_k^\prime x_k)\) are the response residuals. In practice, the component-plus-residual plot can break down as an accurate representation of \(f()\) if there are strong nonlinear relationships between \(x_1\) and the other \(x\)s or if \(y\) is nonlinearly related to another \(x\) that is correlated with \(x_1\).

Fox and Weisberg (2018) extend component-plus-residual plots to more complex regression models, which can, for example, include interactions, by adding partial residuals to predictor effect plots. These graphs also can be applied to linear models fit by 2SLS regression.

crPlots(deq, smooth=list(span=1))

ceresPlots(deq, smooth=list(span=1))

library("effects")

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

plot(predictorEffects(deq, residuals=TRUE), 
     partial.residuals=list(span=1))

deq2 <- update(deq, . ~ I((P - 85)^4/10^5) + D)
crPlots(deq2, smooth=list(span=1))

plot(predictorEffects(deq2, residuals=TRUE), 
     partial.residuals=list(span=1))

Nonconstant Error Variance

Standard least-squares nonconstant variance (“heteroscedasticity”) diagnostics extend straightforwardly to 2SLS regression. We can, for example, plot studentized residuals versus fitted values to discern a tendency for the variability of the former to change (typically to increase) with the level of the latter. For the demand equation in Kmenta’s model,

plot(fitted(deq), rstudent(deq))
abline(h=0)

which seems unproblematic.

A variation of this graph, suggested by Fox (2016), adapts Tukey’s spread-level plot (Tukey 1977) to graph the log of the absolute studentized residuals versus the log of the fitted values, assuming that the latter are positive. If a line fit to the plot has slope \(b\), then a variance-stabilizing power transformation is given by \(y^\lambda = y^{1 - b}\). Thus if \(b > 0\), the suggested transformation is down Tukey’s ladder of powers and roots, with, for example, \(\lambda = 1 - b = 1/2\) representing the square-root transformation, \(\lambda = 1 - b = 0\) the log transformation, and so on. For Kmenta’s model, we have

spreadLevelPlot(deq, smooth=list(span=1))

## 
## Suggested power transformation:  -2.44685

which suggests a slight tendency of spread to increase with level. The transformation \(\lambda = - 2.45\) seems strong, until we notice that the values of \(Q\) are far from 0, and that the ratio of the largest to smallest values \(Q_{\mathrm{max}}/Q_{\mathrm{min}} = 106.23/92.42 = 1.15\) is close to 1, so that \(Q^{-2.45}\) is nearly a linear transformation of \(Q\)—that is, effectively no transformation at all:

with(Kmenta, plot(Q, Q^-2.5))
abline(lm(Q^-2.5 ~ Q, data=Kmenta))

A common score test for nonconstant error variance in least-squares regression, suggested by Breusch and Pagan (1979), is based on the model \[ V(\varepsilon) = g(\gamma_0 + \gamma_1 z_1 + \cdots + \gamma_s z_s) \] where the function \(g()\) is unspecified and the variables \(z_1, \ldots, z_s\) are predictors of the error variance. In the most common application, independently proposed by Cook and Weisberg (1983), there is one \(z\), the fitted values \(\widehat{y}\) from the regression, although it is also common to use the regressors \(x\) from the primary regression as \(z\)s. The test is implemented by regressing the squared standardized residuals \(e_i^2/\widehat{\sigma}^2\) on the \(z\)s, where \(\widehat{\sigma}^2 = \sum e_i^2/n\). The regression sum of squares for this auxiliary regression divided by 2 is then asymptotically distributed as \(\chi^2_s\) under the null hypothesis of constant error variance.

The Breusch-Pagan/Cook-Weisberg test is easily adaptable to 2SLS regression, as implemented by the ncvTest() function in the car package. For Kmenta’s demand equation:

ncvTest(deq)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.2390325, Df = 1, p = 0.62491

ncvTest(deq, var = ~ P + D)

## Non-constant Variance Score Test 
## Variance formula: ~ P + D 
## Chisquare = 0.2392964, Df = 2, p = 0.88723

Here, the first test is against the fitted values and the second more general test is against the explanatory variables in the demand equation; the \(p\)-values for both tests are large, suggesting little evidence against the hypothesis of constant variance.

Remedies for nonconstant variance in 2SLS regression are similar to those in least-squares regression:

• We’ve already suggested that if the error variance increases (or decreases) with the level of the response, and if the response is positive, then we might be able to stabilize the error variance by power-transforming the response.

• If, alternatively, we know the variance of the errors up to a constant of proportionality, then we can use inverse-variance weights for the 2SLS estimator. The ivreg() function supports weighted 2SLS regression, and the diagnostics in the ivreg package work with weighted 2SLS fits (see the next section).

• Finally, we can employ a “sandwich” estimator of the coefficient covariance matrix in 2SLS (or the bootstrap: see, e.g., Davison and Hinkley 1997) to correct standard errors for nonconstant error variance, much as in least-squares regression as proposed by P. J. Huber (1967) and White (1980; also see Long and Ervin 2000).

The ivreg package supports the sandwich() function in the sandwich package (Zeileis 2006). For the Kmenta example, where evidence of nonconstant error variance is slight, the sandwich standard errors are similar to, indeed slightly smaller than, the conventional 2SLS standard errors:

summary(deq, vcov=sandwich::sandwich)

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4305 -1.2432 -0.1895  1.5762  2.4920 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 94.63330    5.14745  18.384 1.18e-12 ***
## P           -0.24356    0.07590  -3.209  0.00515 ** 
## D            0.31399    0.04293   7.315 1.21e-06 ***
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16   142.340 6.43e-11 ***
## Wu-Hausman         1  16    21.898 0.000251 ***
## Sargan             1  NA     2.983 0.084137 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.966 on 17 degrees of freedom
## Multiple R-Squared: 0.7548,  Adjusted R-squared: 0.726 
## Wald test: 34.41 on 2 and 17 DF,  p-value: 1.055e-06

SEs <- round(cbind(sqrt(diag(sandwich::sandwich(deq))), 
                   sqrt(diag(vcov(deq)))), 
             4)
colnames(SEs) <- c("sandwich", "conventional")
SEs

##             sandwich conventional
## (Intercept)   5.1475       7.9208
## P             0.0759       0.0965
## D             0.0429       0.0469

We’ll modify Kmenta’s data to reflect nonconstant error variance, regenerating the data as Kmenta did originally from the reduced-form equations, expressing the endogenous variables \(P\) and \(Q\) as functions of the exogenous variables \(D\), \(F\), and \(A\), and reduced-form errors \(\nu_1\) and \(\nu_2\):

Kmenta2 <- Kmenta[, c("D", "F", "A")]
set.seed(492365) # for reproducibility
Kmenta2 <- within(Kmenta2, {
    EQ <- 75.25 + 0.1125*D + 0.1250*F + 0.225*A
    EP <- 85.00 + 0.7500*D - 0.5000*F - 0.900*A
    d1 <- rnorm(20)
    d2 <- rnorm(20)
    v1 <- 2*d1
    v2 <- -0.5*v1 + d2
    w <- 3*(EQ - min(EQ) + 0.1)/(max(EQ) - min(EQ)) 
    v1 <- v1*w # inducing nonconstant variance
    Q <- EQ + v1
    P <- EP + v2
})

Plotting the sampled reduced-form errors v1 against the expectation of Q shows a clear heteroscedastic pattern:

with(Kmenta2, plot(EQ, v1))

Then refitting the demand equation to the new data set, we get

deq2 <- update(deq, data=Kmenta2)
summary(deq2)

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7805 -1.4920  0.1067  1.9745  5.0629 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 111.42774   11.36644   9.803 2.07e-08 ***
## P            -0.39072    0.13979  -2.795  0.01243 *  
## D             0.28415    0.07965   3.567  0.00237 ** 
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16   122.011 2.06e-10 ***
## Wu-Hausman         1  16    22.883 0.000203 ***
## Sargan             1  NA     0.221 0.638221    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.283 on 17 degrees of freedom
## Multiple R-Squared: 0.521,   Adjusted R-squared: 0.4646 
## Wald test: 6.695 on 2 and 17 DF,  p-value: 0.007172

and the nonconstant error variance is clearly reflected in diagnostics; for example,

spreadLevelPlot(deq2)

## 
## Suggested power transformation:  -22.57328

ncvTest(deq2)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 6.690435, Df = 1, p = 0.0096932

The extreme value of the suggested power transformation of \(Q\) from the spread-level plot, \(\lambda =-23\), reflects (as we noted previously) the fact that \(\max(Q)/\min(Q)\) isn’t much larger than 1.

In our example, the sandwich standard errors are not very different from the conventional standard errors:

SEs2 <- round(cbind(sqrt(diag(sandwich::sandwich(deq2))), 
                   sqrt(diag(vcov(deq2)))), 
             4)
colnames(SEs2) <- c("sandwich", "conventional")
SEs2

##             sandwich conventional
## (Intercept)  13.7782      11.3664
## P             0.1702       0.1398
## D             0.0848       0.0797

As mentioned, bootstrapping provides an alternative to sandwich standard errors as a correction for nonconstant error variance, and the ivreg package supplies an "ivreg" method for the Boot() function in the car package, implementing the case-resampling bootstrap, and returning an object of class "boot" suitable for use with functions in the boot package (Davison and Hinkley 1997; Canty and Ripley 2020). By default, \(R = 999\) bootstrap replications are generated. For example:

set.seed <- 869255 # for reproducibility
b.deq2 <- Boot(deq2)
summary(deq2, vcov.=vcov(b.deq2))

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7805 -1.4920  0.1067  1.9745  5.0629 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 111.42774   16.70085   6.672 3.93e-06 ***
## P            -0.39072    0.19678  -1.986  0.06345 .  
## D             0.28415    0.09762   2.911  0.00974 ** 
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16   122.011 2.06e-10 ***
## Wu-Hausman         1  16    22.883 0.000203 ***
## Sargan             1  NA     0.221 0.638221    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.283 on 17 degrees of freedom
## Multiple R-Squared: 0.521,   Adjusted R-squared: 0.4646 
## Wald test:  4.32 on 2 and 17 DF,  p-value: 0.03042

The bootstrap standard errors are larger than the conventional or sandwich standard errors for this example.

Bootstrap confidence intervals can also be computed from the object returned by Boot(), by default reporting \(BC_a\) (bias-corrected, accelerated) intervals (see the documentation for confint.boot() in the car package):

confint(b.deq2)

## Bootstrap bca confidence intervals
## 
##                   2.5 %      97.5 %
## (Intercept) 64.38447188 137.8541642
## P           -0.61242364   0.2463746
## D            0.05739467   0.4253580

Weighted 2SLS Regression

Suppose that we modify the regression model \(y = X \beta + \varepsilon\) so that now \(\varepsilon \sim N_n(0, \sigma^2 W^{-1})\) where \(W = \mathrm{diag}\{w_i\}\) is an \(n \times n\) diagonal matrix of known inverse-variance weights; that is \(V(\varepsilon_i) = \sigma^2/w_i\). As before, some of the columns of \(X\) may be correlated with the errors \(\varepsilon\), but we have sufficient instrumental variables \(Z\) that are independent of the errors.

Then the weighted 2SLS estimator is \[ b_{\mathrm{W2SLS}} = [X^\top W Z(Z^\top W Z)^{-1} Z^\top W X]^{-1} X^\top W Z (Z^\top W Z)^{-1} Z^\top W y \] Alternatively, we can treat the two stages of 2SLS as weighted least squares (WLS) problems, in each stage minimizing the weighted sum of squared residuals. The ivreg() function computes the weighted 2SLS estimator in this manner.

Phillips’s updating formulas for 2SLS regression could also be modified for the weighted case, but a simpler approach (which is evident in the formula for \(b_{\mathrm{W2SLS}}\) above) is to convert the weighted 2SLS problem into an unweighted problem, by transforming the data to constant variance using \(W^{1/2} = \mathrm{diag}\{\sqrt{w_i}\}\), the Cholesky square root of \(W\). The square root of \(W\) is particularly simple because \(W\) is diagonal. Then in Phillips’s updating formulas, we replace \(y\) with \(y^* = W^{1/2}y\), \(X\) with \(X^* = W^{1/2}X\), and \(Z\) with \(Z^* = W^{1/2}Z\).

For the modified Kmenta2 data, we know that the variances of the errors in the demand equation are inversely proportional to the variable w. This is of course artificial knowledge, reflecting the manner in which the data were constructed. The weighted 2SLS estimator is therefore computed as

deqw <- update(deq, data=Kmenta2, weights=1/w)
summary(deqw)

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta2, weights = 1/w)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.43959 -1.66625 -0.08906  1.81440  3.41694 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 107.88374   10.23415  10.542 7.11e-09 ***
## P            -0.33586    0.12240  -2.744   0.0138 *  
## D             0.26347    0.04405   5.981 1.49e-05 ***
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16   101.172 8.31e-10 ***
## Wu-Hausman         1  16    20.105 0.000376 ***
## Sargan             1  NA     0.087 0.767864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.308 on 17 degrees of freedom
## Multiple R-Squared: 0.7166,  Adjusted R-squared: 0.6833 
## Wald test: 18.79 on 2 and 17 DF,  p-value: 4.95e-05

Plotting studentized residuals against fitted values and testing for nonconstant error variance don’t indicate a heteroscedasticity problem, but there is a relatively large studentized residual of about \(-3\) that stands out somewhat from the other values:

ncvTest(deqw)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 4.21029, Df = 1, p = 0.040179

plot(fitted(deqw), rstudent(deqw))

A Bonferroni outlier test suggests that the studentized residual isn’t unusually large, and once more we’re in the unusual situation of knowing that the model is correct.

outlierTest(deqw)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##       rstudent unadjusted p-value Bonferroni p
## 1937 -3.135343          0.0063887      0.12777

Collinearity Diagnostics

In addition to unusual-data diagnostics, Belsley, Kuh, and Welsch (1980) briefly extend their approach to collinearity diagnostics to 2SLS regression. We believe that this approach, which assimilates collinearity to numerical instability, is flawed, in that it takes into account “collinearity with the intercept.” That is, regressors with values far from 0 have large sums of products with the constant regressor, producing a large standard error of the intercept, and simply reflecting the fact that the intercept extrapolates the fitted regression surface far beyond the range of the data.

Fox and Monette (1992) describe an alternative approach to collinearity diagnostics in linear models fit by least squares based on generalized variance-inflation factors. The implementation of generalized variance inflation factors in the vif() function in the car package, which employs the estimated covariance matrix of the coefficients, applies in general to models with linear predictors, including linear models estimated by 2SLS.

For example, for the demand equation in Kmenta’s model:

sqrt(vif(deq))

##        P        D 
## 1.231124 1.231124

Taking the square-roots of the VIFs puts them on the coefficient standard-error scale. That is, the standard errors of the coefficients of \(P\) and \(D\) are 23% larger than they would be if the estimated coefficients were uncorrelated (which is equivalent to the columns of \(\widehat{X}\) for \(P\) and \(D\) in the second-stage regression being uncorrelated). When, as here, each term in the model has just one coefficient, generalized and ordinary variance-inflation factors coincide. The equality of the VIFs for \(P\) and \(D\) is peculiar to the case of two regressors (beyond the regression constant).

Marginal/conditional plots, produced by the mcPlots() function in the car package, superimpose the added-variable plots on corresponding marginal scatterplots for the regressors. These graphs allow us to visualize the reduction in precision of estimation of each coefficient due to collinearity, which reduces the conditional variation of a regressor relative to its marginal variation. for example, for the demand equation:

mcPlots(deq)

The blue points in each panel represent the marginal scatterplot and the magenta points represent the (partial) added-variable plot, with the arrows showing the relationship between the two sets of points.

Robust Two-Stage Instrumental-Variables Regression

The ivreg() function also implements robust two-stage regression with M and MM estimators (P. Huber and Ronchetti 2009) for both stages, via the rlm() function in the MASS package (Venables and Ripley 2002), by specifying, respectively, the argument type = "M" or type = "MM". The M or MM estimator is applied separately to each first-stage regression and to the second-stage regression, thus down-weighting outliers in all of the regressions. The default for the type argument to ivreg(), type = "OLS", produces the traditional 2SLS estimator.

The details of implementation are largely straightforward. Applying the M or MM estimator to the first stage requires a separate regression for each endogenous regressor, rather than a single multivariate regression, because the robustness weights are generally different for each of the first-stage regressions.

As well, rather than basing the estimated error standard deviation \(\widehat{\sigma}\) on the sum of squared model residuals, we use a robust estimate of \(\sigma\) based on the median absolute residual—the median absolute deviation (or MAD) of the residuals around a “center” of 0. The MAD of the residuals is multiplied by 1.4826 to provide a properly scaled estimate of \(\sigma\) when the error distribution is normal: That is, the MAD estimates the 0.75 quantile of the error distribution, and \(1/\Phi^{-1}(0.75) = 1.4825\), where \(\Phi^{-1}(\cdot)\) is the quantile function of the standard-normal distribution.

To illustrate, we apply a two-stage MM regression to the estimation of the demand equation in the corrupted Kmenta data, and compare the result to the traditional 2SLS estimator applied to the uncorrupted data, to the 2SLS estimator applied to the corrupted data, and to the 2SLS estimator deleting the bad case:

deq.mm <- update(deq1, method = "MM")
summary(deq.mm)

## 
## Call:
## ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1, method = "MM")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -13.22380  -1.56688   0.04599   1.12301   2.55786 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 91.09249   10.62357   8.575 1.40e-07 ***
## P           -0.23742    0.11353  -2.091   0.0518 .  
## D            0.34678    0.05688   6.097 1.19e-05 ***
## 
## Diagnostic tests:
##                  df1 df2 statistic  p-value    
## Weak instruments   2  16    88.025 2.32e-09 ***
## Wu-Hausman         1  16     0.803  0.38339    
## Sargan             1  NA    10.330  0.00131 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.08 on 17 degrees of freedom
## Multiple R-Squared: 0.09864, Adjusted R-squared: -0.007398 
## Wald test: 18.59 on 2 and 17 DF,  p-value: 5.262e-05

compareCoefs(deq, deq1, deq1.20, deq.mm)

## Calls:
## 1: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta)
## 2: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1)
## 3: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1, subset = -20)
## 4: ivreg(formula = Q ~ P + D | D + F + A, data = Kmenta1, method = "MM")
## 
##             Model 1 Model 2 Model 3 Model 4
## (Intercept)   94.63  117.96   92.42   91.09
## SE             7.92   11.64    9.67   10.62
##                                            
## P           -0.2436 -0.4054 -0.2300 -0.2374
## SE           0.0965  0.1417  0.1047  0.1135
##                                            
## D            0.3140  0.2351  0.3233  0.3468
## SE           0.0469  0.0690  0.0527  0.0569
## 

The estimates produced by the robust two-stage estimator are very similar to those for the 2SLS estimator applied to the uncorrupted data and when the bad case is deleted from the corrupted data.

It’s instructive to look at the robustness weights for the first- and second-stage MM estimators:

weights(deq.mm, type = "robustness")

##              P   stage_2
## 1922 0.9744449 0.9843442
## 1923 0.9650580 0.9813862
## 1924 0.9952293 0.8700471
## 1925 0.9998998 0.9618035
## 1926 0.9835828 0.8952875
## 1927 0.9996116 0.9752454
## 1928 0.9668182 0.9525911
## 1929 0.6417885 0.5318460
## 1930 0.8028481 0.9085481
## 1931 0.8905474 0.7734043
## 1932 0.9775792 0.9982440
## 1933 0.9985961 0.9050868
## 1934 0.9691516 0.9241239
## 1935 0.8868046 0.9999572
## 1936 0.7238126 0.8787684
## 1937 0.8387680 0.5293154
## 1938 0.9414371 0.9986454
## 1939 0.5325076 0.6927505
## 1940 0.9996695 0.9783155
## 1941 0.9811313 0.0000000

Thus, 1941 is so discrepant in the second stage that it receives 0 weight in the robust fit.

Most of the diagnostic and other methods for "ivreg" objects also work with robust two-stage regression models. To compute deletion statistics like studentized residuals, dfbeta, and Cook’s \(D\), we use the robustness weights for the M or MM second-stage regression (multiplied by prespecified variance weights, if these are employed) as if they were weights in a 2SLS fit. This is, of course, an approximation: To get exact deletion statistics would require literally omitting each case in turn and recalculating the first- and second-stage robust regressions, an unattractive prospect, because the robust regressions are iterative.

Here, for example, is a graph of second-stage hatvalues, studentized residuals, and Cook’s distances for the 2SMM estimator applied to the corrupted Kmenta data:

influencePlot(deq.mm)

##      StudRes        Hat      CookD
## 1  1.6778245 0.13345298 0.07043221
## 2  0.2892163 0.35576176 0.01023249
## 3 -1.3354194 0.28104141 0.13365515
## 4 -2.7083737 0.09236777 0.06489585
## 5 -0.3801649 0.47603246 0.03387451
## 6 -7.6841103 0.00000000 0.00000000

Because the hatvalue for 1941 in the robust second-stage fit is 0, so is the Cook’s \(D\) for this case, and the “bubble” for 1941 has 0 area (and is therefore invisible). The bad case, however, have a very large (negative) studentized residual of nearly \(-8\).

Concluding Remarks

Careful regression analysis requires methods for looking effectively at the data. Many potential problems can be addressed by examining the data prior to fitting a regression model, decreasing (if not eliminating) the necessity for post-fit diagnostics. No doubt careful data analysts employing 2SLS have always done this. Nevertheless, having methods that allow one to subject a regression model fit by 2SLS to criticism will in at least some cases suggest improvements to the model or perhaps corrections to the data.