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

## ----setup--------------------------------------------------------------------
library(deli)

## ----ratio-of-means-----------------------------------------------------------
set.seed(42)
n <- 200
y1 <- rnorm(n, mean = 4, sd = 1)
y2 <- rnorm(n, mean = 2, sd = 1)

psi <- function(theta) {
  # Row 1: estimating equation for the mean of y1
  mu1 <- y1 - theta[1]
  # Row 2: estimating equation for the mean of y2
  mu2 <- y2 - theta[2]
  # Row 3: estimating equation for the ratio (repeated n times)
  ratio <- rep(theta[1] / theta[2] - theta[3], n)
  # Stack into a 3-by-n matrix. `rbind()` labels each row with the name of the
  # variable it came from, and those labels become the parameter names.
  rbind(mu1, mu2, ratio)
}

m <- m_estimate(stacked_equations = psi, init = c(1, 1, 1))
m@theta

## ----ratio-summary------------------------------------------------------------
summary(m)

## ----naive-ratio--------------------------------------------------------------
mean(y1) / mean(y2)

## ----stacked-odds-ratio-------------------------------------------------------
set.seed(42)
n <- 500
x <- rnorm(n)
pr <- plogis(-0.5 + 0.8 * x)
y <- rbinom(n, 1, pr)
X <- cbind(1, x)

psi <- function(theta) {
  # theta[1:2]: logistic regression coefficients (intercept, slope)
  # theta[3]: odds ratio = exp(theta[2])
  beta <- theta[1:2]
  or <- theta[3]

  # Built-in logistic regression EE (returns a 2-by-n matrix)
  ee_reg <- ee_regression(beta, X = X, y = y, model = "logistic")

  # Custom EE for the odds ratio (deterministic, repeated n times)
  ee_or <- matrix(rep(exp(theta[2]) - or, n), nrow = 1)

  # Stack: 3-by-n matrix
  rbind(ee_reg, ee_or)
}

m <- m_estimate(stacked_equations = psi, init = c(0, 0, 1))
summary(m)

## ----odds-ratio-manual--------------------------------------------------------
exp(m@theta[2])

## ----logistic-fit-------------------------------------------------------------
# Fit the logistic regression only
d <- data.frame(x, y)

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

## ----delta-method-------------------------------------------------------------
# Transform: exponentiate the second coefficient to get the odds ratio
dm_var <- delta_method(m_reg, transform = function(theta) exp(theta[2]))
dm_var

## ----delta-method-ci----------------------------------------------------------
or_est <- exp(m_reg@theta[2])
or_se <- sqrt(dm_var[1, 1])
ci <- or_est + c(-1, 1) * qnorm(0.975) * or_se
round(c(or = unname(or_est), lower = ci[1], upper = ci[2]), 3)

## ----delta-method-raw---------------------------------------------------------
delta_method(
  m_reg@theta,
  transform = function(theta) exp(theta[2]),
  covariance = m_reg@variance
)