rrmix {spagmix} | R Documentation |
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 |
rhotspots |
A |
rsds |
A positive numeric vector of length |
rweights |
A vector of length |
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 |
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 im
ages: 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)