Study subjects

A total of 367,649 children and adolescents (<19 years of age at the time of the accident) were registered as the target population in the first round of TUE according to the latest summary report of the program28. The 59 municipalities at the primary examination were allocated by fiscal year: 2011 for the evacuation zone comprising 13 municipalities near the FNPP, 2012 for 12 municipalities in the middle part of the prefecture, and 2013 for the remaining regions including the most inland part, Aizu, and coastal regions except the evacuation zone (34 municipalities). Since there were late examinees, the first-round examination included 300,473 persons who underwent primary examination from October 2011 to April 2015. When the primary examination using ultrasonography found nodules (size 5.1 mm or more) or cysts (size 20.1 mm or more) or identified the urgent need for confirmatory testing due to clinical reasons, confirmatory examinations were recommended. Confirmatory testing included a detailed ultrasound, blood testing, urine analysis, and fine-needle aspiration cytology (FNAC) if needed. Details of the survey protocol are provided elsewhere16.

The Fukushima Medical University carried out the ultrasound examination programme and compiled the surveillance database which was the data source for this study. The dataset as of 27 February 2017 included 299,908 persons who underwent the primary examination. We excluded 1,523 persons who did not live in Fukushima Prefecture when the earthquake occurred on 11 March 2011. It should be noted that the excluded population included one thyroid cancer case. We further excluded 3,353 persons who did not provide information on residential municipality on the day of the earthquake. Consequently, data on the 295,032 subjects who underwent primary examination were analysed in this study. Among 2,246 examinees who tested positive in the primary examination, 2,051 (91.4%) underwent confirmatory testing. There were 115 cases that were classified as malignant or highly suspicious of malignancy based on FNAC at confirmatory testing. We regarded these cases as diagnosed thyroid cancer cases for this study. These cases were aggregated by municipality for this analysis. According to the summary report of TUE28, the malignancy rate after surgery for the FNAC-based diagnosed thyroid cancer cases was 99.0% for the first round TUE (102 FNAC-based diagnosed cases received surgery, and of those, 101 cases were post-surgically confirmed with thyroid carcinoma).

This study was approved by the Ethics Committee of the Fukushima Medical University (approval no. 1318). Written informed consent for the study was obtained from the parents of every participant. All methods were carried out in accordance with the applicable guidelines and regulations for the use and analysis of the Fukushima Health Management Survey data, managed by the Radiation Medical Science Center, Fukushima Medical University.

Spatial statistical analysis

Since age and sex are the most established determinants of thyroid cancer32, we used the sex and age-standardised prevalence ratio at the municipality level as follows:

$$sp{r}_{i}={o}_{i}/{e}_{i},$$

where i is the index of the municipality, o i is the number of observed cases of thyroid cancer and e i is the sex- and age-adjusted expected number of cases, also defined as:

$${e}_{i}=\sum _{k}po{p}_{ik}p{r}_{k},$$

where k is the index of sex (male and female) and 5-year age group (0–4, 5–9, 10–14, 15–19, and 20 or older), pop ik is the number of primary ultrasound examinees of the kth group in the ith municipality, and pr k is the prefectural marginal prevalence rate of the kth group. The age groups were based on the age of the examinees when they underwent the primary examination. We assume that the geographic distributions of observed cases, {o i }, were independently and randomly generated by the Poisson distribution with the expected number of cases, {e i }:

$${o}_{i} \sim Poisson[{r}_{i}{e}_{i}],$$

where r i is the relative risk of thyroid cancer in the ith municipality. The null hypothesis of our spatial analysis is the situation of geographically homogeneous risk (r i = 1 for all i). In other words, the hypothesis means that the prevalence rates of thyroid cancer are only determined by the age and sex distribution of examinees.

Method 1: detecting locally elevated risk regions

An alternative hypothesis based on the above null hypothesis is that there were regions with higher risks compared to the others. Kulldorff15 proposed an algorithm using a circular scan window which exhaustively searches for elevated risk regions among the set of all possible centre locations and sizes of circle window. The method uses likelihood ratio statistics by comparing risks inside and outside a circular window region. The window region with the highest likelihood ratio is called the ‘most likely cluster’. A Monte Carlo simulation under the assumption of homogeneous risk provides the null distribution of the statistics by which we may test the significance of the most likely cluster. Since the test statistic is defined by only one value in the entire study region, this method avoids the problem of multiple testing.

Tango and Takahashi21 proposed Flexscan to search for irregularly-shaped clusters with elevated risk. In this method, scanning windows are defined as a set of topologically connected spatial units (sharing a border point or line). If the true cluster is non-circular (e.g. linear), the modified spatial scan statistic is likely to have stronger power compared to the traditional circular scanning. Radiation dose is distributed in a non-circular form because radioactive plumes from the FNPP were carried by the wind; areas with high doses stretched from the FNPP in a north-westerly direction and then moved southward in the most northern part of the prefecture. Thus, if the distribution of standardised prevalence ratio is associated with radiation exposure, we expect that the flexibly shaped scan statistic is more likely to detect the elevated risk regions compared to the traditional circular scanning. The results may highlight specific regions with high risk on the map. For the computation, we used Flexscan version 3.1.2. For the technical details on this method, see Supplementary Note.

Method 2: detecting general clustering tendency

A tendency of clustering of prevalence can be considered as a positive spatial autocorrelation of risks. Detection of such a tendency may indicate that the risk of thyroid cancer has a spatial structure reflecting the distributions of some geographic factors. Tango’s MEET22 is considered the method with the highest power for detecting general tendency of spatial clustering33.

The test has two phases. Firstly, it measures the C-index of clustering tendency which requires a scale parameter, which determines “closeness” in geographic space in advance. Although the C-index can be used for statistical testing, the result depends on the choice of scale parameter. Searching for the optimal (most probable) scale of spatial clustering by testing different scale parameter series results in the multiple comparisons/testing. The second phase of MEET involves computing the adjusted p-value of clustering tendency with a predefined set of possible scale parameters and Monte Carlo simulation evaluation, by considering the statistical uncertainty of optimal scale selection. For the technical details, see Supplementary Note. We ran the R script of MEET26 in the R version 3.4 environment (The R Core Development Team, Vienna, Austria).

Method 3: Poisson regression

Ecological regression is another approach to assessing whether regional relative risks are geographically structured, particularly for those which have geographical indicators of suspicious factors of diseases under investigation. Poisson regression is a common choice based on the same assumption of Poisson process for disease clustering tests. A univariate Poisson regression using the log link function is:

$${o}_{i}\sim Poisson[{r}_{i}{e}_{i}],$$

$${\rm{l}}{\rm{n}}({r}_{i})={\beta }_{0}+{\beta }_{1}{x}_{i},$$

where x i is the explanatory variable in the ith municipality. We considered seven municipal indicators as explanatory variables: distance from the FNPP, the estimated external radiation dose, altitude, and several census-based socio-economic indicators.

Using focused tests which are designed to evaluate clustering around a pre-fixed geographic object, we selected the distance from the FNPP and the proportion of estimated external radiation doses of equal to or more than 1 mSv during the first four months after the accident at the FNPP. The distance variable was measured as the Euclidean distance from the FNPP to the town hall of each municipality. The information on estimated effective exposure to radiation from the FNPP was based on individual behavioural data that were collected as part of the Fukushima Health Management Survey. Following the methods of Ohira et al.12,13 who used the same information for defining regions with different radiation dose levels, we defined the variable, external radiation exposure, as a surrogate collective indicator of thyroid equivalent dose.

Related to iodine deficiency is altitude, which is a known geographic risk indicator for thyroid cancer, as it is an indicator of low access to iodine-rich food, such as fish and seaweed18. We employed the altitude of the town hall locations as the representative value for each municipality.

Furthermore, we included 2010 census-based indicators. In several previous studies, rural-urban indices and areal SES indicators were prepared to capture the general characteristics of thyroid cancer geographic distribution19,20,34,35. Since there were no established rural-urban indicators in Japan, we included population density as a simple surrogate of urbanicity36, and the proportion of agriculture, forestry, and fisheries industries workers as a rurality indicator. Regarding areal SES, the rate of unemployment in the entire labour force was included as the most influential component of areal deprivation36. Finally, the proportion of professional workers was used as an index of affluence.

For every explanatory variable, we applied the untransformed values as well as the quartile categories for examining possible non-linear relationships. For the computation of Poisson regression, we used the “glm” function in the R version 3.4 environment. To measure the distance and the cartographic mapping, we used ArcMap 10.5 (Environmental Systems Research Institute, Redlands, CA, USA). Euclidean distances were measured on the basis of map coordinates under rectangular plane (zone 9) projection with the Japanese Geodetic Datum 2011.

Sensitivity analysis of the undiagnosed positive cases on the spatial analysis

To evaluate the uncertainty caused by the 195 undiagnosed positive examinees, we conducted a simulation-based sensitivity analysis for each method as follows:

Step 1: From those who accepted to undergo confirmatory testing, we computed the rates of diagnosed thyroid cancer cases among the positive examinees at primary examination by sex and 5-year age groups for the entire study region. We assumed that the rate was determined by sex and age, with no regional differences in the rates between positive examinees of the primary examination who accepted or declined to undergo confirmatory testing.

Step 2: Applying the binomial probability of being diagnosed cases with the rate computed above for each group of undiagnosed examinees by sex and age groups, we simulated the number of diagnosed cases among undiagnosed examinees for each municipality using random number generation. By adding them to the observed number of cases, we created a hypothetical dataset if every positive examinee underwent confirmatory testing.

Step 3: Using the simulated number of cases, we conducted the same analyses for cluster detection and Poisson regression.

Step 4: We repeated steps 2 and 3 for T times to obtain the distribution of statistical testing results, such as p-values.

We implemented this sensitivity analysis for flexibly shaped scan statistics, MEET, and Poisson regression using R scripts with T = 100.