--- title: "Advanced IV Estimation with ivreg2r" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Advanced IV Estimation with ivreg2r} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction This vignette covers features of `ivreg2r` beyond the basic two-stage least squares workflow of `vignette("introduction")`: alternative estimators, the instrument-validity tests, weak-instrument-robust inference, two-way clustering, and the reduced form. It is organized around the canonical example groups — a Mroz labor-supply equation that carries most of the help file's testing arc, the Griliches wage equation for the help file's orthogonality example and for Baum, Schaffer & Stillman's (2007) weak-instrument and redundancy illustrations, and Klein's consumption function that separates the estimators — and adds a two-way clustered IV example on a cigarette-demand panel. The topics are: - LIML, Fuller, and k-class estimation, including the canonical Klein consumption function - Instrument-validity tests: overidentification, the orthogonality C-statistic, and endogeneity - Weak-instrument-robust inference (the Anderson-Rubin and Stock-Wright statistics) and instrument redundancy - Two-way clustering, including with instrumental variables - Reduced-form regression and degrees-of-freedom adjustments For HAC and AC standard errors, the `l()`/`d()` time-series operators, GMM and the continuously-updated estimator, the panel variance estimators (Kiefer, Driscoll-Kraay, and cluster-plus-kernel), partialling out regressors, and user-supplied weighting matrices, see `vignette("time-series-gmm")`. The four weight types are documented in `vignette("introduction")` and in `?ivreg2`. Five bundled datasets appear below: the Mroz (1987) female labor-supply data (`mroz`), Klein's (1950) Model I macroeconomic series (`klein`), the Griliches (1976) young-men's wage data (`griliches`), the NLS young-women panel (`nlswork`, distributed with Stata's `[XT]` manual), and the Baltagi cigarette-demand panel (`cigar`). ```{r setup, message = FALSE} library(ivreg2r) library(dplyr) data(mroz) data(klein) data(griliches) data(nlswork) ``` ## The Mroz baseline The Mroz (1987) dataset contains `r nrow(mroz)` married women from the 1975 Panel Study of Income Dynamics, of whom `r sum(mroz$inlf == 1)` were working (`inlf == 1`) and have observed wages. Years of education (`educ`) is potentially endogenous in a wage equation because of ability bias. Following the illustrative specification used throughout Stata's `ivreg2` help file and in Baum, Schaffer & Stillman (2007, pp. 471--472), the excluded instruments are `age` and the numbers of children at home under 6 (`kidslt6`) and 6--18 (`kidsge6`). These exclusion restrictions are illustrative and contestable — children plausibly affect wage offers through labor-supply channels — and the help file adopts them precisely because they make the instrument-testing machinery interesting. ```{r mroz-subset} mroz_work <- mroz |> filter(inlf == 1) nrow(mroz_work) ``` The baseline regresses log wages on experience and its square, with education instrumented by the three excluded instruments. Three excluded instruments for one endogenous regressor leave two overidentifying restrictions, which is what makes the model useful for comparing estimators and testing instruments: ```{r mroz-baseline} fit_2sls <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work) summary(fit_2sls) ``` ## LIML, Fuller, and k-class estimation ### LIML Limited-information maximum likelihood (LIML) is an alternative to 2SLS that behaves better under weak instruments: it is approximately median-unbiased and its Wald test suffers less size distortion (Stock & Yogo, 2005). LIML has no finite moments of any order, whatever the degree of overidentification (Mariano & Sawa, 1972), so individual LIML estimates can be extreme in small samples; Fuller's modification, introduced below, restores finite moments. Request it with `method = "liml"`: ```{r liml} fit_liml <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, method = "liml") summary(fit_liml) ``` The output reports `lambda`, the LIML eigenvalue. When `lambda` is close to 1, LIML and 2SLS give similar results, since 2SLS is the k-class member with `k = 1`. A `lambda` far above 1 is evidence against the overidentifying restrictions, not of weak instruments: `N * log(lambda)` is the Anderson-Rubin (1950) overidentification LR statistic reported in the Anderson-Rubin section below, so judge how far above 1 is far by that statistic rather than by eye. Under weak but valid instruments `lambda` stays near 1, and weak identification is diagnosed by the Cragg-Donald and Kleibergen-Paap statistics instead. Here `lambda` is `r round(fit_liml$lambda, 4)`, so LIML barely moves the education coefficient away from 2SLS. ### Fuller modification Fuller's (1977) modified LIML uses `k = lambda - alpha / (N - L)`, where `L` is the number of instruments including the exogenous regressors and `alpha` is the Fuller parameter. The standard choice `fuller = 1` is recommended by Fuller (1977) and by Davidson & MacKinnon (1993, pp. 649--650), as cited in Baum, Schaffer & Stillman (2007, p. 479): ```{r fuller} fit_fuller1 <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, fuller = 1) fit_fuller1$kclass_value # the k-class value, Stata's e(kclass) tidy(fit_fuller1) |> filter(term == "educ") ``` ### k-class estimation The k-class family nests OLS (`k = 0`), 2SLS (`k = 1`), and LIML (`k = lambda`). The `kclass` argument accepts an arbitrary `k`; the documented non-trivial choice is Nagar's (1959) bias-adjusted estimator with `k = 1 + (L - K) / N`, where `L - K` is the number of overidentifying restrictions. Rather than hand-counting `L` and `K`, derive them from the fitted baseline — `rankzz` is the rank of the full instrument matrix and `rank` is the number of estimated coefficients: ```{r kclass} N <- nobs(fit_2sls) L <- fit_2sls$rankzz K <- fit_2sls$rank k_nagar <- 1 + (L - K) / N fit_nagar <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, kclass = k_nagar) c(L = L, K = K, k = k_nagar) tidy(fit_nagar) |> filter(term == "educ") ``` Both the Fuller and Nagar estimators tend to have better finite-sample behavior than 2SLS under weak instruments, though neither is robust to violations of the i.i.d. assumption (Baum, Schaffer & Stillman, 2007, pp. 478--479). ### COVIV With LIML or k-class under robust or cluster-robust standard errors, the sandwich bread uses the k-class projection matrix by default. Setting `coviv = TRUE` computes the covariance from the 2SLS projection (`k = 1`) instead, keeping the LIML residuals in the meat. Under i.i.d. errors, LIML with COVIV coincides with the continuously-updated estimator — the demonstration appears in the Klein suite below: ```{r coviv} fit_coviv <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, method = "liml", vcov = "robust", small = TRUE, coviv = TRUE) all.equal(coef(fit_coviv), coef(fit_liml)) tidy(fit_coviv) |> filter(term == "educ") |> select(term, estimate, std.error) ``` The asserted identity makes the division of labor explicit: COVIV changes only the covariance, so the coefficients are exactly the LIML coefficients, while the displayed standard error comes from the 2SLS-bread sandwich. ### Choosing among the estimators The case for LIML and its relatives is a many-instruments argument. Two-stage least squares is biased toward OLS, and that bias grows with the number of overidentifying restrictions relative to the strength of the instruments; LIML removes its leading term, which is why LIML tolerates weaker instruments before its Wald test becomes unreliable (Stock & Yogo, 2005). When the model is just-identified the argument is moot: with one excluded instrument per endogenous regressor, LIML, Fuller, and 2SLS coincide numerically, and the advantage is entirely a feature of overidentified models. Under genuinely weak instruments LIML's lack of finite moments cuts the other way, since individual estimates can be far from the truth even while the estimator stays median-unbiased. The ranking of LIML over 2SLS is itself established within the i.i.d. Stock-Yogo framework (Baum, Schaffer & Stillman, 2007, pp. 489--490), and switching estimator because a pretest flagged weak instruments carries the usual consequences of pre-testing for the inference that follows. On the Mroz baseline these distinctions barely register, because `lambda` is essentially 1 and all four estimators return nearly the same education coefficient. To see the estimators genuinely diverge we need a small, heavily overidentified problem, which is exactly what the help file's Klein example provides. ### The Klein consumption function Klein's (1950) Model I consumption function regresses consumption on lagged profits and treats current profits and the total wage bill as endogenous, with six excluded instruments: government spending, indirect taxes plus net exports, a linear trend, the government wage bill, lagged capital, and lagged total income. This is the canonical LIML example from Stata's `ivreg2` help file. Two of the terms are lags built on the fly with the `l()` time-series operator against the year variable `yr`; the operator is documented in `vignette("time-series-gmm")`. ```{r klein-formula} klein_form <- consump ~ l(profits, 1) | profits + wagetot | govt + taxnetx + year + wagegovt + capital1 + l(totinc, 1) fit_k_liml <- ivreg2(klein_form, data = klein, tvar = "yr", method = "liml") summary(fit_k_liml) ``` With two endogenous regressors and six excluded instruments the model has four overidentifying restrictions on only `r nobs(fit_k_liml)` observations, and here `lambda` is `r round(fit_k_liml$lambda, 4)`, well above the near-1 value on Mroz. That is the regime in which the choice of estimator matters. Fitting 2SLS, LIML, Fuller(1), and Nagar's k-class on the same equation shows the two endogenous coefficients moving with the estimator: ```{r klein-suite} fit_k_2sls <- ivreg2(klein_form, data = klein, tvar = "yr") fit_k_fuller <- ivreg2(klein_form, data = klein, tvar = "yr", method = "liml", fuller = 1) # Nagar's k = 1 + (L - K)/N = 1 + 4/21, which the help file rounds to 1.19 fit_k_nagar <- ivreg2(klein_form, data = klein, tvar = "yr", method = "kclass", kclass = 1.19) klein_models <- list("2SLS" = fit_k_2sls, LIML = fit_k_liml, "Fuller(1)" = fit_k_fuller, Nagar = fit_k_nagar) bind_rows("2SLS" = tidy(fit_k_2sls), LIML = tidy(fit_k_liml), "Fuller(1)" = tidy(fit_k_fuller), Nagar = tidy(fit_k_nagar), .id = "estimator") |> filter(term %in% c("profits", "wagetot")) |> select(estimator, term, estimate) ``` The estimators genuinely diverge here: the `profits` coefficient changes sign between 2SLS and LIML, and the wage-bill coefficient moves in the second decimal place — differences that are invisible on the strongly identified Mroz equation, and the reason the help file uses Klein to introduce these estimators. ### LIML with COVIV equals CUE The Klein example also demonstrates an equivalence. Under conditional-homoskedasticity (i.i.d.) weighting, LIML with COVIV and the continuously-updated estimator (CUE) are mathematically the same estimator, so they should return the same coefficients. CUE reaches its answer by numerical optimization, so the two agree up to optimizer tolerance rather than exactly: ```{r klein-coviv-cue} fit_k_coviv <- ivreg2(klein_form, data = klein, tvar = "yr", method = "liml", coviv = TRUE) fit_k_cue <- ivreg2(klein_form, data = klein, tvar = "yr", method = "cue") max(abs(coef(fit_k_coviv) - coef(fit_k_cue))) ``` On this tiny (N = `r nobs(fit_k_cue)`) problem the two agree to optimizer tolerance: the small gap printed above is the optimizer's endpoint, not a difference in the estimators, and its exact size depends on where the optimizer stops. CUE generalizes LIML with COVIV to the non-i.i.d. case, where a closed-form eigenvalue solution no longer exists. ## Instrument-validity tests The tests below follow the escalating order of the help file's testing examples: a joint overidentification test, then a test of a subset of instruments (the help file's Griliches example), then a test of a single regressor's endogeneity, then the empty-endogenous form for an included regressor, and finally the LIML overidentification statistics. Every one of these tests prints in `summary()`, and `diagnostics(fit)` returns them all as a tibble with one row per test under a stable `test` key, so pulling a statistic out for prose or a table is a one-line `filter(test == ...)` (see `vignette("introduction")` for the accessor). The `glance()` companion surfaces the headline tests (weak identification, underidentification, overidentification) as columns for model-comparison tables. The chunks below read each statistic from `diagnostics()`, dropping into `glance()` where a column is available and reads better. ### Overidentification With three excluded instruments and one endogenous regressor, the baseline has two overidentifying restrictions. The Sargan (1958) test is a joint test of those restrictions: a rejection signals that some instrument violates the orthogonality conditions but cannot attribute the failure to a specific instrument. ```{r sargan} glance(fit_2sls)[, c("overid_stat", "overid_p")] ``` ### Orthogonality (the C-statistic) The `orthog` argument tests whether a *subset* of instruments is exogenous, conditional on the rest being valid — the difference-in-Sargan C-statistic (Hayashi, 2000). Stata's help file demonstrates it on the Griliches (1976) wage equation, challenging the exogeneity of `age` and `mrt` (marital status) while maintaining `med` (mother's education) and `kww` (a test score), estimated by two-step GMM. The Griliches variables are `lw` (log wage), `s` (years of schooling), `expr` (labor-market experience), `tenure` (job tenure), `rns` (residence in the South), and `smsa` (urban residence); `factor(year)` expands to a set of year dummies, R's spelling of Stata's `i.year`: ```{r orthog} fit_orthog <- ivreg2( lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age + mrt, data = griliches, method = "gmm2s", orthog = c("age", "mrt") ) diagnostics(fit_orthog) |> filter(test == "orthog") ``` The null is that `age` and `mrt` satisfy the exclusion restrictions given that the remaining instruments do. The C-statistic rejects emphatically, so the data are inconsistent with treating `age` and `mrt` as valid exclusions in this equation. The same pair reappears below as Baum, Schaffer & Stillman's weak-instruments illustration. The equation must remain identified after the challenged instruments are dropped, and under robust or cluster-robust standard errors the test is a difference in Hansen J statistics rather than in Sargan statistics. ### Endogeneity of a regressor The `endog` argument tests whether a nominated endogenous regressor is in fact exogenous, comparing the model that instruments it against the model that treats it as exogenous. On the baseline this reproduces Stata's `endog(educ)` example: ```{r endog} fit_endog <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, endog = "educ") diagnostics(fit_endog) |> filter(test == "endogeneity") ``` With a single endogenous regressor this coincides numerically with the endogeneity test that `ivreg2()` already reports for the baseline, so `endog` earns its keep only with two or more endogenous regressors, where it can test a subset: ```{r endog-identity} c(automatic = diagnostics(fit_2sls) |> filter(test == "endogeneity") |> pull(statistic), endog_arg = diagnostics(fit_endog) |> filter(test == "endogeneity") |> pull(statistic)) ``` ### Testing an included regressor: the empty-endogenous form The C-statistic can also challenge an *included* regressor's exogeneity. Write the model with no endogenous regressors — the empty-endogenous form `y ~ exog | 0 | instruments`, Stata's `(= z1 z2)` — so the equation is estimated by OLS while the excluded instruments supply surplus orthogonality conditions, then test the regressor with `orthog`. This is the help file's C-test of `educ` in the Mroz equation: ```{r hols-orthog} fit_hols <- ivreg2( lwage ~ exper + expersq + educ | 0 | age + kidslt6 + kidsge6, data = mroz_work, orthog = "educ" ) summary(fit_hols) ``` Estimation is exactly OLS — the instruments play no role in the coefficients — but the Sargan statistic now tests all surplus conditions, and the C-statistic isolates `educ` by comparing this model against one in which `educ` is instrumented by `age`, `kidslt6`, and `kidsge6`. A small C-statistic means the data do not reject treating `educ` as exogenous, conditional on the three excluded instruments being valid. The same form with `method = "gmm2s"` and `vcov = "robust"` gives Cragg's (1983) heteroskedastic OLS (HOLS) estimator, which exploits the surplus orthogonality conditions for efficiency under heteroskedasticity. Under the default i.i.d. weighting there is no efficiency gain to exploit: the two-step weighting matrix is proportional to `(Z'Z)^{-1}` and the estimates are identical to OLS, in Stata just as here. A cluster or HAC weighting generalizes the efficiency gain to those error structures. ### Anderson-Rubin LIML overidentification Under LIML with i.i.d. errors the summary reports two further overidentification statistics. The Anderson-Rubin (1950) LR statistic is `N * log(lambda)`, and its linearized form is `N * (lambda - 1)`; both are chi-squared with `L - K` degrees of freedom under the null that the overidentifying restrictions hold. When `lambda` is close to 1 the three overidentification tests — Sargan, AR LR, and AR linearized — agree closely: ```{r ar-overid} diagnostics(fit_liml) |> filter(test %in% c("anderson_rubin_overid_lr", "anderson_rubin_overid_lin")) ``` ## Weak-instrument-robust inference ### The Stock-Wright S statistic The Stock-Wright (2000) S statistic is a weak-instrument-robust test of the joint null that the coefficients on the endogenous regressors are zero and the overidentifying restrictions hold. The Anderson-Rubin statistic tests the same null, and under it both are chi-squared with `L1` degrees of freedom, where `L1` is the number of excluded instruments: Anderson-Rubin is the Wald form of the test and S is the LM / GMM-distance form (Baum, Schaffer & Stillman, 2007, pp. 491--492). Unlike Wald tests on 2SLS coefficients, both control size under their maintained assumptions even when instruments are weak. The S statistic is the value of the CUE objective with the exogenous regressors partialled out, and it is computed for every IV model: ```{r stock-wright} diagnostics(fit_2sls) |> filter(test == "stock_wright") ``` ### When instruments are genuinely weak Baum, Schaffer & Stillman (2007, pp. 492--493) illustrate the weak-instrument problem on the Griliches equation, instrumenting `iq` (a measure of ability) with `age` and `mrt` alone, under year dummies and a heteroskedasticity-robust variance: ```{r weak-instruments} weak_form <- lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | age + mrt fit_weak <- ivreg2(weak_form, data = griliches, vcov = "robust") summary(fit_weak) sy10 <- diagnostics(fit_weak) |> filter(test == "sy_iv_size_10") |> pull(statistic) sy10 ``` ```{r weak-values, include = FALSE} kp_f_weak <- diagnostics(fit_weak) |> filter(test == "weak_id_robust") |> pull(statistic) ``` Because the variance is robust, the relevant weak-identification statistic is the Kleibergen-Paap rk Wald F (`r round(kp_f_weak, 2)`), far below the Stock-Yogo 10% maximal-size critical value of `r sy10`. The Stock-Yogo critical values are tabulated for the Cragg-Donald statistic under i.i.d. errors, so under a robust variance they apply to the rk statistic only with caution, or give way to the Staiger-Stock rule of thumb that the first-stage F should exceed 10 (Baum, Schaffer & Stillman, 2007, p. 490). With instruments this weak, the Stock-Wright S and Anderson-Rubin tests remain valid where conventional 2SLS Wald inference rests on a strong-identification approximation that has failed. These robust tests keep correct size but lose power as instruments weaken (Baum, Schaffer & Stillman, 2007, p. 491), and because the null is joint, a rejection can reflect either nonzero coefficients on the endogenous regressors or failure of the overidentifying restrictions — the test cannot say which. ### Redundant instruments Redundancy is a separate question from validity: it asks whether an instrument adds any identifying information once the others are present. An instrument is redundant when dropping it costs no asymptotic efficiency, and on a weak equation it is worth knowing whether the few available instruments are all pulling their weight. Baum, Schaffer & Stillman (2007, p. 492) test `mrt` for redundancy in exactly this weak specification with `redundant`: ```{r redundant} fit_redund <- ivreg2(weak_form, data = griliches, vcov = "robust", redundant = "mrt") diagnostics(fit_redund) |> filter(test == "redundancy") ``` ```{r redund-values, include = FALSE} redund_p <- diagnostics(fit_redund) |> filter(test == "redundancy") |> pull(p_value) ``` The LM statistic is essentially zero (p = `r round(redund_p, 2)`), the published result: `mrt` provides no useful information for identifying the coefficient on `iq`, so identification here rests on `age` alone. Redundancy is about efficiency rather than validity, and it is a conditional concept — an instrument redundant given one set of instruments may cease to be so given another. ### Building the S statistic by hand with `b0` The `b0` argument evaluates the CUE objective at a fixed parameter vector rather than optimizing it, and doing so reconstructs the Stock-Wright S statistic directly. The help file's construction sets the endogenous coefficient (`educ`) to zero and the exogenous coefficients to their OLS values; fixing the exogenous coefficients at OLS is equivalent to partialling them out, so the objective value at that point is the S statistic for the joint null: ```{r b0} ols_fit <- ivreg2(lwage ~ exper + expersq, data = mroz_work) b0 <- c(educ = 0, coef(ols_fit)) fit_b0 <- ivreg2( lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, vcov = "robust", b0 = b0 ) diagnostics(fit_b0) |> filter(test == "overid") # the objective value, Stata's e(j) ``` The identity is easy to confirm against the S statistic that the ordinary robust fit reports automatically: ```{r b0-identity} fit_reg <- ivreg2( lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, vcov = "robust" ) b0_objective <- diagnostics(fit_b0) |> filter(test == "overid") |> pull(statistic) stock_wright_S <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(statistic) all.equal(b0_objective, stock_wright_S) diagnostics(fit_reg) |> filter(test == "stock_wright") ``` ```{r b0-values, include = FALSE} sw_p <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(p_value) sw_df <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(df) ``` One caution: the p-value reported beside the raw objective — the `p_value` on the `overid` row — is calibrated against `L - K` overidentification degrees of freedom, which is the wrong reference distribution for this use. Display the raw objective, as Stata's `di e(j)` does, and take the S-statistic p-value (`r round(sw_p, 4)` here, on L1 = `r sw_df` degrees of freedom) from the regular fit. ### What the package does and does not offer The Anderson-Rubin and Stock-Wright statistics test one particular null — that the endogenous coefficients are zero and the instruments valid. `ivreg2r` reports these joint tests, but it does not test non-zero nulls by re-optimizing the objective at a candidate value (the `b0` mechanism evaluates the objective at a point but is not wired into a search), and it does not invert the Anderson-Rubin test into a confidence region. A reader whose Wald confidence interval is untrustworthy under weak instruments therefore has a size-correct test of the zero null but not, from this package, a ready weak-instrument-robust confidence set. Constructing one means inverting a weak-robust test — the Anderson-Rubin statistic, Kleibergen's LM statistic, or the conditional likelihood-ratio test of Moreira (2003), whose properties are analyzed in Andrews, Moreira & Stock (2006). Stata's `weakiv` command (Finlay, Magnusson & Schaffer) builds these confidence sets for linear IV models with any number of endogenous regressors, and in R the `ivmodel` package computes Anderson-Rubin and conditional-likelihood-ratio confidence intervals for the single-endogenous-regressor case; either accepts the kind of specification used here. ## Two-way clustering ### The nlswork panel The NLS young-women panel is Stata's `[XT]`-manual extract (`webuse nlswork`): person-year observations on women who were 14--24 in 1968 and were interviewed annually through 1988. Panel data of this shape can be correlated within a person over time and within a year across persons, and two-way clustering accounts for both. ```{r nlswork-overview} nrow(nlswork) length(unique(nlswork$idcode)) # persons length(unique(nlswork$year)) # interview years ``` ### One-way versus two-way clustering Compare one-way clustering by person with two-way clustering by person and interview year, the help file's own paired example. Both fits below are ordinary least squares with no instruments — the formula has no `|` groups — and the `clusters` argument behaves identically once the formula gains endogenous and excluded-instrument parts, as the IV example that follows shows. Fitting on the full data frame restricts to the complete cases the formula needs, so no manual subsetting is required: ```{r twoway} fit_1way <- ivreg2( ln_wage ~ grade + age + ttl_exp + tenure, data = nlswork, clusters = ~ idcode ) fit_2way <- ivreg2( ln_wage ~ grade + age + ttl_exp + tenure, data = nlswork, clusters = ~ idcode + year ) ``` Of the `r format(nrow(nlswork), big.mark = ",")` person-year rows, `r format(nobs(fit_1way), big.mark = ",")` have complete data and enter both fits. Compare the standard errors: ```{r twoway-compare} tibble( term = tidy(fit_1way)$term, se_1way = tidy(fit_1way)$std.error, se_2way = tidy(fit_2way)$std.error ) |> mutate(ratio = se_2way / se_1way) ``` Two-way clustering uses the Cameron, Gelbach & Miller (2011) decomposition `V_twoway = V_cluster1 + V_cluster2 - V_intersection`. This estimator needs both cluster dimensions to grow for consistency, so the effective number of clusters governing the asymptotics is the smaller of the two, `min(M1, M2)`. Here there are `r format(fit_2way$n_clusters1, big.mark = ",")` person clusters and only `r fit_2way$n_clusters2` interview-year clusters, so `min(M1, M2)` is `r fit_2way$n_clusters`: the year dimension, not the much larger person dimension, limits how far the two-way asymptotics can be trusted. Two-way clustering helps only when the two variables define distinct, crossed categories, as here, where every person appears in several years and every year holds many persons. That panel-by-time pairing is the typical two-way application. When one category is nested within the other — counties within states, say — the two-way variance reduces to the single-cluster variance on the outer category, and the second dimension adds nothing. ### Two-way clustering with instrumental variables Two-way clustering composes with every estimator, so we pair it with 2SLS on a panel demand equation. The Baltagi cigarette-demand data (`data(cigar)`) are an annual panel of 46 U.S. states over 1963--1992, assembled by Baltagi and Levin (1992) and revisited by Baltagi, Griffin and Xiong (2000). The specification regresses log per-capita cigarette sales on the log real price of cigarettes, with log real per-capita disposable income as an exogenous control. Price is plausibly endogenous, because state demand shocks feed back into prices, so it is instrumented by the minimum real price in adjoining states — a spatial cost-shifter in the Baltagi-Levin tradition — together with the one-year lag of the state's own real price. Clustering on both state and year allows arbitrary within-state serial correlation and arbitrary within-year correlation across states. ```{r two-way-iv} data(cigar) cigar <- cigar |> mutate( lsales = log(sales), lrprice = log(price / cpi), lrndi = log(ndi / cpi), lrpimin = log(pimin / cpi) ) fit_iv_2way <- ivreg2( lsales ~ lrndi | lrprice | lrpimin + l(lrprice, 1), data = cigar, tvar = "year", ivar = "state", clusters = ~ state + year ) summary(fit_iv_2way) ``` The `l(lrprice, 1)` term uses the panel lag operator, which needs the panel structure declared through `tvar` and `ivar`, and which is documented in `vignette("time-series-gmm")`. The summary reports both cluster dimensions — `r fit_iv_2way$n_clusters1` states and `r fit_iv_2way$n_clusters2` years, the first year lost to the lag — and the effective cluster count governing the asymptotics is again the smaller dimension, the `r fit_iv_2way$n_clusters` years. The automatic diagnostics carry through the two-way variance unchanged: the Kleibergen-Paap underidentification and weak-identification statistics and the Hansen J of the single overidentifying restriction are all computed from the two-way cluster-robust covariance. This is a methods illustration rather than a substantive demand study, and the price elasticity should be read in that light. ## Reduced-form regression The reduced form regresses the outcome directly on the full instrument set — exogenous regressors plus excluded instruments — bypassing the first stage. ### Single-equation reduced form Use `reduced_form = "rf"` to store the outcome-equation reduced form alongside the structural estimates: ```{r rf} # Same robust fit as fit_reg above, re-estimated with the reduced form stored fit_rf <- ivreg2( lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, vcov = "robust", reduced_form = "rf" ) fit_rf$reduced_form$coefficients ``` The stored reduced form is not just bookkeeping: its F-statistic on the excluded instruments *is* the Anderson-Rubin F. The AR test is precisely this reduced-form F-test, which is why it stays valid under weak instruments, and the identity holds exactly: ```{r rf-ar-identity} rf_F <- fit_rf$reduced_form$f_stat ar_F <- diagnostics(fit_rf) |> filter(test == "anderson_rubin_f") |> pull(statistic) all.equal(rf_F, ar_F) c(rf_F = rf_F, ar_F = ar_F) ``` ### System reduced form Use `reduced_form = "system"` to estimate the full system — one equation for each endogenous variable plus the outcome — with a cross-equation covariance matrix. The cross-equation covariance is what joint hypotheses spanning the first-stage and outcome reduced forms require, and storing the system reproduces Stata's undocumented `savesfirst` behavior: ```{r rf-system} fit_system <- ivreg2( lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6, data = mroz_work, reduced_form = "system" ) # Equation names fit_system$reduced_form$depvar # Cross-equation VCV dimensions dim(fit_system$reduced_form$vcov) ``` ## Degrees-of-freedom adjustments When fixed effects or other controls have been partialled out before estimation, `dofminus` and `sdofminus` adjust the degrees of freedom without re-estimating: - `dofminus` is subtracted from `N` in the large-sample variance formulas, so `sigma^2 = RSS / (N - dofminus)`. It is the large-sample adjustment (for example, the number of fixed effects), and it is how `xtivreg2, fe` reports correct statistics after within-demeaning consumes one degree of freedom per panel unit. - `sdofminus` is subtracted from the residual degrees of freedom alongside `K` in the small-sample corrections. It is the small-sample adjustment (for example, the number of partialled-out regressors), and it is what Stata's `partial()` sets internally. Most of the time neither option needs to be set by hand. The built-in `partial =` argument — Stata's `partial()` — removes named regressors inside `ivreg2()` and sets `sdofminus` for you. On the Griliches equation, partialling out the year dummies leaves the remaining coefficients exactly where the full regression put them: ```{r partial} fit_full <- ivreg2( lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age, data = griliches, small = TRUE ) fit_partial <- ivreg2( lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age, data = griliches, small = TRUE, partial = "factor(year)" ) shared <- c("s", "expr", "tenure", "rns", "smsa", "iq") all.equal(coef(fit_full)[shared], coef(fit_partial)[shared]) ``` The standard errors match too, because `partial =` set `sdofminus` to the number of partialled-out columns automatically; passing `nopartialsmall = TRUE` (Stata's `nopartialsmall`) would suppress that small-sample adjustment: ```{r partial-se} all.equal( tidy(fit_full) |> filter(term %in% shared) |> pull(std.error), tidy(fit_partial) |> filter(term %in% shared) |> pull(std.error) ) ``` You set `dofminus` and `sdofminus` by hand only when the data were partialled or demeaned *outside* `ivreg2()` — the situation a `xtivreg2, fe` user migrates into, having already within-transformed the panel before fitting. Then `ivreg2()` cannot see how many degrees of freedom were spent, so you declare them: `sdofminus` for data residualized on a set of controls, `dofminus` for the within (fixed-effects) transformation. The `dofminus` case is the recipe `xtivreg2, fe` automates. Demean each variable by its group and tell `ivreg2()` how many group effects were absorbed. Treating the seven Griliches survey years as the groups, subtract each year's mean from every variable and add back the overall mean so the intercept survives --- the within transformation that `xtivreg2, fe` performs --- then declare that the G - 1 absorbed year effects consumed degrees of freedom. The demeaned regression reproduces the dummy-variable fit exactly: shared coefficients, standard errors, sigma, and residual degrees of freedom. ```{r dofminus-fe} gril_within <- griliches |> group_by(year) |> mutate( lw = lw - mean(lw) + mean(griliches$lw), s = s - mean(s) + mean(griliches$s), expr = expr - mean(expr) + mean(griliches$expr), tenure = tenure - mean(tenure) + mean(griliches$tenure) ) |> ungroup() G <- length(unique(griliches$year)) fit_lsdv <- ivreg2(lw ~ s + expr + tenure + factor(year), data = griliches, small = TRUE) fit_fe <- ivreg2(lw ~ s + expr + tenure, data = gril_within, small = TRUE, dofminus = G - 1) slopes <- c("s", "expr", "tenure") all.equal(coef(fit_lsdv)[slopes], coef(fit_fe)[slopes]) all.equal( tidy(fit_lsdv) |> filter(term %in% slopes) |> pull(std.error), tidy(fit_fe) |> filter(term %in% slopes) |> pull(std.error) ) c(sigma_lsdv = fit_lsdv$sigma, sigma_fe = fit_fe$sigma, df_lsdv = fit_lsdv$df.residual, df_fe = fit_fe$df.residual) ``` For a migrating `xtivreg2, fe` user the practical rule is: pass the number of absorbed effects beyond the constant, and use the identity above as the check that you passed the right one. ## Estimator equivalences Two equivalences involving LIML are worth stating here; the authoritative table of variance-estimator equivalences (Driscoll-Kraay, Kiefer, cluster-robust) lives in `vignette("time-series-gmm")`. - **LIML with COVIV = CUE** under conditional-homoskedasticity weighting: demonstrated on the Klein data above, with agreement up to optimizer tolerance. - **LIML = 2SLS when exactly identified**, as noted in the estimator-choice discussion above. ## Publication tables with modelsummary The [`modelsummary`](https://modelsummary.com/) package produces publication-quality regression tables from any model with `tidy()` and `glance()` methods, and it is where the estimator-comparison workflow belongs. `ivreg2r` objects work with it directly. ```{r modelsummary-setup, eval = requireNamespace("modelsummary", quietly = TRUE)} library(modelsummary) ``` ### A single model ```{r ms-basic, eval = requireNamespace("modelsummary", quietly = TRUE)} modelsummary(fit_2sls, output = "markdown") ``` ### Comparing estimators side by side The natural side-by-side comparison is the Klein suite, where the estimator choice visibly matters --- reusing the `klein_models` list built above: ```{r ms-compare, eval = requireNamespace("modelsummary", quietly = TRUE)} modelsummary(klein_models, output = "markdown", gof_map = c("nobs", "r.squared", "adj.r.squared", "sigma")) ``` ### Exponentiated coefficients For log-linear models, `exponentiate = TRUE` reports multiplicative effects (`exp(beta)`). `modelsummary` transforms the standard errors to the exponentiated scale by the delta method, so the displayed standard error is `exp(beta) * se(beta)` and the table is internally consistent on the multiplicative scale. This transformation is `modelsummary`'s own: the package's `tidy(fit, exponentiate = TRUE)` follows the broom convention of exponentiating estimates and confidence bounds while leaving `std.error` on the log scale. ```{r ms-exponentiate, eval = requireNamespace("modelsummary", quietly = TRUE)} modelsummary(fit_2sls, output = "markdown", exponentiate = TRUE) ``` ### Customizing output Rename coefficients and select goodness-of-fit statistics: ```{r ms-custom, eval = requireNamespace("modelsummary", quietly = TRUE)} modelsummary( list("2SLS" = fit_2sls, "LIML" = fit_liml), output = "markdown", coef_rename = c(educ = "Education", exper = "Experience", expersq = "Experience²"), gof_map = c("nobs", "r.squared", "sigma", "weak_id_stat") ) ``` `modelsummary` supports LaTeX, HTML, Word, and gt output; see `?modelsummary::modelsummary` for the full list. ## Stata migration For a Stata-to-R translation table covering the most common options, see the migration section in `vignette("introduction")`; for the full argument list, see `?ivreg2`. Factor variables are handled through R's formula interface: `factor(year)` stands in for Stata's `i.year` dummies, as glossed above. ## References - Anderson, T.W. and Rubin, H. (1949). "Estimation of the Parameters of a Single Equation in a Complete System of Stochastic Equations." *Annals of Mathematical Statistics*, 20(1), 46--63. - Anderson, T.W. and Rubin, H. (1950). "The Asymptotic Properties of Estimates of the Parameters of a Single Equation in a Complete System of Stochastic Equations." *Annals of Mathematical Statistics*, 21(4), 570--582. - Andrews, D.W.K., Moreira, M.J., and Stock, J.H. (2006). "Optimal Two-Sided Invariant Similar Tests for Instrumental Variables Regression." *Econometrica*, 74(3), 715--752. - Baltagi, B.H. and Levin, D. (1992). "Cigarette Taxation: Raising Revenues and Reducing Consumption." *Structural Change and Economic Dynamics*, 3(2), 321--335. - Baltagi, B.H., Griffin, J.M., and Xiong, W. (2000). "To Pool or Not to Pool: Homogeneous Versus Heterogeneous Estimators Applied to Cigarette Demand." *Review of Economics and Statistics*, 82(1), 117--126. - Baum, C.F., Schaffer, M.E., and Stillman, S. (2007). "Enhanced Routines for Instrumental Variables/Generalized Method of Moments Estimation and Testing." *Stata Journal*, 7(4), 465--506. - Cameron, A.C., Gelbach, J.B., and Miller, D.L. (2011). "Robust Inference with Multiway Clustering." *Journal of Business & Economic Statistics*, 29(2), 238--249. - Cragg, J.G. (1983). "More Efficient Estimation in the Presence of Heteroscedasticity of Unknown Form." *Econometrica*, 51(3), 751--763. - Davidson, R. and MacKinnon, J.G. (1993). *Estimation and Inference in Econometrics*. Oxford University Press. - Fuller, W.A. (1977). "Some Properties of a Modification of the Limited Information Estimator." *Econometrica*, 45(4), 939--953. - Griliches, Z. (1976). "Wages of Very Young Men." *Journal of Political Economy*, 84(4, Part 2), S69--S85. - Hayashi, F. (2000). *Econometrics*. Princeton University Press. - Klein, L.R. (1950). *Economic Fluctuations in the United States, 1921--1941*. Cowles Commission Monograph No. 11. Wiley. - Mariano, R.S. and Sawa, T. (1972). "The Exact Finite-Sample Distribution of the Limited-Information Maximum Likelihood Estimator in the Case of Two Included Endogenous Variables." *Journal of the American Statistical Association*, 67(337), 159--163. - Moreira, M.J. (2003). "A Conditional Likelihood Ratio Test for Structural Models." *Econometrica*, 71(4), 1027--1048. - Mroz, T.A. (1987). "The Sensitivity of an Empirical Model of Married Women's Hours of Work to Economic and Statistical Assumptions." *Econometrica*, 55(4), 765--799. - Nagar, A.L. (1959). "The Bias and Moment Matrix of the General k-Class Estimators of the Parameters in Simultaneous Equations." *Econometrica*, 27(4), 575--595. - Sargan, J.D. (1958). "The Estimation of Economic Relationships Using Instrumental Variables." *Econometrica*, 26(3), 393--415. - Stock, J.H. and Wright, J.H. (2000). "GMM with Weak Identification." *Econometrica*, 68(5), 1055--1096. - Stock, J.H. and Yogo, M. (2005). "Testing for Weak Instruments in Linear IV Regression." In D.W.K. Andrews and J.H. Stock (Eds.), *Identification and Inference for Econometric Models: Essays in Honor of Thomas Rothenberg*. Cambridge University Press.