---
title: "Getting Started with deli"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with deli}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

deli is an R package for M-estimation and empirical sandwich variance
estimation. It provides:

1. `m_estimate()`, a one-step fitting function with a formula interface for
   the pre-built estimating equations and a function interface for your own
2. Pre-built estimating equations for common models (means, regression, causal
   inference, survival analysis, and more)
3. Automatic sandwich variance estimation
4. The `MEstimator` and `GMMEstimator` classes, which solve arbitrary
   user-specified estimating equations and hold the results

deli is an R port of the Python
[Delicatessen](https://github.com/pzivich/Delicatessen) library.

## Fitting your first model

The quickest way to fit a model is `m_estimate()`, which takes a formula, a
data frame, and a pre-built estimating equation. It constructs the estimator
and solves it in a single call, much as `lm()` fits a linear model in one
step. Here is a linear regression on the `mtcars` data:

```{r first-fit}
library(deli)

fit <- m_estimate(
  mpg ~ wt + hp,
  data = mtcars,
  .ee = ee_regression,
  model = "linear"
)

coef(fit)
```

The `.ee` argument names the estimating equation, and any further arguments
(such as `model = "linear"`) are passed on to it. The formula and data supply
the response and design matrix, so you do not need to build them by hand.

## Inspecting the fit

Fitted estimators support the standard R accessors. `coef()` returns the point
estimates, `vcov()` returns the sandwich variance-covariance matrix, and
`confint()` returns Wald confidence intervals:

```{r accessors}
vcov(fit)

confint(fit)
```

`summary()` collects the estimates, standard errors, test statistics, and
confidence intervals into a single table:

```{r fit-summary}
summary(fit)
```

deli also provides [broom](https://broom.tidymodels.org/) tidiers, so results
flow into tidyverse pipelines. `tidy()` returns one row per parameter, and
`glance()` returns a one-row model summary:

```{r tidiers}
tidy(fit, conf.int = TRUE)

glance(fit)
```

If you prefer deli's own inference helpers, `confidence_intervals()`,
`z_scores()`, `p_values()`, and `s_values()` each return the corresponding
quantities directly:

```{r inference-helpers}
z_scores(fit)
p_values(fit)
```

## Other regression models

The pre-built regression equation covers several model families. Switch from
linear to logistic regression by changing the `model` argument. Here `vs` is a
binary indicator of engine type:

```{r logistic-fit}
fit_logistic <- m_estimate(
  vs ~ mpg + wt,
  data = mtcars,
  .ee = ee_regression,
  model = "logistic"
)

summary(fit_logistic)
```

The sandwich variance estimator supplies robust standard errors automatically,
without any additional arguments. `ee_regression()` covers linear, logistic,
and Poisson models. See `?ee_glm` for the wider family of generalized linear
models, and the dedicated functions `?ee_robust_regression`,
`?ee_ridge_regression`, `?ee_lasso_regression`, `?ee_elasticnet_regression`,
`?ee_beta_regression`, and `?ee_tobit` for robust, ridge, LASSO, elastic net,
beta, and tobit regression.

## The estimator objects

`m_estimate()` builds and solves an `MEstimator` for you. When a pre-built
estimating equation and the formula interface do not fit your problem, for
example when you supply a custom estimating equation, give `m_estimate()` a
function in place of the formula. The workflow has two steps:

1. Define a `psi` function that returns the estimating equation contributions
2. Pass it to `m_estimate()` with starting values for the parameters

Here is a simple example that estimates a mean:

```{r mean-estimator}
y <- c(1, 2, 3, 1, 4, 5, 3, 2, 6, 7)

# Define the estimating equation
psi <- function(theta) {
  ee_mean(theta, y = y)
}

m <- m_estimate(psi, init = 0)

coef(m)
```

The estimated mean matches `mean(y)`:

```{r mean-check}
mean(y)
```

The result is an `MEstimator` object, the same class the formula interface
returns, so the same accessors work here. The point estimate and variance are
also available on the object as `m@theta` and `m@variance`:

```{r estimator-slots}
m@theta       # Point estimate (the mean)
m@variance    # Sandwich variance estimate
```

## Stacking estimating equations

A key strength of M-estimation is the ability to stack estimating equations.
When you stack equations, the sandwich variance correctly accounts for all
sources of uncertainty. Here is an example that estimates a mean and variance
simultaneously:

```{r mean-variance}
y <- c(1, 2, 3, 1, 4, 5, 3, 2, 6, 7)

psi <- function(theta) {
  ee_mean_variance(theta, y = y)
}

m <- m_estimate(stacked_equations = psi, init = c(0, 0))
m@theta  # c(mean, variance)
```

For more complex stacking with custom equations, see
`vignette("custom-estimating-equations")`.

## Causal inference: inverse probability weighting

deli includes estimating equations for causal inference. Because they combine
several components, they are a natural fit for the function interface. Here is
an example that uses inverse probability weighting (IPW) to estimate an average
treatment effect:

```{r ipw}
set.seed(42)
n <- 500
w <- rbinom(n, 1, 0.5)               # Binary confounder
A <- rbinom(n, 1, plogis(-0.5 + w))  # Treatment depends on w
Y <- 1 + 2 * A + w + rnorm(n)        # Outcome
W <- cbind(1, w)                     # Propensity score design matrix

psi <- function(theta) {
  ee_ipw(theta, y = Y, A = A, W = W)
}

# theta: ACE, E[Y(1)], E[Y(0)], beta0, beta1
m <- m_estimate(stacked_equations = psi, init = c(0, 0, 0, 0, 0))

# ACE (average causal effect) ~ 2
m
```

## Predictions

After fitting a model, use `augment()` for predicted values with confidence
intervals. Give it a data frame of new covariate values as `newdata`, with one
column for each covariate the formula names, and it returns that frame with
`.fitted`, `.se.fit`, `.lower`, and `.upper` beside it:

```{r predictions}
set.seed(42)
n <- 200
x <- rnorm(n)
y <- 1 + 2 * x + rnorm(n)
d <- data.frame(x, y)

m <- m_estimate(y ~ x, data = d, .ee = ee_regression, model = "linear")

# Predict at new values
augment(m, newdata = data.frame(x = seq(-2, 2, by = 1)))
```

Called without `newdata`, `augment()` reports the rows the model was fitted to
and adds a `.resid` column as well.

## Clustered data

For clustered or grouped data, use `aggregate_efuncs()` inside your `psi`
function to get cluster-robust variance estimates:

```{r clustered}
set.seed(42)
n <- 200
n_groups <- 50
group <- rep(1:n_groups, each = n / n_groups)
group_effect <- rnorm(n_groups, sd = 2)
y <- group_effect[group] + rnorm(n)

# Cluster-robust variance
psi <- function(theta) {
  ef <- ee_mean(theta, y = y)
  aggregate_efuncs(ef, group = group)
}

m <- m_estimate(stacked_equations = psi, init = mean(y))

m@theta
m@variance  # Accounts for within-cluster correlation
```

## Exact differentiation

The delta method needs the derivative (Jacobian) of a transform, and the
sandwich variance needs the derivative of the estimating equations. By default
deli computes these with central finite differences (`deriv_method =
"capprox"`), which introduces a small step-size approximation. Passing
`deriv_method = "exact"` switches to forward-mode automatic differentiation,
which evaluates the derivative in closed form. It returns exact derivatives with
no step size to tune. In Python Delicatessen exact differentiation is the
default only for `delta_method`; the estimators and the sandwich variance
default to a forward finite difference.

The `delta_method()` function accepts `deriv_method` directly, and the survival
prediction helpers `survival_predictions()` and `aft_predictions_function()`
forward it to the delta method. The following fits a logistic model and uses
exact differentiation to obtain the variance of the fitted probability at a
covariate pattern:

```{r exact-delta-method}
set.seed(42)
n <- 500
x <- rnorm(n)
pr <- plogis(0.5 + x)
y <- rbinom(n, 1, pr)
d <- data.frame(x, y)

m <- m_estimate(y ~ x, data = d, .ee = ee_regression, model = "logistic")

# Variance of the predicted probability at x = 1, via the delta method.
# The log-odds at the pattern (intercept = 1, x = 1) is theta[1] + theta[2].
transform <- function(theta) inverse_logit(theta[1] + theta[2])
delta_method(m, transform = transform, deriv_method = "exact")
```

That chunk computes the inverse logit twice with two different functions, and
the difference is deliberate. The simulation uses `stats::plogis()`, because it
is ordinary numeric work on plain doubles. The transform uses deli's
`inverse_logit()`, because `plogis()` cannot be differentiated exactly and would
stop the delta method with an error. The two functions return identical numbers;
only one of them survives `deriv_method = "exact"`.

## Next steps
- `vignette("custom-estimating-equations")`: Write your own estimating
  equations and stack them
- `?m_estimate`: the one-step formula and function interface
- `?ee_regression`: linear, logistic, and Poisson regression
- `?ee_glm`: generalized linear models, plus `?ee_ridge_regression`,
  `?ee_lasso_regression`, and the other penalized and robust regression
  functions
- `?ee_gformula`: G-computation / standardization
- `?ee_aipw`: Augmented inverse probability weighting
- `?ee_aft`: Accelerated failure time models
- `?ee_survival_model`: Parametric survival models