---
title: "Regression Models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Regression Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(deli)
```

## Overview

deli provides estimating equations for a wide range of regression models.
Every model in this vignette is fitted the same way: pass a formula, a data
frame, and the estimating equation to `m_estimate()`, which builds the design
matrix, constructs the estimator, and solves it in one call. Arguments belonging
to the estimating equation, such as `model` or `penalty`, go in the same call
and are forwarded to it, while `init`, `subset`, `solver`, and the other solver
controls are arguments of `m_estimate()` itself. The sandwich variance estimator
automatically provides robust standard errors.

The parameters are labeled from the first of three sources that names every one
of them: the names on `init`, the columns of the design matrix, and the row
names the estimating equation writes on its own return. A default `init` carries
no names, so the design columns label the fit, together with any parameter the
equation estimates on top of the coefficients; an explicit `init` is labeled the
same way when its length accounts for one of those two shapes. Where it accounts
for neither, the equation's row names label the fit instead, and where nothing
labels every parameter they are numbered `theta_1` through `theta_p`. Name the
elements of `init` yourself to label the parameters any other way.

`init` defaults to a zero vector with one element per design matrix column. That
is the right length for most estimating equations, but not for the few that
estimate a parameter of their own on top of the coefficients; the gamma and
negative binomial GLMs below are the cases you are likely to meet.

Every model shown here can also be fitted through the function interface, where
you write a `psi` function that supplies the design matrix and response
yourself. That form is what you need for a design the formula notation cannot
express, for a response that is not a vector (`ee_mlogit()`), and for custom or
stacked equations; see `vignette("custom-estimating-equations")`.

## Linear regression

The most basic regression model uses `ee_regression()` with `model = "linear"`:

```{r 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 regression

For binary outcomes, use `model = "logistic"`:

```{r 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 regression

For count data, use `model = "poisson"`:

```{r 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: generalized linear models

`ee_glm()` provides a more flexible interface where you specify the
distribution and link function separately:

```{r 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
```

Available distributions: `"normal"`, `"binomial"`, `"poisson"`, `"gamma"`,
`"negative_binomial"`, `"inverse_gaussian"`, and `"tweedie"`. `"tweedie"` takes
a variance-power `hyperparameter`; see `?ee_glm` for the details. Available
links: `"identity"`, `"log"`, `"logit"`, `"probit"`, `"cauchy"` (alias
`"cauchit"`), `"loglog"`, `"cloglog"`, `"inverse"`, and `"sqrt"`.

`"gamma"` and `"negative_binomial"` estimate one parameter beyond the
regression coefficients: the log of the gamma shape, or the log of the negative
binomial dispersion. The automatic `init` has one element per design matrix
column and so is one element short for these two, which makes an explicit `init`
necessary. It needs no names of its own: the formula interface labels an
unnamed `init` of that length from the model matrix columns and the extra
parameter.

```{r 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
```

The extra parameter is on the log scale, so `exp(m@theta[["log_shape"]])`
recovers the shape, here close to the value of 2 used in the simulation.

## Penalized regression

deli supports several penalized regression methods. These add a penalty term
to the estimating equations.

### Ridge regression

L2 penalty shrinks coefficients toward zero:

```{r 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 regression

L1 penalty produces sparse solutions:

```{r 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
```

### Elastic net

Combines L1 and L2 penalties:

```{r 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 regression

`ee_robust_regression()` replaces the squared loss with a robust loss
function, providing resistance to outliers:

```{r 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)
```

Available loss functions: `"huber"`, `"tukey"`, `"andrew"`, `"hampel"`.

## Weighted regression

All regression EEs support observation weights via the `weights` argument.
Arguments passed through `m_estimate()` are evaluated in the data frame, so the
column name is enough:

```{r 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 after regression

After fitting any regression 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 <- 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)))
```

Called without `newdata`, `augment()` reports the rows the model was fitted to
and adds a `.resid` column as well. For a model with a non-identity link,
`type.predict = "response"` moves `.fitted` and its interval from the linear
predictor to the scale of the response.