264 lines
No EOL
10 KiB
Text
264 lines
No EOL
10 KiB
Text
<<libraries, echo = FALSE, include = FALSE>>=
|
|
require("ggplot2")
|
|
require("ggokabeito")
|
|
require("tidyr")
|
|
require("dplyr")
|
|
require("stringr")
|
|
require("knitr")
|
|
require("pander")
|
|
require("patchwork")
|
|
require("latex2exp")
|
|
@
|
|
|
|
|
|
<<setup, echo = FALSE>>=
|
|
options(dplyr.summarise.inform = FALSE)
|
|
knitr::opts_knit$set(kable.force.latex = TRUE)
|
|
|
|
meanse <- function(x, ...) {
|
|
mean1 <- signif(round(mean(x, na.rm = T), 2), 5) # calculate mean and round
|
|
se1 <- signif(round(sd(x, na.rm = T) / sqrt(sum(!is.na(x))), 2), 2) # std error - round adding zeros
|
|
out <- paste(mean1, "$\\pm$", se1) # paste together mean plus/minus and standard error
|
|
if (str_detect(out, "NA")) {
|
|
out <- "NA"
|
|
} # if missing do not add plusminus
|
|
if (se1 == 0) {
|
|
out <- paste(mean1)
|
|
}
|
|
return(out)
|
|
}
|
|
@
|
|
|
|
<<import-data, echo = FALSE>>=
|
|
filenames <- list.files(
|
|
path = "./data/",
|
|
pattern = "inference_testing_2023-07-*",
|
|
full.names = TRUE
|
|
)
|
|
data_list <- lapply(filenames, readRDS)
|
|
col_id_BICLS <- c(11, 16, 23, 30, 37)
|
|
result_data_frame <- dplyr::bind_rows(data_list)
|
|
|
|
# Compute the preferred model
|
|
result_data_frame <- cbind(result_data_frame, preferred_model = sapply(seq_len(nrow(result_data_frame)), function(n) sub("_BICL", "", names(which.max(result_data_frame[n, col_id_BICLS])))))
|
|
result_data_frame$preferred_model <- factor(result_data_frame$preferred_model, levels = c(
|
|
"sep", "iid", "pi",
|
|
"rho", "pirho"
|
|
))
|
|
@
|
|
|
|
# Efficiency of the inference
|
|
|
|
\paragraph{Simulation settings} For this simulation the data is simulated with
|
|
$M = 2, n_{1}^{m} = 120, n_{2}^{m} = 120, Q_1 = Q_2 = 4$, $\bm{\alpha}, \bm{\pi}$
|
|
and $\bm{\rho}$ are set as follows:
|
|
\begin{align*}
|
|
&&\bm{\alpha} = .25 +
|
|
\begin{pmatrix}
|
|
3 \eps[\alpha] & 2 \eps[\alpha] & \eps[\alpha] & - \eps[\alpha]\\
|
|
2 \eps[\alpha] & 2 \eps[\alpha] & - \eps[\alpha] & \eps[\alpha]\\
|
|
\eps[\alpha] & - \eps[\alpha] & \eps[\alpha] & 2 \eps[\alpha]\\
|
|
- \eps[\alpha] & \eps[\alpha] & 2 \eps[\alpha] & 0
|
|
\end{pmatrix}, \\ \bm{\pi}^1 = \sigma_1
|
|
\begin{pmatrix}
|
|
0.2 & 0.4 & 0.4 & 0
|
|
\end{pmatrix},
|
|
&& \bm{\pi}^2 =
|
|
\begin{pmatrix}
|
|
0.25 & 0.25 & 0.25 & 0.25
|
|
\end{pmatrix}, \\
|
|
\bm{\rho}^1 =
|
|
\begin{pmatrix}
|
|
0.25 & 0.25 & 0.25 & 0.25
|
|
\end{pmatrix}, &&
|
|
\bm{\rho}^2 = \sigma_2
|
|
\begin{pmatrix}
|
|
0 & 0.33 & 0.33 & 0.33
|
|
\end{pmatrix}, &&
|
|
\end{align*}
|
|
with $\eps[\alpha]$ taking nine equally spaced values ranging from 0 to 0.24.
|
|
For each value of $\eps[\alpha]$, 108 datasets ($X_1, X_2$) are simulated,
|
|
resulting in $9 \times 108 = 972$ datasets. More precisely, for each dataset,
|
|
we pick uniformly at random two permutations of $\{ 1, \dots , 4 \}$
|
|
($\sigma_1, \sigma_2$) with the constraint that $\sigma_1(4) \neq \sigma_2(1)$.
|
|
This ensures that each of the two networks have a non-empty block that is empty
|
|
in the other one. Then the networks are simulated with
|
|
$\mathcal{B}$ern-$BiSBM_{120}(4, \bm{\alpha}, \bm{\pi}^m, \bm{\rho}^m)$
|
|
with the previous parameters. Each network has 2 blocks in common and their
|
|
connectivity structures encompass a mix of core-periphery, assortative
|
|
community and disassortative community structures, depending on which 3 of the 4
|
|
blocks are selected for each network. $\eps[\alpha]$ represents the strength of
|
|
these structures, the larger, the easier it is to tell apart one block from
|
|
another.
|
|
The true model of all the simulation is a $\pi\rho\text{-}colBiSBM$.
|
|
|
|
\paragraph{Inference} We want to measure the quality of the
|
|
inference procedure, for this we use the inference described in the section
|
|
\ref{sec:variational-estimation-of-the-parameters}.
|
|
|
|
\paragraph{Quality indicators} To assess the quality of the inference, we will
|
|
use the following indicators:
|
|
\begin{itemize}
|
|
\item First, for each dataset, we put in competition $\pi\text{-}colBiSBM$ with
|
|
$sep\text{-}BiSBM$, $iid\text{-}colBiSBM$, $\rho\text{-}colBiSBM$,
|
|
$\pi\rho\text{-}colBiSBM$
|
|
respectively. To do so, for each dataset, we compute the
|
|
BIC-L of each model $\pi\text{-}colBiSBM$ is preferred to $sep\text{-}BiSBM$
|
|
(resp. $iid\text{-}colBiSBM$, $\rho\text{-}colBiSBM$,
|
|
$\pi\rho\text{-}colBiSBM$) if
|
|
its BIC-L is greater.
|
|
\item When considering $\pi\text{-}colBiSBM$, $\rho\text{-}colBiSBM$,
|
|
$\pi\rho\text{-}colBiSBM$ we compare $\widehat{Q_1}$, $\widehat{Q_2}$ to
|
|
their true values. ($Q_1 = 4$ and $Q_2 = 4$)
|
|
\item Finally, we assess the quality of the node grouping by computing the
|
|
Adjusted Rand Index \parencite[][, ARI = 0 for a random grouping, ARI = 1 for a perfect recovery]{hubertComparingPartitions1985}. For each network, for the
|
|
$\pi\text{-}colBiSBM$, $\rho\text{-}colBiSBM$,
|
|
$\pi\rho\text{-}colBiSBM$ we compare the inferred block memberships to the
|
|
real ones by computing the mean of the ARI per axis over the two networks
|
|
\begin{equation*}
|
|
\overline{\text{ARI}}_d = \frac{1}{2} \text{ARI}\big( \text{ARI}(\widehat{\bm{Z}^1_d},\bm{Z}^1_d) + \text{ARI}(\widehat{\bm{Z}^2_d},\bm{Z}^2_d) \big)
|
|
\end{equation*}
|
|
where $d$ is the dimension or axis (i.e., rows, $d=1$, or columns, $d=2$) of
|
|
the block memberships.
|
|
And we compute the ARI of the whole set of nodes to account for block
|
|
pairing between networks
|
|
\begin{equation*}
|
|
\text{ARI}_d = \text{ARI}\big((\widehat{\bm{Z}^1_d},\widehat{\bm{Z}^2_d}),(\bm{Z}^1_d,\bm{Z}^2_d) \big)
|
|
\end{equation*}
|
|
\end{itemize}
|
|
|
|
All these quality indicators are averaged over the 108 datasets. The results are
|
|
provided in the tables \ref{tab:per_model_sep} to \ref{tab:per_model_pirho}. Each line corresponds to the
|
|
108 datasets for a given value of value of $\eps[\alpha]$.
|
|
|
|
<<inference_table, echo = FALSE>>=
|
|
averaged_data <- result_data_frame %>%
|
|
group_by(epsilon_alpha) %>%
|
|
summarise(across(-preferred_model, list("avrg" = meanse))) %>%
|
|
select(-c(2:10))
|
|
averaged_data <- averaged_data %>%
|
|
select(which(!grepl("*_BICL_*", colnames(averaged_data)),
|
|
arr.ind = TRUE))
|
|
@
|
|
|
|
<<function_per_model, echo = FALSE>>=
|
|
dataframe_per_model <- function(model) {
|
|
averaged_data %>%
|
|
select(epsilon_alpha, starts_with(paste0(model, "_")))
|
|
}
|
|
@
|
|
|
|
\tiny
|
|
<<per_model_table, echo = FALSE, results='asis', message=FALSE, warning = FALSE>>=
|
|
for (model in c("sep", "iid", "pi", "rho", "pirho")) {
|
|
kable_colnames <- c(
|
|
"$\\eps[\\alpha]$", #"BIC-L",
|
|
"$\\overline{\\text{ARI}}_{1}$",
|
|
"$\\overline{\\text{ARI}}_{2}$", "$\\text{ARI}_{1}$", "$\\text{ARI}_{2}$"
|
|
)
|
|
model_name <- model
|
|
if (model != "sep") {
|
|
kable_colnames <- c(
|
|
kable_colnames, "Recovered $Q_1$",
|
|
"Recovered $Q_2$"
|
|
)
|
|
}
|
|
if (model == "pirho") {
|
|
model_name <- "$\\pi\\rho$"
|
|
} else {
|
|
if (model != "iid" && model != "sep") {
|
|
model_name <- paste0("$\\", model, "$")
|
|
} else {
|
|
model_name <- paste0("$", model, "$")
|
|
}
|
|
}
|
|
print(kable(dataframe_per_model(model),
|
|
escape = FALSE,
|
|
booktabs = TRUE,
|
|
digits = 2,
|
|
position = "!h",
|
|
caption = paste0(
|
|
"\\label{tab:per_model_", model,
|
|
"}Quality metrics for ",
|
|
ifelse(model != "sep", paste0(model_name, "$\\text{-}colBiSBM$"),"$sep\\text{-}BiSBM$")
|
|
),
|
|
col.names = kable_colnames
|
|
))
|
|
}
|
|
@
|
|
\normalsize
|
|
|
|
<<proportion-preferred_model, echo = FALSE>>=
|
|
proportion_preferred_data <- result_data_frame %>%
|
|
group_by(epsilon_alpha, preferred_model) %>%
|
|
summarise(n = n()) %>%
|
|
mutate(prop_model = n / sum(n)) %>%
|
|
ungroup() %>%
|
|
select(-n)
|
|
|
|
proportion_preferred_table <- proportion_preferred_data %>%
|
|
pivot_wider(
|
|
names_from = preferred_model,
|
|
values_from = prop_model, values_fill = 0
|
|
)
|
|
|
|
kable(proportion_preferred_table,
|
|
escape = FALSE,
|
|
booktabs = TRUE,
|
|
digits = 2,
|
|
position = "!h",
|
|
caption = "\\label{tab:proportion-preferred-table}Proportions of models selected per \\eps[\\alpha] (data for Figure \\ref{fig:inference-proportion-preferred})",
|
|
col.names = c(
|
|
"\\eps[\\alpha]",
|
|
"$sep\\text{-}BiSBM$",
|
|
"$iid\\text{-}colBiSBM$",
|
|
"$\\pi\\text{-}colBiSBM$",
|
|
"$\\rho\\text{-}colBiSBM$",
|
|
"$\\pi\\rho\\text{-}colBiSBM$"
|
|
),
|
|
align = "rccccc",
|
|
format = "latex"
|
|
)
|
|
@
|
|
|
|
<<proportion_preferred_figure, echo = FALSE>>=
|
|
#| fig.cap="\\label{fig:inference-proportion-preferred}Plot of the proportions of different preferred models in function of \\eps[\\alpha]",
|
|
#| fig.asp = 0.5,
|
|
#| fig.pos = "H",
|
|
#| fig.width = 7,
|
|
#| fig.height = 4,
|
|
#| dpi=300
|
|
|
|
plot <- proportion_preferred_data %>% ggplot() +
|
|
aes(
|
|
x = epsilon_alpha, y = prop_model, color = preferred_model,
|
|
fill = preferred_model
|
|
) +
|
|
guides(
|
|
fill = guide_legend(title = "Preferred Model"),
|
|
color = guide_legend(title = "Preferred Model")
|
|
) +
|
|
scale_x_continuous(breaks = seq(from = 0.0, to = 0.24, by = 0.03)) +
|
|
scale_color_okabe_ito() +
|
|
scale_fill_okabe_ito() +
|
|
xlab(TeX("$\\epsilon_{\\alpha}$")) +
|
|
ylab("Model proportions") +
|
|
geom_col(position = "stack")
|
|
print(plot)
|
|
@
|
|
|
|
\paragraph{Results} For the model comparison, when $\eps[\alpha]$ is small
|
|
($\eps[\alpha]\in[0, .04]$), the simulation model is close to the
|
|
Erd\H{o}s-Reńyi network and it is very hard to find any structure beyond the one
|
|
of a single block on each dimension.
|
|
|
|
On the figure \ref{fig:inference-proportion-preferred} and table
|
|
\ref{tab:proportion-preferred-table} we can see that from
|
|
$\eps[\alpha] = 0.12$ around $70\%$ of the time the $\pi\rho\text{-}colBiSBM$
|
|
model (i.e., the correct one) is selected.
|
|
|
|
An interesting result we can read in the tables is that our models outperform
|
|
the $sep\text{-}BiSBM$ when considering the ARI on the whole set of nodes
|
|
($\text{ARI}_d$). This means that our models are able to recover the block
|
|
pairing \emph{between the networks} in addition to recovering the blocks and
|
|
their parameters. |