---
title: "Overview of using SBIC for network models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(SBICgraph)
library(network) # for visualization 

# to reset par
resetPar <- function() {
    dev.new()
    op <- par(no.readonly = TRUE)
    dev.off()
    op
}
```
The function `comparison` allows for comparison between the true network and the estimated network from the SBIC method.  

First we create a simulated data set using the embedded `simulate` function within SBIC. The function `simulate` generates a data frame, a real network adjacency matrix and a prior network adjacency matrix.  
```{r simulate_data}
p <- 200
m1 <- 100
m2 <- 30
d <- simulate(n=100, p=p, m1=m1, m2=m2)
data<- d$data
real<- d$realnetwork
priori<- d$priornetwork

```
We can visualize the networks  
```{r visualize_networks}
prior_net <- network(priori)
real_net <- network(real)
par(mfrow = c(1,2))
plot(prior_net, main = "Prior network")
plot(real_net, main = "Real network")
par(resetPar())
```

We examine some features of both the prior network and the real network 
```{r examining_networks}
sum(priori[lower.tri(priori)])
sum(priori[lower.tri(priori)])/(p*(p-1)/2)
sum(real[lower.tri(real)])
sum(real[lower.tri(real)])/(p*(p-1)/2)
```

Then we can fit SBIC using one function
```{r fit_models}
lambda<- exp(seq(-10,10, length=30))
# calculating the error rate from the number of edges in the true graph and the number of discordant pairs 
r1 <- m2/m1
r2 <-m2/(p*(p-1)/2-m1)
r <- (r1+r2)/2
model<- sggm(data = data, lambda = lambda, M=priori, prob = r)
```
Comparing the estimated network to the true and prior network. Our comparison function above calcualtes the Positive selection rate (PSR) and the False positive rate (FDR)
```{r}
print("Comparing estimated model with the real network")
comparison(real = real, estimate = model$networkhat)
print("Comparing the prior network with the real network")
comparison(real = real, estimate = priori)
```
We can also compare visualizations 
```{r}
estimated_net <- network(model$networkhat)
par(mfrow = c(1,3))
plot(prior_net, main = "Prior Network")
plot(real_net, main = "Real Network")
plot(estimated_net, main = "Estimated Network")
par(resetPar())
```


The model object also stores all the candidate models generated.  
```{r}
length(model$candidate)
```