Investigation on circular asymmetry of geographical distribution in cancer mortality of Hiroshima atomic bomb survivors based on risk maps: analysis of spatial survival data

While there is a considerable number of studies on the relationship between the risk of disease or death and direct exposure from the atomic bomb in Hiroshima, the risk for indirect exposure caused by residual radioactivity has not yet been fully evaluated. One of the reasons is that risk assessments have utilized estimated radiation doses, but that it is difficult to estimate indirect exposure. To evaluate risks for other causes, including indirect radiation exposure, as well as direct exposure, a statistical method is described here that evaluates risk with respect to individual location at the time of atomic bomb exposure instead of radiation dose. In addition, it is also considered to split the risks into separate risks due to direct exposure and other causes using radiation dose. The proposed method is applied to a cohort study of Hiroshima atomic bomb survivors. The resultant contour map suggests that the region west to the hypocenter has a higher risk compared to other areas. This in turn suggests that there exists an impact on risk that cannot be explained by direct exposure.


Introduction
The risk of disease or death caused by exposure to atomic bomb radiation has been evaluated using estimated radiation doses based on information concerning age, shielding conditions, and distance from the hypocenter under the assumption that the radiation dose decreases with increasing distance from the hypocenter (see, e.g., Preston et al. 2007;Matsuura et al. 1997). For details of the dosimetry system used, see for example the DS02 system (Cullings et al. 2006;Young and Kerr 2005). The corresponding risk analyses focused solely on the risk from direct exposure to the atomic bomb, while the risk from indirect exposure due to residual radioactivity has been not evaluated in previous analyses. This means that the geographical distribution of risk has been structurally restricted to concentric circles under the assumption that the influence of direct exposure essentially depends on the distance from the hypocenter. For example, Peterson et al. (1983) have fitted Cox's proportional hazard models to cancer mortality rates, to investigate circular asymmetry around the hypocenter in Hiroshima and Nagasaki. Gilbert and Ohara (1984) have analyzed data on acute symptoms. They divided the survivors in the Life Span Study (LSS) cohort, registered at the Radiation Effect Research Foundation (RERF), into eight groups according to the survivors' location at the time of atomic bomb exposure relative to the hypocenter and evaluated the relative risk of each octant compared with that for survivors in the octant of east-north-east direction. However, we consider their approach to be not enough to investigate circular asymmetry around the hypocenter, because they evaluated only relative risks for each octant with respect to the location at exposure relative to the hypocenter and did not consider heterogeneity of risk in each octant.
Recently, survivors suspected of having suffered from indirect exposure were reported by Kamada et al. (2006), Kamada and Kawakami (2008), and Tonda et al. (2008) through biological studies and statistical analyses of the incidence of leukemia among the survivors who entered Hiroshima City on August 6, 1945, after the explosion of the atomic bomb. Furthermore, several questionnaire surveys (Uda et al. 1953;Masuda 1989) showed that so-called Black Rain, which might have included radioactivity, fell around the western part of Hiroshima City and the northwest suburbs for several hours just after the explosion. Ohtaki (2011) demonstrated spatial-time distributions of Black Rain using a nonparametric smoothing method applied to data from a questionnaire survey conducted by Hiroshima City in 2008, of about 37,000 inhabitants of Hiroshima and its suburbs that might have experienced Black Rain.
In the present paper, a statistical method is applied to evaluate the risk with respect to individual location at exposure rather than dose and construct a ''risk map,'' that is, a map based on the risk evaluated by location, to visually grasp the geographical distribution of risk without structural restrictions. The risk map allows discussing possible effects of indirect exposure due to ''Black Rain'' and other radioactivity on risk of mortality.

Data
The database of atomic bomb survivors (ABS), registered at the Research Institute for Radiation and Medicine (RI-RBM) at Hiroshima University, was used in the present study. The ABS differs from the LSS of the RERF, because the ABS cohort includes examined survivors residing in Hiroshima Prefecture, and data on health status for survivors also have been cumulatively compiled in the database. The extent of overlap between survivors in the ABS and the LSS was examined by Hayakawa et al. (1994) and Hoshi et al. (1996). Hayakawa et al. (1994) showed that the dose estimates of the ABS were close to those of the LSS among the overlapped subjects. However, it has not been tested how they agree to DS02.
From the ABS, we chose 31,055 subjects for analysis who satisfied the following conditions: (i) being alive and recognized as an atomic bomb survivor as of January 1, 1980 and (ii) having coordinate information on location at the time of atomic bomb exposure (abbreviated in the following as ''location at exposure''). These subjects were followed until December 31, 1997. The endpoint is death from solid cancers (number of deaths: 2,545). Subjects were treated alive at the end of follow-up, in case migration and loss to follow-up for other reasons as censoring (number of subjects: 28,510). Mesh coordinates of 100 m in width were used to define location at exposure [Hoshi et al. (1996)]. Sex, age at atomic bomb exposure (abbreviated in the following as ''age at exposure''), and shielding condition were used as covariates. Estimations of radiation doses were based on Hoshi et al. (1996) and Matsuura et al. (1997). Figure 1 is the scatter plot of location at exposure with the hypocenter as the origin. Gray lines represent the map of Hiroshima city according to the town planning map made between 1925 and 1928. The vertical and horizontal scales are the coordinates in units of kilometers with the origin being the hypocenter. The red and green lines are boundary of heavy and light rainfall area of ''Black Rain'' based on Uda's questionnaire survey (Uda, et al. 1953). Note that upper left region of boundary is rainfall area.

Statistical analyses
Data containing information on location are called ''spatial data.'' Several methods for analyzing spatial data have been proposed, depending on the type of outcome. Geographically weighted regression (GWR), proposed by Fotheringham et al. (2002), corresponds to multiple linear regression analysis of spatial data. GWR is essentially repeated local multiple linear regressions applied to data in the neighborhood of a given location. The GWR approach can be extended to logistic regression for spatial binary data and Poisson regression for spatial count data, but the methodology for spatial survival data, such as those in the study of atomic bomb survivors, still remains to be developed. Recently, Tonda and coworkers (Tonda et al. 2010) proposed a statistical method for spatial data by extending a method proposed for longitudinal data and Satoh et al. 2009). Their approach is applicable not only to spatial continuous and discrete data but also to spatial survival data. In the present paper, a method is developed for estimating the geographical distribution of mortality risk for atomic bomb exposure by extending Cox's proportional hazards model for spatial survival data (Tonda et al. 2010); the resulting method is applied to a cohort study of Hiroshima atomic bomb survivors.

Hazard model with spatially varying coefficients
Consider the proportional hazards model with spatially varying coefficients, which allows the effect of covariates to vary with location. Let (u, v), r, t, sex, and atb denote location at exposure, registered age, attained age, gender (sex = 0 if male, sex = 0 if female), and age at exposure. The proportional hazards function with spatially varying coefficient is then given by where h 0 tjshielding ð Þ is the baseline hazard function dependent on the shielding condition, b l (u, v) is the spatially varying coefficient, and b s and b a are ordinary regression coefficients that are constant with regard to location at exposure. Note that Eq. 1 represents an ordinary Cox model if the spatially varying coefficients are replaced by constant coefficients. Therefore, Eq. 1 represents an extension of a Cox model, and the interpretation of coefficients in Eq. 1 is similar to that with a Cox model. In particular, exp (b l (u, v)) denotes the hazard ratio compared with the location as the reference.
It is assumed that the shape of (b l (u, v)) is in a class described by linear combinations of unknown parameters h and known basis functions We use a polynomial surface basis, which is commonly used in the field of spatial interpolation (Ripley 1981;Venables and Ripley 2002). For example, a quadratic polynomial surface basis is given by In addition, a circular surface basis is expressed by To obtain a smoother shape for the spatially varying coefficient, one can use, for example, a B-spline or a Gaussian basis. Details are given in Satoh et al. (2003), Ruppert et al. (2003), and Konishi and Kitagawa (2010).
. . .; ng, where d i denotes the indicator variable specifying whether subject i is censored or not at time t i , with 1 denoting a failure and 0 denoting censored, the unknown parameters h, b s , and b a can be estimated by maximizing the partial likelihood (Cox 1972(Cox , 1975: where I is the set of indices of failure cases, that is, . . .; ng, and R i is the set of indices of cases who are alive at t i , that is, is expressed byb l ðu; vÞ ¼ĥ 0 xðu; vÞ. Theoretical properties ofb l ðu; vÞ given in Tonda et al. (2010) allows to construct a confidence region for b l (u, v) and test hypotheses about the shape of b l (u, v). In particular, the test of the hypothesis is a test of spatial homogeneity. This test is meaningful because it verifies statistically whether b l (u, v) is spatially varying or not. If the hypothesis of spatial homogeneity is rejected, there exists a regional difference on b l (u, v). Any further discussion of the methodology for confidence regions and tests for b l (u, v) is beyond the scope of this paper; for additional information, see Tonda et al. (2010).

Dose effect model
If dose denotes radiation dose (Hoshi et al. 1996;Matsuura et al. 1997), the risk for direct exposure from atomic bomb radiation was modeled by adding the term b d (t, atb)9 dose to the hazard function in Eq. 1, that is, where h 0 (t) denotes the common baseline hazard. The coefficient b d (t, atb) depends on both of attained age and age at exposure. According to a mathematical model of carcinogenesis, such as the generalized Armitage-Doll model (Ohtaki et al. 1985;Pierce and Mendelsohn 1999;Ohtaki and Niwa 2001;Pierce and Vaeth 2003), we assume the functional structure where k d denotes the effect of radiation dose for survivors exposed at age 0 and k d denotes the influence of age at exposure on sensitivity to radiation dose (the derivation of Eq. 7 is given in the ''Discussion'' section). Equation (7) means that the risk for radiation varies with age at exposure and decreases with increasing attained age. The unknown parameters can be estimated with a slight modification of the partial likelihood in Eq. 4, because 1 þ b d ðt; atbÞ Â dose % exp b d ðt; atbÞ Â dose ð Þis valid for small doses. Note that the hazard function excluding b l (u, v) from Eq. 6 represents the ordinary Cox's hazard model with timedependent variables, given by Equation (8) was used for modeling the relationship between the mortality risk and risk from direct exposure in previous studies [see, e.g., Pierce and Vaeth (2003)].

Results
The proposed method was applied to data from a cohort study of Hiroshima atomic bomb survivors. The method is easy to implement using statistical packages that execute Cox model, such as SAS, SPSS, and R. We used the ''survival'' package version 2.36-2 in R version 2.12.0 [R Development Core Team (2010)]. Table 1 demonstrates the goodness-of-fit among the various models on a q-th order polynomial basis (q = 1, 2, 3, 4) and a circular basis to describe the shape of b l (u, v) in Eq. 1. Statistically comparing five models using Akaike's Information Criteria (AIC; Akaike 1973), we selected the most suitable basis among the five in the manner of variable selection [see Tableman and Kim (2004), chapter 5]. Table 1 suggests the quadratic polynomial model to be optimal. The martingale residuals [Klein and Moeschberger (2003), chapter 11; Therneau and Grambsch (2000), chapter 4] were used to check whether the quadratic polynomial basis adequately describes the geographical distribution of risk. Assessing the presence or absence of spatial trend of residuals by applying generalized additive models [Wood (2006) chapter 4], the hypothesis of spatial homogeneity for residuals in the quadratic polynomial model was not rejected. Therefore, the quadratic basis seems to describe the spatial trend of risk adequately. Figure 2 shows the estimated risk map of mortality based on the quadratic polynomial model, while Table 2 shows the estimated coefficients of b s and b a . In Fig. 2, the contours on the map represent the hazard ratio, exp b l ðu; vÞ ð Þ, for each location compared with the reference location, marked as the blue cross, that is 2 km from the hypocenter toward the east. From Fig. 2, it can be seen that the mortality risk decreases with increasing distance from the hypocenter, but the geographical distribution of the risk map is not concentric: The west area appears to have a higher risk compared with other areas.
In Fig. 3, the decreasing trend of risk with distance from the hypocenter by direction of location at exposure defined by angle from the hypocenter is compared. Angles 178°( about west direction) and 62°(about north-north-east direction) had highest and lowest relative risks, respectively. Figure 3 also suggests that the risk at 2 km from the hypocenter at angle 178°corresponds to the risk at 1,147 m at angle 62°. Figure 4 shows the differences of relative risks by angle from the hypocenter compared with those of angle 62°(about north-north-east direction). Figure 4 suggests that the differences of relative risks become larger with increasing distance from the hypocenter. Figure 5 shows the estimated survival curves at 2 km from the hypocenter at angles 178°and 62°for women with 10 years age at exposure.
Finally, we considered removing the risk for the direct exposure from the risk map in Fig. 2 using the dose effect model. The risk for direct exposure can be drawn as a function of radiation dose. In this case, analysis had to be restricted to those individuals for whom information on radiation dose is available, which slightly decreased the number of subjects. The resulting risk map for the doseeffect model with quadratic basis is given in Fig. 6. Table 3 shows the estimated coefficients except for b l (u, v). In Figs. 7 and 8, the estimated relative risks at 1 Gy with attained age for ages at exposure 10, 20, and 30 years based on Eq. 7 are shown. Figure 7 presents estimated relative risks with adjustment for location at exposure based on the hazard function in Eq. 6, with estimated coefficients given in Table 3, while Fig. 8 presents estimated relative risks without adjustment for location at exposure, based on the ordinary Cox model in Eq. 8, with estimated coefficients given in Table 4. Figure 9 presents the contour map based on estimated direct radiation dose averaged by location at exposure. From Fig. 9, it can be seen that the geographical distribution of direct radiation dose is close to concentric circles. This means that if the risk due to causes other than the direct exposure was negligible compared with that of direct exposure, then the contours in the risk map should be well approximated by concentric circles. If not, however, the risk contours should be far from concentric and circular. According to Table 1 and Figs. 2, 3, 4, and 5, the resultant risk map for a cohort of Hiroshima atomic bomb survivors (Fig. 2) suggests that the quadratic polynomial contours are suitable indeed, but not concentric circles. This suggests that there existed risk factors other than direct radiation exposure.   . 3 Comparison of decreasing trend of relative risks with distance from the hypocenter by angles of location at exposure. The red and blue curves denote the highest and lowest trend, whose angles are 1788 (about west direction) and 648 (about north-north-east direction) Fig. 4 Comparison of increasing trend of relative risks, relative to those at angle 64°(about north-north-east direction), with distance from the hypocenter by angles of location at exposure

Discussion
We also considered removing risks for direct exposure from the risk map in Fig. 2, in order to grasp the geographical distribution of risks for potential causes other than the direct exposure. This purpose was achieved by adding the term b l (t, atb)9 dose to the hazard function. The age dependence of the dose effect has been formulated under the assumption that tumorigenesis requires multi stages of cell variation and that any carcinogenic transitions of cell variation are sensitive to radiation exposure [for details, see Ohtaki and Niwa (2001)]. According to the generalized Armitage-Doll model (Ohtaki et al. 1985;Pierce and Mendelsohn 1999;Ohtaki and Niwa 2001;Pierce and Vaeth 2003), the hazard function in Eq. 1 is modified by where a d denotes the risk due to radiation exposure, a a the relative sensitivity varying with age at exposure, and k the number of mutations required for a normal cell to become malignant. For convenience, the following log-linearization was applied where k d ðk À 1Þa d and k a ðk À 1Þa d a a . Equation (7) is now derived by substituting Eq. 10 into Eq. 9. According to Fig. 6, the resultant risk map, excluding the risks for direct exposure, still has contours skewed toward the west direction. In addition, the test for the hypothesis on spatial homogeneity, formulated by Eq. 5, was rejected (p \ 0.001). These results may provide further evidence of risks for causes other than direct exposure.
As was mentioned in the introduction, several questionnaire surveys showed that Black Rain, which might have included radioactivity, fell around the western part of Hiroshima city and north-west suburbs for several hours just after the explosion. According to the latest results on the geographical distribution of Black Rain (Ohtaki 2011) and Uda's rainfall area described in Fig. 2, the area of rainfall appears roughly similar to the region of high risk in Fig. 2. This similarity suggests that Black Rain might be a possible risk factor accounting for the geographical   and environmental factors that are probably unrelated to radiation exposure due to the atomic bomb. These factors might correlate through association with particular regions, but this will be difficult check.
Note that Peterson et al. (1983) have also studied the circular asymmetry around the hypocenter in Hiroshima and Nagasaki for the LSS cohort of RERF. They divided the survivors into eight groups by the octants according the survivors' location at exposure and fitted a Cox's proportional hazard model. According to their results, the survivors in the west-north-west octant had the highest risk and the relative risk of survivors in the west-north-west compared with those in the east-north-east was about 1.24. As was mentioned in the introduction, their approach suffered from a lack of continuity of risks within groups and between groups. Therefore, they could not grasp any regional spatial trend of risk within and between octants. On the other hand, our results for the ABS cohort of Fig. 7 Estimated relative risks with adjustment for location at exposure based on the hazard function in Eq. 6. Each curve denotes the plot of exp b d ðt; atbÞ ð Þ , which gives relative risks (RR) at 1 Gy with attained age for ages at exposure 10, 20, and 30 years Fig. 8 Estimated relative risks without adjustment for location at exposure, based on the hazard function in Eq. 8. Each curve denotes the plot of exp b d ðt; atbÞ ð Þ , which gives relative risks (abbreviated ''RR'') at 1 Gy with attained age for ages at exposure 10, 20, and 30 years  RIRBM can be used to understand the spatial trend visually. Our result in Fig. 3 is roughly consistent with the areas with higher risks in Peterson et al. (1983). In addition, Fig. 4 shows that the differences in the relative risk among angles of location at exposure become larger with increasing distance from the hypocenter, while Peterson et al. (1983) could only evaluate the relative risks by octants. In this sense, our results are somewhat more valuable than those of Peterson et al. (1983).

Conclusion
The risk map shown in the present work can be interpreted in terms of the radiation dose required to explain the fitted contours for the hazard ratio in Fig. 6. For this, we focused on the hazard ratios at the locations with 2-km distance from the hypocenter in Fig. 6. The relative risk between highest and lowest risk at such locations is about 1.6. This suggests an excess relative risk (ERR) of about 0.6 due to causes other than direct exposure. This value might correspond to quite a large dose (i.e., of more than a Gray) if most of this additional risk is caused by external exposure not yet included in the estimated direct doses, because the ERR per Gray for solid cancer among atomic bomb survivors is on the order of about 0.5 (see, e.g., National Research Council 2006). This is quite unlikely as direct radiation doses where verified experimentally for example by retrospective thermoluminescence measurements on environmental samples (see, e.g., Cullings et al. 2006;Young and Kerr 2005). However, it might be possible that additionally chronic continuous exposure and individual variability caused by internal exposure, that is not included in the current (direct) dose estimates for the atomic bomb survivors, had a large effect on cancer mortality risk among atomic bomb survivors. Unfortunately, data on incorporated radionuclides from fallout are limited, and the effect of any internal exposure requires further clarification. Therefore, the doses corresponding to the contours of risk shown in Fig. 6 also should be an issue in the future. As already mentioned, there might be additional risk factors affecting mortality such as socioeconomic status, life style, and environmental factors that could also explain part of the observed asymmetry, but these factors are difficult to investigate due to limited data available.