In AnthonyEbert/StartNetwork: Calculates KL divergence between mechanistic and statistical network models

knitr::opts_chunk$set(echo = TRUE, cache = TRUE, cache.rebuild = FALSE, autodep = TRUE)

set.seed(1)

library(StartNetwork)
library(parallel)
library(ggplot2)

n = 15
replicates = 400

Triangles 1

mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.5
theta_p <- rep(seq(0.4, 0.6, by = 0.02), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

Triangles 2

mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.4
theta_p <- rep(seq(0.3, 0.5, by = 0.02), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

Triangles 3

mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.3
theta_p <- rep(seq(0.04, 0.4, by = 0.02), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

Triangles 4

n <- 40
replicates <- 4000
mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.3
theta_p <- rep(seq(0.05, 0.35, by = 0.05), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

Triangles 5

n <- 60
replicates <- 20000
mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.3
theta_p <- rep(seq(0.05, 0.35, by = 0.05), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

Triangles 6

n = 15
replicates = 400

mech_net_triangles <- purrr::partial(mech_net_triangles_n, n = !!n)

true_value <- 0.2
theta_p <- rep(seq(0.04, 0.4, by = 0.02), 5)
theta_s <- log(true_value/(1 - true_value))

cl <- parallel::makeCluster(parallel::detectCores())

g <- parallel::parLapply(cl, theta_p, StartNetwork::KL_ss, theta_s = theta_s, replicates = replicates, sorted = TRUE, mech_net = mech_net_triangles, lstat = function(x){length(igraph::triangles(x))/3}, mirror = TRUE, type = c("Bianconi", "Liebenau"))

parallel::stopCluster(cl)

g_tidy <- StartNetwork::tidy_g(g)

ggplot(g_tidy) + aes(x = parameter, y = value, group = interaction(parameter, type), col = type) + geom_boxplot(outlier.shape = NULL, position = "identity") + facet_wrap( ~ key, scales = "free_y") + 
  ggplot2::geom_vline(mapping = aes(xintercept = true_value))

sessionInfo()

AnthonyEbert/StartNetwork documentation built on April 24, 2020, 3:28 a.m.