---
title: "Detailed Function Presentation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Detailed Function Presentation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "ragg_png",
  dpi = 192,
  fig.width = 7,
  fig.height = 4.5,
  out.width = "90%",
  fig.align = "center",
  warning = FALSE,
  message = FALSE
)
```

# PDRobust

This vignette demonstrates validation, prediction, diagnostics, and treatment effect estimation using the bundled data.

``` text
data("BiSample") -> Mapping() -> DataCheck() -> DataStandard() -> prediction / diagnostic / analysis functions
```

## 1. Load the built-in package data

```{r}
library(PDRobust)
data("BiSample", package = "PDRobust")
head(BiSample)
```

## 2. Define roles and analysis settings with `Mapping()`

```{r}
mapping <- Mapping(
  id = "id",
  time = "time",
  treatment = "A",
  survival = "S",
  outcome = "Y",
  baseline_time = 0,
  cutoff_time = 2,
  covariates = c("X1", "X3", "X4", "X5", "X6"),
  interest_vars = c("X1", "X4"),
  y_type = "B"
)

print(mapping)
```

## 3. Validate the raw data with `DataCheck()`

```{r}
check <- DataCheck(BiSample, mapping, strict = FALSE)
names(check)           
```

```{r}
check$valid
check$ready_for_analysis
check$manual_resolution_required
check$can_standardize
```

The itemized report is in `check$checks`; supporting details are in `check$diagnostics`.

```{r}

```

## 4. Standardize the panel with `DataStandard()`

```{r}
pd_data <- DataStandard(BiSample, mapping, drop =TRUE)
head(pd_data)
```

### Imperfect dataset

For imperfect dataset, we have:

```{r}
data("ImperfectConSample", package = "PDRobust")
head(ImperfectConSample)
```

```{r}
con_mapping <- Mapping(
  id = "patient_id",
  time = "visit_month",
  treatment = "treatment",
  survival = "alive_status",
  outcome = "clinical_outcome",
  baseline_time = 0,
  cutoff_time = 12,
  covariates = c("X1", "X2", "X3", "X4", "X5", "X6"),
  interest_vars = c("X1", "X2"),
  y_type = "C"
)

con_check <- DataCheck(ImperfectConSample, con_mapping, strict = FALSE)


```

```{r}
con_check$valid
con_check$ready_for_analysis
con_check$manual_resolution_required
con_check$can_standardize
```

```{r}
con_data <- DataStandard(ImperfectConSample, con_mapping, drop = TRUE)
head(con_data)
```

```{r}
print(dim(ImperfectConSample))
print(dim(con_data))
```

```{r}
names(attributes(con_data))

```

```{r}
attr_standard <- attributes(con_data)
attr_standard$pd_standardization$time_map
head(attr_standard$pd_standardization$id_map)
```

## 5.1 Prediction functions and Diagnostics

```{r}
ps_fo <- A ~ X1 + X3 + X4 + X5 + X6
prin_fo <- S ~ (X1 + X3 + X4 + X5 + X6 ) * A
out_fo <- Y ~ (X1 + X3 + X4 + X5 + X6) * A + S
```

### Propensity score model

```{r}
ps <- PSPred(
  ps_fo = ps_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  mapping = mapping
)
           
head(ps)
```

```{r ps_dgn, fig.alt = "Absolute standardized mean differences before and after propensity-score weighting."}
ps_diagnostic <- PSDiag(data = pd_data,
                        ps_fo = ps_fo)

print(ps_diagnostic)
```

### Principal score model

```{r}
p0 <- PrinPred(
  prin_fo = prin_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  a = 0,
  mapping = mapping
)

head(p0)
```

```{r pps_dgn, fig.alt = "Standardized principal-score balance statistics for the selected covariates."}
principal_diagnostic <- PrinSDiag(
  data = pd_data, 
  ps_fo = ps_fo, 
  prin_fo = prin_fo)

print(principal_diagnostic)
```

### Outcome model

```{r}
mu1 <- OutPred(
  out_fo = out_fo,
  fit_dat = pd_data,
  pred_dat = pd_data,
  a = 1,
  mapping = mapping
)

head(mu1)
```

```{r sa, fig.alt = "Estimated effect-modification coefficients over time at different outcome-noise variance ratios."}
set.seed(12345)
sensitivity <- SA(
  data  = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  ratiovec = c(0.05,0.1, 0.2)
)
print(sensitivity)
```

### Principal-stratum profiling with `QR()`

```{r qr}
principal_profile <- QR(
  data = pd_data,
  prin_fo = prin_fo,
  quantile_level = c(0.25, 0.50, 0.75)
)

print(principal_profile)
principal_profile$data
```

### Treatment-group odds ratios

```{r or_ci, fig.alt = "Cutoff survival odds ratios and confidence intervals within treatment group zero."}
or_control <- ORCI(
  data = pd_data,
  formula = S ~ X1 + X3 + X4,
  a = 0,
  conf_level = 0.95
)

print(or_control)
```

## 5.2 Heterogeneous treatment effect

The five bootstrap replications below are only for a fast demonstration. Substantive standard errors and confidence intervals require more replications and an assessment of their stability. Use `B = 0` for point estimates alone.

```{r htesept, fig.alt = "Time-specific treatment-effect model coefficients and demonstration bootstrap confidence intervals."}
set.seed(12345)
separate_hte <- HTESepT(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  target_time = c(1, 2),
  B = 5,
  conf_level = 0.95,
  max_attempts = NULL,
  verbose = TRUE
)

separate_hte$summary
separate_hte$forest_plot
```

```{r}
head(separate_hte$boot_mat)
```

```{r hteallt, fig.alt = "Pooled treatment-effect model point estimates; bootstrap intervals are not calculated in this example."}
pooled_hte <- HTEAllT(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  B = 0,
  verbose = FALSE
)
pooled_hte$summary
pooled_hte$forest_plot
```