PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Ann Appl Stat. Author manuscript; available in PMC 2010 April 21.
Published in final edited form as:
Ann Appl Stat. 2008 October 8; 3(3): 943–962.
doi:  10.1214/09-AOAS240
PMCID: PMC2857924
NIHMSID: NIHMS128010

ANALYSIS OF MINNESOTA COLON AND RECTUM CANCER POINT PATTERNS WITH SPATIAL AND NONSPATIAL COVARIATE INFORMATION

SHENGDE LIANG and BRADLEY P. CARLIN
Public Health in the Division of Biostatistics School of Public Health University of Minnesota Minneapolis, Minnesota 55455 USA

Abstract

Colon and rectum cancer share many risk factors, and are often tabulated together as “colorectal cancer” in published summaries. However, recent work indicating that exercise, diet, and family history may have differential impacts on the two cancers encourages analyzing them separately, so that corresponding public health interventions can be more efficiently targeted. We analyze colon and rectum cancer data from the Minnesota Cancer Surveillance System from 1998-2002 over the 16-county Twin Cities (Minneapolis-St. Paul) metro and exurban area. The data consist of two marked point patterns, meaning that any statistical model must account for randomness in the observed locations, and expected positive association between the two cancer patterns. Our model extends marked spatial point pattern analysis in the context of a log Guassian Cox process to accommodate spatially referenced covariates (local poverty rate and location within the metro area), individual-level risk factors (patient age and cancer stage), and related interactions. We obtain smoothed maps of marginal log-relative intensity surfaces for colon and rectum cancer, and uncover significant age and stage differences between the two groups. This encourages more aggressive colon cancer screening in the inner Twin Cities and their southern and western exurbs, where our model indicates higher colon cancer relative intensity.

Keywords: Colon cancer, log Guassian Cox process (LGCP), rectum cancer, spatial point process

1. Introduction

1.1. Etiologies of colon and rectum cancer

Traditionally, public health agencies have reported colon and rectum cancers together under the title “colorectal cancer.” Since the turn of the last century, however, an active debate has emerged regarding whether these two cancers really have sufficiently similar etiologies to be aggregated in this way. Some experts have argued that the cancers should be reported and monitored separately, so that public health interventions can be more sensibly and efficiently targeted.

A variety of epidemiological studies have indicated variables that may have a differential impact on colon and rectum cancer. These variables fall under three broad categories. The first is exercise. A very recent study by the Physical Activity Guidelines Advisory Committee (2008) identified 23 publications on this general topic, arising from 12 prospective cohort studies and 8 case-control studies. These studies show a consistent inverse relation between physical activity and colon cancer risk, with this relation being statistically significant for at least one physical activity domain and one sex in 9 of the 12 cohort studies and 5 of the 8 case-control studies. More specifically, the median relative risk (RR) comparing most-versus least-active subjects was 0.7 over all the studies. The advisory committee stated that this overall finding was unlikely to be the result of confounding, since the studies for the most part included relevant covariates, such as body mass index (BMI), smoking, alcohol, diet, screening, menopausal status, and family history of colon cancer. By contrast, the committee found the studies to indicate no apparent relationship between physical activity and rectal cancer risk. Specifically, more than half the studies showed no statistically significant associations, and the median RR over all the studies was 1.0. The fact that moderate-to-vigorous physical activity (say, 30 to 60 minutes per day) may be protective against colon cancer but not rectal cancer suggests that these two cancers should be treated separately in cancer registry reporting and subsequent statistical modeling.

A second broad area that may have a differential impact on the two cancers is diet. Diet has long been suspected as an etiological factor for colorectal cancer; however, studies of individual foods and nutrients have often provided inconsistent results, perhaps due to low statistical power. Flood et al. (2008) address this problem using factor analysis to group dietary variables into three broad groups, and go on to conclude that lower consumption of meat and potatoes, and higher consumption of fruit, vegetables, and fat-reduced foods, are associated with reduced colorectal cancer risk. Wei et al. (2003) use data from two prospective cohort studies (87,733 women from the Nurses' Health Study and 46,632 men from Health Professionals Follow-Up Study) to investigate the effect of dietary variables on colon and rectum cancer separately. In the combined cohort, a variety of variables emerge as significant predictors of elevated colon cancer risk, including intake of beef, pork or lamb as a main dish, intake of processed meat, and alcohol consumption. However, none of these variables emerge as predictors of rectal cancer. Using data from the Iowa Women's Health Study, Folsom and Hong (2005) showed that magnesium and calcium intake were independently associated with significantly lower colon cancer risk, but not rectum cancer risk. The relative risk estimates changed little across baseline subgroups, such as women who did or did not use hormone replacement therapy, or were or were not diabetic. Using data from a different study, Flood et al. (2005) conclude that the protective effect of calcium is present regardless of whether the calcium arises naturally in food, or is delivered through dietary supplements. Finally, Pedersen, Johansen and Gronbaek (2003) observed a dose-response relationship between alcohol and rectal cancer in a Danish cohort of 15,491 men and 13,641 women who did not include wine in their alcohol intake. However, no association between alcohol and colon cancer was found.

The third broad area where the etiologies of colon and rectal cancer appear to differ is family history. For those who reported a family history of colon or rectal cancer, Fuchs et al. (1994) obtained a RR of 1.99 for colon cancer but just 0.86 for rectal cancer, a statistically significant difference based on a simple chi-square test with one degree of freedom. Wei et al. (2003) drew the same conclusion, but using a different dataset and a stepwise polytomous logistic regression procedure.

1.2. MCSS data and problem description

Our specific problem of interest involves the comparison of the spatial distributions of colon and rectum cancer patients in the state of Minnesota. These data are collected by the Minnesota Cancer Surveillance System (MCSS), a program sponsored by the Minnesota Department of Health. The MCSS includes the residential address of essentially every person diagnosed with cancer in Minnesota. Here we consider the subset of patients diagnosed during the period 1998-2002 (an interval chosen partly for its centering around a U.S. Census year, 2000). Figure 1 shows the 7 counties comprising the Twin Cities metro area as those encircled by the dark boundary; also shown are 9 adjacent, exurban counties. Within these 16 counties, we have 6544 individuals for analysis. Figure 1 plots the approximate locations of the cancers after the addition of a random “jitter” to protect patient confidentiality (explaining why some of the cases appear to lie outside of the spatial domain). The physiological adjacency of the colon and the rectum suggests positive dependence in these point patterns; persons with rectum cancer beyond stage 1 (i.e., regional or distant) are at risk for colon cancer due to metastasis. Moreover, the two cancers likely share unmodeled spatially-varying risk factors (such as local health care quality or availability), also implying positive dependence. This may help health care providers or public health policy makers to identify regions of excessive risk requiring intervention (say, a direct mail campaign encouraging more aggressive screening) or other weak links in the health care system.

Fig. 1
Jittered residential locations of colon (light circle) and rectum (dark circle)cancer cases, Twin Cities metro and exurban counties, 1998-2002.

The causes of colon and rectum cancer are unknown. Age is the primary risk factor, with disease incidence increasing significantly after the age of 50. As already mentioned, family medical history may also be helpful in predicting colon and rectum cancer risk, along with several lifestyle factors such as alcohol use, smoking, diet, and exercise. Unfortunately we do not have access to this information for individuals in the MCSS, but the lifestyle factors could reasonably be expected to cluster spatially due to corresponding sociodemographic clustering. We also have census tract-level poverty rates, which should be correlated with these risk factors.

A full analysis of the data in Figure 1 would account for the randomness in the observed locations, their spatial correlation, important covariates (including population density), and any other hierarchical structure in the data (such as the tendency of model residuals to cluster spatially). The output of such an analysis would include maps of the fitted adjusted log-intensity surface, point and interval estimates for important main effects and interactions (e.g., location-age), and perhaps maps of fitted spatial residual surfaces, to help identify spatial covariates still missing from the model.

1.3. Statistical modeling of spatial point patterns

In spatial disease mapping settings, the primary goals are typically to investigate the connections between the disease and (possibly geographically-indexed) covariates, to characterize the spatial variation of the disease occurrence, and to identify areas having elevated disease risk. In such cases, the data are often aggregated to counts within specified areal regions (counties, zip codes, etc.). Indeed, most published statistical analyses to date of data of this type use so-called areal or lattice models; see, for example, Banerjee, Carlin and Gelfand [(2004), Chapter 3 and Section 5.4] for a review. However, if precise geocoded locations of disease cases are available, it is more appealing to study the resulting spatial point pattern using spatial point process modeling. However, such methods are conceptually and computationally more challenging, and are implemented in fewer widely available statistical software programs. Indeed, even when actual geocoded locations are available, a standard computational strategy is to partition the study region and model the counts in each cell of the partition as conditionally independent Poisson observations, obtaining the standard areal model but with an arbitrary partition.

Under a nonhomogeneous Poisson process, the likelihood for the intensity surface generating the locations given the observed locations is well known [see, e.g., Beneš et al. (2002); Diggle (2003); Møller and Waagepetersen (2003)]. We begin with this likelihood, but introduce the following features. First, we accommodate covariate information in a novel way. We envision certain covariates as conditional, that is, we seek to compare point patterns given levels of these covariates. For us, these are cancer type covariates which “mark” the point pattern. We view other patient level characteristics (or risk factors) such as age as nuisance variables for which we wish to adjust. We then model point patterns jointly over geographic space and nuisance covariate space, enabling the notions of both conditional and marginal intensity associated with geographic space. Hence, we obtain an intensity adjusted for these covariates, rather than an intensity which ignores them by not including them in the model. Moreover, we also have purely spatial covariates, some of which are available at areal unit level (say, county-level features), while others may be available at point level (say, distance from a location to the nearest cancer screening facility). Employing spatial information at both scales precludes aggregation of points to counts.

Additionally, we anticipate dependence between the intensity surfaces associated with the two cancers, since, for example, an excess of colon cancer in a portion of the study region may suggest correspondingly high levels of rectum cancer. We capture this dependence using multivariate process realizations for the intensities. Last, working with the above point level likelihood, as well as fairly large numbers of points (e.g., order 103), necessitates approximation to implement the model-fitting.

The analysis of spatial point patterns has a reasonably long history in the literature, initially built using exploratory tools such as distance based methods yielding F functions, G functions, and, perhaps most commonly, Ripley's K function. All are based upon assessing departure from complete spatial randomness (CSR), which is interpreted as a homogeneous Poisson process and for which closed forms for these functions exist. However, no likelihood is specified, and comparison between point patterns is not possible. Another more recent approach involves the use of spatial scan statistics, currently popular in large part due to the SatScan software of Kulldorff (2006). But again, no likelihood is specified so inference is limited to say detection of “hot spots.”

To achieve the foregoing objectives, we instead adopt a model-based focus, and write the intensity of the process as λ(s),where s [set membership] D for some spatial domain D. For a collection of observed cancer case locations si,i = 1,…,n, we work with the likelihood L(λ(s),sD;{si}i=1n) which takes the form eDλ(s)dsΠi=1nλ(si) Often, λ(s) is specified as a parametric function, for example, using a basis representation or a tiled surface. Adding a prior distribution on these parameters, say, θ, yields a posterior distribution p(λ(s; θ)|{si}) for making inferences about the intensity surface.

For us, λ(s) is thought of as a log Guassian process (GP) realization, resulting in the familiar class of Cox processes [Møller and Waagepetersen (2004), page 57]. To specify this prior distribution, we require μ(s), the mean surface, along with σ2 and [var phi], the GP covariance parameters. Below, we express μ(s) in part with a form z′ (s)β, so that the process mean can depend on spatially referenced covariates z(s).

A common class of estimation methods for inhomogeneous spatial point process models avoids full likelihood evaluations by formulating estimating equations [Waagepetersen (2007); Waagepetersen and Guan (2008)]. Guan and Loh (2007) study the distributional properties of the estimation procedure of Waagepetersen (2007), and obtain variance estimates using a thinned block bootstrap procedure. Guan (2006) developed a composite likelihood method based on the second-order intensity function of the underlying process. Diggle and Rowlingson (1994) handle bivariate (case-control) point processes via a conditional likelihood approach to convert the two spatial point process models into an easier-to-fit nonlinear binary regression model. Similarly, Guan, Waagepetersen and Beale (2008) estimate correlation functions via either a consistent nonparametric kernel smoothing estimator, or a parametric conditional likelihood estimator. In all of these approaches, inference on spatial associations and second-order variations proceeds not from the intensity surface, but from pairwise correlation functions and transforms thereof [e.g., the g and K functions in Waagepetersen (2007)]. As such, they do not offer direct attacks on the intensity surface estimation problem, needed for inference regarding the fitted surface itself, its rate of change at any point [as needed for spatial boundary analysis or “wombling”; see Banerjee and Gelfand (2006), and Liang, Banerjee and Carlin (2009)], or model-based comparison of the surfaces for colon and rectum cancer.

As such, we instead adopt a fully Bayesian approach that yields posterior distributions for the intensity surface, or even the spatial residual surface after adjusting for regressors that are allowed to differ for the two cancers. Due to the absence of sufficient covariate (e.g., diet) information in our dataset, we introduce spatially varying random effects which we view as surrogates for these missing covariates. Inference is exact and does not rely upon possibly inappropriate use of infill or increasing-domain asymptotics. However, our more comprehensive approach comes with a price. Specifically, note that if λ(s) is modeled as a random realization of a spatial process, then the likelihood integral is stochastic, precluding explicit evaluation. Indeed, a variety of computational challenges emerge in working with the point-level likelihood in this case: the stochastic integration, the large collection of spatial locations, and a prior specification that is only available through finite dimensional distributions.

Wolpert and Ickstadt (1998) offered one of the first fully Bayesian approaches for spatially nonhomogeneous Poisson process data. Beneš et al. (2002) illustrate one possible Bayesian analysis of a log Guassian Cox process model. While they assume λ(s) constant over grid cells, they do utilize the notion of the population intensity surface, and obtain fitted disease maps under a variety of models (constant, Gaussian kernel, etc.) for this surface. However, they do not consider joint modeling of multiple disease surfaces, nor covariates that are not location-specific (e.g., the age or cancer stage of a case observed at a particular location).

In this paper we employ novel spatial point process approaches that account for both location-specific and nonlocation-specific covariates in the context of multiple dependent point processes to analyze the MCSS dataset. We begin in Section 2 with a brief review of spatial point process modeling. We then go on to present a set of multivariate spatial point process models for investigating the effect of both location-specific and nonlocation-specific covariates, as well as their interactions. Section 3 then gives the results of applying them to our MCSS dataset. Finally, Section 4 discusses our findings and offers directions for future research in this area. Computational challenges in fitting our models are addressed in the supplemental article [Liang, Carlin and Gelfand (2009)].

2. Hierarchical modeling for spatial point processes

2.1. Modeling with spatial covariates

We begin with a brief review of the basics of log Guassian Cox process modeling. Consider a set of random locations which we denote by S={si}i=1n where disease occurrence is observed over a spatial domain D. We model this random set of locations using a nonhomogeneous Poisson process with intensity function λ(s) for all s [set membership] D. Let z(s) be a vector of location-specific covariates corresponding to a disease case observed at s. For us, a key component of z(s) is the indicator of whether the case is in the metro area or not. However, in other contexts, we could envision information such as elevation, climate, exposure to pollutants, and so on to be relevant. We model λ(s) = r(s)π(s), where r(s) is the population density surface at location s. In practice, we may create such a surface using GIS tools and census data, or we may just work with areal unit population counts, letting r(s) = n(A)/|A| if s [set membership] A, where n(A) is the number of persons residing in A and |A| is the area of A. The error introduced by this admittedly crude estimate may be mitigated somewhat by resorting to the non-Bayesian estimating equation alternatives discussed in Section 1.3. Specifically, in this framework one could model the two point processes separately [Waagepetersen (2007)] or jointly [by an extension of Guan (2006)].

Returning to our framework, r(s) serves as an offset and π(s) is interpreted as a population adjusted (or relative) intensity, which we model on the log scale as

π(S)=exp(Z(S)β+w(S)),
(1)

where w(s) is a zero-centered stochastic process, and β is an unknown vector of regression coefficients. If w(s) is taken to be a Gaussian process, then the original point process is called a log Gaussian Cox process [LGCP; Møller and Waagepetersen (2004), page 72]. The likelihood associated with β and wD = {w(s):s [set membership] D} given S takes the form

L(β,wD;S)exp(Dr(S)π(S)dS)×siSr(Si)π(Si).
(2)

Operating formally, a prior on wD along with a prior on β completes the Bayesian specification. Inference proceeds from the posterior which, again formally, is p(β, wD|S) [proportional, variant] L(β, wD; S)p(β)p(wD). Of course, the Gaussian process is only defined through its finite dimensional distributions so that, practically, this posterior is viewed in terms of a finite collection of locations. This motivates discrete approximation of the stochastic integral as we discuss below. One discrete approximation partitions D into a collection of sets (say, Ai, i = 1, 2, … ,m) and creates a Poisson likelihood for the counts given λ(Ai). That is, it models log π(Ai), thus precluding use of point level covariate information. Moreover, since π(Ai) = ∫Ai π(s) ≠ = exp(∫Ai (z(s)′β + w(s)) ds, it is inappropriate to utilize the latter, simpler integration. Indeed, ignoring this inequality can introduce ecological fallacy issues; see, for example, Wakefield and Salway (2001) for a discussion.

We pursue an alternative discrete approximation which still enables us to work at the point level. Suppose we replace ∫Dλ(s) with some choice of numerical integration. For the moment, we allow analytic possibilities as well as Monte Carlo versions, since in either case, we will end up replacing wD with a finite set, say, w={w(sj),j=1,2,,T}. Then we revise (2) to

L(β,w,w(s1),...,w(sn);S)p(w,w(Si),...,w(Sn))p(β).
(3)

Now, we only need to work with an (n + T)-dimensional random variable to handle the w's, hence, their prior is just an (n + T)-dimensional multivariate normal distribution. Note that, in (3), we will require that z(s) be available at each tj; that is, we require the component z(s) surfaces over D. These surfaces are not viewed as random and may be interpolated or tiled, according to the nature of the information for the particular spatial covariate; we assume only that we can assign a value of z for each s [set membership] D.

2.2. Introducing nonspatial covariate information

So far we have indicated how to incorporate covariates that are spatially referenced into the modeling. In our setting, we seek to introduce nonspatial covariates which we think of as being of two types (though the distinction will depend upon the application). One type of covariate provides the “marks” leading to a marked point process model. For us, this covariate is cancer type (colon vs. rectum), and we are interested in whether the two cancer intensity patterns differ.

The second type of covariate we view as an “auxiliary” variable that provides additional information associated with intensity. For us, age and cancer stage are examples of such covariates. Clearly patient age is associated with cancer intensity, but the strength of this association may differ across cancers. We wish to adjust intensity to reflect patient age, analogous to the age standardization used in aggregated areal data settings.

In general, we view these latter covariates as continuous2 and introduce a second argument into the definition of the intensity, yielding a surface in s and v over the product space D × ν (i.e., geographic space by covariate space). We then generalize (1) to

π(s,v)=exp(β0+z(s)β+vα+(vz(s))γ+w(s)).
(4)

where the Kronecker product v [multiply sign in circle] z(s) denotes the set of all the first order multiplicative interaction terms between z(s) and v. When a particular interaction term is not of interest, the corresponding coefficient in γ is set to zero. This expression envisions a conceptual intensity value at each (s, v) combination. The interaction terms between spatial and nonspatial covariates provide the ability to adjust the spatial intensity by individual risk factors. If we fix ν in (4), we can view λ(s, v) = r(s)π(s, v) as a “conditional” intensity at level v. If we integrate over v (see below), we obtain the (cumulative) marginal intensity λ(s) associated with π(s, v).

Now, introducing marks k = 1, 2, …, K, a general additive form for the log relative intensity is

logπk(s,v)=β0k+z(s)βk+vαk+(vz(s))γk+wk(s).
(5)

We can immediately interpret the terms on the right side of (5). The global mark effect is captured with the β0k. Therefore, there is no intercept in z(s) and we have mark-varying coefficients for the spatially-referenced covariates, reflecting the possibility that these covariates can differentially affect the intensity surfaces of the marks. Similarly, we have mark-varying coefficients for the nuisance variables. We also have mark-varying coefficients for the interaction terms, reflecting possibly different effects of the nonspatial covariates over spatial domains. Finally, we allow the spatial random effects to vary with mark, that is, a different Gaussian process realization for each k. Dependence in the wk(s) surfaces may be expected (say, increased intensity at s for one marked outcome encourages increased intensity for another at that s), suggesting the need for a multivariate Gaussian process over the wk. Both separable and nonseparable forms for the associated cross-covariance function are conveniently specified through coregionalization [Gelfand et al. (2004); Banerjee, Carlin and Gelfand (2004), Sections 7.1 and 7.2].

Reduced models of (5) are immediately available, including, for example, wk(s) = w(s), βk = β, and αk = α. Another interesting reduced model obtains by setting γk = 0, leading to

logπk(s,v)=β0k+z(s)βk+vαk+wk(s).
(6)

This separable form enables us to directly study the effect of the marks on spatial intensity. Specifically, the intensity associated with (5) is

λk(s,v)=exp(β0k+vαk)×r(s)exp(z(s)βk+wk(s)).
(7)

We see a factorization into nonspatial nuisance and spatial covariate terms. Presuming the former is integrable over v, the latter, up to a constant, is the “marginal spatial intensity.”

Integration of λk(s, v), based upon (5), can be computed analytically in most cases. When v is categorical, the likelihood integral involves only integration over the spatial domain D. When v is continuous, simple algebra shows

νλk(s,v)dvds=r(s)exp(β0k+z(s)βk+wk(s))×νexp(vαk+(vz(s))γk)dv.

Suppose, for instance, that there is only one component in z(s) and one component in v having range (vl, vu). Provided αk + z(sk ≠ 0, the marginal intensity λk(s) is then

νλk(s,v)dvds=r(s)exp(β0k+βkz(s)+wk(s))×νexp(v(αk+z(s)γk))dv=r(s)exp(β0k+βkz(s)+wk(s))×1αk+z(s)γk[exp(vu(αk+z(s)γk))exp(vl(αk+z(s)γk))].

Turning to the revised likelihood associated with (5), let {ski, vki), i = 1, 2, … ,nk} be the locations and nuisance covariates associated with the nk points having mark k. The likelihood becomes

kexp(Dνλk(s,v)dvds)×kski,vkiλk(ski,vki).
(8)

Using the calculations above, the double integral becomes

Dνλk(s,v)dvds=D(r(s)exp(β0k+βkz(s)+wk(s))×1αk+z(s)γk[exp(vu(αk+z(s)γk))exp(vl(αk+z(s)γk))])ds,

provided that the set {sk + z(sk = 0} has Lebesgue measure zero. Hence, the difficulty in the likelihood evaluation is the same as in (2) and can be treated in the manner described in conjunction with (3). In this regard, note that we bound the components of v in order to integrate explicitly over v. We do not have a stochastic integration with regard to V as we have over D. Of course, sensitivity to the chosen bounds should be investigated.

Last, in the case of k = 2 marks, a common alternative model specification is logistic regression, which views the mark as the response given the locations and covariates. This is conditioning in the opposite order from our model, which views the locations and covariates as random given the marks. In our dataset it seems more natural to compare point patterns for the two different types of cancer, rather than view cancer type as some sort of binary “response” to covariate information.

3. Results

We now present the results of our analysis of the MCSS colon and rectum cancer data. Previous studies suggest that covariates related to a patient's socioeconomic status (SES) may be related to the patient's risk factors through its impact on diet, health care quality, or propensity to seek care. While our dataset lacks any individual-level SES measures, from census data we have several related tract-level variables: percentage of farm population, percentage of rural population, percentage of people with less than high school education, percentage of minority population, and poverty rate. A preliminary population-adjusted nonspatial Poisson regression analysis of our data on these covariates revealed only poverty rate and the metro indicator as significant predictors.

In our initial model, we consider two location-specific covariates: z1(s), the metro area indicator, and z2(s), the poverty rate in the census tract containing s. We also employ two nonlocation-specific covariates: v1, cancer stage [set to 1 if the cancer is diagnosed “late” (regional or distant stage) and 0 otherwise], and v2, the patient's age at diagnosis. The population density r(s) we use for standardization is available at 2000 census tract level, meaning that we assume population density is constant across any tract. The integral of the intensity is approximated by a Monte Carlo sum using a predictive process approximation [Banerjee et al. (2008)]; see the supplemental article by Liang, Carlin and Gelfand (2009) for full details.

The left and middle columns of Figure 2 show maps of the raw mean non-spatially varying covariates (age and proportion diagnosed late), while the right column maps a crude estimate of relative intensity for colon cancer (top row) and rectum cancer (bottom row). Notice these summaries are presented at tract level, even though we have exact (or nearly exact) spatial coordinates here. In the first two columns, tracts containing no cases are simply shaded according to the overall observed mean values for each disease, which are 69.9 and -64.8 for age and 0.618 and 0.555 for proportion diagnosed late for colon and rectum cancer, respectively.

Fig. 2
Minnesota colorectal cancer covariate and response data for colon (top row) and rectum (bottom row) groups: left, tract-specific map of observed mean age; middle, tract-specific map of observed proportion of late diagnosis; right, tract-specific observed ...

None of these four maps show strong spatial patterns, though we do see several areas with higher than average age, late diagnosis fraction, or both. The right column maps the logs of the numbers of cases divided by total number of residents in each tract. These crude maps of the tract-level log relative intensity (unadjusted for any spatial or nonspatial covariates) show somewhat stronger spatial patterns and higher overall rates of colon cancer. The rectum cancer map features an interesting collection of low outlying values in several outer-ring suburban census tracts.

Table 1 breaks down the data by stage and metro/nonmetro area. We see that 38% of colon cancer cases were diagnosed at an early stage, while 44.5% of rectum cancer cases were. In total, colon cancer is nearly three times as prevalent as rectum cancer in both the metro and nonmetro areas. A fact not revealed by the table is that there are 72 individuals who contribute both a colon and a rectum tumor. Since this is only around 1% of the total of 6544 individuals, we do not explicitly model this particular kind of dependence, but rather “lump it in” with the bivariate dependence modeled by ρ.

Table 1
Table of colorectum cancer patients' characteristics in metro and adjacent area of Minnesota.

Figure 3 shows tract-level maps of population density, r(s), and our two location-specific covariates, z1(s) and z2(s). Not surprisingly, the central metro areas are the most populated. The poverty rate is fairly uniform except for high rates in a concentrated portion of the central metro.

Fig. 3
Left, population density by tract; middle, metro/nonmetro area; right, poverty rate by tract .

We now fit our model, using independent Inverse Gamma(2, 0.5) priors for σ12 and σ22, and a Unif(-0.999, 0.999) prior for ρ. The scale of the spatial decay parameter ø is determined by the distance function employed. In this application, we started with a Unif(130, 390) prior for ø, so that the effective range lies between one-fourth and three-fourths of the maximal distance between any two knots. As expected, [var phi] is only weakly identified, so a fairly informative prior is needed for satisfactory MCMC behavior. For simplicity, we simply fix the range parameter at ø = 195, so that the effective range is roughly half of the maximal distance. A random-walk Metropolis-Hastings algorithm is used to draw posterior samples.

Table 2 compares the effective model size and DIC score of three models. It can be seen that the no-random effect model (GLM) is unacceptably bad, and the model with a single set of spatial residuals is not much worse than the bivariate residual model. This suggests that the two sets of residuals are fairly similar, and that ρ is close to 1.

Table 2
Model comparison using effective model size pD and DIC score.

Table 3 shows parameter estimates from some of our models. We parameterize so that the top rows concern the fixed effects for colon cancers, β1, but the second set of rows give the differential effect in the rectum cancer group, β [equivalent] β2 - β1. Thus, any 95% Bayesian confidence intervals that exclude 0 in this part of the table suggest a variable that has a significantly different impact on the two cancers.

Table 3
Parameter estimates for the model with metro indicator and poverty rate as the spatial covariates, and stage and age as individual covariates.

In general, the effects of the non-spatial covariates are fairly similar across the models considered. We find that in the metro area there are relatively fewer cases of both colon and rectum cancer. This is consistent with statewide patterns of colorectal cancer occurrence in Minnesota, where higher age-adjusted rates are often found in nonmetro areas. However, there is no significant change in this relationship in the rectum group relative to the colon group. An interesting and somewhat counterintuitive finding is that poor areas seem to have relatively fewer cases. This appears consistent with the aforementioned finding of Wei et al. (2003) that colon cancer is associated with foods often consumed by relatively more affluent people (beef, pork, or lamb as a main dish, and other processed meat). However, unlike these authors, we find no significant difference in this relationship for rectum cancer.

Turning to the nonlocation-specific covariates, age is significantly associated with increasing colon cancer, but a somewhat surprising relative decrease in rectum cancer. This difference (-0.18) is statistically significant, but not large enough in magnitude to make the overall age effect in the rectum group negative. A look at the data bears this out, with rectum cancers arising in a somewhat younger population; our preliminary Poisson regression also concurs, though here the relative decrease in the rectum group is not significant. Late detection provides another interesting difference between the colon and rectum groups: while there are significantly more cases diagnosed late than early, the effect of late diagnosis is significantly smaller in the rectum group (point estimate -0.26). Thus, public health interventions to encourage screening and early detection of colorectal cancer will have significantly greater impact on prevention for colon than for rectum. The metro-age interaction shows that the effect of age on colon cancer is significantly less pronounced in the metro area; a smaller “age adjustment” to the colon cancer intensity process is needed in the metro area. This effect is largely absent for rectum cancer, but this difference is not quite statistically significant. Finally, the estimate of ρ is very close to 1, indicating very similar spatial residual patterns. This is perhaps a surprisingly strong association, but believable given that these are residual surfaces, which account (at least conceptually) for important missing covariates, which could be spatial (e.g., local screening percentage, other sociodemographic factors) or nonspatial (e.g., the physiological adjacency of the colon and the rectum).

Our Bayesian viewpoint allows us to make probabilistic statements both against and in favor of various null hypotheses of interest. For example, suppose we take 20% as the minimum difference in relative intensity required to conclude a practically meaningful difference between the colon and rectum cancer groups. For predictor i, this amounts to a test of H0i [set membership] [log(0.8), log(1.2) versus Ha : Δi [negated set membership] [log(0.8), log(1.2). The posteriors summarized in the second group of rows in Table 3 enable us to compute posterior odds ratios OR = P(Ha|S)/P (H0|S) for any predictor of interest. In our dataset, the two predictors of greatest substantive interest yield different results. Living in the metro area produces OR = 0.12, or odds of just over 8:1 in favor of no real difference in the colon and rectum However, for late detection we obtain OR = 2.81, or nearly 3:1 odds in favor of a practically meaningful difference (in this case, a relative intensity reduction in the rectum group). Again, this suggests a public health program encouraging more aggressive cancer screening would be sensibly targeted to those living in regions with higher colon cancer relative intensity, since this should lead to a more meaningful reduction in cancer prevalence.

Figure 4 shows maps of the fitted log intensity surfaces both without (left column) and with (right column) the spatial residuals (middle column), for a case at the mean age and diagnosed at an early stage. Without spatial residuals, the two spatial covariates alone predict slightly higher prevalence in the nonmetro areas. However, the residuals (which unlike our spatial covariates are point-level, and are thus summarized using image-contour maps) indicate further reductions are needed in the near southern, southeastern, and northern suburbs, as well as the far north exurban area. This leads to the more mottled fitted patterns in the rightmost column. Note these final two maps in the right column result in spatially smoothed versions of the corresponding maps in the right column of Figure 2, which we recall are something like raw log-relative intensity surfaces. While direct comparison is not really possible since the maps in the right column of Figure 4 are adjusted for both spatial and nonspatial covariates, the overall similarities further confirm the good fit scores achieved by our models. From a practical point of view, when combined with the significant differences in age and late detection between the two cancers found in Table 3, our findings encourage more aggressive colon cancer screening in the inner Twin Cities and the far southern and western exurbs, where the upper right panel of Figure 4 indicates higher colon cancer relative intensity.

Fig. 4
Log-relative intensity surfaces using values at centroid of each census tract at the mean age and assuming an early diagnosis. The top row is for colon cancer and the bottom for rectum cancer. The first column is the log-relative intensity surfaces without ...

4. Discussion

We have offered an analysis of colon and rectum cancer incidence data collected by the Minnesota Cancer Surveillance System during the period 1998-2002. In so doing we extended customary spatial point pattern analysis in the context of a log Guassian Cox process model to accommodate covariates that are spatially referenced, individual-level cancer type marks, and individual-level risk factors that are not of interest in terms of marking. Our approach yields easy-to-interpret fixed effects for testing for equality of epidemiological properties across the two cancers, and fitted maps that can reflect the impacts of the spatially indexed covariates, spatial residuals, or both. These last maps also offer spatially smoothed fitted surfaces reminiscent of those in traditional areal models, but now adjusted for the nonspatially varying covariates, age and cancer stage.

As with many observational data analyses, our findings raise as many questions as they answer. The somewhat counterintuitive negative relationship between tract-level poverty and colorectal cancer shown in Table 3 might be the result of unmodeled confounding between age and poverty: poor areas could very well be significantly younger (especially in the metro, which features a higher proportion of immigrants, who tend to be younger). Since colorectal cancer is so highly associated with age, the apparent beneficial effect of poverty might just be another manifestation of the protective effect of youth. Similarly, modestly negative metro-age interaction may be due to more common use of colorectal screening in the metro area. Such screening methods can reduce colorectal cancer incidence by identifying pre-malignant lesions (polyps) and removing them; failure to screen a population might increase both the number of cases and the ages at which the cancers were diagnosed. Sadly, we currently lack the individual-level income and screening information necessary to precisely address these questions. Moreover, the MCSS database also does not feature information on diet, exercise, or family histories of the patients, the three previously-identified factors most likely to be responsible for any differences between colon and rectum cancer relative hazards. The data collected by MCSS is determined by legislation, and to expand it in any way requires a change in Minnesota state law, attempts at which the Minnesota Department of Health prefers to keep as rare as possible. As a result, future research regarding epidemiological properties of colon and rectum cancer should perhaps focus on obtaining approval and funding for a follow-up questionnaire mailed to all MCSS patients.

Even more fundamentally, an increasing number of authors view the debate over whether colon and rectum cancers have different etiologies as misplaced, arguing that the real distinction is not colon versus rectum, but rather proximal (right, or ascending) versus distal (left, or descending) colon, the latter of which includes the rectum. These authors argue that colorectal cancer is not a single disease, but two distinct diseases with distinct molecular profiles. One of these is more commonly found in the distal colon, and derives from hyperplastic polyps, whose putative successor lesions, serrated adenomas, represent discrete steps along a pathway to cancer [Huang et al. (2004)]. By contrast, the cancers most common in the proximal colon arise from an entirely different molecular pathway [O'Brien et al. (2006)]. Differences in risk factors for these two pathways are not well established but are nonetheless entirely likely. Relatedly, Glebov et al. (2003) found more than 1000 genes expressed differentially in adult ascending versus descending colon. Thus, the real subclassification of interest may not be colon versus rectum or distal colon versus proximal colon, but rather molecular pathology. Of course, this line of thinking encourages a reporting of cancers that is well beyond the capabilities of most U.S. public health reporting (hence intervention) systems, but the idea bears watching.

On the brighter side, recent audits have suggested our MCSS dataset is over 99% complete; that is, due to state reporting requirements, we are aware of essentially every tumor discovered by doctors in Minnesota. However, our methods obviously cannot reflect tumors that are not discovered or otherwise not reported. To the extent that such tumors happen unevenly across the spatial domain, this could lead to bias in our fitted estimates and maps. We do not think differential underreporting is a problem across our current, relatively compact and relatively urban 16-county spatial domain, but datasets that reached further into more remote regions of the state (especially semi-autonomous Native American tribal lands) may well suffer from this problem.

Future work in this area also includes extending our model with more complex interaction terms and perhaps more than two marks (the MCSS database has information on more than 20 cancers), leading to more challenging model-fitting. Another issue to address is the imprecision in the (typically rural) addresses within the point pattern. In some cases error may be simply due to the sensing device (e.g., the GPS unit), while in others it may be due to the practical limits of geocoding: for some of the cancers in our MCSS data, a significant proportion of the geocodes may be based on less than a complete and valid street address (e.g, residence zip + 2, residence zip only, or even the zip of a post office box). A final, perhaps most interesting path for the future lies in space-time point pattern analysis, in order to see evolution of cancer intensities over time. In the case of continuous time, we would now add a time argument to our intensity functions, leading to substantially increased scope for the modeling (e.g., separable versus nonsep-arable models for the space-time intensity). If time is instead viewed as discrete, we might instead extend our framework to temporally dynamic log Guassian Cox process models. Both of these options, while computationally challenging, could pay significant practical dividends.

Supplementary Material

Computational Details

Acknowledgments

The authors are grateful to Drs. Sudipto Banerjee, Andrew Flood, and Aaron Folsom for helpful discussions and Dr. Sally Bushhouse for permitting and facilitating our analysis of the Minnesota Cancer Surveillance System (MCSS) data.

Supported in part by NIH Grant 1-R01-CA95955-01.

Footnotes

2In the case of a discrete valued covariate, any integrals over v in our development are replaced by sums.

Computational Issues (DOI: 10.1214/09-AOAS240SUPP;.pdf). We provide full details of the Monte Carlo algorithms needed to approximate the complex point process likelihoods in the paper. In particular, we flesh out the details of our knot-based predictive process approximation, and give general guidelines for how the knots should be selected in any given application.

REFERENCES

  • Banerjee S, Carlin BP, Gelfand AE. Hierarchical Modeling and Analysis for Spatial Data. Chapman and Hall/CRC Press; Boca Raton, FL: 2004.
  • Banerjee S, Gelfand AE. Bayesian wombling: Curvilinear gradient assessment under spatial process models. J. Amer. Statist. Assoc. 2006;101:1487–1501. MR2279474. [PMC free article] [PubMed]
  • Banerjee S, Gelfand AE, Finley AO, Sang H. Gaussian predictive process models for large spatial datasets. J. Roy. Statist. Soc. Ser. B. 2008;70:825–848. [PMC free article] [PubMed]
  • Beneš V, Bodlák K, Møller J, Waagepetersen R. Research Report R-02-2001. Dept. Mathematical Sciences, Aalborg Univ.; Aalborg: 2002. Bayesian analysis of log Gaussian Cox processes for disease mapping.
  • Diggle PJ. Statistical Analysis of Spatial Point Patterns. 2nd ed. Arnold; London: 2003. MR0743593.
  • Diggle PJ, Rowlingson BS. A conditional approach to point process modelling of elevated risk. J. Roy. Statist. Soc. Ser. A. 1994;157:433–440.
  • Flood A, Peters U, Chatterjee N, Lacey JV, JR., Schairer C, Schatzkin A. Calcium from diet and supplements is associated with reduced risk of colorectal cancer in a prospective cohort of women. Cancer Epidemiology, Biomarkers & Prevention. 2005;14:126–132. [PubMed]
  • Flood A, Rastogi T, Wirfält E, Mitrou PN, Reedy J, Subar AF, Kipnis V, Mouw T, Hollenbeck AR, Leitzmann M, Schatzkin A. Dietary patterns as identified by factor analysis and colorectal cancer among middle-aged Americans. American Journal of Clinical Nutrition. 2008;88:176–184. [PMC free article] [PubMed]
  • Folsom AR, Hong C-P. Magnesium intake and reduced risk of colon cancer in a prospective study of women. American Journal of Epidemiology. 2005;163:232–235. [PubMed]
  • Fuchs C, Giovannucci E, Colditz G, Hunter D, Speizer F, Willett W. A prospective study of family history and the risk of colorectal cancer. New England Journal of Medicine. 1994;331:1669–1674. [PubMed]
  • Gelfand AE, Schmidt AM, Banerjee S, Sirmans CF. Nonstationary multivariate process modelling through spatially varying coregionalization (with discussion) Test. 2004;13:263–312. MR2154003.
  • Glebov OK, Rodriguez LM, Nakahara K, Jenkins J, Cliatt J, Humbyrd C-J, Denobile J, Soballe P, Simon R, Wright G, Lynch P, Patterson S, Lynch H, Gallinger S, Buchbinder A, Gordon G, Hawk E, Kirsch IR. Distinguishing right from left colon by the pattern of gene expression. Cancer Epidemiology, Biomarkers & Prevention. 2003;12:755–762. [PubMed]
  • Guan Y. A composite likelihood approach in fitting spatial point process models. J. Amer. Statist. Assoc. 2006;101:1502–1512. MR2279475.
  • Guan Y, Loh JM. A thinned block bootstrap variance estimation procedure for inhomogeneous spatial point patterns. J. Amer. Statist. Assoc. 2007;102:1377–1386. MR2412555.
  • Guan Y, Waagepetersen R, Beale C. Second-order analysis of inhomogeneous spatial point processes with proportional intensity functions. J. Amer. Statist. Assoc. 2008;103:769–777.
  • Huang CS, O'Brien MJ, Yang S, Farraye F. Hyperplastic polyps, serrated adenomas, and the serrated polyp neoplasia pathway. American Journal of Gastroenterology. 2004;99:2242–2255. [PubMed]
  • Kulldorff M. SaTScan user guide, Version 7.0. 2006 Available at http://www.satscan.org/
  • Liang S, Banerjee S, Carlin BP. Bayesian wombling for spatial point processes. Biometrics. 2009 To appear. [PMC free article] [PubMed]
  • Liang S, Carlin BP, Gelfand AE. Supplement to “Analysis of Minnesota colon and rectum cancer point patterns with spatial and nonspatial covariate information.” 2009 DOI: 10.1214/09-AOAS240SUPP. [PMC free article] [PubMed]
  • Møller J, Waagepetersen RP. An introduction to simulation-based inference for spatial point processes. In: Møller J, editor. Spatial Statistics and Computational Methods. Lecture Notes in Statistics. Vol. 173. Springer; New York: 2003. pp. 143–198. MR2001387.
  • Møller J, Waagepetersen RP. Statistical Inference and Simulation for Spatial Point Processes. Chapman and Hall/CRC Press; Boca Raton, FL: 2004. MR2004226.
  • O'Brien MJ, Yang S, Mack C, Xu H, Huang CS, Mulcahy E, Amorosino M, Farraye FA. Comparison of microsatellite instability, CpG island methylation phenotype, BRAF and KRAS status in serrated polyps and traditional adenomas indicates separate pathways to distinct colorectal carcinoma end points. American Journal of Surgical Pathology. 2006;30:1491–1501. [PubMed]
  • Pedersen A, Johansen C, Gronbaek M. Relations between amount and type of alcohol and colon and rectal cancer in a Danish population based cohort study. Gut. 2003;52:861–867. [PMC free article] [PubMed]
  • PHYSICAL ACTIVITY GUIDELINES ADVISORY COMMITTEE . Physical Activity Guidelines Advisory Committee Report, 2008. U.S. Dept. Health and Human Services; Washington, DC: 2008.
  • Waagepetersen R. An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics. 2007;63:252–258. MR2345595. [PubMed]
  • Waagepetersen R. Estimating functions for inhomogeneous spatial point processes with incomplete covariate data. Biometrika. 2008;95:351–363. MR2422695.
  • Waagepetersen R, Guan Y. Two-step estimation for inhomogeneous spatial point processes. J. Roy. Statist. Soc. Ser. B. 2008 To appear.
  • Wakefield J, Salway R. A statistical framework for ecological and aggregate studies. J. Roy. Statist. Soc. Ser. A. 2001;164:119–137. MR1819026.
  • Wei EK, Giovannucci E, Wu K, Rosner B, Fuchs CS, Willett WC, Colditz GA. Comparison of risk factors for colon and rectal cancer. International Journal of Cancer. 2003;108:433–442. [PMC free article] [PubMed]
  • Wolpert RL, Ickstadt K. Poisson/gamma random field models for spatial statistics. Biometrika. 1998;85:251–269. MR1649114.