Fixed opts bug and polished build/check - Thanks for pointing this ou…

…t ! - This closes #2

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: netresponse
 Type: Package
 Title: Functional Network Analysis
-Version: 1.19.0
-Date: 2015-03-24
+Version: 1.21.11
+Date: 2016-01-09
 Author: Leo Lahti, Olli-Pekka Huovilainen, Antonio Gusmao and Juuso Parkkinen
 Maintainer: Leo Lahti <[email protected]>
 Description: Algorithms for functional network analysis. Includes an
@@ -15,7 +15,7 @@ Depends:
     methods,
     minet,
     mclust,
-    reshape
+    reshape2
 Imports:
     dmt,
     ggplot2,
@@ -30,3 +30,4 @@ BugReports: https://github.com/antagomir/netresponse/issues
 biocViews: CellBiology, Clustering, GeneExpression, Genetics, Network,
     GraphAndNetwork, DifferentialExpression, Microarray,
     Transcription
+RoxygenNote: 5.0.0

NAMESPACE

@@ -1,4 +1,4 @@
-# Generated by roxygen2 (4.1.0): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 export(ICMg.combined.sampler)
 export(ICMg.get.comp.memberships)
@@ -42,13 +42,13 @@ export(model.stats)
 export(order.responses)
 export(pick.model.pairs)
 export(pick.model.parameters)
-export(plot.data)
-export(plot.expression)
-export(plot.response)
-export(plot.responses)
-export(plot.scale)
-export(plot.subnet)
+export(plot_data)
+export(plot_expression)
 export(plot_matrix)
+export(plot_response)
+export(plot_responses)
+export(plot_scale)
+export(plot_subnet)
 export(read.sif)
 export(response.enrichment)
 export(response2sample)
@@ -66,7 +66,6 @@ import(graph)
 import(mclust)
 import(methods)
 import(parallel)
-import(qvalue)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_bar)
@@ -77,7 +76,7 @@ importFrom(ggplot2,geom_line)
 importFrom(ggplot2,geom_point)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggtitle)
-importFrom(ggplot2,opts)
+importFrom(ggplot2,qplot)
 importFrom(ggplot2,theme_bw)
 importFrom(ggplot2,theme_set)
 importFrom(ggplot2,xlab)
@@ -89,5 +88,6 @@ importFrom(igraph,igraph.to.graphNEL)
 importFrom(igraph,subgraph)
 importFrom(minet,build.mim)
 importFrom(plyr,ddply)
-importFrom(reshape,melt)
+importFrom(qvalue,qvalue)
+importFrom(reshape2,melt)
 useDynLib(netresponse)

R/PlotMixture.R

@@ -1,17 +1,13 @@
-#' PlotMixture
-#' 
-#' Arguments:
+#' @title Plot mixtures
+#' @description Plot mixtures.
 #' @param x data vector
 #' @param qofz Mode assignment probabilities for each sample. Samples x modes.
 #' @param binwidth binwidth for histogram
 #' @param xlab.text xlab.text
 #' @param ylab.text ylab.text
 #' @param title.text title.text
-#'
-#' Return:
 #' @return Used for its side-effects
 #' @export
-#'
 #' @author Leo Lahti \email{leo.lahti@@iki.fi}
 #' @references See citation("netresponse") for citation details.
 #' @keywords utilities

R/PlotMixtureUnivariate.R

@@ -1,9 +1,5 @@
-#' PlotMixtureUnivariate
-#' 
-#' Visualize data, centroids and stds for a given univariate
-#' Gaussian mixture model with PCA.
-#'
-#' Arguments:
+#' @title Plot univariate mixtures
+#' @description Visualize data, centroids and stds for a given univariate Gaussian mixture model with PCA.
 #' @param x data vector
 #' @param means mode centroids
 #' @param sds mode standard deviations
@@ -16,12 +12,8 @@
 #' @param density.color Color for density lines
 #' @param cluster.assignments Vector of cluster indices, indicating cluster for each data point
 #' @param ... Further arguments for plot function.
-#'
-#' Return:
 #' @return Used for its side-effects
-#'
 #' @export
-#'
 #' @author Leo Lahti \email{leo.lahti@@iki.fi}
 #' @references See citation("netresponse") for citation details.
 #' @keywords utilities

R/add.ellipse.R

@@ -0,0 +1,28 @@
+#' @title Add ellipse to an existing plot
+#' @description Calculates and plots ellipse corresponding to specified confidence interval in 2-dimensional plot
+#' @usage add.ellipse(centroid, covmat, confidence = 0.95, npoints = 100, col =
+#' "black", ...)
+#' @param centroid Vector with two elements defining the ellipse centroid.
+#' @param covmat Covariance matrix for the investigated data. Only diagonal
+#' covariances supported.
+#' @param confidence Confidence level determining the ellipse borders based on
+#' the covariance matrix.
+#' @param npoints Number of plotting points.
+#' @param col Color.
+#' @param ... Other arguments to be passed.
+#' @return Used for plotting side effects.
+#' @author Leo Lahti \email{leo.lahti@@iki.fi}
+#' @keywords utilities
+#' @export
+#' @examples #add.ellipse(centroid = c(0, 0), covmat = diag(c(1,2)))
+add.ellipse <- function (centroid, covmat, confidence = 0.95, npoints = 100, col = "black", ...) {
+
+  # add ellipse to a plot 
+  el <- ellipse(centroid, covmat, confidence, npoints)
+  points(el, type = "l", col = col, ...)
+
+  el
+}
+
+
+

R/bic.mixture.R

@@ -0,0 +1,27 @@
+#' @title BIC mixture
+#' @description Latent class analysis based on (infinite) Gaussian mixture model. If the input is data matrix, a multivariate model is fitted; if the input is a vector, a univariate model is fitted
+#' @param x samples x features matrix for multivariate analysis, or a vector for univariate analysis 
+#' @param max.modes Maximum number of modes to be checked for mixture model selection
+#' @param bic.threshold BIC threshold which needs to be exceeded before a new mode is added to the mixture.
+#' @param min.modes minimum number of modes
+#' @param ... Further optional arguments to be passed
+#' @return Fitted latent class model (parameters and free energy)
+#' @export
+#' @references See citation("netresponse") 
+#' @author Contact: Leo Lahti \email{leo.lahti@@iki.fi}
+#' @keywords utilities
+bic.mixture <- function (x, max.modes, bic.threshold = 0, min.modes = 1, ...) { 
+
+   # x; max.modes = max.responses; bic.threshold = bic.threshold; min.modes = min.responses
+
+  if (!is.vector(x) && ncol(x) == 1) {x <- x[,1]}	    
+
+  if (is.vector(x)) {
+    bic.mixture.univariate(x, max.modes, bic.threshold, min.modes = min.modes, ...)
+  } else {
+    bic.mixture.multivariate(x, max.modes, bic.threshold, min.modes = min.modes, ...)
+  }
+
+}
+
+

R/bic.mixture.multivariate.R

@@ -0,0 +1,42 @@
+#' @title Multivariate BIC mixture
+#' @description Latent class analysis based on (infinite) Gaussian mixture model. If the input (dat) is data matrix, a multivariate model is fitted. 
+#' @param x matrix (for multivariate analysis)
+#' @param max.modes Maximum number of modes to be checked for mixture model selection
+#' @param bic.threshold BIC threshold which needs to be exceeded before a new mode is added to the mixture.
+#' @param min.modes Minimum number of modes to be checked for mixture model selection
+#' @param ... Further optional arguments to be passed
+#' @return Fitted latent class model (parameters and free energy)
+#' @export
+#' @references See citation("netresponse") 
+#' @author Contact: Leo Lahti \email{leo.lahti@@iki.fi}
+#' @keywords utilities
+bic.mixture.multivariate <- function (x, max.modes, bic.threshold = 0, min.modes = 1, ...) { 
+
+  # x <- mat; max.modes = params$max.responses; bic.threshold = params$bic.threshold
+
+  best.mode <- bic.select.best.mode(x, max.modes, bic.threshold, min.modes) 
+
+  mcl <- Mclust(x, G = best.mode)
+
+  bic <- try(-mclustBIC(x, G = best.mode)[, "VVV"]) 
+  if ( is.na(bic) ) { bic <- Inf } # infinitely bad = Inf
+
+  means <- t(mcl$parameters$mean)
+  vars <- t(apply(mcl$parameters$variance$sigma, 3, function(x){diag(x)}))
+  sds <- sqrt(vars)
+  ws <- as.vector(mcl$parameters$pro)
+  if (is.null(ws)) {ws <- 1} 
+
+  Nparams <- prod(dim(means)) + prod(dim(sds)) + length(ws) 
+
+  # Determine the most likely mode for each sample (-> hard clusters)
+  qofz <- P.r.s(t(x), list(mu = means, sd = sds, w = ws), log = FALSE)
+  rownames(qofz) <- rownames(x)
+  colnames(qofz) <- paste("Mode", 1:ncol(qofz), sep = "-")
+
+  rownames(means) <- rownames(sds) <- names(ws) <- paste("Mode", 1:length(ws), sep = "-")
+  colnames(means) <- colnames(sds) <- colnames(x)
+
+  list(means = means, sds = sds, ws = ws, Nparams = Nparams, free.energy = -mcl$loglik, qofz = qofz, bic = bic)
+
+}

R/bic.mixture.univariate.R

@@ -0,0 +1,47 @@
+#' @title Univariate BIC mixture
+#' @description Latent class analysis based on (infinite) Gaussian mixture
+#' model. If the input (dat) is data matrix, a multivariate model is fitted. If
+#' the input is a vector or a 1-dimensional matrix, a univariate model is
+#' fitted.
+#' @param x  dat vector (for univariate analysis) or a matrix (for multivariate analysis)
+#' @param max.modes Maximum number of modes to be checked for mixture model selection
+#' @param bic.threshold BIC threshold which needs to be exceeded before a new mode is added to the mixture.
+#' @param min.modes minimum number of modes
+#' @param ... Further optional arguments to be passed
+#' @return Fitted latent class model (parameters and free energy)
+#' @author Contact: Leo Lahti \email{leo.lahti@@iki.fi}
+#' @references See citation("netresponse")
+#' @export
+#' @keywords utilities
+bic.mixture.univariate <- function (x, max.modes, bic.threshold = 0, min.modes = 1, ...) { 
+
+  # x <- datamatrix[, node];  max.modes = params$max.responses; bic.threshold = params$bic.threshold
+
+  best.mode <- bic.select.best.mode(x, max.modes, bic.threshold, min.modes = min.modes) 
+  mcl <- Mclust(x, G = best.mode)
+
+  means <- as.vector(mcl$parameters$mean)
+  sds <- as.vector(sqrt(mcl$parameters$variance$sigmasq))
+  if (length(sds) == 1) {sds <- rep(sds, length(means))} 
+  ws <- as.vector(mcl$parameters$pro)
+
+  if (is.null(ws)) {warning("NULL weights, replacing with 1"); ws <- 1} 
+  if (is.null(means)) {warning("NULL means, replacing with 1"); means <- 1} 
+  if (is.null(sds)) {warning("NULL sds, replacing with 1"); sds <- 1} 
+
+  Nparams <- length(means) + length(sds) + length(ws) 
+
+  means <- matrix(means, nrow = length(ws))
+  sds <- matrix(sds, nrow = length(ws))
+
+  # Determine the most likely mode for each sample (-> hard clusters)
+  # save(means, sds, ws, x, file = "~/tmp/tmp.RData")
+  qofz <- P.r.s(matrix(x, nrow = 1), list(mu = means, sd = sds, w = ws), log = FALSE)
+  rownames(qofz) <- names(x)
+
+  names(means) <- names(sds) <- names(ws) <- paste("Mode", 1:length(ws), sep = "-")
+
+  list(means = means, sds = sds, ws = ws, Nparams = Nparams, free.energy = -mcl$loglik, qofz = qofz)
+
+}
+

R/bic.select.best.mode.R

@@ -0,0 +1,65 @@
+#' @title Select best mode with BIC
+#' @description Select optimal number of mixture components by adding components until 
+#' the increase in objective function is below threshold.
+#' @param x  dat vector (for univariate analysis) or a matrix (for multivariate analysis)
+#' @param max.modes Maximum number of modes to be checked for mixture model selection
+#' @param bic.threshold BIC threshold which needs to be exceeded before a new mode is added to the mixture.
+#' @param min.modes Optiomal. Minimum number of modes.
+#' @return Fitted latent class model (parameters and free energy)
+#' @author Contact: Leo Lahti \email{leo.lahti@@iki.fi}
+#' @references See citation("netresponse")
+#' @export
+#' @keywords utilities
+bic.select.best.mode <- function (x, max.modes, bic.threshold, min.modes = 1) {
+
+  # Cost for single mode
+  # BIC : smaller is better
+  # mclustBIC returns the value for -BIC, to be exact
+  nc <- min.modes
+  if (is.vector(x)) { # univariate
+    m <- -mclustBIC(x, G = nc)[, "V"] 
+  } else { # multivariate
+    m <- -mclustBIC(x, G = nc)[, "VVV"] # BIC : smaller is better
+  }
+
+  # ----------------------------------------------------------------
+
+  add.component <- TRUE
+  best.mode <- min.modes
+  if (max.modes == min.modes) {
+    add.component <- FALSE
+  }
+
+  while (add.component && nc < max.modes) {
+
+    nc <- nc + 1
+
+    # BIC : smaller is better
+    if (is.vector(x)) { # univariate
+      m.new <- try(-mclustBIC(x, G = nc)[, "V"]) 
+    } else { # multivariate
+      m.new <- try(-mclustBIC(x, G = nc)[, "VVV"]) 
+    }
+    if ( is.na(m.new) ) { m.new <- Inf } # infinitely bad = Inf
+
+    # FIXME: compressing data with PCA after dimensionality gets otherwise too high?
+    # with around ncol(x) = 30 the mclustBIC is starting to produce NAs
+
+    # FIXME: remove this when code works ok
+    # if (is.na(m.new)) {save(x, nc, file = "m.new.RData")}
+
+    bic.delta <- m.new - m
+
+    if (bic.delta < -bic.threshold) { 
+      best.mode <- nc 
+      m <- m.new
+    } else {
+      add.component <- FALSE
+    }
+  }
+
+  best.mode
+
+}
+
+

R/continuous.responses.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2013 Leo Lahti
+# Copyright (C) 2010-2016 Leo Lahti
 # Contact: Leo Lahti <[email protected]>
 #
 # This program is free software; you can redistribute it and/or modify
@@ -15,19 +15,14 @@
 # "To invent, you need a good imagination and a pile of junk." 
 #     	      	       	      -- Thomas Edison
 
-
-#' Description: Quantify association between modes and continuous variable
-#' 
-#' Arguments:
-#'   @param annotation.vector annotation vector with discrete factor levels, and named by the samples
-#'   @param model NetResponse model object
-#'   @param method method for enrichment calculation
-#'   @param min.size minimum sample size for a response 
-#'   @param data data matrix (samples x features)
-#'
-#' Returns:
-#'   @return List with each element corresponding to one variable and listing the responses according to association strength
-#'            
+#' @title Continuous responses
+#' @description Quantify association between modes and continuous variable
+#' @param annotation.vector annotation vector with discrete factor levels, and named by the samples
+#' @param model NetResponse model object
+#' @param method method for enrichment calculation
+#' @param min.size minimum sample size for a response 
+#' @param data data matrix (samples x features)
+#' @return List with each element corresponding to one variable and listing the responses according to association strength
 #' @author Contact: Leo Lahti \email{leo.lahti@@iki.fi}
 #' @references See citation("netresponse")
 #' @export

R/detect.responses.R

@@ -1,5 +1,3 @@
-
-
   ######################################################################
 
   # Before I put a sketch on paper, the whole idea is worked out
@@ -338,10 +336,10 @@ detect.responses <- function(datamatrix,
   }
 
   # FIXME: if all nodes will be combined (merging.threshold = -Inf), there will be an error. Fix.
-  #'  costs: cost function values at each state
-  #'  moves: indices of groups joined at each state in its columns
-  #'  groupings: groupings at each level of the hierarchy
-  #'  models: compressed representations of the models from each step
+  #  costs: cost function values at each state
+  #  moves: indices of groups joined at each state in its columns
+  #  groupings: groupings at each level of the hierarchy
+  #  models: compressed representations of the models from each step
 
   model <- new("NetResponseModel",
       moves = matrix(move.cost.hist, 3),