### Data

Cases are from the Massachusetts Cancer Registry (MCR), and include all invasive breast cancer cases diagnosed in female residents between 1988 and 1997 (n = 46,333), using the standard definition of invasive breast cancer[18]. Annual reports on MCR data completeness, methods and other issues can be found at the cited website [19]. Each case record included the town, ZIP Code, and census tract of the patient's usual residence at the time of diagnosis, as well as the age at diagnosis, date of diagnosis, race, and stage of disease at diagnosis. Of the 46,333 patients, 43,529 were white females, 1071 black, 6 American Indian, Aleutian, or Eskimo, 115 Chinese, 17 Japanese, 129 other Asian and Pacific Islander, 119 other non-white, and 1347 of unknown race. Race was not included as a covariate in this study due to the large number of unknown. Age was included as a covariate and was divided into fourteen 5-year categories starting with age 15 years up through age 84, with all ages 85 years and above as the fifteenth age group.

### Aggregation unit

Census tracts were used to aggregate cases because they are more uniform than towns or Zip Codes in population size, provide a more sensitive analysis of densely populated areas, and are more homogeneous in their resident characteristics. However, 12.6% (n = 5832) of the cases diagnosed in 1988–1997 could not be assigned a reliable residential census tract because of inaccuracies or omissions in the address information provided to the MCR. In most of these cases, a mailing address had been provided and, even after extensive research, MCR staff could not assign a reliable residential address for the patient at the time of diagnosis. To assign the unassigned cases to census tracts, we compared town and census tract boundaries for 1196 four-digit census tracts. About 90 of the census tract codes included 2-digit suffixes to designate 2, 3, 4, or 5 distinct tracts. In this study, only the first four digits of such tract codes were used, thus combining some separate tracts into one. For example, census tracts 6002.01 and 6002.02 were combined into the 4-digit tract 6002. 1990 Census block numbering areas in Massachusetts are treated here as if they were census tracts. Nearly all of the 1196 four-digit 1990 tracts could be related to town boundaries, with one or more tracts located entirely within a town, or one or more towns entirely within a census tract. Two census tracts overlapped parts of multiple towns and for these it was possible to estimate the proportion of the total tract population living in each town. Unassigned cases were then assigned as follows: for towns completely contained in one tract, cases were assigned to that tract, accounting for 1392 cases. For a town containing two or more census tracts, cases were randomly assigned to tracts proportionate to the town's female population, accounting for 4440 cases. The allocation of cases should therefore be free from systematic error and any error should be localized to a particular town, while the broader state patterns remain correct.

### Statistical analyses

The spatial scan statistic [20] was used to perform purely spatial and space-time analyses. Under the null hypothesis, the incidence of breast cancer follows a Poisson distribution and the probability of a case being diagnosed in a particular location is proportional to the covariate-adjusted population in that location. Tract population data are from the 1990[21] and 2000 Decennial Censuses [22]. This population data is an excellent match to the MCR data since both the Census and the MCR ask patients for their "usual residence." Therefore the numerator and denominator are created in the same way.

Properties that make the spatial scan statistic suitable for geographic surveillance are: 1) it takes into account the uneven geographic distribution of cases and population densities, 2) it does not make assumptions about cluster size or location, 3) it adjusts for multiple testing – a common problem in testing multiple combinations of cluster locations and sizes, 4) it identifies the spatial or space-time locations where the null hypothesis is rejected, and 5) it can detect multiple clusters [10].

Purely spatial analyses were performed first, ignoring time. The maximum spatial cluster size was first set to include up to 50% of the population for both excesses and deficits and then set at 10%, to test for excesses and deficits separately because testing at the 10% level can identify smaller, more defined areas. However, each area identified in an analysis had a likelihood associated with it that was compared to the 9,999 likelihoods from the initial 50% maximum spatial size test, thus accounting the inference for a multiplicity of statistical tests. Space-time analyses were then performed to determine whether the clusters from the purely spatial analysis were long term or temporary. The maximum temporal cluster size was set at 90% and also included purely spatial clusters with a temporal size of 100% for all space-time analyses. This allows for the entire time period all the way down to the smallest amount of time to be considered as the time window, ten years down to one year.

All SaTScan analyses were performed on age-adjusted population counts. To calculate the age-adjusted expected counts, the 1990 and 2000 female population counts were combined into a weighted average of the two based on the years being analyzed, 1988 through 1997. This was done for each age group within each tract. The natural log of this weighted average was entered as the offset variable in a Poisson regression in SAS [23]. There were a few tracts with a zero population for a certain age group; for these, a population of one was entered so a log could be taken. The Poisson regression included age as a class variable and the number of cases within each tract and age group as the dependent variable to calculate the age-adjusted expected counts. The expected counts were aggregated across age for all SaTScan analyses. This method was verified by entering age as a covariate in a SaTScan analysis with the same results. Although SaTScan accurately estimates age-adjusted incidence rates, when other covariates enter the model, SaTScan computes the interaction of each of the 15 age categories with each covariate. In the current study, that would mean estimating 14 × 1 × 1 × 4 = 56 interaction terms because three risk factors were used as covariates. By first computing the age-adjusted counts in SAS, the number of interactions estimated by SaTScan is reduced to 8.

After identifying statistically significant geographic areas and determining whether they were long term or temporary, the next step was to determine whether these areas would change when the model was adjusted for known risk factors. Socioeconomic factors have long been identified as risk factors for breast cancer [9, 24, 25]. Urban/rural status has also been used as a covariate when studying breast cancer [4, 26–29].

An SES index was created following the method of Yost et al [30]. Yost and colleagues included the following variables from the U.S. Decennial Census in a principal component analysis (PCA): percent unemployed, percent working class, percent below the Federal poverty level, median income, median rent, median house value, and percent with at least a high school diploma. We used the same variables from the 1990 Census [21] in a PCA, described in the results section. A PCA factors the matrix of correlations between all pairs of variables into seven uncorrelated components with the hope that a smaller number of these components will account for 70 to 80% of the variation contained in the matrix of correlations among the original seven variables. Instead of using seven highly correlated SES variables as covariates, PCA can reduce the number of covariates to one or two without loss of information. A PCA was also performed on a similar group of variable from the 1990 Census based on what Krieger [24] used to create a SES index, however it did not account for as much of the variance as the variables used by Yost [30].

A variable indicating the urban or rural status of each tract was derived from data available on the Massachusetts Institute for Social and Economic Research (MISER) [31]. An urban area consists of one or more contiguous census blocks with a population density of at least 1,000 persons per square mile, while all other areas are designated as rural [32]. We assigned all tracts entirely within a town the same classification as that town, regardless of their individual populations. Tracts that included several towns were classified as rural since all towns within such tracts were also classified as rural. The capacity of SES and urban/rural status to predict incident cases within census tracts was tested in Poisson regression models, using the GENMOD procedure in SAS [23].

SaTScan [33] created files listing census tracts and the number of the high and low areas they belonged to. These files were brought into Maptitude [34] to create the figures.