Individuals within a family are phenotypically more similar to each other than to a random group of people (Kendler et al., 2015; Polderman et al., 2015). This resemblance has long been recognized and arises from multiple sources, including shared rearing environment, parental influences, cultural values, and genetic inheritance (Abdellaoui et al., 2022; Engzell & Tropf, 2019; Polderman et al., 2015). Quantifying the reasons for similarity in specific traits remains challenging because genetic influences are deeply intertwined with environmental factors. For example, in nuclear families, parent-offspring covariance is caused by direct additive genetic effects of the genotype, half of which is shared between parent and offspring, by environmental influences shared between parent and offspring, and by vertical transmission (VT), in which parental phenotypes causally influence offspring through the rearing environment. Genetic nurture refers to the influence of parental genetic effects—including alleles not transmitted to offspring—on parental traits that in turn affect offspring outcomes through VT (Balbona et al., 2021; Kong et al., 2018). Because parental phenotypes often reflect both genetic and environmental factors, VT also induces a covariance between direct genetic effects and the familial environment (Cloninger, 1980), further increasing parent-offspring covariance. Additional covariance can also arise from assortative mating (AM), which increases the similarity between parents’ and extended family members’ genetic and environmental influences. Disentangling these processes is challenging because they often generate overlapping patterns of phenotypic covariance among family members.

Human genetic researchers have a long history of using structural equation modeling (SEM) to disentangle the sources of familial resemblance (Cloninger, 1980; Keller et al., 2009; Lyu et al., 2025; M. Neale & Cardon, 2013; Wright, 1934). SEM is a statistical framework used to test a hypothesized causal model using observed data. SEM finds the variances and covariances among variables that are implied by a hypothesized data-generating mechanism and compares these to observed data, typically using maximum likelihood methods (Kaplan, 2008). By operationalizing constructs such as genetic and environmental factors as latent variables within an SEM, complex intrafamilial processes can be explicitly modeled. This approach makes it possible to partition variation in traits into genetic and environmental components, to model parental influences such as VT, and to test and account for processes like AM or gene-environment (G-E) covariance (Keller et al., 2009). More broadly, SEM is valuable because it forces explicit specification of causal assumptions, directs attention to effect sizes rather than significance tests, and allows for complex extensions—such as the modeling of recursive relationships or incorporating multivariate data—that would otherwise be mathematically intractable (Balbona et al., 2021; Heath et al., 1985; M. Neale & Cardon, 2013).

Classical twin designs have traditionally been used to partition phenotypic variation into additive and dominant genetic effects, shared environmental influences, and unique environmental influences, but their estimates can be biased when both non-additive genetic and shared environmental influences simultaneously influence the trait (Keller & Coventry, 2005), or when VT, genetic nurture, G-E covariance, or AM are present but not explicitly modeled. Family designs incorporating additional relatives, such as adoption and extended twin family designs, have improved estimation of these effects (Keller et al., 2009; Lyu et al., 2025). Adoption studies can isolate direct genetic effects and VT, but additional biases can be introduced if prenatal effects, selective placement, and/or ongoing contact between biological parents and adoptees are not appropriately accounted for (Horn et al., 1979; Rutter et al., 2001; Shih et al., 2004). Extended twin family designs allow estimation of genetic nurture and AM alongside direct genetic effects, using data from twins’ parents, spouses, and children to account for and mutually estimate these various influences (Heath et al., 1985; Keller et al., 2009; Maes et al., 1997). A hallmark of extended twin family models is the use of non-linear constraints, which describe and constrain recursive relationships between parameters in a way that keeps the overall model internally consistent and identified (e.g., AM and G-E covariance; Keller et al., 2009). This feature illustrates one of the key advantages of using SEM for modeling familial effects, as it allows recursive processes—such as cases where A influences B, which in turn influences A—to be represented explicitly and estimated within a coherent framework. However, extended twin family designs require large samples to estimate latent effects with precision; obtaining such samples is difficult, and the models rely on strong assumptions about the sources of phenotypic covariance; violations of these assumptions can lead to biased estimates. The growing availability of genomic data offers a complementary and potentially less assumption-dependent alternative: by combining parent and offspring phenotypes with measured parental polygenic scores (PGSs)—the sum of trait-associated alleles weighted by their effect sizes from independent genome-wide association studies (GWASs)—it is possible to directly model the relative contributions of genetic sharing and VT (Balbona et al., 2021). This general approach also facilitates genetically informed research in large, publicly available family genomic datasets, which are more widely available than extended twin family data.

A number of models incorporating PGSs into family models have now been developed and tested across a range of traits (McAdams et al., 2023). For example, PGSs were incorporated in twin designs to investigate G-E covariance (Dolan et al., 2021) and to test causal influences between traits (Castro-de-Araujo et al., 2023). Demange et al. (2022) used parental PGSs as measured genetic predictors to estimate the environmental, non-genetic effects of parents’ cognitive and non-cognitive skills on their offspring’s educational outcomes, thereby controlling for genetic transmission. Kong et al. (2018) provided the first empirical evidence for genetic nurture, showing that parental PGSs based on non-transmitted alleles, estimated using haplotype-based PGSs from both parents, predict offspring traits via the environment provided by the parent. Building on this approach, the SEM-PGS model of Balbona et al. (2021) quantifies and accounts for confounding from AM and estimates the total effect of VT, defined as the environmental influence of a parental trait on an offspring trait. Thus, genomic data from large family samples, along with increasingly predictive PGSs, provide information that enables estimation of quantities previously inaccessible or allows them to be estimated in novel ways.

A key limitation of existing PGS-family models is that, with few exceptions—such as the MR-DoC2 model of Castro-de-Araujo et al. (2023), which estimates bidirectional causation—they are univariate, focusing on a single phenotype. The lack of multivariate models creates two challenges. First, it limits the ability to investigate cross-trait effects. In intergenerational transmission, for example, one might ask how parental mental health influences offspring educational attainment, accounting for direct genetic effects, genetic nurture, AM, and G-E covariance of each trait. Univariate approaches cannot address such questions adequately; modeling them as univariate problems, as has sometimes been done (Kong et al., 2018), ignores pleiotropy as well as within- and cross-trait AM and VT, and can therefore produce seriously biased cross-trait estimates, as we show below. Second, ignoring cross-trait effects can also bias within-trait estimates. For instance, when estimating the genetic nurture effect of parental BMI on offspring BMI, failing to account for parental educational attainment—which also affects offspring BMI—may upwardly bias the estimated genetic nurture pathway. Given the well-documented phenotypic, genetic, and environmental correlations between traits (Bulik-Sullivan et al., 2015; Cross-Disorder Group of the Psychiatric Genomics Consortium et al., 2013) and the long tradition of multivariate research in behavioral genetics (Hart et al., 2021; McAdams et al., 2014), it is essential to develop PGS-family models that can estimate and account for cross-trait influences, as univariate models are inherently vulnerable to bias from unmodeled pathways. To address this gap, we introduce the multivariate SEM-PGS model—an extension of our prior univariate SEM-PGS framework (Balbona et al., 2021; Kim et al., 2021)—designed to estimate cross-trait familial influences both within and across generations. Below, we present the structure of the model and introduce novel multivariate path-tracing rules. We then use simulations to demonstrate the importance of accounting for multivariate influences and to evaluate the statistical properties of the model’s estimates under two different approaches to fitting the model.

Multivariate SEM-PGS Model

The multivariate SEM-PGS model jointly estimates and accounts for within- trait and cross-trait effects—including direct genetic effects, VT, genetic nurture, and AM—within a unified structural equation modeling (SEM) framework. For clarity of presentation, we focus on a model involving two traits throughout this paper. While extending the model to include three or more traits is mathematically straightforward, given that expectations are derived using matrices rather than scalar equations, such extensions substantially increase model complexity and the number of parameters to be estimated. This added complexity makes interpretation and communication of results much more challenging. Nevertheless, the qualitative conclusions drawn from the bivariate simulations generalize to higher-dimensional models.

The path diagram of the bivariate model is displayed in Figure 1. The identification of model parameters requires seven vectors of observed components (Table 1): three vectors of phenotypic scores (,,) and four vectors of haplotypic PGSs (,,,), all of length two. To simplify notation, we write [N]T_p/m to indicate any of NT_p , T_p , NT_m , or T_m , where subscripts p and m refer to paternal and maternal haplotypes, respectively, and where NT stands for the haplotypic PGS not transmitted and T for the haplotypic PGS transmitted to offspring. For a bivariate model, there are therefore fourteen observed statistics, forming a 14×14 observed variance-covariance matrix. The definitions of all estimated variance and covariance parameters involved in the multivariate SEM are listed in Table 3, and the definitions of path coefficients are listed in Table 2.

Figure 1

Table 1

Table 2

Table 3

Multivariate Path-Tracing Rules

We derived the expected variance-covariance equations using multivariate pathtracing rules (Vogler, 1985). Multivariate path-tracing in SEM is based on the well-known univariate rules (Cloninger, 1980; Wright, 1934), with some additional rules regarding matrix transposition that are necessary to obtain the correct expectations. The ‘copath’ concept was introduced in univariate SEM to model selection processes such as AM, in which two variables covary without sharing a common antecedent cause (Cloninger, 1980). Formal path-tracing rules for copaths in multivariate SEM have not yet been described and require special considerations. Here, we review the multivariate path-tracing rules and extend them to incorporate copaths. Our focus is on the multivariate case, which builds directly on the univariate rules. Readers unfamiliar with the univariate path-tracing rules are encouraged to review them first, as this background will clarify the logic and interpretation of the multivariate extensions (Balbona et al., 2021).

Multivariate path tracing rules: To aid comprehension of the rules, we first introduce two heuristic definitions: upstream variable and downstream variable. For variables connected by a single-headed arrow, the variable the arrow points to is defined as downstream, and the variable from which the arrow originates is defined as upstream, as this convention is intuitively consistent with how causation flows from one variable to the next. Unlike “exogenous” and “endogenous” terms in SEM, which are qualities of variables themselves, “upstream” and “downstream” are relative terms that refer to how variables relate to one another: a variable may be upstream of one variable and downstream of another. With these definitions in mind, we list below the multivariate path tracing rules for all components in an SEM model:

• Rules for path coefficient matrices: A path coefficient matrix is transposed when traversing from the upstream variable to the downstream variable (with the direction of the arrow) (Wright, 1934). It is untransposed when traversing from the downstream variable to the upstream variable (against the direction of the arrow). If a path coefficient matrix is a full matrix, the upstream variable should be placed along the columns of the matrix, and the downstream variable should be placed along the row of the matrix. For example, f is a full path coefficient matrix that connects the upstream Y_p/m variable to the downstream F_o variable , and thus

  

• Rules for covariance matrices: A covariance matrix is represented by a double-headed arrow in path diagrams.Variables connected by double-headed arrows are not truly upstream or downstream of one another—they occupy the same level in a causal pathway. However, because covariance matrices can sometimes be full rather than symmetric (due to them being a subsection of a larger symmetric variance-covariance matrix), a transposition rule is needed. Thus, as with other matrices, a full covariance matrix is transposed when traversing from the upstream variable to the downstream variable. As long as users apply a consistent convention when arbitrarily designating such variables as upstream or downstream, the path-tracing rules described here will yield the correct variance-covariance expectations. For example, one possible approach, and the one we recommend, is to consistently assign the parental variables as upstream of offspring ones and female spouses upstream of male spouses. As with other path coefficient matrices, upstream variables should be placed along the columns and downstream variables along the rows. For example, if we define the covariance between (upstream) and (downstream) as θ_Tm, then

  

  so when doing the path tracing, we have a direct genetic effect path δk, four paths through the increased genetic covariance from AM 2δg_c + 2ai_c and a path through G-E covariance . Adding up all paths lead us to the correct covariance θ_Tm = δ_k + 2δg_c + 2ai_c + .

• Rules for copath matrices: The copath models multivariate AM and is represented by a straight line, without arrows, connecting two variables (here, Y_p and Y_m). It is a full matrix, and is transposed when traversing from the upstream variable to the downstream variable. As with a covariance matrix, one variable connected by the copath should be designated as the upstream variable and the other as the downstream variable, and this convention should be applied consistently throughout the derivation of all parameter expectations. Here, we adopt the convention that the maternal spouse is upstream of the paternal spouse. As with covariance matrices, the upstream variable should be placed along the columns and the downstream variable along the rows. For example, we define Y_m as the upstream variable and Y_p as the downstream variable, and thus the μ matrix will be

  

  A copath can only be traversed once within a given chain, and a chain must be legitimate prior to traversing the copath. In essence, a copath serves to connect two legitimate chains, thereby creating one longer chain. For example, COV (Y_p,Y_m) = V_YpμV_Ym

• Rules for path-tracing multiplication: Matrices are multiplied in the sequence in which they are encountered while tracing a path. When doing the path tracing, it is crucial to make sure the multiplied matrices are conformant to the matrix multiplication rules. Specifically, the upstream variable of the premultiplied matrix should be the same variable as the downstream variable of the post-multiplied matrix.

Although all matrices involved in the path tracing are of the same n×n dimensions (for example, they are all 2 × 2 matrices in a bivariate model), it is helpful to determine whether each matrix is diagonal, full, or symmetric. Distinguishing among these forms is important for verifying the correctness of mathematical expectations and interpreting their conceptual meaning. A n × n matrix will be full if the designated upstream and downstream variables are not the same variables in the model—for example, when the 2 × 2 covariance matrix is an off-diagonal submatrix of a larger symmetric 4 × 4 covariance matrix that includes all variables. Consider the following 4 × 4 symmetric phenotypic variancecovariance matrix between two paternal and two maternal phenotypic scores:

This matrix can be partitioned into four 2 × 2 covariance submatrices. The two V_y matrices along the diagonal represent the within-person phenotypic variance-covariance matrices between the paternal and two maternal phenotypes, respectively, and are symmetric. In contrast, the two off-diagonal mate submatrices are transpositions of one another and are full matrices, because their own off-diagonal elements correspond to conceptually different constructs. Moreover, given the equation for mate covariance, COV(Y_p,Y_m) = V_YpμV_Ym, and the fact that both V_Y matrices are symmetric, it follows that μ must also be a full matrix. This is intuitive: the assortment between males on trait 1 and females on trait 2 is not necessarily the same as the assortment between males on trait 2 and females on trait 1.

Parameters and Their Interpretations

Using multivariate path-tracing rules, we derived the expectations of all the variance and covariance listed in Table 3 in the Supplement. In this section, we outline the conceptual meaning of key estimated parameters.

Direct additive genetic effects are modeled through the path coefficients δ and a. Both are diagonal matrices, as pleiotropy between the two traits in the base population is captured in the off-diagonal elements of the k and j matrices, while the increase in additive genetic (co)variance induced by AM is quantified in the g, h, and i matrices. The coefficient δ indexes the additive genetic effects of haplotypic PGSs on Y , whereas a indexes residual additive genetic effects not captured by haplotypic PGSs. For almost all traits currently, a » δ. Information for estimating δ is easily estimable from the covariance between the haplotypic PGSs and Y, after correcting for the effects of AM and genetic nurture. Information for estimating a can come from the residual covariance between Y_o and Y_p/m after accounting for all other parameters estimated in the model, though because other factors besides additive direct and indirect genetic effects, VT, and AM can increase parent-offspring covariance, it may be preferable to obtain this parameter from external sources such as RDR regression (Young et al., 2018). The total additive genetic variance of a trait is then derived by tracing all paths that begin with δ and a and return to the phenotype (Y_p/m).

The k matrix represents the variance-covariance of the haplotypic PGSs in the base population (before AM). Similarly, j is defined as the genetic variance of the haplotypic latent genetic score (LGS) in the base population. The values of k depend on the scaling of the haplotypic PGSs. Below are several ways to fix the value of k_ii - the diagonal entries of trait i in the k matrix:

• , when the full PGS is standardized in the base population (which is only possible with simulated data).

• , when the full PGS is standardized in the current generation, which is the approach we recommend in real data.

• , when the haplotypic PGS is standardized in the current generation.

However, unlike k, j is the variance of a latent construct and can take any arbitrary value. The simplest choice is to set j_ii = k_ii and we recommend doing so because this allows the estimates of a and δ to be comparable. k₁₂ and j₁₂ are freely estimated parameters that represent the scaled genetic covariance between the two haplotypic PGSs (k₁₂) or the two haplotypic LGSs (j₁₂) , respectively, after accounting for AM and G-E covariance. In the version of the model described here, the coefficients δ and a for parents are constrained to equal those for offspring. This assumption may be violated when the measured phenotype reflects a different construct, or when its genetic architecture differs across age or cohort (e.g., substance use)—a point we return to in the Discussion. We also assume k₁₂ = j₁₂, implying that the correlation of allelic effects between traits is the same for variants captured and not captured by the PGS. This assumption is probably typically met, but could fail, for instance, if the genetic correlation between traits differs for rare versus common variants, given that common variant effects are currently better captured by PGSs.

The VT effect is modeled by the path coefficient matrix f . VT represents the direct influence of parental phenotypes on offspring phenotypes. Most often this reflects the impact of parental traits on offspring traits through parental behavior, though other pathways (e.g., in utero effects) are also possible. f is a full matrix: the effects of parental trait 1 on offspring trait 2 can be different than the effects of parental trait 2 on offspring trait 1. The estimation of f depends most strongly on the ratio (Balbona et al., 2021). Put simply, when the association of non-transmitted alleles with offspring outcomes grows relative to that of transmitted alleles, this pattern indicates VT, because the nontransmitted alleles can affect offspring only indirectly by shaping parental behavior and environments. By leveraging parental PGSs based on transmitted and non-transmitted alleles in this way, f can be estimated without confounding from shared parent-offspring genes. Embedding this approach within the SEM framework further reduces bias from AM and from genetic variance not captured by the PGS.

Genetic nurture refers to parental genetic effects—including non-transmitted alle- les—that influence offspring outcomes indirectly via the parental environment, a process arising through VT(Balbona et al., 2021, 2022; Kong et al., 2018). Our model quantifies genetic nurture through two mediated pathways, ϕ and ρ. Specifically, ϕ represents the effect of full PGSs (the sum of two haplotypic PGSs) on offspring outcomes as mediated through parental traits, while ρ represents the corresponding effect of full LGSs. These pathways are defined as:

Here, ϕ and ρ represent genetic nurture that stems from founder-generation additive genetic effects—the “pure” genetic nurture effect not distorted by AM. Because f is a full matrix, both ϕ and ρ are also full, allowing asymmetric cross-trait genetic nurture in which the genetic nurture effect of trait 1 on trait 2 can differ from that of trait 2 on trait 1. The ϕ and ρ are defined in a way that is independent of the influence of AM, so they target different estimands compared to other available methods, such as η and regression coefficients in regression-based methods (Kong et al., 2018; McAdams et al., 2023; Young et al., 2022).

As noted above, we model AM using a copath μ. AM refers to nonrandom mating, in which individuals preferentially pair with others who share similar or dissimilar traits at rates exceeding chance (Vandenberg, 1972). The copath was introduced as an elegant way to model AM in path diagrams (Cloninger, 1980), but has rarely been applied in multivariate family models—with only one exception we are aware of (Maes et al., 1997). In this paper, we formally lay out the multivariate copath rules for SEM and use them to model multivariate AM. Copaths provide a path tracing framework that links the expected covariance between parental phenotypes and their haplotypic PGSs to the observed data. In the present model, mate covariance is expressed as , allowing estimation of μ. Empirically, copath estimates correspond to mate phenotypic covariance scaled by phenotypic variance, given by .

AM increases the genetic variance in the population and the genetic covariance between relatives and mates. In our model, this is captured with the g, h, and i parameters. Specifically, g generally represents the increase in the genetic (co)variance of haplotypic PGSs between two traits under AM. We used different subscripts “t” and “c” to denote trans-individual (i.e., between mate) and cis-individual (within-person) haplotypic PGS covariances, respectively. For example, g_t is a full matrix representing the AM-induced covariance between maternal and paternal haplotypic PGSs. g_t is a full matrix for the same reason that μ is full. In contrast, g_c is a symmetric matrix representing (a) the AM-induced increase in the covariance between the two haplotypic PGSs within-person or, equivalently, (b) the AM-induced increase in their variance. Although g_c is the covariance matrix between two distinct variables, T_p/m and NT_p/m , and would be expected to be full according to our previously described rules, it must be symmetric for biological reasons. This symmetry arises from the randomness of haplotype formation during meiosis, as haplotypes are defined relative to the offspring genotype. Analogous to g, h_t and h_c represent increases in latent genetic covariance for trans- and cis- covariances, respectively. Finally, unlike the matrices g_c and h_c, i_c is the full within-individual covariance matrix between and . Because i_c models the covariance between two distinct sets of variables, the matrix is full. To illustrate this, consider an extreme scenario where trait 1 is determined entirely by its PGS and trait 2 entirely by its LGS; here, the term i_c21 will be non-zero while i_c12 would be zero.

G-E covariance rises as a function of VT and AM. We quantify this total G-E covariance using the matrices w and v, which are defined as the covariance between the familial environment () and the combined non-transmitted parental alleles ( + ). Unlike ϕ and ρ, which represent pure genetic nurture pathways, w and v capture the full G-E covariance—both within and across traits—that results from both AM and VT. Because they are a function of f , w and v are inherently full matrices. The choice between ϕ and w as the focal parameter depends on whether the research question pertains to a specific genetic nurture pathway or the total G-E covariance.

In the following sections, we present simulations evaluating the bias and precision of the parameter estimates, and then discuss the simulation results, the model’s assumptions, and recommendations for its application.

Methods

We conducted simulations to evaluate the performance of the Multivariate SEM- PGS model, focusing on two questions. First, we assessed bias by testing whether the model’s median estimates were systematically different from the true simulated parameters. Second, we evaluated the precision of these estimates as a function of sample size and PGS effect size. To address the first question, we simulated evolutionary processes across generations using an R implementation of the GeneEvolve software (Tahmasbi & Keller, 2017). Unlike methods that draw phenotypes directly from a multivariate normal distribution, GeneEvolve uses a forward-time simulation approach by generating a baseline population that evolves according to predefined mating and reproductive rules. Thus, GeneEvolve simulates the underlying causal processes instead of relying on model-implied covariance matrices, thereby providing an independent way to validate our derived equations and to examine how parameter estimates behave when model assumptions deviate from the true causal processes. To address the second question regarding estimate precision, we simulated two traits and their haplotypic PGSs directly from the model-implied multivariate normal distribution. The covariance matrix, Σ, for this distribution was derived from multivariate path-tracing from our model diagram. Because the forward-time simulations already established that the model produced unbiased estimates, we did not need to revisit bias here. Instead, this approach isolates estimation precision by avoiding the additional stochasticity inherent to forward-time simulators, which can inflate the variance of parameter estimates (Chen et al., 2025). All R scripts for this article are available on our GitHub repository (https://github.com/Xuanyu-Lyu/BiSEMPGS).

Forward-Time Simulation Design

We first explain the process for simulating genetic, environmental, and phenotypic scores in the baseline population. Each individual had two phenotypes, which were both standardized to have a total variance of 1. Heritability of each trait (the proportion of phenotypic variance (V_Y) attributable to genetic factors (V_G) and genetic effects in the base population were simulated by drawing m = 50 causal variants (CVs) from a binomial distribution with minor allele frequencies uniformly distributed as p ~ U(0.1, 0.5). The effect sizes of the CVs were drawn from a multivariate normal distribution ~ N([0,0],), where r_g represents pleiotropy (the genetic correlation) in the base population. We then computed the haplotypic PGS and haplotypic LGS from these variants, ensuring that each score had a variance of 1 (thereby using the first method of standardized PGSs where ) and that their genetic correlation aligned with the specified r_g in first generation. These scores were subsequently scaled using the δ and a matrices to achieve the desired proportions and , thereby generating the PGSs and LGSs for the base population (note that these were by definition uncorrelated in the base population). We then simulated residual environmental scores (ϵ) to account for residuals in the phenotypic variance unexplained by V_G . Phenotypic variance was scaled to be 1 in the based population and therefore the distribution of ϵ was ~ N([0,0], I — diag(V_G)). Phenotypic scores for each individual were obtained by summing their genetic and environmental scores.

After generating the base population, the GeneEvolve algorithm paired males and females such that their phenotypic correlation matched the specified mate correlation matrix (r_mate). Each pair of mates produced a number of offspring that followed a Poisson distribution with a mean of 2. We simulated genetic scores of offspring (including both PGS and LGS) using randomly sampled CVs from the mother and separately from the father. Family environmental scores were created for offspring by , reflecting the VT process. We created phenotypic scores of offspring by summing genetic scores, family environmental scores, and residual environmental scores (ϵ). The algorithms for mating and reproduction iterated through several generations. We used PGSs and phenotypic scores from the final generation as inputs for model fitting.

A key challenge in forward-time simulations is modeling cross-trait AM and its genetic consequences (Border & Malik, 2023). To impose the target mate correlation matrix (r_mate) across two traits, we first split the generation’s phenotypic data into males and females and constructed a four-dimensional template of correlated scores, X_sim = (m₁ , m₂, f₁, f₂), by drawing from N (0, Σ). The off-diagonal symmetric entries of Σ were set to the empirical correlations between the two traits in that generation, while the 2-by-2 block linking males and females was defined by the user-specified r_mate . We then computed Euclidean distance matrices between observed, standardized phenotypes {Y₁, Y₂} and their corresponding template vectors for males and females, and greedily paired individuals using a duplicate-removal/matching routine that minimized total distance, thereby aligning the realized cross-partner correlation structure with Σ. This template-guided matching roughly approach achieved both within-trait and cross-trait AM patterns specified by r_mate while preserving the observed phenotype distributions within each sex.Nevertheless, this approach consistently produced r_mate values slightly (e.g., < .01) below the target correlation in each generation. This discrepancy resulted in small deviations in the estimated μ and g_c values from their true population values (see details in Results). To address this limitation, we also fit the model using data generated directly from a covariance matrix derived from mathematical expectations.

Multivariate Normal Simulations

We used multivariate normal simulations to assess the precision of parameter estimates. These simulations were performed by sampling PGSs and phenotypic scores from the model-implied covariance matrix derived from Figure 1. This task was not straightforward because, while some parameters can be pre-defined, others are nonlinear functions of multiple parameters. Such dependencies arise from the recursive interrelationships among variables: for example, AM and VT alter G-E covariance across generations, which in turn modifies genetic variances, further feeding back into G-E covariance. Therefore, to find the model-implied covariance, we implemented an iterative algorithm that specified baseline values and then simulated how parameters evolved across generations under the evolutionary processes represented in our model. Our models assumed equilibrium—that is, a population state in which phenotypic and genetic variances remain constant across generations due to a balance between variance-increasing forces (AM and VT) and variance-reducing forces such as recombination (Bulmer, 1971). Thus, parameter values were retained once they stabilized across generations, indicating equilibrium. Specifically, in the multivariate SEM-PGS model (Table 4), f, r_mate, V_G, V_E, , and were fixed across generations, while all other parameters were algebraic functions of these pre-defined values and obtained from the iterative procedure. Using these equilibrium values, we then computed the full model-implied covariance matrix.

Table 4

We conducted both the bias analysis using forward-time simulation and the precision analysis using multivariate normal simulations under a range of conditions detailed in Table 4. To assess model precision, we systematically varied the sample size of trios (N ∈ {4k, 8k, 16k, 32k, 48k, 64k}) and the variance explained by the PGS for the first trait () across five levels. The variance explained by the second trait’s PGS (), the VT parameter (f), and the additive genetic variance in the base generation (V_G) were held constant across all conditions. Both the forward-time simulations and the analytical derivation of the expected equilibrium covariance were iterated for 20 generations to ensure that parameters influenced by AM and VT reached equilibrium. For each parameter combination, we generated 500 replicate datasets. The simulation code is available at https://github.com/Xuanyu-Lyu/BiSEMPGS.

Model Fitting and Estimation Analyses

We fit models using OpenMx with the NPSOL optimizer to address the non-linear constraints in the model (M. C. Neale et al., 2016). We used the median estimate across 500 models run on forward-time simulated data to evaluate bias, and the median absolute deviation (MAD) on 500 models run on multivariate normally distributed data to evaluate precision. We report median and MAD as primary statistics because they are robust to outliers, which complex SEM models with non-linear constraints can occasionally produce, even with large sample sizes. Based on our experience, approximately 3-4% of models yielded outlier estimates when N_trio = 8k, decreasing to around 1% when N_trio = 64k. Because the number of estimates in the bivariate SEM-PGS model is very large, we limited those presented to keep the paper manageable and accessible. Specifically, we report two sets of results: eight within-trait estimates for trait 1 and eight cross-trait estimates between traits 1 and 2 (Figure 2). We show the bias of parameter estimates as a function of N_trio and their precision as a function of both N_trio and .

Figure 2

The estimation of the latent genetic path coefficient a is one of the most challenging aspects of the SEM-PGS model, as the only information to estimate it comes from the covariance between and after accounting for the other estimated factors in the model. In empirical studies, researchers can choose to fix a by calculating genetic variance of the base population (V_G-Base) using the RDR or sibling regression approach (Visscher et al., 2006; Young et al., 2018) in their own dataset. To evaluate the effect of fixing a on the precision of other parameter estimates, we fitted each simulated MVN dataset in two ways: (1) with the latent genetic path coefficient a fixed to its true value, and (2) with a freely estimated. As a follow-up, we performed a sensitivity analysis investigating the bias introduced by fixing a to incorrect values. To do this, we fit the bivariate SEM-PGS model where a₁₁ was fixed across a range of ±0. 125 from its true value (≈ 0.7746) in increments of 0.025, with N_trio = 32,000 and = 4%.

Finally, we investigated the effect on parameter estimates when a univariate model is applied to a truly bivariate data-generating phenotypic architecture. We simulated two scenarios (Table 5): (1) a phenotypic architecture characterized by strong within-trait and weak cross-trait effects; and (2) a phenotypic architecture dominated by weak within-trait and strong cross-trait effects. In both scenarios, the genetic correlation in the founder generation was set to zero to simplify interpretation. To quantify the bias, we fit a univariate SEM-PGS model on trait 1. For each condition, we generated 500 replicate datasets with N_trio = 64,000.

Table 5

Results

Unbiasedness of the Parameters

Figures 2 and 3 show the median and MAD of Bivariate SEM-PGS model estimates from forward-time simulation as a function of N_trios . The results indicate that parameters were unbiased or nearly unbiased across all sample sizes, with the bootstrap tests suggesting some of them significantly deviated from their true value (see Supplement for p-values). For the within trait estimates (Figure 2), very small deviations in the median values of V_Y , a, μ, and g_c from their true values were likely due to the imperfect simulation of bivariate AM and genetic effects in the forward-time simulation algorithms. Similarly, the crosstrait estimates (Figure 3) were unbiased or nearly so, with slight deviations observed in μ and V_Y . All significant biases were small relative to the sampling error of the estimates, accounting for < 10% of the total variability in estimated parameters, regardless of sample size (see Appendix I). Similar patterns were observed across all other conditions of N_trios and . Using the alternative multivariate normal simulation approach with covariances based on the model math, no estimates were biased (no estimate median was statistically different from the true parameter values; Figures S1 and S2 in the Supplement).

Figure 3

Estimation Precision

In Figure 4, we show the change in MAD of within-trait effect estimates of trait 1 (the [1,1] entry of each parameter estimate matrix) as a function of N_tri_o when fitting the model with fixed a compared to estimating a, given = .04. As expected, larger sample sizes led to better precision for all the within-trait estimates. For most within- trait estimates, the MAD showed a sharp decline as N_trio increased from 4k to 32k and decreased more gradually thereafter. Fixing a elements to their true values resulted in more precise latent parameter estimates across the board except for g_c. Among the within- trait parameters, the latent G-E covariance v11 benefited most from fixing a, because it is a direct function of a₁₁ and f₁₁ . For other parameters, the precision gains from fixing a were modest compared to those achieved by simply increasing N_trio . The pattern of changes for trait 1 within-trait estimates ([1,1] entry) was qualitatively similar to that for trait 2 ([2,2] entry).

Figure 4

Figure 5 shows how within-trait parameter precision changes with at fixed N = 32k . Larger values yielded much more precise estimates of f and v, because more predictive PGSs provide more reliable information about both direct and indirect genetic effects, thereby reducing sampling error of estimates that depend on these effects. In contrast, estimates of μ, δ, w, and g_c showed little dependence on . The small fluctuations observed were minor relative to the scale of the y-axes and likely reflected sampling error. Consistent with Figure 4, fixing a improved the precision of all parameters except g_c, with the largest benefit for v when PGSs are weak. For example, when = .16, the gap in SEs for v₁₁ between fixed and estimated a was smaller than when =.01, illustrating that fixing a is most valuable when PGSs explain little variance. Patterns for trait 2 mirrored those for trait 1.

Figure 5

Figure 6 illustrates the MAD of cross-trait parameter estimates as a function of N_trio . Consistent with the within-trait findings, larger sample sizes led to improved estimation precision for all parameters, with the most pronounced gains occurring for sample sizes up to 32K trios. Unlike in the within-trait analysis, fixing the direct genetic effect (a) benefited all cross-trait estimates, though the magnitude of the improvement varied by parameter. This benefit was particularly substantial for the VT (f) and G-E covariance (w) parameters at smaller sample sizes. For these same parameters (f₂₁ and w₂₁), precision increased more rapidly with sample size when a was freely estimated, but the fixed-a model maintained lower MAD across all sample sizes. For the remaining parameters, increasing sample size had a comparable effect on precision regardless of whether a was fixed or estimated.

Figure 6

Figure 7 also reveals that fixing a consistently improved the precision of cross-trait estimates, but the effect of depended on whether a was fixed or freely estimated. Higher led to more precise estimates only when a was freely estimated; when a was fixed, the precision of most cross-trait estimates showed little dependence on the PGS’s predictive power. The only exception was the cross-trait latent G-E covariance parameter, v21 , which still benefited from higher . Overall, these results suggest that fixing a provides more precise estimates of cross-trait effects, particularly when PGSs are weak—a dynamic not observed for within-trait estimates.

Figure 7

Parameter Bias When Fixing a at Incorrect Values

Figures 8 and 9 show the results of a sensitivity analysis on within-trait and crosstrait parameter estimates, respectively, when a₁₁ was fixed to incorrect values. Overall, the VT (f) and G-E covariance (v) estimates were affected to varying degrees, whereas the estimates for μ, δ, V_Y, and g_c were minimally impacted, with their variation dominated by stochastic noise. Among the affected parameters, f₁₁, w₁₁, and w₂₁ were influenced to a smaller extent than the cross-trait parameters f₂₁, v₁₁, and v₂₁, with the latent GE covariance parameters (v₁₁ and v₂₁) most affected. Additionally, both the within- and between-trait familial variances (V_F11 and V_F12) exhibited negative relationships with the fixed value of a₁₁ , as both components are functions of f₁₁ and f₂₁ (Figure S4).

Figure 8

Figure 9

Estimate Bias Using a Univariate Model When the Phenotypic Architecture Is Truly Bivariate

Figure 10 illustrates the bias that arises when a univariate model for trait 1 is fit when the true data-generating process has weak within-trait but strong cross-trait VT and AM effects (Table 5, Condition 2). This model misspecification resulted in significant bias for nearly all parameters, with the exceptions of δ , a , and V_Y . Most notably, the within-trait VT parameter, f₁₁ , was severely overestimated: despite a true value of zero, its estimates were centered around 0.25. At the same time, the familial variance estimate for trait 1, V_F11 , was substantially underestimated (median = 0.2396 vs. the true value of 0.3665). These two results—an overestimate of f₁₁ alongside an underestimate of V_F11—are not contradictory. The true V_F11 arises entirely from cross-trait VT originating from trait 2. When both traits are modeled jointly, the bivariate model correctly identifies no VT from trait 1 to trait 1 but substantial VT from trait 2 to trait 1. In contrast, the univariate model misattributes the cross-trait influence as a spurious within-trait VT effect (f₁₁) while simultaneously underestimating the total familial variance for trait 1.

Figure 10

A related practice is to study cross-trait VT or genetic nurture by fitting a univariate model that pairs different traits across generations—for example, parental education and offspring health (Kong et al., 2018). However, this approach also fails to account for pleiotropy or for both within- and cross-trait AM and VT, and therefore yields biased estimates. Just as unmodeled AM can masquerade as VT in a univariate setting, strong cross-trait AM can masquerade as cross-trait VT when only a single trait is modeled in each generation—for example, parental substance use and offspring educational attainment.

Discussion

In this study, we introduced the multivariate SEM-PGS model. This represents the first multivariate framework to simultaneously estimate the effects of both within- and cross-trait AM, genetic nurture, VT, G-E covariance, and direct genetic effects. To support the development of this model, we described a multivariate path-tracing rule for copaths, complementing the multivariate path-tracing rules introduced by Vogler (1985). Through Monte-Carlo simulations and model fitting with OpenMx, we validated the accuracy of these path-tracing rules, demonstrating that both within- and cross-trait estimates were unbiased. Larger sample sizes and more powerful PGSs () improved estimation precision, with the steepest gains between 4k and 32k trios. Fixing a to its population value enhanced estimate precision for most estimates, with cross-trait ones benefiting more than within-trait ones and latent factors such as VT (f) and G-E covariance (w and v) benefiting more than estimates based on observed variables such as (δ and g_c). However, fixing a to incorrect values introduced bias for these latent parameters. Furthermore, we showed that analyzing data using a univariate model when the data-generating process was actually multivariate can introduce substantial bias to parameter estimates, especially when cross-trait effects are strong. In the following sections, we discuss the strengths and assumptions of the multivariate SEM-PGS model, interpret the results in greater detail, and offer best-practice guidance for its application.

A Novel Tool for Gene-Environment Interplay

By integrating parental and offspring phenotypes with transmitted and nontransmitted PGSs for multiple traits in a SEM framework, the multivariate SEM-PGS model introduces a novel approach to studying genetic and environmental processes within families. Unlike traditional family designs that rely on assumptions about the sources of phenotypic resemblance, this framework uses directly observable relationships between phenotypes and transmitted and non-transmitted PGSs to provide empirical evidence of VT and AM. Moreover, recent family-based models that utilize PGSs have estimated direct and indirect genetic effects for single traits (Balbona et al., 2021; Kong et al., 2018), but have not properly accounted for cross-trait influences. For instance, parental traits may often influence different offspring traits (i.e., cross-trait VT, such as parental education influencing offspring BMI) and parents may assort across different traits (e.g., across different psychiatric disorders). In this manuscript, we have demonstrated that within-trait estimates may also reflect unmodeled cross-trait effects, and not modeling these cross-trait influences, or modeling them in a univariate framework, can lead to incorrect inference.

Through forward-time simulations, we confirmed that the multivariate SEM-PGS model provides unbiased estimates of direct genetic effects (a and δ), VT effects (f), genetic nurture effects (ϕ and ρ), G-E covariance w and v, AM effects (μ), and other parameters when its assumptions are met. The forward-time simulation approach also provided an independent check on our mathematical derivations and path-tracing rules for copaths. Furthermore, we found that the precision of latent parameter estimates improved with both N_trio and ; we recommend N_trio > 32K and > 0.01 as rough lower bounds for reliable estimation, although these trade off, so smaller sample sizes can be used as grows and vice-versa. Nevertheless, this highlights two limitations of the current model. First, it requires large family-based datasets with genomic information; at present, only a few resources, such as the Norwegian Mother, Father, and Child Cohort Study (MoBa), meet these criteria. However, the number of such datasets is expected to grow in the coming decade, given the distinct advantages of family-based genomic data over individual-level data (Davies et al., 2024). Second, the model is currently limited to traits with sufficiently high . While this threshold is already met for many traits, the expansion of large-scale GWASs will likely increase the range of traits for which the model can be applied.

Furthermore, our analysis demonstrates that fitting a univariate SEM-PGS model to data with an underlying bivariate causal structure introduces bias in multiple estimates, particularly when there are strong cross-trait effects. For instance, a significant indirect genetic effect estimate for f in a univariate analysis of depressive symptoms could be an artifact from an unmodeled trait that is genetically correlated with depression. While these indirect effects are themselves real, the parameters must be interpreted with the caveat that they may represent integrated effects from multiple pathways. This same caution applies to the interpretation of multivariate SEM-PGS results, as it is practically infeasible to account for all potential inter-generational transmission pathways with a limited number of measured phenotypes. By the same token, our results suggest that attempting to estimate cross-trait effects using a univariate model can also lead to biased estimates, and such models should be interpreted with caution.

Building on our previous univariate model (Balbona et al., 2021), the present framework provides a more refined parameterization of genetic nurture. In the univariate setting, the G-E covariance (w) and genetic nurture (the effect of parental genes on offspring mediated through parental phenotype) are conceptually distinct but numerically identical. This equivalence arises because the same set of genes that influence the parental phenotype, and thereby affect the offspring phenotype through VT, also influence the offspring phenotype directly. This equivalence does not hold once cross-trait effects are introduced. For example, when parental trait 1 influences offspring trait 2 through VT, the genes that shape parental trait 1 are not identical to the genes directly influencing offspring trait 2. In such cases, genetic nurture still contributes to cross-trait G-E covariance, but it is no longer reducible to either of the within-trait G-E covariance terms or even to the cross-trait covariance term itself. In a multivariate context, therefore, G-E covariance and genetic nurture are related but not synonymous.

To capture this distinction, we introduce two parameters in the multivariate SEM-PGS model, ϕ and ρ, that explicitly represent genetic nurture pathways. These parameters separate genetic nurture effects arising from both within-trait and cross-trait influences, and they are defined so as to be independent of AM. This framework clarifies the interpretation of genetic nurture in multivariate models and differentiates our approach from earlier measures such as η in Kong et al. (2018) and the regression-based α estimates summarized in Young et al. (2022), both of which can be conflated with AM.

Model Assumptions

Like all statistical models, the current model relies on several assumptions regarding genetic effects, AM, and environmental influences to provide accurate estimates. First, as currently implemented when estimating a within the model itself rather than using external estimates to fix it, the model assumes that the residual parent-offspring covariance after

accounting for VT, AM, and direct PGS effects is purely due to additive genetic factors, and uses this residual covariance to estimate a. However, in truth, this residual covariance is unlikely to be purely genetic in origin; for example, shared environmental influences can also contribute. As a result, estimating a directly within the model can lead to biased estimates of other quanities. Unbiased estimates of a can instead be obtained from external methods such as Relatedness Disequilibrium Regression (RDR) (Young et al., 2018) or sibling-based regression (Visscher et al., 2006), including multivariate extensions of these approaches. For example, in trio datasets with parental genotypes and offspring phenotypes, RDR can be used to estimate founder-generation genetic variance, from which a can be derived by imposing an additional constraint on the total additive genetic variance in the base population: VA_Base = (a²+δ²)(a²+δ²+V_ϵ). Fixing a in this way reduces the standard errors of other parameters, enabling analyses in smaller samples and with weaker PGSs. However, for this approach to provide valid standard errors of estimates, the uncertainty in the external estimate of a should be incorporated into the SEM-PGS model. One strategy (a parametric bootstrap approach)is to draw a repeatedly from a normal distribution with mean equal to its point estimate and standard deviation equal to its standard error, refitting the model at each draw. This propagates uncertainty from the independent method into the SEM-PGS framework. Alternatively, this same principle could be implemented in a Bayesian framework by placing a prior on a informed by its external estimate. Our sensitivity analyses suggest that small misspecifications of a have limited impact on most estimates, but cross-trait VT and G-E covariance parameters are more vulnerable. It is important to note that, if parental phenotypes are included in the model, an additional latent factor (e.g., a shared environmental factor across both generations) should also be specified to absorb any residual parent-offspring covariance, ensuring that the model fits well. Alternatively, it is possible to simply omit parental phenotypes altogether when fixing a and still achieve unbiased estimates of all other parameters. Overall, future SEM-PGS models should aim to better estimate a internally (e.g., by incorporating sibling or extended relative information), but in the meantime, we believe that drawing on external approaches such as RDR will generally yield nearly unbiased estimates.

Second, the model assumes that genetic effects are consistent between parents and offspring. Specifically, it posits that δ_p/m = δ_o and a_p/m = a_o . The model also implicitly assumes a genetic correlation of unity between parental and offspring phenotypes. This assumption may be reasonable for traits such as height or weight, but may be violated in practice when social traits are measured at different ages or when the factors influencing the phenotype differ across cohorts. To address this limitation, we are working on extending the model by using informatio on additional relatives to allow parental genetic effects to differ from offspring genetic effects.

Third, the model assumes equivalence in genetic correlation between latent and observed genetic variables (k₁₂ = j₁₂). This assumption is generally plausible, since measured causal alleles are likely to resemble unmeasured (latent) alleles in their cross-trait genetic correlations. However, this assumption might be violated if common variants (more likely captured by PGSs) and rare variants differ in their contributions to cross-trait genetic correlations.

Fourth, the model assumes the absence of dominance or epistasis effects. This assumption is common in most statistical genetic models, as estimating dominance or epistasis at the level of PGSs is challenging and rarely attempted. However, since dominance and epistasis genetic effects are statistically orthogonal components to additive genetic effects, ignoring them is only likely to inflate V_ϵ, which alters the interpretation of V_ϵ but should not bias or otherwise alter interpretations of other parameter estimates.

Fifth, in its present form, it assumes that VT effects are identical for both father and mother (f_p = f_m). This assumption may not hold true for certain traits where VT from mothers differs from fathers. It would be simple to reparameterize the model to allow for sex-different VT paths, although doing so would reduce estimates’ precision.

Finally, the SEM-PGS model makes two assumptions about AM. One is that assortment is primarily phenotypic, even though in practice it may also arise through social homogamy or genetic homogamy (Robinson et al., 2017). Because these mechanisms are not mutually exclusive, misclassifying the source of AM can bias parameter estimates by leaving residual spousal similarity unaccounted for. For example, treating social homogamy as phenotypic assortment could lead to underestimation of environmental influences such as V_F and w. In our simulations, however, misspecifying social or genetic homogamy as phenotypic assortment produced only limited bias in the final parameter estimates. Furthermore, a major advantage of incorporating genomic data is that the type of AM can be tested empirically before fitting the full bivariate SEM-PGS model. This is done by comparing the observed spousal phenotypic correlation (g) to the correlation implied by haplotypic PGSs: if the latter exceeds the observed correlation, it suggests genetic homogamy, whereas if the observed correlation is larger, it implies social homogamy. The SEM-PGS framework has sufficient data to estimate these correlations and can be extended to model each contribution directly by adding latent AM factors. We have previously implemented such modifications in the Cascade model, which is based on extended twin family data, and adapting them here would be straightforward (Keller et al., 2009). Alternatively, a shared environmental factor between parents could be added to the model to soak up residual spousal covariance not captured by that implied from the spousal PGS covariance.

Another assumption is that AM has reached equilibrium, meaning that the level of assortment in the parental generation matches that in the offspring generation. Again, an advantage of using genomic data is that this assumption can be evaluated by comparing the within-trait increase in the covariance of within-individual haplotypic PGSs (g_c), which reflects AM in the grandparental generation and before, to the covariance of cross-individual haplotypic PGSs (.5()), which reflects AM in the parental generation. When disequilibrium is present, the mathematical expectations of the parameters can be modified accordingly, and we provide a derivation in the supplement for when AM is only present in the parental generation.

Conclusion

In summary, we introduced the multivariate SEM-PGS model, a flexible statistical framework that integrates within- and cross-trait direct genetic effects, genetic nurture, G-E covariance, VT, and AM into a single structural equation modeling approach. A key innovation is the development of multivariate path-tracing rules, which enable rigorous modeling of cross-trait AM. Simulations demonstrated that the model yields unbiased estimates and increasing precision with larger trio samples and more predictive polygenic scores, while a complementary strategy of borrowing information from external methods highlights its adaptability. We provide practical suggestions for fitting the model in the Supplement, but emphasize that the exemplar specification presented here should be adapted to the unique properties of different datasets. Although the framework rests on assumptions and requires relatively large family-genomic samples, it offers a general foundation that can be tailored for different data and questions. By providing tools to disentangle complex intergenerational pathways, the multivariate SEM-PGS model advances the study of how genetic and environmental influences combine to shape human traits and lays the groundwork for future extensions to more nuanced and biologically informed models.

Supplement

Parameter Expectations

This section provides a step-by-step guide to deriving the expected values of all estimated parameters. Throughout, symbols in the equations correspond to the matrices defined in Tables 2 and 3. To simplify notation, we write [N]T_p/m to indicate any of NT_p, T_p, NT_m, or T_m , where subscripts p and m refer to paternal and maternal haplotypes, respectively, and where NT stands for the haplotypic PGS not transmitted and T for the haplotypic PGS transmitted to offspring. The “L” in LNT and LT stands for “latent” (as opposed to the observed haplotypic PGSs). Finally, subscripts “LO” in, e.g., itLO stand for “latent to observed". For each parameter, we first show the relevant path-tracing chain (given in parentheses for convenience). For example, Y_p/m → [N]T_p/m denotes the path tracing starts with Y_p/m and finishes at [N]T_p/m , which means Y_p/m is the downstream variable and [N]T_p/m is the upstream variable. We then present the corresponding mathematical derivation. This approach is intended to make the logic of each step transparent, so that readers can see how expectations arise from simple path-tracing rules and so that the model can be altered depending on the data and questions at hand.

Assuming Equilibrium AM

We define the Ω parameters as “shortcuts” or shorthand covariances to facilitate the derivation of other parameter expectations. Specifically, they represent the covariance between haplotypic PGSs and parental phenotypes. We show both non-transposed and transposed forms, as each is useful for deriving subsequent expectations.

Similar to Ω, Γs are shortcut covariance between LGS and parental phenotypes.

With Ω and Γ defined, we can derive other expectations. AM induces directional covariance between the effects of all causal variants. We quantify this as g for the increase in genetic (co)variance for the haplotypic PGSs, h for the increase in genetic (co)variance for the haplotypic LGS, and i for the covariance between these two constructs. The t and c subscripts stand for “trans” (across mates) and “cis” (within-person) respectively.

Unlike g and h, the covariance between the the haplotypic PGSs and the haplotypic LGS, i, is a full matrix (the covariance between the PGS of trait 1 and the LGS of trait 2 is conceptually different than the covariance between the PGS of trait 2 and the LGS of trait 1) and so requires special attention regarding when it is transposed. For i_t , we adopt the convention (arbitrarily but consistently with μ in the model) that the paternal variable (L[N]T_p and [N]T_p) is downstream and the maternal (L[N]T_m and [N]T_m) variable is upstream, and denote the order in the i subscripts. For i_c , we adopt a convention that LGSs are the downstream and PGSs are the upstream.

The gene-environment (G-E) covariance matrices w and v are full because they are derived from f, which is itself a full matrix as a consequence of vertical transmission (VT). Conceptually, the covariance between family environment for trait 2 and PGS/LGS for trait 1 can differ from the covariance between family environment for trait 1 and PGS/LGS for trait 2.

Here, we note a minor correction to the equations for θ_NT and θ_LNT presented in the supplement of Balbona et al. (2021). In the previous paper, the covariance terms w and v were given a scalar coefficient of 2, whereas the correct coefficient is 1, as reflected in our equations above.

Assuming Disequilibrium AM

The following equations are parameter expectations in a model that has only one generationof AM, where AM occurs only in the parental generation but not before. Note that other types of disequilibrium AM are also possible, and we believe the model could be adjusted to incorporate these.

Supplementary Figures

Figure S1

Figure S2

Figure S3

Figure S4

Suggested Approach for Fitting the multivariate SEM-PGS Model

Given the complexity of the Multivariate SEM-PGS model, we have outlined a step-by-step procedure to guide researchers in its application. Below is a concise and clear workflow to ensure accurate implementation and interpretation of the model.

1. Prepare the genomic data. Begin by performing quality control (QC) to identify Mendelian errors and verify that the parent-offspring relatedness is approximately 0.5, ensuring accurate kinship specification. Next, separate the PGS of offspring and parents into four distinct haplotypic PGSs.

2. Compute for all traits. Ensure that the sample size (N) is sufficient relative to the observed by consulting supplementary tables or conducting simulations.

3. Scale the haplotypic PGSs. We recommend dividing each haplotypic PGS by the variance of the sum of the two haplotypic PGSs within an individual. With this approach, the scaling factors k and j are set equal to and , respectively. This parameterization is preferable to simply fixing the variance of each haplotypic PGS to 1/2, because it allows the haplotypic variances, and not just their covariances, to carry information about AM.

4. Determine the estimation of parameter a. Based on the trait and whether the residual parent-offspring covariance is likely to be due only to additive genetics, decide whether to estimate a within the model or fix it using external methods. In general, we recommend estimating a outside the model. An unbiased estimation of V_A-Base can be achieved through approaches such as Relatedness Disequilibrium Regression (RDR) (Young et al., 2018) or sibling-based methods (Visscher et al., 2006). Then, set VA_Base = (a² + δ²)(a² + δ² + V_E).

5. Check the PGS Predictability. Compare (Parental PGS Predictibility) with (Offspring PGS Predictibility) by running a regression Y ~ PGS + Covariates for both parental and offspring traits. If is substantially different from , consider modeling distinct values of a and δ for parents and offspring. Different δ values for parents and offspring can be modeled directly, as they come from distinct information. However, the model presented in Figure 1 can not identify both a_parent and a_offspring at the same time. Therefore, fixing at least one of a_parent or a_offspring using the approaches mentioned in the last bullet point is required for model identification.

6. Evaluate assortative mating equilibrium. Verify whether g_c (genetic covariance) approximates . If not, this suggests that AM has not reached equilibrium. Refer to the equations provided in Appendix I, which account for deviations from AM equilibrium.

7. Evaluate the mechanism of assortative mating. Compare the observed average g with its expected value under phenotypic AM. If inconsistencies are found, introduce latent AM factors to represent different AM mechanisms and conduct model comparisons to identify the best fit (Keller et al., 2009).

8. Constrain or free parental effects. Determine whether f_p (paternal effect) and f_m (maternal effect) can be constrained to equality based on statistical or theoretical grounds. If necessary, estimate these parameters independently.

9. Fit the model using OpenMx. Implement the model using an OpenMx script that aligns with the decisions made in previous steps. Utilize the MxTryhard function to explore multiple starting values, adjust feasibility and optimality tolerance settings, and increase the number of trials.

10. Evaluate model fit. Examine the OpenMx output for the following:

    1. Confirm that the OpenMx Status Code is 1 or 0, indicating successful convergence.

    2. Ensure that the estimated V_Y (variance of the outcome) aligns with the observed V_Y in the covariance matrix. If there are discrepancies, interpret parameter estimates with caution, and/or consider modifying the model to account for the discrepancies.

    3. Check for any parameter estimates that deviate greatly from the typical range of parameter values in an SEM model.

    If issues arise, adjust starting values or MxTryhard parameters and rerun the model. Persistent problems may require revisiting earlier steps. For additional support, researchers are encouraged to open an issue on the paper’s GitHub page (https://github.com/Xuanyu-Lyu/BiSEMPGS) or contact the corresponding authors via email.

11. Report and interpret results. Present the final estimates and provide a detailed interpretation of the findings, ensuring alignment with the model’s theoretical framework.

Appendix I

Descriptive statistics are provided for all parameters for the condition = .04 values and all sample sizes. The p-value from the Wilcoxon rank-sum test (W-R test) is used to compare the sample median with the true value of the parameter. Systematic variance refers to the proportion of total variance that can be attributed to the systematic bias arising from the imperfect simulation of assortative mating.

Effect size = .04 with Freely Estimated a

Data availability

The current manuscript is a method development study. All the data used in this study were simulated using computer software. All simulation code is available at https://github.com/Xuanyu-Lyu/BiSEMPGS.

Additional information

Funding

National Institutes of Health (MH130448)

• Matthew C Keller

• Xuanyu Lyu

National Institutes of Health (K99AA032540)

• Tong Chen