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

## ----first-fit----------------------------------------------------------------
library(deli)

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

coef(fit)

## ----accessors----------------------------------------------------------------
vcov(fit)

confint(fit)

## ----fit-summary--------------------------------------------------------------
summary(fit)

## ----tidiers------------------------------------------------------------------
tidy(fit, conf.int = TRUE)

glance(fit)

## ----inference-helpers--------------------------------------------------------
z_scores(fit)
p_values(fit)

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

summary(fit_logistic)

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

## ----mean-check---------------------------------------------------------------
mean(y)

## ----estimator-slots----------------------------------------------------------
m@theta       # Point estimate (the mean)
m@variance    # Sandwich variance estimate

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

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

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