rrmix {spagmix}

Spatial relative risk surface generation

Description

Generates an appropriately scaled spatial (bivariate) relative risk surface using a supplied control density and N isotropic Gaussian-style hotspots.

Usage

rrmix(g, rhotspots, rsds, rweights, rbase = 1, log = TRUE)

Arguments

g

A pixel image representing the control density; this will be internally rescaled to integrate to 1 if it does not already do so.

rhotspots

A 2 \times N matrix giving the centers of the N peaks and troughs in the relative risk density.

rsds

A positive numeric vector of length N giving the isotropic standard deviations for each relative Gaussian peak or trough.

rweights

A vector of length N giving relative weightings for each peak (positive weights) or trough (negative).

rbase

The base level of the relative risk surface (default is 1). The peaks and troughs will be added or subtracted from this base level prior to normalisation.

log

A logical value. If TRUE (default), the relative risk surface is returned logged.

Details

A useful tool for the comparison of two estimated density functions on the same spatial region W \subset R^2 is the relative risk function, r, (Bithell, 1990; 1991; Kelsall and Diggle, 1995), defined simply as a density-ratio:

r(x) = f(x) / g(x); x \in W.

Various methods have been developed to improve estimation of r, most commonly with a motivation in geographical epidemiology, where the ‘numerator’ density f pertains to the observed disease cases and the ‘denominator’ density g reflects the distribution of the at-risk controls (Kelsall and Diggle, 1995; Hazelton and Davies, 2009; Davies and Hazelton, 2010). To test newly developed methodology, simulations based on known relative risk scenarios are usually necessary. This function allows the user to design such scenarios, as used in Hazelton and Davies (2009), Davies and Hazelton (2010), and Davies (2013) for example.

This function calculates a relative risk surface based on N Gaussian-style ‘bumps’ added and subtracted from a base level of rbase, with the peaks and troughs centered at the coordinates given by rhotspots with relative weights of rweights and isotropic standard deviations of rsds. The risk surface r is computed as

r(x) \propto rbase + \sum_{i=1}^{N} rweights[i]*exp(-0.5*rsds[i]^(-2)*||x-rhotspots[,i]||^2)

where || . || denotes Euclidean norm. Because f and g are both densities, the risk surface as defined above must then be rescaled with respect to the supplied control density g (argument g) to ensure that

\int_W r(x)g(x) dx = 1

This is automatically performed inside the function. The case density that gives rise to the designed r is then easily recovered because f = r * g. By default, the function returns the log-relative risk surface \log r = \log f - \log g alongside the case and control densities.

Value

An object of class rrim. This is a solist of three pixel images: f as the case density, g the control density (a copy of the argument of the same name, integrating to 1), and r as the (log) relative risk surface.

Author(s)

A.K. Redmond and T.M. Davies

References

Bithell, J.F. (1990), An application of density estimation to geographical epidemiology, Statistics in Medicine, 9, 691-701.

Bithell, J.F. (1991), Estimation of relative risk functions, Statistics in Medicine, 10, 1745-1751.

Davies, T.M. (2013), Jointly optimal bandwidth selection for the planar kernel-smoothed density-ratio, Spatial and Spatio-temporal Epidemiology, 5, 51-65.

Davies, T.M. and Hazelton, M.L. (2010), Adaptive kernel estimation of spatial relative risk, Statistics in Medicine, 29(23) 2423-2437.

Kelsall, J.E. and Diggle, P.J. (1995a), Kernel estimation of relative risk, Bernoulli, 1, 3-16.

Examples

set.seed(1)
gg <- rgmix(3,window=toywin,S=matrix(c(0.08^2,0,0,0.1^2),nrow=2),p0=0.2)

rho <- rrmix(g=gg,
             rhotspots=cbind(c(0.8,0.3),c(0.4,0.4),c(0.6,0.5),c(0.3,0.5)),
             rsds=c(0.005,0.025,0.01,0.025),
             rweights=c(3,2,10,5)*10)

rho.sample <- rrpoint(c(400,800),rho,toywin)


oldpar <- par(mfrow=c(2,2))
plot(rho$g,main="control density")
plot(rho$f,main="case density")
plot(rho$r,main="log relative risk surface")
plot(rho.sample$controls,main="sample data")
points(rho.sample$cases,col=2)
legend("topright",col=2:1,legend=c("cases","controls"),pch=1)
par(oldpar)