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

## ----message = FALSE----------------------------------------------------------
library(oalasso)

set.seed(20260702)
n <- 300
X <- matrix(rnorm(n * 6), n, 6, dimnames = list(NULL, paste0("x", 1:6)))
A <- rbinom(n, 1, plogis(X %*% c(1, 1, 0, 0, 1, 1)))
Y <- as.numeric(2 * A + X %*% c(0.6, 0.6, 0.6, 0.6, 0, 0) + rnorm(n))
sim <- data.frame(A, Y, X)

fit <- oal(A ~ x1 + x2 + x3 + x4 + x5 + x6, data = sim, outcome = ~ Y)

## the path records, per grid point: lambda2, delta, lambda.n = n^delta,
## the paired gamma = 2*(2 - delta + 1), the exact glmnet s actually used,
## the wAMD, and the selected-set size
fit$path[, c("delta", "lambda.n", "gamma", "s", "wamd", "n.selected")]
all.equal(fit$path$gamma, 2 * (2 - fit$path$delta + 1))

## ----eval = FALSE-------------------------------------------------------------
# oal(A ~ ., data = dat, outcome = ~ y,
#     lambda = c(-10, -5, -2, -1, -0.75, -0.5, -0.25, 0.25, 0.49) + 0.5)

## ----eval = FALSE-------------------------------------------------------------
# ## Reference construction: what oal() computes internally
# f  <- A ~ x1 + x2 + x3 + x4 + x5 + x6
# X  <- model.matrix(f, sim)[, -1, drop = FALSE]
# Xs <- scale(X, TRUE, TRUE)                     # the package's fixed protocol
# b  <- coef(lm(sim$Y ~ sim$A + Xs))[-(1:2)]     # drop intercept and treatment
# names(b) <- colnames(Xs)
# fit.hook <- oal(f, data = sim, outcome.coef = b)   # same selections as outcome = ~ Y
#                                                    # (provenance label differs)
# 
# ## Ridge outcome model for p close to n
# set.seed(1)                                    # cv.glmnet folds are random
# cv <- glmnet::cv.glmnet(x = cbind(A = sim$A, Xs), y = sim$Y, alpha = 0,
#                         penalty.factor = c(0, rep(1, ncol(Xs))))
# b  <- as.numeric(coef(cv, s = "lambda.min"))[-(1:2)]
# names(b) <- colnames(Xs)
# fit.ridge <- oal(f, data = sim, outcome.coef = b)

## ----eval = FALSE-------------------------------------------------------------
# ## EXPERIMENTAL — screening use only, per the ten points above. Requires `party`.
# library(party)
# 
# f  <- A ~ x1 + x2 + x3 + x4 + x5 + x6
# X  <- model.matrix(f, sim)[, -1, drop = FALSE]
# Xs <- scale(X, TRUE, TRUE)                 # same standardization as oal()
# df <- data.frame(Y = sim$Y, Xs, check.names = FALSE)
# 
# set.seed(20260702)                         # forests are stochastic
# cf <- cforest(reformulate(colnames(Xs), response = "Y"), data = df,
#               controls = cforest_unbiased(ntree = 1000,
#                                           mtry  = ceiling(sqrt(ncol(Xs)))))
# vi <- varimp(cf, conditional = TRUE)       # conditional permutation importance
# vi <- pmax(vi[colnames(Xs)], 0)            # negative importance -> 0 = hard exclusion
# 
# fit.imp <- oal(f, data = sim, outcome.coef = vi)
# fit.imp    # provenance line flags the user-supplied weights; screening use only