This package aims to enable users to build and interpret multivariate machine learning models harnessing the tidyverse (tidy model syntax in particular). This package builds off ideas from Gradient Forests (Ellis et al 2012), ecological genomic approaches (Fitzpatrick and Keller, 2014) and multi-response stacking algorithms (Xing et al 2019).

This package can be of use for any multi-response machine learning problem, but was designed to handle data common to community ecology (site by species data) and ecological genomics (individual or population by SNP loci).

Installation

Install the stable version of the package:

#install.packages("devtools")
#devtools::install_github('nfj1380/mrIML')
library(mrIML)

Quick start

mrIML is designed to be used with a single function call or to be used in an ad-hoc fashion via individual function calls. In the following section we give an overview of the simple use case. For more on using each function see the function documentation. The core functions for both regression and classification are: mrIMLpredicts, mrIMLperformance, and mrInteractions,for plotting and visualization mrVip, mrFlashlight, andplot_vi. Estimating the interactions alone can be substantially computationally demanding depending on the number of outcomes you want to test. The first step to using the package is to load it as follows.

Model component

Now all the data is loaded and ready to go we can formulate the model using tidymodel syntax. In this case we have binary data (SNP presence/absence at each loci) but the data could also be counts or continuous (the set_model argument would be “regression” instead of “classification”). The user can specify any model from the tidymodel universe as ‘model 1’ (see https://www.tidymodels.org/find/ for details). However, we have done most of our testing on random forests (rf), xgr boost and glms (generalized linear models). Here we will specify a random forest classification model as the model applied to each response.

model1 <- 
  rand_forest(trees = 100, mode = "classification", mtry = tune(), min_n = tune()) %>% #100 trees are set for brevity
  set_engine("randomForest")

mrIMLpredicts

This function represents the core functionality of the package and includes results reporting, plotting and optional saving. It requires a data frame of X t( the snp data for example) and Y represented by the covariates or features.

For example, for a regression problem lets look at the data from Fitzpatrick and Keller (2015) from mrIML. Columns 1-2 represent geographical coordinates, 3-9 environmental data for each location, 10-14 are spatial eigenvectors (MEMs) and the remaining columns are proportion of the poulations with a each reference allelle.

data <- gfData[1:20]
head(data)
#>         x     y bio_1 bio_10 bio_15 bio_18 bio_2 bio_7 elevation     MEM.1      MEM.2      MEM.3
#> 1 -101.57 49.14    31    182     63    189   131   489       498  24.55053   7.444415  -9.018635
#> 2 -147.89 69.10  -107     85     60     89   110   507       400 -14.18000 -22.945314 -11.845836
#> 3 -148.51 63.39   -52     90     63    196   102   399      1040 -15.77654 -27.420568 -13.038813
#> 4  -98.35 54.51   -12    159     52    204   103   510       213  20.73482   8.715802 -19.187171
#> 5 -112.79 53.33    18    150     66    234   119   427       778  18.51313  -3.488470  24.300653
#> 6 -146.17 64.59   -43    135     65    154   109   493       460 -17.18802 -31.070507 -13.633361
#>         MEM.4 REFERENCE_X80928 REFERENCE_X172196 REFERENCE_X172210 REFERENCE_X172625
#> 1  0.06879434        0.4000000        0.03333333         0.2692308         0.4230769
#> 2 21.59719769        0.3611111        0.00000000         0.2222222         0.4705882
#> 3 14.65433486        0.4333333        0.00000000         0.1666667         0.1666667
#> 4 -5.13771830        0.3846154        0.07142857         0.4230769         0.4615385
#> 5  7.73910268        0.4285714        0.00000000         0.4000000         0.4000000
#> 6  5.78938201        0.3333333        0.00000000         0.1785714         0.3846154
#>   REFERENCE_X177996 REFERENCE_X180626 REFERENCE_X181299
#> 1        0.00000000        0.03846154        0.06666667
#> 2        0.00000000        0.05555556        0.00000000
#> 3        0.03571429        0.16666667        0.00000000
#> 4        0.10714286        0.00000000        0.00000000
#> 5        0.06666667        0.19230769        0.06666667
#> 6        0.00000000        0.08333333        0.00000000

For a classiciation problem, the data is similar except are response (by definition) is binary (i.e., is there a muutation at that locus).

# Define set of features
FeaturesnoNA<-Features[complete.cases(Features), ]
Y <- FeaturesnoNA #for simplicity
# Define set the outcomes of interst
fData <- filterRareCommon (Responsedata, lower=0.4, higher=0.7) 
X <- fData
head(X)

yhats <- mrIMLpredicts(X=X,Y=Y, model1=model1, balance_data='no', mod='classification', parallel = FALSE)

#save(yhats, file='logreg_model')
ModelPerf <- mrIMLperformance(yhats, model1, X=X) #
ModelPerf[[2]] #this measures performance across all loci.
#> [1] 0.5689655

Plotting

VI <- mrVip(yhats, Y=Y) 
plot_vi(VI=VI,  X=X,Y=Y, modelPerf=ModelPerf, cutoff= 0, plot.pca='yes') #the cutoff reduces the number of individual models printed in the second plot. 

#> Press [enter] to plot individual variable importance summaries

#> Press [enter] to plot the importance PCA plot

Effect of a feature on genetic change

We also wrap some flashlight functionality to visualize the marginal (i.e. partial dependencies) or conditional (accumulated local effects) effect of a feature on genetic change. Partial dependencies take longer to calculate and are more sensitive to correlated features

flashlightObj <- mrFlashlight(yhats, X, Y, response = "multi", model='classification')
#plot prediction scatter for all responses. Gets busy with 
plot(light_scatter(flashlightObj, v = "Forest", type = "predicted"))


#profileData_pd <- light_profile(flashlightObj,  v = "Grassland")
#mrProfileplot(profileData_pd , sdthresh =0.05) #sdthresh removes responses from the first plot that do not vary with the feature
profileData_ale <- light_profile(flashlightObj, v = "Grassland", type = "ale") #acumulated local effects
mrProfileplot(profileData_ale , sdthresh =0.01)

#the second plot is the cumulative turnover function

Interacting predictors or features

Finally, we can assess how features interact overall to shape genetic change. Be warned this is memory intensive. Future updates to this package will enable users to visualize these interactions and explore them in more detail using 2D ALE plots for example.

#interactions <-mrInteractions(yhats, X, Y,  mod='classification') #this is computationally intensive so multicores are needed. If stopped prematurely - have to reload things
#mrPlot_interactions(interactions, X,Y, top_ranking = 2, top_response=2)

References

Xing, L, Lesperance, ML and Zhang, X (2020). Simultaneous prediction of multiple outcomes using revised stacking algorithms. Bioinformatics, 36, 65-72.

Fitzpatrick, M.C. & Keller, S.R. (2015) Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecology Letters 18, 1–16.

Ellis, N., Smith, S.J. and Pitcher, C.R. (2012), Gradient forests: calculating importance gradients on physical predictors. Ecology, 93: 156-168. doi:10.1890/11-0252.1

Session info

devtools::session_info()
#> - Session info -----------------------------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.6.3 (2020-02-29)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  English_United States.1252  
#>  ctype    English_United States.1252  
#>  tz       America/New_York            
#>  date     2021-05-14                  
#> 
#> - Packages ---------------------------------------------------------------------------------------
#>  package         * version    date       lib source                        
#>  assertthat        0.2.1      2019-03-21 [1] CRAN (R 3.6.3)                
#>  backports         1.2.1      2020-12-09 [1] CRAN (R 3.6.3)                
#>  BBmisc            1.11       2017-03-10 [1] CRAN (R 3.6.3)                
#>  broom           * 0.7.6      2021-04-05 [1] CRAN (R 3.6.3)                
#>  cachem            1.0.4      2021-02-13 [1] CRAN (R 3.6.3)                
#>  callr             3.7.0      2021-04-20 [1] CRAN (R 3.6.3)                
#>  caret           * 6.0-86     2020-03-20 [1] CRAN (R 3.6.3)                
#>  cellranger        1.1.0      2016-07-27 [1] CRAN (R 3.6.3)                
#>  checkmate         2.0.0      2020-02-06 [1] CRAN (R 3.6.3)                
#>  class             7.3-15     2019-01-01 [2] CRAN (R 3.6.3)                
#>  cli               2.5.0      2021-04-26 [1] CRAN (R 3.6.3)                
#>  cluster           2.1.0      2019-06-19 [2] CRAN (R 3.6.3)                
#>  codetools         0.2-16     2018-12-24 [2] CRAN (R 3.6.3)                
#>  colorspace        2.0-0      2020-11-11 [1] CRAN (R 3.6.3)                
#>  cowplot           1.1.1      2020-12-30 [1] CRAN (R 3.6.3)                
#>  crayon            1.4.1      2021-02-08 [1] CRAN (R 3.6.3)                
#>  data.table        1.14.0     2021-02-21 [1] CRAN (R 3.6.3)                
#>  DBI               1.1.1      2021-01-15 [1] CRAN (R 3.6.3)                
#>  dbplyr            2.1.1      2021-04-06 [1] CRAN (R 3.6.3)                
#>  desc              1.3.0      2021-03-05 [1] CRAN (R 3.6.3)                
#>  devtools          2.4.0      2021-04-07 [1] CRAN (R 3.6.3)                
#>  dials           * 0.0.9      2020-09-16 [1] CRAN (R 3.6.3)                
#>  DiceDesign        1.9        2021-02-13 [1] CRAN (R 3.6.3)                
#>  digest            0.6.27     2020-10-24 [1] CRAN (R 3.6.3)                
#>  doParallel      * 1.0.16     2020-10-16 [1] CRAN (R 3.6.3)                
#>  dplyr           * 1.0.5      2021-03-05 [1] CRAN (R 3.6.3)                
#>  ellipsis          0.3.1      2020-05-15 [1] CRAN (R 3.6.3)                
#>  evaluate          0.14       2019-05-28 [1] CRAN (R 3.6.3)                
#>  extrafont         0.17       2014-12-08 [1] CRAN (R 3.6.2)                
#>  extrafontdb       1.0        2012-06-11 [1] CRAN (R 3.6.0)                
#>  fansi             0.4.2      2021-01-15 [1] CRAN (R 3.6.3)                
#>  farver            2.1.0      2021-02-28 [1] CRAN (R 3.6.3)                
#>  fastmap           1.1.0      2021-01-25 [1] CRAN (R 3.6.3)                
#>  fastmatch         1.1-0      2017-01-28 [1] CRAN (R 3.6.0)                
#>  flashlight      * 0.8.0      2021-04-21 [1] CRAN (R 3.6.3)                
#>  FNN               1.1.3      2019-02-15 [1] CRAN (R 3.6.3)                
#>  forcats         * 0.5.1      2021-01-27 [1] CRAN (R 3.6.3)                
#>  foreach         * 1.5.1      2020-10-15 [1] CRAN (R 3.6.3)                
#>  fs                1.5.0      2020-07-31 [1] CRAN (R 3.6.3)                
#>  furrr             0.2.2      2021-01-29 [1] CRAN (R 3.6.3)                
#>  future            1.21.0     2020-12-10 [1] CRAN (R 3.6.3)                
#>  gbm             * 2.1.8      2020-07-15 [1] CRAN (R 3.6.3)                
#>  gdtools           0.2.3      2021-01-06 [1] CRAN (R 3.6.3)                
#>  generics          0.1.0      2020-10-31 [1] CRAN (R 3.6.3)                
#>  ggplot2         * 3.3.3      2020-12-30 [1] CRAN (R 3.6.3)                
#>  ggrepel         * 0.9.1      2021-01-15 [1] CRAN (R 3.6.3)                
#>  globals           0.14.0     2020-11-22 [1] CRAN (R 3.6.3)                
#>  glue              1.4.2      2020-08-27 [1] CRAN (R 3.6.3)                
#>  gower             0.2.2      2020-06-23 [1] CRAN (R 3.6.3)                
#>  GPfit             1.0-8      2019-02-08 [1] CRAN (R 3.6.3)                
#>  gridExtra         2.3        2017-09-09 [1] CRAN (R 3.6.3)                
#>  gtable            0.3.0      2019-03-25 [1] CRAN (R 3.6.3)                
#>  hardhat           0.1.5      2020-11-09 [1] CRAN (R 3.6.3)                
#>  haven             2.4.0      2021-04-14 [1] CRAN (R 3.6.3)                
#>  highr             0.9        2021-04-16 [1] CRAN (R 3.6.3)                
#>  hms               1.0.0      2021-01-13 [1] CRAN (R 3.6.3)                
#>  hrbrthemes      * 0.8.0      2020-03-06 [1] CRAN (R 3.6.3)                
#>  htmltools         0.5.1.1    2021-01-22 [1] CRAN (R 3.6.3)                
#>  httr              1.4.2      2020-07-20 [1] CRAN (R 3.6.3)                
#>  infer           * 0.5.4      2021-01-13 [1] CRAN (R 3.6.3)                
#>  ipred             0.9-11     2021-03-12 [1] CRAN (R 3.6.3)                
#>  iterators       * 1.0.13     2020-10-15 [1] CRAN (R 3.6.3)                
#>  janitor         * 2.1.0      2021-01-05 [1] CRAN (R 3.6.3)                
#>  jsonlite          1.7.2      2020-12-09 [1] CRAN (R 3.6.3)                
#>  knitr             1.32       2021-04-14 [1] CRAN (R 3.6.3)                
#>  labeling          0.4.2      2020-10-20 [1] CRAN (R 3.6.3)                
#>  lattice         * 0.20-38    2018-11-04 [2] CRAN (R 3.6.3)                
#>  lava              1.6.9      2021-03-11 [1] CRAN (R 3.6.3)                
#>  lhs               1.1.1      2020-10-05 [1] CRAN (R 3.6.3)                
#>  lifecycle         1.0.0      2021-02-15 [1] CRAN (R 3.6.3)                
#>  listenv           0.8.0      2019-12-05 [1] CRAN (R 3.6.3)                
#>  lubridate         1.7.10     2021-02-26 [1] CRAN (R 3.6.3)                
#>  magrittr          2.0.1      2020-11-17 [1] CRAN (R 3.6.3)                
#>  MASS              7.3-51.5   2019-12-20 [2] CRAN (R 3.6.3)                
#>  Matrix            1.2-18     2019-11-27 [2] CRAN (R 3.6.3)                
#>  memoise           2.0.0      2021-01-26 [1] CRAN (R 3.6.3)                
#>  MetricsWeighted   0.5.2      2021-04-16 [1] CRAN (R 3.6.3)                
#>  mgcv              1.8-31     2019-11-09 [2] CRAN (R 3.6.3)                
#>  mlr               2.19.0     2021-02-22 [1] CRAN (R 3.6.3)                
#>  modeldata       * 0.1.0      2020-10-22 [1] CRAN (R 3.6.3)                
#>  ModelMetrics      1.2.2.2    2020-03-17 [1] CRAN (R 3.6.3)                
#>  modelr            0.1.8      2020-05-19 [1] CRAN (R 3.6.3)                
#>  mrIML           * 1.0.0      2021-05-13 [1] Github (nfj1380/mrIML@7d971b9)
#>  munsell           0.5.0      2018-06-12 [1] CRAN (R 3.6.3)                
#>  nlme              3.1-144    2020-02-06 [2] CRAN (R 3.6.3)                
#>  nnet              7.3-12     2016-02-02 [2] CRAN (R 3.6.3)                
#>  parallelly        1.24.0     2021-03-14 [1] CRAN (R 3.6.3)                
#>  parallelMap       1.5.0      2020-03-26 [1] CRAN (R 3.6.3)                
#>  ParamHelpers      1.14       2020-03-24 [1] CRAN (R 3.6.3)                
#>  parsnip         * 0.1.5      2021-01-19 [1] CRAN (R 3.6.3)                
#>  pdp             * 0.7.0      2018-08-27 [1] CRAN (R 3.6.3)                
#>  permute         * 0.9-5      2019-03-12 [1] CRAN (R 3.6.3)                
#>  pillar            1.6.0      2021-04-13 [1] CRAN (R 3.6.3)                
#>  pkgbuild          1.2.0      2020-12-15 [1] CRAN (R 3.6.3)                
#>  pkgconfig         2.0.3      2019-09-22 [1] CRAN (R 3.6.3)                
#>  pkgload           1.2.1      2021-04-06 [1] CRAN (R 3.6.3)                
#>  plyr              1.8.6      2020-03-03 [1] CRAN (R 3.6.3)                
#>  prettyunits       1.1.1      2020-01-24 [1] CRAN (R 3.6.3)                
#>  pROC              1.17.0.1   2021-01-13 [1] CRAN (R 3.6.3)                
#>  processx          3.5.0      2021-03-23 [1] CRAN (R 3.6.3)                
#>  prodlim           2019.11.13 2019-11-17 [1] CRAN (R 3.6.3)                
#>  ps                1.6.0      2021-02-28 [1] CRAN (R 3.6.3)                
#>  purrr           * 0.3.4      2020-04-17 [1] CRAN (R 3.6.3)                
#>  R6                2.5.0      2020-10-28 [1] CRAN (R 3.6.3)                
#>  randomForest    * 4.6-14     2018-03-25 [1] CRAN (R 3.6.3)                
#>  RANN              2.6.1      2019-01-08 [1] CRAN (R 3.6.3)                
#>  Rcpp              1.0.6      2021-01-15 [1] CRAN (R 3.6.3)                
#>  readr           * 1.4.0      2020-10-05 [1] CRAN (R 3.6.3)                
#>  readxl            1.3.1      2019-03-13 [1] CRAN (R 3.6.3)                
#>  recipes         * 0.1.16     2021-04-16 [1] CRAN (R 3.6.3)                
#>  remotes           2.3.0      2021-04-01 [1] CRAN (R 3.6.3)                
#>  reprex            2.0.0      2021-04-02 [1] CRAN (R 3.6.3)                
#>  reshape2          1.4.4      2020-04-09 [1] CRAN (R 3.6.3)                
#>  rlang           * 0.4.10     2020-12-30 [1] CRAN (R 3.6.3)                
#>  rmarkdown         2.7        2021-02-19 [1] CRAN (R 3.6.3)                
#>  ROSE              0.0-3      2014-07-15 [1] CRAN (R 3.6.3)                
#>  rpart             4.1-15     2019-04-12 [2] CRAN (R 3.6.3)                
#>  rpart.plot        3.0.9      2020-09-17 [1] CRAN (R 3.6.3)                
#>  rprojroot         2.0.2      2020-11-15 [1] CRAN (R 3.6.3)                
#>  rsample         * 0.0.9      2021-02-17 [1] CRAN (R 3.6.3)                
#>  rstudioapi        0.13       2020-11-12 [1] CRAN (R 3.6.3)                
#>  Rttf2pt1          1.3.8      2020-01-10 [1] CRAN (R 3.6.2)                
#>  rvest             1.0.0      2021-03-09 [1] CRAN (R 3.6.3)                
#>  scales          * 1.1.1      2020-05-11 [1] CRAN (R 3.6.3)                
#>  sessioninfo       1.1.1      2018-11-05 [1] CRAN (R 3.6.3)                
#>  snakecase         0.11.0     2019-05-25 [1] CRAN (R 3.6.3)                
#>  stringi           1.5.3      2020-09-09 [1] CRAN (R 3.6.3)                
#>  stringr         * 1.4.0      2019-02-10 [1] CRAN (R 3.6.3)                
#>  survival          3.1-8      2019-12-03 [2] CRAN (R 3.6.3)                
#>  systemfonts       1.0.1      2021-02-09 [1] CRAN (R 3.6.3)                
#>  testthat          3.0.2      2021-02-14 [1] CRAN (R 3.6.3)                
#>  themis          * 0.1.3      2020-11-12 [1] CRAN (R 3.6.3)                
#>  tibble          * 3.1.0      2021-02-25 [1] CRAN (R 3.6.3)                
#>  tidymodels      * 0.1.3      2021-04-19 [1] CRAN (R 3.6.3)                
#>  tidyr           * 1.1.3      2021-03-03 [1] CRAN (R 3.6.3)                
#>  tidyselect        1.1.0      2020-05-11 [1] CRAN (R 3.6.3)                
#>  tidyverse       * 1.3.1      2021-04-15 [1] CRAN (R 3.6.3)                
#>  timeDate          3043.102   2018-02-21 [1] CRAN (R 3.6.3)                
#>  tune            * 0.1.4      2021-04-20 [1] CRAN (R 3.6.3)                
#>  unbalanced        2.0        2015-06-26 [1] CRAN (R 3.6.3)                
#>  usethis           2.0.1      2021-02-10 [1] CRAN (R 3.6.3)                
#>  utf8              1.2.1      2021-03-12 [1] CRAN (R 3.6.3)                
#>  vctrs           * 0.3.6      2020-12-17 [1] CRAN (R 3.6.3)                
#>  vegan           * 2.5-7      2020-11-28 [1] CRAN (R 3.6.3)                
#>  vip             * 0.3.2      2020-12-17 [1] CRAN (R 3.6.3)                
#>  viridis         * 0.6.0      2021-04-15 [1] CRAN (R 3.6.3)                
#>  viridisLite     * 0.4.0      2021-04-13 [1] CRAN (R 3.6.3)                
#>  withr             2.4.2      2021-04-18 [1] CRAN (R 3.6.3)                
#>  workflows       * 0.2.2      2021-03-10 [1] CRAN (R 3.6.3)                
#>  workflowsets    * 0.0.2      2021-04-16 [1] CRAN (R 3.6.3)                
#>  xfun              0.22       2021-03-11 [1] CRAN (R 3.6.3)                
#>  xgboost         * 1.3.2.1    2021-01-18 [1] CRAN (R 3.6.3)                
#>  xml2              1.3.2      2020-04-23 [1] CRAN (R 3.6.3)                
#>  yaml              2.2.1      2020-02-01 [1] CRAN (R 3.6.3)                
#>  yardstick       * 0.0.8      2021-03-28 [1] CRAN (R 3.6.3)                
#> 
#> [1] C:/Users/gmachad/Documents/R/win-library/3.6
#> [2] C:/Program Files/R/R-3.6.3/library