Introduction: The MXM R package for Feature Selection

The MXM R Package, short for the latin 'Mens ex Machina' ( Mind from the Machine ), is a collection of utility functions for feature selection, cross validation and Bayesian Networks. The package supports conditional independence tests for various combinations of target and predictor variables (continuous, categorical). MXM offers many feature selection algorithms focused on providing one or more minimal feature subsets, refered also as variable signatures, that can be used to improve the performance of downstream analysis tasks such as regression and classification, by excluding irrelevant and redundant variables. In this tutorial we will learn how to use the MMPC algorithm.

Overview of the MMPC algorithm

MMPC stands for Max-Min Parents and Children, a constraint based feature selection algorithm, as first described by Brown, Tsamardinos and Aliferis, (2004). Parents and Children refers to the fact that the algorithm identifies the parents and children of the variable of interest ( target ), assuming a Bayesian Network for all observed variables. What it will not recover is the spouses of the children , and for this the FBED algorithmcan be applied The FBED algorithm is also available in the MXM R package and can essentially recover the Markov Blanket of the variable of interest. For simplicity, we will use a dataset with fewer variables, but the algorithms perform especially well in datasets with very large feature space, such as in biomedical datasets (eg. millions of SNPs as variables in GWAS studies, thousands of genes in NGS omics datasets, etc).

Selecting the appropriate Conditional Independence Test

At its core, the MMPC algorithm performs multiple conditional independance tests, and progressively excludes irrelevant and/or redundant variables. The final variables that have “survived” through all those elimination stages, are the MMPC output signature.

The selection of the appropriate conditional independence test is a crucial decision for the validity and success of downstream statistical analysis and machine learning tasks. Currently the MXM R package supports nummerous tests for different combinations of target ( dependent ) and predictor ( independent ) variables. A detailed summary table to guide you through the selection of the most suitable test is included in MXM's reference manual ( “CondInditional independence tests” )

The MMPC() function: Required and optional arguments

There are 3 mandatory arguments for the MXM::MMPC() function: 1) an object with the target variable, 2) one with the complete dataset but with the target variable removed and 3) a conditional indepence test, selected by the reference manual table mentioned above. The dataset has to have the instances (N) as rows and the features (f) as columns. For the target and dataset it is recommended to also retain colnames and rownames information.

Several hyper-parameters are also provided as optional arguments in the MXM::MMPC() function. In the following block, the function along with the most important hyperparameters are presented.

# Overview the MXM::`MMPC()` function 

mod <- MXM::MMPC(
       target,           # The target variable vector
       dataset,          # The dataset with the target column removed
       max_k = 3,        # The maximum size of the conditioning set to use 
       threshold = 0.05, # level of alpha for statistical  significance
       test = 'testIndFisher',   
       ini = NULL,        # if TRUE, the calculated univariate associations
                          # are stored for runtime efficiency in subsequent 
                          # MMPC runs with diferent hyper-parameters.
       hash =  TRUE,      # if TRUE, the calculated statistics are stored.
       hashObject = NULL, # the mmpcobject from a previous run
       ncores = 1,        # number of cores for parallel execution. 
                          # Recommended for thousands of variables.
       backward = TRUE)   # If TRUE, the backward phase
                          # (or symmetry correction) is implemented.
                          # Falsely included variables,
                          # in the MMPC output signature are removed. 



Applying the MMPC algorithm on the UCI wine dataset

For this tutorial we will use the UCI wine dataset. The dataset contains the results of a chemical analysis performed on 3 different types of wines and includes 12 quality related characteristics plus the information of the wine class as the first attribute (Type: 1,2 or 3). More information about the wine dataset is available at the UCI repository.

The following block installs the MXM package, takes care of package dependencies, downloads and then cleans the dataset for the subsequent steps. The categorical variable (class information) is omitted for this example and we retain only the numerical variables (continuous and count data). We will then apply the MMPC algorithm to acquire a minimal, highly relevant subset of variables that can be used to best model the "NonFlavanoids" content.

PACKAGE DEPENDENCIES and UCI WINE DATASET LOADING:

# 0. INSTALL and LOAD the MXM R Package:
#install.packages('MXM', dependencies = TRUE )
library(MXM)

# 1. DOWNLOAD the wine dataset from UCI:
URL  <- "ftp://ftp.ics.uci.edu/pub/machine-learning-databases/wine/wine.data"
wine <-  read.csv(URL, header = FALSE) 


# 2. SET variables' names as header:
colnames(wine) <- c('Type', 'Alcohol', 'Malic', 'Ash', 
                    'Alcalinity', 'Magnesium', 'Phenols', 
                    'Flavanoids', 'Nonflavanoids', 'Proanthocyanins',
                    'Color', 'Hue', 'Dilution', 'Proline')

# 3. REMOVE the 1st attribute, which is the class information:
wine <- wine[,-1] 

# 4. PREVIEW UCI's wine dataset:
head(wine, 2)
##   Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids
## 1   14.23  1.71 2.43       15.6       127    2.80       3.06          0.28
## 2   13.20  1.78 2.14       11.2       100    2.65       2.76          0.26
##   Proanthocyanins Color  Hue Dilution Proline
## 1            2.29  5.64 1.04     3.92    1065
## 2            1.28  4.38 1.05     3.40    1050
# The header should include the wine attributes sans the class labels, 
# in the following order:

# Alcohol | Malic | Ash | Alcalinity | Magnesium | Phenols | Flavanoids 
# Nonflavanoids | Proanthocyanins | Color | Hue | Dilution | Proline 

Exploratory Data Analysis: Inspecting the wine dataset

# 5. CHECK for missing or non-numeric values in the dataframe:

sum(is.na(wine))
## [1] 0
sum(is.nan(as.matrix(wine)))  #if 0, then No NAs, none NaNs, good to go!
## [1] 0

Even if the dataset contains missing values, they will automatically be replaced by the current variable (column) mean value with an appropriate warning to the user after the execution. For optimal results, it is advised to use a more sophisticated imputation method according to your data's needs before running MMPC().

# 6. CHECK `wine` object's data type, dimensions:
str(wine)
## 'data.frame':    178 obs. of  13 variables:
##  $ Alcohol        : num  14.2 13.2 13.2 14.4 13.2 ...
##  $ Malic          : num  1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...
##  $ Ash            : num  2.43 2.14 2.67 2.5 2.87 2.45 2.45 2.61 2.17 2.27 ...
##  $ Alcalinity     : num  15.6 11.2 18.6 16.8 21 15.2 14.6 17.6 14 16 ...
##  $ Magnesium      : int  127 100 101 113 118 112 96 121 97 98 ...
##  $ Phenols        : num  2.8 2.65 2.8 3.85 2.8 3.27 2.5 2.6 2.8 2.98 ...
##  $ Flavanoids     : num  3.06 2.76 3.24 3.49 2.69 3.39 2.52 2.51 2.98 3.15 ...
##  $ Nonflavanoids  : num  0.28 0.26 0.3 0.24 0.39 0.34 0.3 0.31 0.29 0.22 ...
##  $ Proanthocyanins: num  2.29 1.28 2.81 2.18 1.82 1.97 1.98 1.25 1.98 1.85 ...
##  $ Color          : num  5.64 4.38 5.68 7.8 4.32 6.75 5.25 5.05 5.2 7.22 ...
##  $ Hue            : num  1.04 1.05 1.03 0.86 1.04 1.05 1.02 1.06 1.08 1.01 ...
##  $ Dilution       : num  3.92 3.4 3.17 3.45 2.93 2.85 3.58 3.58 2.85 3.55 ...
##  $ Proline        : int  1065 1050 1185 1480 735 1450 1290 1295 1045 1045 ...
# The output should be a datarame: 
#'data.frame':  178 obs. of  13 variables

Preparing the MMPC() input objects: 'target', 'dataset'

And now, we will tailor the target ( non-flavanoids content variable ) and the complete dataset objects to the the MMPC()'s function needs. Both will be converted to matrices, using the built in function as.matrix() and we will make sure that the dataset matrix is given as instances (N) by features (f). After selecting the target variable, create a matrix that includes only the remaining variables. This would be the dataset input for the MMPC() function. This is necessary to assure that the signature does not include the target variable.

# 0. Exclude target variable column
targetVariable <- wine$Nonflavanoids
targetVariable <- NULL


# 1. Convert dataframe to matrix:
wine_dataset <- as.matrix(wine[, -8])
wine_dataset[, 12] <- as.numeric(wine_dataset[, 12])
head(wine_dataset, 2)
##      Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids
## [1,]   14.23  1.71 2.43       15.6       127    2.80       3.06
## [2,]   13.20  1.78 2.14       11.2       100    2.65       2.76
##      Proanthocyanins Color  Hue Dilution Proline
## [1,]            2.29  5.64 1.04     3.92    1065
## [2,]            1.28  4.38 1.05     3.40    1050

We check once more the dimension of our dataset object. We expect it to be in the form Instances x Features.

# 2. Check dimensions of the wine_dataset
# REMINDER: We need it as N x f    // N for instances, f or features

dim(wine_dataset)
## [1] 178  12
# The output should be 178 x 12, 
#178 instances and 12 features; if so, we're good to go

If desired, you can allocate the target variable in a new data structure, merely for the purpose of keeping functions' input easy to read.

# 3. Select the target variable (`Nonflavanoids`) and store as a matrix:
target_NonFlav <- as.vector(wine$Nonflavanoids)
str(target_NonFlav,2)
##  num [1:178] 0.28 0.26 0.3 0.24 0.39 0.34 0.3 0.31 0.29 0.22 ...

After the brief EDA (Exploratory Data Analysis) of the UCI wine dataset, it's time to actually apply MMPC. For a comprehensive overview of the hyper-parameter options type ?MMPC in your Rstudio console. Here, we have set hash = TRUE, for faster runtimes in subsequent runs of the MMPC with different hyper-parameters. We have selected the testIndFisher test which is the appropriate for our data, that we have retrieved from the Conditional Independece Tests cheatsheet from the MXM reference manual. The backward = TRUE was also used, for acquiring a minimal signature. All other parameters were left to default or the first run.

First run of MMPC()

# MMPC on the wine dataset: 

library('MXM')
mmpcobject_wine_NonFlav <- MXM::MMPC( target  = target_NonFlav,            
                                       dataset = wine_dataset,            
                                       max_k = 3,          
                                       threshold = 0.05,                                         
                                       test = 'testIndFisher',   
                                       ini = NULL,                                                
                                       hash =  TRUE,      
                                       hashObject = NULL,                                        
                                       ncores = 1,         
                                       backward = TRUE)   

MMPC() output: The mmpcobject and the feature signature

Let's now explore the output of the MMPC function, the mmpcobject. The generic function summary() can be used to display the indices of the features that are included in the signature, and also contains information about the MMPC() run, such as the selected hyperparameters, execution runtime and also statistics about the distribution of the p-values from the performed tests. For more specific information, you can access the mmpcobject fields with the @ operator, the operator for accessing S4 class objects in R . This is essentially the same as using the dollar operator $ for accessing R5 class objects ' slots, typically used with dataframes for example. The fields contain the output results of the MMPC() run and also two lists that can be re-used for subsequent MMPC runs for computation time efficiency.

Below, we can see the two objects that facilitate the computation time efficiency in the MMPC() re-runs. After running MMPC_ with some hyper-parameters you might want to run the algorithm again with different hyper-parameters (max_k for example). To avoid calculating the univariate associations (first step of MMPC) again, you can take the list univ from the first run and provide it as input to the argument ini in the subsequent runs the algorithm. This can speed up the second run (and subequent runs of course) up to 50%, which is crucial if you are handling datasets with a very high number of features.

# Cache of the stats calculated in the MMPC run
str(mmpcobject_wine_NonFlav@hashObject) 
## List of 2
##  $ stat_hash  :Class 'Hash' <environment: 0x0000000025855b40> 
##  $ pvalue_hash:Class 'Hash' <environment: 0x0000000025855788>
# a list with the univariate associations
str(mmpcobject_wine_NonFlav@univ)        
## List of 2
##  $ stat  : num [1:12] -2.08 3.99 2.49 5.01 -3.47 ...
##  $ pvalue: num [1:12] -3.24 -9.25 -4.3 -13.56 -7.32 ...

Let's also save the runtime of the first MMPC run in a variable, to compare later.

execution_time_1st_MMPC_run <- mmpcobject_wine_NonFlav@runtime
execution_time_1st_MMPC_run
##    user  system elapsed 
##    0.01    0.00    0.02

MMPC() subsequent run: Efficiency with the ini and hashObject arguments

Now, we can run MMPC() again to check how we can use the hashObject and univ lists to maximize runtime efficiency when multiple runs are requirted.

# MMPC on the wine dataset: 

library('MXM')
mmpcobject_2nd_run <- MXM::MMPC(target  = target_NonFlav,            
                                 dataset = wine_dataset  ,            
                           # it was set to 3 in the 1st run 
                           max_k = 5, 
                           # it was set to 0.05 in the 1st run
                           threshold = 0.01, 
                           test = 'testIndFisher',
                           #the cached univariate tests
                           ini = mmpcobject_wine_NonFlav@univ,
                           # cached stats, p-values
                           hashObject = mmpcobject_wine_NonFlav@hashObject)  

We used the stored stats, p-values and univariate tests performed in the first run for avoiding redundant calculations. Our example dataset is very small to highlight the impact of such an implementation in the algorithm, but let's compare the runtimes for first and second run of MMPC().

execution_time_2nd_MMPC_run <- mmpcobject_2nd_run@runtime

execution_time_1st_MMPC_run
##    user  system elapsed 
##    0.01    0.00    0.02
execution_time_2nd_MMPC_run
##    user  system elapsed 
##       0       0       0

Even in our small wine dataset, the difference in runtime is impressive. MMPC is designed for, and actually shines in high feature space datasets, such as those in the domains of computer vision and -omics approaches in Life Sciences.

Grid search for hyperparameter tuning: the mmpc.path() function

In the spirit of automation and performance efficiency, it would be an impossible task to manually search for the optimal set of hyper-parameters. Thus, the MXM package supports a grid search function, named mmpc.path. The function returns an object that includes matrices with the following information that can be used for selecting the optimal configuration.

bic: matrix with the BIC values of the final fitted model based on the selected variables identified by each combination of the hyper-parameters.
size:   matrix with the legnth of the selected variables identified by each configuration of MMPC.
variables: A list containing the variables from each configuration of MMPC
runtime: The run time of the function

The desired hyperparameters to be checked can be given as vectors in the relative argument.

# Grid Search for MMPC hyper-parameter tuning 

library('MXM')
mmpcGridSearch <- MXM::mmpc.path(target  = target_NonFlav,            
                  dataset = wine_dataset,            
                  max_ks = c(3,4,5,6),  # a vector of k to try
                  alphas = NULL,   # a vector of thresholds; 
                                        # If NULL, 0.1, 0.05 and 0.01 
                                        # will be tested.
                  test = 'testIndFisher',   
                  ncores = 1)

We can acces the mmpc.path() objects using the dollar $ operator and then use the generic which function to retrieve the lowest value for BIC and the repsective hyper-parameters.

BIC_results <- as.data.frame(mmpcGridSearch$bic)
head(BIC_results, 4)
##              max_k=6   max_k=5   max_k=4   max_k=3
## alpha=0.1  -288.7666 -288.7666 -288.7666 -288.7666
## alpha=0.05 -288.7666 -288.7666 -288.7666 -288.7666
## alpha=0.01 -286.5281 -286.5281 -286.5281 -286.5281
# We can retrieve the indices of the minimum BIC values:
which(BIC_results == min(BIC_results), arr.ind = TRUE)
##            row col
## alpha=0.1    1   1
## alpha=0.05   2   1
## alpha=0.1    1   2
## alpha=0.05   2   2
## alpha=0.1    1   3
## alpha=0.05   2   3
## alpha=0.1    1   4
## alpha=0.05   2   4

Above we observe that the highest alpha level, 0.1, is the one with the lowest BIC. However, since the difference is miniscule in the BIC change, we will select a lower alpha level, to avoid including false positives in our model.

We can display also the size of the selected signatures, which is the number variables the signature of ech configuration contains.

size_of_signature_results <- as.data.frame(mmpcGridSearch$size)
head(size_of_signature_results, 4)
##            max_k=6 max_k=5 max_k=4 max_k=3
## alpha=0.1        4       4       4       4
## alpha=0.05       4       4       4       4
## alpha=0.01       2       2       2       2
# We can retrieve the indices of the maximum subset:
which(size_of_signature_results == max(size_of_signature_results), arr.ind = TRUE)
##            row col
## alpha=0.1    1   1
## alpha=0.05   2   1
## alpha=0.1    1   2
## alpha=0.05   2   2
## alpha=0.1    1   3
## alpha=0.05   2   3
## alpha=0.1    1   4
## alpha=0.05   2   4

We can also preview, the indices of the actual variables that were included in the signature in each configuration:

head(mmpcGridSearch$variables, 4)
## $`alpha=0.1 & max_k=6`
## [1]  4  5  7 11
## 
## $`alpha=0.1 & max_k=5`
## [1]  4  5  7 11
## 
## $`alpha=0.1 & max_k=4`
## [1]  4  5  7 11
## 
## $`alpha=0.1 & max_k=3`
## [1]  4  5  7 11

Let's now get back at our initial MMPC() run, and explore the output summary:

summary(mmpcobject_wine_NonFlav)
##     Length      Class       Mode 
##          1 MMPCoutput         S4

We can retrieve the selected MMPC signature, by choosing selectedVarsOrder from the MMPC output object; this way, the features are printed based on highest statistical significance.

mmpcobject_wine_NonFlav@selectedVarsOrder
## [1]  7  4  5 11
# The signature should include the variables with indices 7, 4, 5

We can retrieve the names of the features in the signature by selecting the colnames with the above indices from the dataset provided as input in the MMPC() function.

colnames(wine_dataset)[7]
## [1] "Flavanoids"
colnames(wine_dataset)[4]
## [1] "Alcalinity"
colnames(wine_dataset)[5]
## [1] "Magnesium"

The MMPC() signature, highlights the "Flavanoids", "Alcalinity" and " Magnesium" content as variables of high importance in relation to our selected target variable "NonFlavanoids" content.

We can actually create a regression model, using the above signature variables and check how the model performs. We will call the mmpc.model function and use the mmpcObject from the MMPC() run above, as input. For a more detailed overview the the function mmpc.model(), you can type"?mmpc.model" in your rstudio console.

# MODEL ESTIMATES USING MMPC'S FEATURE SUBSET AS PrEDICTORs 
mmpcmodel_wine_NonFlav<- mmpc.model(
                                   target = target_NonFlav, 
                                   dataset = wine_dataset,
                                   wei = NULL, 
                                   mmpcObject = mmpcobject_wine_NonFlav, 
                                   test = 'testIndFisher')

summary(mmpcmodel_wine_NonFlav) ; 
##           Length Class  Mode   
## mod       12     lm     list   
## signature  5     -none- numeric
mmpcmodel_wine_NonFlav$ypografi
## NULL

The MMPC algorithm follows a forward-backward filter approach for feature selection in order to provide a minimal, highly-predictive feature subset of a high dimensional dataset. The max_k hyper-parameter dictates the maximum number of variables as a conditioning set to use in the conditioning independence test. In our example, 12 model where built in the process, amongst them also the one with the final selected features as variables.

Above, the ypografi variable denotes the indices of the MMPC selected features and also the BIC (Bayesian Information Criterion) of the final model. If you are not familiar with BIC as a model selection criterion, you can type ?BIC in your Rstudio console for a brief introduction. As a rough estimate, when comparing models of the same Y as target variable, the model with the lowest BIC is prefered. Below, we will retrieve the beta Coefficients and the intercept for the selected model, so that we can write the actual formula for the regression model, for the "NonFlavanoids" as target variable.

mmpcmodel_wine_NonFlav$mod
## 
## Call:
## lm(formula = target ~ ., data = as.data.frame(dataset[, signature]), 
##     weights = wei)
## 
## Coefficients:
## (Intercept)   Flavanoids   Alcalinity    Magnesium     Dilution  
##    0.549466    -0.029591     0.007240    -0.001543    -0.043972

Let the "NonFlavanoids" target variable, be Y, then the formula for the model can be written as follows:

Y = -0.087194 + 0.032057 Alcalinity -0.006664 Magnesium -0.236471 Flavanoids

Optimized performance tips and tricks

Permutation option

The MXM package offers a permutation version of the statistical tests, which is recommended for small sample size datasets. If you are working with very large datasets, in order to save computational time, there is a trick to avoind doing all permutations. As soon as the number of times the permuted test statistic is more than the observed test statistic is more than 50 (if threshold = 0.05 and R = 999), the p-value has exceeded the signifiance level (threshold value) and hence the predictor variable is not significant. There is no need to continue do the extra permutations, as a decision has already been made.

The function perm.MMPC()has the same arguments as theMMPC()function, with an additional option calledR which is the number of permutations.

# TESTS WITH PERMUTATIONS:  

library('MXM')
permutation_MMPC_model <- MXM::perm.mmpc (target_NonFlav,
                                          wine_dataset,  
                                          R = 999, # Number of permutations 
                                          max_k      = 3, 
                                          threshold  = 0.05, 
                                          test       = 'permFisher', 
                                          ncores     = 1,
                                          backward   = FALSE)

Choosing the appropriate max_k for the conditioning set

As a rule of thumb, when the sample size is rather small eg. < 100, the default max_k = 3 is recommended. In small feature space datasets, as for example in the case of the wine dataset the default option is also recommended. When working with high dimensional datasets (many hundreds to thousands of features), the following approximation can be used:

max_k = floor(N/10) where N is the number of instances.

For example, in a gene expression dataset, where N = 100 and features = 50000, max_k = 100/10 = 10, could be used. For more information about this heuristic there is an interesting topic in stackoverflow http://bit.ly/stack_ovrflw_max_predictors_in_multi_regression.

Parallel computing: Choose number of cores option

You can use more than one cores to speed up execution time. However, this only pays off when working with large datasets. For small datasets, like in our example ncore = 1 is recommended.

You are now ready to apply the MMPC() algorithm and explore your very own dataset! The MXM is under intensive development and will be regularly updated with new functionalities. For questions, algorithm requests and suggestions, do not hesitate to contact us at mtsagris@uoc.gr.

References

Lagani, V., Athineou, G., Farcomeni, A., Tsagris, M. & Tsamardinos, I. Feature Selection with the R Package MXM: Discovering Statistically Equivalent Feature Subsets. J. Stat. Softw. 80, 1-25 (2017).

Borboudakis, G. & Tsamardinos, I. Forward-Backward Selection with Early Dropping. (2017).

Tsamardinos, Brown and Aliferis (2006). The max-min hill-climbing Bayesian network structure learning algorithm. Machine learning, 65(1), 31-78.

Brown, L. E., Tsamardinos, I., & Aliferis, C. F. (2004). A novel algorithm for scalable and accurate Bayesian network learning. Medinfo, 711-715.

Tsamardinos, I., Aliferis, C. F., & Statnikov, A. (2003). Time and sample efficient discovery of Markov blankets and direct causal relations. In Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining (pp. 673-678). ACM.

Statnikov, A. R., Tsamardinos, I., & Aliferis, C. F. (2003). An Algorithm For Generation of Large Bayesian Networks.

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                  LC_CTYPE=Greek_Greece.1253   
## [3] LC_MONETARY=Greek_Greece.1253 LC_NUMERIC=C                 
## [5] LC_TIME=Greek_Greece.1253    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_0.8.3 MXM_1.4.5  
## 
## loaded via a namespace (and not attached):
##  [1] foreach_1.4.4       jsonlite_1.6        splines_3.6.0      
##  [4] R.utils_2.8.0       Formula_1.2-3       ucminf_1.1-4       
##  [7] assertthat_0.2.1    latticeExtra_0.6-28 slam_0.1-45        
## [10] numDeriv_2016.8-1   pillar_1.4.0        backports_1.1.4    
## [13] lattice_0.20-38     quantreg_5.38       glue_1.3.1         
## [16] digest_0.6.18       RColorBrewer_1.1-2  checkmate_1.9.4    
## [19] minqa_1.2.4         colorspace_1.4-1    htmltools_0.3.6    
## [22] Matrix_1.2-17       R.oo_1.22.0         Rfast2_0.0.4       
## [25] pkgconfig_2.0.2     SparseM_1.77        purrr_0.3.2        
## [28] relations_0.6-8     scales_1.0.0        lme4_1.1-21        
## [31] MatrixModels_0.4-1  htmlTable_1.13.2    tibble_2.1.1       
## [34] ggplot2_3.2.0       geepack_1.2-1       nnet_7.3-12        
## [37] lazyeval_0.2.2      survival_2.44-1.1   magrittr_1.5       
## [40] crayon_1.3.4        evaluate_0.13       ordinal_2019.4-25  
## [43] R.methodsS3_1.7.1   doParallel_1.0.14   R.cache_0.13.0     
## [46] nlme_3.1-139        MASS_7.3-51.4       foreign_0.8-71     
## [49] R.rsp_0.43.1        RcppZiggurat_0.1.5  tools_3.6.0        
## [52] coxme_2.2-10        data.table_1.12.2   energy_1.7-5       
## [55] stringr_1.4.0       munsell_0.5.0       cluster_2.0.8      
## [58] Rfast_1.9.6         compiler_3.6.0      rlang_0.4.0        
## [61] grid_3.6.0          nloptr_1.2.1        iterators_1.0.10   
## [64] rstudioapi_0.10     bigmemory.sri_0.1.3 htmlwidgets_1.3    
## [67] visNetwork_2.0.6    base64enc_0.1-3     boot_1.3-22        
## [70] codetools_0.2-16    gtable_0.3.0        sets_1.0-18        
## [73] markdown_0.9        bigmemory_4.5.33    R6_2.4.0           
## [76] gridExtra_2.3       knitr_1.22          bdsmatrix_1.3-3    
## [79] zeallot_0.1.0       Hmisc_4.2-0         stringi_1.4.3      
## [82] parallel_3.6.0      Rcpp_1.0.1          vctrs_0.2.0        
## [85] rpart_4.1-15        acepack_1.4.1       tidyselect_0.2.5   
## [88] xfun_0.7