Fixed graphs and degreeDistributions

.Rhistory

@@ -1,58 +1,3 @@
-library(ggplot2)
-library(GeoStratR)
-data("Bios")
-plot(Bios)
-library(raster)
-library(ggplot2)
-library(GeoStratR)
-data("Bios")
-plot(Bios)
-breaks <- 3
-Temp <- Bios[[1]]
-Prec <- Bios[[4]]
-data("Bios")
-breaks <- 3
-Temp <- Bios[[1]]
-Prec <- Bios[[4]]
-cmat <- rg_biv_cmat(breaks, style = 1)
-legend <- rg_biv_get_legend(cmat, xlab = 'Temp', ylab = 'Prec')
-xy <- rg_biv_create_raster(Temp, Prec, breaks)
-map <- rg_biv_plot_raster(xy, cmat, border = st_border_proj, xlab = 'Gain', ylab = 'Loss', limits = c(10.7, 12, 46.2, 46.8))
-xy
-plot(xy)
-rg_biv_plot_raster(xy, cmat)
-rg_biv_plot_raster
-rg_biv_plot_raster
-rg_biv_plot_raster
-bivraster = xy
-crs = "+init=epsg:4326"
-bivraster %>% raster::projectRaster(crs = crs) %>%
-raster::as.data.frame(xy = TRUE)
-bivraster %>% raster::projectRaster(crs = crs)
-bivraster
-crs = "+init=epsg:4326"
-bivraster %>% raster::projectRaster(crs = crs)
-bivraster %>% #raster::projectRaster(crs = crs) %>%
-raster::as.data.frame(xy = TRUE) %>% tibble::as_tibble() %>%
-dplyr::rename(BivValue = 3) %>% tidyr::pivot_longer(names_to = "Variable",
-values_to = "bivVal", cols = BivValue)
-r_df <- bivraster %>% #raster::projectRaster(crs = crs) %>%
-raster::as.data.frame(xy = TRUE) %>% tibble::as_tibble() %>%
-dplyr::rename(BivValue = 3) %>% tidyr::pivot_longer(names_to = "Variable",
-values_to = "bivVal", cols = BivValue)
-crs
-bivraster %>% raster::projectRaster(crs = crs)
-?raster::projectRaster
-projectRaster(Bios[[1]], crs = crs)
-projectRaster(Bios[[1]], crs = "+init=epsg:4326")
-devtools::check(remote = T)
-devtools::check(remote = T)
-library(NetworkExtinction)
-library(NetworkExtinction)
-install.packages("pkgdown")
-install.packages("ragg")
-install.packages("ragg")
-library(NetworkExtinction)
 library(NetworkExtinction)
 install.packages("cranlogs")
 ?cranlogs::cran_downloads()
@@ -510,3 +455,58 @@ A <- CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)
 class(History)
 class(History$sims)
 class(NullHyp$sims)
+library(NetworkExtinction)
+NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)
+NullHyp$sims
+library(ggplot2)
+library(scales) # for muted
+NullHyp$sims |> colnames()
+library(ggplot2)
+library(scales) # for muted
+ggplot(NullHyp$sims, aes(x = NumExt, y = AccSecExt_mean)) +
+geom_ribbon(aes(ymax = Upper, ymin = 0, group = Lower))
+NullHyp$sims
+library(ggplot2)
+library(scales) # for muted
+ggplot(NullHyp$sims, aes(x = NumExt, y = AccSecExt_mean)) +
+geom_ribbon(aes(ymax = Upper, ymin = 0, group = Lower)) +
+geom_col(aes(nsim))
+nsim
+nsim
+NullHyp$sims$nsim
+max(NullHyp$sims$nsim)
+max(NullHyp$sims$nsim)/2
+NullHyp$sims
+library(ggplot2)
+library(scales) # for muted
+ggplot(NullHyp$sims, aes(x = NumExt, y = AccSecExt_mean)) +
+geom_ribbon(aes(ymax = Upper, ymin = 0, group = Lower)) +
+scale_fill_gradient2(position="bottom" , low = "blue", mid = muted("blue"), high = "red",
+midpoint = median(max(NullHyp$sims$nsim)/2))
+library(NetworkExtinction)
+library(NetworkExtinction)
+## Not run:
+data("Less_Connected")
+History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected")
+NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)
+CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)
+library(NetworkExtinction)
+## Not run:
+data("Less_Connected")
+History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected")
+NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)
+CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)
+CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)
+library(NetworkExtinction)
+## Not run:
+data("Less_Connected")
+History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected")
+NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)
+CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)
+library(NetworkExtinction)
+library(NetworkExtinction)
+## Not run:
+data("Less_Connected")
+History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected")
+NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)
+CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)

.Rproj.user/shared/notebooks/paths

@@ -1 +1,10 @@
-C:/Users/erike/Documents/[Research] Active Projects/[R Package] NetworkExtinction/NetworkExtinction_Github/README.Rmd="59A99C17"
+/home/au687614/Documents/NetworkExtinction/.github/workflows/pkgdown.yaml="49E1A421"
+/home/au687614/Documents/NetworkExtinction/README.Rmd="A12DEE6D"
+/home/au687614/Documents/NetworkExtinction/tests/testthat.R="FD93C6CD"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-CompareExtinctions.R="ECCC3341"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-DegreeDistribution.R="B9B6621D"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-ExtinctionOrder.R="37B1C277"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-ExtinctionPlot.R="497ADBD2"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-RandomExtinctions.R="9AD4FAF8"
+/home/au687614/Documents/NetworkExtinction/tests/testthat/test-SimulateExtinctions.R="AB0B02E9"
+/home/au687614/Downloads/Degree distributonV3.R="00835093"

.github/workflows/pkgdown.yaml

@@ -8,7 +8,7 @@ name: pkgdown
 
 jobs:
   pkgdown:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-20.04
     env:
       RSPM: https://packagemanager.rstudio.com/cran/__linux__/bionic/latest
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}

NAMESPACE

@@ -65,7 +65,9 @@ importFrom(network,network.size)
 importFrom(parallel,clusterExport)
 importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
+importFrom(patchwork,plot_annotation)
 importFrom(patchwork,wrap_plots)
+importFrom(purrr,discard)
 importFrom(purrr,map)
 importFrom(purrr,reduce)
 importFrom(rlang,sym)

R/Extintions.R

@@ -588,6 +588,9 @@ ExtinctionOrder <- function(Network, Order, NetworkType = "Trophic", clust.metho
 #' @importFrom utils setTxtProgressBar
 #' @importFrom utils txtProgressBar
 #' @importFrom patchwork wrap_plots
+#' @importFrom patchwork plot_annotation
+
+
 #' @author Derek Corcoran <[email protected]>
 #' @author M. Isidora Ávila-Thieme <[email protected]>
 #' @author Erik Kusch <[email protected]>
@@ -658,7 +661,9 @@ RandomExtinctions <- function(Network, nsim = 10,
   if(plot == TRUE){
     g <- ggplot(sims, aes(x = NumExt, y = AccSecExt_mean)) + geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = scales::muted("red")) + geom_line() + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() + ggplot2::theme(axis.title.x = ggplot2::element_blank())
     h <- ggplot(sims, aes(x = NumExt, y = nsim/max(nsim))) + geom_line() + theme_bw() + labs(y = "Prop", x = "Primary extinctions")
-    I <- patchwork::wrap_plots(g, h, ncol = 1, guides = "keep", heights = c(3,1))
+    I <- patchwork::wrap_plots(g, h, ncol = 1, guides = "keep", heights = c(3,1)) + patchwork::plot_annotation(tag_levels = 'A')
+
+
     print(I)
   }
 
@@ -703,6 +708,9 @@ RandomExtinctions <- function(Network, nsim = 10,
 #' @importFrom ggplot2 element_blank
 #' @importFrom ggplot2 theme
 #' @importFrom patchwork wrap_plots
+#' @importFrom patchwork plot_annotation
+
+
 #' @importFrom scales muted
 #' @author Derek Corcoran <[email protected]>
 #' @author M. Isidora Ávila-Thieme <[email protected]>
@@ -714,7 +722,7 @@ CompareExtinctions <- function(Nullmodel, Hypothesis){
     g <- ggplot(Nullmodel$sims, aes(x = NumExt, y = AccSecExt_mean)) + geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = scales::muted("red")) + geom_line(aes(color = "blue")) + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() + ggplot2::theme(axis.title.x = ggplot2::element_blank())
     g <- g + geom_point(data = Hypothesis$sims, aes(y = AccSecExt), color = "black") + geom_line(data = Hypothesis$sims, aes(y = AccSecExt, color = "black")) + scale_color_manual(name = "Comparison",values =c("black", "blue"), label = c("Observed","Null hypothesis"))
     h <- ggplot(Nullmodel$sims, aes(x = NumExt, y = nsim/max(nsim))) + geom_line() + theme_bw() + labs(y = "Prop", x = "Primary extinctions")
-    I <- patchwork::wrap_plots(g, h, ncol = 1, guides = "keep", heights = c(3,1))
+    I <- patchwork::wrap_plots(g, h, ncol = 1, guides = "keep", heights = c(3,1)) + patchwork::plot_annotation(tag_levels = 'A')
     I
     return(I)
   }

R/Powerlaw.R

@@ -47,6 +47,7 @@
 #'@importFrom MASS fitdistr
 #'@importFrom purrr map
 #'@importFrom purrr reduce
+#'@importFrom purrr discard
 #'@importFrom stats ks.test
 #'@importFrom stats glm
 #'@importFrom stats logLik
@@ -69,64 +70,109 @@ DegreeDistribution <- function(Network, scale = "arithmetic"){
 
 
   #exponential model nls
-  exp.model <- nls(Cumulative~exp(K*lambda+ c),start= list(lambda=0.1, c = 0), data = For.Graph)
-  For.Graph$Exp <- predict(exp.model)
-  Summs.exp <- glance(exp.model)
-  Summs.exp$model <- "Exp"
-  Summs.exp$Normal.Resid <- ifelse(tidy(ks.test(augment(exp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
-  Summs.exp$family <- "Exponential"
-  Summs.exp$AICcNorm <- logLik(MASS::fitdistr(augment(exp.model)$.resid, "normal"))[1]
-  Summs.exp$AICcNorm <- (4 - 2*Summs.exp$AICcNorm) + (12/(nrow(augment(exp.model)) - 1))
-  Params.exp <- tidy(exp.model)
-  Params.exp$model <- "Exp"
-  #exponential model Log
+  tryCatch({
+    exp.model <- nls(Cumulative~exp(K*lambda+ c),start= list(lambda=0.1, c = 0), data = For.Graph)
+    For.Graph$Exp <- predict(exp.model)
+    Summs.exp <- glance(exp.model)
+    Summs.exp$model <- "Exp"
+    Summs.exp$Normal.Resid <- ifelse(tidy(ks.test(augment(exp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
+    Summs.exp$family <- "Exponential"
+    Summs.exp$AICcNorm <- logLik(MASS::fitdistr(augment(exp.model)$.resid, "normal"))[1]
+    Summs.exp$AICcNorm <- (4 - 2*Summs.exp$AICcNorm) + (12/(nrow(augment(exp.model)) - 1))
+    Params.exp <- tidy(exp.model)
+    Params.exp$model <- "Exp"
+
+  }, error = function(e) {
+    # display error message when there is an error
+    message("Could not fit Exponential model: ", e$message)
+  })
+  if(!exists("Summs.exp")) {
+    Summs.exp <- NULL
+  }
+  if(!exists("Params.exp")) {
+    Params.exp <- NULL
+  }
 
-  power <- filter(For.Graph, K != 0 & Cumulative != 0)
-  logexp.model <- glm(LogCum ~ K, data = power)
-  power$LogExp <- exp(predict(logexp.model))
-  Summs.logexp <- glance(logexp.model)
-  Summs.logexp$model <- "LogExp"
-  Summs.logexp$Normal.Resid <- ifelse(tidy(ks.test(augment(logexp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
-  Summs.logexp$family <- "Exponential"
-  Summs.logexp$AICcNorm <- logLik(MASS::fitdistr(augment(logexp.model)$.resid, "normal"))[1]
-  Summs.logexp$AICcNorm <- (4 - 2*Summs.logexp$AICcNorm) + (12/(nrow(augment(logexp.model)) - 1))
-  Params.logexp <- tidy(logexp.model)
-  Params.logexp$model <- "LogExp"
+  #exponential model Log
+  tryCatch({
+    power <- filter(For.Graph, K != 0 & Cumulative != 0)
+    logexp.model <- glm(LogCum ~ K, data = power)
+    power$LogExp <- exp(predict(logexp.model))
+    Summs.logexp <- glance(logexp.model)
+    Summs.logexp$model <- "LogExp"
+    Summs.logexp$Normal.Resid <- ifelse(tidy(ks.test(augment(logexp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
+    Summs.logexp$family <- "Exponential"
+    Summs.logexp$AICcNorm <- logLik(MASS::fitdistr(augment(logexp.model)$.resid, "normal"))[1]
+    Summs.logexp$AICcNorm <- (4 - 2*Summs.logexp$AICcNorm) + (12/(nrow(augment(logexp.model)) - 1))
+    Params.logexp <- tidy(logexp.model)
+    Params.logexp$model <- "LogExp"}, error = function(e) {
+      # display error message when there is an error
+      message("Could not fit Logexponential model: ", e$message)
+    })
+  if(!exists("Summs.logexp")) {
+    Summs.logexp <- NULL
+  }
+  if(!exists("Params.logexp")) {
+    Params.logexp <- NULL
+  }
 
   #logpowerlaw
-
-  logpower.model <- glm(LogCum ~ I(log(K)), data = power)
-  power$LogPower <- exp(predict(logpower.model))
-  For.Graph <- full_join(For.Graph, power)
-  Summs.logpower <- glance(logpower.model)
-  Summs.logpower$model <- "LogPower"
-  Summs.logpower$Normal.Resid <- ifelse(tidy(ks.test(augment(logpower.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
-  Summs.logpower$family <- "PowerLaw"
-  Summs.logpower$AICcNorm <- logLik(MASS::fitdistr(augment(logpower.model)$.resid, "normal"))[1]
-  Summs.logpower$AICcNorm <- (4 - 2*Summs.logpower$AICcNorm) + (12/(nrow(augment(logpower.model)) - 1))
-  Params.logpower <- tidy(logpower.model)
-  Params.logpower$model <- "LogPower"
+  tryCatch({
+    logpower.model <- glm(LogCum ~ I(log(K)), data = power)
+    power$LogPower <- exp(predict(logpower.model))
+    For.Graph <- full_join(For.Graph, power)
+    Summs.logpower <- glance(logpower.model)
+    Summs.logpower$model <- "LogPower"
+    Summs.logpower$Normal.Resid <- ifelse(tidy(ks.test(augment(logpower.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
+    Summs.logpower$family <- "PowerLaw"
+    Summs.logpower$AICcNorm <- logLik(MASS::fitdistr(augment(logpower.model)$.resid, "normal"))[1]
+    Summs.logpower$AICcNorm <- (4 - 2*Summs.logpower$AICcNorm) + (12/(nrow(augment(logpower.model)) - 1))
+    Params.logpower <- tidy(logpower.model)
+    Params.logpower$model <- "LogPower"}, error = function(e) {
+      # display error message when there is an error
+      message("Could not fit Logpower model: ", e$message)
+    })
+  if(!exists("Summs.logpower")) {
+    Summs.logpower <- NULL
+  }
+  if(!exists("Params.logpower")) {
+    Params.logpower <- NULL
+  }
 
   #powerlaw
+  tryCatch({
+    power.model <- nls(Cumulative~a*K^y, start= list(y=0, a = 1), data = power)
+    power$Power <- predict(power.model)
+    For.Graph <- full_join(For.Graph, power)
+    Summs.power <- glance(power.model)
+    Summs.power$model <- "Power"
+    Summs.power$Normal.Resid <- ifelse(tidy(ks.test(augment(power.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
+    Summs.power$family <- "PowerLaw"
+    Summs.power$AICcNorm <- logLik(MASS::fitdistr(augment(power.model)$.resid, "normal"))[1]
+    Summs.power$AICcNorm <- (4 - 2*Summs.power$AICcNorm) + (12/(nrow(augment(power.model)) - 1))
+    Params.power <- tidy(power.model)
+    Params.power$model <- "Power"
+
+  }, error = function(e) {
+    # display error message when there is an error
+    message("Could not fit Power model: ", e$message)
+  })
+  if(!exists("Summs.power")) {
+    Summs.power <- NULL
+  }
+  if(!exists("Params.power")) {
+    Params.power <- NULL
+  }
 
-  powerlaw.model <- nls(Cumulative~a*K^y, start= list(y=0, a = 1), data = power)
-  power$Power <- predict(powerlaw.model)
-  For.Graph <- full_join(For.Graph, power)
-  Summs.power <- glance(powerlaw.model)
-  Summs.power$model <- "Power"
-  Summs.power$Normal.Resid  <- ifelse(tidy(ks.test(augment(powerlaw.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")
-  Summs.power$family <- "PowerLaw"
-  Summs.power$AICcNorm <- logLik(MASS::fitdistr(augment(powerlaw.model)$.resid, "normal"))[1]
-  Summs.power$AICcNorm <- (4 - 2*Summs.power$AICcNorm) + (12/(nrow(augment(powerlaw.model)) - 1))
-  Params.power <- tidy(powerlaw.model)
-  Params.power$model <- "Power"
 
   #all together
-  Summs <- full_join(Summs.exp, Summs.power)
-  Summs <- full_join(Summs, Summs.logexp)
-  Summs <- full_join(Summs, Summs.logpower) %>% select(logLik, AIC, BIC, model, Normal.Resid, family, AICcNorm)
-  Summs <- arrange(Summs, Normal.Resid, AIC) %>%  dplyr::select(logLik, AIC, BIC, model, Normal.Resid, family)
-  params <- bind_rows(Params.logpower, Params.power, Params.logexp, Params.exp) %>%
+  Summs <- list(Summs.exp, Summs.power, Summs.logexp, Summs.logpower) |>
+    purrr::discard(is.null) |>
+    purrr::reduce(full_join) |>
+    select(logLik, AIC, BIC, model, Normal.Resid, family, AICcNorm) |>      arrange(Normal.Resid, AIC) |>
+    dplyr::select(logLik, AIC, BIC, model, Normal.Resid, family)
+  params <- list(Params.logpower, Params.power, Params.logexp, Params.exp) |> purrr::discard(is.null) |>
+    purrr::reduce(bind_rows) %>%
     dplyr::filter(model %in% Summs$model) %>%
     mutate(term = case_when(term == "y" ~ "Beta",
                             term == "a" ~ "c",