Nothing
inst/doc/topological-analysis.R
In robber: Using Block Model to Estimate the Robustness of Ecological Network
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
warning = FALSE,
message = FALSE,
fig.width = 6,
fig.height = 4
)
## ----setup--------------------------------------------------------------------
library(robber)
library(ggplot2)
library(tidyr)
library(dplyr)
library(forcats)
library(purrr)
library(igraph)
theme_set(theme_bw(base_size = 15))
## -----------------------------------------------------------------------------
data("hostparasite", package = "robber")
data("seeddispersal", package = "robber")
data("pollination", package = "robber")
## -----------------------------------------------------------------------------
print(paste0("Host-Parasite: ", nrow(hostparasite), " ",
ncol(hostparasite), " ", mean(hostparasite)), sep = ",")
print(paste0("Seed Dispersal: ", nrow(seeddispersal), " ",
ncol(seeddispersal), " ", mean(seeddispersal)), sep = ",")
print(paste0("Pollination: ", nrow(pollination), " ",
ncol(pollination), " ", mean(pollination)), sep = ",")
## -----------------------------------------------------------------------------
(rob_emp_unif <- robustness_emp(hostparasite))
plot(rob_emp_unif)
## -----------------------------------------------------------------------------
(lbm_param <- get_lbm_param(hostparasite, ncores = 1L))
## -----------------------------------------------------------------------------
(rob_lbm_unif <- do.call(robustness_lbm, lbm_param))
## -----------------------------------------------------------------------------
rob_emp_unif$auc
rob_lbm_unif$auc
## -----------------------------------------------------------------------------
plot(rob_emp_unif) +
plot(rob_lbm_unif, add = TRUE)
## -----------------------------------------------------------------------------
rob_fun <- lapply(
X = seq_len(10),
FUN = function(i) {
A <- simulate_lbm(con = lbm_param$con,
pi = lbm_param$pi,
rho = lbm_param$rho,
nr = lbm_param$nr,
nc = lbm_param$nc)$A
return(list(unif = robustness_emp(A),
dec = robustness_emp(A, ext_seq = "decreasing"),
inc = robustness_emp(A, ext_seq = "increasing")))
}
)
## -----------------------------------------------------------------------------
pu <- plot(rob_emp_unif) +
plot(rob_lbm_unif, add = TRUE) +
plot(rob_fun[[1]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[2]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[3]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[4]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[5]]$unif, add = TRUE, lty = "dotted", col = "red") +
ggtitle("Uniform")
pu
## -----------------------------------------------------------------------------
rob_emp_dec <- robustness_emp(hostparasite, ext_seq = "decreasing")
rob_emp_inc <- robustness_emp(hostparasite, ext_seq = "increasing")
rob_emp_dec$auc
rob_emp_inc$auc
## -----------------------------------------------------------------------------
rob_emp_dec_lin <- robustness_emp(hostparasite, ext_seq = "decreasing",
method = "linear")
rob_emp_inc_lin <- robustness_emp(hostparasite, ext_seq = "increasing",
method = "linear")
rob_emp_dec_lin$auc
rob_emp_inc_lin$auc
## -----------------------------------------------------------------------------
rob_lbm_dec <- do.call(robustness_lbm, c(lbm_param, ext_seq = "decreasing"))
rob_lbm_inc <- do.call(robustness_lbm, c(lbm_param, ext_seq = "increasing"))
rob_lbm_dec$auc
rob_lbm_inc$auc
## -----------------------------------------------------------------------------
plot(rob_lbm_unif, col = "black", lty = 1) +
plot(rob_emp_unif, add = TRUE, col = "blue", lty = 2) +
plot(rob_emp_dec, add = TRUE, lty = 3, col = "blue") +
plot(rob_emp_dec_lin, add = TRUE, lty = 4, col = "blue") +
plot(rob_emp_inc, add = TRUE, lty = 3, col = "blue") +
plot(rob_emp_inc_lin, add = TRUE, lty = 4, col = "blue") +
plot(rob_lbm_dec, add = TRUE, lty = 1, col = "black") +
plot(rob_lbm_inc, add = TRUE, lty = 1, col = "black") +
ggtitle("Robustness of hostparasite")
## -----------------------------------------------------------------------------
pl_param <- get_lbm_param(pollination, ncores = 1L)
sd_param <- get_lbm_param(seeddispersal, ncores = 1L)
## -----------------------------------------------------------------------------
pl_param
## -----------------------------------------------------------------------------
sd_param
## -----------------------------------------------------------------------------
pl_lbm_unif <- do.call(robustness_lbm, pl_param)
pl_lbm_dec <- do.call(robustness_lbm, c(pl_param, ext_seq = "decreasing"))
pl_lbm_inc <- do.call(robustness_lbm, c(pl_param, ext_seq = "increasing"))
sd_lbm_unif <- do.call(robustness_lbm, sd_param)
sd_lbm_dec <- do.call(robustness_lbm, c(sd_param, ext_seq = "decreasing"))
sd_lbm_inc <- do.call(robustness_lbm, c(sd_param, ext_seq = "increasing"))
## ----echo=FALSE---------------------------------------------------------------
df <- data.frame(
Extinction = c("Uniform", "Increasing", "Decreasing"),
hostparasite = c(rob_lbm_unif$auc, rob_lbm_inc$auc, rob_lbm_dec$auc),
pollination = c(pl_lbm_unif$auc, pl_lbm_inc$auc, pl_lbm_dec$auc),
seeddispersal = c(sd_lbm_unif$auc, sd_lbm_inc$auc, sd_lbm_dec$auc))
knitr::kable(df)
## -----------------------------------------------------------------------------
comp_unif <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51)
comp_inc <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51, ext_seq = "increasing")
comp_dec <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51, ext_seq = "decreasing")
## ----echo = FALSE-------------------------------------------------------------
dfc <- data.frame(
Extinction = c("Uniform", "Increasing", "Decreasing"),
hostparasite = c(comp_unif[[1]], comp_inc[[1]], comp_dec[[1]]),
pollination = c(comp_unif[[2]], comp_inc[[2]], comp_dec[[2]]),
seeddispersal = c(comp_unif[[3]], comp_inc[[3]], comp_dec[[3]]))
knitr::kable(dfc)
## ----set_constant-------------------------------------------------------------
dens <- .0156
nr <- 100
nc <- 100
rob_er <- auc_robustness_lbm(matrix(dens, 1,1), 1, 1, nr, nc)
rob_er
## ----echo=FALSE---------------------------------------------------------------
mod <- matrix(c("a", 1, 1, "a"), 2, 2)
cp <- matrix(c("a", "a", "a", 1), 2, 2)
## ----echo=FALSE---------------------------------------------------------------
knitr::kable(mod)
## ----echo=FALSE---------------------------------------------------------------
knitr::kable(cp)
## ----parameters---------------------------------------------------------------
pi <- c(1/4, 3/4)
## ----analyse_topo, echo = FALSE-----------------------------------------------
if (! file.exists("res_topo.rds")) {
robust_topology <- tibble::tibble()
eps <- c(1/seq(8, 1.5, by = -.5), seq(1, 8, by = .5))
for(i in seq(19)) {
rho <- c(i*.05, (20-i)*.05)
list_con_mod <- lapply( eps,
function(j) {
list(con = matrix(c(j, 1,
1, j), 2, 2),
pi = pi,
rho = rho)})
rob_mod <- purrr::map_dfc(
.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
.f = function(x)
unlist(compare_robustness(list_param = list_con_mod,
dens = dens,
new_nr = nr,
new_nc = nc,
ext_seq = x))) %>%
dplyr::mutate(Topology = "Modular",
rho = rho[1],
j = seq_along(eps))
list_con_nest <- lapply( eps,
function(j) {
list(con = matrix(c(j, j,
j, 1), 2, 2),
pi = pi,
rho = rho)})
rob_nest <- purrr::map_dfc(.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
.f = function(x)
unlist(compare_robustness(list_param = list_con_nest,
dens = dens,
new_nr = nr,
new_nc = nc,
ext_seq = x))) %>%
dplyr::mutate(Topology = "Core-Periphery",
rho = rho[1],
j = seq_along(eps))
robust_topology <- dplyr::bind_rows(robust_topology, rob_mod, rob_nest)
saveRDS(robust_topology, "res_topo.rds")
}
} else {
robust_topology <- readRDS("res_topo.rds")
}
## ----eval=FALSE---------------------------------------------------------------
# robust_topology <- tibble::tibble()
# eps <- c(1/seq(8, 1.5, by = -.5), seq(1, 8, by = .5))
# for(i in seq(19)) {
# rho <- c(i*.05, (20-i)*.05)
# list_con_mod <- lapply( eps,
# function(j) {
# list(con = matrix(c(j, 1,
# 1, j), 2, 2),
# pi = pi,
# rho = rho)})
# rob_mod <- purrr::map_dfc(
# .x = purrr::set_names(c("uniform", "increasing", "decreasing")),
# .f = function(x)
# unlist(compare_robustness(list_param = list_con_mod,
# dens = dens,
# new_nr = nr,
# new_nc = nc,
# ext_seq = x))) %>%
# dplyr::mutate(Topology = "Modular",
# rho = rho[1],
# j = seq_along(eps))
# list_con_nest <- lapply( eps,
# function(j) {
# list(con = matrix(c(j, j,
# j, 1), 2, 2),
# pi = pi,
# rho = rho)})
# rob_nest <- purrr::map_dfc(.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
# .f = function(x)
# unlist(compare_robustness(list_param = list_con_nest,
# dens = dens,
# new_nr = nr,
# new_nc = nc,
# ext_seq = x))) %>%
# dplyr::mutate(Topology = "Core-Periphery",
# rho = rho[1],
# j = seq_along(eps))
# robust_topology <- dplyr::bind_rows(robust_topology, rob_mod, rob_nest)
# }
## ----print_topo---------------------------------------------------------------
prt <- robust_topology %>%
# rename(Uniform = uniform, Increasing = increasing, Decreasing = decreasing) +
tidyr::pivot_longer(cols = c("decreasing","uniform", "increasing"),
names_to = "Extinction",
values_to = "Robustness") %>%
mutate(Extinction = forcats::as_factor(
case_when(Extinction == "decreasing" ~ "Decreasing",
Extinction == "uniform" ~ "Uniform",
Extinction == "increasing" ~ "Increasing"))) %>%
ggplot(ggplot2::aes(x = j, y = rho)) +
geom_tile(ggplot2::aes(fill = Robustness)) +
facet_grid(Topology ~ Extinction) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = rob_er,
guide = guide_colorbar(title = "R")) +
annotate(x = 15, y = .5, geom = "point", col = "black", size = 3) +
annotate(x = 15, y = .55, geom = "text", label = "ER",
col = "black", size = 3) +
scale_x_continuous(breaks = c(0, 7, 15, 22, 30),
labels = c(1/8, 1/4, 1, 4, 8)) +
coord_fixed(ratio = 30) +
theme_bw(base_size = 15, base_line_size = 0)
# df_shape <- data.frame(
# x = c(1, 1, 29, 29, 1, 14),
# y = c(.05, .05, .05, .5, .95, .05),
# button = c("a", "d", "b","e", "c", "f"),
# Topology = rep(c("Core-Periphery", "Modular"), 3),
# Extinction = c(rep("Decreasing", 2), rep("Increasing", 2), rep("Uniform", 2)))
df_shape <- data.frame(
x = c(1 , 29 , 29, 29, 1, 1),
y = c(.85, .05, .35, .95, .6, .4),
button = c(15L, 17L, 18L, 0L, 2L, 5L),#,
Topology = rep(c("Core-Periphery", "Modular"), each = 3))
# Extinction = c(rep("Decreasing", 2), rep("Increasing", 2), rep("Uniform", 2)))
prt + geom_point(data = df_shape, mapping = aes(x = x, y = y, shape = button),
size = 4, show.legend = FALSE)+
scale_shape_identity()
## ----compute_gnp--------------------------------------------------------------
if (! file.exists("res_gnp.rds")) {
n_iter <- 300
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
nr <- 100
nc <- 100
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp <- bind_rows(res)
n_iter <- 300
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
nr <- 100
nc <- 100
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod <- bind_rows(res)
res_gnp <- bind_rows(res_cp, res_mod, res_er)
saveRDS(res_gnp, "res_gnp.rds")
} else {
res_gnp <- readRDS("res_gnp.rds")
}
## ----compute_gnm--------------------------------------------------------------
if (! file.exists("res_gnm.rds")) {
n_iter <- 300
nr <- 100
nc <- 100
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp_gnm <- bind_rows(res)
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er_gnm <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod_gnm <- bind_rows(res)
res_gnm <- bind_rows(res_cp_gnm, res_mod_gnm, res_er_gnm)
saveRDS(res_gnm, "res_gnm.rds")
} else {
res_gnm <- readRDS("res_gnm.rds")
}
## ----compute_fixed_block------------------------------------------------------
if (! file.exists("res_block.rds")) {
n_iter <- 300
nr <- 100
nc <- 100
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
simz <- function() {
A <- diag(0, 100)
for (k in seq_along(pi)) {
for (q in seq_along(rho)) {
A[Z == k, W == q] <- igraph::as_incidence_matrix(
igraph::sample_bipartite(sum(Z == k), sum(W == q),
type = "gnp", p = con[k,q]))
}
}
return(A)
}
Z <- rep(seq(2), times = 100*pi)
W <- rep(seq(2), times = 100*rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp_block <- bind_rows(res)
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
Z <- rep(seq(1), times = 100*pi)
W <- rep(seq(1), times = 100*rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er_block <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
W <- rep(seq(2), times = 100*rho)
Z <- rep(seq(2), times = 100*pi)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod_block <- bind_rows(res)
res_block <- bind_rows(res_cp_block, res_mod_block, res_er_block)
saveRDS(res_block, "res_block.rds")
} else {
res_block <- readRDS("res_block.rds")
}
## -----------------------------------------------------------------------------
res_gnp %>%
mutate(Model = "Canonical") %>%
pivot_longer(cols = - c(type, Model), names_to = "Method", values_to = "Robustness") %>%
# rename(Topology = type) %>%
mutate(
Topology = case_when(
type == "cp" ~ "Core-Periphery",
type == "er" ~ "ER",
type == "mod" ~ "Modular"
)) %>%
filter(Method == "mc") %>%
ggplot(aes(x = Robustness, fill = Topology, lty = Topology)) +
geom_density(alpha = .3) +
# ggridges::geom_density_ridges(aes(y = Model), alpha = .3) +
# geom_density(alpha = .2) +
ylab(label = "") +
# facet_wrap(~ Model) +
theme_minimal(base_size = 15)
Try the robber package in your browser
Any scripts or data that you put into this service are public.
robber documentation built on May 29, 2024, 5:48 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
warning = FALSE,
message = FALSE,
fig.width = 6,
fig.height = 4
)
## ----setup--------------------------------------------------------------------
library(robber)
library(ggplot2)
library(tidyr)
library(dplyr)
library(forcats)
library(purrr)
library(igraph)
theme_set(theme_bw(base_size = 15))
## -----------------------------------------------------------------------------
data("hostparasite", package = "robber")
data("seeddispersal", package = "robber")
data("pollination", package = "robber")
## -----------------------------------------------------------------------------
print(paste0("Host-Parasite: ", nrow(hostparasite), " ",
ncol(hostparasite), " ", mean(hostparasite)), sep = ",")
print(paste0("Seed Dispersal: ", nrow(seeddispersal), " ",
ncol(seeddispersal), " ", mean(seeddispersal)), sep = ",")
print(paste0("Pollination: ", nrow(pollination), " ",
ncol(pollination), " ", mean(pollination)), sep = ",")
## -----------------------------------------------------------------------------
(rob_emp_unif <- robustness_emp(hostparasite))
plot(rob_emp_unif)
## -----------------------------------------------------------------------------
(lbm_param <- get_lbm_param(hostparasite, ncores = 1L))
## -----------------------------------------------------------------------------
(rob_lbm_unif <- do.call(robustness_lbm, lbm_param))
## -----------------------------------------------------------------------------
rob_emp_unif$auc
rob_lbm_unif$auc
## -----------------------------------------------------------------------------
plot(rob_emp_unif) +
plot(rob_lbm_unif, add = TRUE)
## -----------------------------------------------------------------------------
rob_fun <- lapply(
X = seq_len(10),
FUN = function(i) {
A <- simulate_lbm(con = lbm_param$con,
pi = lbm_param$pi,
rho = lbm_param$rho,
nr = lbm_param$nr,
nc = lbm_param$nc)$A
return(list(unif = robustness_emp(A),
dec = robustness_emp(A, ext_seq = "decreasing"),
inc = robustness_emp(A, ext_seq = "increasing")))
}
)
## -----------------------------------------------------------------------------
pu <- plot(rob_emp_unif) +
plot(rob_lbm_unif, add = TRUE) +
plot(rob_fun[[1]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[2]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[3]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[4]]$unif, add = TRUE, lty = "dotted", col = "red") +
plot(rob_fun[[5]]$unif, add = TRUE, lty = "dotted", col = "red") +
ggtitle("Uniform")
pu
## -----------------------------------------------------------------------------
rob_emp_dec <- robustness_emp(hostparasite, ext_seq = "decreasing")
rob_emp_inc <- robustness_emp(hostparasite, ext_seq = "increasing")
rob_emp_dec$auc
rob_emp_inc$auc
## -----------------------------------------------------------------------------
rob_emp_dec_lin <- robustness_emp(hostparasite, ext_seq = "decreasing",
method = "linear")
rob_emp_inc_lin <- robustness_emp(hostparasite, ext_seq = "increasing",
method = "linear")
rob_emp_dec_lin$auc
rob_emp_inc_lin$auc
## -----------------------------------------------------------------------------
rob_lbm_dec <- do.call(robustness_lbm, c(lbm_param, ext_seq = "decreasing"))
rob_lbm_inc <- do.call(robustness_lbm, c(lbm_param, ext_seq = "increasing"))
rob_lbm_dec$auc
rob_lbm_inc$auc
## -----------------------------------------------------------------------------
plot(rob_lbm_unif, col = "black", lty = 1) +
plot(rob_emp_unif, add = TRUE, col = "blue", lty = 2) +
plot(rob_emp_dec, add = TRUE, lty = 3, col = "blue") +
plot(rob_emp_dec_lin, add = TRUE, lty = 4, col = "blue") +
plot(rob_emp_inc, add = TRUE, lty = 3, col = "blue") +
plot(rob_emp_inc_lin, add = TRUE, lty = 4, col = "blue") +
plot(rob_lbm_dec, add = TRUE, lty = 1, col = "black") +
plot(rob_lbm_inc, add = TRUE, lty = 1, col = "black") +
ggtitle("Robustness of hostparasite")
## -----------------------------------------------------------------------------
pl_param <- get_lbm_param(pollination, ncores = 1L)
sd_param <- get_lbm_param(seeddispersal, ncores = 1L)
## -----------------------------------------------------------------------------
pl_param
## -----------------------------------------------------------------------------
sd_param
## -----------------------------------------------------------------------------
pl_lbm_unif <- do.call(robustness_lbm, pl_param)
pl_lbm_dec <- do.call(robustness_lbm, c(pl_param, ext_seq = "decreasing"))
pl_lbm_inc <- do.call(robustness_lbm, c(pl_param, ext_seq = "increasing"))
sd_lbm_unif <- do.call(robustness_lbm, sd_param)
sd_lbm_dec <- do.call(robustness_lbm, c(sd_param, ext_seq = "decreasing"))
sd_lbm_inc <- do.call(robustness_lbm, c(sd_param, ext_seq = "increasing"))
## ----echo=FALSE---------------------------------------------------------------
df <- data.frame(
Extinction = c("Uniform", "Increasing", "Decreasing"),
hostparasite = c(rob_lbm_unif$auc, rob_lbm_inc$auc, rob_lbm_dec$auc),
pollination = c(pl_lbm_unif$auc, pl_lbm_inc$auc, pl_lbm_dec$auc),
seeddispersal = c(sd_lbm_unif$auc, sd_lbm_inc$auc, sd_lbm_dec$auc))
knitr::kable(df)
## -----------------------------------------------------------------------------
comp_unif <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51)
comp_inc <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51, ext_seq = "increasing")
comp_dec <- compare_robustness(list(lbm_param, pl_param, sd_param),
dens = .148, new_nr = 36, new_nc = 51, ext_seq = "decreasing")
## ----echo = FALSE-------------------------------------------------------------
dfc <- data.frame(
Extinction = c("Uniform", "Increasing", "Decreasing"),
hostparasite = c(comp_unif[[1]], comp_inc[[1]], comp_dec[[1]]),
pollination = c(comp_unif[[2]], comp_inc[[2]], comp_dec[[2]]),
seeddispersal = c(comp_unif[[3]], comp_inc[[3]], comp_dec[[3]]))
knitr::kable(dfc)
## ----set_constant-------------------------------------------------------------
dens <- .0156
nr <- 100
nc <- 100
rob_er <- auc_robustness_lbm(matrix(dens, 1,1), 1, 1, nr, nc)
rob_er
## ----echo=FALSE---------------------------------------------------------------
mod <- matrix(c("a", 1, 1, "a"), 2, 2)
cp <- matrix(c("a", "a", "a", 1), 2, 2)
## ----echo=FALSE---------------------------------------------------------------
knitr::kable(mod)
## ----echo=FALSE---------------------------------------------------------------
knitr::kable(cp)
## ----parameters---------------------------------------------------------------
pi <- c(1/4, 3/4)
## ----analyse_topo, echo = FALSE-----------------------------------------------
if (! file.exists("res_topo.rds")) {
robust_topology <- tibble::tibble()
eps <- c(1/seq(8, 1.5, by = -.5), seq(1, 8, by = .5))
for(i in seq(19)) {
rho <- c(i*.05, (20-i)*.05)
list_con_mod <- lapply( eps,
function(j) {
list(con = matrix(c(j, 1,
1, j), 2, 2),
pi = pi,
rho = rho)})
rob_mod <- purrr::map_dfc(
.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
.f = function(x)
unlist(compare_robustness(list_param = list_con_mod,
dens = dens,
new_nr = nr,
new_nc = nc,
ext_seq = x))) %>%
dplyr::mutate(Topology = "Modular",
rho = rho[1],
j = seq_along(eps))
list_con_nest <- lapply( eps,
function(j) {
list(con = matrix(c(j, j,
j, 1), 2, 2),
pi = pi,
rho = rho)})
rob_nest <- purrr::map_dfc(.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
.f = function(x)
unlist(compare_robustness(list_param = list_con_nest,
dens = dens,
new_nr = nr,
new_nc = nc,
ext_seq = x))) %>%
dplyr::mutate(Topology = "Core-Periphery",
rho = rho[1],
j = seq_along(eps))
robust_topology <- dplyr::bind_rows(robust_topology, rob_mod, rob_nest)
saveRDS(robust_topology, "res_topo.rds")
}
} else {
robust_topology <- readRDS("res_topo.rds")
}
## ----eval=FALSE---------------------------------------------------------------
# robust_topology <- tibble::tibble()
# eps <- c(1/seq(8, 1.5, by = -.5), seq(1, 8, by = .5))
# for(i in seq(19)) {
# rho <- c(i*.05, (20-i)*.05)
# list_con_mod <- lapply( eps,
# function(j) {
# list(con = matrix(c(j, 1,
# 1, j), 2, 2),
# pi = pi,
# rho = rho)})
# rob_mod <- purrr::map_dfc(
# .x = purrr::set_names(c("uniform", "increasing", "decreasing")),
# .f = function(x)
# unlist(compare_robustness(list_param = list_con_mod,
# dens = dens,
# new_nr = nr,
# new_nc = nc,
# ext_seq = x))) %>%
# dplyr::mutate(Topology = "Modular",
# rho = rho[1],
# j = seq_along(eps))
# list_con_nest <- lapply( eps,
# function(j) {
# list(con = matrix(c(j, j,
# j, 1), 2, 2),
# pi = pi,
# rho = rho)})
# rob_nest <- purrr::map_dfc(.x = purrr::set_names(c("uniform", "increasing", "decreasing")),
# .f = function(x)
# unlist(compare_robustness(list_param = list_con_nest,
# dens = dens,
# new_nr = nr,
# new_nc = nc,
# ext_seq = x))) %>%
# dplyr::mutate(Topology = "Core-Periphery",
# rho = rho[1],
# j = seq_along(eps))
# robust_topology <- dplyr::bind_rows(robust_topology, rob_mod, rob_nest)
# }
## ----print_topo---------------------------------------------------------------
prt <- robust_topology %>%
# rename(Uniform = uniform, Increasing = increasing, Decreasing = decreasing) +
tidyr::pivot_longer(cols = c("decreasing","uniform", "increasing"),
names_to = "Extinction",
values_to = "Robustness") %>%
mutate(Extinction = forcats::as_factor(
case_when(Extinction == "decreasing" ~ "Decreasing",
Extinction == "uniform" ~ "Uniform",
Extinction == "increasing" ~ "Increasing"))) %>%
ggplot(ggplot2::aes(x = j, y = rho)) +
geom_tile(ggplot2::aes(fill = Robustness)) +
facet_grid(Topology ~ Extinction) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = rob_er,
guide = guide_colorbar(title = "R")) +
annotate(x = 15, y = .5, geom = "point", col = "black", size = 3) +
annotate(x = 15, y = .55, geom = "text", label = "ER",
col = "black", size = 3) +
scale_x_continuous(breaks = c(0, 7, 15, 22, 30),
labels = c(1/8, 1/4, 1, 4, 8)) +
coord_fixed(ratio = 30) +
theme_bw(base_size = 15, base_line_size = 0)
# df_shape <- data.frame(
# x = c(1, 1, 29, 29, 1, 14),
# y = c(.05, .05, .05, .5, .95, .05),
# button = c("a", "d", "b","e", "c", "f"),
# Topology = rep(c("Core-Periphery", "Modular"), 3),
# Extinction = c(rep("Decreasing", 2), rep("Increasing", 2), rep("Uniform", 2)))
df_shape <- data.frame(
x = c(1 , 29 , 29, 29, 1, 1),
y = c(.85, .05, .35, .95, .6, .4),
button = c(15L, 17L, 18L, 0L, 2L, 5L),#,
Topology = rep(c("Core-Periphery", "Modular"), each = 3))
# Extinction = c(rep("Decreasing", 2), rep("Increasing", 2), rep("Uniform", 2)))
prt + geom_point(data = df_shape, mapping = aes(x = x, y = y, shape = button),
size = 4, show.legend = FALSE)+
scale_shape_identity()
## ----compute_gnp--------------------------------------------------------------
if (! file.exists("res_gnp.rds")) {
n_iter <- 300
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
nr <- 100
nc <- 100
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp <- bind_rows(res)
n_iter <- 300
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
nr <- 100
nc <- 100
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con)
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod <- bind_rows(res)
res_gnp <- bind_rows(res_cp, res_mod, res_er)
saveRDS(res_gnp, "res_gnp.rds")
} else {
res_gnp <- readRDS("res_gnp.rds")
}
## ----compute_gnm--------------------------------------------------------------
if (! file.exists("res_gnm.rds")) {
n_iter <- 300
nr <- 100
nc <- 100
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp_gnm <- bind_rows(res)
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er_gnm <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
pi <- nr * pi
rho <- nc * rho
con <- round(outer(pi, rho) * con)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simulate_lbm(nr = nr, nc = nc, pi = pi, rho = rho, con = con, method = "gnm")
auc_mc <- robustness_emp(G$A, nb_iter = n_iter)$auc
par <- get_lbm_param(G$A)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod_gnm <- bind_rows(res)
res_gnm <- bind_rows(res_cp_gnm, res_mod_gnm, res_er_gnm)
saveRDS(res_gnm, "res_gnm.rds")
} else {
res_gnm <- readRDS("res_gnm.rds")
}
## ----compute_fixed_block------------------------------------------------------
if (! file.exists("res_block.rds")) {
n_iter <- 300
nr <- 100
nc <- 100
pi <- c(.1,.9)
rho <- c(.2,.8)
con <- matrix(c(.8, .4, .4, .1), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
simz <- function() {
A <- diag(0, 100)
for (k in seq_along(pi)) {
for (q in seq_along(rho)) {
A[Z == k, W == q] <- igraph::as_incidence_matrix(
igraph::sample_bipartite(sum(Z == k), sum(W == q),
type = "gnp", p = con[k,q]))
}
}
return(A)
}
Z <- rep(seq(2), times = 100*pi)
W <- rep(seq(2), times = 100*rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "cp", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_cp_block <- bind_rows(res)
pi <- 1
rho <- 1
con <- matrix(0.1, 1, 1)
Z <- rep(seq(1), times = 100*pi)
W <- rep(seq(1), times = 100*rho)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "er", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_er_block <- bind_rows(res)
pi <- c(.2, .8)
rho <- c(.9, .1)
con <- matrix(c(.6, .1, .2, .2), 2, 2)
con <- 0.1*con/as.vector(pi%*%con%*%rho)
W <- rep(seq(2), times = 100*rho)
Z <- rep(seq(2), times = 100*pi)
res <- pbmcapply::pbmclapply(
X = seq(500),
FUN = function(i) {
G <- simz()
auc_mc <- robustness_emp(G, nb_iter = n_iter)$auc
par <- get_lbm_param(G)
auc_lbm <- robustness_lbm(par$con, par$pi, par$rho, nr, nc)$auc
return(tibble(type = "mod", mc = auc_mc, lbm = auc_lbm ))
}, mc.cores = 8
)
res_mod_block <- bind_rows(res)
res_block <- bind_rows(res_cp_block, res_mod_block, res_er_block)
saveRDS(res_block, "res_block.rds")
} else {
res_block <- readRDS("res_block.rds")
}
## -----------------------------------------------------------------------------
res_gnp %>%
mutate(Model = "Canonical") %>%
pivot_longer(cols = - c(type, Model), names_to = "Method", values_to = "Robustness") %>%
# rename(Topology = type) %>%
mutate(
Topology = case_when(
type == "cp" ~ "Core-Periphery",
type == "er" ~ "ER",
type == "mod" ~ "Modular"
)) %>%
filter(Method == "mc") %>%
ggplot(aes(x = Robustness, fill = Topology, lty = Topology)) +
geom_density(alpha = .3) +
# ggridges::geom_density_ridges(aes(y = Model), alpha = .3) +
# geom_density(alpha = .2) +
ylab(label = "") +
# facet_wrap(~ Model) +
theme_minimal(base_size = 15)
Try the robber package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.