-
Notifications
You must be signed in to change notification settings - Fork 2
Hypothesis testing
Shai Pilosof edited this page Jul 21, 2020
·
1 revision
The function run_infomap_monolayer
can use shuffling algorithms built into vegan
. To use this, we need to set signif=T
and provide a shuffling method to shuff_method
. The shuffling methods are the ones detailed in ?vegan::commsim
.
network_object <- create_monolayer_object(memmott1999, bipartite = T, directed = F, group_names = c('A','P'))
# Run with shuffling
infomap_object <- run_infomap_monolayer(network_object, infomap_executable='Infomap',
flow_model = 'undirected',
silent=T, trials=20, two_level=T, seed=123,
signif = T, shuff_method = 'r00', nsim = 50)
# Plot histograms
plots <- plot_signif(infomap_object, plotit = T)
plot_grid(
plots$L_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20)),
plots$m_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20))
)
Another way is to provide shuff_method
with a list containing already shuffled networks. You can for example use the function shuffle_infomap
to produce shuffled networks with vegan
. This is necessary on the shuffling algorithm from ?vegan::commsim
needs additional arguments such as burning.
# Or can shuffle like this, if additional arguments are needed for the shuffling algorithm
shuffled <- shuffle_infomap(network_object, shuff_method='curveball', nsim=50, burnin=2000)
infomap_object <- run_infomap_monolayer(network_object, infomap_executable='Infomap',
flow_model = 'undirected',
silent=T, trials=20, two_level=T, seed=123,
signif = T, shuff_method = shuffled, nsim = 50)
plots <- plot_signif(infomap_object, plotit = F)
plot_grid(
plots$L_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20)),
plots$m_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20))
)
You can also create your own shuffling algorithm and a list of shuffled link lists. As in this example.
tur2016_altitude <- tur2016 %>% filter(altitude==1800) %>% select(from=donor,to=receptor,weight=no.grains)
tur_network <- create_monolayer_object(tur2016_altitude, directed = T, bipartite = F)
# A dedicated function to shuffle the networks, considering the flow.
shuffle_tur_networks <- function(x){ # x is a network object
m <- x$mat
# Assign off-diagona values
off_diagonal_lower <- m[lower.tri(m, diag = FALSE)]
off_diagonal_upper <- m[upper.tri(m, diag = FALSE)]
out <- matrix(0, nrow(m), ncol(m), dimnames = list(rownames(m), colnames(m)))
out[lower.tri(out, diag = FALSE)] <- sample(off_diagonal_lower, size = length(off_diagonal_lower), replace = F)
out[upper.tri(out, diag = FALSE)] <- sample(off_diagonal_upper, size = length(off_diagonal_upper), replace = F)
# Re-assign the diagonal (left intact)
diag(out) <- diag(m)
# Sanity checks
t1 <- setequal(out[upper.tri(out, diag = FALSE)], m[upper.tri(m, diag = FALSE)]) #Should be TRUE
t2 <- setequal(out[lower.tri(out, diag = FALSE)], m[lower.tri(m, diag = FALSE)]) #Should be TRUE
t3 <- identical(out[upper.tri(out, diag = FALSE)], m[upper.tri(m, diag = FALSE)]) #Should be FALSE
t4 <- identical(out[lower.tri(out, diag = FALSE)], m[lower.tri(m, diag = FALSE)]) #Should be FALSE
t5 <- all(table(m)==table(out))# Should be TRUE because it includes all the values, including diagonal
if (any(c(t1,t2,!t3,!t4,t5)==F)){stop('One or more sanity checks failed')}
return(out)
}
# Create a list with shuffled link lists
nsim <- 500
shuffled <- NULL
for (i in 1:nsim){
print(i)
x <- shuffle_tur_networks(tur_network) #Shuffle the network
x <- create_monolayer_object(x,directed = T,bipartite = F) # Create a monolayer object
shuffled[[i]] <- create_infomap_linklist(x)$edge_list_infomap #Create a link-list
}
# Use the shuffled link lists.
tur_signif <- run_infomap_monolayer(tur_network, infomap_executable='Infomap',
flow_model = 'directed',
silent=T,
trials=100, two_level=T, seed=200952, ...='-k --markov-time 50',
signif = T, shuff_method = shuffled)
print(tur_signif$pvalue)
plots <- plot_signif(tur_signif, plotit = F)
pdf('/Users/shai/Dropbox (BGU)/Apps/Overleaf/A dynamical perspective on community detection in ecological networks/figures/null_model_example.pdf',12,8)
plot_grid(
plots$L_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20)),
plots$m_plot+
theme_bw()+
theme(legend.position='none',
axis.text = element_text(size=20),
axis.title = element_text(size=20))
)
dev.off()