# Introduction

This is a short introduction and test for graphical models with stability selection. We first load the relevant packages. We use package QUIC to fit graphical networks:

``````library(stabs)
library(huge) # Need package huge for generating test data
``````
``````## Warning: package 'huge' was built under R version 3.4.1
``````
``````library(QUIC)
library(igraph)
``````

## Generate some test data using huge

``````N <- 200
set.seed(10010)
dat.hubs <- huge.generator(n=N, d=40, graph="hub")
set.seed(10010)
dat.cluster <- huge.generator(n=N, d=40, graph="cluster")
set.seed(10010)
dat.rand <- huge.generator(n=N, d=40, graph="random")
``````

Next, we can visualize the test data

``````plot(dat.hubs)
``````

``````plot(dat.cluster)
``````

``````plot(dat.rand)
``````

# Stability selection

To run stability selection on a graphical model, we can simply use the dedicated `fitfun` called `quic.graphical_model`. Note that we don't supply a `y` argument to stabsel.

``````s.hubs <- stabsel(x = dat.hubs\$data, fitfun = quic.graphical_model,
cutoff = 0.75, PFER = 10)
s.cluster <- stabsel(x = dat.cluster\$data, fitfun = quic.graphical_model,
cutoff = 0.75, PFER = 10)
s.rand <- stabsel(x = dat.rand\$data, fitfun = quic.graphical_model,
cutoff = 0.75, PFER = 10)
``````

## Plot comparisons

In order to make comparisons possible, we use the following user defined plot function

``````plot_comparison <- function(stabsout, orig) {
## display comparison of original graph and stabs estimation
j <- orig\$omega * 0
orig.graph <- graph.adjacency(orig\$theta != 0, mode = "max", diag = FALSE)
ut <- upper.tri(j)
j[ut][stabsout\$selected] <- 1
stabs.graph <- graph.adjacency(j != 0, mode = "max", diag = FALSE)
layout <- layout.fruchterman.reingold(orig.graph)

par(mfrow = c(1,2))
plot(orig.graph, layout = layout, edge.color = "gray50", vertex.color = 'red',
main = "Real graph",  vertex.size = 3, vertex.label = NA)
plot(stabs.graph, layout = layout, edge.color = "gray50", vertex.color = 'red',
main = "Stabs estimated graph", vertex.size = 3, vertex.label = NA)
}
``````

Now, we compare the three graphs:

``````plot_comparison(s.hubs, dat.hubs)
``````

``````plot_comparison(s.cluster, dat.cluster)
``````

``````plot_comparison(s.rand, dat.rand)
``````

As one can see, the original graphs are resembled by the estimated graphs rather well in this example.

# User specified `fitfun` for graphical models

For general instruction on how to create a `fitfun` please refer to

``````vignette("Using_stabs", package = "stabs")
``````

For an example for grahical model see

``````quic.graphical_model
``````
``````## function (x, y, q, ...)
## {
##     if (!requireNamespace("QUIC"))
##         stop("Package ", sQuote("QUIC"), " is required but not available")
##     extraargs <- list(...)
##     empirical.cov <- cov(x)
##     lams <- extraargs\$lams
##     if (is.null(lams)) {
##         max.cov <- max(abs(empirical.cov[upper.tri(empirical.cov)]))
##         lams <- getLamPath(max.cov, max.cov * 0.05, len = 40)
##     }
##     est <- QUIC::QUIC(empirical.cov, rho = 1, path = lams, msg = 0)
##     ut <- upper.tri(empirical.cov)
##     qvals <- sapply(1:length(lams), function(idx) {
##         m <- est\$X[, , idx]
##         sum(m[ut] != 0)
##     })
##     qq <- qvals >= q
##     if (!any(qq)) {
##         stop("Didn't reach the required number of variables. Try supplying lambda manually")
##     }
##     lamidx <- which.max(qvals >= q)
##     M <- est\$X[, , lamidx][ut]
##     selected <- (M != 0)
##     s <- sapply(1:lamidx, function(idx) {
##         m <- est\$X[, , idx][ut] != 0
##         return(m)
##     })
##     colnames(s) <- as.character(1:ncol(s))
##     return(list(selected = selected, path = s))
## }
## <bytecode: 0x00000000107a1970>
## <environment: namespace:stabs>
## attr(,"class")
## [1] "function"        "graphical_model"
``````

A speciality of graphical models is that the function has an additional class `"graphical_model"`. You can set this class immediately after having defined your `fitfun`, say, `my.graphical_model` as follows:

``````class(my.graphical_model) <- c(class(my.graphical_model), "graphical_model")
``````

This will avoid that stabsel expects a `y` argument and will change naming conventions for the resulting selections.