Introduction to R package copulaSim

What is the usage of the R package copulaSim?

This package serves to simulate jointly distributed patient-level data from historical data based on the copula invariance property.

Why do we need to perform virtual patient simulation?

To consistently optimize clinical trial designs and data analysis methods through trial simulation, we need to simulate multivariate mixed-type virtual patient data independent of designs and analysis methods under evaluation. To make the outcome of optimization more realistic, we should utilize relevant empirical patient-level data when it is available.

The challenges faced when simulating small data

When simulating small empirical data, the underlying marginal distributions and their dependence structure cannot be understood or verified thoroughly due to the limited sample size.

Solution: copula invariance property

To resolve this issue, we use the copula invariance property to generate the joint distribution without making a strong parametric assumption. The theoretical background is addressed below.

The demonstration of R Package copulaSim

Generate empirical data:

library(mvtnorm)
arm1 <- rmvnorm(n = 30, mean  = rep(10, 5), sigma = diag(5) + 0.5)
test_data <- as.data.frame(cbind(1:30, rep(1, 30), arm1))
colnames(test_data) <- c("id","arm",paste0("time_", 1:5))
knitr::kable((test_data), "simple")

Load package copulaSim

library(copulaSim)

Use function copula.sim to generate ONE simulated dataset for the emperical data

• As a means of avoiding extreme sampling results when performing one simulated dataset, it is advised to use “energy” or “ball” in the argument validation.type. The purpose of doing this is to perform data validation, which ensures the similarity between empirical data and the simulated data based on 2-sample test.

• Below shows the 2-sample test result, which indicates that the joint distribution of the simulated data is not significantly different from the empirical data.

## Generate 1 simulated dataset
simu_S1 <- copula.sim(data.input = test_data[,-c(1,2)], 
                      id.vec = test_data$id, 
                      arm.vec = test_data$arm,
                      n.patient = 50 , 
                      n.simulation = 1, 
                      seed = 2022, 
                      validation.type = "energy",
                      verbose = TRUE)

## Simulate 1th Dataset

## p.value for energy test: 0.7410
## Compelete simulating 1th Dataset in 0.112 seconds

## Obtain the simulated long-form dataset
simu_S1$data.simul

## # A tibble: 250 × 6
##       id   arm col.num col.name data.sim sim.id
##    <dbl> <dbl>   <dbl> <chr>       <dbl>  <int>
##  1     1     1       1 time_1      11.1       1
##  2     2     1       1 time_1       9.20      1
##  3     3     1       1 time_1      11.9       1
##  4     4     1       1 time_1      11.8       1
##  5     5     1       1 time_1       9.25      1
##  6     6     1       1 time_1       9.76      1
##  7     7     1       1 time_1      11.9       1
##  8     8     1       1 time_1      10.3       1
##  9     9     1       1 time_1       8.95      1
## 10    10     1       1 time_1      10.0       1
## # ℹ 240 more rows

• Also, from the marginal density plot, we can see that the copula simulated data is close to the empirical data.

library(dplyr)
## Obtain the empirical long-form dataset
empir <- simu_S1$data.input.long %>% mutate(cate = "empirical_n30") %>% rename(data = data.input)

## Produce the marginal density plot
simul <- simu_S1$data.simul %>%  mutate(cate = "copulaSim_n50") %>% 
         rename(data = data.sim) %>% select(-sim.id)  
library(ggplot2)
rbind(empir, simul) %>% filter(grepl('time', col.name)) %>%
  ggplot(aes(x = data, color = cate, fill = cate)) +
  facet_wrap(.~col.name, ncol = 5) +
  geom_density(alpha = 0.001, size = 1)

Use function extract.data.sim to convert simulated data into a list of wide-form matrices

## Converting the long-form simulated dataset to wide-form 
simu.wide <- extract.data.sim(simu_S1)
simu.wide

## $`sim.id=1`
## $`sim.id=1`$`arm=1`
##          time_1    time_2    time_3    time_4    time_5
##  [1,] 11.127030 11.325799 11.786916 11.769028 11.623775
##  [2,]  9.203379 11.054563 10.525130  8.630098  7.353626
##  [3,] 11.924163 11.812299 11.844487 11.806753 11.479405
##  [4,] 11.782184 10.671019 10.598892 10.565174 11.773225
##  [5,]  9.249597  7.703238 10.030171  9.672915  9.954657
##  [6,]  9.759826  8.532776  8.634592  9.728963 10.550438
##  [7,] 11.911844 10.918633  9.594576 10.703994 11.603200
##  [8,] 10.339377  8.527943  7.389786  9.167303  9.444960
##  [9,]  8.946892  9.539930  9.783363 10.031107  8.893265
## [10,] 10.045914  9.962378  8.532143  9.540387 10.551221
## [11,]  9.447876  7.729853 10.610530  8.381324 11.550052
## [12,] 10.865246  9.468578  9.129310  7.169852  7.720598
## [13,]  9.782136  9.642103 10.582774  8.811649  9.475999
## [14,] 11.996090 10.984104 11.788005 11.203946 11.218175
## [15,] 10.856761  9.613014  8.545528  8.221973  9.967647
## [16,]  9.631585  8.174877  7.988731  8.384340  9.780465
## [17,]  8.644228  7.449185  7.805120  7.101433  7.161039
## [18,]  8.530870  6.813254  8.090133  6.821707  7.152442
## [19,] 10.847863  8.523417 10.776287  8.367029  9.145951
## [20,]  8.588792  8.877201  8.746274  7.471506  7.839983
## [21,] 10.892390  9.528601  9.626109 10.812180 11.521295
## [22,]  7.780555  8.272926  9.424159 10.009823  9.137528
## [23,]  9.620048  9.649411 10.659039 10.540234  9.789219
## [24,] 11.280750 11.113207  9.125195  9.698404  9.349368
## [25,]  9.280459  9.578634 10.774475  9.675663  9.957034
## [26,] 10.228388  9.165494  7.380182  9.037418  9.663526
## [27,] 10.447642  9.574762  7.648439  9.224875  9.117408
## [28,] 10.851224  9.479127  9.923657 10.944821 11.484996
## [29,] 11.201100  9.474327 10.215350 10.769165 11.431480
## [30,]  8.636720  9.412743 10.196494 10.731829  9.289992
## [31,] 11.477822 10.999094 10.766194  9.674396  8.355948
## [32,] 10.889466 10.946969  8.568303  9.553431  9.448822
## [33,] 11.014796  8.605659  8.644384  9.451497  9.915779
## [34,] 10.131326 10.399864  9.998449 10.658874 11.623009
## [35,]  9.475226  9.483848 11.888333 10.664428 10.451769
## [36,] 11.615231  8.664343  9.563163  9.172283  9.924852
## [37,] 10.244199  8.414233  9.694495  8.964739  8.594206
## [38,]  9.064099  8.563456  9.671788  9.601551 11.528291
## [39,] 10.853139  9.539907  9.949911  9.098377  9.970234
## [40,]  9.398315  8.594853 10.614556  8.729805  9.867850
## [41,] 10.307585 10.550908 10.244333 10.015723  9.823887
## [42,]  8.875791  8.889099  9.953758  9.121014 10.525687
## [43,] 10.249062 11.108764  8.608275  9.028071  7.824461
## [44,]  9.580994  7.864856 10.176663  9.357681  9.236152
## [45,]  8.990349  6.761218  9.516359  8.519611  9.488778
## [46,]  8.020852  7.913411  9.681183  9.523703  9.634955
## [47,] 11.895511 10.747081  9.939754  9.188391  8.705108
## [48,] 10.120452 10.693661  8.089758  9.423969  9.116244
## [49,]  9.938313  9.540008  8.057750  8.765834  8.798658
## [50,] 11.132536  9.567156 10.739242  9.667163  9.475772

Use function compare.copula.sim to perform the comparison between empirical data and multiple simulated datasets.

• To increase the data variability, we suggest that the data validation option should be turned off when producing multiple simulated datasets. That is, the argument for validation.type should be set to “none”.

## Generate 100 simulated datasets
simu_S100 <- copula.sim(data.input = test_data[,-c(1,2)], 
                         id.vec = test_data$id, 
                         arm.vec = test_data$arm,
                         n.patient = 50 , 
                         n.simulation = 100, 
                         seed = 2022, 
                         validation.type = "none",
                         verbose = FALSE)

## Compare the marginal mean via the function compare.copula.sim 
compare <- compare.copula.sim(simu_S100)
knitr::kable(compare$mean.comparison, "simple")

Employ function new.arm.copula.sim to simulate new multivariate datasets with shifted mean vector from existing empirical data.

• This function assumes that the dependence structure of the new-arm data and that of the empirical data are the same.

## Generate Empirical Data 
##  Assume that the single-arm, 3-dimensional empirical data follows multivariate normal data
arm1 <- rmvnorm(n = 80, mean  = c(10,10.5,11), sigma = diag(3) + 0.5)
test_data2 <- as.data.frame(cbind(1:80, rep(1,80), arm1))
colnames(test_data2) <- c("id", "arm", paste0("time_", 1:3))

## Generate 1 simulated datasets with one empirical arm and two new-arms
## The mean difference between empirical arm and 
## (i) the 1st new arm is assumed to be 2.5, 2.55, and 2.6 at each time point
## (ii) the 2nd new arm is assumed to be 4.5, 4.55, and 4.6 at each time point

newARM <- new.arm.copula.sim(data.input = test_data2[,-c(1,2)], 
                             id.vec = test_data2$id,
                             arm.vec = test_data2$arm, 
                             n.patient = 100 , 
                             n.simulation = 1, 
                             seed = 2022,
                             shift.vec.list = list(c(2.5,2.55,2.6), c(4.5,4.55,4.6)),
                             verbose = FALSE)

## Obtain the simulated long-form dataset
newARM$data.simul

## # A tibble: 900 × 6
##       id   arm col.num col.name data.sim sim.id
##    <dbl> <dbl>   <dbl> <chr>       <dbl>  <int>
##  1     1     1       1 time_1       9.63      1
##  2     2     1       1 time_1       9.65      1
##  3     3     1       1 time_1       8.24      1
##  4     4     1       1 time_1      12.1       1
##  5     5     1       1 time_1       9.12      1
##  6     6     1       1 time_1      10.0       1
##  7     7     1       1 time_1       8.56      1
##  8     8     1       1 time_1      10.6       1
##  9     9     1       1 time_1      10.5       1
## 10    10     1       1 time_1       8.48      1
## # ℹ 890 more rows

## Verify the mean difference
newARM$data.simul %>%  
  group_by(.data$arm, .data$col.num) %>%
  summarise(N = n(), Mean = mean(.data$data.sim), SD = sd(.data$data.sim))

## # A tibble: 9 × 5
## # Groups:   arm [3]
##     arm col.num     N  Mean    SD
##   <dbl>   <dbl> <int> <dbl> <dbl>
## 1     1       1   100  9.71 1.01 
## 2     1       2   100 10.2  1.26 
## 3     1       3   100 10.8  1.23 
## 4     2       1   100 12.0  0.975
## 5     2       2   100 12.7  1.30 
## 6     2       3   100 13.3  1.14 
## 7     3       1   100 14.0  0.975
## 8     3       2   100 14.7  1.30 
## 9     3       3   100 15.3  1.14

Acknowledgement

This research project and the development of the R package are supported by AbbVie Experiential Internship Program. I am also grateful to Dr. Xuemin Gu, Dr. Jenny Jiao, and Dr. Jane Zhang at the Eyecare Clinical Statistics Team for valuable comments on this work.