## ---- echo=FALSE-------------------------------------------------------------- knitr::opts_chunk$set( cache = FALSE, fig.width = 9, message = FALSE, eval=FALSE, warning = FALSE) ## ----deps, eval=TRUE---------------------------------------------------------- library(ggplot2) # library(igraph) library(colourvalues) library(GGally) library(network) library(sna) library(ape) library(dplyr) library(philentropy) library(cluster) library(miaSim) ## ----rpower------------------------------------------------------------------- # # FIXME: probably R already has some function for this; use that to simplify the code and # # avoid unnecessary duplication # rpower <- function(n, alpha, norm = FALSE) { # u <- runif(n, min = 0, max = 1 - .Machine$double.eps^0.5) # power <- (1 - u)^(1 / (1 - alpha)) # if (norm) { # return(power / mean(power)) # } else { # return(power) # } # } ## ----gibs--------------------------------------------------------------------- # params <- data.frame( # alpha = c(7, 3, 2, 1.6, 1.01), # t_end = c(20, 10, 5, 2, 1), # t_step = c(0.02, 0.01, 0.005, 0.002, 0.001) # ) # set.seed(42) # n_species <- 100 # # # Define local communities; # # we use here a smaller number in order to speed up demo calculations # # local_communities <- seq_len(100) # Too slow! # local_communities <- c(1, 5, 10, 20) # Faster ## ----store-------------------------------------------------------------------- # H <- list() # df <- list() # plot_hist <- list() # N <- list() # interactions_custom <- list() # A <- list() # g <- list() # g_plot <- list() # local_species_pool <- list() # local_A <- list() ## ----misc, warning=FALSE, eval=FALSE------------------------------------------ # set.seed(42) # for (row in seq_len(nrow(params))) { # # for each power-law distribution # print(paste("row =", row)) # # H[[row]] <- rpower(n = n_species, alpha = params$alpha[row], norm = TRUE) # df[[row]] <- data.frame(H[[row]]) # plot_hist[[row]] <- ggplot(df[[row]], aes(H[[row]])) + # ggtitle(paste("alpha =", params$alpha[[row]])) + # geom_histogram(fill = "steelblue", color = "grey20", bins = 10) + # scale_y_log10(breaks = c(1, 10, 100), limits = c(1, 100)) + # theme_bw() + # theme(plot.title = element_text(hjust = 0.5)) # # # plot_hist[[row]] # print(plot_hist[[row]]) # # N[[row]] <- matrix(rnorm(n = n_species^2, mean = 0, sd = 1), nrow = n_species) # diag(N[[row]]) <- 0 # interactions_custom[[row]] <- N[[row]] %*% diag(H[[row]]) # A[[row]] <- randomA( # n_species = 100, # diagonal = -1, # scale_off_diagonal = 0.07, # connectance = 1, # interactions = interactions_custom[[row]] # ) # g[[row]] <- graph_from_adjacency_matrix(A[[row]], weighted = TRUE, diag = FALSE) # print(paste("max edge weight =", max(E(g[[row]])$weight))) # E(g[[row]])$color <- head(colour_values(c(abs(E(g[[row]])$weight), 0), alpha = 0, include_alpha = TRUE), -1) # # g_plot[[row]] <- ggnet2( # net = g[[row]], # mode = "circle", # node.color = "black", # node.size = 1, # edge.color = "color", edge.alpha = 0.25 # ) + # labs(title=paste("alpha =", params$alpha[[row]])) + # theme(plot.title = element_text(hjust = 0.5)) # # print(g_plot[[row]]) # # local_species_pool[[row]] <- list() # local_A[[row]] <- list() # # for (local_community in local_communities) { # local_species_pool[[row]][[local_community]] <- sample(x = 100, size = 80) # local_A[[row]][[local_community]] <- A[[row]][local_species_pool[[row]][[local_community]], local_species_pool[[row]][[local_community]]] # } # } ## ----sis, eval=FALSE---------------------------------------------------------- # SIS <- list() # group_sis <- list() # for (row in seq_len(nrow(params))) { # SIS[[row]] <- head(order(H[[row]], decreasing = TRUE), 1) # group_sis[[row]] <- list() # for (local_community in local_communities) { # if (paste0("sp", SIS[[row]]) %in% rownames(local_A[[row]][[local_community]])) { # group_sis[[row]][[local_community]] <- "with SIS" # } else { # group_sis[[row]][[local_community]] <- "without SIS" # } # } # group_sis[[row]] <- unlist(group_sis[[row]]) # } ## ----glv, eval=FALSE---------------------------------------------------------- # simulation_GLV <- list() # otu_table <- list() # jsd <- list() # PCoA <- list() # PCoA_coord <- list() # bestk <- list() # PAM <- list() # SI <- list() # groups <- list() # pcoa_plot <- list() # dfs_to_bind <- list() # # scenarios <- c("original", "various growth rates", "SIS grows faster", "non-SIS grows faster") # growthRates <- vector(mode = "list", length = 4) # growthRates[[1]] <- rep(1, 80) # growthRates[[2]] <- NULL # # growthRates for scenario 3 and 4 are assigned in the following loop # customPalette <- c("#d95319", "#0072bd", "#edb120", "#7e2f8e") ## ----scenario, eval=FALSE----------------------------------------------------- # scenario <- 1 # print(paste0("scenario: ", scenarios[scenario])) # # ### make new lists #### # simulation_GLV[[scenario]] <- list() # otu_table[[scenario]] <- list() # jsd[[scenario]] <- list() # PCoA[[scenario]] <- list() # PCoA_coord[[scenario]] <- list() # bestk[[scenario]] <- list() # PAM[[scenario]] <- list() # SI[[scenario]] <- list() # groups[[scenario]] <- list() # pcoa_plot[[scenario]] <- list() # # ### run simulations #### # for (row in seq_len(nrow(params))) { # simulation_GLV[[scenario]][[row]] <- list() # otu_table[[scenario]][[row]] <- data.frame(matrix(nrow = 0, ncol = 100)) # colnames(otu_table[[scenario]][[row]]) <- paste0("sp", 1:100) # dfs_to_bind[[scenario]] <- list() # # #### overwrite growth rates #### # # for (local_community in local_communities) { # for (local_community in local_communities) { # if (scenario == 3) { # growthRates[[scenario]] <- 9.9 * local_species_pool[[row]][[local_community]] == SIS[[row]] + 0.1 # } # if (scenario == 4) { # growthRates[[scenario]] <- -9.9 * local_species_pool[[row]][[local_community]] == SIS[[row]] + 10 # } # # simulation_GLV[[scenario]][[row]][[local_community]] <- simulateGLV( # n_species = 80, # names_species = paste0("sp", local_species_pool[[row]][[local_community]]), # A = local_A[[row]][[local_community]], # x0 = rep(1, 80), # growth_rates = growthRates[[scenario]], # t_end = params$t_end[row], # t_step = params$t_step[row], # stochastic = FALSE, # norm = TRUE # ) # # # makePlot(simulation_GLV[[scenario]][[row]][[local_community]]$matrix) # dfs_to_bind[[scenario]] <- append(dfs_to_bind[[scenario]], list(simulation_GLV[[scenario]][[row]][[local_community]]$matrix[time = 1000, ])) # } # # ### form otu tables #### # otu_table[[scenario]][[row]] <- bind_rows(dfs_to_bind[[scenario]]) # otu_table[[scenario]][[row]] <- subset(otu_table[[scenario]][[row]], select = -time) # otu_table[[scenario]][[row]][is.na(otu_table[[scenario]][[row]])] <- 0 # # ### calculate jensen-shannon distance and plot PCoA #### # jsd[[scenario]][[row]] <- JSD(as.matrix(otu_table[[scenario]][[row]])) # PCoA[[scenario]][[row]] <- ape::pcoa(as.dist(jsd[[scenario]][[row]])) # PCoA_coord[[scenario]][[row]] <- PCoA[[scenario]][[row]]$vectors[, 1:2] # colnames(PCoA_coord[[scenario]][[row]]) <- c("PCo1", "PCo2") # # bestk[[scenario]][[row]] <- 2 # PAM[[scenario]][[row]] <- pam(PCoA[[scenario]][[row]]$vectors, k = 2) # SI[[scenario]][[row]] <- PAM[[scenario]][[row]]$silinfo$avg.width # # # Too slow, let us limit k to 2..3 in the example (earlier 3:10) # for (k in 2:3) { # PAMtemp <- pam(PCoA[[scenario]][[row]]$vectors, k = k) # if (PAMtemp$silinfo$avg.width > SI[[scenario]][[row]]) { # SI[[scenario]][[row]] <- PAMtemp$silinfo$avg.width # bestk[[scenario]][[row]] <- k # PAM[[scenario]][[row]] <- PAMtemp # } # } # # groups[[scenario]][[row]] <- factor(PAM[[scenario]][[row]]$clustering) # # dataframe_pcoa <- as.data.frame(PCoA_coord[[scenario]][[row]]) # dataframe_pcoa$group <- groups[[scenario]][[row]] # p <- ggplot( # dataframe_pcoa, # aes(x = PCo1, y = PCo2, color = group) # ) + # labs(title = paste("scenario", scenario, scenarios[scenario], ";", "alpha =", params$alpha[[row]])) + # geom_point() + # theme_bw() + # scale_color_manual(values = customPalette) + # theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) # # # Store the plot # pcoa_plot[[scenario]][[row]] <- p # # Show the plot # print(pcoa_plot[[scenario]][[row]]) # } # ## ----pcoa, eval=FALSE--------------------------------------------------------- # pcoa_plot_sis <- list() # for (row in seq_len(nrow(params))) { # pcoa_plot_sis[[row]] <- ggplot(as.data.frame(PCoA_coord[[1]][[row]]), aes(x = PCo1, y = PCo2)) + # geom_point(aes(colour = group_sis[[row]])) + # labs(title=paste("alpha =", params$alpha[[row]])) + # theme_bw() + # theme(plot.title = element_text(hjust = 0.5)) # print(pcoa_plot_sis[[row]]) # pcoa_validation_ plots # }