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

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

## ----linear-------------------------------------------------------------------
set.seed(42)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
y <- 1 + 2 * x1 - 0.5 * x2 + rnorm(n)
d <- data.frame(x1, x2, y)

m <- m_estimate(y ~ x1 + x2, data = d, .ee = ee_regression, model = "linear")
summary(m)

## ----logistic-----------------------------------------------------------------
set.seed(42)
n <- 500
x <- rnorm(n)
y <- rbinom(n, 1, plogis(0.5 + x))
d <- data.frame(x, y)

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

## ----poisson------------------------------------------------------------------
set.seed(42)
n <- 500
x <- rnorm(n)
y <- rpois(n, lambda = exp(0.5 + 0.3 * x))
d <- data.frame(x, y)

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

## ----glm-poisson--------------------------------------------------------------
set.seed(42)
n <- 500
x <- rnorm(n)
y <- rpois(n, lambda = exp(0.5 + 0.3 * x))
d <- data.frame(x, y)

m <- m_estimate(
  y ~ x,
  data = d,
  .ee = ee_glm,
  distribution = "poisson",
  link = "log"
)
m@theta

## ----glm-gamma----------------------------------------------------------------
set.seed(42)
n <- 500
x <- rnorm(n)
mu <- exp(0.5 + 0.3 * x)
y <- rgamma(n, shape = 2, scale = mu / 2)
d <- data.frame(x, y)

m <- m_estimate(
  y ~ x,
  data = d,
  .ee = ee_glm,
  distribution = "gamma",
  link = "log",
  init = c(0, 0, 0)
)
m@theta

## ----ridge--------------------------------------------------------------------
set.seed(42)
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1 + 0.5 * x1 + 0.3 * x2 + rnorm(n)
d <- data.frame(x1, x2, y)

m <- m_estimate(
  y ~ x1 + x2,
  data = d,
  .ee = ee_ridge_regression,
  model = "linear",
  penalty = 0.5
)

summary(m)

## ----lasso--------------------------------------------------------------------
# Not differentiable, so the sandwich variance should not be trusted here
m <- m_estimate(
  y ~ x1 + x2,
  data = d,
  .ee = ee_lasso_regression,
  model = "linear",
  penalty = 0.1
)
m@theta

## ----elasticnet---------------------------------------------------------------
# The L1 half is not differentiable, so again distrust the sandwich variance
m <- m_estimate(
  y ~ x1 + x2,
  data = d,
  .ee = ee_elasticnet_regression,
  model = "linear",
  penalty = 0.1,
  ratio = 0.5
)
m@theta

## ----robust-------------------------------------------------------------------
set.seed(42)
n <- 200
x <- rnorm(n)
y <- 1 + 2 * x + rnorm(n)
# Add some outliers
y[1:5] <- y[1:5] + 20
d <- data.frame(x, y)

# The Huber loss is convex, so its estimating function has a single root. What
# makes the seed necessary is that the Huber psi is bounded: far from the
# solution every residual is past the tuning constant k, every contribution
# saturates at k, and the estimating function is constant with a Jacobian of
# exactly zero. Starting from zero lands in that flat region, where the solver
# has no slope to follow, so start from a least-squares fit.
start <- coef(lm(y ~ x, data = d))

# Huber loss with k = 1.345
m <- m_estimate(
  y ~ x,
  data = d,
  .ee = ee_robust_regression,
  model = "linear",
  loss = "huber",
  k = 1.345,
  init = start
)

# Compare with OLS (affected by outliers)
m_ols <- m_estimate(y ~ x, data = d, .ee = ee_regression, model = "linear")

rbind(robust = m@theta, ols = m_ols@theta)

## ----weighted-----------------------------------------------------------------
set.seed(42)
n <- 200
x <- rnorm(n)
y <- 1 + 2 * x + rnorm(n)
w <- runif(n, 0.5, 1.5)
d <- data.frame(x, y, w)

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

## ----predictions--------------------------------------------------------------
set.seed(42)
n <- 300
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 covariate values
augment(m, newdata = data.frame(x = seq(-2, 2, length.out = 5)))