Package {PDRobust}

PDRobust: Principal-stratification treatment-effect estimation

Description

PDRobust estimates treatment effects for longitudinal outcomes that may be unavailable after death. The effects apply to the always-survivor principal stratum: subjects who would survive through the selected cutoff time under either treatment. Mapping() identifies the variables and analysis times, and DataStandard() prepares the data for the remaining functions.

Details

PSPred(), PrinPred(), and OutPred() fit their models each time they are called and return numeric predictions. HTESepT() estimates effects separately at user-selected times, whereas HTEAllT() estimates a joint trajectory over every observed time from baseline through cutoff.

Treatment coding

The implemented estimator uses treatment 1 as the survival-favorable arm: potential survival satisfies S^1 \ge S^0 at cutoff. Its always-survivor principal score is therefore the survival probability under treatment 0. If the survival-favorable arm is coded as 0 in the raw data, recode the raw treatment as 1 - A before mapping and standardizing. To report the original contrast, negate the package estimate and transform an interval ⁠[lower, upper]⁠ to ⁠[-upper, -lower]⁠. Mapping() does not infer or reverse treatment coding.

Interpretation and assumptions

The target population comprises subjects who would survive through the selected cutoff under either treatment. The effect model describes outcome differences at earlier analysis times within that fixed population. Continuous-outcome models describe effects on a linear scale. Binary-outcome models transform the model's linear predictor with 2 * plogis(eta) - 1; their coefficients are not odds ratios.

Causal interpretation requires consistency, no interference, adequate treatment and survival overlap, treatment ignorability conditional on the measured covariates, the stated survival monotonicity, and principal ignorability. Checking the data and covariate balance cannot establish these assumptions. Triple robustness means that, under the method's assumptions and regularity conditions, at least two of the propensity score, principal score, and outcome mean models must be correctly specified. It does not guarantee unbiased results in every finite sample. Limiting extreme probabilities can also affect the estimates.

SA() examines sensitivity by adding random noise to the outcome. It does not vary the principal-ignorability assumption. See the method-and-coding vignette for treatment coding and other implementation limits.

Author(s)

Maintainer: Huan Wang whhuan42@gmail.com

Authors:

• Huan Wang whhuan42@gmail.com

• Yilin Zhang

References

Zhang, Y., Shardell, M., Falvey, J., McCoy, R., Stuart, E., and Chen, C. (2026). A Novel Tool for Evaluating Effect Modification in Older Adults with ADRD Using Medicare Claims. arXiv:2608.06654. doi:10.48550/arXiv.2608.06654.

See Also

Useful links:

• https://github.com/whhuan/PD_Robust

• https://whhuan.github.io/PD_Robust/

• Report bugs at https://github.com/whhuan/PD_Robust/issues

Binary longitudinal example data

Description

Simulated long-format data for illustrating binary-outcome analyses. Potential survival is generated with S^1 \ge S^0, matching the package's treatment-1 survival-favorable convention. The simulation variables are not observed counterfactual information available in a real study.

Usage

BiSample

Format

A simulated long-format data frame with 1,200 rows (400 subjects at three visits) and 16 variables:

id

Subject identifier.

time

Analysis time.

Pi

Simulated probability of treatment 1 conditional on baseline covariates, stored to three decimal places.

S1, S0

Simulated potential survival indicators under treatment 1 and 0, respectively; 1 denotes alive and 0 denotes dead.

Y1, Y0

Simulated binary potential outcomes under treatment 1 and 0, respectively. These simulation variables are retained for illustration; package analyses use the observed outcome Y.

X1, X2, X3

Continuous baseline covariates.

X4, X5, X6

Binary baseline covariates.

A

Binary treatment indicator.

S

Binary survival or intermediate-status indicator.

Y

Binary outcome, structurally missing after death.

Source

Simulated for package examples.

Examples

data("BiSample", package = "PDRobust")
head(BiSample)

Check whether longitudinal data are ready for analysis

Description

Checks the columns, values, visit structure, and analysis settings specified by mapping. Every observed time from baseline through the cutoff is treated as an analysis time, and the input data are left unchanged.

Usage

DataCheck(data, mapping, strict = FALSE)

Arguments

Value

A pd_data_check list with the following components:

valid

TRUE when no check classified as an error fails. Some warnings about encoding or ordering may still prevent analysis.

ready_for_analysis

TRUE when the data pass every check required for analysis.

manual_resolution_required

TRUE when a failed check requires the user to correct the data before standardization.

can_standardize

TRUE when no problem requires manual correction. Standardization can still fail if rows must be removed but drop = FALSE, or if removal leaves no observations or only one treatment group.

checks

A data frame with one row per performed check, including the result, its importance, details, and a recommended action.

settings

A list containing the validated mapping.

diagnostics

Detailed row indices, subject identifiers, and summary tables for the performed checks. Missing columns or empty input cause an early return with only the checks possible at that stage.

Numeric summaries intended for display are rounded to three decimals; counts, row indices, identifiers, and logical flags retain their types.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
check <- DataCheck(BiSample, map)
check$ready_for_analysis

Prepare longitudinal data for PDRobust analyses

Description

Converts supported binary values to 0 and 1, replaces subject IDs and analysis times with consecutive integers, sorts the data by subject and time, and stores the information needed by other PDRobust functions.

Usage

DataStandard(data, mapping, drop = FALSE)

Arguments

Value

A data frame inheriting from pd_data, retaining the input column names and additional unmapped columns. Rows outside the mapped time window are removed; retained rows are sorted by recoded ID and time. Attributes are:

pd_mapping

The mapping updated to use the standardized baseline and cutoff times. Column names remain unchanged.

pd_original_mapping

The mapping supplied for the raw data.

pd_check

The final pd_data_check report, including attrition. A warning is issued if the returned data are not ready for analysis, for example if deletion removes one treatment group.

pd_standardization

A list containing time_map, id_map, attrition, and initial_check. The maps link raw values to their standardized values; the initial check describes the input data.

The function checks the data both before and after preparation. The stored reports describe this call and are not recalculated if the returned data are later edited or subsetted; see pd_methods. The ID map has one row per retained subject, so its size increases with the number of subjects. Analysis values retain full precision; only summaries shown to users and percentages describing removed data are rounded.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
attr(pd_dat, "pd_mapping")

Joint estimation of heterogeneous treatment effects across time

Description

Performs joint longitudinal analysis of heterogeneous treatment effects across all time points with time included as a covariate.

Usage

HTEAllT(
  data,
  ps_fo,
  prin_fo,
  out_fo,
  B,
  conf_level = 0.95,
  max_attempts = NULL,
  verbose = TRUE,
  progress_callback = NULL
)

Arguments

Details

HTEAllT() combines predictions from the propensity score, principal score, and outcome mean models to form the data used for effect estimation. It then estimates one treatment-effect trajectory across all observed times and, when requested, uses bootstrap resampling to calculate confidence intervals.

Value

A pd_hte_pooled object containing the jointly estimated trajectory, the analysis times, confidence intervals when B > 0, model-checking information, and a summary of successful and failed bootstrap attempts. time_effect_estimable indicates whether the data allowed a time effect to be included. Displayed estimates are rounded to three decimal places; boot_mat stores the unrounded bootstrap coefficients.

Numerical safeguards

Propensity scores are limited to ⁠[0.01, 0.99]⁠. Their product with estimated survival probabilities under treatment 1 is limited to ⁠[0.005, 0.995]⁠ when effects are estimated. These limits prevent division by probabilities very close to zero, but they can affect the estimates and do not demonstrate adequate overlap or validate the causal assumptions. Review the returned model information, especially when few subjects remain at risk.

Treatment coding

The implemented estimator uses treatment 1 as the survival-favorable arm: potential survival satisfies S^1 \ge S^0 at cutoff. Its always-survivor principal score is therefore the survival probability under treatment 0. If the survival-favorable arm is coded as 0 in the raw data, recode the raw treatment as 1 - A before mapping and standardizing. To report the original contrast, negate the package estimate and transform an interval ⁠[lower, upper]⁠ to ⁠[-upper, -lower]⁠. Mapping() does not infer or reverse treatment coding.

See Also

PDRobust-package, pd_methods

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
fit <- HTEAllT(
  pd_dat,
  A ~ X1 + X2 + X4,
  S ~ X1 + X2 + X4 + A + time,
  Y ~ X1 + X2 + A,
  B = 0
)
fit$summary

Separate estimation of heterogeneous treatment effects by time point

Description

Performs separate analyses of heterogeneous treatment effects at each selected time point.

Usage

HTESepT(
  data,
  ps_fo,
  prin_fo,
  out_fo,
  target_time,
  B,
  conf_level = 0.95,
  max_attempts = NULL,
  verbose = TRUE,
  progress_callback = NULL
)

Arguments

Details

HTESepT() combines predictions from the propensity score, principal score, and outcome mean models to form the data used for effect estimation. It fits a separate treatment-effect model at each selected time and, when requested, uses bootstrap resampling to calculate confidence intervals.

Value

A pd_hte_timevarying object containing the estimate for each requested time, confidence intervals when B > 0, model-checking information, and a summary of successful and failed bootstrap attempts. Displayed estimates are rounded to three decimal places; boot_mat stores the unrounded bootstrap coefficients.

Numerical safeguards

Propensity scores are limited to ⁠[0.01, 0.99]⁠. Their product with estimated survival probabilities under treatment 1 is limited to ⁠[0.005, 0.995]⁠ when effects are estimated. These limits prevent division by probabilities very close to zero, but they can affect the estimates and do not demonstrate adequate overlap or validate the causal assumptions. Review the returned model information, especially when few subjects remain at risk.

Treatment coding

The implemented estimator uses treatment 1 as the survival-favorable arm: potential survival satisfies S^1 \ge S^0 at cutoff. Its always-survivor principal score is therefore the survival probability under treatment 0. If the survival-favorable arm is coded as 0 in the raw data, recode the raw treatment as 1 - A before mapping and standardizing. To report the original contrast, negate the package estimate and transform an interval ⁠[lower, upper]⁠ to ⁠[-upper, -lower]⁠. Mapping() does not infer or reverse treatment coding.

See Also

PDRobust-package, pd_methods

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
fit <- HTESepT(
  pd_dat,
  A ~ X1 + X2 + X4,
  S ~ X1 + X2 + X4 + A + time,
  Y ~ X1 + X2 + A,
  target_time = c(0, 2), B = 0
)
fit$summary

Imperfect Continuous Longitudinal Example Data

Description

A deliberately imperfect continuous-outcome longitudinal data set derived from a simulated continuous-outcome panel. The data mimic common issues encountered in raw clinical data exports while remaining recoverable using DataCheck and DataStandard with drop = TRUE.

Usage

ImperfectConSample

Format

A data frame with 599 rows and 11 variables in long format, with one row per recorded subject and visit:

patient_id

Noncanonical character subject identifier.

visit_month

Character-encoded visit time in months.

treatment

Character-encoded binary treatment assignment.

alive_status

Character-encoded binary survival or intermediate status.

X1, X2, X3

Continuous baseline covariates.

X4, X5, X6

Binary baseline covariates.

clinical_outcome

Continuous longitudinal clinical outcome.

Details

The data include nonstandard subject identifiers, character-encoded visit times and binary variables, unsorted records, an incomplete longitudinal record, missing required covariate values, a missing outcome among survivors, and a record with a missing subject identifier. Structural outcome missingness for records with alive_status = 0 is retained.

Source

Simulated for package examples.

Examples

data("ImperfectConSample", package = "PDRobust")
head(ImperfectConSample)

Identify variables and analysis times for PDRobust

Description

Records which columns contain the subject ID, time, treatment, survival, outcome, and covariates, together with the analysis time range and outcome type. Other package functions use this information to interpret the data consistently.

Usage

Mapping(
  id,
  time,
  treatment,
  survival,
  outcome,
  baseline_time,
  cutoff_time,
  covariates,
  interest_vars,
  y_type
)

Arguments

Details

All ten arguments are required. target_time is specified separately when calling HTESepT() because it selects time points for that analysis rather than describing the data.

Value

A pd_mapping object that can be supplied to DataCheck() and DataStandard().

Examples

map <- Mapping(
  id = "id", time = "time", treatment = "A",
  survival = "S", outcome = "Y",
  baseline_time = 3,
  cutoff_time = 9,
  covariates = c("X1", "X2", "X4"),
  interest_vars = c("X1", "X2"),
  y_type = "C"
)
map

Estimate covariate associations with survival at the cutoff time

Description

Estimates odds ratios with confidence intervals for associations between covariates and survival at the cutoff time within a selected treatment group.

Usage

ORCI(data, formula, a, conf_level = 0.95)

Arguments

Details

ORCI() fits the supplied logistic regression model using observations from treatment group a at the cutoff time. It reports an odds ratio and Wald confidence interval for every non-intercept coefficient that can be estimated. Covariates are not selected according to statistical significance.

Value

An odds_ratios object containing odds-ratio estimates and confidence intervals, the fitted model, model-checking information, and a forest plot. Reported estimates are rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
result <- ORCI(
  pd_dat, S ~ X1 + X2 + X4, a = 0
)
result$forestplotdat

Estimate outcome predictions

Description

Fits the outcome mean model and predicts each row of pred_dat under treatment a and survival status 1. The model is fitted again each time the function is called.

Usage

OutPred(out_fo, fit_dat, pred_dat, a, mapping, ...)

Arguments

Value

A numeric vector of predicted outcome means, one for each row of pred_dat, rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
mu1 <- OutPred(Y ~ X1 + X2 + A + S, pd_dat, pd_dat, a = 1, mapping = map)
head(mu1)

Evaluate how well a propensity score model performs

Description

Calculates the standardized mean difference (SMD) for each covariate before and after propensity score weighting.

Usage

PSDiag(data, ps_fo)

Arguments

Details

PSDiag() fits the propensity score model using baseline observations and uses inverse-probability-of-treatment weighting to make the treatment groups more comparable. It limits estimated probabilities to ⁠[0.01, 0.99]⁠ to avoid extremely large weights. A smaller absolute SMD after weighting indicates better balance for that covariate.

Value

A PSDiag object containing SMDs before and after weighting, the estimated propensity scores and weights, and a balance plot. SMDs are rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
result <- PSDiag(pd_dat, A ~ X1 + X2 + X4)
result$smd_after

Estimate propensity scores

Description

Fits a logistic model using baseline observations from fit_dat and returns each row's estimated probability of receiving treatment 1 in pred_dat. The model is fitted again each time the function is called.

Usage

PSPred(ps_fo, fit_dat, pred_dat, mapping, ...)

Arguments

Value

A numeric vector of propensity scores, one for each row of pred_dat, rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
ps <- PSPred(A ~ X1 + X2 + X4, pd_dat, pd_dat, map)
head(ps)

Estimate cumulative principal scores

Description

Fits the principal score model and returns each row's estimated probability of surviving from baseline through its observed time under treatment a. All observed times from baseline through cutoff are used, and the model is fitted again each time the function is called.

Usage

PrinPred(prin_fo, fit_dat, pred_dat, a, mapping, ...)

Arguments

Details

When the data contain multiple times, each post-baseline observation is used to model the next survival step only if the subject was alive at the previous observed time. If the data contain only one observed time, all complete observations at that time are used.

Value

A numeric vector of cumulative survival probabilities, one for each row of pred_dat, rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
score0 <- PrinPred(
  S ~ X1 + X2 + X4 + A + time,
  pd_dat, pd_dat, a = 0, mapping = map
)
head(score0)

Evaluate covariate balance for the principal score model

Description

Calculates a standardized balance statistic for each numeric covariate at the cutoff time after accounting for treatment assignment and estimated survival. Values nearer zero indicate better balance between the weighted treatment groups.

Usage

PrinSDiag(data, ps_fo, prin_fo)

Arguments

Details

The function fits both the propensity score and principal score models. It uses all observed times from baseline through cutoff to estimate cumulative survival probabilities, limits propensity scores to ⁠[0.01, 0.99]⁠, and then calculates the balance statistics at the cutoff time.

Value

A PrinSDiag object containing the standardized balance statistics, estimated probabilities, and a diagnostic plot. Balance statistics are rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
result <- PrinSDiag(
  pd_dat,
  A ~ X1 + X2 + X4,
  S ~ X1 + X2 + X4 + A + time
)
result$statistics

Summary statistics of covariates within the always-survivor principal stratum

Description

Estimates the user-specified quantile for continuous covariates and the mean of covariates for subjects within the always-survivor principal stratum.

Usage

QR(data, prin_fo, quantile_level = 0.5)

Arguments

Details

QR() uses estimated survival probabilities under treatment 0 to weight the numeric variables listed in interest_vars. It reports a weighted mean for every variable. For variables with more than two observed values, it also reports the requested weighted quantiles. Variables with no more than two observed values are treated as binary and receive a mean but no quantile.

Value

A QR object containing the weighted means, requested quantiles, variable-type indicators, and weights. Reported means and quantiles are rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
result <- QR(
  pd_dat,
  S ~ X1 + X2 + X4 + A + time,
  quantile_level = c(0.25, 0.5, 0.75)
)
result$mean

Sensitivity analysis of outcome mean model misspecification

Description

Examines how time-specific heterogeneous treatment effect estimates change when random noise is added to the outcome. The amount of noise is determined from the observed outcome variance at each analysis time, and subjects from both treatment groups at the cutoff contribute to effect estimation.

Usage

SA(data, ps_fo, prin_fo, out_fo, ratiovec = c(0, 0.05, 0.1))

Arguments

Details

For each observed analysis time and each value in ratiovec, SA() adds mean-zero random noise whose variance equals that value multiplied by the observed outcome variance. It then re-estimates the heterogeneous treatment effect for that time.

For continuous outcomes, the perturbed outcomes are used both to refit the outcome mean model and to estimate the treatment effect. For binary outcomes, the outcome mean model is fitted to the original binary outcomes, while the perturbed outcomes are used only during effect estimation. Binary treatment effects use the same bounded scale as HTESepT().

The three nuisance models are fitted again as needed for each scenario. The results show sensitivity to this particular form of random outcome noise; they do not cover every possible model error or violation of the causal assumptions. Set an R random seed before calling SA() to reproduce the same noise.

Value

An SA object containing estimates for every analysis time and noise level, model-checking information, warnings, and plots. Displayed estimates are rounded to three decimal places.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
pd_dat <- DataStandard(BiSample, map)
set.seed(20260912)
result <- SA(
  pd_dat,
  A ~ X1 + X2 + X4,
  S ~ X1 + X2 + X4 + A + time,
  Y ~ X1 + X2 + A,
  ratiovec = c(0, 0.05)
)
head(result$data)

Print, plot, and subset PDRobust results

Description

Provides standard ways to print analysis summaries, draw stored plots, and select rows or columns from data prepared by DataStandard(). Plotting a result does not refit its model.

Usage

## S3 method for class 'pd_mapping'
print(x, ...)

## S3 method for class 'pd_data_check'
print(x, ...)

## S3 method for class 'pd_data'
x[...]

## S3 method for class 'pd_hte_timevarying'
print(x, ...)

## S3 method for class 'pd_hte_pooled'
print(x, ...)

## S3 method for class 'PSDiag'
print(x, ...)

## S3 method for class 'PrinSDiag'
print(x, ...)

## S3 method for class 'odds_ratios'
print(x, ...)

## S3 method for class 'QR'
print(x, ...)

## S3 method for class 'SA'
print(x, ...)

## S3 method for class 'pd_hte_timevarying'
plot(x, ...)

## S3 method for class 'pd_hte_pooled'
plot(x, ...)

## S3 method for class 'PSDiag'
plot(x, ...)

## S3 method for class 'PrinSDiag'
plot(x, ...)

## S3 method for class 'odds_ratios'
plot(x, ...)

Arguments

Details

Subsetting preserves the stored information but does not check the data again. Before analyzing subsetted or edited data, validate them because removing rows or columns can break the required longitudinal structure. QR() has a print method but no package-specific plot method.

Value

Print methods show the main result and invisibly return x. When a result contains a plot, printing also draws it; SA draws all stored sensitivity plots. Plot methods invisibly return the stored ggplot object. Subsetting returns the selected data and preserves its PDRobust mapping and preparation information when the result remains a data frame.

Examples

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", "X4"),
  interest_vars = c("X1", "X2"), y_type = "B"
)
print(map)
print(DataCheck(BiSample, map))
prepared <- DataStandard(BiSample, map)
prepared[1:3, ]
diagnostic <- PSDiag(prepared, A ~ X1 + X2 + X4)
print(diagnostic)
p <- plot(diagnostic)