Following standard practice (Bogdan et al., 2018; Purcell et al., 2009), we first conducted a series of preliminary analyses to select the PRS p-value threshold that provided the best fit and maximized effect sizes. In these models, the log of HED frequency was regressed onto age, PRS, marital status, the two-way interactions, and the three-way interaction between age, PRS, and marital status. The selected PRS was standardized and used in all subsequent analyses. We then fit GLMMs with marital status and the selected PRS separately to examine the main effects of the predictors on HED. In each model, marital status or PRS and its interaction with age were included as predictors for HED. We then fit and interpreted the full model with marital status, the selected PRS, the two-way interactions, and the three-way interaction between age, PRS, and marital status. We used the GLIMMIX procedure in SAS and maximum likelihood estimator based on Laplace approximation for parameter estimation.