# Advanced IV Estimation with ivreg2r

## 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`).

    library(ivreg2r)
    library(dplyr)
    data(mroz)
    data(klein)
    data(griliches)
    data(nlswork)

## The Mroz baseline

The Mroz (1987) dataset contains 753 married women from the 1975 Panel Study of Income Dynamics, of whom 428 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.

    mroz_work <- mroz |> filter(inlf == 1)
    nrow(mroz_work)
    #> [1] 428

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:

    fit_2sls <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                       data = mroz_work)
    summary(fit_2sls)
    #>
    #> 2SLS Estimation
    #>
    #> Call:
    #> ivreg2(formula = lwage ~ exper + expersq | educ | age + kidslt6 +
    #>     kidsge6, data = mroz_work)
    #>
    #> Observations: 428
    #> VCV type:     Classical (iid)
    #>
    #> Coefficients:
    #>               Estimate Std. Error z value Pr(>|z|)
    #> (Intercept) -0.3848718  1.0115512  -0.380  0.70359
    #> educ         0.0964002  0.0814278   1.184  0.23646
    #> exper        0.0421930  0.0138831   3.039  0.00237 **
    #> expersq     -0.0008323  0.0004204  -1.980  0.04773 *
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      0.1556
    #> Adj. R-squared: 0.1496
    #> Wald chi2(3):  7.5 (p = 0.0001)
    #> Root MSE:       0.6638
    #>
    #> Underidentification test (Anderson canon. corr. LM statistic):
    #>   Chi-sq(3) = 12.82 (p = 0.0051)
    #>
    #> Weak identification test:
    #>   Cragg-Donald Wald F:           4.34
    #>   Stock-Yogo critical values (IV size):
    #>      10%  maximal IV size       22.30
    #>      15%  maximal IV size       12.83
    #>      20%  maximal IV size       9.54
    #>      25%  maximal IV size       7.80
    #>   Stock-Yogo critical values (IV relative bias):
    #>      5%  maximal IV relative bias   13.91
    #>      10%  maximal IV relative bias   9.08
    #>      20%  maximal IV relative bias   6.46
    #>      30%  maximal IV relative bias   5.39
    #>
    #> Overidentification test (Sargan):
    #>   Chi-sq(2) = 0.70 (p = 0.7042)
    #>
    #> Weak-instrument-robust inference:
    #>   H0: B1=0 and orthogonality conditions are valid
    #>   Anderson-Rubin Wald F(3,422) = 0.61 (p = 0.6076)
    #>   Anderson-Rubin Wald Chi-sq(3) = 1.86 (p = 0.6016)
    #>   Stock-Wright LM S Chi-sq(3) = 1.85 (p = 0.6033)
    #>
    #> Endogeneity test:
    #>   Chi-sq(1) = 0.02 (p = 0.8899)
    #>   Tested: educ
    #>
    #> First-stage diagnostics:
    #>   Endogenous        F-stat   p-value  Partial R2  Shea PR2      SW F      AP F
    #>   educ                4.34    0.0050      0.0299    0.0299      4.34      4.34
    #>
    #> Instrumented:          educ
    #> Included instruments:  exper, expersq
    #> Excluded instruments:  age, kidslt6, kidsge6

## 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"`:

    fit_liml <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                       data = mroz_work, method = "liml")
    summary(fit_liml)
    #>
    #> LIML Estimation
    #>
    #> Call:
    #> ivreg2(formula = lwage ~ exper + expersq | educ | age + kidslt6 +
    #>     kidsge6, data = mroz_work, method = "liml")
    #>
    #> Observations: 428
    #> VCV type:     Classical (iid)
    #> lambda:       1.001642
    #> kclass:       1.001642
    #>
    #> Coefficients:
    #>               Estimate Std. Error z value Pr(>|z|)
    #> (Intercept) -0.3769294  1.0394246  -0.363  0.71688
    #> educ         0.0957581  0.0836906   1.144  0.25254
    #> exper        0.0422292  0.0139270   3.032  0.00243 **
    #> expersq     -0.0008335  0.0004220  -1.975  0.04827 *
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      0.1555
    #> Adj. R-squared: 0.1495
    #> Wald chi2(3):  7.5 (p = 0.0001)
    #> Root MSE:       0.6638
    #>
    #> Underidentification test (Anderson canon. corr. LM statistic):
    #>   Chi-sq(3) = 12.82 (p = 0.0051)
    #>
    #> Weak identification test:
    #>   Cragg-Donald Wald F:           4.34
    #>   Stock-Yogo critical values (LIML size):
    #>      10%  maximal LIML size      6.46
    #>      15%  maximal LIML size      4.36
    #>      20%  maximal LIML size      3.69
    #>      25%  maximal LIML size      3.32
    #>
    #> Overidentification test (Sargan):
    #>   Chi-sq(2) = 0.70 (p = 0.7042)
    #>
    #> Anderson-Rubin overidentification:
    #>   LR Chi-sq(2) = 0.702 (p = 0.7040)
    #>   Linearized Chi-sq(2) = 0.703 (p = 0.7038)
    #>
    #> Weak-instrument-robust inference:
    #>   H0: B1=0 and orthogonality conditions are valid
    #>   Anderson-Rubin Wald F(3,422) = 0.61 (p = 0.6076)
    #>   Anderson-Rubin Wald Chi-sq(3) = 1.86 (p = 0.6016)
    #>   Stock-Wright LM S Chi-sq(3) = 1.85 (p = 0.6033)
    #>
    #> Endogeneity test:
    #>   Chi-sq(1) = 0.02 (p = 0.8899)
    #>   Tested: educ
    #>
    #> First-stage diagnostics:
    #>   Endogenous        F-stat   p-value  Partial R2  Shea PR2      SW F      AP F
    #>   educ                4.34    0.0050      0.0299    0.0299      4.34      4.34
    #>
    #> Instrumented:          educ
    #> Included instruments:  exper, expersq
    #> Excluded instruments:  age, kidslt6, kidsge6

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 1.0016, 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):

    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)
    #> [1] 0.9992719
    tidy(fit_fuller1) |> filter(term == "educ")
    #> # A tibble: 1 × 7
    #>   term  estimate std.error statistic p.value conf.low conf.high
    #>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
    #> 1 educ    0.0967    0.0805      1.20   0.230  -0.0611     0.254

### 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:

    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)
    #>        L        K        k
    #> 6.000000 4.000000 1.004673
    tidy(fit_nagar) |> filter(term == "educ")
    #> # A tibble: 1 × 7
    #>   term  estimate std.error statistic p.value conf.low conf.high
    #>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
    #> 1 educ    0.0944    0.0884      1.07   0.286  -0.0789     0.268

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:

    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))
    #> [1] TRUE
    tidy(fit_coviv) |> filter(term == "educ") |> select(term, estimate, std.error)
    #> # A tibble: 1 × 3
    #>   term  estimate std.error
    #>   <chr>    <dbl>     <dbl>
    #> 1 educ    0.0958    0.0869

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")`.

    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)
    #>
    #> LIML Estimation
    #>
    #> Call:
    #> ivreg2(formula = klein_form, data = klein, method = "liml", tvar = "yr")
    #>
    #> Observations: 21
    #> VCV type:     Classical (iid)
    #> lambda:       1.498746
    #> kclass:       1.498746
    #>
    #> Coefficients:
    #>               Estimate Std. Error z value Pr(>|z|)
    #> (Intercept)   17.14766    1.84030   9.318   <2e-16 ***
    #> profits       -0.22251    0.20175  -1.103   0.2701
    #> wagetot        0.82256    0.05538  14.853   <2e-16 ***
    #> l(profits, 1)  0.39603    0.17360   2.281   0.0225 *
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      0.9566
    #> Adj. R-squared: 0.9489
    #> Wald chi2(3):  118.4 (p = 0.0000)
    #> Root MSE:       1.3953
    #>
    #> Underidentification test (Anderson canon. corr. LM statistic):
    #>   Chi-sq(5) = 12.01 (p = 0.0347)
    #>
    #> Weak identification test:
    #>   Cragg-Donald Wald F:           2.89
    #>   Stock-Yogo critical values (LIML size):
    #>      10%  maximal LIML size      4.06
    #>      15%  maximal LIML size      2.95
    #>      20%  maximal LIML size      2.63
    #>      25%  maximal LIML size      2.46
    #>
    #> Overidentification test (Sargan):
    #>   Chi-sq(4) = 6.99 (p = 0.1365)
    #>
    #> Anderson-Rubin overidentification:
    #>   LR Chi-sq(4) = 8.497 (p = 0.0750)
    #>   Linearized Chi-sq(4) = 10.474 (p = 0.0332)
    #>
    #> Weak-instrument-robust inference:
    #>   H0: B1=0 and orthogonality conditions are valid
    #>   Anderson-Rubin Wald F(6,13) = 18.01 (p = 0.0000)
    #>   Anderson-Rubin Wald Chi-sq(6) = 174.61 (p < 2.2e-16)
    #>   Stock-Wright LM S Chi-sq(6) = 18.75 (p = 0.0046)
    #>
    #> Endogeneity test:
    #>   Chi-sq(2) = 8.98 (p = 0.0112)
    #>   Tested: profits, wagetot
    #>
    #> First-stage diagnostics:
    #>   Endogenous        F-stat   p-value  Partial R2  Shea PR2      SW F      AP F
    #>   profits             2.92    0.0497      0.5742    0.5926      3.65      3.14
    #>   wagetot            38.92    0.0000      0.9473    0.9777    100.64     41.88
    #>
    #> Instrumented:          profits, wagetot
    #> Included instruments:  l(profits, 1)
    #> Excluded instruments:  govt, taxnetx, year, wagegovt, capital1, l(totinc, 1)

With two endogenous regressors and six excluded instruments the model has four overidentifying restrictions on only 21 observations, and here `lambda` is 1.4987, 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:

    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)
    #> # A tibble: 8 × 3
    #>   estimator term    estimate
    #>   <chr>     <chr>      <dbl>
    #> 1 2SLS      profits   0.0173
    #> 2 2SLS      wagetot   0.810
    #> 3 LIML      profits  -0.223
    #> 4 LIML      wagetot   0.823
    #> 5 Fuller(1) profits  -0.169
    #> 6 Fuller(1) wagetot   0.820
    #> 7 Nagar     profits  -0.0497
    #> 8 Nagar     wagetot   0.814

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:

    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)))
    #> [1] 3.937378e-07

On this tiny (N = 21) 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.

    glance(fit_2sls)[, c("overid_stat", "overid_p")]
    #> # A tibble: 1 × 2
    #>   overid_stat overid_p
    #>         <dbl>    <dbl>
    #> 1       0.702    0.704

### 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`:

    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")
    #> # A tibble: 1 × 8
    #>   test   test_name  statistic    df   df2  p_value tested_vars note
    #>   <chr>  <chr>          <dbl> <int> <int>    <dbl> <chr>       <chr>
    #> 1 orthog C (orthog)      86.6     2    NA 1.55e-19 age, mrt    <NA>

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:

    fit_endog <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                        data = mroz_work, endog = "educ")
    diagnostics(fit_endog) |> filter(test == "endogeneity")
    #> # A tibble: 1 × 8
    #>   test        test_name   statistic    df   df2 p_value tested_vars note
    #>   <chr>       <chr>           <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 endogeneity Endogeneity    0.0191     1    NA   0.890 educ        <NA>

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:

    c(automatic = diagnostics(fit_2sls) |> filter(test == "endogeneity") |> pull(statistic),
      endog_arg = diagnostics(fit_endog) |> filter(test == "endogeneity") |> pull(statistic))
    #>  automatic  endog_arg
    #> 0.01914714 0.01914714

### 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:

    fit_hols <- ivreg2(
      lwage ~ exper + expersq + educ | 0 | age + kidslt6 + kidsge6,
      data = mroz_work, orthog = "educ"
    )
    summary(fit_hols)
    #>
    #> OLS Estimation
    #>
    #> Call:
    #> ivreg2(formula = lwage ~ exper + expersq + educ | 0 | age + kidslt6 +
    #>     kidsge6, data = mroz_work, orthog = "educ")
    #>
    #> Observations: 428
    #> VCV type:     Classical (iid)
    #>
    #> Coefficients:
    #>               Estimate Std. Error z value Pr(>|z|)
    #> (Intercept) -0.5220407  0.1977017  -2.641  0.00828 **
    #> exper        0.0415665  0.0131135   3.170  0.00153 **
    #> expersq     -0.0008112  0.0003914  -2.073  0.03822 *
    #> educ         0.1074896  0.0140802   7.634 2.27e-14 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      0.1568
    #> Adj. R-squared: 0.1509
    #> Wald chi2(3):  26.3 (p = 0.0000)
    #> Root MSE:       0.6633
    #>
    #> Overidentification test (Sargan):
    #>   Chi-sq(3) = 0.72 (p = 0.8681)
    #>
    #> Orthogonality test (C-stat):
    #>   Chi-sq(1) = 0.02 (p = 0.8899)
    #>   Tested: educ
    #>
    #> Included instruments:  exper, expersq, educ
    #> Excluded instruments:  age, kidslt6, kidsge6

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:

    diagnostics(fit_liml) |>
      filter(test %in% c("anderson_rubin_overid_lr", "anderson_rubin_overid_lin"))
    #> # A tibble: 2 × 8
    #>   test                 test_name statistic    df   df2 p_value tested_vars note
    #>   <chr>                <chr>         <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 anderson_rubin_over… Anderson…     0.702     2    NA   0.704 <NA>        <NA>
    #> 2 anderson_rubin_over… Anderson…     0.703     2    NA   0.704 <NA>        <NA>

## 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:

    diagnostics(fit_2sls) |> filter(test == "stock_wright")
    #> # A tibble: 1 × 8
    #>   test         test_name         statistic    df   df2 p_value tested_vars note
    #>   <chr>        <chr>                 <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 stock_wright Stock-Wright LM S      1.85     3    NA   0.603 <NA>        <NA>

### 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:

    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)
    #>
    #> 2SLS Estimation
    #>
    #> Call:
    #> ivreg2(formula = weak_form, data = griliches, vcov = "robust")
    #>
    #> Observations: 758
    #> VCV type:     Robust
    #>
    #> Coefficients:
    #>                 Estimate Std. Error z value Pr(>|z|)
    #> (Intercept)    10.550965   2.781762   3.793 0.000149 ***
    #> iq             -0.094890   0.041890  -2.265 0.023500 *
    #> s               0.339712   0.118327   2.871 0.004092 **
    #> expr           -0.006604   0.029255  -0.226 0.821406
    #> tenure          0.084885   0.030668   2.768 0.005643 **
    #> rns            -0.376939   0.155997  -2.416 0.015678 *
    #> smsa            0.218119   0.103112   2.115 0.034399 *
    #> factor(year)67  0.007775   0.166325   0.047 0.962717
    #> factor(year)68  0.037799   0.152359   0.248 0.804061
    #> factor(year)69  0.334703   0.163799   2.043 0.041016 *
    #> factor(year)70  0.628643   0.246846   2.547 0.010875 *
    #> factor(year)71  0.444610   0.186188   2.388 0.016942 *
    #> factor(year)73  0.439027   0.166866   2.631 0.008513 **
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      -6.4195
    #> Adj. R-squared: -6.5390
    #> Wald chi2(12):  4.4 (p = 0.0000)
    #> Root MSE:       1.1676
    #>
    #> Underidentification test (Kleibergen-Paap rk LM statistic):
    #>   Chi-sq(2) = 5.90 (p = 0.0524)
    #>
    #> Weak identification test:
    #>   Cragg-Donald Wald F:           2.72
    #>   Kleibergen-Paap rk Wald F:     2.93
    #>   (Stock-Yogo critical values are for iid errors)
    #>   Stock-Yogo critical values (IV size):
    #>      10%  maximal IV size       19.93
    #>      15%  maximal IV size       11.59
    #>      20%  maximal IV size       8.75
    #>      25%  maximal IV size       7.25
    #>
    #> Overidentification test (Hansen J):
    #>   Chi-sq(1) = 1.56 (p = 0.2111)
    #>
    #> Weak-instrument-robust inference:
    #>   H0: B1=0 and orthogonality conditions are valid
    #>   Anderson-Rubin Wald F(2,744) = 46.95 (p < 2.2e-16)
    #>   Anderson-Rubin Wald Chi-sq(2) = 95.66 (p < 2.2e-16)
    #>   Stock-Wright LM S Chi-sq(2) = 71.33 (p = 0.0000)
    #>
    #> Endogeneity test:
    #>   Chi-sq(1) = 56.06 (p = 0.0000)
    #>   Tested: iq
    #>
    #> First-stage diagnostics:
    #>   Endogenous        F-stat   p-value  Partial R2  Shea PR2      SW F      AP F
    #>   iq                  2.93    0.0539      0.0073    0.0073      2.93      2.93
    #>
    #> Instrumented:          iq
    #> Included instruments:  s, expr, tenure, rns, smsa, factor(year)67, factor(year)68, factor(year)69, factor(year)70, factor(year)71, factor(year)73
    #> Excluded instruments:  age, mrt

    sy10 <- diagnostics(fit_weak) |> filter(test == "sy_iv_size_10") |> pull(statistic)
    sy10
    #> [1] 19.93

Because the variance is robust, the relevant weak-identification statistic is the Kleibergen-Paap rk Wald F (2.93), far below the Stock-Yogo 10% maximal-size critical value of 19.93. 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`:

    fit_redund <- ivreg2(weak_form, data = griliches, vcov = "robust",
                         redundant = "mrt")
    diagnostics(fit_redund) |> filter(test == "redundancy")
    #> # A tibble: 1 × 8
    #>   test       test_name           statistic    df   df2 p_value tested_vars note
    #>   <chr>      <chr>                   <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 redundancy Redundancy test (L…   0.00176     1    NA   0.967 mrt         <NA>

The LM statistic is essentially zero (p = 0.97), 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:

    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)
    #> # A tibble: 1 × 8
    #>   test   test_name statistic    df   df2 p_value tested_vars note
    #>   <chr>  <chr>         <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 overid Hansen J       1.66     2    NA   0.436 <NA>        <NA>

The identity is easy to confirm against the S statistic that the ordinary robust fit reports automatically:

    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)
    #> [1] TRUE

    diagnostics(fit_reg) |> filter(test == "stock_wright")
    #> # A tibble: 1 × 8
    #>   test         test_name         statistic    df   df2 p_value tested_vars note
    #>   <chr>        <chr>                 <dbl> <int> <int>   <dbl> <chr>       <chr>
    #> 1 stock_wright Stock-Wright LM S      1.66     3    NA   0.645 <NA>        <NA>

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 (0.6454 here, on L1 = 3 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.

    nrow(nlswork)
    #> [1] 28534
    length(unique(nlswork$idcode))  # persons
    #> [1] 4711
    length(unique(nlswork$year))    # interview years
    #> [1] 15

### 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:

    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 28,534 person-year rows, 28,099 have complete data and enter both fits. Compare the standard errors:

    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)
    #> # A tibble: 5 × 4
    #>   term         se_1way se_2way ratio
    #>   <chr>          <dbl>   <dbl> <dbl>
    #> 1 (Intercept) 0.0338   0.0431   1.28
    #> 2 grade       0.00216  0.00265  1.23
    #> 3 age         0.000945 0.00164  1.74
    #> 4 ttl_exp     0.00184  0.00268  1.46
    #> 5 tenure      0.00163  0.00305  1.87

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 4,697 person clusters and only 15 interview-year clusters, so `min(M1, M2)` is 15: 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.

    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)
    #>
    #> 2SLS Estimation
    #>
    #> Call:
    #> ivreg2(formula = lsales ~ lrndi | lrprice | lrpimin + l(lrprice,
    #>     1), data = cigar, clusters = ~state + year, tvar = "year",
    #>     ivar = "state")
    #>
    #> Observations: 1,334
    #> VCV type:     Cluster-robust
    #> Clusters:     46 (state), 29 (year)
    #>
    #> Coefficients:
    #>             Estimate Std. Error z value Pr(>|z|)
    #> (Intercept)  3.42880    0.34460   9.950  < 2e-16 ***
    #> lrprice     -0.84886    0.11442  -7.419 1.18e-13 ***
    #> lrndi        0.27957    0.07535   3.710 0.000207 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #> ---
    #> R-squared:      0.3271
    #> Adj. R-squared: 0.3261
    #> Wald chi2(2):  27.7 (p = 0.0000)
    #> Root MSE:       0.1840
    #>
    #> Underidentification test (Kleibergen-Paap rk LM statistic):
    #>   Chi-sq(2) = 12.14 (p = 0.0023)
    #>
    #> Weak identification test:
    #>   Cragg-Donald Wald F:           4766.55
    #>   Kleibergen-Paap rk Wald F:     275.63
    #>   (Stock-Yogo critical values are for iid errors)
    #>   Stock-Yogo critical values (IV size):
    #>      10%  maximal IV size       19.93
    #>      15%  maximal IV size       11.59
    #>      20%  maximal IV size       8.75
    #>      25%  maximal IV size       7.25
    #>
    #> Overidentification test (Hansen J):
    #>   Chi-sq(1) = 0.94 (p = 0.3325)
    #>
    #> Weak-instrument-robust inference:
    #>   H0: B1=0 and orthogonality conditions are valid
    #>   Anderson-Rubin Wald F(2,28) = 27.10 (p = 0.0000)
    #>   Anderson-Rubin Wald Chi-sq(2) = 56.26 (p = 0.0000)
    #>   Stock-Wright LM S Chi-sq(2) = 12.43 (p = 0.0020)
    #>
    #> Endogeneity test:
    #>   Chi-sq(1) = 0.09 (p = 0.7636)
    #>   Tested: lrprice
    #>
    #> First-stage diagnostics:
    #>   Endogenous        F-stat   p-value  Partial R2  Shea PR2      SW F      AP F
    #>   lrprice           275.63    0.0000      0.8776    0.8776    275.63    275.63
    #>
    #> Instrumented:          lrprice
    #> Included instruments:  lrndi
    #> Excluded instruments:  lrpimin, l(lrprice, 1)

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 — 46 states and 29 years, the first year lost to the lag — and the effective cluster count governing the asymptotics is again the smaller dimension, the 29 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:

    # 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
    #>   (Intercept)         exper       expersq           age       kidslt6
    #>  0.9940565569  0.0452258622 -0.0009716342 -0.0027608640  0.0175153014
    #>       kidsge6
    #> -0.0385237413

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:

    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)
    #> [1] TRUE
    c(rf_F = rf_F, ar_F = ar_F)
    #>      rf_F      ar_F
    #> 0.5795624 0.5795624

### 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:

    fit_system <- ivreg2(
      lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
      data = mroz_work, reduced_form = "system"
    )

    # Equation names
    fit_system$reduced_form$depvar
    #> [1] "lwage" "educ"

    # Cross-equation VCV dimensions
    dim(fit_system$reduced_form$vcov)
    #> [1] 12 12

## 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:

    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])
    #> [1] TRUE

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:

    all.equal(
      tidy(fit_full)    |> filter(term %in% shared) |> pull(std.error),
      tidy(fit_partial) |> filter(term %in% shared) |> pull(std.error)
    )
    #> [1] TRUE

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.

    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])
    #> [1] TRUE
    all.equal(
      tidy(fit_lsdv) |> filter(term %in% slopes) |> pull(std.error),
      tidy(fit_fe)   |> filter(term %in% slopes) |> pull(std.error)
    )
    #> [1] TRUE
    c(sigma_lsdv = fit_lsdv$sigma, sigma_fe = fit_fe$sigma,
      df_lsdv = fit_lsdv$df.residual, df_fe = fit_fe$df.residual)
    #> sigma_lsdv   sigma_fe    df_lsdv      df_fe
    #>   0.337452   0.337452 748.000000 748.000000

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` 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.

    library(modelsummary)

### A single model

    modelsummary(fit_2sls, output = "markdown")

(Intercept) | -0.385
| (1.012)
educ | 0.096
| (0.081)
exper | 0.042
| (0.014)
expersq | -0.001
| (0.000)
Num.Obs. | 428
R2 | 0.156
R2 Adj. | 0.150
vcov_type | iid
weak_id_stat | 4.34207086242851
underid_stat | 12.8158231068423
underid_p | 0.00505230978693376
overid_stat | 0.701512431708445
overid_p | 0.704155394854523

### 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:

    modelsummary(klein_models, output = "markdown",
                 gof_map = c("nobs", "r.squared", "adj.r.squared", "sigma"))

| 2SLS | LIML | Fuller(1) | Nagar
--- | --- | --- | --- | ---
(Intercept) | 16.555 | 17.148 | 17.008 | 16.711
| (1.321) | (1.840) | (1.702) | (1.437)
profits | 0.017 | -0.223 | -0.169 | -0.050
| (0.118) | (0.202) | (0.180) | (0.137)
wagetot | 0.810 | 0.823 | 0.820 | 0.814
| (0.040) | (0.055) | (0.051) | (0.044)
l(profits, 1) | 0.216 | 0.396 | 0.355 | 0.266
| (0.107) | (0.174) | (0.156) | (0.122)
Num.Obs. | 21 | 21 | 21 | 21
R2 | 0.977 | 0.957 | 0.963 | 0.973
R2 Adj. | 0.973 | 0.949 | 0.956 | 0.968
Sigma | 1.022 | 1.395 | 1.296 | 1.105

### 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.

    modelsummary(fit_2sls, output = "markdown", exponentiate = TRUE)

(Intercept) | 0.681
| (0.688)
educ | 1.101
| (0.090)
exper | 1.043
| (0.014)
expersq | 0.999
| (0.000)
Num.Obs. | 428
R2 | 0.156
R2 Adj. | 0.150
vcov_type | iid
weak_id_stat | 4.34207086242851
underid_stat | 12.8158231068423
underid_p | 0.00505230978693376
overid_stat | 0.701512431708445
overid_p | 0.704155394854523

### Customizing output

Rename coefficients and select goodness-of-fit statistics:

    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")
    )

| 2SLS | LIML
--- | --- | ---
(Intercept) | -0.385 | -0.377
| (1.012) | (1.039)
Education | 0.096 | 0.096
| (0.081) | (0.084)
Experience | 0.042 | 0.042
| (0.014) | (0.014)
Experience² | -0.001 | -0.001
| (0.000) | (0.000)
Num.Obs. | 428 | 428
R2 | 0.156 | 0.155
Sigma | 0.664 | 0.664
weak_id_stat | 4 | 4

`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.