Study setting

The study involves major asbestos-user countries in which historical asbestos consumption and mesothelioma mortality data are available. Males are a target population of this study, as the main gender involved in occupational exposure and related health effects. Being based on country specific consumption data, the study addressed asbestos exposure as a whole, without distinguishing it between occupational and environmental sources, using a longitudinal approach by calculating exposure history with a pre-defined lag period with a further limitation on the working age population. Next, we modelled the relationship between past asbestos consumption and mesothelioma mortality to estimate country-specific exposure-response, lag-response and age-response functions.

Asbestos consumption data

Data on asbestos consumption were derived from worldwide asbestos supply and consumption estimations provided by the United States Geological Surveys (USGS) technical reports [6], which include historical country-specific data of national production, import, and export based on information extracted from a variety of sources. Asbestos apparent consumption has been defined as domestic production plus the difference between imports and exports.

As a first step, we selected the major user countries of raw asbestos according to international trade figures, considering all countries with more than two million of raw asbestos consumption in the period 1933–2012. For descriptive trend analysis, we have divided this period into four 20-year sub-periods (1933–1952; 1953–1972; 1973–1992; 1993–2012), analysing the temporal distribution of asbestos consumption, comparatively by country.

Mesothelioma mortality data

Malignant mesothelioma mortality has been retrieved from the WHO mortality database [9]. For each selected country, we have considered the number of males who died from mesothelioma (ICD-10: C.45) in a 20-year period (1998–2017). Russia, China, and India reported no mesothelioma death in the WHO mortality database archives, while data from Brazil were only partially available (not national data). Thailand presented 31 male deaths for mesothelioma in the entire period. For all other selected countries, we calculated the percentage variation in the four, five-year periods. As for the modelling of the association between consumption and mesothelioma mortality, we used yearly data from 1920 to 2018. For modelling the exposure, lag, and age-response functions, mortality data available in five-year age classes were interpolated obtaining occurrences for each yearly age.

Male population data

Male population data during the period of 1950–2017 were retrieved from Global World Statistics 2022, which is based on data provided by World Bank and United Nation Organization [10]. Data are provided in five age classes (0–19; 20–39; 40–59; 60–79; 80+) and were interpolated at an yearly resolution in the interval 0–100 years of age. The total population, needed to calculate per capita asbestos consumption, was interpolated up to the year 1920 using historical data.

Asbestos exposure

To estimate the lifetime exposure to asbestos of the population, we assigned the asbestos consumption per capita by year during the period of 1920–2018 to the whole population regardless of age. We calculated the related exposure history with a pre-defined lag period, setting lag 0 equal to the exposure that occurred during the last year. Therefore, assuming a maximum lag of 60 years, we generated a matrix of exposure histories of the population for each age and year that was used to define temporal dependencies using the model presented in the next section. We assumed that asbestos consumption before 1920 was equal to that estimated for this year. Finally, we assumed that occupational exposure to asbestos was limited to the working age (15–65 years old). The matrix of historical asbestos exposure was then merged by age and year with mesothelioma mortality to fit the statistical model. Modelling was restricted to the countries with reliable mortality data (12 of the 18 original countries). China, Russia and India were discarded due to the lack of mesothelioma mortality data in the WHO database. Thailand and Brazil were discarded due to incompleteness in mortality data, and Mexico for lack of consistency in the asbestos exposure range with respect to other countries. The model was applied to male subjects of age 25–89, considered as those likely involved with mesothelioma mortality due to past asbestos exposure.

Statistical model

We extended the modelling approach used in a study carried out in Italy [7] to determine the association between past asbestos consumption and mesothelioma male mortality, applying it in a multi-country context using a two stages analysis. We have modelled only mesothelioma deaths in men in non-linear models, according to the less straightness of the association in females. In the first stage, we fitted a quasi-Poisson generalized linear model with age/year-specific mesothelioma deaths to model the association between past asbestos exposure and mesothelioma mortality rates in each country.

The model is defined in the equation below:

Where the outcome Y_t corresponds to yearly counts of mesothelioma deaths assumed to follow a Poisson distribution with over dispersion, the function f(E_t;θ) specifies the association with the asbestos exposure at year t, the function $g\left(ag{e}_t;\beta \right)$M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ g\left( {ag{e_t};\beta } \right) \] \end{document} specifies the association with age of the subjects at year t, and log(pop) is the corresponding log-transformed population at year t used as an offset.

The exposure–response associations have been modelled using distributed lag non-linear models (DLNM), a statistical framework to represent complex exposure–lag–response risk surfaces [11]. The function f(E_t;θ) is defined by a bi-dimensional cross-basis term, using a natural spline function with two degrees of freedom and equally spaced knots for the exposure–response dimension, and a natural spline function with three degrees of freedom to model the lagged-response dimension, accounting for 60 years of lag. The age-response function $g\left(ag{e}_t;\beta \right)$M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ g\left( {ag{e_t};\beta } \right) \] \end{document} was represented by a one-dimensional natural spline function with a knot placed at age 65. In the second stage, the cross-basis coefficients obtained for each country were then used to obtain an overall pooled exposure-response function using mixed-effects meta-analysis [12]. We also estimated the best linear unbiased estimators (BLUPs) to obtain country-specific exposure response functions. The meta-analysis was also carried out for both lag-response and age-response functions defined at specific exposure and lag values, respectively. A sensitivity analysis about the dependence of the meta-analysed exposure, lag, and age-response curves, using different degrees of freedom of the natural splines functions, has been carried out.

Trend analysis of asbestos consumption and mesothelioma mortality

According to the criteria defined above, we selected 18 countries with raw asbestos cumulative consumptions higher than two million tons in the period 1933–2012. Eight countries in Europe (Belgium, France, Germany, Italy, Poland, Russia, Spain, and the United Kingdom), five in Asia (China, India, Japan, South Korea, and Thailand), four in North and South America (Brazil, Canada, Mexico, and the USA), together with Australia. According to USGS technical reports, these countries represented 79% of total asbestos consumption in the world during the study period [6].

Table 1 shows the asbestos consumption in each country during four, twenty-year periods and the percent of variation with respect to the previous one. The huge increase in asbestos consumption in the period 1953–1972 is noteworthy, with respect to the previous 20-year period in all countries. In the subsequent 20 years (1973–1992), some countries registered a reduction in asbestos consumption: Canada and United States of America (–11% and –58, respectively), Belgium, France, and the United Kingdom in Europe (–21%, –5% and –56%, respectively), and Australia (–23%). In 1993–2012, all countries presented a decrease in asbestos consumption from the previous period, with the exception of China, India, and Thailand, which consistently still increased the use of asbestos (161%, 127% and 64%, respectively).

Table 2 shows the number of persons who died from malignant mesothelioma in major asbestos-user countries during four periods of five years between 1998 and 2017. Data about China, India, and Russia are missing, as data were not available in the WHO mortality database archives. Thailand is not present in Table 2 due to the inconsistency of reported cases (31 male deaths due to mesothelioma in the period). Most of the selected countries presented an increase in the number of cases by five-year period, with a decreasing trend in percentage variation (Belgium, Canada, France, Germany, Italy, Japan, Mexico, and South Korea). The USA is the only country for which a decrease in the number of deceased subjects is documented in the last period with respect to the previous one (13,056 deaths for mesothelioma in 2008–2012, and 12,659 in 2013–2017). These findings support the positive effects of asbestos ban (or restrictions policies) on mesothelioma mortality trend. Brazil, which included data only for a region, presented an irregular trend with a first decrease between 1998–2002 and 2003–2007 (558 and 366 deaths, respectively) and a clear increasing trend subsequently (412 cases in 2008–2012 and 508 in the subsequent period).

Figure 1 tabulates the association between asbestos consumption (1933–2012) and mesothelioma deaths (1998–2017), showing the distribution of countries. Russia, China, and India have been tabulated with no mesothelioma cases registered. A clear linear relationship can be observed between total consumption and total deaths for mesothelioma.

Results of the statistical models

Country-specific and pooled meta-analytical exposure, lag, and age-response relationships are shown in Figure 2. High heterogeneity is observed among the countries. The overall pooled exposure-response function (Figure 2b) exhibits an increasing risk for asbestos consumption level. Heterogeneity is also detected among the country-specific lag-response functions (Figure 2c). The overall pooled lag-exposure response (latency) for an exposure of 3.9 tons for 10,000 inhabitants shows an interesting pattern. The lag-specific RR grew almost linearly for 20 years since exposure, followed by a drop which reflects the behaviour of some countries. The age-response relationship (Figure 2a) shows a slow increase until age 50 and a sharp increase thereafter with a large scatter for each countries. Figure SM1 displays the full exposure-lag-response through the bi-dimensional representation of the predicted risk over a grid of exposures and lags, obtained from pooled meta-analysis of country-specific results. Yearly country-specific model estimated mesothelioma deaths in men and cases registered in the WHO archive were compared to obtain a goodness of fit of the country model. Cumulative observed and predicted cases, mean squared error (MSE), and Spearman correlation coefficient of yearly data are presented in Table 3 for the years of data availability. The country-specific models show a good ability to reproduce the total observed malignant mesothelioma cases with a small bias. Low MSE values and a good correlation with observed values are also identified.

The exposure, lag, and age response curves obtained in the main analysis were tested for sensitivity to different values of the parameters used in the natural spline functions included in the cross-basis terms of the model. We tested different degrees of freedom for exposure and lag functions, and different number of knots for the age function. Results are shown in Supplemental Material (Figure SM2). In summary, the number of degree of freedom of the exposure function introduces minimal differences both in the exposure and the lag response curves. Similarly, minimal modifications are estimated when the number of degrees of freedom are increased in the lag function. Higher variations in the age-response curves are instead estimated as far as the number of knots in the natural spline function for age is increased. However, such differences are within the 95% confidence interval of the main analysis results.