In ha0ye/portalDS: Compute Dynamic Stability for Rodents at Portal

LOCAL <- identical(Sys.getenv("LOCAL"), "true")
knitr::opts_chunk$set(
  purl = LOCAL,
  collapse = TRUE,
  comment = "#>",
  fig.width = 6
)

Introduction

This document demonstrates the dynamic stability analysis, using the Maizuru Bay fish community as an example dataset [@Ushio_2018].

First, we load the packages that we are going to use, and set the seed (for the random number generation associated with creating surrogate data):

library(portalDS)
library(dplyr)
library(ggplot2)

set.seed(42)

The Data

First, we load in the data and take a look at the time series for the fish populations.

data(maizuru_block)

maizuru_palette <- viridis::viridis(NCOL(maizuru_block) - 1)
names(maizuru_palette) <- sort(names(maizuru_block[, -1]))

x_axis_time <- scale_x_date(
  limits = as.Date(c("2005-03-01", "2011-01-01")),
  date_breaks = "1 year", date_labels = "%Y", expand = c(0.01, 0)
)


plot_time_series(maizuru_block) +
  x_axis_time +
  scale_color_manual(values = maizuru_palette) +
  guides(color = FALSE)

Dynamic Stability Calculations

Perform the dynamic stability analysis for the maizuru system. The steps involved in the full analysis are:

1a. Run simplex projection on each time series to identify the optimal embedding dimension 1b. Generate surrogate time series (here, using the twin-surrogate method)

simplex_results <- compute_simplex(maizuru_block,
  E_list = 1:24,
  surrogate_method = "twin",
  num_surr = 200,
  surr_params = list(
    T_period = 24,
    quantile_vec = c(
      0.875, 0.88, 0.89, 0.90, 0.91, 0.92,
      0.93, 0.94, 0.95, 0.85, 0.84, 0.83,
      0.82, 0.81, 0.80, 0.96
    )
  ), 
  id_var = "censusdate"
)

1. Run ccm on each pairwise interaction, including the surrogate data, and using the embedding dimension identified previously. Since this takes a while to run, we've commented out the actual code, and instead load a pre-generated results file.

# future::plan(future.callr::callr)
# ccm_results <- compute_ccm(simplex_results,
#                            lib_sizes = c(6, 12, 24, 40, 80, 140, 220, 320, NROW(maizuru_block)),
#                            num_samples = 200)

data_path <- system.file("extdata", "maizuru_ccm_results.RDS",
  package = "portalDS", mustWork = TRUE
)
ccm_results <- readRDS(data_path)

1. identify the significant interactions by comparing the CCM results for the real time series against the comparable calculations for the surrogate data

ccm_links <- compute_ccm_links(ccm_results)

1. run S-map models for each time series, using the appropriate number of lags, and including the important interacting variables (we rescale the time series to each have mean = 0, variance = 1)

smap_coeffs <- compute_smap_coeffs(maizuru_block, 
                                   ccm_links,
                                   rescale = TRUE, 
                                   id_var = "censusdate")

1. extract out the s-map coefficients from the models and assemble matrices for the system

smap_matrices <- compute_smap_matrices(smap_coeffs,
                                       ccm_links)

1. perform eigen-decomposition on the s-map coefficient matrices

eigen_decomp <- compute_eigen_decomp(smap_matrices)
eigenvalues <- eigen_decomp$values
eigenvectors <- eigen_decomp$vectors

Figures

Interaction Network

First, we look at the interaction network. plot_network has a default layout and palette, but we can specify both. Here, we re-use the viridis color scheme used to plot the time series.

maizuru_network <- plot_network(ccm_links, palette = maizuru_palette)

Note that the output of plot_network is a list. The plot element is the ggraph object, which we can plot:

print(maizuru_network$plot)

The remaining components allow us to plot a subset of the network, while keeping the same layout of nodes and the associated colors.

halichoeres_links <- ccm_links %>%
  filter(stringr::str_detect(lib_column, "Halichoeres") |
    stringr::str_detect(target_column, "Halichoeres"))

halichoeres_network <- plot_network(halichoeres_links,
  palette = maizuru_palette,
  existing_graph = maizuru_network$graph
)

print(halichoeres_network$plot)

Eigenvalues

Dynamic stability [@Ushio_2018] is defined as the time-varying property that is the largest (absolute value) eigenvalue. With the results of our analysis, we can plot this as a time series, applying the same x-axis transformation as for the population time series earlier:

eigenvalues %>%
  plot_eigenvalues() +
  x_axis_time

Note that plot_eigenvalues has a number of additional parameters, which allow us to plot multiple eigenvalues simultaneously, or identify when the dominant eigenvalue is complex (red points):

eigenvalues %>%
  plot_eigenvalues(num_values = 2, highlight_complex = TRUE) +
  x_axis_time

Eigenvectors

Just as with the eigenvalues, we can plot the dominant eigenvector(s) as it changes in time. Note that there are only 13 dimensions here, because we drop the 2 species that share no links with the others.

eigenvectors %>%
  plot_eigenvectors(num_values = 2) +
  scale_color_manual(values = maizuru_palette) +
  x_axis_time

References

ha0ye/portalDS documentation built on March 28, 2020, 1:21 p.m.