---
title: "Inference after a glmnet fit: carving for logistic regression"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Inference after a glmnet fit}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      eval = requireNamespace("glmnet", quietly = TRUE))
```

This vignette shows **carving** in logistic regression. We select variables with glmnet on a
random half of the data, then do inference with all of it, accounting for selection.
`glmnet_inference` does this from the fitted `glmnet` object. It uses
`glmnet_problem_from_fit`, which returns the LASSO problem glmnet solved in the original
coordinates: the solution, the gradient and Hessian of the smooth part, the penalty weights
and the bounds. For inference without held-out data see
`vignette("no_randomization", package = "lassoinf")`.

## Carving as randomization

glmnet divides the loss by the number of observations, and we use this per-observation scale
throughout. On the subsample of size $m = \pi n$, glmnet's problem gives the LASSO solution
$\hat{\beta}$, the gradient $G$ at $\hat{\beta}$, and the information $\hat{Q}$ at the
one-step relaxed fit on the selected model $E$ (see `?glmnet_problem`). `LassoInference`
linearizes the selection with $\hat{Q}$ and selects with

$$
Z + \omega = -G + \hat{Q}\hat{\beta}.
$$

The full-data statistic uses the same matrix and the full-data gradient $G_{full}$, at one
Newton step for the full-data loss on $E$ from $\hat{\beta}$,
$\tilde{\beta} = \hat{\beta} - H_{E,E}^{-1} G_{full,E}$:

$$
Z = -G_{full}(\tilde{\beta}) + \hat{Q}\tilde{\beta}.
$$

Both are $\hat{Q}$ times approximations to the unpenalized MLE, on the subsample and on all
the data. With $I$ the per-observation Fisher information, the subsample MLE minus the
full-data MLE is asymptotically independent of the full-data MLE, with variance
$I^{-1}(1/m - 1/n)$, and $\hat{Q} \approx I$. So

$$
\Sigma = \text{Var}(Z) \approx \hat{Q} / n, \qquad
\text{Var}(\omega) \approx \frac{1 - \pi}{\pi}\, \Sigma,
$$

and selecting on a subsample is randomized selection with `scalar_noise` $= (1 - \pi) / \pi$,
which is 1 for $\pi = 1/2$. `glmnet_inference` with `selection_rows` does all of this.

## Data

```{r}
library(lassoinf)
set.seed(0)

n <- 600; p <- 30
beta <- c(0.6, -0.6, 0.4, rep(0, p - 3))
make_data <- function(n) {
  X <- matrix(rnorm(n * p), n, p)
  y <- rbinom(n, 1, plogis(drop(X %*% beta)))
  list(X = X, y = y)
}
d <- make_data(n)
```

## Select on half of the data with glmnet

We fit a path that ends at a fixed $\lambda$, with the default standardization and an
unpenalized intercept. The coordinates are $(\text{intercept}, \beta)$, so index 0 is the
intercept. `glmnet_inference` warns if the fit does not satisfy the KKT conditions closely
enough, so a small `thresh` matters.

```{r}
pi_selection <- 0.5
lam <- 0.06
lambda_path <- lam * exp(seq(log(10), 0, length.out = 30))

select <- function(d) {
  rows <- sample(nrow(d$X), round(pi_selection * nrow(d$X)))
  fit <- glmnet::glmnet(d$X[rows, ], d$y[rows], family = "binomial",
                        lambda = lambda_path, control = list(thresh = 1e-14))
  list(fit = fit, rows = rows)
}
sel <- select(d)
which(coef(sel$fit, s = lam)[-1] != 0)
```

## Inference

```{r}
li <- glmnet_inference(sel$fit, d$X, d$y, s = lam, selection_rows = sel$rows, level = 0.9)
li$scalar_noise
summ <- li$summary()
summ$variable <- ifelse(summ$index == 0, "(Intercept)", paste0("x", summ$index))
summ
```

The estimate for the selected model is $W Z_E$, with $W = \hat{Q}_{E,E}^{-1}$, and its target is
$W (\hat{Q}\beta^*)_E$, with $\beta^* = (0, \beta)$ the true coefficients including the
intercept. When $E$ contains all the true signals, the target is $\beta^*_E$. Otherwise it is
the projection of $\beta^*$ onto the coordinates in $E$, in the metric given by the information
matrix $\hat{Q}$.

```{r}
target <- function(li) drop(li$W %*% (li$Q_hat %*% c(0, beta))[li$E + 1])
summ$truth <- target(li)
summ$cover <- summ$lower_conf <= summ$truth & summ$truth <= summ$upper_conf
summ[, c("variable", "lower_conf", "upper_conf", "p_value", "truth", "cover")]
```

For comparison, `li$"_naive"` holds Wald intervals that use all the data and ignore
selection. Data splitting refits the selected model on the half of the data not used for
selection:

```{r}
splitting <- function(d, rows, li, level = 0.9) {
  held_out <- setdiff(seq_len(nrow(d$X)), rows)
  X1 <- cbind(1, d$X)[held_out, li$E + 1, drop = FALSE]
  fit <- glm.fit(X1, d$y[held_out], family = binomial())
  est <- unname(coef(fit))
  mu <- fit$fitted.values
  sd <- sqrt(diag(solve(crossprod(X1, X1 * (mu * (1 - mu))))))
  q <- qnorm(1 - (1 - level) / 2)
  data.frame(index = li$E, estimate = est, lower_conf = est - q * sd, upper_conf = est + q * sd)
}
li$`_naive`
splitting(d, sel$rows, li)
```

The held-out refit targets the coefficients of the best-fitting logistic model on $X_E$. That
is the same as the target above when $E$ contains all the true signals.

## Coverage over repeated experiments

```{r}
one_replication <- function() {
  d <- make_data(n)
  sel <- select(d)
  li <- glmnet_inference(sel$fit, d$X, d$y, s = lam, selection_rows = sel$rows, level = 0.9)
  summ <- li$summary()
  keep <- summ$index != 0   # report the selected variables, not the intercept
  truth <- target(li)[keep]
  covers <- function(df) (df$lower_conf <= truth & truth <= df$upper_conf)
  lengths <- function(df) df$upper_conf - df$lower_conf
  tabs <- list(carving = summ[keep, ], naive = li$`_naive`[keep, ],
               splitting = splitting(d, sel$rows, li)[keep, ])
  data.frame(method = rep(names(tabs), each = sum(keep)),
             cover = unlist(lapply(tabs, covers)), length = unlist(lapply(tabs, lengths)),
             signal = rep(beta[summ$index[keep]] != 0, length(tabs)))
}

set.seed(1)
results <- do.call(rbind, replicate(100, one_replication(), simplify = FALSE))
data.frame(coverage = tapply(results$cover, results$method, mean),
           median_length = tapply(results$length, results$method, median))
```

Carving and data splitting both have close to the nominal 90% coverage. Carving also uses
the half of the data that selected the model, so its intervals are shorter. The naive
intervals ignore selection and undercover, almost entirely on the selected null variables:

```{r}
with(results, tapply(cover, list(method, signal = signal), mean))
```