---
title: "Instrumental Variables Estimation with ivreg2r"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Instrumental Variables Estimation with ivreg2r}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction

`ivreg2r` provides comprehensive instrumental variables and GMM estimation with automatic diagnostic tests. It supports 2SLS, LIML, Fuller, k-class, two-step efficient GMM, and continuously-updated (CUE) estimators with classical, robust, cluster-robust, HAC, and Driscoll-Kraay standard errors. It is an R implementation inspired by Stata's `ivreg2` (Baum, Schaffer & Stillman), designed for applied researchers --- especially those migrating from Stata who expect rich diagnostics reported alongside their IV estimates.

This vignette covers the core workflow: OLS, 2SLS, diagnostics, robust and cluster-robust inference, weights, and broom integration. For LIML, Fuller, and k-class estimation, weak-instrument-robust inference, and model-comparison workflows, see `vignette("advanced-iv")`. For HAC/AC standard errors, GMM, CUE, partialling, and panel VCE options, see `vignette("time-series-gmm")`.

We'll use the Card (1995) dataset throughout this vignette. Card uses geographic proximity to a four-year college as an instrument for years of education in a log-wage equation, exploiting the idea that growing up near a college reduces the cost of attending and therefore shifts educational attainment. The exclusion restriction --- that proximity affects wages only through schooling --- is an identifying *assumption*, not something the data can fully verify: proximity could in principle be correlated with unobserved determinants of wages, a concern Card (1995) discusses at length.

```{r setup, message = FALSE}
library(ivreg2r)
library(dplyr)
data(card)
```

## OLS baseline

Give `ivreg2()` an ordinary regression formula with no `|` separators --- the same `y ~ x1 + x2` you would pass to `lm()` --- and it fits by ordinary least squares. In Stata this is the OLS special case of the same command: `ivreg2 lwage educ exper ...` with no `(educ = ...)` block. We use the control set from Card (1995, Table 2), which is also the specification in Example 15.4 of Wooldridge (2020): years of schooling, experience and its square, race, urban residence (`smsa`), current region (`south`), and the 1966 residence indicators (`smsa66` together with the census-region dummies). The Card data ship the census-region indicators already expanded as separate dummies, `reg661`--`reg669` with `reg661` the omitted base, so we list `reg662`--`reg669` directly:

```{r ols}
fit_ols <- ivreg2(
  lwage ~ educ + exper + expersq + black + smsa + south + smsa66 +
    reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + reg669,
  data = card
)
summary(fit_ols)
```

The point estimates are identical to `lm()`, and the model uses the same estimation infrastructure as the IV models below, so switching between OLS and 2SLS requires only a formula change. The default *reporting* follows Stata's conventions, however: z statistics rather than t statistics, standard errors without the finite-sample N - K correction, and Root MSE computed as sqrt(RSS/N). Setting `small = TRUE` reproduces `lm()`'s standard errors, t statistics, and sigma exactly. (For *weighted* models, coefficients and `small = TRUE` standard errors still match `lm()`, but sigma is defined differently --- see Weighted estimation below.)

## Two-stage least squares

To estimate by two-stage least squares, mark the endogenous regressor and its instruments with the `|` separator. Where Stata writes `ivreg2 lwage <controls> (educ = nearc4)`, `ivreg2r` writes `lwage ~ <controls> | educ | nearc4`: the exogenous regressors first, then the endogenous regressor after the first `|`, then the excluded instruments after the second. Each variable is listed once, and the exogenous regressors are used as instruments automatically.

Here we instrument `educ` with `nearc4` (grew up near a four-year college). The robust, first-stage, and weighted fits below all reuse this same specification, so we write the formula out once and store it:

```{r tsls}
iv_formula <- lwage ~ exper + expersq + black + smsa + south + smsa66 +
  reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + reg669 |
  educ | nearc4

fit_iv <- ivreg2(iv_formula, data = card)
summary(fit_iv)
```

All `r format(nobs(fit_iv), big.mark = ",")` observations are used, and the IV estimate of the return to schooling (`r round(coef(fit_iv)["educ"], 4)`, versus `r round(coef(fit_ols)["educ"], 4)` by OLS) replicates the published estimates for this specification (Card, 1995, Table 2; Wooldridge, 2020, Example 15.4).

The formula has three groups, separated by the two `|` bars, mapping onto the pieces of the Stata command:

| Group | Contents | In `ivreg2 lwage ... (educ = nearc4)` |
|-------|----------|----------------------------------------|
| Before the first `|` | Dependent variable, then exogenous regressors | `lwage`, then the controls |
| Between the two `|` | Endogenous regressor(s) | `educ` |
| After the second `|` | Excluded instruments | `nearc4` |

The full instrument set is the union of the exogenous regressors and the excluded instruments; you never list the exogenous regressors twice.

## Understanding the diagnostics

The `summary()` output reports several diagnostic tests automatically. This section walks through each one, using the values from the fit above.

```{r diag-values, include = FALSE}
diag_iv <- diagnostics(fit_iv)
cd_f <- diag_iv |> filter(test == "weak_id") |> pull(statistic)
sy10 <- diag_iv |> filter(test == "sy_iv_size_10") |> pull(statistic)
sy15 <- diag_iv |> filter(test == "sy_iv_size_15") |> pull(statistic)
underid_stat <- diag_iv |> filter(test == "underid") |> pull(statistic)
underid_p <- diag_iv |> filter(test == "underid") |> pull(p_value)
endog_stat <- diag_iv |> filter(test == "endogeneity") |> pull(statistic)
endog_p <- diag_iv |> filter(test == "endogeneity") |> pull(p_value)
ar_chi2_p <- diag_iv |> filter(test == "anderson_rubin_chi2") |> pull(p_value)
educ_z_p <- tidy(fit_iv) |> filter(term == "educ") |> pull(p.value)
```

### Weak identification

The fit above reports a **Cragg-Donald Wald F** of `r round(cd_f, 2)`. This statistic tests whether the excluded instruments are *weak* --- that is, only weakly correlated with the endogenous regressors. Weak instruments lead to biased 2SLS estimates and unreliable standard errors.

**Stock-Yogo critical values** provide formal thresholds for the F statistic. Two types exist:

- **IV size**: the critical value such that the *maximal* rejection rate of a nominal 5% Wald test on the endogenous regressors' coefficients is no more than 10%, 15%, 20%, or 25%. "Maximal" means the worst case over all possible true parameter values. For one endogenous variable and one instrument, the 10% maximal-size critical value is `r sy10`.
- **IV relative bias**: the critical value such that the *maximum* relative bias of 2SLS (relative to OLS bias) is no more than 5%, 10%, 20%, or 30%. Stock and Yogo tabulate these only when the number of overidentifying restrictions is at least two: the measure relies on the finite-sample moments of the 2SLS estimator, which exist only up to the degree of overidentification (Stock & Yogo, 2005; Baum, Schaffer & Stillman, 2007, p. 490). That is why neither this just-identified fit nor the overidentified fit below displays them: both have fewer than two overidentifying restrictions.

Here the Cragg-Donald F of `r round(cd_f, 2)` falls short of the 10% maximal-size critical value of `r sy10` but clears the 15% value of `r sy15` --- a genuinely borderline case. A 5% Wald test on `educ` could reject as often as 15% of the time under the null, so the conventional t-statistic should be read with some suspicion. This is exactly the situation the weak-instrument-robust Anderson-Rubin test (below) is designed for.

One caution: the Stock-Yogo critical values are tabulated for the Cragg-Donald statistic under i.i.d. errors. Under a robust, cluster, or HAC VCE the reported statistic is instead the Kleibergen-Paap rk Wald F, and Baum, Schaffer & Stillman (2007, p. 490) recommend applying the i.i.d. values to it only with caution, or falling back on the Staiger & Stock (1997) rule that F should exceed 10. See `?ivreg2` and `vignette("advanced-iv")` for the fuller discussion.

### Underidentification

The fit above reports an Anderson LM statistic of `r round(underid_stat, 2)` (p = `r format.pval(underid_p, digits = 2, eps = 1e-4)`). The underidentification test asks whether the excluded instruments have *any* correlation with the endogenous regressors. Formally, the null hypothesis is that the matrix of reduced-form coefficients on the excluded instruments has rank K1 - 1, where K1 is the number of endogenous regressors --- i.e., that the equation is underidentified. With a single endogenous regressor, as here, this reduces to the intuitive statement that the instruments are irrelevant. Rejection is good news: it means the model is identified.

With *several* endogenous regressors the test covers the whole system, failing to reject if any one of them is unidentified, so the per-regressor first-stage statistics below are the finer diagnostic. Under i.i.d. errors this is the Anderson (1951) canonical-correlations LM test; under heteroskedasticity or clustering it becomes the Kleibergen-Paap (2006) rk LM statistic.

Underidentification is distinct from weak identification. The underidentification test asks "is there *any* correlation at all?" while the weak identification test asks "is there *enough* correlation for reliable inference?" Here the underidentification test rejects decisively even though the instrument is borderline-weak by the Stock-Yogo standard --- a common pattern, and a reminder that rejecting underidentification is not by itself evidence of strong instruments.

### Overidentification

When there are more excluded instruments than endogenous regressors, the model is overidentified and we can test whether the instruments are valid (i.e., uncorrelated with the structural error). Adding proximity to a *two*-year college (`nearc2`) as a second instrument gives the specification of Computer Exercise C5, Chapter 15, in Wooldridge (2020):

```{r overid}
fit_overid <- ivreg2(
  lwage ~ exper + expersq + black + smsa + south + smsa66 +
    reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + reg669 |
    educ | nearc2 + nearc4,
  data = card
)
summary(fit_overid)
```

```{r overid-values, include = FALSE}
sargan <- diagnostics(fit_overid) |> filter(test == "overid")
```

The Sargan (1958) test (for iid models) or Hansen (1982) J test (for robust/cluster models) tests the null hypothesis that all instruments are valid. Here the Sargan statistic is `r round(sargan$statistic, 2)` (p = `r round(sargan$p_value, 4)`), so the single overidentifying restriction is not rejected. A rejection would have indicated that the instruments are not satisfying the required orthogonality conditions --- either because an instrument is not truly exogenous, or because it is incorrectly excluded from the regression. The test is a joint test of instrument validity *and* correct model specification (Baum, Schaffer & Stillman, 2003, pp. 16--17).

The test exists only when the model is overidentified; with exact identification there are no testable restrictions. It has low power when instruments are weak, and no power at all against violations that leave every instrument consistent for the same wrong parameter value, so a failure to reject does not by itself establish validity (see `?ivreg2` for the fuller treatment).

To test a *subset* of instruments while maintaining the rest, use the `orthog` argument --- see `vignette("advanced-iv")` for details.

### Endogeneity test (C-statistic)

The endogeneity test (also called the C-statistic or difference-in-Sargan test) asks whether the regressors treated as endogenous are actually endogenous. The null hypothesis is that they are exogenous. For the just-identified fit above, the statistic is `r round(endog_stat, 2)` (p = `r round(endog_p, 4)`): we cannot reject that `educ` is exogenous. In that case OLS may be preferred because it is more efficient --- this fit is a live example of the trade-off, since the IV point estimate is similar to OLS but far less precise.

The test maintains instrument validity under both hypotheses, so it is uninformative about endogeneity if the instruments are themselves invalid, and it shares the low power under weak instruments noted above. Using its result to choose between OLS and IV is a pre-test whose selection step can badly distort the size of subsequent inference (Guggenberger, 2010); see `?ivreg2` and Hayashi (2000, pp. 218--221, 232--234) for the underlying theory.

Stata's `ivreg2` computes this test only when `endog()` is requested, whereas `ivreg2r` reports it automatically for all endogenous regressors; use the `endog` argument to test a specific subset (e.g., `endog = "educ"`).

### Anderson-Rubin test

The Anderson-Rubin (1949) test provides inference on the endogenous regressors that remains valid even when instruments are weak. Unlike the standard Wald or t-test on 2SLS coefficients --- which can be severely distorted by weak instruments --- the AR test controls size whether instruments are strong or weak, under the model's maintained error assumptions.

The null hypothesis is that the coefficients on the endogenous regressors are jointly zero *and* the overidentifying restrictions hold. When the model is overidentified, rejection could indicate either nonzero coefficients or invalid overidentifying restrictions (or both). In the just-identified case, the overidentifying restrictions are trivially satisfied and the test is purely about the coefficients.

This serves as a useful cross-check: if the standard t-test says an endogenous regressor is significant but the AR test does not reject, weak instruments may be distorting inference. In the running example the two tests agree but the gap is visible: the conventional z-test on `educ` gives p = `r round(educ_z_p, 4)`, while the Anderson-Rubin chi-squared test gives p = `r round(ar_chi2_p, 4)` --- still a rejection at the 5% level, but a weaker one, consistent with the borderline first stage.

Robustness to weak instruments comes at a price: as instruments weaken, the power of the AR test declines --- Baum, Schaffer & Stillman (2007, p. 491) describe this as an intuitively appealing feature, since weak instruments *should* make confident inference harder. In overidentified models the AR chi-squared statistic has L1 degrees of freedom (L1 = number of excluded instruments), so power also falls as instruments are added without strengthening identification. More powerful weak-instrument-robust procedures of the conditional likelihood ratio type exist for those cases; `vignette("advanced-iv")` discusses what this package does and does not offer and names software that constructs weak-instrument-robust confidence sets.

The Stock-Wright (2000) S statistic tests the same joint null hypothesis using a GMM-distance formulation. It is also automatically computed. See `vignette("advanced-iv")` for more on weak-instrument-robust inference.

## First-stage diagnostics

The `summary()` output reports a table of first-stage diagnostics for each endogenous variable. Fitting with `first_stage = TRUE` additionally stores the underlying first-stage regressions, which `first_stage()` returns as a named list --- one `ivreg2_first_stage` object per endogenous regressor, each supporting `summary()`, `coef()`, `vcov()`, `tidy()`, and `glance()`:

```{r first-stage}
fit_fs <- ivreg2(iv_formula, data = card, first_stage = TRUE)
fs <- first_stage(fit_fs)
summary(fs$educ)
```

The headline number is the partial F-statistic of the excluded instruments, the standard measure of instrument relevance. With one endogenous regressor and one instrument it equals the Cragg-Donald Wald F reported above (`r round(glance(fs$educ)$f_stat, 2)` here), so the borderline weak-identification reading carries over directly. The `tidy()` method returns the same fit as a coefficient tibble for use in tables:

```{r first-stage-tidy}
tidy(fs$educ)
```

The remaining columns are per-regressor diagnostics that become informative only with several endogenous regressors: Shea's (1997) partial R-squared adjusts for intercorrelation among the instruments, while the Sanderson-Windmeijer (2016) and Angrist-Pischke (2009) conditional F statistics assess each endogenous regressor's identification given the others. With a single endogenous regressor the two conditional F statistics collapse to the ordinary first-stage F, and Shea's partial R-squared to the ordinary partial R-squared. See `?ivreg2` for their definitions and `?first_stage` for the stored per-regressor regressions.

## Robust and clustered inference

### Heteroskedasticity-robust standard errors

Use `vcov = "robust"` for White (1980) heteroskedasticity-robust standard errors. Add `small = TRUE` for a finite-sample correction (equivalent to Stata's `robust small`):

```{r robust}
fit_robust <- ivreg2(iv_formula, data = card, vcov = "robust", small = TRUE)
summary(fit_robust)
```

```{r robust-values, include = FALSE}
diag_rob <- diagnostics(fit_robust)
kp_f_rob <- diag_rob |> filter(test == "weak_id_robust") |> pull(statistic)
cd_f_rob <- diag_rob |> filter(test == "weak_id") |> pull(statistic)
```

Notice that the diagnostics automatically adapt. The underidentification test is now the Kleibergen-Paap rk LM statistic instead of the Anderson LM statistic. In the weak-identification block, the Kleibergen-Paap rk Wald F (`r round(kp_f_rob, 2)` here) is reported *alongside* the Cragg-Donald Wald F (`r round(cd_f_rob, 2)`), which remains in the output for reference --- only the rk statistic is valid under non-i.i.d. errors. In overidentified models the Hansen J replaces the Sargan statistic (this fit is exactly identified, so the line reads "(equation exactly identified)").

The `small = TRUE` option applies finite-sample corrections throughout: t-statistics and F-tests replace z-statistics and chi-squared tests.

### Cluster-robust standard errors

When observations are correlated within groups (e.g., repeated observations on the same individual in panel data), use the `clusters` argument. Panel data is the natural home of clustering, so we switch here to the NLS young women's panel (`nlswork`) and regress log wages on grade, age, total labor-market experience, and job tenure, clustering on the individual (`idcode`). The example is ordinary least squares with no instruments; the `clusters` argument behaves identically once the formula gains endogenous and excluded-instrument parts, and the clustered-IV demonstrations live in the other two vignettes.

```{r cluster}
data(nlswork)
fit_cluster <- ivreg2(
  ln_wage ~ grade + age + ttl_exp + tenure,
  data = nlswork, clusters = ~ idcode
)
summary(fit_cluster)
```

Supplying `clusters` selects the cluster-robust VCE automatically, without setting `vcov`. To see how much clustering changes the standard error on `grade`, we compare it against a classical (i.i.d.) fit on the same data. broom's `tidy()` returns a fit's coefficient table as a tibble, standard errors and all; stacking the tidied i.i.d. and clustered fits and keeping the `grade` row puts the two standard errors side by side:

```{r cluster-values}
fit_iid <- ivreg2(ln_wage ~ grade + age + ttl_exp + tenure, data = nlswork)

grade_se <- bind_rows(iid = tidy(fit_iid), cluster = tidy(fit_cluster), .id = "vce") |>
  filter(term == "grade")
grade_se
```

Fitting on the full `r format(nobs(fit_cluster), big.mark = ",")` observations, the cluster-robust standard error on `grade` is `r formatC(grade_se$std.error[grade_se$vce == "cluster"], format = "f", digits = 4)`, versus `r formatC(grade_se$std.error[grade_se$vce == "iid"], format = "f", digits = 4)` under the i.i.d. VCE. Clustering roughly doubles the standard error because wages are serially correlated within person across the panel's multiple survey years, a dependence that classical and heteroskedasticity-robust standard errors both ignore.

Cluster-robust inference is asymptotic in the number of clusters rather than observations: with only a handful of clusters the standard errors can be biased downward even when N is large, and `small = TRUE` only partially compensates. Here the `r format(fit_cluster$n_clusters, big.mark = ",")` clusters, one per individual, put us well inside the asymptotic regime. Researchers clustering on a dozen states or regions should consider wild-cluster bootstrap or other few-cluster methods.

When `small = TRUE` is added, `ivreg2r` applies the standard cluster finite-sample correction factor, `(N-1)/(N-K) * M/(M-1)`, where M is the number of clusters, and reports t- and F-statistics instead of z- and chi-squared statistics.

For two-way clustering (e.g., by firm and year in a panel), use `clusters = ~ var1 + var2`. See `vignette("advanced-iv")` for details.

### VCE option summary

| `vcov` | `clusters` | `small` | Equivalent Stata command |
|--------|-----------|---------|------------------------|
| `"iid"` | — | `FALSE` | `ivreg2 ...` |
| `"iid"` | — | `TRUE` | `ivreg2 ..., small` |
| `"robust"` | — | `FALSE` | `ivreg2 ..., robust` |
| `"robust"` | — | `TRUE` | `ivreg2 ..., robust small` |
| — | `~ var` | `FALSE` | `ivreg2 ..., cluster(var)` |
| — | `~ var` | `TRUE` | `ivreg2 ..., cluster(var) small` |

When `clusters` is supplied, VCE is automatically cluster-robust regardless of the `vcov` argument.

For HAC, AC, Kiefer, and Driscoll-Kraay VCE options, see `vignette("time-series-gmm")`.

## Weighted estimation

`ivreg2r` supports Stata's four weight types through `weight_type`: sampling weights (`"pweight"`), analytic weights (`"aweight"`, the default), frequency weights (`"fweight"`), and importance weights (`"iweight"`). All four are passed with the `weights` argument. We work through two of them --- sampling weights on the genuine survey weight the Card data ship with, and analytic weights on a grouped-data example --- then summarize how the remaining two differ.

### Sampling weights (pweight)

The Card data include a genuine survey weight: `weight` is the NLS sampling weight, reflecting each respondent's probability of selection into the survey. Probability/sampling weights (equivalent to Stata's `[pw=varname]`) are specified with `weights` plus `weight_type = "pweight"`:

```{r weights}
fit_pw <- ivreg2(iv_formula, data = card, weights = weight, weight_type = "pweight",
                 vcov = "robust")
summary(fit_pw)
```

Because pweights describe the sampling design rather than the error variance, valid standard errors must be robust: `ivreg2r` forces a robust VCE whenever `weight_type = "pweight"`, regardless of the `vcov` setting, exactly as Stata does for `[pw=]`. We pass `vcov = "robust"` explicitly above; leaving `vcov` at its default `"iid"` would trigger the same forced override, with a warning noting it. The weighted point estimate answers a population-level question --- what is the return to schooling in the population the NLS sample was drawn from --- whereas the unweighted fit describes the sample.

### Analytic weights (aweight)

An analytic weight is inversely proportional to an observation's variance. The textbook case is grouped data: each row is a cell mean, weighted by the number of individuals that went into the cell, and the defining property is that a weighted regression on the cell means reproduces the coefficients of the individual-level regression. We can watch that identity hold. Collapse the individuals to the cells of three binary covariates, taking each cell's mean log wage and its size --- `group_by() |> summarize()` is Stata's `collapse (mean) lwage (count) n = lwage, by(black smsa south)`:

```{r aweight}
cells <- card |>
  group_by(black, smsa, south) |>
  summarize(lwage = mean(lwage), n = n(), .groups = "drop")

micro   <- ivreg2(lwage ~ black + smsa + south, data = card)
grouped <- ivreg2(lwage ~ black + smsa + south, data = cells,
                  weights = n, weight_type = "aweight")

all.equal(coef(micro), coef(grouped))
```

The cell means are computed in double precision, so nothing is lost in the collapse, and the two coefficient vectors agree to numerical tolerance: `aweight` on the cell means recovers exactly what the individual-level regression would have given.

### The other weight types

The remaining two types are passed the same way --- `weights =` plus `weight_type =` --- and differ only in what the weight *means*:

- **Frequency weights** (`weight_type = "fweight"`, Stata's `[fw=]`) record how many identical observations a row stands for, so a row with weight 3 counts as three copies and the sample size becomes the sum of the weights: `ivreg2(y ~ x, data = d, weights = count, weight_type = "fweight")`. They obey an identity of the same kind as analytic weights: an fweighted fit reproduces the fit on the fully expanded rows, standard errors and sample size included. The [validation article](https://restatr.com/ivreg2r/articles/validation.html) checks that identity against Stata.
- **Importance weights** (`weight_type = "iweight"`, Stata's `[iw=]`) attach an arbitrary multiplier to each observation without a probabilistic interpretation. Like frequency weights they are not normalized, so N becomes the sum of the weights, but they are restricted to the classical i.i.d. VCE and are incompatible with robust, cluster, HAC, and GMM estimation. They exist mainly for building weighted computations on top of `ivreg2r`.

See `?ivreg2` for the precise rules.

For aweights and pweights the weights are normalized internally to sum to N, following Stata, so coefficients, standard errors, and test statistics do not depend on the scale of the weights. One consequence is that for these two types sigma (the Root MSE) differs from `lm()`'s by a factor of `sqrt(N / sum(w))`, because `lm()` uses the raw weights.

## Beyond 2SLS

This vignette focuses on 2SLS, but `ivreg2r` also supports LIML, Fuller, and k-class estimation via the `method`, `fuller`, and `kclass` arguments. LIML has smaller approximate bias than 2SLS when there are many or weak instruments. Fuller's (1977) modified LIML further improves on LIML by having finite moments. These finite-sample advantages are established under i.i.d. errors; neither estimator is robust to violations of that assumption (Baum, Schaffer & Stillman, 2007, pp. 478--479). See `vignette("advanced-iv")` for examples and comparisons.

## Tidyverse workflows

`ivreg2` objects integrate with the tidyverse through `tidy()`, `glance()`, and `augment()`.

### Coefficient table

`tidy()` returns a tibble of coefficient estimates with standard errors, test statistics, p-values, and confidence intervals:

```{r tidy}
tidy(fit_iv)
```

For the intervals alone, `confint(fit_iv)` returns the same bounds as a matrix; they are based on the t-distribution when `small = TRUE` was used at estimation time and on the standard normal otherwise.

### Model summary

`glance()` returns a single-row tibble with fit statistics and the headline diagnostic tests (weak identification, underidentification, overidentification) --- a compact set chosen so that model-comparison tables render sensibly:

```{r glance}
glance(fit_iv)
```

### Diagnostic tests as data

`diagnostics()` returns the full set of specification tests that `summary()` prints --- including the Stock-Yogo critical values --- as a tibble with one row per test:

```{r diagnostics}
diagnostics(fit_iv)
```

The `test` column holds stable keys, so pulling one test out for a table or a threshold comparison is a one-line `filter(test == "weak_id")`. Tests you request through `endog`, `orthog`, or `redundant` --- Stata's `endog()`, `orthog()`, and `redundant()` options --- appear as additional rows.

### Augmented data

`augment()` returns the original data with fitted values and residuals appended:

```{r augment}
augment(fit_iv) |> head()
```

## Stata migration quick reference

#### Estimation basics and accessors

| Task | Stata | R (`ivreg2r`) |
|------|-------|---------------|
| Basic IV | `ivreg2 y x (endo = z)` | `ivreg2(y ~ x | endo | z, data = d)` |
| Extra instruments, no endogenous regressor | `ivreg2 y x (=z1 z2)` | `ivreg2(y ~ x | 0 | z1 + z2, data = d)` (the `0` marks an empty endogenous group) |
| Lag / difference | `L.x` / `D.x` | `l(x, 1)` / `d(x, 1)` (requires `tvar`, plus `ivar` for panels) |
| Coefficients | `e(b)` | `coef(fit)` |
| VCV matrix | `e(V)` | `vcov(fit)` |
| Residuals | `predict resid, resid` | `residuals(fit)` |
| Fitted values | `predict yhat` | `fitted(fit)` |
| Prediction SE | `predict s, stdp` | `predict(fit, se.fit = TRUE)` |
| Confidence intervals | — | `confint(fit)` |
| Tidy coef table | — | `tidy(fit)` |
| Summary statistics | `ereturn list` | `glance(fit)` |
| Diagnostics | Displayed after estimation | Displayed by `summary(fit)`; as a tibble via `diagnostics(fit)` |

#### Estimators

| Task | Stata | R (`ivreg2r`) |
|------|-------|---------------|
| LIML | `ivreg2 ..., liml` | `ivreg2(..., method = "liml")` |
| Fuller(1) | `ivreg2 ..., fuller(1)` | `ivreg2(..., fuller = 1)` |
| k-class | `ivreg2 ..., kclass(0.5)` | `ivreg2(..., kclass = 0.5)` |
| COVIV | `ivreg2 ..., liml coviv` | `ivreg2(..., method = "liml", coviv = TRUE)` |
| GMM2S | `ivreg2 ..., gmm2s` | `ivreg2(..., method = "gmm2s")` |
| CUE | `ivreg2 ..., cue` | `ivreg2(..., method = "cue")` |
| Evaluate CUE at fixed b | `ivreg2 ..., b0(bmat)` | `ivreg2(..., b0 = c(...))` |
| User W matrix | `ivreg2 ..., wmatrix(name)` | `ivreg2(..., wmatrix = W)` |
| User S matrix | `ivreg2 ..., smatrix(name)` | `ivreg2(..., smatrix = S)` |
| Moment covariance / weight matrix | `e(S)` / `e(W)` | `fit$S` / `fit$W` |

#### Variance estimation and clustering

| Task | Stata | R (`ivreg2r`) |
|------|-------|---------------|
| Robust SEs | `ivreg2 ..., robust` | `ivreg2(..., vcov = "robust")` |
| Robust + small | `ivreg2 ..., robust small` | `ivreg2(..., vcov = "robust", small = TRUE)` |
| Cluster SEs | `ivreg2 ..., cluster(g)` | `ivreg2(..., clusters = ~ g)` |
| Cluster + small | `ivreg2 ..., cluster(g) small` | `ivreg2(..., clusters = ~ g, small = TRUE)` |
| Two-way cluster | `ivreg2 ..., cluster(g1 g2)` | `ivreg2(..., clusters = ~ g1 + g2)` |
| AC kernel | `ivreg2 ..., kernel(bar) bw(3)` | `ivreg2(..., kernel = "bartlett", bw = 3, tvar = "year")` |
| HAC kernel | `ivreg2 ..., robust kernel(bar) bw(3)` | `ivreg2(..., vcov = "robust", kernel = "bartlett", bw = 3, tvar = "year")` |
| Auto bandwidth | `ivreg2 ..., kernel(bar) bw(auto)` | `ivreg2(..., kernel = "bartlett", bw = "auto", tvar = "year")` |
| Driscoll-Kraay | `ivreg2 ..., dkraay(3)` | `ivreg2(..., dkraay = 3, tvar = "year", ivar = "id")` |
| Kiefer VCE | `ivreg2 ..., kiefer` | `ivreg2(..., kiefer = TRUE, tvar = "year", ivar = "id")` |
| Stock-Watson VCE | `ivreg2 ..., sw` | `ivreg2(..., sw = TRUE, ivar = "id")` |
| Center moments | `ivreg2 ..., center` | `ivreg2(..., center = TRUE)` |

#### Weights, tests, and other options

| Task | Stata | R (`ivreg2r`) |
|------|-------|---------------|
| Analytic weights | `ivreg2 ... [aw=w]` | `ivreg2(..., weights = w)` |
| Frequency weights | `ivreg2 ... [fw=w]` | `ivreg2(..., weights = w, weight_type = "fweight")` |
| Probability weights | `ivreg2 ... [pw=w]` | `ivreg2(..., weights = w, weight_type = "pweight")` |
| Importance weights | `ivreg2 ... [iw=w]` | `ivreg2(..., weights = w, weight_type = "iweight")` |
| Orthogonality test | `ivreg2 ..., orthog(z1)` | `ivreg2(..., orthog = c("z1"))` |
| Endogeneity test | `ivreg2 ..., endog(z1)` | `ivreg2(..., endog = c("z1"))` |
| Redundancy test | `ivreg2 ..., redundant(z1)` | `ivreg2(..., redundant = c("z1"))` |
| DoF adjustment | `ivreg2 ..., dofminus(5)` | `ivreg2(..., dofminus = 5)` |
| Small-sample DoF | `ivreg2 ..., sdofminus(3)` | `ivreg2(..., sdofminus = 3)` |
| Partialling | `ivreg2 ..., partial(x1)` | `ivreg2(..., partial = c("x1"))` |
| Partial constant / all | `ivreg2 ..., partial(_cons)` | `ivreg2(..., partial = "_cons")` or `"_all"` |
| No partial DoF adjustment | `ivreg2 ..., nopartialsmall` | `ivreg2(..., nopartialsmall = TRUE)` |
| Save first stage | `ivreg2 ..., savefirst` | `ivreg2(..., first_stage = TRUE)`, then `first_stage(fit)` |
| Save reduced form | `ivreg2 ..., saverf` | `ivreg2(..., reduced_form = "rf")` (`rf` displays, `saverf` stores) |
| Save first-stage system | `ivreg2 ..., savesfirst` | `ivreg2(..., reduced_form = "system")` (`savesfirst` is undocumented in Stata) |
| Suppress ID stats | `ivreg2 ..., noid` | `ivreg2(..., noid = TRUE)` |

See `vignette("advanced-iv")` for LIML, k-class, two-way clustering, and model comparison examples. See `vignette("time-series-gmm")` for HAC/AC, GMM, CUE, partialling, and panel VCE examples.

### Intentional differences from Stata

| Area | Stata `ivreg2` | `ivreg2r` | Rationale |
|------|----------------|-----------|-----------|
| Formula | `ivreg2 y x (endo = z)` | `y ~ x | endo | z` | R convention; each variable listed once |
| Numerical method | Cross-products via Mata | QR-based regression via `lm.fit` | Better conditioning for core estimation |
| Diagnostics | Some computed at post-estimation | All computed at estimation time | Returned as a tibble by `diagnostics()`; headline tests also in `glance()` |
| Factor variables | `i.` prefix syntax | R `model.matrix()` / contrasts | R convention |
| Time-series | `tsset` declares time/panel | `tvar=` / `ivar=` arguments | R has no `tsset`; explicit args |
| CUE optimizer | Mata `optimize()` | R `optim()` (BFGS + Nelder-Mead) | Different optimizer, same objective |
| `redundant` with no endogenous regressors | Silently ignored | Warns that the test is skipped | An explicitly requested test should not disappear silently |
| Missing values in a cluster variable | Silently drops the affected observations from the estimation sample | Raises an error asking the user to drop or handle those rows explicitly | The estimation sample should never change silently |

## References

- Anderson, T.W. (1951). "Estimating Linear Restrictions on Regression Coefficients for Multivariate Normal Distributions." *Annals of Mathematical Statistics*, 22(3), 327--351.
- 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.
- Angrist, J.D. and Pischke, J.-S. (2009). *Mostly Harmless Econometrics*. Princeton University Press.
- Baum, C.F., Schaffer, M.E., and Stillman, S. (2003). "Instrumental Variables and GMM: Estimation and Testing." *Stata Journal*, 3(1), 1--31.
- 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.
- Card, D. (1995). "Using Geographic Variation in College Proximity to Estimate the Return to Schooling." In L.N. Christofides, E.K. Grant, and R. Swidinsky (Eds.), *Aspects of Labour Market Behaviour: Essays in Honour of John Vanderkamp*. University of Toronto Press.
- Cragg, J.G. and Donald, S.G. (1993). "Testing Identifiability and Specification in Instrumental Variable Models." *Econometric Theory*, 9(2), 222--240.
- Fuller, W.A. (1977). "Some Properties of a Modification of the Limited Information Estimator." *Econometrica*, 45(4), 939--953.
- Guggenberger, P. (2010). "The Impact of a Hausman Pretest on the Asymptotic Size of a Hypothesis Test." *Econometric Theory*, 26(2), 369--382.
- Hansen, L.P. (1982). "Large Sample Properties of Generalized Method of Moments Estimators." *Econometrica*, 50(4), 1029--1054.
- Hayashi, F. (2000). *Econometrics*. Princeton University Press.
- Kleibergen, F. and Paap, R. (2006). "Generalized Reduced Rank Tests Using the Singular Value Decomposition." *Journal of Econometrics*, 133(1), 97--126.
- Sanderson, E. and Windmeijer, F. (2016). "A Weak Instrument F-Test in Linear IV Models with Multiple Endogenous Variables." *Journal of Econometrics*, 190(2), 212--221.
- Sargan, J.D. (1958). "The Estimation of Economic Relationships Using Instrumental Variables." *Econometrica*, 26(3), 393--415.
- Shea, J. (1997). "Instrument Relevance in Multivariate Linear Models: A Simple Measure." *Review of Economics and Statistics*, 79(2), 348--352.
- Staiger, D. and Stock, J.H. (1997). "Instrumental Variables Regression with Weak Instruments." *Econometrica*, 65(3), 557--586.
- 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.
- White, H. (1980). "A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity." *Econometrica*, 48(4), 817--838.
- Wooldridge, J.M. (2020). *Introductory Econometrics: A Modern Approach*, 7th edition. Cengage.