Applying Latent Variable Models to Estimate Cumulative Exposure Burden to Chemical Mixtures and Identify Latent Exposure Subgroups: A Critical Review and Future Directions.

There are many open questions in chemical mixtures research; below, we detail a few questions and posit how latent variable methods could be used to address them. Liu et al. [ 29 ] introduced the innovative use of item response theory (IRT) to create summary exposure burden scores for chemical mixtures, specifically focusing on PFAS mixtures. They demonstrated that in IRT, the concentrations of the individual PFAS biomarkers, as well as exposure patterns to the PFAS mixture, collectively influence the latent exposure burden score through data-driven, nonlinear functions that are unrelated to specific health outcomes. The authors also highlighted that IRT facilitates comparisons of exposure burden across studies even if they do not measure the exact same set of chemicals. They showed that IRT provides a straightforward means to include exposure biomarkers with low detection frequencies, and can help mitigate exposure measurement errors. To promote harmonization and future meta-analyses, they provided a PFAS burden score calculator for 2017–2018 based on recent US population reference ranges. This calculator enables researchers to calculate PFAS burden scores on a common scale, even if the studies did not measure exactly the same set of PFAS biomarkers, and found that findings would remain robust to associations with health outcomes. When quantifying cumulative exposure burden in diverse populations, in which different people may be exposed to different sets of chemicals due to heterogeneity in exposure sources and patterns, it is important to consider whether a single exposure burden scoring algorithm is appropriate for the entire population. Customized scoring algorithms, which could account for systematic differences in exposure sources within the population, may be needed. If these characteristics are known a priori (for example, exposure burden burden is hypothesized to differ across sex), then multiple-group IRT may be used. However, if researchers hypothesize that complex, unknown combinations of socio-demographic, dietary and behavioral characteristics underlie different exposure patterns to the chemical mixture such that the groups are not known in advance, then mixture IRT (MixIRT) approaches may be used. MixIRT combines IRT and LCA, by identifying latent subpopulations characterized by different scoring algorithms, using anchor item(s) to set a common burden scale across latent subpopulations, and estimating MixIRT exposure burden based on a participant’s weighted likelihood of belonging to each latent subpopulation. Liu et al. [ 10 ] demonstrated the use of MixIRT to develop a customized exposure burden scoring algorithm, to ensure that PFAS burden scores can be equitably compared across population subgroups. They found that using MixIRT, they could detect that Asian Americans had significantly higher PFAS burden compared with non-Hispanic Whites, but this disparity was sometimes masked when using summed PFAS concentrations as the exposure metric. The researchers suggest that it may be important to account for sources of exposure variation when estimating burden scores, to ensure findings are fair and informative for all people. Another potential use of latent variable modeling is to determine how many burden scores are needed. Accordingly, latent variable models for multidimensional [ 109 ], bifactor [ 110 ] and hierarchical [ 111 ] data structures are available. In the literature on latent variables, the problem is often phrased as determining the dimensionality of the data. Statistical methods to determine the number of dimensions to retain is a major topic in exploratory factor analysis. Auerswald and Moshagen [ 112 ] provide an overview and comparison of these methods using simulated data. Rather than discussing the mechanics of these procedures, we will instead focus on some practical considerations when determining the number of latent variables to use for a given application. First, it is recommended to use multiple methods to find converging evidence for the appropriate number of factors to use [ 112 ]. Second, it may be informative to fit, compare, and interpret multiple models with different numbers of factors [ 112 , 113 ]. For instance, one may wish to use a unidimensional latent variable models to estimate an overall score, even if multiple subfactors are identified, especially if the subfactors are correlated with each other [ 114 , 115 ]. Finally, it is important to also consider whether the models make theoretical sense [ 116 ], as burden scores based on an implausible models may not validly represent true latent burdens. For example, Chen et al. [ 22 ] developed phthalate burden subscores using multi-dimensional IRT. They compared model fits of three theoretically driven IRT models and found that a correlated factors model representing low molecular weight phthalate burden, high molecular weight phthalate burden and DEHP phthalate burden provided the optimal fit. They showed that associations of phthalate burden and insulin resistance were consistent with findings using molar sums, but the IRT approach also enabled researchers to account for potential missingness in some phthalate metabolites and allowed researchers to account for measurement error when verifying consistency of associations. There are limitations, however, to IRT approaches. Specifically, they do not allow researchers to predict the effect of chemicals not included in the original sample. However, because IRT can be used to develop an overall burden score, it may be particularly useful for studying an entire class of chemicals, such as PFAS, in which some PFAS have been phased out of production and replaced by other PFAS that exert similar health effects (known as regrettable substitution). While the concentrations of individual chemical biomarkers may vary over time due to regrettable substitution, the overall exposure burden to the chemical class may remain high. Further, IRT may be more appropriate for chemical classes that are known toxins, as an underlying assumption is that the probability of having a higher concentration of a biomarker increases monotonically with exposure burden to the chemical class. In the case of some mixtures, such as metal mixtures in which some metals are toxicants and others are nutrients, or an individual metal could be both a nutrient and toxicant depending on exposure concentration, simple IRT methods may not be appropriate as the scale and interpretation of the latent trait would be unclear. However, it may be possible to develop latent variable models that represent such scenarios. There is limited but emerging research applying such latent structures to polytomous data in educational contexts [ 117 , 118 ]. Data harmonization is needed in consortia and other research, to achieve a sample size to detect smaller effect sizes or to make inferences across a heterogeneous population. IRT provides tools to place test scores onto a common scale, using the common items as anchors, so that they can be compared across students even if they did not take the same test version [ 119 ]. In educational testing, in order to protect test integrity, different groups of students are often administered different versions of a test that include a set of common items. In the environmental health setting, cohorts may measure slightly different sets of chemical mixtures, with a common set of overlapping chemical biomarkers. Using IRT, we can standardize measures of mixture burden across studies, so that the burden score can be compared, even if the studies did not measure exactly the same set of chemicals. Existing approaches for studying mixtures are limited to analyzing only the common set of chemicals assayed in all studies. IRT could make full use of the exposure data, even if some-what different sets of chemicals are measured in different studies. In latent variable modeling, especially in the context of IRT models, data harmonization is also known as linking [ 120 ], and its major goal is to place latent variables on the same scale for all groups. In other words, appropriate harmonization or linking is necessary for latent variable scores that originally came from different datasets to be comparable to each other. For example, if different sets of chemicals are measured for two different groups, the latent burden scores are not necessarily directly comparable. Generally speaking, for scores to be comparable, there must be some common element across the datasets [ 121 ]. One strategy is to assume that both groups have the same latent burden distribution, known as the equivalent groups design [ 120 ]. However, this assumption is often unrealistic, and may only be an appropriate strategy if the datasets represent random samples from the same population. In such a case, the latent variable models can be estimated separately, assuming the same (e.g., standard normal) latent trait distribution for both groups and the resulting latent variables will automatically be on the same scale. Another strategy is to have one or more chemicals in common across the two datasets. Here, the chemicals in common are known as the anchor (which functions similarly to the anchor discussed above for DIF detection). In this design, it may be possible to estimate a common latent variable model for all groups [ 122 ] or estimate the models separately and transform one model to the scale of the other [ 122 ]. A third scenario may arise if a common chemical is not available across the datasets one wishes to combine. Normally, this lack of commonality would prevent the fair comparison of latent variable scores. However, it may be possible to harmonize these datasets if another dataset can be identified that includes some chemical measurements that belong to each of the original datasets. This strategy is known as an external anchor and similar statistical techniques may be used as are available for the anchor items strategy [ 123 ]. An important caveat to all these methods is that they assume that the appropriate latent variable model is the same for all groups. In other words, these procedures assume that no chemicals exhibit differential item functioning, and the presence of differential item functioning can negatively affect the accuracy of the linking procedure [ 124 ]. Therefore, great care should be taken to select an appropriate set of anchors that are unlikely to function differently in different groups.

The discipline of psychometrics, or the measurement of latent (unobservable) constructs underpinning psychological and educational assessment, has primarily been used in psychology to study constructs such as ability [ 11 ], depression [ 12 ], and anxiety [ 13 ]. Although there are a small number of studies using latent variable models in chemical mixtures research [ 14 – 17 ], we believe that there is great potential to answer new and increasingly sophisticated research questions using these models. The derivation of latent constructs can be conceptualized and operationalized in a variety of ways, and we briefly describe a few methods commonly deployed in psychometrics. Table 1 provides an overview of latent variable methods and the research questions they address with respect to our focus on quantifying exposure burden and assessing latent subpopulations. Factor analysis (FA) [ 18 ] follows from the insight that a set of observed variables can be summarized more efficiently as a smaller set of latent variables. Since Spearman’s discovery, FA has become a widely used approach in the social and behavioral sciences. FA can be either exploratory (EFA) or confirmatory (CFA), depending on the constraints specified. EFA typically has no or very few a priori specifications on the number of factors or the pattern of relationships between the factors and observed variables. Further, it can be used as an exploratory or descriptive analysis to determine whether all observed variables in the model accurately represent the latent construct, or whether some observed variables do not meaningfully contribute to the explanatory variance in the overall model. CFA is often conducted after an EFA and typically includes more researcher-imposed constraints, such as specifying the number of factors to fit and the relationship of observed variables to each other, and model specification is often guided by theory. Structural equation modeling (SEM) is a class of techniques that enables researchers to simultaneously define latent variables, and also study the relationship of those latent variables with other latent or observed variables. SEM consists of a measurement portion, which creates latent variables from the observed variables, and a structural portion, which assesses the relationships of these latent variables and other covariates [ 14 , 19 ], and it can be extended by including group variables so that parameters for each group are different [ 14 ]. Both SEM and FA have close relations with item response theory (IRT), which models the probability of a categorical response as a function of one or more underlying unobserved traits, often denoted as theta ( θ ). The prevailing idea is that individual item scores result from an interaction between qualities of the item (e.g., item difficulty) and qualities of the person (e.g., the ability of the person), permitting the estimation of individual scores on a common latent continuum and precise standard errors of those estimated scores. Since its conception, IRT has been used in a variety of fields including education testing, personality assessment, patient-reported outcomes measurement, clinical assessment, and epigenetic data [ 10 , 20 – 35 ]. Applications of IRT methodology include incorporating multidimensionality (the idea that an item assesses not one but multiple constructs) and differential item functioning (the idea that subgroups may respond asymmetrically to the same item). An important distinction is that not all unsupervised learning methods are latent variable models. In general, unsupervised learning methods are concerned with predicting the value or distribution of an unobserved variable [ 36 ]. Latent variable models are a subset of unsupervised learning methods in which the unobserved variable to be predicted is (in theory) causally related to the observed variables used to infer it [ 37 ]. Even though latent variable models are built on the assumptions of causality, simply fitting a latent variable model is not evidence of causality (see [ 37 , 38 ]). Unsupervised learning methods that are not latent variable models include principal components analysis [ 39 ] and k-means clustering [ 40 ] whereas latent variable models include factor analysis, item response theory, structural equation modeling, and more [ 41 ]. Latent variable models often serve similar goals (e.g., dimension reduction) and use similar computational strategies (e.g., eigenvalue decomposition) as the broader class of unsupervised learning methods. However, latent variable models are usually more scientifically plausible because they more faithfully represent the causal mechanisms theorized to underlie the data and inherently account for the probabilistic relationship between the observed and latent variables. Therefore, we advocate the use of latent variable models where possible (e.g., using FA instead of PCA). Latent class analysis (LCA) [ 42 – 45 ] is a family of methods used to identify previously undefined subgroups, called latent classes , based on a set of observed variables. The observed variables in LCA are categorical and often binary. LCA has enjoyed increasing popularity in recent years and across areas such as phenotype identification [ 46 ]. Given a well-developed theoretical and statistical justification for how many distinct groups can best represent the data [ 44 ], LCA can uncover hidden groups or subpopulations underlying the data. Similarly, latent profile analysis (LPA) attempts to estimate an underlying latent variable of group membership, using continuous variables rather than categorical variables [ 47 ]. Given that a central question in environmental health measurement is appropriate quantification of exposure burden during critical windows of susceptibility, another approach may be latent change point analysis, which can identify unusual sudden changes in an individual’s score pattern compared to their baseline score pattern [ 48 , 49 ]. We have briefly discussed commonly used latent variable modeling methods in the field of behavioral and social sciences, though we note they are only a small subset of analytic techniques that the discipline of psychometrics can offer to novel analysis of environmental health outcomes. We now discuss the prior work using latent variable modeling in environmental epidemiology.

As aforementioned, LCA and LPA are latent variable modeling methods that have been used extensively in education, psychology, and sociological research [ 42 ]. While they have also been used in clinical and epidemiologic research settings to characterize underlying, mutually exclusive subgroups (e.g., endotypes, phenotypes) based on symptomatology or disease biomarkers [ 62 – 65 ], their use with environmental chemical exposure measures is limited. Importantly, while several studies have used LCA or LPA to identify subgroups of alcohol, tobacco, or illicit drug use [ 66 – 71 ], we do not discuss them here as the purpose of those studies was often to understand use behaviors and their risk factors as opposed to underlying chemical exposures resulting from such behaviors. To the best of our knowledge, only three prior studies have leveraged a standard LCA approach in the context of characterizing environmental chemical exposure mixtures [ 72 – 74 ]; none have used LPA. One study focused on exposures to a single chemical class (polycyclic aromatic hydrocarbons [PAHs]) among all non-pregnant 2013–2014 National Health and Nutrition Examination Survey (NHANES) participants, and relations of the resulting latent PAH classes with body mass index (BMI) [ 73 ]. The same authors conducted an LCA of 47 urinary biomarkers of exposure to six chemical classes (metals, pesticides, environmental phenols, phthalates, polycyclic aromatic hydrocarbons, and volatile organic compounds), among child and adolescent (ages 6–19 years) participants in the 2011–2012 cycle of NHANES, and evaluated associations with lymphocyte and neutrophil cell counts [ 74 ]. Most recently, LCA was applied to 15 urinary biomarkers of phthalates and environmental phenols, among pregnant women from LIFECODES, a prospective birth cohort, and the resulting classes were validated against those identified among pregnant women from NHANES cycles 2005–2006 and 2007–2008 [ 72 ]. Additionally, authors estimated associations of the latent classes with urinary biomarkers of oxidative stress. While these studies introduce and demonstrate the utility of LCA for expressing complex, high-order interactions among environmental exposure mixtures as readily interpretable, mutually exclusive exposure mixture subgroups, there are important limitations for future studies to overcome. All three of these studies dichotomized exposure biomarker concentrations based on percentiles (e.g., median or 75th percentile) as opposed to using the continuous concentrations directly with an LPA approach [ 75 ]. As such, resulting latent classes may not optimally characterize the underlying exposure mixture subgroups given the loss of information through dichotomization as well as sensitivity of findings to chosen cut-offs. Additionally, none of these studies accounted for the uncertainty (i.e., measurement error) in latent class assignment when conducting latent class regression, which is important for obtaining unbiased estimates of associations with latent classes, especially when the LCA/LPA does not have a high entropy (i.e., high certainty in class assignment) [ 76 ]. There are also alternative parameterizations of the standard LCA/LPA model, or extensions thereof, which have yet to be utilized for the study of environmental exposure mixtures. For example, whereas all aforementioned studies used an unsupervised approach (i.e., latent classes estimated independent of outcomes), it is possible to estimate latent classes and latent class membership jointly with a distal outcome, approximating a supervised approach. Further, few studies have yet extended LCA or LPA to identify subgroups of exposure mixtures based on longitudinal exposure measures, either as trajectory subgroups (i.e., repeated-measures LCA/LPA) [ 77 ] or by latent transition analysis (LTA) [ 42 ]. To date, only one study has utilized a repeated-measures LCA/LPA in the context of longitudinal chemical mixtures. Kuiper et al. used this approach to characterize longitudinal, serum PFAS concentrations measured at multiple timepoints from delivery through midchildhood, among children enrolled in the Health Outcomes and Measures of the Environment (HOME) Study cohort, finding two latent subgroups: “higher PFAS” and “lower PFAS” (higher or lower relative to the other latent group [ 78 ]). Other authors have extended the broad concept of underlying, latent categorical variables to develop methodologies that more readily accommodate high-dimensional chemical exposure biomarkers. Zhang et al. [ 79 ] developed a latent variable modeling approach that jointly estimated latent profiles of 62 semicontinuous polychlorinated biphenyl (PCB) congeners measured in serum and associations with endometriosis risk. Specifically, their approach utilized a Bayesian framework and random effects to allow for the identification of the most toxic PCB with respect to endometriosis [ 79 ]. Similarly, a recent study leveraged this same Bayesian latent variable modeling framework to jointly assess latent profiles of serum PCB biomarkers with couple’s infertility risk, estimating separate latent profiles for each member of the dyad [ 80 ]. Finally, four studies have applied a latent unknown clustering integrating multiomics data (LUCID) approach [ 81 ] to identify clusters of individuals based on exposure biomarker and high-dimensional metabolomics data and evaluated associations of cluster membership with health outcomes [ 82 – 85 ]. While not yet routinely utilized by researchers, these methods demonstrate the powerful generalizability and flexibility of latent variable models to identify underlying exposure subgroups.

Advancements in latent variable modeling and psychometrics can be used to address contemporary questions in mixtures research, particularly around the understudied topics of quantifying cumulative exposure burden to mixtures and identifying hidden subpopulations with distinct exposure patterns to mixtures. We believe there is great potential to answer new and increasingly sophisticated research questions using these models.

The above sections demonstrate some ways in which latent variable models have been used for environmental mixtures research. However, we believe that the full strength of these models has not yet been realized. For example, latent variable models can be used to identify ways in which the appropriate measurement model for a particular biomarker differs across groups. In literature on latent variable models, these techniques are known as differential item functioning (DIF) [ 86 ] or measurement invariance analysis [ 87 ]. The presence of DIF suggests that for individuals with the same level of the latent trait (e.g., exposure burden), participants from different subgroups have a different probability of having a certain concentration level of that biomarker. There may be differences in exposure sources that affect mixture burden for different socio-demographic subgroups (e.g., sex, race/ethnicity, socioeconomic status). For example, females are more likely to use personal care products that contain certain phthalates; individuals identifying as Black are less likely to use sunscreen and thus have lower oxybenzone exposure [ 88 ]. When analyzing DIF (or measurement invariance), it is important to distinguish overall group differences (i.e., population differences in latent burden scores) from differences that are idiosyncratic to a particular biomarker/item. The former is known as differential impact , and if it exists, it is appropriate to use the same latent variable model for all groups. The latter is known as DIF, and if it exists, different latent variable models should be used for that biomarker for different groups. While environmental epidemiology studies have recognized the potential for measurement invariance/DIF, they have not always tested for it in a rigorous manner. For example, in a study on endocrine-disrupting compounds (EDCs) and thyroid hormones, the authors hypothesized that the loadings of EDCs on the latent exposure variable and the hormonal loads, as well as the binding capacity of thyroid hormones could be different by sex. Thus, they provided separate SEMs for males and females [ 17 ], but did not perform statistical tests to confirm meaningful DIF and whether separate models were needed. Many methods exist for identifying DIF, some of which are not specific to latent variable models [ 89 ], but many of which are based on a latent variable approach. Within the latent variable framework, one popular method is the likelihood ratio test (LRT) [ 90 ]. When using the LRT to test a specific biomarker for DIF, two multiple-group latent variable models are fit. In one model, the parameters for that biomarker are fixed to be the same for all groups, and in the other model, the parameters for that biomarker are freely estimated. The relative fit of the two models is then compared by taking the difference in the model deviances. Under the null hypothesis that the two models fit equally well, this difference in model deviances follows a chi-square distribution with degrees of freedom equal to the number of parameters associated with that biomarker. Therefore, a significant p -value indicates statistically significant DIF that the biomarker does not function in the same way in both groups. In the LRT, it is necessary to fix the parameters for at least one biomarker, known as the anchor , to be the same for all groups. That is, at least one biomarker must be identified that does not exhibit DIF. Various strategies exist for identifying an appropriate anchor [ 91 , 92 ]. In addition, the more the biomarkers that are used as anchors, the greater the statistical power to detect DIF, but at the risk of an inflated Type I error rate [ 93 ]. If possible, it is recommended to use 10–20% of the biomarkers as anchors [ 94 ]. The previously described methods for detecting DIF items are most often used when attempting to identify individual items that function differently across groups. Another major approach focuses on establishing no differential functioning for the entire set of items simultaneously. This approach involves sequentially testing various measurement invariance models that make increasingly strict assumptions [ 95 – 97 ]. Conventions vary, but typically begin by fitting a configural invariance model in which separate latent variable models are fit to each group. In the configural invariance model, the same latent variable model (i.e., the same pattern of which biomarkers represent which latent variables) is fit separately to each group, but there are no assumptions that the model parameters are the same across groups (however, some restrictions are necessary for model identification, such as setting the latent mean and standard deviation to 0 and 1 for each group). Model-data fit can then be evaluated using a suite of overall fit indices [ 98 , 99 ], and good fit for the configural invariance model indicates that the biomarkers measure the same latent variable(s) for all groups. If good model-data fit is not found for the configural invariance model, scores from the different groups represent qualitatively different latent variables. In this situation, different sets or patterns of biomarkers would need to be identified if the researcher hoped to measure the same latent variable in all groups. In contrast, if the configural invariance model fits well, the metric invariance model is fit next. The metric invariance model constrains the discrimination (or factor loading) parameters to have the same values across groups. The absolute fit of the metric can be assessed using overall fit indices and by comparing the relative fit of the metric invariance model to the configural invariance model [ 98 , 100 ]. If the metric invariance model does not fit well, the cause of the invariance (i.e., the particular biomarkers that require different discriminations for different groups) can be investigated, and a partially invariant model may be identified [ 101 ]. Notably, the metric invariance model, like the configural invariance model, requires model identification constraints such as fixing the latent mean and standard deviation to 0 and 1 for all groups. This identification requirement prevents scores from different groups from being compared directly, even if a well-fitting (partial) metric invariance model is identified. Therefore, once this well-fitting (partial) metric invariance model is found, a third scalar invariance model may be fit that modifies the previous model by also constraining threshold or intercept parameters to be the same values across groups. The same statistical techniques mentioned above can be used to compare the scalar invariance model to the metric invariance model and to identify a partially invariant model if needed. An important feature of the scalar invariance model is that the identification restrictions (e.g., setting the latent mean and standard deviation to 0 and 1) only is required for one group, and the latent mean and standard deviation can be freely estimated for the remaining groups. Therefore, once a well-fitting (partial) scalar invariance model is identified, latent burden scores across groups can be directly compared to each other. A final level of invariance is strict invariance and requires that the residual variance for each biomarker is the same across groups. Strict invariance does not affect the comparability of the latent burden scores and rarely holds in practice. Therefore, many researchers omit this final step of the measurement invariance analysis [ 96 ]. Although the sequential measurement invariance models and the LRT for DIF employ similar models, they approach the analysis in distinct ways. Namely, the LRT focuses on identifying differential functioning for individual biomarkers whereas the sequential approach proceeds by testing different types of invariance for all biomarkers simultaneously. In the sequential approach, tools are available to identify partial invariance, and those biomarkers that exhibit non invariance are those that exhibit DIF. In practice, the choice of method may be informed by the purpose of the research and researcher expectations about the likely prevalence of DIF. When DIF is detected for a particular biomarker, there are several strategies for how to proceed with analyses. For example, it may be appropriate to use different measurement models with different parameters for each group [ 102 ]. However, note that scores from different groups will only be comparable to each other in this case if the model parameters come from a multiple-group model with at least one item that has the same item parameters (i.e., the same measurement model) for all groups. Other strategies include using more complicated multidimensional models [ 103 ] or omitting the differentially functioning biomarkers. Liu and Rogers [ 102 ] provide some general advice on when to use each strategy, but generally recommend against dropping differentially functioning items. Another limitation of existing work is that measurements are often taken to be known with certainty. However, latent variable models typically indicate that there is non-negligible measurement error associated with the scores. If not explicitly accounted for, these measurement errors can bias correlations with other variables and reduce statistical power [ 104 ]. Although there are a small number of studies [ 15 , 16 , 105 , 106 ] that have used structural equation models to properly account for measurement error, there is opportunity for these methods to become more widespread. A potential solution to quantify and thus account for measurement error is plausible value imputation, which uses the observed score distributions to produce multiple datasets. These datasets can then be aggregated to produce statistics at the sample level, accounting for the variance and facilitating calculation of standard errors. Such techniques [ 107 ], while popular in wide-scale education assessment and patient-reported outcomes research [ 108 ], have yet to be explored in the context of exposure burden research.

Structural equation modeling (SEM) has been used to assess the association of biochemical exposures and health outcomes [ 16 , 17 , 19 , 50 – 58 ]. In the measurement portion of SEMs, latent variables, either unidimensional or multidimensional, have been created from one or more classes of chemical mixtures to reflect an overall exposure [ 16 , 17 , 19 , 51 – 53 , 55 – 57 ]. Unidimensional models have been used to reflect an overall exposure: Grandjean et al., used SEM to calculate a latent variable that reflects an overall exposure to perfluorinated compounds (PFCs) [ 19 ], Mogensen et al. calculated an overall measure of exposure to perfluorinated alkylate substances (PFASs) [ 16 ] from PFOS, PFOA, and PFHxS concentrations; Grandjean et al. used PFOS, PFOS, PFHxS, PFNA, and PFDA to estimate an overall effect of PFAS mixtures [ 51 ]. Wang et al. used 10 urinary OH-PAHs metabolites to measure an overall polycyclic aromatic hydrocarbons (PAHs) exposure [ 57 ]. Application of SEMs has also been used to generate latent variables for air pollution chemicals. Jaafari et al. and Tu et al. applied SEMs to calculate latent variables for air pollution using CO, NO 2 , O 3 , SO 2 , and particulate matter (PM 1 , PM 2.5 , and PM 10 ) in their analyses [ 53 , 56 ]. There are also examples using SEMs to generate latent exposure variables from more than one class of chemicals, e.g., Przybyla et al. created a latent variable reflecting exposure to endocrine-disrupting chemicals (EDCs) from seven phthalates, three phenols, and perchlorate [ 17 ]. Multidimensional models have also been used to create latent variables for multiple classes of chemical exposures, or for time-varying exposure data. For example, Shook-Sa et al. developed a structural equation model to assess the effects of volatile organic compounds (VOCs) and tobacco use on adolescents’ respiratory function, which included a latent tobacco exposure variable and two latent variables for VOC exposures (xylene/ethylbenzene and toluene) [ 54 ]. Buncher et al. suggested that SEMs could be more informative in longitudinal environmental research when data are measured at different time points [ 59 ]. In a study assessing the impacts of prenatal and postnatal exposure to PCBs on children’s antibody response to immunizations, Heilmann et al. first fit SEMs for prenatal and postnatal exposures separately and then developed a model including both prenatal and postnatal latent exposures, in which prenatal exposure affected postnatal exposure and both latent variables were associated with outcomes [ 52 ]. Mogensen et al. included two latent PFAS exposure variables from chemical concentrations measured at age 5 and age 7 in their SEM investigating the impacts of PFAS chemicals on children’s antibody responses, and they also estimated a joint latent PFAS exposure variable taking the measurements at two time points into account [ 16 ]. Bayesian inference for SEMs has also been used: when studying the impacts of traffic related air pollution exposure on health outcomes, Baja et al. developed SEMs in a Bayesian framework to account for air pollution measured at multiple time points [ 60 , 61 ]. SEMs can be used to assess and adjust for imprecision of chemical exposure when exposure data are measured using different methods. For example, when estimating the latent mercury exposure, Grandjean et al. used four measures of children’s mercury exposure across different tissue (measured in cord blood, umbilical cord tissue, maternal full-length hair, and maternal short-segment hair). They accounted for measurement imprecision for maternal hair and cord measurements, by allowing the measurement errors to be correlated from these two types of tissue [ 50 ]. They also used structural equation modeling approach to adjust for imprecision in their research determining the effects of PFCs as well as PFASs exposure on antibody responses to vaccinations in children [ 16 , 19 , 51 ]. Welch et al. also applied SEMs in a longitudinal study in which exposure variables had been measured differently. In their study, they used SEMs to assess the association of concentrations of diphtheria and tetanus antibodies and latent variables of exposure to metals (arsenic, lead, and manganese) measured in two places (in drinking water and in blood) at three time points (pregnancy, toddlerhood, and childhood) [ 58 ].

Environmental health and risk assessment research increasingly emphasizes understanding the health impacts of environmental mixtures, as compared with previous interest in single environmental exposures. This shift acknowledges that individuals are consistently exposed to multiple environmental agents simultaneously [ 1 ]. Subsequently, there has been a proliferation of supervised statistical methods that address exposure mixtures [ 2 ], which model the impacts of environmental mixtures in association with a pre-specified health outcome. Statistical methods have been developed to address the potential independent, joint, interactive, nonlinear or overall effects of a mixture on a pre-specified health outcome [ 3 – 7 ], with the method selection dependent on the research question of interest [ 8 , 9 ]. However, because these methods simultaneously model the relationship between the chemical mixtures and the health outcome, findings may vary across related health outcomes even when they relate to the same underlying physiological system, leading to inconsistent results on how to mitigate the impacts of mixtures on health. An attractive alternative, therefore, is the use of unsupervised methods to characterize and estimate burden to exposure mixtures. Using an unsupervised approach would allow for transportability of exposure mixture characterization across studies with different health outcomes of interest. Latent variable models, a suite of methods foundational in psychometric research, provide several unsupervised approaches that have infrequently been used in other fields, including exposure science and environmental epidemiology. Latent variable models may be used to characterize and/or quantify exposure to complex chemical mixtures, particularly in the context of cumulative exposure burden or patterns of exposure mixtures. Cumulative or total exposure burden could be represented by a single summary metric (or set of summary metrics) to quantify exposure to all components of a mixture, analogous to using a polygenic risk score to summarize genetic risks, or an IQ score to summarize cognitive ability. This would also be useful for biomonitoring and risk assessment, when the interest lies in identifying individuals with high total exposure burdens. Researchers may also be interested in comparing burden across race/ethnicity groups and socioeconomic strata, and geographical regions. For example, Liu et al. recently applied mixture item response theory (MixIRT) to characterize cumulative exposure burden to per- and polyfluoroalkyl substances (PFAS) across self-identified racial and ethnic groups in the National Health and Nutrition and Examination Survey (NHANES). They found higher cumulative PFAS exposure burden among Asian Americans, a novel finding that was masked when more traditional approaches for summarizing cumulative exposure were used [ 10 ]. A summary metric of exposure burden would allow researchers to compare the relative “contribution” of exposure burden to chemical mixtures across a range of health outcomes, while maintaining a consistent metric of exposure for the chemical mixture. A distinct but related approach may involve assessing patterns of exposure to chemical mixtures (cross-sectionally or longitudinally), in order to determine how exposures co-occur, and to identify hidden subpopulations with differing exposure burden and patterns of exposure. Such an approach may be particularly useful in succinctly characterizing complex, higher-order interactions among components of the chemical mixture, and especially useful in policy-oriented research given the readily interpretable characterization of exposure patterns (e.g., “high exposure group,” “low exposure group,” etc.). Given the broad utility of latent variable models and infrequent application in understanding complex chemical mixtures, the focus of this review is on how latent variable methods can be used to address two specific, understudied questions in environmental epidemiology: (1) quantifying cumulative exposure burden to mixtures and (2) identifying hidden subpopulations with distinct exposure patterns to mixtures. We will provide an overview of latent variable methods and review prior applications, along with summarizing limitations of existing work and exploring future directions for how advanced latent variable methods could be used to address contemporary mixtures questions, such as creating an exposure burden scale, and addressing harmonization of exposure burden across studies. Our goal is to demonstrate that transdisciplinary applications of novel latent variable methods from the psychometrics literature could be used to address pressing questions in environmental epidemiology.