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

## ----setup--------------------------------------------------------------------
library(subsampling)

## -----------------------------------------------------------------------------
set.seed(2)
N <- 5000
Z1 <- rbinom(N, 1, 0.04)
Z2 <- rbinom(N, 1, 0.07)
X1 <- rnorm(N)
X2 <- rnorm(N)
eta <- 0.5 + 0.5 * Z1 + 0.5 * Z2 + 0.5 * X1 + 0.5 * X2
Y <- rbinom(N, 1, plogis(eta))

data <- data.frame(Y, Z1, Z2, X1, X2)
formula <- Y ~ .
rareFeature.index <- 1:2
n.plt <- 300
n.ssp <- 700

## -----------------------------------------------------------------------------
fit_bl <- ssp.glm.rF(
  formula = formula,
  data = data,
  n.plt = n.plt,
  n.ssp = n.ssp,
  family = "quasibinomial",
  criterion = "BL-Uni",
  rareFeature.index = rareFeature.index
)

summary(fit_bl)

## -----------------------------------------------------------------------------
fit_rlopt <- ssp.glm.rF(
  formula = formula,
  data = data,
  n.plt = n.plt,
  n.ssp = n.ssp,
  family = "quasibinomial",
  criterion = "R-Lopt",
  balance.X.plt = TRUE,
  rareFeature.index = c("Z1", "Z2")
)

summary(fit_rlopt)

## -----------------------------------------------------------------------------
fit_auto <- ssp.glm.rF(
  formula = formula,
  data = data,
  n.plt = n.plt,
  n.ssp = n.ssp,
  family = "quasibinomial",
  criterion = "BL-Uni",
  rareFeature.index = NULL,
  rareThreshold = 0.09
)

fit_auto$rareFeature.index

## -----------------------------------------------------------------------------
fit_y_balanced <- ssp.glm.rF(
  formula = formula,
  data = data,
  n.plt = n.plt,
  n.ssp = n.ssp,
  family = "quasibinomial",
  criterion = "BL-Uni",
  balance.Y.ssp = TRUE,
  rareFeature.index = c("Z1", "Z2")
)

c(
  full_Y_rate = mean(data$Y),
  subsample_Y_rate = fit_y_balanced$Y.proportion.ssp
)

## -----------------------------------------------------------------------------
fit_estimated <- ssp.glm.rF(
  formula = formula,
  data = data,
  n.plt = n.plt,
  n.ssp = n.ssp,
  family = "quasibinomial",
  criterion = "R-Lopt",
  balance.X.plt = TRUE,
  rareFeature.index = c("Z1", "Z2"),
  control = list(poi.method = "estimated", b = 2),
  record.stage.time = TRUE
)

fit_estimated$stage.time

## -----------------------------------------------------------------------------
set.seed(3)
N_g <- 3000
Z <- rbinom(N_g, 1, 0.05)
X <- rnorm(N_g)
y <- 1 + 0.6 * Z + 0.2 * X + rnorm(N_g)
gaussian_data <- data.frame(y, Z, X)

fit_gaussian <- ssp.glm.rF(
  y ~ .,
  data = gaussian_data,
  n.plt = 200,
  n.ssp = 500,
  family = "gaussian",
  criterion = "BL-Uni",
  rareFeature.index = "Z"
)

summary(fit_gaussian)