Worked Example Using the IST Dataset

Introduction

This R Markdown document illustrates example usage of the RESIDE package using the IST dataset.

Setup

Load the RESIDE package and set a seed for reproducibility and store the folder directory for export / import.

# Load the Library
library(RESIDE)
# Set the seed
set.seed(1234)
# Store the folder path used for import / export
folder_path <- tempdir()

Summarise original data

Select the variables of interest and summarise. Selection of variables is optional, if the variables are known. They may not be known until the marginal distributions have been received.

# Select variables of interest from the IST dataset.
IST_original <- IST |> dplyr::select(
  AGE, # AGE at Randomisation
  SEX, # SEX M/F
  RATRIAL, # Atrial Fibrillation Y/N at Randomisation 
  # (not coded for 984 patients in the pilot phase)
  RSBP, # Systolic Blood Pressure at Randomisation
  STRK14 # Indicator of Any Stroke at 14 days
)

# Convert the character variables to factors (to allow for summary)
IST_original <- IST_original |> dplyr::mutate_if(is.character, factor)

# Produce a summary of the variables
summary(IST_original)
#>       AGE        SEX       RATRIAL        RSBP           STRK14       
#>  Min.   :16.00   F: 9028    :  984   Min.   : 70.0   Min.   :0.00000  
#>  1st Qu.:65.00   M:10407   N:15282   1st Qu.:140.0   1st Qu.:0.00000  
#>  Median :73.00             Y: 3169   Median :160.0   Median :0.00000  
#>  Mean   :71.72                       Mean   :160.2   Mean   :0.04152  
#>  3rd Qu.:80.00                       3rd Qu.:180.0   3rd Qu.:0.00000  
#>  Max.   :99.00                       Max.   :295.0   Max.   :1.00000

Fit a cox model on the original data

# Load survival and dplyr libraries
library(survival) # For Cox PH model
library(dplyr) # For data manipulation

# Stroke event is measured at 14 days, so set this for patients
IST_original$DAY <- 14

# Illustrate the 984 missing values
sum(IST_original$RATRIAL == "")
#> [1] 984

# Remove the missing values
IST_original <- IST_original[!IST_original$RATRIAL == "",]

# Drop the factor name for the missing values
IST_original$RATRIAL <- droplevels(IST_original$RATRIAL)

# Summarise the variable to show there are no longer missing values
summary(IST_original$RATRIAL)
#>     N     Y 
#> 15282  3169

# Fit a Cox PH model
cox.ph <- coxph(Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, data = IST_original) 

# Output the summary of the Cox PH Model
cox.ph
#> Call:
#> coxph(formula = Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, 
#>     data = IST_original)
#> 
#>              coef exp(coef) se(coef)     z      p
#> AGE      0.005237  1.005251 0.003387 1.546 0.1220
#> SEXM     0.077489  1.080570 0.074388 1.042 0.2976
#> RATRIALY 0.231692  1.260732 0.091818 2.523 0.0116
#> RSBP     0.000994  1.000995 0.001310 0.759 0.4479
#> 
#> Likelihood ratio test=11.8  on 4 df, p=0.01893
#> n= 18451, number of events= 764

Get the marginal distributions

Use the get_marginal_distributions() function to get the marginal distributions, additionally selecting which variables using the variables parameter.

# Get the Marginal Distributions for the selected variables
marginals <- get_marginal_distributions(
  IST,
  variables = c(
    "AGE",
    "SEX",
    "RATRIAL",
    "RSBP",
    "STRK14"
  )
)

Export the marginal distributions

Export the marginal distributions using the export_marginal_distributions() function, using the force parameter to override any existing files.

# Export the Marginal Distributions
export_marginal_distributions(marginals,
                              folder_path = folder_path,
                              force = TRUE)
#> Exporting  Categorical variables to:  C:\Users\Ryan\AppData\Local\Temp\RtmpMRdU12/categorical_variables.csv
#> Exporting  Binary variables to:  C:\Users\Ryan\AppData\Local\Temp\RtmpMRdU12/binary_variables.csv
#> Exporting  Continuous variables to:  C:\Users\Ryan\AppData\Local\Temp\RtmpMRdU12/continuous_variables.csv
#> Exporting  Summary to:  C:\Users\Ryan\AppData\Local\Temp\RtmpMRdU12/summary.csv

Reimport marginal distributions

Import the exported marginal distributions using the import_marginal_distributions() function

# Import the Marginal Distributions
imported_marginals <- import_marginal_distributions(folder_path = folder_path)

Synthesise data from marginal distributions (without correlations)

Synthesise data from the imported marginals using the synthesise_data function.

# Synthesise a dataset from the imported Marginal Distributions (without correlations)
sim_df <- synthesise_data(imported_marginals)

Summarise the synthesised data

Summarise the simulated data

# Convert any Character variables to Factors
sim_df <- sim_df |> dplyr::mutate_if(is.character, factor)
# Summarise the synthesised data
summary(sim_df)
#>        id        SEX       RATRIAL       STRK14             AGE       
#>  Min.   :    1   F: 9086    :  955   Min.   :0.00000   Min.   :16.00  
#>  1st Qu.: 4860   M:10349   N:15327   1st Qu.:0.00000   1st Qu.:65.00  
#>  Median : 9718             Y: 3153   Median :0.00000   Median :73.00  
#>  Mean   : 9718                       Mean   :0.04281   Mean   :71.55  
#>  3rd Qu.:14576                       3rd Qu.:0.00000   3rd Qu.:80.00  
#>  Max.   :19435                       Max.   :1.00000   Max.   :99.00  
#>       RSBP      
#>  Min.   : 72.0  
#>  1st Qu.:141.0  
#>  Median :160.0  
#>  Mean   :159.9  
#>  3rd Qu.:176.0  
#>  Max.   :295.0

Fit a cox model on the simulated data

Fit the same cox model as earlier except this time on the simulated data.


# As before the events are measured at day 14
sim_df$DAY <- 14

# Show that the missing observations are in the data
sum(sim_df$RATRIAL == "")
#> [1] 955

# Remove the missing observations
sim_df <- sim_df[!sim_df$RATRIAL == "",]

# Remove the missing factor name
sim_df$RATRIAL <- droplevels(sim_df$RATRIAL)

# Show that there are no missing observations
summary(sim_df$RATRIAL)
#>     N     Y 
#> 15327  3153

# Fit the model on the synthesised data
cox.ph.sim <- coxph(Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, data = sim_df) 

# Show a summary of the model
cox.ph.sim
#> Call:
#> coxph(formula = Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, 
#>     data = sim_df)
#> 
#>                coef  exp(coef)   se(coef)      z     p
#> AGE       0.0012645  1.0012653  0.0030791  0.411 0.681
#> SEXM      0.0232422  1.0235145  0.0719414  0.323 0.747
#> RATRIALY -0.0924565  0.9116889  0.0982948 -0.941 0.347
#> RSBP     -0.0005778  0.9994223  0.0012911 -0.448 0.654
#> 
#> Likelihood ratio test=1.4  on 4 df, p=0.8445
#> n= 18480, number of events= 778

Synthesise data with correlations

Synthesise data from the imported marginals with correlations, using the correlations parameter to specify a list of assumed correlations created with the correlation() function. Categorical variables are correlated using a single category, specified with factor_name.x or factor_name.y.

# Synthesise data specifying assumed correlations
sim_df_cor <- synthesise_data(
  imported_marginals,
  correlations = list(
    # Patients without atrial fibrillation are less likely to have a stroke
    correlation("RATRIAL", "STRK14", -0.2, factor_name.x = "N"),
    # Older patients have a higher systolic blood pressure
    correlation("AGE", "RSBP", 0.3)
  )
)

Summarise the synthesised data

Summarise the synthesised data (with correlations)

# Convert to Factors from Character variables
sim_df_cor <- sim_df_cor |> dplyr::mutate_if(is.character, factor)
# Summarise the synthesised dataset
summary(sim_df_cor)
#>        id        SEX       RATRIAL       STRK14             AGE      
#>  Min.   :    1   F: 9103    :  923   Min.   :0.00000   Min.   :18.0  
#>  1st Qu.: 4860   M:10332   N:15256   1st Qu.:0.00000   1st Qu.:65.0  
#>  Median : 9718             Y: 3256   Median :0.00000   Median :73.0  
#>  Mean   : 9718                       Mean   :0.03885   Mean   :71.5  
#>  3rd Qu.:14576                       3rd Qu.:0.00000   3rd Qu.:80.0  
#>  Max.   :19435                       Max.   :1.00000   Max.   :99.0  
#>       RSBP      
#>  Min.   : 70.0  
#>  1st Qu.:141.0  
#>  Median :160.0  
#>  Mean   :159.9  
#>  3rd Qu.:176.0  
#>  Max.   :278.0

Fit a cox model on the simulated data with correlations

Using the synthesised data (with correlations) fit a Cox PH model, with the same parameters as earlier.


# Again events are measured at 14 days
sim_df_cor$DAY <- 14

# Again check that the missing values where added
sum(sim_df_cor$RATRIAL == "")
#> [1] 923

# Again remove the missing values
sim_df_cor <- sim_df_cor[!sim_df_cor$RATRIAL == "",]

# Again drop the missing factor
sim_df_cor$RATRIAL <- droplevels(sim_df_cor$RATRIAL)

# Show there are no missing values
summary(sim_df_cor$RATRIAL)
#>     N     Y 
#> 15256  3256

# Fit the model on the synthesised data (with correlations)
cox.ph.sim.cor <- coxph(Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, data = sim_df_cor)

# Show a summary of the model
cox.ph.sim.cor
#> Call:
#> coxph(formula = Surv(DAY, STRK14) ~ AGE + SEX + RATRIAL + RSBP, 
#>     data = sim_df_cor)
#> 
#>                coef  exp(coef)   se(coef)      z        p
#> AGE      -0.0042511  0.9957579  0.0033247 -1.279    0.201
#> SEXM      0.0324843  1.0330176  0.0765993  0.424    0.672
#> RATRIALY  0.6712446  1.9566711  0.0840131  7.990 1.35e-15
#> RSBP     -0.0009317  0.9990687  0.0014465 -0.644    0.520
#> 
#> Likelihood ratio test=60.51  on 4 df, p=2.266e-12
#> n= 18512, number of events= 686