In diazrenata/conditionalsads: What the Package Does (One Line, Title Case)

library(conditionalsads)

This analysis will also source functions from meteR and feasiblesads.

1. Generate toy data

Generate toy dataset that captures 4 scenarios:

S and N are the same for control and manipulated communities, and both are drawn from the statistical constraint.
S and N change with manipulation, but both are drawn from their respective constraints.
S and N remain the same, but the manipulated community is drawn from a near-uniform distribution rather than the statistical constraint.
S and N change, and the manipulated community is drawn from a near-uniform distribution.

Within each scenario there are nsamples communities in each treatment group. The control communities always have ctrl_s species and ctrl_n individuals. For scenarios 2 and 4, the treatment communities have approximately 80% of the species and individuals as the control communities.

toy_fs = make_toy_data(constraint = "FS", nsamples = 5, ctrl_s = 5, ctrl_n = 100)
str(toy_fs)
head(toy_fs)

2. Generate samples from constraints

For each sample, generate nsamples samples from the feasible set and from the METE distribution.

These will be the sample_sads for comparison.

nsamples <- 10


get_both_samples <- function(this_vect, nsamples, inpar = FALSE){
  library(conditionalsads)
  this_rad = as.integer(this_vect[3:length(this_vect)])
  this_rad = na.omit(this_rad)
  s = length(this_rad)
  n = sum(this_rad)
  this_fs <- sample_feasibleset(s = s, n = n, nsamples, inpar)
  this_mete <- sample_METE(s = s, n= n, nsamples)

  these_constraint_samples <- list(this_fs, this_mete)

  return(these_constraint_samples)
}

constraint_samples <- apply(toy_fs, MARGIN = 1, FUN = get_both_samples,
                            nsamples = nsamples, inpar = TRUE)

3. Calculate test statistics

Empirical R2 and sample R2 compared to most-likely constraint vector.
Empirical and sample evenness (Evar and Simpson's) and skew.

fs_r2 <- list()
fs_kl <- list()
mete_r2 <- list()
mete_kl <- list()

fs_evar <- list()
fs_simp <- list()
fs_skew <- list()

mete_evar <- list()
mete_simp <- list()
mete_skew <- list()

poilogs <- list()


for(i in 1:nrow(toy_fs)) {

 this_rad = as.integer(toy_fs[i, 3:ncol(toy_fs)])
  this_rad = na.omit(this_rad)
  s = length(this_rad)
  n = sum(this_rad)

  fs_constraint = get_fs_ct(s, n, nsamples = nsamples, newsamples = F, oldsamples = constraint_samples[[i]][[1]])
  mete_constraint  = get_mete_prediction(s, n)

  fs_r2[[i]] <- get_stat_list(empirical_rad = this_rad,
                              sampled_rads = constraint_samples[[i]][[1]],
                              stat = "r2", constraint = fs_constraint)
  fs_kl[[i]] <- get_stat_list(empirical_rad = this_rad,
                              sampled_rads = constraint_samples[[i]][[1]],
                              stat = "kl_div", constraint = fs_constraint)


  mete_r2[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[1]],
                                stat = "r2", constraint = mete_constraint)


  mete_kl[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[1]],
                                stat = "kl_div", constraint = mete_constraint)

  fs_evar[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[1]],
                                stat = "evar")
  fs_simp[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[1]],
                                stat = "simpson")
  fs_skew[[i]] <-  get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[1]],
                                stat = "rad_skew")

  mete_evar[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[2]],
                                stat = "evar")
  mete_simp[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[2]],
                                stat = "simpson")
  mete_skew[[i]] <- get_stat_list(empirical_rad = this_rad,
                                sampled_rads = constraint_samples[[i]][[2]],
                                stat = "rad_skew")

  poilogs[[i]] <- rad_poilog_cs(this_rad)

  rm(s)
  rm(n)
}

4. Get quantiles of empirical test statistics vs. constraints

fs_r2_quantile <- vapply(fs_r2, FUN = test_quantile, FUN.VALUE = 1)
fs_kl_quantile <- vapply(fs_kl, FUN = test_quantile, FUN.VALUE =1)
fs_evar_quantile <- vapply(fs_evar, FUN = test_quantile, FUN.VALUE = 1)
fs_simp_quantile <- vapply(fs_simp, FUN = test_quantile, FUN.VALUE = 1)
fs_skew_quantile <- vapply(fs_skew, FUN = test_quantile, FUN.VALUE = 1)


mete_r2_quantile <- vapply(mete_r2, FUN = test_quantile, FUN.VALUE = 1)
mete_kl_quantile <- vapply(mete_kl, FUN = test_quantile, FUN.VALUE =1)
mete_evar_quantile  <- vapply(mete_evar, FUN = test_quantile, FUN.VALUE = 1)
mete_simp_quantile  <- vapply(mete_simp, FUN = test_quantile, FUN.VALUE = 1)
mete_skew_quantile  <- vapply(mete_skew, FUN = test_quantile, FUN.VALUE = 1)


# Poilog pars

pull_poilog <- function(pl_list, item) {
  return(pl_list[[item]])
}

poilog_expmu <- vapply(poilogs, FUN = pull_poilog, FUN.VALUE = 1, item = 1)
poilog_sig <- vapply(poilogs, FUN = pull_poilog, FUN.VALUE = 1, item = 2)

toy_fs <- cbind(toy_fs, fs_r2_quantile,fs_kl_quantile, fs_evar_quantile, fs_simp_quantile, fs_skew_quantile, mete_r2_quantile, mete_kl_quantile, mete_evar_quantile, mete_simp_quantile, mete_skew_quantile, poilog_expmu,
                poilog_sig)

if(FALSE) {
  write.csv(toy_fs, "toy_fs_done250.csv")
}

5. Compare quantiles across experimental treatments

library(ggplot2)


for(i in 1:4) {

  this_h <- levels(toy_fs$h)[i]

  this_data <- toy_fs %>%
    dplyr::filter(h == this_h)

  r2_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = fs_r2_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$fs_r2_quantile)) +
    ggtitle(paste0(this_h, " r2")) +
    theme_bw()

  kl_plot <- ggplot(data = this_data, aes(x = trtmnt, y = fs_kl_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$fs_kl_quantile)) +
    ggtitle(paste0(this_h, " kl")) +
    theme_bw()

  evar_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = fs_evar_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$fs_evar_quantile)) +
    ggtitle(paste0(this_h, " evar")) +
    theme_bw()

  simp_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = fs_simp_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$fs_simp_quantile)) +
    ggtitle(paste0(this_h, " Simpson evenness")) +
    theme_bw()

  skew_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = fs_skew_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$fs_skew_quantile)) +
    ggtitle(paste0(this_h, " skewness")) +
    theme_bw()


print(gridExtra::grid.arrange(r2_plot, kl_plot, evar_plot, simp_plot, skew_plot,
                          nrow = 3))

}

for(i in 1:4) {

  this_h <- levels(toy_fs$h)[i]

  this_data <- toy_fs %>%
    dplyr::filter(h == this_h)

  r2_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = mete_r2_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$mete_r2_quantile)) +
    ggtitle(paste0(this_h, " r2")) +
    theme_bw()

  kl_plot <- ggplot(data = this_data, aes(x = trtmnt, y = mete_kl_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$mete_kl_quantile)) +
    ggtitle(paste0(this_h, " kl")) +
    theme_bw()

  evar_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = mete_evar_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$mete_evar_quantile)) +
    ggtitle(paste0(this_h, " evar")) +
    theme_bw()

  simp_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = mete_simp_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$mete_simp_quantile)) +
    ggtitle(paste0(this_h, " Simpson evenness")) +
    theme_bw()

  skew_plot <-  ggplot(data = this_data, aes(x = trtmnt, y = mete_skew_quantile)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$mete_skew_quantile)) +
    ggtitle(paste0(this_h, " skewness")) +
    theme_bw()

    poilog_mu_plot <- ggplot(data = this_data, aes(x = trtmnt, y = poilog_expmu)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$poilog_expmu)) +
    ggtitle(paste0(this_h, " poilog_expmu")) +
    theme_bw()

   poilog_sig_plot <- ggplot(data = this_data, aes(x = trtmnt, y = poilog_sig)) +
    geom_jitter(aes(x = this_data$trtmnt, y = this_data$poilog_sig)) +
    ggtitle(paste0(this_h, " poilog_sig")) +
    theme_bw()


print(gridExtra::grid.arrange(r2_plot, kl_plot, evar_plot, simp_plot, skew_plot,
                          nrow = 3))
print(gridExtra::grid.arrange(poilog_mu_plot, poilog_sig_plot,
                          nrow = 1))
}