---
title: "Diagnostics, profiling, and sensitivity analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Diagnostics, profiling, and sensitivity analysis}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 80
---

```{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
)
```

## Introduction

The complete workflow is illustrated as follows. This article focuses on the
package's diagnostics, profiling and sensitivity analysis, including five
functions `PSDiag()`, `PrinSDiag()`, `SA()`, `OR()` and `ORCI()`.

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

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

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

pd_data <- DataStandard(BiSample, map)
```

```{r}
head(pd_data)
```

```{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
```

## Propensity score, covariate balance

To assess the adequacy of the propensity score model specification, `PSDiag()`
evaluate covariate balance before and after weighting. Before weighting,
standardized mean differences(SMD) are calculated using the original pooled
standard-deviation denominator. After ordinary inverse probability of treatment
weighting, the denominator is calculated using the corresponding weighted
effective sample sizes.

The argument `data` specifies the standardized dataset used for the diagnostic
analysis, whereas `ps_fo` specifies the propensity score model formula. The
returned object contains the following components:

```{r}
ps_diag <- PSDiag(data = pd_data, 
                  ps_fo = ps_fo)
names(ps_diag)
```

The two primary outputs are data frames containing the standardized mean
differences before and after weighting and the corresponding diagnostic plot.

A suitably specified propensity score model should improve covariate balance
after weighting. Accordingly, the absolute SMD **should generally decrease
toward zero, with values below 0.1** commonly regarded as indicating acceptable
residual imbalance. However, satisfactory covariate balance **does not** by
itself establish that the propensity score model is correctly specified or
eliminate the possibility of unmeasured confounding.

```{r, fig.alt = "Absolute standardized mean differences before and after propensity-score weighting."}
print(ps_diag)
```

Additional components provide further information about the diagnostic
procedure. For example, `ps_diag$weights` contains the calculated inverse
probability weights, and `ps_diag$weight_type` identifies the weighting method
used.

```{r}
ps_diag$weight_type
```

## Principal score, covariate-specific balance statistic

`PrinSDiag()` computes a standardized statistic obtained by comparing the
weighted covariates contribution of surviving subjects across both treatment
groups, using both propensity score model and principal score model.

The argument `data` specifies the standardized dataset used for the diagnostic
analysis, `ps_fo` specifies the propensity score model formula and `prin_fo`
specifies the principal score model. The returned object contains the following
components:

```{r}
prin_diag <- PrinSDiag(data = pd_data, 
                       ps_fo = ps_fo, 
                       prin_fo = prin_fo)
names(prin_diag)
```

The two primary outputs are data frames containing the standardized statistics
and the corresponding diagnostic plot. Values close to zero indicate the
residual discrepancies; values between -1.96 to 1.96 are also acceptable. Values
outside that range warrant further examination of the principal score model and
propensity score model.

```{r, fig.alt = "Standardized principal-score balance statistics for the selected covariates."}
print(prin_diag)
```

The returned object also provides additional diagnostic components, such as the
cumulative principal scores under treatment levels 0 and 1.

```{r}
head(prin_diag$p0)
head(prin_diag$p1)
```

## Outcome-noise sensitivity analysis

`SA()` evaluates the sensitivity of the estimated heterogeneous treatment
effects to additional unexplained variation in the outcome. This function
estimates how the estimated effect-modification coefficients change when random
outcome noise is introduced.

The argument `data` specifies the standardized dataset used in the sensitivity
analysis. The arguments `ps_fo`, `prin_fo`, and `out_fo` specify the propensity
score, principal score, and conditional outcome model formulas, respectively.
The argument `ratiovec` specifies the noise level and it should be a list,
defaulted by `c(0, 0.05, 0.10)` The returned object contains the following
components:

```{r}
set.seed(20160878)
sa <- SA(
  data = pd_data,
  ps_fo = ps_fo,
  prin_fo = prin_fo,
  out_fo = out_fo,
  ratiovec = c(0, 0.05, 0.1)
)
names(sa)
```

The two primary outputs are data frame containing the estimated coefficients
across different perturbation level and the corresponding curves for the
covariates of interest.

Greater similarity in the magnitude and direction of estimated coeffcients
across increasing noise levels indicates greater robustness of estimated
heterogeneous treatment effects.

This comparison concerns the specified random outcome-noise perturbations. It
does not test the causal identifying assumptions or implement the
principal-ignorability sensitivity parameter in the methodological paper.

```{r, fig.alt = "Estimated effect-modification coefficients over time at different outcome-noise variance ratios."}
print(sa)
```

The returned object also includes supplementary diagnostic information. For
example, `variance_by_time` records the empirical outcome variance used to scale
the perturbation at each analysis time, `convergence` summarizes the
estimating-equation solution for each scenario, `model_diagnostics` contains
model-fitting diagnostics, and `warnings` records consolidated warnings
generated during the analysis.

```{r}
sa$variance_by_time
sa$warnings
```

## Principal-stratum summaries

`QR()` characterizes the distribution of selected covariates within the
estimated "always-survivor" principal stratum at the cutoff time. This function
calculate the cumulative principal score under treatment level 0 and uses these
values as subject-specific weights. When level 0 represents the unexposed or
reference condition, subjects with a higher estimated probability of surviving
under that condition receive greater weight in the principal-stratum summaries.

The argument `data` specifies the standard dataset used for principal stratum
summaries, `prin_fo` specifies the principal model formula, and `quantile_level`
specifies one or more quantile levels, and defaluts to 0.5. It can also be a
numeric vector such as `c(0.5, 0.95)`.

```{r}
profile <- QR(
  data = pd_data,
  prin_fo = prin_fo,
  quantile_level = c(0.5, 0.95)
)

names(profile)
```

The primary outputs are principal-score-weighted means and quantiles. For each
selected variables, `QR()` calculates weighted means at the cutoff time. For
non-binary variables, it additionally calculates the weighted quantile using an
intercept-only weighted quantile regression model.

```{r}
print(profile)
```

Additional components include `profile$binary` , which identifies variables
treated as binary, and `profile$weights`, which contains the full-precision
cumulative principal-score weights used in the calculations.

```{r}
head(profile$weights)
```

## Treatment-specific survial odds ratios

`ORCI()` estimates the association between selected covariates and survival odds
within a specified treatment group. This function fits a logistic regression
model and exponentiates the non-intercept coefficients to obtain odd ratios and
Wald confidence interval.

The argument `data` specifies the standardized dataset used for this function.
The argument `formula` specifies the logistic regression model, with the mapped
survival variable as the response variable on the left-hand side and covariates
on the right-hand side. The argument `a` specifies the treatment group used for
this analysis. The argument `conf_level` specifies confidence level and defaults
to 0.95.

The returned object contains the following components:

```{r}
or_fo <- S ~ X1 + X2 + X4
or0 <- ORCI(data = pd_data, 
            formula = or_fo,
            a = 0,
            conf_level = 0.95)

names(or0)
```

The two primary outputs are estimated odds ratios and their Wald confidence
interval, together with the corresponding plot. An odds ratio equal to 1
indicates no estimated association with survival on the odds scale.

For a continuous variable, an odd ratio greater than 1 indicates that higher
covariate values are associated with higher survival odds at cutoff time point,
conditional on the other covariates in the model. For a categorical variable, an
add ratio compares the specified category with its reference category.

The confidence intervals are based on the fitted logistic regression
coefficients and their estimated standard errors. A confidence interval that
excludes 1 provides evidence of an association at the corresponding confidence
level, whereas an interval containing 1 indicates that the direction of
association remains uncertain. Very wide intervals may indicate limited
information, sparse outcome events, poor covariate overlap, or unstable model
estimation within the selected group.

```{r, fig.alt = "Cutoff survival odds ratios and confidence intervals within treatment group zero."}
print(or0)
```

The returned object also contains supplementary information, such as
or0\$model_diagnostics, which stores the full-precision fitted logistic
regression model.

```{r}
or0$model_diagnostics
```