03_analysis_real_dataset.Rmd

---
title: "Stochastic Block Model Prior with Ordering Constraints for Gaussian Graphical Models"
author:
    - Alessandro Colombi (Supervisor)^[a.colombi10@campus.unimib.it]
    - Teo Bucci^[teo.bucci@mail.polimi.it]
    - Filippo Cipriani^[filippo.cipriani@mail.polimi.it]
    - Filippo Pagella^[filippo.pagella@mail.polimi.it]
    - Flavia Petruso^[flavia.petruso@mail.polimi.it]
    - Andrea Puricelli^[andrea3.puricelli@mail.polimi.it]
    - Giulio Venturini^[giulio.venturini@mail.polimi.it]
output:
    pdf_document:
        toc: true
        toc_depth: 3
        number_section: true
    html_document:
        toc: true
        toc_float: true
        number_sections: true
#date: "2023-01-17"
editor_options:
    chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align='center')
```

# Simulations

## Load the necessary packages

```{r, include = FALSE}
shhh = function(lib_name){ # It's a library, so shhh!
    suppressWarnings(suppressMessages(require(lib_name, character.only = TRUE)))
}
shhh("tidyverse")
shhh("ACutils")
shhh("mvtnorm")
shhh("salso")
shhh("FGM")
shhh("gmp")
shhh("mcclust")
shhh("mcclust.ext")
shhh("logr")
shhh("tidygraph")
shhh("ggraph")
shhh("igraph")
```

```{r, include = FALSE}
paths = c(
    "src/utility_functions.R",
    "src/bulky_functions.R",
    "src/data_generation.R"
)

for(p in paths){
    path = file.path(p)
    if(file.exists(path)){
        source(path)
    } else {
        cat("File",path,"was not found in directory, please check.")
    }
}
```

## Load data

```{r}
suono_centered <- load(here::here("dataset","Suono_centered.Rdat"))
suono_normalized <- load(here::here("dataset","Suono_Quantile_Normalized.Rdat"))

n = nrow(suono_cent)
p = ncol(suono_quant_norm)
```

## Set options for the simulation

```{r}
options = set_options(sigma_prior_0=0.5,
                      sigma_prior_parameters=list("a"=1,"b"=1,"c"=1,"d"=1),
                      theta_prior_0=1,
                      theta_prior_parameters=list("c"=1,"d"=1),
                      rho0=p, # start with one group
                      weights_a0=rep(1,p-1),
                      weights_d0=rep(1,p-1),
                      alpha_target=0.234,
                      beta_mu=0.25,
                      beta_sig2=1/16,
                      d=3,
                      alpha_add=0.5,
                      adaptation_step=1/(p*1000),
                      update_sigma_prior=TRUE,
                      update_theta_prior=TRUE,
                      update_weights=TRUE,
                      update_partition=TRUE,
                      update_graph=TRUE,
                      perform_shuffle=TRUE)
```


## Run the simulation

```{r}
res <- Gibbs_sampler(
    data = suono_cent,
    niter = 8000,
    nburn = 2000,
    thin = 1,
    options = options,
    seed = 123456,
    print = FALSE
)

res_norm <- Gibbs_sampler(
    data = suono_quant_norm,
    niter = 8000,
    nburn = 2000,
    thin = 1,
    options = options,
    seed = 123456,
    print = FALSE
)
```

```{r}
# create output directory if not present
dir.create(file.path("output", "real_dataset"), showWarnings = FALSE, recursive = TRUE)

# save an object to a file
filename_data = file.path("output","real_dataset","realtest_cent.rds")
saveRDS(res, file = filename_data)

filename_data = file.path("output","real_dataset","realtest_norm.rds")
saveRDS(res_norm, file = filename_data)
```


# Posterior analysis (Case: Normalized)

Restore the object

```{r}
filename_data = file.path("output","real_dataset","realtest_norm.rds")
res = readRDS(file = filename_data)
```

## Partition

Recomputing the partition in other forms and the number of groups

```{r}
rho = res$rho
r = do.call(rbind, lapply(res$rho, rho_to_r))
z = do.call(rbind, lapply(res$rho, rho_to_z))
num_clusters = do.call(rbind, lapply(res$rho, length))
num_clusters = as.vector(num_clusters)
```

### Acceptance frequency

```{r}
mean(res$accepted)
```


### Barplot of changepoints

```{r}
bar_heights = colSums(r)
barplot(
    bar_heights,
    names = seq_along(bar_heights),
    border = "NA",
    space = 0,
    yaxt = "n",
    main="Changepoint frequency distribution",
    #col = color,
    cex.names=.6,
    las=2
)
```

### Evolution of the number of clusters

```{r}
plot(
    x = seq_along(num_clusters),
    y = num_clusters,
    type = "n",
    xlab = "Iterations",
    ylab = "Number of groups",
    main = "Number of groups - Traceplot"
)
lines(x = seq_along(num_clusters), y = num_clusters)
```

```{r}
barplot(
    prop.table(table(num_clusters)),
    xlab = "Number of groups",
    ylab = "Relative Frequency",
    main = paste(
        "Number of groups - Relative Frequency\n",
        "Last:",
        tail(num_clusters, n = 1),
        "- Mean:",
        round(mean(num_clusters), 2)
    )
)
```

### Retrieving best partition using VI on visited ones (order is guaranteed here)

Here we are satisfied with finding the optimal partition only in the set of those visited, not in all the possible ones. I expect it could work even worse. But at least it guarantees to find an admissible one.
I would say that it is the implementation of formula (13) of the Corradin-Danese paper (https://doi.org/10.1016/j.ijar.2021.12.019).

```{r}
library("Rcpp")
library("RcppArmadillo")
sourceCpp("src/wade.cpp")

# compute VI
sim_matrix <- salso::psm(z)
dists <- VI_LB(z, psm_mat = sim_matrix)

# select best partition (among the visited ones)
best_partition_index = which.min(dists)
rho_est = rho[[best_partition_index]]
z_est = z[best_partition_index,]

# VI loss
dists[best_partition_index]

# select best partition
unname(z_est)
```

## Graph

Extract last plinks

```{r}
last_plinks = tail(res$G, n=1)[[1]]
```

Criterion 1 to select the threshold (should not work very well) and assign final graph

```{r}
threshold = 0.5
G_est <- matrix(0,p,p)
G_est[which(last_plinks>threshold)] = 1
```

Criterion 2 to select the threshold

```{r}
bfdr_select = BFDR_selection(last_plinks, tol = seq(0.1, 1, by = 0.001))
```

Inspect the threshold and assign final graph

```{r}
bfdr_select$best_treshold
G_est = bfdr_select$best_truncated_graph
```

### Plot estimated matrices


```{r}
ACutils::ACheatmap(
    last_plinks,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated plinks matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    G_est,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Graph",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    tail(res$K,n=1)[[1]],
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Precision matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)
```

### Graph visualization

Plot against the original

```{r}
g2 <- graph.adjacency(G_est)
edges2 <- get.edgelist(g2)
edges2 <- cbind(edges2, rep("estimated", nrow(edges2)))
edges <- as.data.frame(edges2)
names(edges) = c("from", "to", "graph")
nodes = data.frame(
    vertices = 1:p,
    clust_est = as.factor(z_est)
)
# nodes
g = graph_from_data_frame(edges, directed = FALSE, nodes)
lay = create_layout(g, layout = "linear", circular = TRUE)

output_plot <- ggraph(lay) +
    geom_edge_arc(edge_colour = "grey") +
    geom_node_point(aes(color = clust_est), size = 4) +
    geom_node_text(aes(label = name), repel = TRUE)
print(output_plot)
```




# Posterior analysis (Case: Centered)

Restore the object

```{r}
filename_data = file.path("output","real_dataset","realtest_cent.rds")
res = readRDS(file = filename_data)
```

## Partition

Recomputing the partition in other forms and the number of groups

```{r}
rho = res$rho
r = do.call(rbind, lapply(res$rho, rho_to_r))
z = do.call(rbind, lapply(res$rho, rho_to_z))
num_clusters = do.call(rbind, lapply(res$rho, length))
num_clusters = as.vector(num_clusters)
```

### Acceptance frequency

```{r}
mean(res$accepted)
```


### Barplot of changepoints

```{r}
bar_heights = colSums(r)
barplot(
    bar_heights,
    names = seq_along(bar_heights),
    border = "NA",
    space = 0,
    yaxt = "n",
    main="Changepoint frequency distribution",
    #col = color,
    cex.names=.6,
    las=2
)
```

### Evolution of the number of clusters

```{r}
plot(
    x = seq_along(num_clusters),
    y = num_clusters,
    type = "n",
    xlab = "Iterations",
    ylab = "Number of groups",
    main = "Number of groups - Traceplot"
)
lines(x = seq_along(num_clusters), y = num_clusters)
```

```{r}
barplot(
    prop.table(table(num_clusters)),
    xlab = "Number of groups",
    ylab = "Relative Frequency",
    main = paste(
        "Number of groups - Relative Frequency\n",
        "Last:",
        tail(num_clusters, n = 1),
        "- Mean:",
        round(mean(num_clusters), 2)
    )
)
```

### Retrieving best partition using VI on visited ones (order is guaranteed here)

Here we are satisfied with finding the optimal partition only in the set of those visited, not in all the possible ones. I expect it could work even worse. But at least it guarantees to find an admissible one.
I would say that it is the implementation of formula (13) of the Corradin-Danese paper (https://doi.org/10.1016/j.ijar.2021.12.019).

```{r}
library("Rcpp")
library("RcppArmadillo")
sourceCpp("src/wade.cpp")

# compute VI
sim_matrix <- salso::psm(z)
dists <- VI_LB(z, psm_mat = sim_matrix)

# select best partition (among the visited ones)
best_partition_index = which.min(dists)
rho_est = rho[[best_partition_index]]
z_est = z[best_partition_index,]

# VI loss
dists[best_partition_index]

# select best partition
unname(z_est)
```

## Graph

Extract last plinks

```{r}
last_plinks = tail(res$G, n=1)[[1]]
```

Criterion 1 to select the threshold (should not work very well) and assign final graph

```{r}
threshold = 0.5
G_est <- matrix(0,p,p)
G_est[which(last_plinks>threshold)] = 1
```

Criterion 2 to select the threshold

```{r}
bfdr_select = BFDR_selection(last_plinks, tol = seq(0.1, 1, by = 0.001))
```

Inspect the threshold and assign final graph

```{r}
bfdr_select$best_treshold
G_est = bfdr_select$best_truncated_graph
```

### Plot estimated matrices


```{r}
ACutils::ACheatmap(
    last_plinks,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated plinks matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    G_est,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Graph",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    tail(res$K,n=1)[[1]],
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Precision matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)
```

### Graph visualization

Plot against the original

```{r}
g2 <- graph.adjacency(G_est)
edges2 <- get.edgelist(g2)
edges2 <- cbind(edges2, rep("estimated", nrow(edges2)))
edges <- as.data.frame(edges2)
names(edges) = c("from", "to", "graph")
nodes = data.frame(
    vertices = 1:p,
    clust_est = as.factor(z_est)
)
# nodes
g = graph_from_data_frame(edges, directed = FALSE, nodes)
lay = create_layout(g, layout = "linear", circular = TRUE)

output_plot <- ggraph(lay) +
    geom_edge_arc(edge_colour = "grey") +
    geom_node_point(aes(color = clust_est), size = 4) +
    geom_node_text(aes(label = name), repel = TRUE)
print(output_plot)
```