Principal Component Pursuit for Pattern Identification in Environmental Mixtures

We present PCP as a robust method for dimensionality reduction and pattern identification [19]. Given an exposure data matrix $X_{n,p}$, where $n$ is the number of observations and $p$ is the number of pollutants, PCP seeks to express $X$ as a superposition of two matrices: a low-rank matrix $L_{n,p}$ where $r=\mathrm{rank}(L)\ll\min(n,p)$, and a sparse matrix $S_{n,p}$ where most entries are zero. Because $L$ is of rank $r\ll p$, its rank can correspond to underlying patterns in exposure, such as specific sources or certain behaviors. $L$ is still defined in terms of the original variables, i.e., patterns are not directly estimated. PCP may be paired with various matrix factorization techniques (e.g., SVD, PCA, factor analysis, or NMF) to extract chemical loadings and individual scores. The rank of $L$ or the number or location of nonzero entries in $S$ do not need to be a priori defined.

We incorporated two PCP extensions that suit features of environmental mixtures data. First, Zhang et al. [20] recently proposed $\sqrt{PCP}$ with a noise-independent universal choice of regularization parameters. Previous formulations of PCP required knowledge of the true noise level to determine the appropriate parameters [23, 24, 22]. This is problematic in environmental mixtures where we cannot know or accurately estimate the underlying noise level, and it would leave the researcher with the subjective task of tuning parameters on a per-dataset basis. Zhang et al. provide a more practical approach to pattern recognition in environmental mixtures.

As first proposed, PCP minimizes a weighted combination of the nuclear norm of $L$, $\|L\|_{*}$, and the $\ell^{1}$ norm of $S$, $\|S\|_{\ell^{1}}$ [25, 19, 23]. Notably, this formulation has the desirable quality of convexity, meaning that every local optimum is a global optimum. This, together with the particular structure of the $\ell^{1}$ and nuclear norms guarantees that the resulting optimization problem can be solved efficiently [26, 27, 28]. In practice, however, the nuclear norm assumes a stronger low-rank structure (i.e., slowly decaying singular values) than what is the case in many real-life environmental mixtures (e.g., POPs or air pollution). To address unsatisfactory performance with the nuclear norm, we replaced it with a rank-$r$ projection. While the nuclear norm is convex, the rank-$r$ projection is not. However, closely related non-convex formulations are accompanied by theoretical guarantees of equivalent performance with the convex implementation [21, 22]. Combining a non-convex rank projection with $\sqrt{PCP}$, we solve the following optimization problem:

where $X$ denotes the original data matrix. The two parameters, $\lambda$ and $\mu$, are not tuned by the researcher; instead, they are each set using single universal values, $\lambda=1/\sqrt{n}$ from Candès et al. [19], and $\mu=\sqrt{p/2}$ from Zhang et al. [20] which have been shown theoretically to yield near-optimal estimation performance. The indicator function $\mathbbm{1}_{\mathrm{rank}(L)\leq r}$ constrains $L$ to be of rank $\leq r$; the $\ell^{1}$ norm $\|S\|_{1}$ is the sum of the absolute values of the entries of $S$ and encourages $S$ to be sparse; the final term is the error between the predicted and the observed values, which favors a solution that is close to the original data.

We simulated 100 exposure matrices for all combinations of two mixture sizes, three noise structures, and three detection proportions (1800 total). We generated datasets of 500 observations each, $\left\{\mathbf{x}_{i}\right\}_{i=1}^{500},$ where $\mathbf{x}_{i}=\left(x_{i,1},\ldots,x_{i,p}\right)^{\mathrm{T}}$ presents an exposure profile with $p$ mixture components. We specified $r=4$ underlying patterns and investigated two mixture sizes ($p=16$ and $p=48$). We first simulated chemical loadings ($r\times p$) to represent realistic environmental patterns where some chemicals were distinct to a single pattern and some chemicals appeared in multiple patterns. Each pattern included $\dfrac{p}{8}$ chemicals that loaded distinctly and $\dfrac{p}{4}$ chemicals that overlapped with a second pattern. Distinct chemicals were given a loading of 1 on the single pattern on which they loaded and a loading of 0 for the remaining patterns. One third of the chemicals appeared in only one pattern; two thirds of the chemicals appeared in two patterns. This design corresponds to multiple environmental sources giving rise to the chemicals in the mixture. Overlapping chemicals were drawn from a Dirichlet distribution so that their loadings would sum to 1 over all patterns. Of the four loadings across the four patterns for each chemical, two were drawn from $\operatorname{Dir}(\alpha_{1}=1,\alpha_{2}=1)$ and two were set to zero. This introduced variability into the overlapping chemical loadings (Figure 1(a)).

We next generated individual scores ($n\times r$). We drew scores independently from $\log\mathcal{N}(\mu=1,\sigma=1)$. We created the simulated data from matrix products of individual scores by chemical loadings with added noise, replacing negative values with zero. We generated noise in one of three ways, (1) low Gaussian noise ($\mathcal{N}(0,1)$), (2) high Gaussian noise ($\mathcal{N}(0,5)$), (3) or low Gaussian noise with high sparse events. Figure 1(b) shows an example simulated correlation matrix. Finally, we designated a quantile (25^th, 50^th, or 75^th) and set all values below the threshold as $<$LOD.

For pattern recognition in an environmental mixture with varying detection limits across chemicals, we chose a mixture of dioxins, furans, and polychlorinated biphenyls (PCBs) measured in U.S. adults from the 2001–2002 NHANES cycle. NHANES inclusion criteria have been reported previously [29]. For the chosen cycle, 11,039 participants were interviewed. One third of participants aged 12 years and older were eligible for environmental chemical analysis. We removed individuals below 18 years of age or without any POP measurements, resulting in a final study sample of 1,000. Eighteen PCBs, seven dioxins, and nine furans were measured. Exposure assessment of POPs in NHANES has been described previously [30, 31]. Of the POPs measured, 21 detected in at least 50% of all samples were included in our analyses. All POP values were lipid-adjusted by the U.S. Centers for Disease Control and Prevention (CDC) [32].

We determined the appropriate rank for PCP-LOD and the number of components to retain from PCA in the same manner for all experiments. For PCA, we a priori defined our component retention criterion as the first $k$ components that explained $\geq 80\%$ of the variance in the data, as seen previously in environmental mixtures applications (e.g., Gibson et al. [4]). While it is possible to perform cross-validation on PCA [33, 34], it is not a common practice in applied environmental health research. For PCP-LOD we used the default parameters for $\lambda$ and $\mu$ and cross-validated to select the rank of the $L$ matrix. We set an initial grid of rank values from 1 to 10 for all scenarios. We performed this cross-validation approach on a single representative dataset for each combination of simulated mixture sizes ($p$ = 16 and 48), proportions $<$ LOD (25%, 50%, and 75%), and noise structures (low, high, and sparse) and for the POP mixture.

To cross-validate PCP-LOD on a single dataset, $X$, we repeated the following steps 100 times for each rank $r\in[1,10]$. (1) We randomly corrupted 20% of the mixture $X$ as missing (i.e., set the value to NA) to serve as a held-out test set, denoted $X_{\Omega}$, yielding the corrupted matrix $\tilde{X}$. (2) We ran PCP-LOD on $\tilde{X}$ to obtain $\hat{L}$ and $\hat{S}$. (3) We recorded the relative recovery error of $\hat{L}_{\Omega}+\hat{S}_{\Omega}$ compared with the observed data $X_{\Omega}$ in the held-out set, calculated via the Frobenius norm, $||X_{\Omega}-\hat{L}_{\Omega}-\hat{S}_{\Omega}||_{F}\,/\,||X_{\Omega}||_{F}$. Finally, for each rank, we aggregated the average relative recovery error across 100 runs and chose the optimal rank, $\hat{r}$, as that with the lowest mean relative recovery error on the held-out set. We subsequently ran PCP-LOD on the full dataset $X$ with the selected rank $\hat{r}$.

We ran PCP-LOD and PCA on all simulated datasets. We compared PCP-LOD and PCA to assess their relative performance when faced with large proportions of non-detectable observations. For PCP-LOD we estimated the rank of $\hat{L}$, the sparsity of $\hat{S}$, and their relative change to assess stability of the solution across increasing proportions of data $<$ LOD. Because the sparse matrix may contain non-zero values so close to zero as to be considered zero, we set a threshold above which to regard values as legitimate extreme exposures. We evaluated sparse events two standard deviations of the model residuals ($\hat{S}+\hat{\varepsilon}$), per chemical, from zero, i.e., $2\times\sqrt{\mathrm{Var}(X_{p}^{\mathrm{obs}}-\hat{L}_{p}^{\mathrm{obs}})}$, where $\mathrm{obs}$ indicates values above the LOD in the simulated data.

For both PCP-LOD and PCA, we calculated relative predictive error as the ratio of the error to the truth in terms of their Frobenius norm: $\left\lVert\mathrm{Truth}-\mathrm{Predicted}\right\rVert_{F}/\left\lVert\mathrm{Truth}\right\rVert_{F}$. For PCP-LOD we interpreted $\hat{L}$ as the predicted values, and for PCA we constructed predicted values as the product of the score matrix (i.e., the coordinates of the rotated data on the principal components) by the rotation matrix (i.e., right eigenvectors), truncated at the chosen rank. We defined the ‘truth’ as simulated values before noise or sparse events were added. Finally, we assessed the stability of the identified patterns using the relative prediction error of the singular value decomposition (SVD).

Prior to the application to the NHANES POP mixture, we examined distributional plots and descriptive statistics for all variables. We scaled all exposure concentrations by their standard deviations to make variances comparable across chemicals. The solution, thus, cannot be influenced by high-variance pollutants. We used PCP-LOD to separate unique events from underlying patterns. Following PCP, we extracted individual scores and pattern loadings from $\hat{L}$ using SVD. We compared scores, loadings, and overall relative error with those obtained from PCA. We present unique events and interpret observed patterns. All analyses were conducted using R version 4.0.4 [35].

We ran PCP-LOD and PCA on all simulated datasets. PCP-LOD had lower relative prediction error across the majority of mixture sizes ($p$ = 16 and 48), proportions $<$ LOD (25%, 50%, and 75%), and noise structures (low, high, and sparse). PCP-LOD outperformed PCA on all simulations with low noise, simulations with high noise with up to 50% $<$ LOD, and simulations with low noise and added sparse events with up to 50% $<$ LOD (Figure 2). Figures 2 and 3 present simulations where $p$ = 16; corresponding figures where $p$ = 48 are included in Supplemental Figures S1 and S2.

PCP-LOD was more affected by the proportion of data $<$ LOD, which can be seen in the larger step size between box plots in Figure 2. The decline in PCP-LOD predictive accuracy as the proportion of values $<$ LOD increased appears because of poorer performance on values $<$ LOD in high noise scenarios (Figure 3). Relative prediction error for values $>$ LOD was approximately constant for PCP-LOD and PCA. Supplemental Tables S1, S2, and S3 contain the median and inter-quartile range (IQR) of relative error for predicted values overall and stratified by LOD.

Next, we assessed the stability of the identified patterns using the SVD of the simulated data before noise or sparse events were added and compared this with the SVD of the $\hat{L}$ matrix and of PCA results. Figure 4 depicts the relative prediction error comparing the left eigenvectors (comparable to scaled individual scores) of the PCP-LOD and PCA solutions with those of the simulated ‘truth.’ PCP-LOD’s median relative prediction error is generally lower than PCA’s for the larger mixture size and higher than PCA’s for the smaller mixture size. However, these patterns appear quite stable over increasing proportions of data $<$ LOD for both methods. PCP-LOD solutions achieved lower relative prediction error on chemical loadings (i.e., right eigenvectors) across all simulations (Supplemental Figure S3).

Across PCP-LOD solutions, between 2% and 10% of $\hat{S}$ entries were non-sparse. We found decreasing sparsity as the proportion $<$ LOD increased, with 3% (IQR: 2%, 4%), 6% (IQR: 4%, 7%), and 7% (IQR: 3%, 8%) unique events, on average, found in simulations with 25%, 50%, and 75% $<$ LOD, respectively. For simulations that included sparse events in the noise structure, PCP-LOD correctly included 69% (IQR: 67%, 71%), 70% (IQR: 68%, 72%), and 65% (IQR: 62%, 67%) of sparse values in the $\hat{S}$ matrix, on average for simulations with 25%, 50%, and 75% $<$ LOD, respectively.

Thirty-four POPs were measured in the NHANES 2001–2002 cycle. Detection frequency is presented in Figure 5. Fourteen PCBs, four furans, and three dioxins were detected in $>$ 50% of samples. Exposure levels of POPs were all positively correlated (Figure 6(a)).

We applied PCP-LOD to identify underlying patterns of POP exposure and extreme exposure events that were not explained by these patterns without making a priori assumptions concerning the number of patterns or sparse events. PCP-LOD returned a low-rank matrix of rank three, which corresponds with three patterns of POP exposure in the $\hat{L}$ matrix. Figure 6(b) depicts $\hat{L}$’s correlation matrix along side the correlation matrix of the raw data (Figure 6(a)). By removing sparse events and residual noise, PCP-LOD increased the correlations between POPs. To characterize underlying patterns, we extracted principal components from the low-rank matrix using SVD.

The three components distinguished by PCP-LOD included one component of overall POP exposure, a component that separated dioxins and furans from PCBs, and a third component that separated higher molecular weight PCBs from lower molecular weight PCBs (Figure 7). The first component explained 79.4% of the variance in the low-rank matrix, the second explained 14.6%, and the third explained 6.0%.

PCP-LOD partitioned the variation that was unexplained by the low-rank structure into a sparse matrix of large outlying values and the remaining residuals. The $\hat{S}$ matrix contained mostly zero values, with 5.7% of entries being non-sparse. Sparse observations were generally weakly correlated, with the absolute value of r $<$ 0.15 for 70% of Spearman correlations between sparse chemical exposure events. Table 1 summarizes the number of individuals with uniquely high $\left(>2\times\sqrt{\mathrm{Var}(X_{p}^{\mathrm{obs}}-L_{p}^{\mathrm{obs}})}\right)$ or low $\left(<-2\times\sqrt{\mathrm{Var}(X_{p}^{\mathrm{obs}}-L_{p}^{\mathrm{obs}})}\right)$ exposure events. Figure 8 describes participant-specific sparse events. Most participants had no extreme exposures (44%) or only extremely low exposures (18%). Twenty-two percent had one high unique event on a single chemical, and 16% had between two and six high exposures across 21 chemicals left unexplained by the identified patterns.

PCA conducted on the POP mixture chose three components that explained $\geq$ 80% of the variance and returned loadings and scores much the same as those from $\hat{L}$ (results not shown). Using the three chosen components, PCA’s relative prediction error on values $>$ LOD was 0.30, similar to PCP-LOD’s relative error of 0.32 when comparing only $\hat{L}$ with the original data. However, when including $\hat{S}$ in the solution ($\hat{L}+\hat{S}$), PCP-LOD’s relative error on values $>$ LOD was 0.07. This is more comparable to the PCA solution when including all 21 components, 0.06, which does not accomplish any dimension reduction.