From fe2ba5186d587fb8ac5ef75a33c2989e0de5f5c4 Mon Sep 17 00:00:00 2001 From: Adrian Baddeley Date: Sun, 6 Oct 2024 15:53:28 +0800 Subject: [PATCH] Diffusion smoothing --- DESCRIPTION | 4 +- NAMESPACE | 12 ++++ NEWS | 20 +++++- R/SmoothHeat.R | 78 +++++++++++++++++++++ R/blurHeat.R | 137 +++++++++++++++++++++++++++++++++++++ R/densityHeat.ppp.R | 4 ++ R/relriskHeat.R | 108 +++++++++++++++++++++++++++++ inst/doc/packagesizes.txt | 2 +- inst/info/packagesizes.txt | 2 +- man/SmoothHeat.Rd | 39 +++++++++++ man/SmoothHeat.ppp.Rd | 58 ++++++++++++++++ man/blurHeat.Rd | 83 ++++++++++++++++++++++ man/bw.relriskHeatppp.Rd | 99 +++++++++++++++++++++++++++ man/relriskHeat.Rd | 103 ++++++++++++++++++++++++++++ 14 files changed, 744 insertions(+), 5 deletions(-) create mode 100644 R/SmoothHeat.R create mode 100644 R/blurHeat.R create mode 100644 R/relriskHeat.R create mode 100644 man/SmoothHeat.Rd create mode 100644 man/SmoothHeat.ppp.Rd create mode 100644 man/blurHeat.Rd create mode 100644 man/bw.relriskHeatppp.Rd create mode 100644 man/relriskHeat.Rd diff --git a/DESCRIPTION b/DESCRIPTION index f1edcc2..7716fdd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: spatstat.explore -Version: 3.3-2.001 -Date: 2024-09-06 +Version: 3.3-2.002 +Date: 2024-10-06 Title: Exploratory Data Analysis for the 'spatstat' Family Authors@R: c(person("Adrian", "Baddeley", role = c("aut", "cre", "cph"), diff --git a/NAMESPACE b/NAMESPACE index 87431d3..8b97f7f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -66,6 +66,8 @@ export("bind.ratfv") export("bits.envelope") export("bits.test") export("blur") +export("blurHeat") +export("blurHeat.im") export("boyce") export("bw.abram.ppp") export("bw.CvL") @@ -78,6 +80,7 @@ export("bw.pcf") export("bw.ppl") export("bw.pplHeat") export("bw.relrisk") +export("bw.relriskHeatppp") export("bw.relrisk.ppp") export("bw.scott") export("bw.scott.iso") @@ -424,6 +427,8 @@ export("rectcontact") export("RelevantDeviation") export("reload.or.compute") export("relrisk") +export("relriskHeat") +export("relriskHeat.ppp") export("relrisk.ppp") export("rename.fv") export("resolve.2D.kernel") @@ -482,6 +487,9 @@ export("smoothcrossEngine") export("Smoothfun") export("Smoothfun.ppp") export("Smooth.fv") +export("SmoothHeat") +export("SmoothHeat.im") +export("SmoothHeat.ppp") export("Smooth.im") export("smoothpointsEngine") export("Smooth.ppp") @@ -568,6 +576,7 @@ S3method("as.data.frame", "fv") S3method("as.tess", "quadrattest") S3method("auc", "ppp") S3method("berman.test", "ppp") +S3method("blurHeat", "im") S3method("bw.abram", "ppp") S3method("bw.relrisk", "ppp") S3method("cbind", "fv") @@ -671,6 +680,7 @@ S3method("quadrat.test", "quadratcount") S3method("quadrat.test", "splitppp") S3method("range", "ssf") S3method("[", "rat") +S3method("relriskHeat", "ppp") S3method("relrisk", "ppp") S3method("resolve.lambdacross", "ppp") S3method("resolve.lambda", "ppp") @@ -691,6 +701,8 @@ S3method("shift", "quadrattest") S3method("simulate", "rhohat") S3method("Smoothfun", "ppp") S3method("Smooth", "fv") +S3method("SmoothHeat", "im") +S3method("SmoothHeat", "ppp") S3method("Smooth", "im") S3method("Smooth", "ppp") S3method("Smooth", "solist") diff --git a/NEWS b/NEWS index a59f1bb..005202e 100644 --- a/NEWS +++ b/NEWS @@ -1,10 +1,28 @@ - CHANGES IN spatstat.explore VERSION 3.3-2.001 + CHANGES IN spatstat.explore VERSION 3.3-2.002 OVERVIEW + o relative risk estimation using diffusion. + + o smoothing using diffusion. + o Tweaks to bandwidth selection. +NEW FUNCTIONS + + o relriskHeat, relriskHeat.ppp + Relative risk estimation using diffusion. + + o blurHeat, blurHeat.im + Image smoothing using diffusion. + + o SmoothHeat, SmoothHeat.ppp + Smoothing numerical values observed at points, using diffusion. + + o bw.relriskHeatppp + Bandwidth selection for relriskHeat.ppp + SIGNIFICANT USER-VISIBLE CHANGES o bw.ppl diff --git a/R/SmoothHeat.R b/R/SmoothHeat.R new file mode 100644 index 0000000..cf58cf2 --- /dev/null +++ b/R/SmoothHeat.R @@ -0,0 +1,78 @@ +#' +#' SmoothHeat.R +#' +#' Nadaraya-Watson style smooth regression using diffusion +#' +#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit +#' +#' $Revision: 1.3 $ $Date: 2024/10/06 01:26:29 $ + +SmoothHeat <- function(X, ...) { + UseMethod("SmoothHeat") +} + +SmoothHeat.im <- function(X, sigma, ...) { + blurHeat(X, sigma, ...) +} + +SmoothHeat.ppp <- function(X, sigma, ..., weights=NULL) { + stopifnot(is.ppp(X)) + stopifnot(is.marked(X)) + marx <- marks(X) + if(!is.vector(marx)) stop("Marks of X should be a numeric vector") + marx <- as.numeric(marx) + if(is.null(weights)) { + numwt <- marx + denwt <- NULL + } else { + check.nvector(weights, npoints(X), oneok=TRUE) + if(length(weights) == 1) weights <- rep(weights, npoints(X)) + numwt <- marx * weights + denwt <- weights + } + Y <- unmark(X) + numer <- densityHeat(Y, sigma, weights=numwt, ...) + denom <- densityHeat(Y, sigma, weights=denwt, ...) + return(numer/denom) +} + +bw.SmoothHeatppp <- function(X, ..., weights=NULL, + srange=NULL, ns=16, sigma=NULL, + leaveoneout=TRUE, verbose=TRUE) { + stopifnot(is.ppp(X)) + stopifnot(is.marked(X)) + marx <- marks(X) + if(!is.vector(marx)) stop("Marks of X should be a numeric vector") + marx <- as.numeric(marx) + if(is.null(weights)) { + numwt <- marx + denwt <- NULL + } else { + check.nvector(weights, npoints(X), oneok=TRUE) + if(length(weights) == 1) weights <- rep(weights, npoints(X)) + numwt <- marx * weights + denwt <- weights + } + #' compute weighted and unweighted intensity estimates + U <- unmark(X) + aNumer <- HeatEstimates.ppp(U, ..., weights=numwt, + srange=srange, ns=ns, sigma=sigma, + leaveoneout=leaveoneout, verbose=verbose) + aDenom <- HeatEstimates.ppp(U, ..., weights=denwt, + srange=srange, ns=ns, sigma=sigma, + leaveoneout=leaveoneout, verbose=verbose) + h <- aDenom$h + hname <- aDenom$hname + #' compute smoother + zhat <- aNumer$lambda/aDenom$lambda + #' compute least squares cross-validation criterion + zobs <- matrix(marx, nrow(zhat), ncol(zhat), byrow=TRUE) + CV <- rowSums((zhat - zobs)^2) + iopt <- which.min(CV) + result <- bw.optim(CV, h, iopt, + criterion="Least squares cross-validation", + hname=hname, + unitname=unitname(X)) + return(result) +} + diff --git a/R/blurHeat.R b/R/blurHeat.R new file mode 100644 index 0000000..00fe83b --- /dev/null +++ b/R/blurHeat.R @@ -0,0 +1,137 @@ +#' +#' blurHeat.R +#' +#' Image blurring by diffusion +#' +#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit +#' +#' Licence: GNU Public Licence >= 2 +#' +#' $Revision: 1.3 $ $Date: 2024/10/06 02:28:55 $ + +blurHeat <- function(X, ...) { + UseMethod("blurHeat") +} + +blurHeat.im <- function(X, sigma, ..., connect=8, + symmetric=FALSE, k=1, show=FALSE) { + Y <- as.im(X) + check.1.integer(k) + stopifnot(k >= 1) + if(!(connect %in% c(4,8))) + stop("connectivity must be 4 or 8") + if(is.im(sigma)) { + # ensure Y and sigma are on the same grid + A <- harmonise(Y=Y, sigma=sigma) + Y <- A$Y + sigma <- A$sigma + } else if(is.function(sigma)) { + sigma <- as.im(sigma, as.owin(Y)) + } else check.1.real(sigma) + #' initial state + v <- as.matrix(Y) + u <- as.vector(v) + #' symmetric random walk? + if(symmetric) { + asprat <- with(Y, ystep/xstep) + if(abs(asprat-1) > 0.01) + warning(paste("Symmetric random walk on a non-square grid", + paren(paste("aspect ratio", asprat))), + call.=FALSE) + } + #' determine appropriate jump probabilities & time step + pmax <- 1/(connect+1) # maximum permitted jump probability + xstep <- Y$xstep + ystep <- Y$ystep + minstep <- min(xstep, ystep) + if(symmetric) { + #' all permissible transitions have the same probability 'pjump'. + #' Determine Nstep, and dt=sigma^2/Nstep, such that + #' Nstep >= 16 and M * pjump * minstep^2 = dt + M <- if(connect == 4) 2 else 6 + Nstep <- max(16, ceiling(max(sigma)^2/(M * pmax * minstep^2))) + sn <- (sigma^2)/Nstep + px <- py <- pxy <- sn/(M * minstep^2) + } else { + #' px is the probability of jumping 1 step to the right + #' py is the probability of jumping 1 step up + #' if connect=4, horizontal and vertical jumps are exclusive. + #' if connect=8, horizontal and vertical increments are independent + #' Determine Nstep, and dt = sigma^2/Nstep, such that + #' Nstep >= 16 and 2 * pmax * minstep^2 = dt + Nstep <- max(16, ceiling(max(sigma)^2/(2 * pmax * minstep^2))) + sn <- (sigma^2)/Nstep + px <- sn/(2 * xstep^2) + py <- sn/(2 * ystep^2) + if(max(px) > pmax) stop("Internal error: px exceeds pmax") + if(max(py) > pmax) stop("Internal error: py exceeds pmax") + if(connect == 8) pxy <- px * py + } + #' construct adjacency matrices + dimv <- dim(v) + my <- gridadjacencymatrix(dimv, across=FALSE, down=TRUE, diagonal=FALSE) + mx <- gridadjacencymatrix(dimv, across=TRUE, down=FALSE, diagonal=FALSE) + if(connect == 8) + mxy <- gridadjacencymatrix(dimv, across=FALSE, down=FALSE, diagonal=TRUE) + #' restrict to window + if(anyNA(u)) { + ok <- !is.na(u) + u <- u[ok] + mx <- mx[ok,ok,drop=FALSE] + my <- my[ok,ok,drop=FALSE] + if(connect == 8) + mxy <- mxy[ok,ok,drop=FALSE] + if(is.im(sigma)) { + px <- px[ok] + py <- py[ok] + if(connect == 8) + pxy <- pxy[ok] + } + } else ok <- TRUE + #' construct iteration matrix + if(connect == 4) { + A <- px * mx + py * my + } else { + A <- px * (1 - 2 * py) * mx + py * (1 - 2 * px) * my + pxy * mxy + } + #' debug + stopifnot(min(rowSums(A)) >= 0) + stopifnot(max(rowSums(A)) <= 1) + #' + diag(A) <- 1 - rowSums(A) + #' k-step transition probabilities + if(k > 1) { + Ak <- A + for(j in 2:k) Ak <- Ak %*% A + } else Ak <- A + k <- as.integer(k) + Nstep <- as.integer(Nstep) + Nblock <- Nstep/k + Nrump <- Nstep - Nblock * k + #' run + U <- u + Z <- Y + if(!show) { + for(istep in 1:Nblock) U <- U %*% Ak + } else { + opa <- par(ask=FALSE) + each <- max(1, round(Nblock/60)) + for(istep in 1:Nblock) { + U <- U %*% Ak + if(istep %% each == 0) { + Z[] <- as.vector(U) + f <- sqrt(istep/Nstep) + main <- if(is.im(sigma)) paste(signif(f, 3), "* sigma") else + paste("sigma =", signif(f * sigma, 3)) + plot(Z, main=main) + Sys.sleep(0.4) + } + } + par(opa) + } + if(Nrump > 0) for(istep in 1:Nrump) U <- U %*% A + #' pack up + Z[] <- as.vector(U) + return(Z) +} + diff --git a/R/densityHeat.ppp.R b/R/densityHeat.ppp.R index a64a080..7083ddb 100644 --- a/R/densityHeat.ppp.R +++ b/R/densityHeat.ppp.R @@ -3,6 +3,10 @@ #' #' Diffusion estimator of density/intensity #' +#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit +#' +#' Licence: GNU Public Licence >= 2 +#' densityHeat <- function(x, sigma, ...) { UseMethod("densityHeat") diff --git a/R/relriskHeat.R b/R/relriskHeat.R new file mode 100644 index 0000000..087174c --- /dev/null +++ b/R/relriskHeat.R @@ -0,0 +1,108 @@ +#' +#' relriskHeat.R +#' +#' Relative risk/conditional probability using diffusion smoothing +#' +#' Copyright (C) 2018-2024 Adrian Baddeley, Tilman Davies and Suman Rakshit +#' +#' $Revision: 1.7 $ $Date: 2024/10/06 02:55:11 $ +#' +#' + +relriskHeat <- function(X,...) { + UseMethod("relriskHeat") +} + +relriskHeat.ppp <- function(X,..., sigmaX=NULL, weights=NULL){ + stopifnot(is.ppp(X)) + stopifnot(is.multitype(X)) + nX <- npoints(X) + + Y <- split(X) + ntypes <- length(Y) + if(ntypes == 1) + stop("Data contains only one type of points") + + type <- marks(X) + + if(length(sigmaX)) { + check.nvector(sigmaX, nX) + sigmaX <- split(sigmaX, type) + } else sigmaX <- rep(list(NULL), ntypes) + + if(length(weights)) { + check.nvector(weights, nX) + weights <- split(weights, type) + } else weights <- rep(list(NULL), ntypes) + + Deach <- mapply(densityHeat, x=Y, sigmaX=sigmaX, weights=weights, + MoreArgs=list(...), SIMPLIFY=FALSE) + + Dall <- Reduce("+", Deach) + + probs <- solapply(Deach, "/", e2=Dall) + + return(probs) +} + +bw.relriskHeatppp <- function(X, ..., + method=c("likelihood", "leastsquares"), + weights=NULL, + srange=NULL, ns=16, sigma=NULL, + leaveoneout=TRUE, verbose=TRUE) { + stopifnot(is.ppp(X)) + stopifnot(is.multitype(X)) + method <- match.arg(method) + U <- unmark(X) + Y <- split(X) + Denominator <- HeatEstimates.ppp(U, ..., weights=weights, + srange=srange, ns=ns, sigma=sigma, + leaveoneout=leaveoneout, verbose=verbose) + h <- Denominator$h + hname <- Denominator$hname + # extract denominator value for each sigma (row) and each data point (col) + lambda.denom <- Denominator$lambda + #' + if(is.null(weights)) { + Numerators <- lapply(Y, HeatEstimates.ppp, ..., + srange=srange, ns=ns, sigma=sigma, + leaveoneout=leaveoneout, verbose=verbose) + } else { + check.nvector(weights, npoints(X), oneok=TRUE) + if(length(weights) == 1) weights <- rep(weights, npoints(X)) + wsplit <- split(weights, marks(X)) + Numerators <- mapply(HeatEstimates.ppp, + X=Y, weights=wsplit, + MoreArgs = list(..., + srange=srange, ns=ns, sigma=sigma, + leaveoneout=leaveoneout, + verbose=verbose), + SIMPLIFY=FALSE) + } + #' extract estimates of numerator terms + lamlist <- lapply(Numerators, getElement, name="lambda") + #' reassemble into original position + #' (tried to do this with 'unsplit' but it's too messy) + opos <- split(seq_len(npoints(X)), marks(X)) + lambda.numer <- matrix(, nrow=nrow(lambda.denom), ncol=ncol(lambda.denom)) + for(k in seq_along(opos)) { + if(length(opos.k <- opos[[k]])) + lambda.numer[ , opos.k] <- lamlist[[k]] + } + #' compute predicted probability of observations + phat <- lambda.numer/lambda.denom + #' compute cross-validation criterion + switch(method, + likelihood = { + CV <- -rowMeans(log(phat)) + cname <- "Likelihood cross-validation" + }, + leastsquares = { + CV <- rowMeans((1 - phat)^2) + cname <- "Least squares cross-validation" + }) + result <- bw.optim(CV, h, criterion=cname, hname=hname, + unitname=unitname(X)) + return(result) +} + diff --git a/inst/doc/packagesizes.txt b/inst/doc/packagesizes.txt index fdf0e3c..2d50e05 100755 --- a/inst/doc/packagesizes.txt +++ b/inst/doc/packagesizes.txt @@ -11,4 +11,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines "2024-03-21" "3.2-7" 241 535 0 32149 6365 "2024-07-09" "3.3-1" 230 517 0 31571 6365 "2024-08-20" "3.3-2" 230 517 0 31571 6364 -"2024-09-06" "3.3-2.001" 230 517 0 31576 6364 +"2024-10-06" "3.3-2.002" 235 525 0 31903 6364 diff --git a/inst/info/packagesizes.txt b/inst/info/packagesizes.txt index fdf0e3c..2d50e05 100755 --- a/inst/info/packagesizes.txt +++ b/inst/info/packagesizes.txt @@ -11,4 +11,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines "2024-03-21" "3.2-7" 241 535 0 32149 6365 "2024-07-09" "3.3-1" 230 517 0 31571 6365 "2024-08-20" "3.3-2" 230 517 0 31571 6364 -"2024-09-06" "3.3-2.001" 230 517 0 31576 6364 +"2024-10-06" "3.3-2.002" 235 525 0 31903 6364 diff --git a/man/SmoothHeat.Rd b/man/SmoothHeat.Rd new file mode 100644 index 0000000..cfaf788 --- /dev/null +++ b/man/SmoothHeat.Rd @@ -0,0 +1,39 @@ +\name{SmoothHeat} +\alias{SmoothHeat} +\title{Spatial Smoothing of Data by Diffusion} +\description{ + Generic function to perform spatial smoothing of spatial data + by diffusion. +} +\usage{ + SmoothHeat(X, \dots) +} +\arguments{ + \item{X}{Some kind of spatial data} + \item{\dots}{Arguments passed to methods.} +} +\details{ + This generic function calls an appropriate method + to perform spatial smoothing on the spatial dataset \code{X} + using diffusion. + + Methods for this function include + \itemize{ + \item \code{\link[spatstat.explore]{SmoothHeat.ppp}} for point patterns + \item \code{\link[spatstat.explore]{SmoothHeat.im}} for pixel images. + } +} +\seealso{ + \code{\link[spatstat.explore]{SmoothHeat.ppp}}, + \code{\link[spatstat.explore]{SmoothHeat.im}}. +} +\value{ + An object containing smoothed values of the input data, + in an appropriate format. See the documentation for the methods. +} +\author{ + \adrian. +} +\keyword{spatial} +\keyword{methods} +\keyword{smooth} diff --git a/man/SmoothHeat.ppp.Rd b/man/SmoothHeat.ppp.Rd new file mode 100644 index 0000000..7c9b0cf --- /dev/null +++ b/man/SmoothHeat.ppp.Rd @@ -0,0 +1,58 @@ +\name{SmoothHeat.ppp} +\alias{SmoothHeat.ppp} +\title{Spatial Smoothing of Observations using Diffusion Estimate of Density} +\description{ + Performs spatial smoothing of numeric values observed + at a set of irregular locations, using the diffusion estimate + of the density. +} +\usage{ +\method{SmoothHeat}{ppp}(X, sigma, \dots, weights=NULL) +} +\arguments{ + \item{X}{ + Point pattern (object of class \code{"ppp"}) + with numeric marks. + } + \item{sigma}{ + Smoothing bandwidth. A single number giving the equivalent + standard deviation of the smoother. + } + \item{\dots}{ + Arguments passed to \code{\link[spatstat.explore]{densityHeat}} + controlling the estimation of each marginal intensity, + or passed to \code{\link[spatstat.geom]{pixellate.ppp}} + controlling the pixel resolution. + } + \item{weights}{Optional numeric vector of weights associated with each + data point. + } +} +\details{ + This is the analogue of the Nadaraya-Watson smoother, using the + diffusion smoothing estimation procedure (Baddeley et al, 2022). + The numerator and denominator of the Nadaraya-Watson smoother are + calculated using \code{\link[spatstat.explore]{densityHeat.ppp}}. +} +\value{ + Pixel image (object of class \code{"im"}) giving the smoothed + mark value. +} +\seealso{ + \code{\link[spatstat.explore]{Smooth.ppp}} for the usual kernel-based + smoother (the Nadaraya-Watson smoother) + and \code{\link[spatstat.explore]{densityHeat}} for the diffusion estimate of density. +} +\author{ + \adrian, \tilman and Suman Rakshit. +} +\examples{ + plot(SmoothHeat(longleaf, 10)) +} +\references{ + Baddeley, A., Davies, T., Rakshit, S., Nair, G. and McSwiggan, G. (2022) + Diffusion smoothing for spatial point patterns. + \emph{Statistical Science} \bold{37}, 123--142. +} +\keyword{spatial} +\keyword{smooth} diff --git a/man/blurHeat.Rd b/man/blurHeat.Rd new file mode 100644 index 0000000..777b846 --- /dev/null +++ b/man/blurHeat.Rd @@ -0,0 +1,83 @@ +\name{blurHeat} +\alias{blurHeat} +\alias{blurHeat.im} +\alias{SmoothHeat.im} +\title{ + Diffusion Blur +} +\description{ + Blur a Pixel Image by Applying Diffusion +} +\usage{ +blurHeat(X, \dots) + +\method{blurHeat}{im}(X, sigma, \dots, + connect = 8, symmetric = FALSE, k= 1, show = FALSE) + +\method{SmoothHeat}{im}(X, sigma, \dots) + +} +\arguments{ + \item{X}{ + Pixel image (object of class \code{"im"}). + } + \item{sigma}{ + Smoothing bandwidth. A numeric value, a pixel image or + a \code{function(x,y)}. + } + \item{\dots}{ + Ignored by \code{blurHeat.im}. + } + \item{connect}{ + Grid connectivity: either 4 or 8. + } + \item{symmetric}{ + Logical value indicating whether to \emph{force} the algorithm + to use a symmetric random walk. + } + \item{k}{ + Integer. Calculations will be performed by repeatedly multiplying + the current state by the \code{k}-step transition matrix. + } + \item{show}{ + Logical value indicating whether to plot successive iterations. + } +} +\details{ + The function \code{blurHeat} is generic. + + This help file documents the method \code{blurHeat.im} for pixel images + (objects of class \code{"im"}). This is currently equivalent + to \code{SmoothHeat.im}, which is also documented here. + + If \code{sigma} is a numeric value, then + the classical time-dependent heat equation is solved + up to time \code{t = sigma^2} starting with the initial + condition given by the image \code{X}. This has the effect + of blurring the input image \code{X}. + + If \code{sigma} is a function or a pixel image, then + it is treated as a spatially-variable diffusion rate, + and the corresponding heat equation is solved. + + This command can be used to calculate the expected value + of the diffusion estimator of intensity (\code{\link[spatstat.explore]{densityHeat}}) + when the true intensity is known. +} +\value{ + A pixel image on the same raster as \code{X}. +} +\author{ + \adrian. +} +\seealso{ + \code{\link[spatstat.explore]{densityHeat}}, + \code{\link[spatstat.explore]{blur}}. +} +\examples{ + Z <- as.im(function(x,y) { sin(10*x) + sin(9*y) }, letterR) + ZZ <- blurHeat(Z, 0.2) + plot(solist(original=Z, blurred=ZZ), main="") +} +\keyword{spatial} +\keyword{math} diff --git a/man/bw.relriskHeatppp.Rd b/man/bw.relriskHeatppp.Rd new file mode 100644 index 0000000..3be49f2 --- /dev/null +++ b/man/bw.relriskHeatppp.Rd @@ -0,0 +1,99 @@ +\name{bw.relriskHeatppp} +\alias{bw.relriskHeatppp} +\title{ + Bandwidth Selection for Relative Risk using Diffusion +} +\description{ + Performs data-based bandwidth selection for + the diffusion estimate of relative risk \code{\link[spatstat.explore]{relriskHeat.ppp}} + using either likelihood cross-validation or least squares +} +\usage{ +bw.relriskHeatppp(X, \dots, method = c("likelihood", "leastsquares"), + weights = NULL, srange = NULL, ns = 16, sigma = NULL, + leaveoneout = TRUE, verbose = TRUE) +} +\arguments{ + \item{X}{ + A multitype point pattern (object of class \code{"ppp"}). + } + \item{\dots}{ + Arguments passed to \code{\link[spatstat.explore]{relriskHeat.ppp}}. + } + \item{method}{ + Character string specifying the cross-validation method. + Partially matched to \code{"likelihood"} for binary likelihood + cross-validation or \code{"leastsquares"} for least squares + cross-validation. + } + \item{weights}{ + Optional numeric vector of weights associated with each point of \code{X}. + } + \item{srange}{ + Numeric vector of length 2 specifying a range of bandwidths to be + considered. + } + \item{ns}{ + Integer. Number of candidate bandwidths to be considered. + } + \item{sigma}{ + Maximum smoothing bandwidth. + A numeric value, or a pixel image, or a \code{function(x,y)}. + Alternatively a numeric vector containing a sequence of + candidate bandwidths. + } + \item{leaveoneout}{ + Logical value specifying whether intensity values at data points + should be estimated using the leave-one-out rule. + } + \item{verbose}{ + Logical value specifying whether to print progress reports. + } +} +\details{ + This algorithm selects the optimal global bandwidth for + kernel estimation of relative risk for the dataset \code{X} + using diffusion smoothing \code{\link[spatstat.explore]{relriskHeat}}. + + If \code{sigma} is a numeric value, the algorithm finds the + optimal bandwidth \code{tau <= sigma}. + + If \code{sigma} is a pixel image or function, the algorithm + finds the optimal fraction \code{0 < f <= 1} such that + smoothing with \code{f * sigma} would be optimal. +} +\value{ + A numerical value giving the selected bandwidth + (if \code{sigma} was a numeric value) + or the selected fraction of the maximum bandwidth + (if \code{sigma} was a pixel image or function). + The result also belongs to the class \code{"bw.optim"} which can be + plotted. +} +\author{ + \adrian, \tilman and Suman Rakshit. +} +\seealso{ + \code{\link[spatstat.explore]{relriskHeat.ppp}} +} +\examples{ + ## bovine tuberculosis data + X <- subset(btb, select=spoligotype) + if(interactive()) { + smax <- 40 + ns <- 16 + dimyx <- NULL + } else { + ## reduce data and resolution to speed up + X <- X[c(TRUE, rep(FALSE, 7))] + smax <- 9 + ns <- 8 + dimyx <- 32 + } + b <- bw.relriskHeatppp(X, sigma=smax, ns=ns, dimyx=dimyx) + b + plot(b) +} +\keyword{spatial} +\keyword{nonparametric} + diff --git a/man/relriskHeat.Rd b/man/relriskHeat.Rd new file mode 100644 index 0000000..d2847fb --- /dev/null +++ b/man/relriskHeat.Rd @@ -0,0 +1,103 @@ +\name{relriskHeat} +\alias{relriskHeat} +\alias{relriskHeat.ppp} +\title{ + Diffusion Estimate of Conditional Probabilities +} +\description{ + Computes the conditional probability estimator of relative risk + based on a multitype point pattern using the diffusion estimate + of the type-specific intensities. +} +\usage{ +relriskHeat(X, \dots) + +\method{relriskHeat}{ppp}(X, \dots, sigmaX=NULL, weights=NULL) +} +\arguments{ + \item{X}{ + A multitype point pattern (object of class \code{"ppp"}). + } + \item{\dots}{ + Arguments passed to \code{\link[spatstat.explore]{densityHeat}} + controlling the estimation of each marginal intensity. + } + \item{sigmaX}{ + Optional. + Numeric vector of bandwidths, one associated with each data point in + \code{X}. + } + \item{weights}{ + Optional numeric vector of weights associated with each point of + \code{X}. + } +} +\details{ + The function \code{relriskHeat} is generic. This file documents the + method \code{relriskHeat.ppp} for spatial point patterns (objects of + class \code{"ppp"}). + + This function estimates the spatially-varying conditional probability + that a random point (given that it is present) will belong to + a given type. + + The algorithm separates \code{X} into + the sub-patterns consisting of points of each type. + It then applies \code{\link[spatstat.explore]{densityHeat}} to each sub-pattern, + using the same bandwidth and smoothing regimen for each sub-pattern, + as specified by the arguments \code{\dots}. + + If \code{weights} is specified, it should be a numeric vector + of length equal to the number of points in \code{X}, so that + \code{weights[i]} is the weight for data point \code{X[i]}. + + Similarly when performing lagged-arrival smoothing, + the argument \code{sigmaX} must be a numeric vector of the same length + as the number of points in \code{X}, and thus contain the + point-specific bandwidths in the order corresponding to each of these + points regardless of mark. +} +\value{ + A named list (of class \code{\link[spatstat.geom]{solist}}) + containing pixel \code{\link[spatstat.geom]{im}}ages, + giving the estimated conditional probability surfaces for each type. +} +\seealso{ + \code{\link[spatstat.explore]{relrisk.ppp}} for the + traditional convolution-based kernel estimator of + conditional probability surfaces, + and the function \code{risk} in the \pkg{sparr} package for the + density-ratio-based estimator. +} +\references{ + Agarwal, N. and Aluru, N.R. (2010) + A data-driven stochastic collocation approach for + uncertainty quantification in MEMS. + \emph{International Journal for Numerical Methods in Engineering} + \bold{83}, 575--597. + + Baddeley, A., Davies, T., Rakshit, S., Nair, G. and McSwiggan, G. (2022) + Diffusion smoothing for spatial point patterns. + \emph{Statistical Science} \bold{37}, 123--142. + + Barry, R.P. and McIntyre, J. (2011) + Estimating animal densities and home range in regions with irregular + boundaries and holes: a lattice-based alternative to the kernel + density estimator. \emph{Ecological Modelling} \bold{222}, 1666--1672. + + Botev, Z.I. and Grotowski, J.F. and Kroese, D.P. (2010) + Kernel density estimation via diffusion. + \emph{Annals of Statistics} \bold{38}, 2916--2957. +} +\author{ + \adrian and \tilman. +} +\examples{ + ## bovine tuberculosis data + X <- subset(btb, select=spoligotype) + plot(X) + P <- relriskHeat(X,sigma=9) + plot(P) +} +\keyword{spatial} +\keyword{smooth}