- Email: [email protected]
A method for spatial analysis of risk in a population-based case-control study
International Journal of Hygiene and Environmental Health
Int. J. Hyg. Environ. Health 205, 115 ± 120 (2002) ¹ Urban & Fischer Verlag http: // www.urbanfischer.de/journals/intjhyg
A method for spatial analysis of risk in a population-based casecontrol study Vero¬nica Vieiraa, Thomas Webstera, Ann Aschengraub, David Ozonoffa a b
Department of Environmental Health, Boston University School of Public Health, Boston, MA, USA Department of Epidemiology and Biostatistics, Boston University School of Public Health, Boston, MA, USA
Received September 22, 2001 ¥ Accepted November 7, 2001
Abstract Spatial distributions of disease occurrence and risk have traditionally served as a tool for identifying exposures of public health concern. Current software for manipulating geographic data (GIS) now allows many kinds of analyses not feasible before. This paper presents a method for producing a ™picture∫ of disease risks based on residential history data from a population based case-control study. We illustrate the method using geographically coded data on individual-level risk factors, such as age and smoking, from a cancer study of the Upper Cape Cod region of Massachusetts for 1983 ± 1986. Our results show the association between lung cancer and residential location as an indicator of exposure. Crude and adjusted odds ratios were estimated by adaptive rate stabilization and mapped using kriging as an interpolation method to examine the risk of lung cancer in the region. The crude and adjusted surfaces for various smoothing parameters were compared to identify areas of increased lung cancer not explained by standard risk factors. Such spatial patterns of disease risk may provide clues to exposures of importance or confirm associations with previously suspected exposures. Key words: Disease mapping ± kriging ± smoothing ± lung cancer ± Cape Cod
Introduction In 1988, an unusually high incidence of cancers in the upper region of Cape Cod, Massachusetts (U. S. A.) prompted a series of epidemiological studies to investigate possible environmental risk factors (™exposures∫) associated with the region, among them proximity to the Massachusetts Military Reservation (MMR), pesticide applications to the cranberry bogs in the area, a large electric power plant, and tetrachloroethyelene-contaminated drinking water (Aschengrau et al. 1992, 1993,
1996, 1998; Ozonoff et al. 1994; Paulu et al. 1999). Positive associations were observed, but none completely explained the excess incidence of cancer. These studies also provide an invaluable data set for spatial analysis. Here we present methodology for depicting risk in the region as a function of location, using lung cancer for purposes of illustrating the method. The ultimate objective of identifying areas of higher risk is to produce new hypotheses that might explain why cancer rates are elevated in
Corresponding author: Thomas Webster, Dept. of Environmental Health, Boston University School of Public Health, Talbot 2E, 715 Albany St., Boston, MA 02118, USA. Phone: 617 638 4620, Fax: 617 638 4857, E-mail: [email protected]
1438-4639/02/205-115 $ 15.00/0
116
V. Vieira et al.
the upper Cape Cod region compared to the rest of the state. Results for breast cancer will be presented elsewhere (Paulu et al., submitted 2001; Kustron, paper in preparation 2001).
Materials and methods Study population Begun in response to community concern, the Upper Cape Cancer Incidence Study was a population-based casecontrol study designed to evaluate the association between nine different kinds of cancer (breast, lung, colorectal, liver, brain, pancreas, bladder, kidney cancers and leukemia) and environmental exposures. During the period 1983 ± 1986, the Massachusetts Cancer Registry recorded 243 incident cases of lung cancer among permanent residents of Barnstable, Bourne, Falmouth, Mashpee and Sandwich, collectively referred to as the Upper Cape (Figure 1). Controls were chosen to represent the underlying population that gave rise to the cases. Selection criteria required controls to be permanent residents of the same towns during 1983 ± 86, and were frequency matched to cases on age, gender, and vital status. Because many of the cases were elderly or deceased, three different sources of controls were used: (1) Random Digit Dialing (RDD) identified living controls less than 65 years of age; (2) Health Care Finance Agency lists (HCFA files for Medicare) for the population older than 65 years of age, where RDD is usually very inefficient; and (3) death certificates for controls who had died from 1983 onward. The residences of the resulting controls provide an estimate of the underlying population distribution. Subjects or their next-of-kin completed an extensive questionnaire, which provided information on demographics (age, sex, marital status, education), a 40-year residential history, potential confounders, and exposure to risk factors. For lung cancer, potential confounders
Fig. 1.
Location of study area.
included age, smoking of cigarettes, pipe, and cigars, living with a smoker, and occupational exposure to lung carcinogens (work involving arsenic, asbestos, chromium, coal tar pitch). ™Index years∫ were randomly assigned to controls in a distribution similar to that of case diagnosis years. An index year for controls was the counter-part to year of diagnosis for cases, and was needed to estimate length and time of exposure for controls in a fashion comparable to the cases. 1206 controls were used in the final analysis. For a detailed description of the methods, see Aschengrau et al. 1993. Geographical Information System (GIS) All subject residential addresses in the Upper Cape area over the 40-year period were included in the spatial analysis. Detailed records were constructed for location (latitude/longitude), length and calendar years of residency. All addresses where residency time began after diagnosis date for cases or index date for controls were excluded. The final analysis included 367 residential locations for cases and 1800 for controls, for a total of 2167 points. The original analyses of these data referenced each subject on a paper USGS 1 : 24,000 scale map using small paper dots marked with identification numbers. This allowed certain very crude types of geographic information to be used in epidemiological analyses (distance from a hazardous waste site, for example), but made more sophisticated analyses cumbersome or infeasible. The advent of affordable and usable Geographic Information Systems in recent years now enables easy and complex queries from geocoded databases, thus relaxing this constraint. We digitized the locations of the subject residences from the paper database into Universal Transverse Mercator (UTM) coordinates and linked them to the individual's epidemiologic record. Digitizing was done without knowledge of case status and the final data were scrutinized for accuracy and a sub-sample subjected to repeat digitization as a cross-check. Our digitized data and
Method for spatial risk analysis
117
previously collected information were compared to identify which residences were missed or duplicated in the digitization process. The GIS allows us to map the coordinates of subjects and link the location to additional individual and environmental information already available from a variety of sources. Figure 2 shows the distribution of lung cancer cases and controls in the study area. (To preserve confidentiality, this figure was created by randomly placing residences within a 1.2 km grid that includes the actual location. Actual locations were used in the analysis.) Adaptive rate stabilization (ARS) Because the precision of rate estimates depends, in part, on the size of the underlying population at risk, mapping small areas across a large region produces estimates of varying precision, depending upon the corresponding population density. This makes the visual comparison induced by a map sometimes misleading, as depicted extremes might merely represent regions of unusual instability of the estimates. To produce estimates that are more comparable in this sense, we devised a method we call adaptive rate stabilization (Paulu et al., submitted 2001). Adaptive rate stabilization produces maps of smoothed odds ratios estimated from case-control data. Each subject serves as the center of a circular subregion whose radius is expanded to include k controls. For example, when k 10, the radius of a circle is expanded from each subject (case or control) until it encompasses 10 controls, as shown in Figure 3. The case-control odds of the subregion equals the number of cases within the circle divided by k. The odds for the subregion is then divided by the odds for the entire region, yielding a crude odds ratio (we used the
Fig. 2.
Distribution of lung cancer cases and controls.
Fig. 3.
Example of adaptive rate stabilization method.
entire Upper Cape as the reference population). Thus, for k 10, a subregion with 3 cases has an odds of 0.3. Dividing by 367/1800, the case-control odds for the fivetown area, produces an odds ratio of 1.47, an estimated 47% increased rate of lung cancer in the circle area. Adaptive rate stabilization is essentially a type of nearest neighbor smoothing, with the degree of smoothing determined by the number of controls k in a circle. Using a constant number of controls for each denominator keeps the precision more uniform throughout the study area by taking population into account, i. e., smaller circles are
118
V. Vieira et al.
used where data points are dense, and larger circles where points are scarce. As k increases, the surface of odds ratios flattens. The ultimate limit, of course, would be to include in each circle all 1800 control locations (the entire population), producing a completely flat surface. At the other extreme is a circle so small it encompasses only the focal point, producing an odds ratio of either zero (if the focal point is a control) or infinity (if the focal point is a case), a ™totally spiky∫ surface. Intermediate sized circles (intermediate levels of k > 0) produce intermediate levels of smoothing. To control for confounding, an adjusted odds ratio was calculated for each point using multiple logistic regression with confounders included as covariates. Exposure was treated as a binary variable depending upon whether or not a subject lived within the circle (Paulu et al., submitted 2001). Since the goal of this paper is to illustrate the method, no induction period was assumed for the current analysis, although this can be done relatively easily. Calculations were repeated with k values of 10, 20, 30, 40, and 50. The resulting crude and adjusted odds ratios were plotted in histograms and q-q normal plots to verify normality, an assumption needed for kriging, the interpolation method used to draw the maps (Webster and Oliver, 1994). Kriging Adaptive rate stabilization calculates a crude or adjusted odds ratio at each subject residence location, but mapping methods require values between the available points as well. This is accomplished with kriging, a method that interpolates continuous surfaces using a sophisticated averaging of nearby points. Kriging is preferred to conventional methods such as inverse distance weighting because weights are assigned so that kriging variance is
Fig. 4.
minimized and estimates are unbiased (Webster and Oliver, 2001). Furthermore, weights are not sensitive to clustering of points due to population density. In kriging, a description of how values vary with increasing separation is used to describe the degree of spatial dependence. The degree of similarity of values of nearby points is used to determine how much weight to give in a weighted average of the points. When spatial dependency is present, a plot of the model (the semivariogram) shows lower variance (i. e., similar values) among points at small separation distances that increases with increasing distances. The distance where the variance ceases to increase is the distance beyond which spatial structure no longer exists between points. This distance serves as the size of the search window used for interpolation. The semivariogram model that best fits the data then determines the appropriate weights for kriging. We employed ™spherical semivariograms∫ to provide the weights used in modeling the kriged surface of crude and adjusted odds ratios for each smoothing parameter.
Results Maps of crude and adjusted odds ratios for lung cancer in the Upper Cape were produced for k values (number of controls) of 10, 20, 30, 40, and 50. Spatial variability, denoted by the range of odds ratios over the surface, decreased as the smoothing parameter k increased, as expected (Figure 4). Figure 5 shows surfaces of crude and adjusted lung cancer odds ratios for k 50. The crude and
Effect of increased smoothing on lung cancer risk map (k 50 smoother than k 10).
Method for spatial risk analysis
Fig. 5.
119
Effect of adjusting for lung cancer risk factors.
adjusted maps are relatively similar. In the crude odds ratios map, an area of elevated lung cancer risk is identified in the northwestern corner of the region. However, the area becomes of lower risk in the adjusted map while a region to the east becomes of higher risk. Other prominent regions of increased risk identified in the adjusted map correspond to urban areas of Falmouth and Barnstable.
Discussion While local disease mapping (™cluster∫) investigations are often demanded by concerned communities, they have not been popular with epidemiologists. Critics argue that they often combine unrelated diseases, apply arbitrary boundaries, contain insufficient numbers of subjects, and ignore population density, latency, and confounding. The combination of a population-based case-control study with adaptive rate stabilization addresses many of these limitations and improves upon the use of GIS as a tool for describing disease risk. Because of the large data set from the Upper Cape study, a sufficient number of lung cancer cases allowed for spatial analysis without having to include other cancer types. There was no need for arbitrary boundaries because all residences were geocoded in a way independent of jurisdictional boundaries. While not a focus of the current study, the availability of a 40year residential history makes investigation of
latency effects and time-series analysis possible as well. In a crude analysis, the elevated rates in some areas may be due to geographic concentration of known risk factors, i. e., ™spatial confounding.∫ For example, lung cancer might be elevated due to geographical concentration of smokers. Figure 5 adjusts for several known lung cancer risk factors, addressing this problem. The region of elevated lung cancer risk in the northwestern part of the study area is particularly interesting. It does not disappear after adjustment for known risk factors (although its location shifts slightly). In an earlier non-spatial analysis (Ozonoff et al., 1994), the authors found a modest but unstable increase in lung cancer risk among subjects who ever lived within 3 kilometers of gun and mortar training sites at the Massachusetts Military Reservation (Compare Figures 1 and 5). The MMR is a large area in the middle of the study region where excess propellant was burned after artillery practice, leading to potential exposure of the surrounding population to airborne carcinogens. The map locates the position of the hotspot that may account for this association (no conclusions should be drawn until latency is considered). Our map suggests an analysis based on more detailed exposure information of residents in this location and a comparison group. The approach used here has some limitations. Ordinary kriging assumes that variance does not change with direction (isotropy). However, analysis of the semivariograms showed some directional
120
V. Vieira et al.
variation in the odds ratios for lung cancer. This suggests that a global (™nondirectional∫) variogram, commonly used in spatial epidemiology, is not the most appropriate model. Rather, a local variogram may be required (Kustron, paper in preparation 2001). Our current method also does not give any measure of the precision of the odds ratios, such as surfaces displaying confidence intervals. In addition, the reference populations in the crude and adjusted analyses are not precisely the same. For the crude OR map, we use the entire Upper Cape Cod region; the adjusted OR map does not include those who fall within the circle. With over 1800 control locations, a difference of 30 ± 50 should not be noticeable in the odds ratios. To confirm this, we carried out a series of logistic regressions with no covariates where k 30, 40, 50. The resulting crude maps showed no difference from those constructed using ARS. Hence, the comparison made in Figure 5 reflects the effect of adjustment, not a shifting reference population. We are currently exploring the use of generalized additive models (GAMs) to simultaneously perform spatial smoothing, adjustment for covariates, and generation of estimates of precision (Kelsall and Diggle, 1998). GIS is a tool that not only allows visualization but also permits the use of sophisticated analytical methods for investigating the role of environmental exposures in disease risk. By overlaying the odds ratio surface with maps of environmental exposures, new causal hypotheses may be suggested or suspected ones confirmed. Of particular interest are the combined effects of two or more exposures on human health. All possible combinations of environmental exposures are too numerous to test statistically, but maps can potentially reveal a concurrence of joint exposures and disease risk that can be investigated further.
Acknowledgements. This work was supported by the Superfund Basic Research Grant 1P42ES 07381.
References Aschengrau, A., Ozonoff, D.: Upper Cape Cancer Incidence Study. Final Report. Massachusetts Dept. of Public Health, Boston 1992. Aschengrau, A., Ozonoff, D., Paulu, C., Coogan, P., Vezina, R., Heeren, T., Zhang, Y.: Cancer risk and tetrachloroethylene-contaminated drinking water in Massachusetts. Arch. Environ. Health 48, 284 ± 292 (1993). Aschengrau, A., Ozonoff, D., Coogan, P., Vezina, R., Heeren, T., Zhang, Y.: Cancer risk and residential proximity to cranberry cultivation in Massachusetts. Amer. J. Public Health 86, 1289 ± 1296 (1996). Aschengrau, A., Paulu, C., Ozonoff, D.: Tetrachloroethylene-contaminated drinking water and the risk of breast cancer. Environ. Health Perspect. 106, 947 ± 953 (1998). Kelsall, J., Diggle, P.: Spatial variation in risk of disease: a nonparametric binary regression approach. Appl. Statist. 47, 559 ± 573 (1998). Ozonoff, D., Aschengrau, A., Coogan, P.: Cancer in the vicinity of a Department of Defense superfund site in Massachusetts. Toxicol. Indust. Health 10, 119 ± 141 (1994). Paulu, C., Aschengrau, A., Ozonoff, D.: Tetrachloroethylene-contaminated drinking water and the risk of colon-rectum, lung, and other cancers. Environ. Health Perspect. 107, 265 ± 271 (1999). Webster, R., Oliver, M.: Kriging the local risk of a rare disease from a register of diagnoses. Geographical Analysis. 26, 168 ± 184 (1994). Webster, R., Oliver, M.: Geostatistics for Environmental Scientists. John Wiley & Sons, Ltd., New York 2001.
Int. J. Hyg. Environ. Health 205, 115 ± 120 (2002) ¹ Urban & Fischer Verlag http: // www.urbanfischer.de/journals/intjhyg
A method for spatial analysis of risk in a population-based casecontrol study Vero¬nica Vieiraa, Thomas Webstera, Ann Aschengraub, David Ozonoffa a b
Department of Environmental Health, Boston University School of Public Health, Boston, MA, USA Department of Epidemiology and Biostatistics, Boston University School of Public Health, Boston, MA, USA
Received September 22, 2001 ¥ Accepted November 7, 2001
Abstract Spatial distributions of disease occurrence and risk have traditionally served as a tool for identifying exposures of public health concern. Current software for manipulating geographic data (GIS) now allows many kinds of analyses not feasible before. This paper presents a method for producing a ™picture∫ of disease risks based on residential history data from a population based case-control study. We illustrate the method using geographically coded data on individual-level risk factors, such as age and smoking, from a cancer study of the Upper Cape Cod region of Massachusetts for 1983 ± 1986. Our results show the association between lung cancer and residential location as an indicator of exposure. Crude and adjusted odds ratios were estimated by adaptive rate stabilization and mapped using kriging as an interpolation method to examine the risk of lung cancer in the region. The crude and adjusted surfaces for various smoothing parameters were compared to identify areas of increased lung cancer not explained by standard risk factors. Such spatial patterns of disease risk may provide clues to exposures of importance or confirm associations with previously suspected exposures. Key words: Disease mapping ± kriging ± smoothing ± lung cancer ± Cape Cod
Introduction In 1988, an unusually high incidence of cancers in the upper region of Cape Cod, Massachusetts (U. S. A.) prompted a series of epidemiological studies to investigate possible environmental risk factors (™exposures∫) associated with the region, among them proximity to the Massachusetts Military Reservation (MMR), pesticide applications to the cranberry bogs in the area, a large electric power plant, and tetrachloroethyelene-contaminated drinking water (Aschengrau et al. 1992, 1993,
1996, 1998; Ozonoff et al. 1994; Paulu et al. 1999). Positive associations were observed, but none completely explained the excess incidence of cancer. These studies also provide an invaluable data set for spatial analysis. Here we present methodology for depicting risk in the region as a function of location, using lung cancer for purposes of illustrating the method. The ultimate objective of identifying areas of higher risk is to produce new hypotheses that might explain why cancer rates are elevated in
Corresponding author: Thomas Webster, Dept. of Environmental Health, Boston University School of Public Health, Talbot 2E, 715 Albany St., Boston, MA 02118, USA. Phone: 617 638 4620, Fax: 617 638 4857, E-mail: [email protected]
1438-4639/02/205-115 $ 15.00/0
116
V. Vieira et al.
the upper Cape Cod region compared to the rest of the state. Results for breast cancer will be presented elsewhere (Paulu et al., submitted 2001; Kustron, paper in preparation 2001).
Materials and methods Study population Begun in response to community concern, the Upper Cape Cancer Incidence Study was a population-based casecontrol study designed to evaluate the association between nine different kinds of cancer (breast, lung, colorectal, liver, brain, pancreas, bladder, kidney cancers and leukemia) and environmental exposures. During the period 1983 ± 1986, the Massachusetts Cancer Registry recorded 243 incident cases of lung cancer among permanent residents of Barnstable, Bourne, Falmouth, Mashpee and Sandwich, collectively referred to as the Upper Cape (Figure 1). Controls were chosen to represent the underlying population that gave rise to the cases. Selection criteria required controls to be permanent residents of the same towns during 1983 ± 86, and were frequency matched to cases on age, gender, and vital status. Because many of the cases were elderly or deceased, three different sources of controls were used: (1) Random Digit Dialing (RDD) identified living controls less than 65 years of age; (2) Health Care Finance Agency lists (HCFA files for Medicare) for the population older than 65 years of age, where RDD is usually very inefficient; and (3) death certificates for controls who had died from 1983 onward. The residences of the resulting controls provide an estimate of the underlying population distribution. Subjects or their next-of-kin completed an extensive questionnaire, which provided information on demographics (age, sex, marital status, education), a 40-year residential history, potential confounders, and exposure to risk factors. For lung cancer, potential confounders
Fig. 1.
Location of study area.
included age, smoking of cigarettes, pipe, and cigars, living with a smoker, and occupational exposure to lung carcinogens (work involving arsenic, asbestos, chromium, coal tar pitch). ™Index years∫ were randomly assigned to controls in a distribution similar to that of case diagnosis years. An index year for controls was the counter-part to year of diagnosis for cases, and was needed to estimate length and time of exposure for controls in a fashion comparable to the cases. 1206 controls were used in the final analysis. For a detailed description of the methods, see Aschengrau et al. 1993. Geographical Information System (GIS) All subject residential addresses in the Upper Cape area over the 40-year period were included in the spatial analysis. Detailed records were constructed for location (latitude/longitude), length and calendar years of residency. All addresses where residency time began after diagnosis date for cases or index date for controls were excluded. The final analysis included 367 residential locations for cases and 1800 for controls, for a total of 2167 points. The original analyses of these data referenced each subject on a paper USGS 1 : 24,000 scale map using small paper dots marked with identification numbers. This allowed certain very crude types of geographic information to be used in epidemiological analyses (distance from a hazardous waste site, for example), but made more sophisticated analyses cumbersome or infeasible. The advent of affordable and usable Geographic Information Systems in recent years now enables easy and complex queries from geocoded databases, thus relaxing this constraint. We digitized the locations of the subject residences from the paper database into Universal Transverse Mercator (UTM) coordinates and linked them to the individual's epidemiologic record. Digitizing was done without knowledge of case status and the final data were scrutinized for accuracy and a sub-sample subjected to repeat digitization as a cross-check. Our digitized data and
Method for spatial risk analysis
117
previously collected information were compared to identify which residences were missed or duplicated in the digitization process. The GIS allows us to map the coordinates of subjects and link the location to additional individual and environmental information already available from a variety of sources. Figure 2 shows the distribution of lung cancer cases and controls in the study area. (To preserve confidentiality, this figure was created by randomly placing residences within a 1.2 km grid that includes the actual location. Actual locations were used in the analysis.) Adaptive rate stabilization (ARS) Because the precision of rate estimates depends, in part, on the size of the underlying population at risk, mapping small areas across a large region produces estimates of varying precision, depending upon the corresponding population density. This makes the visual comparison induced by a map sometimes misleading, as depicted extremes might merely represent regions of unusual instability of the estimates. To produce estimates that are more comparable in this sense, we devised a method we call adaptive rate stabilization (Paulu et al., submitted 2001). Adaptive rate stabilization produces maps of smoothed odds ratios estimated from case-control data. Each subject serves as the center of a circular subregion whose radius is expanded to include k controls. For example, when k 10, the radius of a circle is expanded from each subject (case or control) until it encompasses 10 controls, as shown in Figure 3. The case-control odds of the subregion equals the number of cases within the circle divided by k. The odds for the subregion is then divided by the odds for the entire region, yielding a crude odds ratio (we used the
Fig. 2.
Distribution of lung cancer cases and controls.
Fig. 3.
Example of adaptive rate stabilization method.
entire Upper Cape as the reference population). Thus, for k 10, a subregion with 3 cases has an odds of 0.3. Dividing by 367/1800, the case-control odds for the fivetown area, produces an odds ratio of 1.47, an estimated 47% increased rate of lung cancer in the circle area. Adaptive rate stabilization is essentially a type of nearest neighbor smoothing, with the degree of smoothing determined by the number of controls k in a circle. Using a constant number of controls for each denominator keeps the precision more uniform throughout the study area by taking population into account, i. e., smaller circles are
118
V. Vieira et al.
used where data points are dense, and larger circles where points are scarce. As k increases, the surface of odds ratios flattens. The ultimate limit, of course, would be to include in each circle all 1800 control locations (the entire population), producing a completely flat surface. At the other extreme is a circle so small it encompasses only the focal point, producing an odds ratio of either zero (if the focal point is a control) or infinity (if the focal point is a case), a ™totally spiky∫ surface. Intermediate sized circles (intermediate levels of k > 0) produce intermediate levels of smoothing. To control for confounding, an adjusted odds ratio was calculated for each point using multiple logistic regression with confounders included as covariates. Exposure was treated as a binary variable depending upon whether or not a subject lived within the circle (Paulu et al., submitted 2001). Since the goal of this paper is to illustrate the method, no induction period was assumed for the current analysis, although this can be done relatively easily. Calculations were repeated with k values of 10, 20, 30, 40, and 50. The resulting crude and adjusted odds ratios were plotted in histograms and q-q normal plots to verify normality, an assumption needed for kriging, the interpolation method used to draw the maps (Webster and Oliver, 1994). Kriging Adaptive rate stabilization calculates a crude or adjusted odds ratio at each subject residence location, but mapping methods require values between the available points as well. This is accomplished with kriging, a method that interpolates continuous surfaces using a sophisticated averaging of nearby points. Kriging is preferred to conventional methods such as inverse distance weighting because weights are assigned so that kriging variance is
Fig. 4.
minimized and estimates are unbiased (Webster and Oliver, 2001). Furthermore, weights are not sensitive to clustering of points due to population density. In kriging, a description of how values vary with increasing separation is used to describe the degree of spatial dependence. The degree of similarity of values of nearby points is used to determine how much weight to give in a weighted average of the points. When spatial dependency is present, a plot of the model (the semivariogram) shows lower variance (i. e., similar values) among points at small separation distances that increases with increasing distances. The distance where the variance ceases to increase is the distance beyond which spatial structure no longer exists between points. This distance serves as the size of the search window used for interpolation. The semivariogram model that best fits the data then determines the appropriate weights for kriging. We employed ™spherical semivariograms∫ to provide the weights used in modeling the kriged surface of crude and adjusted odds ratios for each smoothing parameter.
Results Maps of crude and adjusted odds ratios for lung cancer in the Upper Cape were produced for k values (number of controls) of 10, 20, 30, 40, and 50. Spatial variability, denoted by the range of odds ratios over the surface, decreased as the smoothing parameter k increased, as expected (Figure 4). Figure 5 shows surfaces of crude and adjusted lung cancer odds ratios for k 50. The crude and
Effect of increased smoothing on lung cancer risk map (k 50 smoother than k 10).
Method for spatial risk analysis
Fig. 5.
119
Effect of adjusting for lung cancer risk factors.
adjusted maps are relatively similar. In the crude odds ratios map, an area of elevated lung cancer risk is identified in the northwestern corner of the region. However, the area becomes of lower risk in the adjusted map while a region to the east becomes of higher risk. Other prominent regions of increased risk identified in the adjusted map correspond to urban areas of Falmouth and Barnstable.
Discussion While local disease mapping (™cluster∫) investigations are often demanded by concerned communities, they have not been popular with epidemiologists. Critics argue that they often combine unrelated diseases, apply arbitrary boundaries, contain insufficient numbers of subjects, and ignore population density, latency, and confounding. The combination of a population-based case-control study with adaptive rate stabilization addresses many of these limitations and improves upon the use of GIS as a tool for describing disease risk. Because of the large data set from the Upper Cape study, a sufficient number of lung cancer cases allowed for spatial analysis without having to include other cancer types. There was no need for arbitrary boundaries because all residences were geocoded in a way independent of jurisdictional boundaries. While not a focus of the current study, the availability of a 40year residential history makes investigation of
latency effects and time-series analysis possible as well. In a crude analysis, the elevated rates in some areas may be due to geographic concentration of known risk factors, i. e., ™spatial confounding.∫ For example, lung cancer might be elevated due to geographical concentration of smokers. Figure 5 adjusts for several known lung cancer risk factors, addressing this problem. The region of elevated lung cancer risk in the northwestern part of the study area is particularly interesting. It does not disappear after adjustment for known risk factors (although its location shifts slightly). In an earlier non-spatial analysis (Ozonoff et al., 1994), the authors found a modest but unstable increase in lung cancer risk among subjects who ever lived within 3 kilometers of gun and mortar training sites at the Massachusetts Military Reservation (Compare Figures 1 and 5). The MMR is a large area in the middle of the study region where excess propellant was burned after artillery practice, leading to potential exposure of the surrounding population to airborne carcinogens. The map locates the position of the hotspot that may account for this association (no conclusions should be drawn until latency is considered). Our map suggests an analysis based on more detailed exposure information of residents in this location and a comparison group. The approach used here has some limitations. Ordinary kriging assumes that variance does not change with direction (isotropy). However, analysis of the semivariograms showed some directional
120
V. Vieira et al.
variation in the odds ratios for lung cancer. This suggests that a global (™nondirectional∫) variogram, commonly used in spatial epidemiology, is not the most appropriate model. Rather, a local variogram may be required (Kustron, paper in preparation 2001). Our current method also does not give any measure of the precision of the odds ratios, such as surfaces displaying confidence intervals. In addition, the reference populations in the crude and adjusted analyses are not precisely the same. For the crude OR map, we use the entire Upper Cape Cod region; the adjusted OR map does not include those who fall within the circle. With over 1800 control locations, a difference of 30 ± 50 should not be noticeable in the odds ratios. To confirm this, we carried out a series of logistic regressions with no covariates where k 30, 40, 50. The resulting crude maps showed no difference from those constructed using ARS. Hence, the comparison made in Figure 5 reflects the effect of adjustment, not a shifting reference population. We are currently exploring the use of generalized additive models (GAMs) to simultaneously perform spatial smoothing, adjustment for covariates, and generation of estimates of precision (Kelsall and Diggle, 1998). GIS is a tool that not only allows visualization but also permits the use of sophisticated analytical methods for investigating the role of environmental exposures in disease risk. By overlaying the odds ratio surface with maps of environmental exposures, new causal hypotheses may be suggested or suspected ones confirmed. Of particular interest are the combined effects of two or more exposures on human health. All possible combinations of environmental exposures are too numerous to test statistically, but maps can potentially reveal a concurrence of joint exposures and disease risk that can be investigated further.
Acknowledgements. This work was supported by the Superfund Basic Research Grant 1P42ES 07381.
References Aschengrau, A., Ozonoff, D.: Upper Cape Cancer Incidence Study. Final Report. Massachusetts Dept. of Public Health, Boston 1992. Aschengrau, A., Ozonoff, D., Paulu, C., Coogan, P., Vezina, R., Heeren, T., Zhang, Y.: Cancer risk and tetrachloroethylene-contaminated drinking water in Massachusetts. Arch. Environ. Health 48, 284 ± 292 (1993). Aschengrau, A., Ozonoff, D., Coogan, P., Vezina, R., Heeren, T., Zhang, Y.: Cancer risk and residential proximity to cranberry cultivation in Massachusetts. Amer. J. Public Health 86, 1289 ± 1296 (1996). Aschengrau, A., Paulu, C., Ozonoff, D.: Tetrachloroethylene-contaminated drinking water and the risk of breast cancer. Environ. Health Perspect. 106, 947 ± 953 (1998). Kelsall, J., Diggle, P.: Spatial variation in risk of disease: a nonparametric binary regression approach. Appl. Statist. 47, 559 ± 573 (1998). Ozonoff, D., Aschengrau, A., Coogan, P.: Cancer in the vicinity of a Department of Defense superfund site in Massachusetts. Toxicol. Indust. Health 10, 119 ± 141 (1994). Paulu, C., Aschengrau, A., Ozonoff, D.: Tetrachloroethylene-contaminated drinking water and the risk of colon-rectum, lung, and other cancers. Environ. Health Perspect. 107, 265 ± 271 (1999). Webster, R., Oliver, M.: Kriging the local risk of a rare disease from a register of diagnoses. Geographical Analysis. 26, 168 ± 184 (1994). Webster, R., Oliver, M.: Geostatistics for Environmental Scientists. John Wiley & Sons, Ltd., New York 2001.