---
title: "Getting started with estimatr"
output: rmarkdown::html_vignette
bibliography: estimatr.bib
link-citations: yes
vignette: >
  %\VignetteIndexEntry{Getting started with estimatr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE)
options(digits = 4, width = 90)
```

estimatr provides a small set of estimators for design-based inference, with the variance estimators and degrees-of-freedom corrections that social scientists actually use, and it fits them quickly. Base R and most packages default to classical standard errors, which are rarely if ever justified by the research design; getting robust or cluster-robust standard errors usually means fitting a model in one package and correcting it in another. Here the correction is an argument directly in the estimation function.

The package supports six estimators:

| function | what it estimates |
|---|---|
| `lm_robust()` | linear regression with heteroskedasticity-robust or cluster-robust standard errors |
| `lm_lin()` | a treatment effect adjusted for pre-treatment covariates, following @lin2013 |
| `iv_robust()` | two-stage least squares |
| `difference_in_means()` | a difference in means, with the variance that matches the randomization |
| `horvitz_thompson()` | an average treatment effect by inverse probability weighting |
| `lh_robust()` | a linear combination of coefficients, or several tested jointly |

This vignette shows each of them on one running example. `vignette("mathematical-notes")` gives the definitions and the citations; `vignette("estimatr2.0")` is for readers coming from estimatr 1.x, and lists what changed.

## The example data

One hundred units in twenty clusters of five, with a pre-treatment covariate `X`, a sampling weight `W`, ten blocks formed from `X`, and two potential outcomes. The same experiment gets run three ways (complete, clustered, and blocked random assignment), and a fourth time with noncompliance for `iv_robust()`. Each assignment is made in the section that uses it. The assignments come from randomizr, whose declarations `horvitz_thompson()` reads later on.

```{r, warning = FALSE}
library(estimatr)
library(randomizr)
library(dplyr)
```

```{r}
set.seed(343)
N <- 100

dat <- tibble(
  X = runif(N),
  cluster = rep(letters[1:20], each = 5),
  block = cut(X, breaks = quantile(X, seq(0, 1, 0.1)), labels = FALSE, include.lowest = TRUE),
  W = runif(N, 0.5, 1.5), # sampling weights, used further down
  Y_Z_0 = rnorm(N, mean = X) + rep(rnorm(20, sd = 0.5), each = 5), # control potential outcome, with a shock shared within each cluster
  Y_Z_1 = Y_Z_0 + 0.35 # treatment potential outcome (constant effects)
)

print(dat, n = 5)
```

## `lm_robust()`

Usage follows `lm()`. The difference is the default standard error: HC2 rather than classical. Start with complete random assignment of half the units, and reveal the potential outcome that corresponds to each unit's assignment.

```{r}
dat <- dat |>
  mutate(
    Z = complete_ra(N, m = N / 2),
    Y = if_else(Z == 1, Y_Z_1, Y_Z_0)
  )

fit <- lm_robust(Y ~ Z + X, data = dat)
fit
```

`summary()` adds the R-squared and the F statistic; `tidy()` returns the coefficient table as a tibble, which is usually what you want next.

```{r}
tidy(fit)
```

### Clustered standard errors

Now assign treatment by cluster instead: ten of the twenty clusters, with every unit in a cluster sharing its cluster's assignment.

```{r}
dat <- dat |>
  mutate(
    Z_cl = cluster_ra(clusters = cluster, m = 10),
    Y_cl = if_else(Z_cl == 1, Y_Z_1, Y_Z_0)
  )
```

Pass the name of the cluster variable. The default becomes CR2, the cluster analogue of HC2, using the small-sample correction of @pustejovskytipton2018.

```{r}
# Ignoring the clustered assignment understates the standard error
fit_naive <-
  lm_robust(Y_cl ~ Z_cl + X, data = dat)
tidy(fit_naive)

# Accounting for it
fit_cluster_aware <-
  lm_robust(Y_cl ~ Z_cl + X, data = dat, clusters = cluster)
tidy(fit_cluster_aware)
```

Notice the degrees of freedom as well as the standard error. Under CR2 they are computed per coefficient by a Satterthwaite approximation, so they are not the residual degrees of freedom and they differ across the terms.

### Choosing a different variance estimator

Without clusters, `se_type` takes `"classical"`, `"HC0"`, `"HC1"`, `"HC2"`, and `"HC3"`; with clusters, `"CR0"` and `"CR2"`. `"stata"` is accepted in both cases and reproduces what Stata's `robust` or `cluster` option gives, which is HC1 without clusters and CR0 with Stata's two small-sample corrections with them.

```{r}
fit_HC2 <- lm_robust(Y ~ Z + X, data = dat, se_type = "HC2")

fit_classical <- lm_robust(Y ~ Z + X, data = dat, se_type = "classical")

fit_stata <- lm_robust(Y ~ Z + X, data = dat, se_type = "stata")

# The standard error on Z under each
c(
  HC2 = fit_HC2$std.error[["Z"]],
  classical = fit_classical$std.error[["Z"]],
  stata = fit_stata$std.error[["Z"]]
)
```

### Fixed effects

We can demonstrate how to include fixed effects in an lm_robust call with a block-randomized experiment:

```{r}
dat <- dat |>
  mutate(
    Z_bl = block_ra(blocks = block, prob = 0.5),
    Y_bl = if_else(Z_bl == 1, Y_Z_1, Y_Z_0)
  )
```

`fixed_effects` takes a right-sided formula and absorbs those groups rather than adding them as dummy columns. The estimates and standard errors are identical to writing the dummies out. Because `lm_robust` absorbs rather than building the dummy matrix, the computation is very fast, even with many fixed effects.

```{r}
lm_robust(Y_bl ~ Z_bl, data = dat, fixed_effects = ~ block)

# The same numbers, using the slower dummy matrix way.
lm_robust(Y_bl ~ Z_bl + factor(block), data = dat)
```

The absorbed levels are kept in the `felevels` element of the result. One caution: `fixed_effects` combined with `clusters` defaults to `se_type = "CR0"`, because CR2 is the one estimator that still has to expand the dummies, and it warns once per session to say so. Asking for `se_type = "CR2"` by name gets it. The full menu by case is in `?lm_robust`.

### Weights

`weights` takes a column of the data, and the fit is weighted least squares with the matching robust variance.

```{r}
weighted_fit <- lm_robust(Y ~ Z + X, data = dat, weights = W)
tidy(weighted_fit)
```

Weighted HC2 and HC3 differ from Stata's, because the two use different definitions of the hat matrix. `vignette("mathematical-notes")` sets out the difference and why estimatr follows the R convention.

## `lm_lin()`

Adjusting for pre-treatment covariates in an experimental context can increase precision, but @freedman2008 showed that doing it with an ordinary regression can bias the estimated treatment effect and can in some cases *reduce* precision. @lin2013 proposed a fix for the precision problem: center every covariate, interact each with treatment, and regress the outcome on treatment, the centered covariates, and the interactions.

`lm_lin()` does the pre-processing for you. It takes the treatment in `formula` and the covariates in a right-sided `covariates` formula; everything else is `lm_robust()`'s.

```{r}
lin_fit <- lm_lin(Y ~ Z, covariates = ~ X, data = dat)
tidy(lin_fit)
```

The `Z` row is the estimate of the average treatment effect, since when the covariates are at their centered means of zero, the interaction terms drop out. The centers used are returned in `scaled_center`, and centering happens after any function in the formula is evaluated, so `~ log(X)` centers the log.

```{r}
lin_fit$scaled_center
```

## `iv_robust()`

Instrumental variables estimation comes up most often in experiments with noncompliance. The random assignment `Z` says who was encouraged to take the treatment, and `D` records who actually took it. Here 60 percent of units are compliers, who take the treatment if and only if assigned; the rest never take it, whatever their assignment (one-sided noncompliance). Outcomes follow treatment receipt, not assignment.

```{r}
dat <- dat |>
  mutate(
    Y_D_0 = Y_Z_0, # potential outcome if treatment is not received
    Y_D_1 = Y_D_0 + 0.35, # potential outcome if treatment is received
    complier = rbinom(N, size = 1, prob = 0.6),
    D = Z * complier,
    Y_nc = if_else(D == 1, Y_D_1, Y_D_0)
  )
```

Because `Z` is randomized, it is an instrument for `D`, and two-stage least squares estimates the complier average causal effect [@gerbergreen2012], which is 0.35 here because the effect is constant. `iv_robust()` takes the same `se_type` menu as `lm_robust()`. The instruments go after a vertical bar, and a covariate goes on both sides of it.

```{r}
iv_fit <- iv_robust(Y_nc ~ D + X | Z + X, data = dat)
summary(iv_fit)
```

`diagnostics = TRUE` adds the weak-instrument, Wu-Hausman, and overidentification tests; the last is Sargan's under `se_type = "classical"` and Wooldridge's robust score test otherwise. `residuals()` returns the structural residuals, $y - X\widehat{\beta}$, rather than the second-stage ones.

## `difference_in_means()`

A difference in means is simple until the design is blocked or clustered, at which point the point estimate has to average over blocks and the variance has to reflect how units were actually assigned. `difference_in_means()` picks the estimator that matches the design and tells you which one it picked.

```{r}
dim_fit <- difference_in_means(Y ~ Z, data = dat)
dim_fit
dim_fit$design
```

With no blocks or clusters, the variance and degrees of freedom are what `t.test()` computes. Pass `clusters` and the design changes, and so does the estimator:

```{r}
dim_cl <- difference_in_means(Y_cl ~ Z_cl, data = dat, clusters = cluster)
dim_cl$design
tidy(dim_cl)
```

`blocks` works the same way, and blocks need not all be the same shape. Blocks with at least two units in each arm carry their own variance; blocks with a singleton arm have their variance estimated across such blocks, following @pashleymiratrix2021. Check `design` rather than assuming, and see `?difference_in_means` for the full classification and the two designs it refuses.

```{r}
# The blocked version of the experiment
dim_bl <- difference_in_means(Y_bl ~ Z_bl, data = dat, blocks = block)
dim_bl$design
tidy(dim_bl)

# R's sleep data, ten patients each measured under two drugs, is a
# matched-pairs design, and is recognised as such
difference_in_means(extra ~ group, data = sleep, blocks = ID)$design
```

## `horvitz_thompson()`

When you know how treatment was assigned, you can estimate the average treatment effect without a model, by weighting each observed outcome by the inverse of the probability of the condition it was observed under. The estimator is unbiased, and it handles designs that regression handles awkwardly: clusters of unequal size, per-unit probabilities that differ, assignment schemes with dependence across units.

`condition_prs` is the one argument that describes the design, and what you pass decides which variance you get. A named vector of condition probabilities gives the conservative bound, valid for any design and tight only for Bernoulli assignment:

```{r}
horvitz_thompson(Y ~ Z, data = dat, condition_prs = c("0" = 0.5, "1" = 0.5))
```

Passing the randomization declaration instead is what buys the design-aware variance, because a declaration already carries the block structure, the cluster structure, the per-unit probabilities and whether assignment was simple or complete. Clustered and blocked designs need no further arguments: describe the design once, in the declaration.

```{r}
# declare_ra is from the randomizr package
decl <- declare_ra(N = N, prob = 0.5, simple = FALSE)
horvitz_thompson(Y ~ Z, data = dat, condition_prs = decl)
```

`?horvitz_thompson` covers the rest: per-unit probability matrices, designs supplied as a permutation matrix, and contrasting two arms of a multi-arm design.

## `lh_robust()`

To test a linear combination of coefficients, or several at once, `lh_robust()` fits the model and applies `car::linearHypothesis()` to it, keeping the robust variance and the degrees of freedom of the fit.

```{r}
lh_fit <- lh_robust(Y ~ Z + X, data = dat, linear_hypothesis = "Z + 2*X = 0")
tidy(lh_fit)
```

Several restrictions give a joint Wald test as well.

```{r}
joint <- lh_robust(Y ~ Z + X, data = dat, linear_hypothesis = c("Z = 0", "X = 0"))
joint$joint_hypothesis
```

## Working with the output

Every estimator returns an object supporting the S3 methods you would expect: `summary()`, `print()`, `tidy()`, `glance()`, `coef()`, `confint()`, `vcov()`, `nobs()`, and, for the regression estimators, `predict()`, `residuals()` and `update()`.

```{r}
glance(fit)
confint(fit)
```

Regression tables work through the usual packages. `texreg` and `modelsummary` both take `lm_robust` objects directly.

```{r, results = "asis"}
texreg::htmlreg(list(fit, lin_fit), include.ci = FALSE, caption = "")
```

`emmeans` also has methods for `lm_robust` fits, for estimated marginal means and contrasts, computed with the fit's robust variance.

```{r, eval = requireNamespace("emmeans", quietly = TRUE)}
emmeans::emmeans(fit, "Z", at = list(Z = 0:1))
```

## In a simulation

Design-based work means fitting the same model many times over, which is what estimatr is built to do quickly. [DeclareDesign](https://declaredesign.org) handles the repetition: declare the model, the inquiry, the assignment, and the estimator once, and `diagnose_design()` re-runs the whole experiment and reports how the estimator does against the truth.

```{r, eval = requireNamespace("DeclareDesign", quietly = TRUE) && requireNamespace("future.apply", quietly = TRUE)}
library(DeclareDesign)

declaration_lin <-
  declare_model(N = 100, 
                X = runif(N), 
                U = rnorm(N, mean = X), 
                potential_outcomes(Y ~ 0.35 * Z + U)) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(Z = complete_ra(N)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, covariates = ~ X, 
                    .method = lm_lin, inquiry = "ATE")

set.seed(343)
diagnose_design(declaration_lin, sims = 500)
```

Adding a second `declare_estimator()` step, say `declare_estimator(Y ~ Z, .method = lm_robust, inquiry = "ATE", label = "unadjusted")`, is how you compare two estimators on one design.

## Where to go next

- `vignette("mathematical-notes")` for the definition of every estimator and variance above, each with its citation and with a live check that estimatr computes what the definition says.
- `vignette("estimatr2.0")` if you are porting code from estimatr 1.x.
- The [Performance](https://declaredesign.org/r/estimatr/articles/performance.html) page for what each estimator costs at scale.
- The [tidyverse](https://declaredesign.org/r/estimatr/articles/estimatr-in-the-tidyverse.html) and [regression table](https://declaredesign.org/r/estimatr/articles/regression-tables.html) pages for what to do with a fit once you have one.

# References