3.3 Spatial variation in risk

This section could arguably be renamed “intensity ratios”, as the main idea is to look at the quotient of nonparamteric estimates of the intensity of the process describing the cases and the process describing the controls. However, the theoretical framework for the methods in this section derives from nonparametric density estimation (i.e. the estimation of probability density functions), so we stick with the term “density ratios”.

This method first uses kernel smoothing for separate estimation of $\lambda (x)$ and ${\lambda }_0(x)$, then uses these estimates in the expression $\lambda (x)=\alpha {\lambda }_0(x)\rho (x)$ to estimate $\rho (x)$.

The basic problem addressed by the kernel smoothing method is to estimate an unknown probability density function from a observed sample of data. This problem is very familiar if one is prepared to assume a parametric form for the density function, $f(x)$ say, as then all that is required is to estimate the parameters of $f(x)$, for example using maximum likelihood estimation.

It often occurs, especially in estimating unknown intensity functions, that it is not practical to assume a particular parametric form. In that case, non-parametric or semi-parametric methods are required.

Consider a sample $x_1,\mathrm{\dots },x_n$ from an unknown univariate density function $f(x)$. A simple but rather unsatisfactory method is to treat $f(x)$ as a discrete probability distribution, resulting in the empirical distribution function:

A natural improvement is based on histogram methods. Note that the bandwidth (i.e. width of each bar) has an impact on the appearance of the estimate of $f(x)$.

Kernel smoothing improves further on these simple estimates by treating the underlying distribution as having a continuous, smooth density function. Specifically, we estimate $f(x)$ by

$$\widehat{f}(x;h)=h^{-1}\sum _{i=1}^{n}W\{(x-x_i)/h\}.$$

Here, $W(t)$ is a symmetric (even) function, non-increasing in $|t|$, known as the kernel, and $h>0$ is a scalar which determines the amount of smoothing. Each data-point $x_i$ contributes to the estimate of $f(x)$ at each possible value of $x$, with the largest contributions coming from data-points close to $x$.

Possible choices for $W(\cdot )$ include the standard Gaussian density function, the uniform density function on a specified interval, and the Epanechnikov kernel:

Another way of looking at the effect of changing $h$ appears in the graphs in Figure 3.5, for which the data are a random sample of size $n=25$ from the standard Gaussian distribution.

One method for bandwidth determination is least-squares cross-validation, see Wand and Jones (1994). We would like to choose $\widehat{f}$ to minimise the integrated mean square error

$$R(\widehat{f},f)=𝔼\left[\int {(\widehat{f}(x;h)-f(x))}^{2}dx\right].$$

We cannot do this directly, as $f(x)$ is unknown, but instead we can estimate the part of this function that depends on $h$,

$$J(h)=\int {[\widehat{f}(x;h)]}^{2}dx-2\int \widehat{f}(x;h)f(x)dx,$$

by

$$\widehat{J}(h)=\int {[\widehat{f}(x;h)]}^{2}dx-\frac{2}{n}\sum _{i=1}^n{\widehat{f}}_{(-i)}(x;h),$$

where

$${\widehat{f}}_{(-i)}(x;h)=\frac{n}{(n-1)}{h}^{-1}\sum _{j\ne i}^{n}W\{(x-x_j)/h\}.$$

is the so-called ‘jackknife’ estimate of $f(x)$, obtained using all except the $i$th data-point. Then choose $h$ to minimise $\widehat{J}(h)$. Stone’s theorem (Theorem 3.1 below) guarantees (under certain conditions) that, in large samples, the estimate $\widehat{h}$ obtained using cross-validation converges to the value that minimises $J(h)$. Further technical details are provided by Scott and Terrell (1987).

In a spatial setting, the method is similar. The $x_i$ are now case or control locations (treated separately), and $W(\cdot )$ is a circularly symmetric bivariate density function. The same value of $h$ is usually used to estimate both $\lambda (x)$ and ${\lambda }_0(x)$.

Since $\lambda (x)=\alpha {\lambda }_0(x)\rho (x)$, we have

$$\widehat{\rho }(x)=\frac{\widehat{\lambda }(x)}{\widehat{\alpha }{\widehat{\lambda }}_0(x)}$$

.

An estimate of the probability that a given individual at location $x$ is a case is

$$\widehat{p}(x)=\frac{\widehat{\lambda }(x)}{\widehat{\lambda }(x)+{\widehat{\lambda }}_0(x)}$$

.

• •

  For $n$ cases and $m=N-n$ controls, let $Y_i=1/0$ for case/control respectively at location $x_i$.

• •

  Define weights,

  $$w_i(x)=\frac{W\{(x-x_i)/h\}}{{\sum }_{j=1}^{N}W\{(x-x_j)/h\}}$$

• •

  Kernel estimator is

  $$\widehat{p}(x)=\sum _{i=1}^N{w}_i(x)Y_i$$

  and so

  $$\widehat{\rho }(x)=\frac{\widehat{p}(x)}{\widehat{\alpha }\{1-\widehat{p}(x)\}}.$$

As in the kernel smoothing method, the amount of smoothing is determined primarily by the bandwidth, $h$. Kelsall and Diggle (1998) recommend cross-validation to choose $h$, similar to the cross-validation method described above. For binary regression, with no explanatory variables, this is defined as follows:

• •

  for each $h$, let ${\widehat{p}}_{(-i)}(x_i)$ be the estimate of $p(x_i)$ using all data except $Y_i$;

• •

  choose $h$ to maximise the cross-validated likelihood, or equivalently minimise $CV(h)$ (minus the cross-validated log-likelihood), defined as

  $$CV(h)=-\sum _{i=1}^n\log {\widehat{p}}_{(-i)}(x_i)-\sum _{i=n+1}^N\log \{1-{\widehat{p}}_{(-i)}(x_i)\}$$

Jarner et al. (2002) develop the modifications needed to analyse matched case-control data. Even though in general the matched design is not recommended when spatial variation is of scientific interest, if a matched design is used the method of analysis must respect the design.

If explanatory variables are included, the generalized additive model (below) can be used, which uses a similar idea by incorporating a cross-validated likelihood step in order to estimate residual risk after estimating the effects of covariates on $p(x)$.

Generalized additive models (GAMs) are an extension of generalized linear models (GLMs) that allow a more flexible functional relationship between covariates and the outcome of interest. In GLMs, this relationship is typically assumed to be linear, or even if not, must be specified in advance. In much the same way as the kernel smoothing methods described above, GAMs allow the form of this relationship to be estimated as part of the model-fitting procedure.

The model assumption is

$$\log \left(\frac{p(x)}{1-p(x)}\right)=u{(x)}^{\prime }\beta +g(x)$$

where $u(x)$ is vector of known risk factors and $g(x)$ is a function that models smooth residual spatial variation.

In a spatial setting, $x$ is a geographical location, but terms equivalent to $g(x)$ could also be included to investigate the nature of measured covariates, both in spatial and non-spatial regression analyses.

We have assumed a logit link, $\log \{p(x)/[1-p(x)]\}$, canonical for the Bernoulli distribution. As for GLMs, other link functions could also be specified, and $g(x)$ could be replaced by ${\sum }_{i=1}^p{g}_i(x_i)$ for greater generality.

We have $\widehat{\alpha }\widehat{\rho }(x)=\exp \{u{(x)}^{\prime }\widehat{\beta }+\widehat{g}(x)\}$, and the log-likelihood is of the form

$$\sum _{i=1}^n\log p(x_i)+\sum _{i=n+1}^N\log \{1-p(x_i)\}.$$

Details of the fitting algorithm are available in Hastie and Tibshirani (1990), and include a kernel smoothing step for $\widehat{g}(x)$ within iteratively weighted least squares. The method is implemented in the R library mgcv.