inst/doc/ELAplus-intro.R

In ELAplus: Energy Landscape Analysis for Ecological Communities

## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
opts_chunk$set(
  fig.align = "center",
  out.extra = 'style="max-width:100%; overflow-x:auto; white-space: nowrap;"',
  tidy = FALSE
)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  message  = FALSE,
  warning  = FALSE
)

## -----------------------------------------------------------------------------
# load package
library(ELAplus)
# load test data set
# species abundance table
data("baseabtable")
head(baseabtable)

# metadata (environmental factors)
data("basemetadata")
head(basemetadata)

## -----------------------------------------------------------------------------
# load package
library(ELAplus)
# Generate random parameters
hb_params <- hb.paramgen(24, ne = 2) # the number of environmental factors
h.act <- hb_params[[1]]
j.act <- hb_params[[2]]
g.act <- hb_params[[3]]

# Simulate community composition data
rand_data <- HeatBath(100, 512, h.act, j.act, g = g.act)
rand_mat   <- rand_data[[1]]
rand_enmat <- rand_data[[2]]


## -----------------------------------------------------------------------------
formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata)

## ----eval=FALSE---------------------------------------------------------------
# formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata,
#                         normalize = 1, parameters = c(0.01, 0.01, 0.99))
# #> Processed 256 samples.
# #> Relative abundance threshold = 0.01
# #> Occurrence threshold (lower) = 0.01
# #> Occurrence threshold (upper) = 0.99
# #> Selected  16  out of  16 species.
# ocvecs <- formatted[[1]]  # community composition binary matrix
# abvecs <- formatted[[2]]  # community abundance matrix (not binarized)
# envecs <- formatted[[3]]  # environmental variables (optional)

## ----echo=FALSE, out.width="80%"----------------------------------------------
knitr::include_graphics("figures/ParameterMatrix.png")

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/EPMEM_abstract2.png")

## ----eval=FALSE---------------------------------------------------------------
# bpresult <- Findbp(ocvecs, rep = 2, threads = 2,
#                    cvmode="10fold",nrepeat = 10, we=c(0.001),totalit=4000,
#                    lmd=c(1e-05,1e-03,0.01),maxlrs = c(0.005),
#                    tol=0.03, runadamW=TRUE,sparse=TRUE,fastfitting=TRUE)
# #Start hyper-parameter search:
# #There are 300 tasks
# #Finished 300 tasks
# #SA: elapsed time 15.08 sec
# 
# bp <- bpresult[[1]]  # typically best bpresult
# bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_")))
# allresults <- bpresult[[2]]

## ----plotSAtest, eval=FALSE---------------------------------------------------
# plotSAtest(allresults)
# # "id" corresponds to hyperparameter: lmd_we_maxlr

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/plotSAtest.png")

## ----eval=FALSE---------------------------------------------------------------
# sa <- runSA(
#   ocvecs,
#   enmat   = NULL,
#   rep     = 32,
#   getall  = FALSE,
#   lambda  = bpm[[1]],
#   we      = bpm[[2]],
#   maxlr = bpm[[3]],
#   totalit = bpm[[4]]
# )
# #> Start parameter fitting:
# #> .SA: elapsed time 0.13 sec

## ----eval=FALSE---------------------------------------------------------------
# ela <- ELA(
#   sa,
#   env            = NULL,
#   threads        = 1,
#   SS.itr         = 20000,
#   FindingTip.itr = 10000,
#   reporting      = TRUE
# )
# #> Start ELA:
# #> 7 stable states were found.
# #> Checking 21 tipping points.
# #> converting...
# #> ELA: elapsed time 2.19 sec

## ----eval=FALSE---------------------------------------------------------------
# ela[[1]] # stable state IDs
# #> [1] "09x" "01t" "EWB" "EY8" "1uV" "1mL" "17s"
# ela[[2]] # stable state energies
# #> [1] -11.099418 -10.297931  -9.544050  -8.961634  -8.005040  -6.440385  -6.269140

## ----eval=FALSE---------------------------------------------------------------
# elap <- ELPruning(ela, th=0.1, type="relative")
# #> Start pruning:
# #> *...
# #> ELPruning: elapsed time 0.05 sec

## ----eval=FALSE---------------------------------------------------------------
# elap[[1]][[1]] # updated stable state IDs
# #> [1] "09x" "EWB" "EY8" "1uV"
# elap[[1]][[2]] # updated stable state energies
# #> [1] -11.099418  -9.544050  -8.961634  -8.005040

## ----eval=FALSE---------------------------------------------------------------
# elap[[2]]
# #>  ss.before.pruning ss.after.pruning
# #> 1               09x              09x
# #> 2               01t              09x
# #> 3               EWB              EWB
# #> 4               EY8              EY8
# #> 5               1uV              1uV
# #> 6               1mL              09x
# #> 7               17s              09x

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/deltaE_scheme.png")

## ----eval=FALSE---------------------------------------------------------------
# elap <- ELPruning(ela, th=0.2, type="relative")
# #> Start pruning:
# #> *...
# #> ELPruning: elapsed time 0.05 sec
# elap[[1]][[1]] # updated stable state IDs
# #> [1] "09x" "EWB" "1uV"

## ----showDG, eval=FALSE-------------------------------------------------------
# # Before pruning
# showDG(ela, ocvecs, "example")

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/DGexample.png")

## ----showDG_pruned, eval=FALSE------------------------------------------------
# # After pruning
# showDG(elap[[1]], ocvecs, "example")

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/DGexample_pruned.png")

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/deltaE_hist.png")

## ----showDG_pvalue, eval=FALSE------------------------------------------------
# bs_params <- bootstrap_SA(ocvecs, rep = 16, threads=8,　bootstrap.it = 1000,
#                           lambda  = bpm[[1]],we = bpm[[2]],
#                           maxlr = bpm[[3]], totalit = bpm[[4]])
# 
# bootstrap <- bootstrap_ELA(bs_params,ocvecs,ela=ela,threads = 1)
# ela_rep <- bootstrap[[1]]
# bootstrap_ene <- bootstrap[[2]]
# showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene, min_pval = 0.01,
#        scale_labels = c(0.05,0.1,0.2,0.5))

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/DGexample_pval.png")

## ----showDG_IQR_pruned, eval=FALSE--------------------------------------------
# elap_rep <- ELPruning(ela, type = "bootstrap",bootstrap_ene=bootstrap_ene,lower_qr = 0.25)
# showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene, min_pval = 0.01,
#        scale_labels = c(0.05,0.1,0.2,0.5))

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/DGexample_pval_pruned.png")

## ----showDG_manual, eval=FALSE------------------------------------------------
# grobj_to_plot <- showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene,getGraphObj = TRUE)
# grobj_point <- grobj_to_plot[grobj_to_plot$point_str != "",]
# g <- ggplot(
#   grobj_to_plot,
#   aes(x=.data$nodes2xposi, y=.data$energy, label=.data$point_str)) +
#   xlab("") + geom_text(hjust=0.75, vjust=2, aes(fontface=2)) + geom_path() +
#   geom_point(
#     data = grobj_point,
#     mapping = aes(
#       x = .data$nodes2xposi,
#       y = .data$energy,
#       color = .data$energy,
#       shape = .data$type,
#     ),size = 4
#   ) +
#   scale_color_viridis_c(option = "plasma",alpha=0.8,
#                         direction=-1,na.value = "grey") +
#   coord_cartesian(xlim=c(0.75,7.5), ylim=c(-5.8,-11.5))
# plot(g)

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/DG_manual.png")

## ----ShowIntrGraph, eval=FALSE------------------------------------------------
# showIntrGraph(
#   elap[[1]],
#   sa,
#   th = 0.01,
#   annot_adj = c(0.75, 1.00)
# )

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/IntrGraph.png")

## ----eval=FALSE---------------------------------------------------------------
# bpresult <- Findbp(ocvecs, enmat = envecs, rep = 2, threads = 2,
#                    cvmode="10fold",nrepeat = 10, we=c(0.001),totalit=4000,
#                    lmd=c(1e-05,1e-03,0.01),maxlrs = c(0.005),
#                    tol=0.03, runadamW=TRUE,sparse=TRUE,fastfitting=TRUE)
# bp <- bpresult[[1]]
# bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_")))
# 
# sa <- runSA(
#   ocvecs,
#   enmat   = envecs,
#   rep     = 16,
#   getall  = FALSE,
#   lambda = bpm[[1]],
#   we = bpm[[2]],
#   maxlr = bpm[[3]],
#   totalit = bpm[[4]]
# )
# 

## ----eval=FALSE---------------------------------------------------------------
# # This process takes few minutes with a small number of threads, and here pre-computed results are shown.
# gela <- GradELA(
#   sa      = sa,
#   eid     = "factor.1",
#   enmat   = envecs,
#   env     = NULL,
#   range   = c(-1, 1),
#   steps   = 32,
#   th = 0.05,
#   threads = 8,
#   pruning_type = "relative",
#   bs_params = bs_params,
#   ocmat = ocvecs)

## ----eval=FALSE---------------------------------------------------------------
# # This process takes few minutes with a small number of threads, and here pre-computed results are shown.
# bs_params <- bootstrap_SA(ocvecs, enmat = envecs,rep = 16, threads=8,bootstrap.it = 1000,
#                           lambda  = bpm[[1]],we = bpm[[2]],
#                           maxlr = bpm[[3]], totalit = bpm[[4]])
# 
# gela <- GradELA(
#   sa      = sa,
#   eid     = "factor.1",
#   enmat   = envecs,
#   env     = NULL,
#   range   = c(-1, 1),
#   steps   = 32,
#   lower_qr = 0.25,
#   threads = 8,
#   pruning_type = "bootstrap",
#   bs_params = bs_params,
#   ocmat = ocvecs)

## ----eval=FALSE---------------------------------------------------------------
# gela[[1]][[1]][[1]][[1]] # stable state IDs at 1st environmental parameter
# # [1] "0Ht" "EWB" "1uV"
# gela[[1]][[20]][[1]][[1]] # stable state IDs at 20th environmental parameter
# # [1] "09x" "EWB"

## ----showSSD, eval=FALSE, fig.width=8, fig.height=6, message=FALSE------------
# # GradELA result (ELPruning with type="relative").
# showSSD(gela)

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/SSD.png")

## ----showSSD_bootstrap, eval=FALSE, fig.width=8, fig.height=6, message=FALSE----
# # GradELA result (ELPruning with type="bootstrap").
# showSSD(gela)

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/SSD_pruned.png")

## ----showSSD_ggplot, eval=FALSE, fig.width=8, fig.height=6, message=FALSE-----
# ssd_df <- showSSD_ggplot(gela,plot = FALSE,getSSDobj = TRUE)
# g <- ggplot(ssd_df, aes(x = env, y = Energy, color = SS_ID)) +
#             geom_point(shape = 19) +
#             geom_line(aes(group = SS_ID))
# 

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/SSD_ggplot.png")

## ----GELA3D-demo, eval=FALSE, fig.width=8, fig.height=6, message=FALSE--------
# # This process takes few minutes, and here pre-computed result is drawn.
# gelaobj <- GELAObj(gela, sa=sa,threads=2)
# showGELA3D(gelaobj,theta = 30,phi = 30)

## ----echo=FALSE, out.width="95%"----------------------------------------------
knitr::include_graphics("figures/GELA3D.png")