Modgo simulation"
In modgo: Mock Data Generation

knitr::opts_chunk$set(fig.width = 8,
                      fig.height = 8,
                      out.height = "80%",
                      out.width = "80%",
                      dpi = 300)
counter <- 0

library(survival)

1. Default modgo simulation

For illustration, we selected the Cleveland Clinic Heart Disease Data set from the University of California in Irvine (UCI) machine learning data repository [@Dua2019]. Below, we are using eleven variables, five of which are continuous, four are dichotomous, and two categorical variables.

library(modgo)
data("Cleveland", package = "modgo")

# Specifying dichotomous and ordinal categorical variables
binary_variables <- c("Sex","HighFastBloodSugar","CAD","ExInducedAngina")
categorical_variables <- c("Chestpaintype","RestingECG")
nrep <- 500
plot_variables <- c("Age", "STDepression", binary_variables[c(1,3)], categorical_variables)

In this section, we run modgo with its default settings. For modgo to produce results that mimic the original data set efficiently, user needs to specify dichotomous and ordinal categorical variables. Variables will be considered as continuous, otherwise. All modgo runs in this and the following sections will produce 500 data sets with the specification nrep = 500; the default is 100.

Figure 1 shows the correlation plots for the default modgo run, and Figure 2 displays the distribution plots for the original data set and one simulated data set. The default displayed simulated data set is the first one. Moreover, for all the plots a set of variables are used.

test <- modgo(data = Cleveland,
              bin_variables = binary_variables,
              categ_variables = categorical_variables,
              nrep = nrep)

corr_plots(test, variables = plot_variables)
counter <- counter  + 1

distr_plots(test, variables = plot_variables)
counter <- counter  + 1

2. Expansions

2.1 Selection by thresholds of variables

modgo provides an option so that only subjects (instances) are simulated that fulfill a specific requirement. In the simplest case (Section 2.1), the user can specify an upper or a lower boundary, or an interval for a variable. The use may alternatively specify a combination of variables and thresholds.

Three steps are required when subjects need to fulfill a specific selection criterion for a continuous variable. First, the name of the variable needs to be specified, for which the threshold needs to be set. Second, the left and right boundaries need to be specified. Third, a data frame with three columns is defined with Column 1: variable name of threshold variable, Column 2: left boundary, i.e., lower bound, Column 3: right boundary, i.e., upper bound. Finally, the data frame is imported using the thresh_var argument. In the example, all subjects have to be at least 66 years old. The selection variable therefore is Age with left threshold 65 and right threshold infinity NA.

If the percentage of samples fulfilling the indicated threshold requirements are less than 10% of the simulated samples, modgo stops to avoid excessive computation time. However, users can force thresh_force = TRUE the requested simulation to be run.

Figure 3 shows the correlation plot for this illustration. Substantial differences between the original and the simulated correlation plots can be observed for the RestingECG and several other variables. Figure 4 displays the corresponding distribution plot. The age distribution is shifted as expected. Furthermore, the distribution of subjects with coronary artery disease (CAD = 1) is higher in the simulated than the original data set.

Variables <- c("Age")
thresh_left <- c(65)
thresh_right <- c(NA)
thresholds <- data.frame(Variables, thresh_left, thresh_right)

print(as.matrix(thresholds))

test_thresh <-  modgo(data = Cleveland,
                      bin_variables = binary_variables,
                      categ_variables = categorical_variables,
                      thresh_var = thresholds,
                      nrep = nrep,
                      thresh_force = TRUE)

corr_plots(test_thresh, variables = plot_variables)
counter <- counter  + 1

distr_plots(test_thresh, variables = plot_variables)
counter <- counter  + 1

2.2 Perturbation analysis - Unchanged variance

For continuous variables, modgo provides the option to add a normally distributed noise with mean 0 and variance $\sigma_{p}^2$. With this perturbation, the variance of the perturbed variable is identical to the variance of the original variable. This option permits the generation of values from continuous variables, which were not observed in the original data set.

To specify which variables are to be perturbed and to which degree, i.e., percentage, the user needs to provide modgo with a named vector of the percentages and with the corresponding variables names as the names of the vector.

Similar to the previous examples, Figure 5 shows the correlation plots for the expansion to perturbations, and Figure 6 displays the distribution plots. Figure 6 shows that the distribution of both resting blood pressure and cholesterol change substantially due to the perturbation.

#Create named vector
perturb_vector <- c(0.9,0.7)
names(perturb_vector) <- c("RestingBP","Cholsterol")

test_pertru <-  modgo(data = Cleveland,
                      bin_variables = binary_variables,
                      categ_variables = categorical_variables,
                      pertr_vec = perturb_vector,
                      nrep = nrep)

corr_plots(test_pertru, variables = c(plot_variables, names(perturb_vector)))
counter <- counter  + 1

distr_plots(test_pertru, variables = c(plot_variables, names(perturb_vector)))
counter <- counter  + 1

3. Generalized lambda modgo simulation

Another feature provided by modgo is the ability to simulate a data set using the Generalized Lambdas Distribution (GLD) method. This method allows users to generate values that are not present in the original data set, as the default method only uses values from the original data. The GLD method is based on the GLDEX package [@SteveSU]. More information on Generalized Lambdas Distributions can be found in Fitting Statistical Distributions [@Zaven_Karian].

test_GLD <- modgo(data = Cleveland,
                  bin_variables = binary_variables,
                  categ_variables = categorical_variables,
                  generalized_mode = TRUE,
                  nrep = nrep)

corr_plots(test_GLD, variables = plot_variables)
counter <- counter  + 1

distr_plots(test_GLD, variables = plot_variables)
counter <- counter  + 1

GLDEX provides three basic models for calculating the four Lambdas for each distribution. These models are called rmfmkl (default model), rprs, and star [@SteveSU]. They can also be combined for a bi-modal estimation. We give you the option to specify your desired model(or a combination of models) for each variable in the dataset. In the next step, we will show you how to specify the desired GLD models.

Variables <- c("Age","STDepression")
Model <- c("rprs", "star-rmfmkl")
model_matrix <- cbind(Variables,
                      Model)

test_GLD_define_model <- modgo(data = Cleveland,
                  bin_variables = binary_variables,
                  categ_variables = categorical_variables,
                  generalized_mode = TRUE,
                  generalized_mode_model = model_matrix,
                  nrep = nrep)

corr_plots(test_GLD_define_model, variables = plot_variables)
counter <- counter  + 1

distr_plots(test_GLD_define_model, variables = plot_variables)
counter <- counter  + 1

By examining the distribution plots shown above, we can observe that the GLD method has the capability to generate values for certain variables, such as STDepression, that do not exist in the original data set and can also be extremely high. Additionally, there is an alternative option to compute Generalized Lambdas independently of the modgo function and subsequently utilize them as an input for a subsequent modgo run.

gener_lambdas_matrix <- generalizedMatrix(data = Cleveland,
                                          generalized_mode_model = model_matrix,
                                          bin_variables = binary_variables) 
test_GLD_define_model_set_lambdas <- modgo(data = Cleveland,
                  bin_variables = binary_variables,
                  categ_variables = categorical_variables,
                  generalized_mode = TRUE,
                  generalized_mode_lmbds = gener_lambdas_matrix,
                  nrep = nrep)

Lastly, an intriguing feature provided by modgo in conjunction with the GLD method is the capability to simulate a data set without requiring the original data. To execute modgo without a data set, the user must provide the following:
1) Correlation matrix of the data set
2) Generalized matrix of the data set
3) Sample size of the simulated data set
4) Variable names
Below, we present an example of this case.

# Necessary arguments
gener_lambdas_matrix <- generalizedMatrix(data = Cleveland,
                                          generalized_mode_model = model_matrix,
                                          bin_variables = binary_variables)
sigma <- cor(Cleveland)
variables_names <- colnames(sigma)
sample_size <- 100

test_GLD_no_data_set <- modgo(data = NULL,
                              variables = variables_names,
                              bin_variables = binary_variables,
                              categ_variables = categorical_variables,
                              sigma = sigma,
                              generalized_mode = TRUE,
                              generalized_mode_lmbds = gener_lambdas_matrix,
                              n_samples = sample_size,
                              nrep = nrep)

corr_plots(test_GLD_no_data_set, variables = plot_variables)
counter <- counter  + 1

4. Survival example

To demonstrate the simulation of survival variables, we chose the cancer data set from the survival package [@survival-package]. This data set contains 167 samples and 10 variables. In order to set up modgo_survival(), the user must specify a status variable and a time variable, in addition to the other arguments of modgo.

# cancer prepare
data("cancer", package = "survival")

cancer <- na.omit(cancer)
cancer$sex <- cancer$sex - 1
cancer$status <- cancer$status - 1

time_var_cancer <- "time"
status_var_cancer <- "status"
bin_var_cancer <- c("status", "sex")
cat_var_list_cancer <- c("ph.ecog")

plot_variables_surv <- colnames(cancer)[1:6]

The modgo_survival function divides the data set into two separate data sets based on the status variable, and then it simulates each data set individually using the Generalized Lambdas Distribution method. The user can specify which GLDEX model should be used for each data set, with the default being "rprs".

# Survival run
test_surv <- modgo_survival(data = cancer,
               surv_method = 1,
               bin_variables = bin_var_cancer,
               categ_variables = cat_var_list_cancer,
               event_variable = status_var_cancer,
               time_variable = time_var_cancer,
               generalized_mode_model_no_event = "rmfmkl",
               generalized_mode_model_event = "rprs")

corr_plots(test_surv, variables = plot_variables_surv)
counter <- counter  + 1

distr_plots(test_surv, variables = plot_variables_surv)
counter <- counter  + 1

In the following plot, the surv_fit() curves from the survival package for both original and simulated data are depicted.

data_set_info <- c(rep("Original", dim(test_surv$original_data)[1]),
                   rep("Simulated", dim(test_surv$simulated_data[[1]])[1]))
combine_data_set <- rbind(test_surv$original_data,
                          test_surv$simulated_data[[1]])
combine_data_set <- cbind(combine_data_set,
                          data_set_info)
fit <- survfit(Surv(time, status) ~ data_set_info,
               data=combine_data_set)
plot(fit,
     fun = "F",
     col=1:2)

legend(700, 1,
       c("Original data set", "Simulated data set"),
       lty=c(1,1),
       col=c(1,2),
       bty='n',
       lwd=2)