---
title: "Integrative Bayesian Multi-Platform Regression with IntegMultiReg"
author: "Sinian Zhang, Thierry Chekouo"
output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Integrative Bayesian Multi-Platform Regression with IntegMultiReg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.2
)
set.seed(1)
```

# Overview

`IntegMultiReg` implements the *integrative multi-regression* (IMR) model of
Chekouo, Stingo, Doecke and Do (2017) and extends it from time-to-event
outcomes to continuous (Gaussian) and binary (probit) outcomes.  The central
idea is to integrate **all** available subjects across several molecular
platforms, even when many subjects are measured on only a subset of the
platforms.

Given $K$ platforms, subjects are partitioned into the $2^K-1$ non-empty
**availability subgroups** of a $K$-set Venn diagram.  One regression model is
built per subgroup, and information is shared across availability subgroups
through

* **non-local (product moment) priors** on the regression coefficients, which
  sharply penalise small effects and yield parsimonious models, and
* a **Markov random field (MRF) prior** on the variable-selection indicators,
  which encourages the same biomarkers to be selected across availability
  subgroups.

```{r load}
library(IntegMultiReg)
data("simIMR", package = "IntegMultiReg")
sapply(simIMR$platforms, dim)
```

The bundled `simIMR` data set has three platforms (`genomic`, `proteomic`,
`metabolomic`) measured on overlapping but partially missing sets of subjects,
clinical covariates for everyone, and three outcome types generated from the
same latent signal.  The truly associated features are recorded in
`simIMR$truth`.

# Real TCGA-KIRC example

The package also includes `kircIMR`, a reduced real-data example derived from
public UCSC Xena TCGA kidney renal clear cell carcinoma (KIRC) sampleMap files,
not from controlled-access TCGA/GDC files.  It mirrors the case-study structure
of Chekouo et al. (2017): mRNA expression, miRNA expression and DNA methylation
measured on overlapping patient sets, clinical covariates, and a right-censored
survival outcome.  Patient IDs are package-internal labels such as `KIRC001`;
the package does not distribute TCGA barcodes, and users should not attempt
participant re-identification or linkage to external resources.

```{r load-kirc}
data("kircIMR", package = "IntegMultiReg")
sapply(kircIMR$platforms, dim)
kircIMR$model_subgroup_sizes
```

The original Biometrics analysis used a much larger screened panel
(776 mRNA, 91 miRNA and 729 methylation features) and long MCMC chains.  The
package object keeps a smaller Cox-screened panel (50, 30 and 50 features) so
that examples remain lightweight while retaining the same scientific data
structure.  A short survival fit follows the same call pattern:

```{r kirc-fit, eval = FALSE}
kirc_fit <- imr(
  platform_data_list = kircIMR$platforms,
  outcome            = kircIMR$outcome.survival,
  cov                = kircIMR$covariates,
  type_outcome       = "right.censored",
  nu                 = c(-4, -3, -4),
  sample_mcmc        = c(4000, 1000),
  ssize              = 30,
  seed               = 1
)
summary(kirc_fit)
plot_top_features(kirc_fit, top = 12)
```

# Fitting the model

The workhorse is `imr()`.  Here we fit a binary (probit) outcome with a short
MCMC run for illustration; in practice use longer chains.

```{r fit}
fit <- imr(
  platform_data_list = simIMR$platforms,
  outcome            = simIMR$outcome.binary,
  cov                = simIMR$covariates,
  type_outcome       = "binary",
  nu                 = c(-4, -3, -4),
  sample_mcmc        = c(1500, 500),
  ssize              = 30,
  seed               = 1
)
fit
```

The four modelled availability subgroups correspond to the availability patterns
`011` (genomic + proteomic), `111` (all three), `101` (genomic + metabolomic)
and `100` (metabolomic only); the bitstring digits run from the first platform
(right) to the last (left).  `plot_subgroup_sizes()` shows how the sample splits
across them.

```{r plot-sizes, fig.height = 3.6}
plot_subgroup_sizes(fit)
```

# Inspecting selected biomarkers

`summary()` ranks, per platform, the features whose marginal posterior
inclusion probability (mPIP) exceeds a threshold.

```{r summary}
summary(fit, threshold = 0.5)
```

The posterior inclusion probabilities themselves are available through `coef()`
(a list of subgroup $\times$ feature matrices) and can be visualised either as
per-platform heatmaps or as a single ranked bar chart of the most strongly
supported features across all platforms.

```{r plot-selection}
plot(fit, type = "selection", platform = 1)
```

```{r plot-top, fig.height = 4}
plot_top_features(fit, top = 8)
```

The estimated MRF interaction parameters, which measure how strongly each pair
of subgroups shares biomarkers, are shown with `type = "theta"`, and the
log-posterior trace with `type = "trace"`.

```{r plot-trace}
plot(fit, type = "trace")
```

# Predicting new subjects

`predict()` routes new subjects to the appropriate availability subgroup,
standardises their features with the training scaling, and performs Bayesian
model averaging.  When the new platforms are supplied in the same order as at
training, `platform_names` can be omitted.  Binary outcomes are returned on the
probability scale.

```{r predict}
new_x <- simIMR$platforms$genomic[1:20, ]
new_p <- simIMR$platforms$proteomic[1:20, ]
pred <- predict(fit, newdata = list(new_x, new_p),
                covariates = simIMR$covariates)
head(pred[["model:011"]])
```

# Assessing predictive performance

`cv_imr()` runs repeated $K$-fold splits using the fitted MCMC samples and
reports the accuracy measure appropriate to the outcome type (AUC for binary,
the concordance index for survival, mean squared error for continuous), overall
and within each availability subgroup.  The MCMC sampler is not re-run inside
each fold; to compare IMR with BMS, fit a separate `method = "BMS"` object.

```{r cv}
cv <- cv_imr(fit, k = 5, rounds = 3)
attr(cv, "metric")
round(colMeans(cv$total_cindex), 3)
```

# Other outcome types

The same interface handles continuous and right-censored outcomes; only
`type_outcome` and the `outcome` data frame change.

```{r other-types, eval = FALSE}
# continuous (Gaussian) outcome
fit_c <- imr(
  simIMR$platforms, simIMR$outcome.continuous, cov = simIMR$covariates,
  type_outcome = "continuous", nu = c(-4, -3, -4),
  sample_mcmc = c(1500, 500), ssize = 30, seed = 1)

# right-censored survival outcome
fit_s <- imr(
  simIMR$platforms, simIMR$outcome.survival, cov = simIMR$covariates,
  type_outcome = "right.censored", nu = c(-4, -3, -4),
  sample_mcmc = c(1500, 500), ssize = 30, seed = 1)
```

# References

Chekouo T, Stingo FC, Doecke JD, Do K-A (2017). "A Bayesian Integrative
Approach for Multi-Platform Genomic Data: A Kidney Cancer Case Study."
*Biometrics*, **73**(2), 615--624. <doi:10.1111/biom.12587>