Design and participants
We used data from the US National Health and Nutrition Examination Survey (NHANES). The NHANES is an ongoing survey conducted by the Centers for Disease Control (CDC) and Prevention that uses a representative sample of non-institutionalized civilians in the US; selected by a complex, multistage, stratified, clustered probability design. Information on participants was collected by interviews and in personal physical examinations. The interview includes background information such as socio-demographic, dietary and health-related questions. The examination component consists of medical and physiological measurements, as well as laboratory tests. The National Center for Health Statistics Ethics Review Board approved all NHANES protocols, and all survey participants completed a consent form. The detailed protocol on NHANES methodology and data collection is available on https://www.cdc.gov/nchs/nhanes/index.htm. For this study, we initially included adults aged ≥20 years enrolled in the NHANES-III between 1988 and 1994, with data on blood lead and urinary cadmium concentrations (n = 16,040). Exposure and covariate data from NHANES-III were then linked to the National Death Index mortality data.
Mortality data
A full description of mortality linkage methods is available from the National Center for Health Statistics (NCHS) (18). Briefly, the de-identified and anonymized data of the NHANES III participants were linked to National Death Index mortality files based on 12 identifiers for each participant (eg, Social Security number, sex, and date of birth) with a probabilistic matching algorithm to determine mortality status. The NCHS public-use linked mortality file provides mortality follow-up data from the date of NHANES III survey participation until December 31, 2015 (1988–2015). Participants with no matched death record at this date were assumed alive during the entire follow-up period. In a validation study using mortality-linked data from the first NHANES study (NHANES-I; 1971–75), 96% of deceased participants and 99% of those still alive were classified correctly (19). The underlying cause of death were recorded in the public-use linked mortality files using the following ICD-10 codes: Cardiovascular diseases including heart diseases (I00-I09, I11, I13, I20-I51) and cerebrovascular diseases (I60-I69) and malignant neoplasms (C00-C97). From the 16,040 participants with baseline data, 1,729 had missing data on mortality and other covariates. The final study sample included 14,311.
Measurements of blood lead and urinary cadmium
Blood and urine samples were collected during the medical examination. The laboratory methods for the processing of these samples are described in detail elsewhere (20). Briefly, the blood and urine specimens were frozen (−30°C and −20°C; respectively), stored, and shipped for analysis to the Division of Laboratory Sciences, National Center for Environmental Health at the Centers for Disease Control and Prevention in Atlanta, Georgia (USA). Lead (µg/dL) concentration was measured in whole blood using inductively coupled plasma mass spectrometry. Urinary cadmium was measured by graphite furnace atomic absorption with Zeeman background correction using the CDC and Prevention modification (20) of the method proposed by Pruszkowska et al. (21). Specimens were analyzed in duplicate, and the average of the two measurements was reported. The detection limit were 1.0 μg/dL (0.048 μmol/L) and 0.03 μg/L for the blood lead and the urinary cadmium; respectively. For study participants who had concentrations of blood lead below the level of detection (n= 1217; 8.5%), values were imputed using LOD/√2 [0.7 μg/dL]. Urinary creatinine measured using the Jaffe reaction with a Beckman Synchron AS/ASTRA Clinical Analyzer (Beckman Instruments, Inc., Brea, CA), was used to account for urine dilution.
Covariates
Baseline covariates were collected when individuals participated in a household interview and demographic information—including sex (male/female), age (continuous; years), ethnicity (Mexican-American, other Hispanic, not Hispanic), poverty to income ratio (categorized in tertiles), the number of years of education attended and completed (continuous; years), area of residence (metro and non-metro counties), and smoking status (current, former and never) was obtained. Information on body-mass index ([BMI] continuous; kg/m2), physical activity (None, One to 14 times, 15 or more times; per month) and overall dietary quality indexes (continuous) was obtained during the medical examination. Dietary intake was collected using a 24-hour dietary recall. We derived the diet quality indexes as measured by the Healthy Eating Index 2015 (HEI-2015) (22) and the adapted dietary inflammatory index (23), from the daily intakes of foods/beverages, energy, and nutrients of the 24-hour dietary recall. A complete description of the development of these scores is described in Supplementary Methods.
Statistical analysis
Complete data on exposures, covariates, and mortality were available for 14,311 participants. We log-transformed (base 2) Blood lead and urinary cadmium concentrations to reduce the influence of outliers and descriptive and bivariate analyses are reported as geometric means and geometric standard errors (SE) by population characteristics. We used the parametric g-computation to estimate risk ratios (RR) and risk differences (RD) of all-cause and specific causes of mortality under hypothetical interventions. In 1986, Robins introduced the g-methods, a class of causal inference techniques that allows building outcome prediction model based on observed quantities, and then predicts potential outcomes under potential hypothetical intervention (24, 25). In recent years, there have been substantial advances in the application of this method, and have been used for instance to evaluate hypothetical interventions on sources of lead exposure on BLLs (26). The parametric g-formula is a generalization of the standardization method and allows to flexibly simulate and estimate survival curves to visualize time-specific effect estimates of any form of hypothetical intervention. A more detailed discussion of this method is presented elsewhere (27, 28). Briefly, we first fitted a pooled logistic regression conditional on covariates and follow-up time, after arranging the data into a person-time structure. The discrete-times hazards of all-cause and specific causes of mortality for each 2-fold increase in the baseline metals concentrations (log2-transformed to reduce skewness) were then estimated. This estimated risk was used to predict the mortality under two hypothetical interventions specifying thresholds for metals concentrations. We then compared the estimated risk of mortality under the intervention 1: had all participants were assigned a high concentration (e.g. the 95th percentiles values of the metals distributions) with the estimated risk of mortality under the intervention 2: had all participants were assigned a low concentration (e.g. the 5th percentiles values of the metals distributions). This approach assumed a linear association between metals concentration and death; to check the assumption, we then used multivariable restricted cubic splines with three knots placed at the 10th, 50th, and 90th percentiles of each metal concentrations distribution to provide a graphical presentation (29). Splines allowed us to test whether there was any departure from linearity.
Finally, we also considered interventions comparing quartile groups of lead and cadmium concentrations. We categorized metal concentrations into quartiles and estimated the discrete-times hazards of all-cause and specific causes of mortality for each quartile with the first quartile group, the lowest metal concentration, as the reference category. We then compared the estimated risk of mortality under the intervention 1: had all participants belonged to the bottom quartile group with the estimated risk of mortality under the intervention 2: had all participants belonged to the low quartile group.
There is no know safety levels for the blood and urinary concentrations of these metals in adults. We therefore chose the interventions listed above based on previous epidemiologic analyses, as discussed elsewhere (14). Models were adjusted for age, sex, ethnicity, poverty index, educational level, area of residence, smoking status, BMI, physical activity, diet quality evaluated by the healthy eating index and metals concentrations (mutual adjustment). The selection of potential confounders was done a priori. We also included product terms between the metals concentrations and time in all models to account for the time-varying risk. All analyses were weighted by the provided sample weights to account for the unequal probabilities of inclusion and response rates. We estimated 95% confidence intervals (CIs) of the RR and the RD using non-parametric bootstrap (M=200) and used the 2.5th and 97.5th percentiles as the lower and upper confidence interval limits, respectively samples.
We additionally investigated age (< 50 years and ≥ 50 years)- and sex-specific estimates in stratified analyses as previous studies reported potential effect modification of the associations between metals concentrations and mortality by sex and age (12, 14). We also evaluated effect modification by age and sex using Cochran Q tests.
Sensitivity Analyses
We ran sensitivity analyses by performing unweighted models that not account for the NHANES survey unequal probabilities of inclusion and response rates because the weighted method is inefficient analysis due to the large variability in assigned weights (30). The unweighted analysis yields correct estimates when models are adjusted for the auxiliary variables used to define the weights (i.e., age, sex, and ethnicity) (30).
Statistical analyses were performed using R version 4.0.4 and the Statistical Analysis System software, version 9.4 (SAS Institute, Cary, NC, USA).