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.8470
## Compelete simulating 1th Dataset in 0.127 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      10.2       1
##  2     2     1       1 time_1       8.61      1
##  3     3     1       1 time_1      11.1       1
##  4     4     1       1 time_1      11.8       1
##  5     5     1       1 time_1       9.31      1
##  6     6     1       1 time_1       9.93      1
##  7     7     1       1 time_1      12.1       1
##  8     8     1       1 time_1      10.8       1
##  9     9     1       1 time_1       8.61      1
## 10    10     1       1 time_1       9.82      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,] 10.151302 11.262896 11.179125 12.336607 10.755550
##  [2,]  8.611930 11.315514 11.047152  8.561220  7.506208
##  [3,] 11.126120 11.577182 11.281983 12.611838 10.208543
##  [4,] 11.752735  9.988710 10.556013 10.078319 13.065829
##  [5,]  9.310427  7.443341 10.430644  9.985622  9.999425
##  [6,]  9.928467  7.545930  9.074725  9.999663 10.257611
##  [7,] 12.061901 10.388941 10.011796 10.701000 10.952909
##  [8,] 10.837329  7.520669  8.066801  9.329762  9.596770
##  [9,]  8.607893 10.033156 10.347485 10.462500  8.356089
## [10,]  9.817994 10.042937  8.784328  9.639264 10.249570
## [11,]  9.831326  7.420905 10.871118  8.206978 10.934924
## [12,] 10.904633  9.363302 10.021266  7.767827  7.793121
## [13,]  9.475139 10.077049 11.001252  8.484578  9.309441
## [14,] 12.363950 10.386280 11.160094 12.081534 10.244775
## [15,] 10.813434  9.566561  8.773451  7.975407 10.113657
## [16,]  9.974961  7.466789  8.384725  8.277364  9.975148
## [17,]  8.747835  7.514959  8.552011  7.887507  7.529547
## [18,]  8.859644  7.266393  8.915675  7.742814  7.527039
## [19,] 11.118052  7.519730 11.051283  8.193885  9.045561
## [20,]  8.512230  9.403146  9.703764  7.945986  7.822941
## [21,] 10.818910  9.376858 10.082516 11.843499 10.787450
## [22,]  7.643678  7.916587 10.130554 10.696183  9.084092
## [23,]  9.159416 10.128242 11.049562 10.787312  9.481023
## [24,] 10.807697 11.202243  9.856331  9.814826  9.190226
## [25,]  8.913029 10.024658 11.057402  9.948391  9.672249
## [26,] 10.435270  8.926810  8.063619  8.617271  9.650825
## [27,] 10.312298  9.567646  8.193035  9.308108  9.216468
## [28,] 10.634116  9.316347 10.243001 12.050085 10.461273
## [29,] 11.117572  9.189028 10.431482 11.631533 10.356882
## [30,]  7.854881  9.570265 10.887115 12.013592  9.008375
## [31,] 11.116546 11.076438 11.068651  9.777916  8.080101
## [32,] 10.293340 11.113292  9.020730  9.642551  9.338996
## [33,] 11.483485  7.651354  9.091693  9.571358 10.031431
## [34,]  9.706704 10.131359 10.335925 10.764851 10.965803
## [35,]  9.129314  9.537628 11.426575 10.870881 10.005800
## [36,] 12.330935  7.578337 10.026179  8.736298 10.020030
## [37,] 10.706147  7.515011 10.182795  8.639284  8.696482
## [38,]  9.058985  8.831833 10.134090  9.762597 10.906351
## [39,] 10.765721  9.489975 10.308114  8.594825 10.020310
## [40,]  9.320699  8.915188 11.009409  8.509022  9.679105
## [41,]  9.834014 10.425242 10.769724 10.065575  9.604299
## [42,]  8.714811  9.325008 10.332542  8.938965 10.238226
## [43,]  9.645307 11.255328  9.315267  8.630444  7.805637
## [44,]  9.847089  7.450173 10.678122  9.537576  9.211284
## [45,]  9.302916  7.188245 10.072721  8.559050  9.575030
## [46,]  7.885152  7.571727 10.232449  9.856905  9.451851
## [47,] 12.035023 10.357907 10.347732  8.733867  8.229622
## [48,]  9.605216 10.550296  8.620983  9.578498  9.035317
## [49,]  9.849406  9.572772  8.574798  8.482307  8.961818
## [50,] 10.949633  9.534782 11.050449  9.751648  9.270638

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.