Simulation of data

[szymekdr]

The function is now ready, it's in R/simulate_data - just source file simBACE_new.R.

What you can setup:

simBACE <- function(
    response_type = "gaussian",  # type of the response
    predictor_types = c("gaussian", "gaussian"),  # type of each predictor
    var_names = NULL,  # custom names, can be left NULL - they will be x1, x2, ...
    beta_matrix = NULL,  # matrix defining dependencies between predictors (if NULL - generated automatically)
    beta_resp = NULL,  # betas for y~x1+x2+... dependencies (can be left NULL)
    beta_sparsity = 0.7,  # how sparse betas should be for predictors (0 = all predictors depend on each other, 1 = they are independent)
    ix_matrix = NULL,  # interactions between predictors specified using "Tukey-like" notation, see function description
    beta_ix = NULL,  # betas for interaction terms
    intercepts = NULL,  # intercepts for variables if not zero
    birth = 0.8,  # birth rate for tree generation
    death = 0.4,  # death rate for tree generation
    phylo_signal = NULL,  # vector of phylogenetic signals for variables y, x1, x2, ...
    n_cases = 200,  # total number of obs to generate
    n_species = 75,  # number of species on the tree
    missingness = NULL,  # vector of missigness scores for each variable, 0 = no missing data, y is the first variable
    sigmas = NULL,  # SD for random effects of non-phylogenetic species effect and residuals, if NULL - generated automatically
    rr = FALSE,  # if TRUE - allow for random slopes in gaussian and/or poisson predictors
    rr_form = NULL) # names of predictors where random slope is possible, separately for each random effect

[szymekdr]

What you get out of it: list of objects with data and parameters:

List of 4
 $ data          :'data.frame': 150 obs. of  6 variables:
  ..$ species: chr [1:150] "eHZDN" "REEEZ" "CxOyx" "goBMd" ...
  ..$ y      : num [1:150] 0.1175 4.8218 4.2323 -0.0728 -0.5029 ...
  ..$ x1     : num [1:150] 1.84 NA 2.45 NA 2.32 ...
  ..$ x2     : int [1:150] 0 NA NA 0 0 NA 0 0 0 NA ...
  ..$ x3     : chr [1:150] NA "B" "B" NA ...
  ..$ x4     : num [1:150] -3.316 -4.102 -5.196 -0.725 -4.446 ...
 $ tree          :List of 4
  ..$ edge       : int [1:94, 1:2] 49 50 55 61 62 63 77 77 63 70 ...
  ..$ edge.length: num [1:94] 0.00076 1.3353 1.57791 0.10763 0.16093 ...
  ..$ tip.label  : chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  ..$ Nnode      : int 47
  ..- attr(*, "class")= chr "phylo"
  ..- attr(*, "order")= chr "cladewise"
 $ params        :List of 17
  ..$ response_type   : chr "gaussian"
  ..$ predictor_types : chr [1:4] "gaussian" "binary" "multinomial3" "gaussian"
  ..$ var_names       : chr [1:5] "y" "x1" "x2" "x3" ...
  ..$ beta_matrix     : num [1:4, 1:4] 0 0.0626 0 0.4136 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ beta_resp       : num [1:4] -0.263 -0.333 0.386 -0.404
  ..$ intercepts      :List of 2
  .. ..$ predictors: num [1:4] 0 0 0 0
  .. ..$ response  : num 0
  ..$ phylo_signal    : num [1:5] 0.8 0.5 0.2 0.2 0.7
  ..$ sigmas          :List of 2
  .. ..$ sigma_species : num [1:5] 1 1 1 1 1
  .. ..$ sigma_residual: num [1:5] 1 1 1 1 1
  ..$ n_cases         : num 150
  ..$ n_species       : num 50
  ..$ n_species_actual: int 48
  ..$ birth           : num 0.8
  ..$ death           : num 0.4
  ..$ missingness     : num [1:5] 0 0.3 0.4 0.5 0.2
  ..$ rr              : logi FALSE
  ..$ rr_form         : NULL
  ..$ design          :List of 2
  .. ..$ ns       : Named num [1:4] 1 1 2 1
  .. .. ..- attr(*, "names")= chr [1:4] "gaussian" "binary" "multinomial3" "gaussian"
  .. ..$ sim_order: Named int [1:4] 1 2 3 4
  .. .. ..- attr(*, "names")= chr [1:4] "x1" "x2" "x3" "x4"
 $ random_effects:List of 4
  ..$ u_species:List of 5
  .. ..$ y : num [1:48] -0.548 -0.444 -1.554 0.371 2.194 ...
  .. ..$ x1: num [1:48] 0.5566 -2.159 1.9682 -0.0591 -2.356 ...
  .. ..$ x2: num [1:48] 1.224 -0.599 1.586 1.265 0.41 ...
  .. ..$ x3: num [1:48] -0.463 -0.426 0.226 -2.245 -0.996 ...
  .. ..$ x4: num [1:48] 0.129 -0.7956 0.2753 -0.0276 -0.9274 ...
  ..$ u_phylo  :List of 5
  .. ..$ y : Named num [1:48] -0.00973 -2.84198 -3.25706 -1.50313 2.74144 ...
  .. .. ..- attr(*, "names")= chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  .. ..$ x1: Named num [1:48] 0.43 0.78 0.566 1.779 0.201 ...
  .. .. ..- attr(*, "names")= chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  .. ..$ x2: Named num [1:48] -0.806 0.193 0.373 0.238 0.214 ...
  .. .. ..- attr(*, "names")= chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  .. ..$ x3: Named num [1:48] -1.3357 -0.436 -0.3631 -0.0673 0.9932 ...
  .. .. ..- attr(*, "names")= chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  .. ..$ x4: Named num [1:48] -2.61 1.22 -2 -1.09 -1.89 ...
  .. .. ..- attr(*, "names")= chr [1:48] "0L21U" "51eFU" "5nQQT" "7U6SR" ...
  ..$ u_slopes : list()
  ..$ residuals:List of 5
  .. ..$ y : num [1:150] -0.154 1.727 0.653 1.777 0.147 ...
  .. ..$ x1: num [1:150] 1.315 -0.828 0.835 0.511 0.532 ...
  .. ..$ x2: num [1:150] 0.9608 1.0049 -0.5292 0.3414 0.0153 ...
  .. ..$ x3: num [1:150] 0.0475 0.3536 0.4154 -0.6016 -0.1808 ...
  .. ..$ x4: num [1:150] 1.583 0.407 -0.344 -0.116 1.75 ...

[szymekdr]

Function print_simBACE_summary() gives a structured summary of simulated data:

r$> print_simBACE_summary(simdata2)
=== simBACE Simulation Summary ===

Response type: gaussian 
Predictor types: gaussian, binary, multinomial3, gaussian 
Variable names: y, x1, x2, x3, x4 

Sample size: 150 observations
Species: 48 (requested: 50 )
Tree tips: 48 

Phylogenetic signal: 0.8, 0.5, 0.2, 0.2, 0.7 
Missingness: 0 %, 30 %, 40 %, 50 %, 20 % 

Beta matrix:
          x1         x2        x3 x4
x1 0.0000000  0.0000000 0.0000000  0
x2 0.0626379  0.0000000 0.0000000  0
x3 0.0000000 -0.3269397 0.0000000  0
x4 0.4135717 -0.2585901 0.2630108  0


Beta response coefficients:
[1] -0.2631214 -0.3334791  0.3863068 -0.4043886

Random slopes: NO

--- Data preview ---
  species           y         x1 x2   x3         x4
1   eHZDN  0.11752327  1.8373702  0 <NA> -3.3161371
2   REEEZ  4.82175939         NA NA    B -4.1015317
3   CxOyx  4.23233938  2.4483857 NA    B -5.1958789
4   goBMd -0.07278349         NA  0 <NA> -0.7245938
5   qqH6v -0.50290770  2.3159161  0    A -4.4459088
6   hoK1u -4.00282560 -0.2204237 NA <NA>  1.7056161

--- Data summary ---
   species                y                 x1                x2        
 Length:150         Min.   :-6.9974   Min.   :-5.3018   Min.   :0.0000  
 Class :character   1st Qu.:-3.2301   1st Qu.:-1.2951   1st Qu.:0.0000  
 Mode  :character   Median :-1.1683   Median : 0.3149   Median :0.0000  
                    Mean   :-0.9538   Mean   : 0.2723   Mean   :0.3556  
                    3rd Qu.: 0.6742   3rd Qu.: 1.8480   3rd Qu.:1.0000  
                    Max.   : 7.7736   Max.   : 4.5284   Max.   :1.0000  
                                      NA's   :45        NA's   :60      
      x3                  x4         
 Length:150         Min.   :-6.6282  
 Class :character   1st Qu.:-2.7246  
 Mode  :character   Median :-0.9744  
                    Mean   :-0.6309  
                    3rd Qu.: 2.1738  
                    Max.   : 5.4822  
                    NA's   :30       

*** NOTE ***
Variables 'x2' are categorical but represented as integers.
Consider converting to factor if needed for analysis.

[szymekdr]

3 wrappers can be run quickly to generate simple Gaussian, binary or Poisson datasets: simBACE_binary(), simBACE_gaussian(), simBACE_poisson() - they run without any parameters.