## ----setup, include=FALSE-----------------------------------------------------
library(knitr)

# Global options
opts_chunk$set(fig.path="fig/")

## ----NetResponse1, warning=FALSE, message=FALSE-------------------------------
library(netresponse)

# Generate simulated data
res <- generate.toydata(Dim = 3, Nc = 3, Ns = 200, sd0 = 3, rgam.shape = 1, rgam.scale = 1)

D <- res$data
component.means <- res$means
component.sds   <- res$sds
sample2comp     <- res$sample2comp

# Use fully connected network
network <- matrix(rep(1, 9), nrow = 3) 

# Fit NetResponse model
# Various network formats are supported, see help(detect.responses) for
# details. With large data sets, consider the 'speedup' option.
set.seed(4243)
res <- detect.responses(D, network, mixture.method = "vdp", pca.basis = TRUE)

# List subnets (each is a list of nodes)
subnet.id <- names(get.subnets(res))[[1]]

## ----NetResponse2, fig.width=6, fig.height=6, warning=FALSE, message=FALSE, fig.show="hide", eval=FALSE----
#  library(ggplot2)
#  vis <- plot_responses(res, subnet.id, plot_mode = "pca")

## ----NetResponse2b, fig.width=6, fig.height=5, warning=FALSE, message=FALSE, eval=FALSE----
#  # Modify the resulting ggplot2 object to enhance visualization
#  p <- vis$p # Pick the ggplot2 object from results
#  p <- p + geom_point(size = 3) # Modify point size
#  print(p)

## ----NetResponse3, fig.width=8, fig.height=8, warning=FALSE, message=FALSE, eval=FALSE----
#  vis <- plot_responses(res, subnet.id, plot_mode = "network")

## ----NetResponse4, fig.width=8, fig.height=8, warning=FALSE, message=FALSE, eval=FALSE----
#  vis <- plot_responses(res, subnet.id, plot_mode = "heatmap")

## ----NetResponse5, fig.width=8, fig.height=8, warning=FALSE, message=FALSE, eval=FALSE----
#  vis <- plot_responses(res, subnet.id, plot_mode = "boxplot_data")

## ----NetResponse7, fig.width=8, fig.height=8, warning=FALSE, message=FALSE, eval=FALSE----
#  plot_scale(vis$breaks, vis$palette, two.sided = TRUE)

## ----NetResponse8, warning=FALSE, message=FALSE-------------------------------
subnet.id <- 'Subnet-1'

# Sample - response probabilities (soft cluster assignment)
response.probs <- sample2response(res, subnet.id)
tail(round(response.probs, 6))

# Sample - response hard assignments
hard.clusters <- response2sample(res, subnet.id)
print(hard.clusters)

## ----NetResponse9, warning=FALSE, message=FALSE-------------------------------
params <- get.model.parameters(res, subnet.id) 
names(params)

## ----vdp, warning=FALSE, message=FALSE----------------------------------------
# Generate 2-dimensional simulated data with 3 clusters
res <- generate.toydata(Dim = 2, Nc = 3, Ns = 200, sd0 = 3, rgam.shape = 1, rgam.scale = 1)

D <- res$data
real.means <- res$means
real.sds   <- res$sds
real.sample2comp     <- res$sample2comp

# Infinite Gaussian mixture model with       
# Variational Dirichlet Process approximation       
mixt <- vdp.mixt( D )
            
# Centroids of the detected Gaussian components       
estimated.means <- mixt$posterior$centroids

# The colors denote the known clusters
# The blue ball denotes the original (known) cluster centroids and
# the triangle denotes the estimated centroids
plot(D, col = real.sample2comp, pch = 1)
points(real.means, col = "blue", pch = 16, cex = 2)
points(estimated.means, col = "blue", pch = 17, cex = 2)
            
# Hard mixture component assignment for each sample
estimated.sample2comp <- apply(mixt$posterior$qOFz, 1, which.max)

# Compare known and estimated mixture components
# (note that cluster indices may have switched due to unidentifiability)
# nearly all samples have one-to-one match between the real and estimated 
# clusters
head(table(estimated.sample2comp, real.sample2comp))

## ----cite, warning=FALSE, message=FALSE---------------------------------------
citation("netresponse")

## ----version, warning=FALSE, message=FALSE------------------------------------
sessionInfo()