---
title: "ICONIC: Causal Model Selection with Genetic Instruments and Negative Controls"
author: "Sean T. Bresnahan and Charis Xiong"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
    toc_float: true
vignette: >
  %\VignetteIndexEntry{ICONIC Walkthrough}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
subtitle: Model selection, sensitivity analysis, and prospective design for observational omics
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  fig.width = 7, fig.height = 4.5, comment = "#>"
)

set.seed(1)
library(iconic)
has_torch <- check_torch_setup()
```

# Introduction

This vignette walks through the complete `iconic` workflow on a simulated dataset: assembling an `iconic_data` object, diagnosing estimator eligibility, estimating natural direct and indirect effects, stress-testing the result, and planning a prospective study.

# Installation

`iconic` is being submitted to Bioconductor. Once accepted, install the
release version with:

```{r install-bioc, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("iconic")
```

The development version is available from GitHub:

```{r install-github, eval=FALSE}
BiocManager::install("sbresnahan/iconic")
```


# What this package is for

`iconic` is a workbench for causal inference in observational omics studies
where unmeasured confounding biases ordinary regression. The motivating
application is estimating the causal effect of an exposure (e.g. gestational
diabetes mellitus) on a panel of molecular outcomes (e.g. placental
transcriptome), using genetic instruments and negative controls to identify
the effect.

The package fits eight estimators of the natural direct and indirect effects
(NDE / NIE), diagnoses which are valid for *your* data, stress-tests them
against confounding and pleiotropy violations, and recommends the one most
likely to be unbiased. A generative texture model (hybrid GAN + copula) lets
the sensitivity analysis mirror the marginal and joint structure of your
cohort instead of relying on generic simulation defaults.

This vignette walks through the core workflow on simulated data. The
statistical theory -- the eight estimators, the negative-control assumptions
(A1, A2, A2'), the completeness condition, the sensitivity surface, and the
model-selection criterion -- is developed in full in the accompanying
manuscript and technical supplement.

# Quick start

The workflow is six steps:

1. **Standardize** your data with `iconic_data()`.
2. **Diagnose** which estimators are eligible with `iconic_diagnose()`.
3. **Estimate** the NDE/NIE with all eligible estimators via `iconic_estimate()`.
4. **Stress-test** robustness with `iconic_sensitivity()`.
5. **Recommend** a single estimator with `iconic_recommend()`.
6. **Prospect** for future studies with `iconic_prospect()`.

We illustrate on a simulated mediation panel with a genetic instrument (`G`),
an exposure (`X`), a mediator panel (`M`), an outcome panel (`Y`), and a panel
of negative controls (`W`).

```{r workflow-data}
set.seed(1)
dat <- generate_toy_data(n = 200, n_features = 10, seed = 1)

# Assemble an iconic_data object. Y and M are features x samples matrices;
# X is a length-n vector; G is an n x n_snps matrix; W is an n x n_ncs matrix.
idat <- iconic_data(
  X = dat$X, Y = dat$Y, M = dat$M,
  G = dat$G, W = dat$W,
  covariates = dat$synthetic_data
)
```

# Step 1: Diagnose

`iconic_diagnose()` reports instrument strength (partial F for `G` and `Gm`),
negative-control validity screens (A1, A2, A2'), completeness, and an
eligibility table for all eight estimators.

```{r workflow-diagnose}
diag <- iconic_diagnose(idat)
diag$eligibility
```

Each estimator is marked eligible only when its required assumptions are
met. For example, `IV2SLS` requires a valid instrument (F >= 10), while
`PGC` additionally requires negative-control completeness. The
instrument-variable estimators (`IV2SLS`, `IV2SLS2`) are identified by the
instrument(s) alone, so they remain eligible even when no negative controls
are supplied; `W`, when present, is used as an optional augmentation.

# Step 2: Estimate

`iconic_estimate()` fits every eligible estimator and returns per-feature
NDE, NIE, standard error, and p-value.

```{r workflow-estimate}
est <- iconic_estimate(idat)
head(est)
```

The `se_method` argument controls standard errors (`"delta"` by default;
`"bootstrap"` for confidence intervals) and the `p_value_method` argument
controls the null-distribution p-value (`"sobel"` or `"jt-comp"` for higher
power under sparse signals).

# Step 3: Train a texture model once

The sensitivity and prospective analyses use a generative texture model to
simulate data that mirrors your cohort's marginal and joint structure. Train
it **once** and attach it to the `iconic_data` object so every downstream
function reuses it without retraining.

For a runnable example we train a texture model briefly here; in practice you
would train for more epochs on your own cohort with `train_gan_on_real_data()`
(see the *Texture Model* vignette).

```{r workflow-gan, eval=has_torch}
# Train a texture model once (few epochs for the vignette).
input <- load_real_input_data(
  X_matrix = matrix(dat$X, nrow = 1),
  Y_matrix = t(dat$Y), M_matrix = t(dat$M),
  W_matrix = t(dat$W), covariates_df = dat$synthetic_data
)

gan <- train_gan_on_real_data(
  input$gan_training_data,
  feature_correlations = input$feature_correlations,
  feature_texture = input$feature_texture,
  epochs = 5, seed = 1, verbose = FALSE
)

# Attach it: iconic_sensitivity() and iconic_prospect() will now reuse it.
idat$trained_gan <- gan
```

```{r workflow-no-torch, eval=!has_torch, echo=FALSE, results="asis"}
cat(
"The sensitivity, recommendation, and prospective chunks in this vignette ",
"are skipped because torch is not available. Install torch with ",
"`install.packages(\"torch\")` and run `torch::install_torch()` to run them."
)
```

The texture model requires `torch`; install it with
`install.packages("torch")` before training.

# Step 4: Sensitivity analysis

`iconic_sensitivity()` sweeps the exogeneity of the instrument
(`rho_G1`, `rho_G2`) and reports how each estimator's bias degrades. It
reuses the attached GAN, so no retraining happens here.

```{r workflow-sensitivity, eval=has_torch}
sens <- iconic_sensitivity(idat, n_iter = 3)
sens
```

The output reports the maximum absolute bias across the sensitivity surface
for each estimator, the tipping points (first `rho_G2` at `rho_G1 = 0` where
`|bias| > 0.10`), and the texture source. Use `confounding = "inferred"` to
calibrate the held-fixed confounding parameters from the data:

```{r workflow-confounding}
conf <- infer_confounding(idat, diag, est)
conf
```

`infer_confounding()` estimates the confounding strength, mediator-outcome
confounding, negative-control coverage (`omega_1`, `omega_2`), and the number
of latent confounders (`k`) from the gap between estimators and the
negative-control residual variance. Pass `confounding = "inferred"` to
`iconic_sensitivity()` to use these values instead of defaults.

# Step 5: Recommend

`iconic_recommend()` combines eligibility, point estimates, and the
sensitivity surface into a single ranked recommendation.

```{r workflow-recommend, eval=has_torch}
rec <- iconic_recommend(idat, diag, est, sens)
rec
```

The recommended estimator is the eligible estimator with the best
per-estimand robustness across the sensitivity surface (NDE and NIE are
ranked separately). The full ranking and the ineligible estimators (with
reasons) are printed below the recommendation.

# Step 6: Prospective analysis

`iconic_prospect()` answers the design question: *if you don't yet have a
strong instrument, how strong does it need to be?* It sweeps the instrument
strength (`gamma_G`) and reports how estimates converge, then simulates a
prospective study at a target strength.

```{r workflow-prospect, eval=has_torch}
pros <- iconic_prospect(idat, n_iter = 3)
pros
```

The output reports the F-statistics at each strength level and recommends
which estimator would be best if you collected an instrument of that quality.
It also reuses the attached GAN.

# Survival outcomes

`iconic` supports time-to-event outcomes via Cox log-hazard-ratio and RMST
(restricted mean survival time) effect scales. RMST is collapsible and
admits an exact NDE/NIE decomposition; the log-HR scale is reported for
familiarity.

```{r workflow-survival}
set.seed(1)
sdat <- generate_toy_data(
  n = 100, n_features = 3,
  outcome_type = "survival", seed = 1
)

# For survival, pass surv_time + surv_event instead of Y.
sdata <- iconic_data(
  X = sdat$X,
  surv_time = sdat$surv_time, surv_event = sdat$surv_event,
  G = sdat$G, W = sdat$W,
  outcome_type = "survival",
  covariates = sdat$synthetic_data
)

sest <- iconic_estimate(sdata)
head(sest)
```

# Putting it together

The full pipeline, end to end:

```r
# 1. Standardize data (optionally attach a pre-trained GAN).
dat <- generate_toy_data(n = 200, n_features = 10, seed = 1)
idat <- iconic_data(X = dat$X, Y = dat$Y, M = dat$M,
                    G = dat$G, W = dat$W, covariates = dat$synthetic_data)

# 2. Train a texture model once and attach it.
input <- load_real_input_data(
  X_matrix = matrix(dat$X, nrow = 1), Y_matrix = t(dat$Y),
  M_matrix = t(dat$M), W_matrix = t(dat$W),
  covariates_df = dat$synthetic_data
)
idat$trained_gan <- train_gan_on_real_data(
  input$gan_training_data,
  feature_correlations = input$feature_correlations,
  feature_texture = input$feature_texture,
  epochs = 100, seed = 1
)

# 3. Diagnose, estimate, stress-test, recommend.
diag <- iconic_diagnose(idat)
est <- iconic_estimate(idat)
sens <- iconic_sensitivity(idat, confounding = "inferred")
rec <- iconic_recommend(idat, diag, est, sens)

# 4. Prospect for future studies (reuses the attached GAN).
pros <- iconic_prospect(idat, confounding = "inferred")
```

# Further reading

- **Manuscript and technical supplement** -- the eight estimators, the
  negative-control assumptions, the completeness condition, the sensitivity
  surface, and the model-selection criterion.
- **Function reference** -- `help(package = "iconic")` for the full API,
  including `gan_sensitivity()`, `pleiotropy_sweep()`,
  `mediation_sensitivity()`, and the negative-control diagnostics
  (`nc_validity_screen()`, `nc_independence_check()`,
  `nc_completeness_check()`).

# Session information

```{r sessionInfo}
sessionInfo()
```