#' Summary functions for one pattern.
#'
#' Uses the function \code{\link{summ_func_random_labelling}}. Eats the
#' number of permutations the user has given and replaces it with zero.
#'
#' @seealso summ_func_random_labelling
#' @inheritParams summ_func_random_labelling
#' @param nsim Ignored.
#' @return A matrix with dimensions: summ_func, r.
#' @importFrom abind adrop
#' @export
summ_func <- function(..., nsim = 0L) {
res <- summ_func_random_labelling(..., nsim = 0L)
# There are no simulations.
res[['a']] <- adrop(res[['a']], drop = 1L)
res
}
#' Summary functions for random labelling
#'
#' Calculates summary functions for a pattern and its simulations under the
#' hypothesis of random labelling.
#'
#' @param pattern A \code{\link[spatstat]{ppp}} object as the simple marked
#' point pattern to be analysed. The marks need to be in the form of a
#' numeric vector. The window has to have the type "rectangle".
#' @param edge_correction The name of the edge correction to be used. Options
#' are 'translate' and 'none'.
#' @param mtf_name A vector of mark test function names. "1" stands for the
#' unmarked K-function. Accepted values are '1', 'm', 'mm', 'gamma',
#' 'gammaAbs' and 'morAbs'.
#' @param nsim The number of permutations.
#' @param r_max A positive scalar value representing the maximum radius that
#' should be considered. r_vec overrides r_max. By default, r_max is NULL
#' and will get a sensible default.
#' @param r_vec A monotonically increasing vector of non-negative r-values
#' to act as the endpoints of the bins for the K_f-functions. r_vec
#' overrides r_max. The bins are exclusive on the left and inclusive on
#' the right. If the first vector element has value zero, it will be
#' regarded as the collapsed bin [0, 0], and the next bin will start from
#' and exclude 0.
#' @param do_besags_L A boolean describing whether Besag's L-function should
#' also be returned where available.
#' @param use_biased_lambda2 A logical scalar on whether to use the biased
#' or the unbiased (in the Poisson case) estimate of the intensity
#' squared.
#' @param method The name of the method to create simulations under the null
#' hypothesis. 'permute' results in permutations of the marks. Using
#' 'sample' will sample the marks from the empirical mark distribution
#' with replacement. 'permute' is the default.
#' @param ... Currently unused.
#' @return A list containing the radius vector (r) and an array (a)
#' containing the summary function estimates for the original pattern and
#' the randomly labelled patterns with dimensions orig_and_perm,
#' summ_func, r. The estimates for the given marked point pattern are on
#' the row named "original".
#' @importFrom abind abind
#' @import spatstat
#' @export
#' @examples
#' require(spatstat)
#' pp <- spruces
#' res <- summ_func_random_labelling(pp)
summ_func_random_labelling <-
function(pattern, edge_correction = 'translate',
mtf_name = c('1', 'm', 'mm', 'gamma', 'gammaAbs', 'morAbs'),
nsim = 999L, r_max = NULL, r_vec = NULL, do_besags_L = TRUE,
method = 'permute', use_biased_lambda2 = FALSE, ...) {
check_pattern(pattern)
check_n_perm(nsim)
# 'translate' edge correction, only rectangle and disc windows are allowed.
# For 'none' correction, also other windows can be used.
if(edge_correction != 'none') {
if ( !( spatstat::is.rectangle(pattern[['window']]) || is.disc(pattern[['window']])) ) {
stop('The window must have type \"rectangle\" or be a disc specified as \"polygonal\".')
}
}
# Handle things related to distances of points.
radius_l <- consider_radius(pattern, r_max, r_vec)
r_vec <- radius_l[['r_vec']]
n_bin <- radius_l[['n_bin']]
nearby_arr_idx <- radius_l[['nearby_arr_idx']]
bin_idx <- radius_l[['bin_idx']]
# Handle other things related to the unmarked point pattern.
one_per_lambda2 <-
one_per_lambda_squared(pattern,
use_biased_lambda2 = use_biased_lambda2)
edge_corr <- do_edge_correction(pattern, edge_correction,
nearby_arr_idx)
orig_marks <- pattern[['marks']]
# Create mark test functions and the corresponding weights for pairs.
mtf_l <- create_mtfs_and_weights(orig_marks, mtf_name, edge_corr,
one_per_lambda2)
orig_mtf_func_l <- mtf_l[['mtf_func_l']]
orig_weight_m <- mtf_l[['weight_m']]
# Pick relevant marks.
orig_mark_l <- marks_within_radius(orig_marks, nearby_arr_idx)
orig_mark1 <- orig_mark_l[['mark1']]
orig_mark2 <- orig_mark_l[['mark2']]
# Calculate summary function estimates.
# Dims: r, mtf
orig_summ_func_m <- K_f(orig_mark1, orig_mark2, bin_idx, orig_weight_m,
orig_mtf_func_l, mtf_name, n_bin)
# Change dims into: orig_and_perm, summ_func, r
orig_summ_func_a <- t(orig_summ_func_m)
dim(orig_summ_func_a) <- c(1L, dim(orig_summ_func_a))
dimnames(orig_summ_func_a) <- c(list(orig_and_perm = c('original')),
rev(dimnames(orig_summ_func_m)))
# FIXME: K_1 gets calculated again and again. It is enough to do it for
# the original marks if it is asked for by the user at all. Then remove
# '1' from mtf_name, orig_weight_m, orig_mtf_func_l and add K_1 and
# possibly L_1 back later for the simulations.
if (nsim > 0L) {
# Simulate from the null hypothesis and estimate summary functions.
if (method %in% 'permute') {
# Result dimensions: r, summ_func, permutation.
sim_summ_func_a <-
replicate(nsim, {
perm_marks <- sample(orig_marks, replace = FALSE)
perm_mark_l <- marks_within_radius(perm_marks,
nearby_arr_idx)
perm_mark1_vec <- perm_mark_l[['mark1']]
perm_mark2_vec <- perm_mark_l[['mark2']]
perm_summ_func_m <- K_f(perm_mark1_vec,
perm_mark2_vec, bin_idx,
orig_weight_m,
orig_mtf_func_l, mtf_name,
n_bin)
})
} else if (method %in% 'sample') {
# Result dimensions: r, summ_func, permutation.
sim_summ_func_a <-
replicate(nsim, {
perm_marks <- sample(orig_marks, replace = TRUE)
perm_mtf_l <-
create_mtfs_and_weights(perm_marks,
mtf_name,
edge_corr,
one_per_lambda2)
perm_mtf_func_l <- perm_mtf_l[['mtf_func_l']]
perm_weight_m <- perm_mtf_l[['weight_m']]
perm_mark_l <- marks_within_radius(perm_marks,
nearby_arr_idx)
perm_mark1_vec <- perm_mark_l[['mark1']]
perm_mark2_vec <- perm_mark_l[['mark2']]
perm_summ_func_m <- K_f(perm_mark1_vec,
perm_mark2_vec, bin_idx,
perm_weight_m,
perm_mtf_func_l, mtf_name,
n_bin)
})
} else {
stop('method was not recognized.')
}
sim_summ_func_a <-
array(sim_summ_func_a,
dim=c(n_bin, dim(orig_summ_func_a)[2], nsim),
dimnames=list(r=NULL,
summ_func=dimnames(orig_summ_func_a)[[2]],
orig_and_perm=NULL))
sim_summ_func_a <- aperm(sim_summ_func_a, c(3L, 2L, 1L))
all_summ_func_a <- abind(orig_summ_func_a, sim_summ_func_a,
along = 1L)
names(dimnames(all_summ_func_a)) <- c('orig_and_perm', 'summ_func',
'r')
} else {
all_summ_func_a <- orig_summ_func_a
}
if (do_besags_L) {
is_any_mark_neg <- any(orig_marks < 0)
is_any_mark_pos <- any(orig_marks > 0)
L_summ_a <- all_besags_L(all_summ_func_a,
is_any_mark_neg = is_any_mark_neg,
is_any_mark_pos = is_any_mark_pos)
all_summ_func_a <- abind(all_summ_func_a, L_summ_a, along = 2L)
names(dimnames(all_summ_func_a)) <- c('orig_and_perm', 'summ_func',
'r')
}
res <- list(r = r_vec, a = all_summ_func_a, call = match.call())
class(res) <- 'all_summ_func_rl'
res
}
plot.all_summ_func_rl <- function(x, mtf_name = NULL, L = TRUE, nrow = NULL,
ncol = NULL, ...) {
# FIXME
# Make a data frame according to the arguments, i.e. choose which
# functions to plot.
# Plot using ggplot.
}
#' Checks if a window object is a circle
#'
#' @param x An object to be checked.
#' @param digits The accuracy to check whether the points at the polygonial boundary
#' are on the arch of a circle.
#' @param ... Ignored.
#' @importFrom spatstat is.owin
is.disc <- function(x, digits=6, ...) {
if(!is.owin(x)) stop("x is not an owin object.")
if(x[['type']] != 'polygonal') return(FALSE)
if(diff(x[['xrange']]) != diff(x[['yrange']])) return(FALSE)
r <- diff(x[['xrange']])/2
x0 <- mean(x[['xrange']])
y0 <- mean(x[['yrange']])
if ( sum( round((x[['bdry']][[1]][['x']]-x0)^2 + (x[['bdry']][[1]][['y']]-y0)^2, digits=digits) == round(r^2, digits=digits) ) != length(x[['bdry']][[1]][['x']]) ) return(FALSE)
return(TRUE)
}
#' Finds the radius and centre of a circle defined by a polygonial boundary
#'
#' @param x An owin object.
#' @param ... Ignored.
disc_param <- function(x, ...) {
if(!is.disc(x)) stop("x is not a disc.")
R <- diff(x[['xrange']])/2
x0 <- mean(x[['xrange']])
y0 <- mean(x[['yrange']])
list(R = R, centre = c(x0, y0))
}
#' Checks that the given pattern is valid.
#'
#' @param pattern A pattern to be checked.
#'
#' @importFrom spatstat is.ppp
#' @importFrom spatstat is.empty.ppp
#' @importFrom spatstat is.polygonal
#' @importFrom spatstat is.rectangle
check_pattern <- function(pattern) {
if (length(pattern) < 1L) {
stop('Pattern was not given.')
}
if (!spatstat::is.ppp(pattern)) {
stop('The pattern must be a ppp object.')
}
if (spatstat::is.empty.ppp(pattern)) {
stop('The pattern must not be empty.')
}
marks <- pattern[['marks']]
if (length(marks) < 1L) {
stop('The pattern must have marks.')
}
if (pattern[['markformat']] != 'vector') {
stop('The markformat must be vector.')
}
if (!is.numeric(marks)) {
stop('The marks must be of type numeric or integer.')
}
if ( !(spatstat::is.polygonal(pattern[['window']]) || spatstat::is.rectangle(pattern[['window']])) ) {
stop('The window must be specified as \"polygonal\" or \"rectangle\".')
}
# if ( !( is.rectangle(pattern[['window']]) || is.disc(pattern[['window']])) ) {
# stop('The window must have type \"rectangle\" or be a disc specified as \"polygonal\".')
# }
if (any(duplicated(matrix(c(pattern[['x']], pattern[['y']]), ncol = 2L),
MARGIN = 1L))) {
stop('The pattern must be simple i.e. without duplicate points ',
'location-wise.')
}
}
#' Check that the given number of permutations makes sense.
#'
#' @param nsim The value/object to be checked.
check_n_perm <- function(nsim) {
if (length(nsim) != 1L || !is.finite(nsim) || nsim < 0L) {
stop('Number of permutations has to be a scalar, finite, ',
'non-negative number.')
}
}
#' Estimate 1 / lambda ^ 2.
#'
#' Notice that this estimator is biased.
#'
#' @param pattern A \code{\link[spatstat]{ppp}} object.
#' @param use_biased_lambda2 A logical. Use the unbiased estimator or not.
#' @import spatstat
one_per_lambda_squared <- function(pattern, use_biased_lambda2) {
if (length(use_biased_lambda2) != 1L ||
!is.logical(use_biased_lambda2) ||
!is.finite(use_biased_lambda2)) {
stop('use_biased_lambda2 must be either TRUE or FALSE.')
}
area <- spatstat::area.owin(pattern[['window']])
n_point <- pattern[['n']]
if (use_biased_lambda2) {
one_per_lambda2 <- area * area / (n_point * n_point)
} else {
one_per_lambda2 <- area * area / (n_point * (n_point - 1L))
}
one_per_lambda2
}
#' Besag's L-transform.
#'
#' @param k A vector of values to be transformed.
besags_L <- function(k) {
sqrt(1 / pi * k)
}
#' Channel to Besag's L.
#'
#' Pick only those K_f that have a valid Besag's L-transform and return the
#' transformed summary functions.
#'
#' @param all_k_a An array containing all K_f for all patterns. Dimensions:
#' orig_and_perm, summ_func, r.
#' @param is_any_mark_neg A logical.
#' @param is_any_mark_pos A logical.
#' @param sep_char A character used to separate L_char and mtf_name.
#' @param L_char A character used for the tranformed function.
#' @return An array of dimensions: orig_and_perm, summ_func, r. The size of
#' the second dimension depends on how many K_f functions had a valid
#' L-transform.
all_besags_L <- function(all_k_a, is_any_mark_neg, is_any_mark_pos,
sep_char = '_', L_char = 'L') {
K_f_name <- dimnames(all_k_a)[['summ_func']]
# First row of sapply result is all 'K', second row mtf.
mtf_name <- sapply(strsplit(K_f_name, sep_char, fixed=TRUE),
identity)[2, , drop = TRUE]
valid_mtf_name <- besags_L_valid(mtf_name,
is_any_mark_neg = is_any_mark_neg,
is_any_mark_pos = is_any_mark_pos)
match.idx <- match(valid_mtf_name, mtf_name)
if (length(match.idx) < 1L) {
all_l_a <- NULL
} else {
all_l_a <- besags_L(all_k_a[, match.idx, , drop = FALSE])
L_f_name <- paste(L_char, mtf_name, sep = sep_char)
dimnames(all_l_a)[2] <- list(summ_func = L_f_name)
}
all_l_a
}
#' Calculates the cumulated K_f-estimates.
#'
#' @param ... Arguments to be passed to rho_f.
K_f <- function(...) {
rho_f_m <- rho_f(...)
# as.matrix is needed if rho_f_m has only one r-value.
n_r <- nrow(rho_f_m)
K_f_m <- as.matrix(apply(rho_f_m, 2, cumsum), nrow = n_r)
if (n_r == 1L) {
# apply flips the matrix over if there is only one r-value.
K_f_m <- t(K_f_m)
}
colnames(K_f_m) <- paste('K_', colnames(rho_f_m), sep = '')
names(dimnames(K_f_m)) <- c('r', 'summ_func')
K_f_m
}
#' Weigh and sum into bins.
#'
#' Weigh the value of the f-function for each pair, sum together results
#' of the same bin and place the sum in the correct bin.
#'
#' @param f_value A vector of values from the mark test function
#' corresponding to each pair of points.
#' @param weight A vector of weights for each pair of points.
#' @param idx_vec A vector of bin indices to guide each pair of points into
#' the correct bin.
#' @param n_bin The amount of bins that the returned vector should have.
#' @return A rho_f vector with length of n_bin.
#' @useDynLib marksummary
#' @import Rcpp
weigh_and_bin <- function(f_value, weight, idx_vec, n_bin) {
.Call('weighAndBin', f_value, weight, idx_vec, n_bin,
PACKAGE = 'marksummary')
}
#' Calculates the uncumulated rho_f-estimates
#'
#' @param mark1 The marks of the first points of pairs in a vector.
#' @param mark2 The marks of the second points of pairs in a vector.
#' @param bin_idx A vector of indices which tell which radius bin each point
#' pair belongs to.
#' @param weight_m A weight matrix for point pairs.
#' @param mtf_func_l A list of mark test functions all of which take two
#' real marks and return a real value.
#' @param mtf_name The names of the mark test functions. Matches mtf_func_l
#' in order.
#' @param n_bin The number of bins in the radius vector.
#' @return A matrix of rho_f values. Dimensions r, mtf.
rho_f <- function(mark1, mark2, bin_idx, weight_m, mtf_func_l, mtf_name,
n_bin) {
val_m <- as.matrix(vapply(seq_along(mtf_func_l),
function(mtf_idx) {
mtf_func <- mtf_func_l[[mtf_idx]]
f_value <- mtf_func(mark1, mark2)
one_rho_f <-
weigh_and_bin(f_value,
weight_m[, mtf_idx, drop = TRUE],
bin_idx, n_bin)
one_rho_f
}, numeric(n_bin)))
# vapply acts differently depending on whether the result is a vector or
# a scalar.
if (n_bin == 1L) {
val_m <- t(val_m)
}
colnames(val_m) <- mtf_name
val_m
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.