Libraries

library(tidyverse)
library(foreign)
library(stats)
library(viridis)
library(survey)
library(plotly)
library(kableExtra)
library(gridExtra)
library(sjPlot)

Data Import

First, we read in NHANES cycles 2003 - 2004, 2005 - 2006, 2007 - 2008, 2009 - 2010, 2011 - 2012, 2013 - 2014, and 2015-2016. As of June 2020, phthalate data for 2017 - 2018 was not yet available. The NHANES datasets used are Demographics (DEMO) and Phthalates and Plasticizers Metabolites - Urine (PHTHTE). For the 2015-2016 cycle we also need the Albumin & Creatinine (ALB_CR_I) file since creatinine data were removed from the phthalates files during this cycle (more on creatinine below).

# reading in 2003 - 2004 data
demo_03_04 = read.xport("./data/phthalates/2003_2004_DEMO_C.XPT")
phthalate_03_04 = read.xport("./data/phthalates/2003_2004_PHTHTE_C.XPT.txt")

# reading in 2005 - 2006 data
demo_05_06 = read.xport("./data/phthalates/2005_2006_DEMO_D.XPT")
phthalate_05_06 = read.xport("./data/phthalates/2005_2006_PHTHTE_D.XPT.txt")

# reading in 2007 - 2008 data
demo_07_08 = read.xport("./data/phthalates/2007_2008_DEMO_E.XPT")
phthalate_07_08 = read.xport("./data/phthalates/2007_2008_PHTHTE_E.XPT.txt")

# reading in 2009 - 2010 data
demo_09_10 = read.xport("./data/phthalates/2009_2010_DEMO_F.XPT")
phthalate_09_10 = read.xport("./data/phthalates/2009_2010_PHTHTE_F.XPT.txt")

# reading in 2011 - 2012 data
demo_11_12 = read.xport("./data/phthalates/2011_2012_DEMO_G.XPT.txt")
phthalate_11_12 = read.xport("./data/phthalates/2011_2012_PHTHTE_G.XPT.txt")

# reading in 2013 - 2014 data
demo_13_14 = read.xport("./data/phthalates/2013_2014_DEMO_H.XPT.txt")
phthalate_13_14 = read.xport("./data/phthalates/2013_2014_PHTHTE_H.XPT.txt")

# reading in 2015 - 2016 data (note change in creatinine source file for this cycle)
demo_15_16 = read.xport("./data/phthalates/2015_2016_DEMO_I.XPT.txt")
phthalate_15_16 = read.xport("./data/phthalates/2015_2016_PHTHTE_I.XPT.txt")
creat_15_16 = read.xport("./data/phthalates/2015_2016_ALB_CR_I.XPT.txt")

Data Tidy

Next we’ll bind the data files for each cycle using left-joins, merging on the unique identifier SEQN.

data_03_04 = 
  left_join(demo_03_04, phthalate_03_04, by = "SEQN") %>% 
  mutate(cycle = "2003-2004")

data_05_06 = 
  left_join(demo_05_06, phthalate_05_06, by = "SEQN") %>% 
  mutate(cycle = "2005-2006")

data_07_08 = 
  left_join(demo_07_08, phthalate_07_08, by = "SEQN") %>% 
  mutate(cycle = "2007-2008")

data_09_10 = 
  left_join(demo_09_10, phthalate_09_10, by = "SEQN") %>% 
  mutate(cycle = "2009-2010")

data_11_12 = 
  left_join(demo_11_12, phthalate_11_12, by = "SEQN") %>% 
  mutate(cycle = "2011-2012")

data_13_14 = 
  left_join(demo_13_14, phthalate_13_14, by = "SEQN") %>% 
  mutate(cycle = "2013-2014")

data_15_16 = 
  left_join(demo_15_16, phthalate_15_16, by = "SEQN") %>% 
  left_join(creat_15_16, by = "SEQN") %>% 
  mutate(cycle = "2015-2016") 

all_data = 
  bind_rows(data_03_04, data_05_06, data_07_08, data_09_10, data_11_12, data_13_14, data_15_16) 

# select variables of interest, drop every observation without biomarker data, filter out women aged 20-44, and consolidate household and family income variables
all_data = 
  all_data %>% 
  select(SEQN, cycle, RIAGENDR, RIDAGEYR, RIDRETH1, DMDEDUC2, DMDMARTL, INDHHINC, INDFMINC, INDHHIN2, INDFMIN2, WTMEC2YR, SDMVPSU, SDMVSTRA, URXUCR, URXECP, URXMBP, URXMC1, URXMEP, URXMHH, URXMHP, URXMIB, URXMOH, URXMZP, URXMNP, URXCOP, URXMNP, URXCNP) %>%
  drop_na(URXMEP) %>% 
  rename(gender = RIAGENDR, age = RIDAGEYR, ethnicity = RIDRETH1, education = DMDEDUC2, marital_status = DMDMARTL, creatinine = URXUCR, mecpp = URXECP, mnbp = URXMBP, mcpp = URXMC1, mep = URXMEP, mehhp = URXMHH, mehp = URXMHP, mibp = URXMIB, meohp = URXMOH, mbzp = URXMZP, mcnp = URXCNP, mcop = URXCOP, minp = URXMNP) %>% 
  filter(gender == 2, age %in% (20:44)) %>% 
  mutate(
    household_income = if_else(cycle %in% c("2003-2004", "2005-2006"), INDHHINC, INDHHIN2), 
    family_income = if_else(cycle %in% c("2003-2004", "2005-2006"), INDFMINC, INDFMIN2)
    ) %>% 
  select(-c(INDHHINC, INDHHIN2, INDFMINC, INDFMIN2))

all_data_c = 
  all_data %>% 
  mutate(mecpp_c = 100*mecpp/creatinine, 
         mnbp_c = 100*mnbp/creatinine,
         mcpp_c = 100*mcpp/creatinine,
         mep_c = 100*mep/creatinine, 
         mehhp_c = 100*mehhp/creatinine,
         mehp_c = 100*mehp/creatinine, 
         mibp_c = 100*mibp/creatinine, 
         meohp_c = 100*meohp/creatinine, 
         mbzp_c = 100*mbzp/creatinine, 
         mcnp_c = 100*mcnp/creatinine,
         mcop_c = 100*mcop/creatinine, 
         minp_c = 100*minp/creatinine)

A brief historical interlude

In 1976, Congress passed the Toxic Substances Control Act, which sought to measure and regulate industrial usage of chemicals deemed to pose “unreasonable risk” to humans and the environment.

The EPA was tasked with administration of the act and since passage, an inventory of about 84,000 chemicals has been established. In short, the inventory is problematic. First, it’s only updated every four years. Second, it gives exemptions to smaller manufacturers and an ever-growing list of exempt chemicals. Finally (and most importantly), quantities are often reported in extremely wide ranges (“1,000,000 - <20,000,000 lbs” was an entry I saw in the 2016 data…) (Zota et al, 2014, Sciences et al).

Basically, I once assumed that some governing body has a solid grasp on the quantities and types of chemicals used in the US. I no longer assume this.

In another facet of this regulatory picture, we have the Consumer Product Safety Improvement Act of 2008. The act focuses on toxic substances in childrens’ toys and banned presence of certain phthalates (check out this disturbing pile of dolls seized by US Customs in 2013 due to excess phthalate levels).

One might think this would provide some clues on why MEP declined, but the timing doesn’t work out- the act went into effect in 2009 and our trend started in 2004/2005. Additionally, DEP is not included in the scope, it’s hard to ascribe the act as the root cause. It is possible, however, that the industry responded to pressure from advocacy groups and consumers, or in anticipation of further bans, and undertook formulation changes outside of (and prior to) passage of the act.

Back to the numbers

Up until now, we’ve explored the temporal trends of common phthalates and did a bit of unsuccessful (but fun?) digging through the exciting world of toxic chemical legislative history. Now we will backtrack and summarize the average levels, as posed by Question 1.

We proceed to summarizing the raw and creatinine-adjusted values from the 2015-2016 NHANES cycle. Mean and median values for both measures are reported in Table 1 below. Both are shown because there is a high degree of right-skew in the data. To illustrate this, here is a quick histogram of unadjusted MEP levels.

all_data %>% 
  ggplot(aes(x = mep)) +
  geom_histogram(binwidth = 500) +
  labs(x = "Unadjusted MEP (ng/mL)", 
       title = "Figure 2: MEP Histogram")

This right-skew is predictable - most people (thankfully) have very low levels of urinary MEP. As such, the median is a better measure of central tendency than the mean. We will look at both since some of the literature I reference below uses mean values and keeping the means will allow for comparison.

# calculating means and medians for phthalate values, not adjusted for creatinine
means_raw = 
all_data_c %>% 
  filter(cycle == "2015-2016") %>% 
  select(-c(gender:creatinine), -SEQN, -c(mecpp_c:minp_c)) %>% 
  pivot_longer(c(mecpp:mcnp), names_to = "chemical", values_to = "value") %>% 
  drop_na() %>% 
  mutate(chemical = recode(chemical, 
                           "mecpp"= "Mono-2-ethyl-5-carboxypentyl phthalate (MECPP)", 
                           "mnbp"= "Mono-n-butyl phthalate (MnBP)",
                           "mcpp"= "Mono-(3-carboxypropyl) phthalate (MCPP)",
                           "mep"= "Mono-ethyl phthalate (MEP)",
                           "mehhp"= "Mono-(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP)",
                           "mehp"= "Mono-(2-ethyl)-hexyl phthalate (MEHP)",
                           "mibp"= "Mono-isobutyl phthalate (MiBP)",
                           "meohp"= "Mono-(2-ethyl-5-oxohexyl) phthalate (MEOHP)", 
                           "mbzp"= "Mono-benzyl phthalate (MBzP)", 
                           "mcnp"= "Mono(carboxyisononyl) phthalate (MCNP)", 
                           "mcop"= "Mono(carboxyisoctyl) phthalate (MCOP)",
                           "minp"= "Mono-isononyl phthalate (MiNP)")) %>%
  group_by(chemical) %>% 
  summarize(mean = mean(value), median = median(value)) %>% 
  mutate(id = row_number())

# Calculating means and medians for adjusted values
means_c_15_16 = 
  means_c %>% 
  filter(cycle == "2015-2016") %>% 
  mutate(id = row_number()) %>% 
  select(id, chemical_c, mean_c, median_c)

# joining adjusted and un-adjusted into one table
left_join(means_c_15_16, means_raw, by = "id") %>% 
  select(c("cycle", "chemical", "mean", "median", "mean_c", "median_c")) %>% 
  rename("Chemical" = "chemical", "Adjusted Mean" = "mean_c", "Adjusted Median" = "median_c", "Mean" = "mean", "Median" = "median") %>% 
  knitr::kable(digits = 1, caption = "Table 1: Phthalate Concentration in Women Ages 20-44, per NHANES 2015-2016 Cycle.") %>% 
  kable_styling(c("striped", "bordered")) %>%
  add_header_above(c(" " = 2, "Unadjusted (ng/mL)" = 2, "Creatinine-Adjusted (ug/g creatinine)" = 2))

So what do these values tell us? Are they safe? And what does “safe” even mean in this context?

Figuring out the answers gets us into a gray area. Some (myself included) would argue that no amount of plasticizers should be present in the human body. However, given how frequently we interact with plastics, this is probably not a realistic goal.

Alternatively, we can draw certain conclusions about threshold values above which detrimental health effects were observed in past studies, and compare those values to what we saw in Table 1. But this, too, is tricky. For one thing, there’s no way to ethically study the effects of phthalates in humans using the gold standard of study design - the randomized controlled trial (RCT). You can’t just pick a group of pregnant women, split them up, and force one group to eat large amounts of phthalates.

This leaves us with two options: animal studies and observational human studies. The former has the widely-acknowledged limitation that humans and lab animals (usually rats) are not, in fact, biologically equivalent. The latter has analytic limitations, the major one being unmeasured confounding.

Confounding is a huge consideration in health/medical research and many epidemiologists spend their entire careers studying its effects and how it distorts measured relationships between exposures and outcomes (you can find a nice primer here). For our purposes it’s sufficient to say that observational studies of phthalates might tell us, for instance, that women with high levels of urinary MEP gave birth to smaller babies, but we cannot say that phthalates were the cause unless we control for all other potential causes of small babies (randomization accomplishes this, hence the usage of RCTs in humans). Both options have pretty big limitations, but if we combine the knowledge from both, we can perhaps glean conclusions to help us make sense of Table 1.

Here is what I found:

• Several studies have established adverse effects in rats, as summarized by Lyche et al (Lyche et al, 2009). The outcomes were primarily related to reproductive effects and included sperm count reduction in males (NTP-CERHR Expert Panel Update on the Reproductive and Developmental Toxicity of Di(2-ethylhexyl) phthalate and delayed puberty in females Grande et al, 2006. Looking at the doses, however, these effects were observed in the range of 100 - 2000 mg chemical/kg animal per day. If we take a number in this range, say 375 mg/kg animal/day, we get a dose of roughly 150 mg per rat, per day (the average adult rat weighs about 400 grams. If we extrapolate this to humans, using an average weight of 68 kilograms (150 lbs) for women, we get that there are about 170 rats per human (not literally of course), meaning that the human equivalent of 150 mg is 25,500 mg, or 25 grams of phthalate (almost 1 oz) per day. This is pretty huge and I highly doubt that anyone is ingesting this much on a daily basis.

• Human studies on reproductive outcomes are limited, but Polanska et al (Polanska et al, 2016) found an association between MEP and pregnancy duration in a cohort of Polish women. The other outcomes studied were weight. length, head and chest circumference of the baby and no significant associations were found between these outcomes and any other phthalate. In this study, the median creatinine-adjusted MEP was 22.7 g/g creatinine, which is on par with the value calculated here, (33.2 g/g creatinine). So barring unmeasured confounding in this study, we would conclude that MEP levels in US women are still too high and pose risk for adverse reproductive outcomes.

• In a cohort of US women who gave birth to males, Swan et al (Swan et al, 2005) found an association between the boys’ anogenital distance and mothers’ exposure to MEP, MnBP, MBzP, and MiBP, with MnBP having the biggest effect. In this study, the boys with the shortest anogenital distance (i.e. the most reproductive impairment) had median concentrations of MEP, MnBP, MBzP, and MiBP of 225 ng/mL, 24.5 ng/mL, 16.1 ng/mL, and 4.8 ng/mL, respectively (the authors did not adjust for creatinine). Aside from MiBP, these values are much higher than the median values computed above.