`%notin%` <- function(a, b) !(a %in% b)
# get package names
get_package_names <- function() {
gsub("package:", "", search()[grep("package:", search())])
}
#' Is the relative error between two numbers is less than the given tolerance?
#'
#' Given two consecutive numbers in a sequence, return \code{TRUE} if the
#' relative change is positive but less than the given tolerance.
#' @param current A scalar; the most recent value of the sequence.
#' @param previous A scalar; the second most recent value of the sequence, or a
#' reference value.
#' @param tol The tolerance, a positive scalar near zero.
#' @keywords internal
has_converged <- function(current, previous, tol = 0.01) {
rel_change <- (current - previous) / abs(previous)
rel_change > 0 && rel_change < tol
}
#' Flip a matrix upside down
#' @param x A matrix
#' @return A matrix, \code{x} with rows reversed.
#' @keywords internal
flipud <- function(x) {
x[rev(seq_len(nrow(x))), , drop = FALSE]
}
#' Convert a long data frame to a wide matrix.
#'
#' Convert a long data frame to a wide matrix, with time along the row dimension
#' and locations along the column dimension. Values in the matrix could be e.g.
#' the observed counts or the population.
#' @param df A data frame with at least 3 columns.
#' @param time_col Integer or string that specifies the time column.
#' @param location_col Integer or string that specifies the location column.
#' @param value_col Integer or string that specifies the value column.
#' @return A matrix with time on rows and locations on columns.
#' @importFrom dplyr select
#' @importFrom tidyr spread_
#' @importFrom magrittr %>%
#' @export
df_to_matrix <- function(df, time_col = 1, location_col = 2, value_col = 3) {
if (is.numeric(time_col)) {
time_col <- names(df)[time_col]
}
if (is.numeric(location_col)) {
location_col <- names(df)[location_col]
}
if (is.numeric(value_col)) {
value_col <- names(df)[value_col]
}
x <- df %>%
as.data.frame %>%
select(time_col, location_col, value_col) %>%
spread_(key_col = location_col, value_col = value_col)
times <- x[, 1]
x <- as.matrix(x[, -1])
attributes(x)$dimnames[[1]] <- times
return(x)
}
#' Convert a matrix to a data frame.
#'
#' Convert a matrix to a data frame with columns time, location, and one more
#' containing the elements of the input matrix.
#' @param mat A matrix.
#' @param name The name of the third column in the output matrix.
#' @param locations If not \code{NULL}, a vector with the names of the
#' locations.
#' @param times If not \code{NULL}, a vector with the time points. If
#' \code{NULL}, the matrix is assumed to be ordered with time point 1 in the
#' first row.
#' @return A matrix with columns \code{time, location, name}, where \code{name}
#' is specified in the input.
#' @keywords internal
matrix_to_df <- function(mat, name, locations = NULL, times = NULL) {
if (!is.null(locations) && length(locations) != ncol(mat)) {
stop("The number of locations must be equal to col(mat)")
}
if (!is.null(times) && length(times) != nrow(mat)) {
stop("The number of times must be equal to nrow(mat)")
}
if (is.null(locations)) {
locations <- seq_len(ncol(mat))
}
if (is.null(times)) {
times <- seq_len(nrow(mat))
}
df <- data.frame(time = rep(times, ncol(mat)),
location = rep(locations, each = nrow(mat)),
V1 = as.vector(mat))
names(df)[3] <- name
return(df)
}
#' Run a scan statistic analysis.
#'
#' Run a scan statistic analysis with the given scan statistic and arguments.
#' @param scanstat A scan statistic function.
#' @param args A named list of arguments to be passed to \code{scanstat}.
#' @param gumbel Logical: should a Gumbel P-value be calculated? Default is
#' \code{FALSE}.
#' @return A list with components
#' \describe{
#' \item{observed}{The table of observed statistics.}
#' \item{simulated}{The table of simulated statistics.}
#' \item{MC_pvalue}{The Monte Carlo P-value of the scan statistic.}
#' \item{Gumbel_pvalue}{The Gumbel P-value of the scan statistic.}
#' }
#' @keywords internal
run_scan <- function(scanstat, args, gumbel = FALSE) {
scan <- do.call(scanstat, args)
# Extract the most likely cluster (MLC)
scan$observed %<>% arrange(-score)
MLC <- scan$observed[1, ]
# Get P-values
gumbel_pvalue <- NULL
MC_pvalue <- NULL
if (nrow(scan$simulated) > 0) {
if (gumbel) {
gumbel_pvalue <- gumbel_pvalue(MLC$score, scan$simulated$score,
method = "ML")$pvalue
}
MC_pvalue <- mc_pvalue(MLC$score, scan$simulated$score)
}
return(list(observed = scan$observed,
replicates = scan$simulated,
MC_pvalue = MC_pvalue,
Gumbel_pvalue = gumbel_pvalue))
}
# Clean up when package is unloaded.
.onUnload <- function (libpath) {
library.dynam.unload("scanstatistics", libpath)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.