Add stability_score and signal to noise ratio plot functions

2 files changed, 65 insertions, 0 deletions

diff --git a/R/plot_snr.R b/R/plot_snr.R
new file mode 100644
index 0000000..3fb1f9c
--- /dev/null
+++ b/R/plot_snr.R
@@ -0,0 +1,39 @@
+#' Plot Signal to Noise Ratio
+#'
+#' Compute the stability scores for each entry and plot it
+#' against the estimate correlation coefficients.
+#'
+#' @param X The data matrix to operate on
+#' @param method Required to select a correlation method
+#'
+#' @return ggplot object
+plot_snr <- function(
+	X,
+	method = c('rmad', 'pearson', 'spearman', 'kendall')
+) {
+	method <- match.arg(method)
+	# Compute reference correlation
+	R <- if(method == 'rmad'){
+		rmad(X)
+	}else{
+		cor(X, method = method)
+	}
+
+	W <- stability_score(X)
+	
+	idx <- which(upper.tri(W), arr.ind = TRUE)
+	w <- W[upper.tri(W)]
+	r <- R[upper.tri(R)]
+
+	df <- data.frame(r = r, w = w)
+
+	ggplot2::ggplot(data = df, ggplot2::aes(x = r, y = w)) +
+	ggplot2::geom_point(alpha = 0.5) +
+	ggplot2::geom_smooth(se = FALSE) +
+	ggplot2::labs(
+		title = paste('Weighted vs', method, 'correlation'),
+		x = paste(method, 'correlation'),
+		y = 'Weight'
+	) +
+	ggplot2::theme_classic()
+}
diff --git a/R/stability_score.R b/R/stability_score.R
new file mode 100644
index 0000000..7841214
--- /dev/null
+++ b/R/stability_score.R
@@ -0,0 +1,26 @@
+#' Stability Score
+#'
+#' Computing entry wise stability score for selection
+#'
+#' @param X The data matrix to operate on
+#' @param K Number of k-folds (default 25)
+#' @param subsample_fraction Hold out fraction (default 0.7)
+#'
+#' @return Score array
+#' @export
+stability_score <- function(X, K = 25, subsample_fraction = 0.7){
+	n <- nrow(X)
+	p <- ncol(X)
+	rho <- vector('list', K)
+	for(r in seq_len(K)){
+		idx <- sample(n, size = floor(subsample_fraction * n))
+		rho[[r]] <- cor(X[idx, ])
+	}
+	R <- simplify2array(rho)
+	theta <- acos(pmax(pmin(R, 1), -1))
+	delta <- theta - pi/2
+	delta_sd <- apply(delta, c(1, 2), sd)
+	delta_sd <- pmax(delta_sd, quantile(delta_sd[upper.tri(delta_sd)], 0.05))
+	score <- 1 / delta_sd
+	score / median(score[upper.tri(score)])
+}