Genomic SEM in umx (umxGSEM)

Purpose

This vignette is a verification stub for genomic structural equation modeling (GSEM) in umx via umxGSEM().

It documents that a one-factor model on the built-in LDSC object Psych_LDSC recovers the same DWLS point estimates and standard errors (to printed precision) as GenomicSEM’s commonfactor() on the same input.

Rails, not vibes. umx builds OpenMx summary WLS with modern
mxData(numObs = N, observedStats = list(cov = S, useWeight = W, asymCov = V))
after reordering \(V\) into OpenMx residual order (var_*, poly_*_*).
Legacy OpenMx type = "acov" is not used or accepted.

Reference: Grotzinger et al. (2019). Genomic structural equation modeling. Nature Human Behaviour, 3, 513–525. doi:10.1038/s41562-019-0566-x.

Data: `Psych_LDSC`

Five traits from multivariate LDSC: SCZ, BIP, MDD, EA, INSOM.

data(Psych_LDSC, package = "umx")
str(Psych_LDSC, max.level = 1)
# names: V, S, I, N, m
# S is 5×5 genetic covariance; V is 15×15 sampling covariance of vech(S)
colnames(Psych_LDSC$S)

Slot	Role
`S`	Genetic covariance (heritabilities on diagonal)
`V`	Sampling covariance of non-redundant elements of `S`
`I`	LDSC intercepts
`N`	Effective sample sizes for each `S` element
`m`	SNP count used in LDSC

Fit a common factor with `umxGSEM`

Unit-variance identification (std.lv = TRUE by default). Path-style ~= is accepted as lavaan =~.

m1 = umxGSEM(
  model      = "g ~= SCZ + BIP + MDD + EA + INSOM",
  covstruc   = Psych_LDSC,
  estimation = "DWLS"
)
umxSummary(m1)

Optional explicit lavaan syntax (same model):

m1b = umxGSEM(
  model = "
    g =~ NA*SCZ + BIP + MDD + EA + INSOM
    g ~~ 1*g
  ",
  covstruc   = Psych_LDSC,
  estimation = "DWLS",
  std.lv     = FALSE
)

Inspect summary WLS data (modern names only — not acov):

tmx_show(m1, matrices = c("data.S", "data.V", "data.W"), report = "markdown")

Decimal match: umxGSEM vs GenomicSEM `commonfactor`

Verified against GenomicSEM commonfactor(Psych_LDSC, estimation = "DWLS") on the same Psych_LDSC object (unstandardized loadings and residual variances; unit variance on the factor).

Factor loadings

Path	GenomicSEM	umxGSEM	SE (both)
g → SCZ	0.4314	0.4314	0.0408
g → BIP	0.4235	0.4235	0.0397
g → MDD	0.1310	0.1310	0.0148
g → EA	0.0087	0.0087	0.0058
g → INSOM	0.0165	0.0165	0.0055

Residual variances

Residual	GenomicSEM	umxGSEM
SCZ	0.0583	0.0583
BIP	0.0834	0.0834
MDD	0.0836	0.0836
EA	0.1123	0.1123
INSOM	0.0457	0.0457

Standardized loadings (GenomicSEM / same model structure) are approximately SCZ 0.87, BIP 0.83, MDD 0.41, EA 0.03, INSOM 0.08 — SCZ and BIP dominate the genetic factor; EA and INSOM barely load (as expected from S).

Fit indices (do not require decimal match)

Omnibus \(\chi^2\), CFI, and related indices can differ across packages because of how \(r'Wr\) is scaled, how the independence model is built, and dummy numObs. Parameters and SEs are the cross-check. For GSEM model comparison in umx, prefer change in CFI / SRMR and nested tests via umxCompare() (see package docs).

Package	Example omnibus (illustrative)
umxGSEM	\(\chi^2(5) \approx 388.5\); fit \(r'Wr \approx 407.3\) with `numObs = 1`
GenomicSEM `commonfactor`	\(\chi^2(5) \approx 504\); CFI \(\approx 0.45\); SRMR \(\approx 0.16\)

Implementation notes (for maintainers)

DWLS weights: \(W = \mathrm{diag}(1/\mathrm{diag}(V))\) after reordering \(V\) into OpenMx residual order (all variances, then free covariances).
OpenMx data: observedStats = list(cov = S, useWeight = W, asymCov = V) — never legacy type = "acov".
Identification: genetic scale prefers unit factor variance (free first loading), not unit loading.
Numerics: elements of LDSC \(V\) are often \(10^{-4}\)–\(10^{-6}\); weights are large; tiny print digits can hide real structure (tmx_show auto-bumps digits for small matrices).
Local minima: DWLS can sit in a secondary basin; umxGSEM defaults to tryHard = "yes" and eigen-based starts.

Re-run verification (maintainers): fit with eval = TRUE above; optionally compare to GenomicSEM commonfactor when that stack is available. No umx dependency on GenomicSEM or lavaan GSEM helpers is required for users.

SNP summary data for multivariate GWAS (pipeline roadmap)

Adding SNPs (GenomicSEM’s sumstats → commonfactorGWAS / userGWAS) needs per-trait GWAS files plus a 1000 Genomes MAF reference. Built-in Psych_LDSC / Anthro_LDSC are LDSC-only (S, V, …); they do not ship SNP betas.

How this maps to our LDSC objects

umx object	Traits	Closest workshop object	SNP tutorial files
`Psych_LDSC`	SCZ, BIP, MDD, EA, INSOM	`PRAC_COV` (5-trait user model)	Workshop has no bundled EA/INSOM SNP subsets for the 5-trait path model
(3-trait subset of Psych)	SCZ, BIP, MDD only	`PSYCH_COV`	Yes — tutorial multivariate GWAS uses these three
`Anthro_LDSC`	BMI, WHR, CO, Waist, Hip, Height, IHC, BL, BW	`Anthro_LDSC.RData` in workshop materials	Full anthropometric GWAS from original papers / consortia (not in the 1k SNP toy pack)

Workshop practical (IBG / Grotzinger materials):

Toy (subset) SNPs — for smoke tests only

From the practical (order must match ldsc / sumstats):

SCZ_subset.txt
BIP_subset.txt
MDD_subset.txt
reference.1000G.subset.txt    # drastically reduced 1000G ref for MAF

Notes from the tutorial authors:

Subsets are on the order of ~1k SNPs (munge/ldsc demo only).
For multivariate GWAS they still pair SNP files with a full-SNP premade PSYCH_COV LDSC object — i.e. subset SNPs are for sumstats/commonfactorGWAS smoke tests, not for rebuilding LDSC from scratch with only 1k SNPs.
Full reference: see wiki; full 1000G MAF file is not the .subset file.

Example (GenomicSEM naming; umx will mirror later):

files       = c("SCZ_subset.txt", "BIP_subset.txt", "MDD_subset.txt")
ref         = "reference.1000G.subset.txt"
trait.names = c("SCZ", "BIP", "MDD")
se.logit    = c(TRUE, TRUE, TRUE)
# GenomicSEM::sumstats(files, ref, trait.names, se.logit = se.logit)
# commonfactorGWAS(covstruc = PSYCH_COV, SNPs = p_sumstats, estimation = "DWLS")

For umx, the matching 3-trait genetic cov structure is the SCZ/BIP/MDD block of Psych_LDSC (or a dedicated 3-trait LDSC when we ship one).

Full-fat SNPs — production

Download full GWAS summary statistics for each trait from the original consortia (PGC for SCZ/BIP/MDD, SSGAC for EA, etc.). Do not use “top hits only” dumps.
Munge to LDSC/HapMap3 format (SNP, A1, A2, effect, P, …) — in umx: umxGSEM_dl_RefList() + umxGSEM_munge() (OpenMx C++ munge). HapMap3 list: Zenodo w_hm3.snplist (see ?umxGSEM_dl_RefList).
LD scores for European-ancestry LDSC: eur_w_ld_chr/ (Bulik-Sullivan / LDSC standard bundle; also linked from GenomicSEM install materials).
Run multivariable LDSC → covstruc (S, V, I, N, m) — same shape as Psych_LDSC.
Prepare SNP file for SEM-GWAS (sumstats-equivalent): align alleles to a 1000 Genomes Phase 3 reference used for SNP variance / MAF.
Full ref download (GenomicSEM wiki):
https://utexas.app.box.com/s/vkd36n197m8klbaio3yzoxsee6sxo11v
(reference.1000G.maf.0.005.txt style; not the workshop .subset.)
Trait order in munge / ldsc / sumstats / SEM-GWAS must stay identical.
Binary traits: know whether SEs are logistic (se.logit = TRUE is common even when effects are ORs). Continuous traits: OLS = TRUE or standardized betas as appropriate.

Minimum columns before munge: rsID, A1 (effect allele), A2, effect (beta or OR), P. Sample size \(N\) (or effective \(N\)) is required for scaling.

umx API for SNPs (implemented)

# Toy files under inst/developer/GenomicSEM/ (or workshop Box pack)
dir = system.file("developer", "GenomicSEM", package = "umx")
snps = umxGSEM_sumstats(
  files = file.path(dir, c("SCZ_subset.txt", "BIP_subset.txt", "MDD_subset.txt")),
  ref   = file.path(dir, "reference.1000G.subset.txt"),
  trait.names = c("SCZ", "BIP", "MDD"),
  se.logit = TRUE
)

# Match SNP traits: pass full LDSC object + traits= (do not V[1:6,1:6] — wrong pairs)
gwas = umxGSEM_GWAS(covstruc = Psych_LDSC, SNPs = snps,
  traits = c("SCZ", "BIP", "MDD"), maxSNPs = 50)
head(gwas)

umxGSEM() — one structural model (string or existing umxRAM/mxModel).
umxGSEM_sumstats() — allele-align + continuous-scale betas for multivariate GWAS.
umxGSEM_GWAS() — per-SNP expand S/V, DWLS fit, return SNP table (est = factor ~ SNP).

Session

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] umx_4.80.0     OpenMx_2.22.11 rmarkdown_2.31
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.60           bslib_0.11.0       
#>  [4] ggplot2_4.0.3       htmlwidgets_1.6.4   visNetwork_2.1.4   
#>  [7] lattice_0.22-9      numDeriv_2016.8-1.1 quadprog_1.5-8     
#> [10] vctrs_0.7.3         tools_4.6.1         generics_0.1.4     
#> [13] stats4_4.6.1        curl_7.1.0          parallel_4.6.1     
#> [16] tibble_3.3.1        xts_0.14.2          pkgconfig_2.0.3    
#> [19] Matrix_1.7-6        MuMIn_1.48.19       RColorBrewer_1.1-3 
#> [22] S7_0.2.2            RcppParallel_6.0.0  lifecycle_1.0.5    
#> [25] compiler_4.6.1      farver_2.1.2        stringr_1.6.0      
#> [28] textshaping_1.0.5   mnormt_2.1.2        htmltools_0.5.9    
#> [31] sys_3.4.3           buildtools_1.0.0    sass_0.4.10        
#> [34] yaml_2.3.12         pillar_1.11.1       jquerylib_0.1.4    
#> [37] MASS_7.3-66         admisc_0.40         cachem_1.1.0       
#> [40] nlme_3.1-170        lavaan_0.7-2        polycor_0.8-2      
#> [43] tidyselect_1.2.1    digest_0.6.39       stringi_1.8.7      
#> [46] dplyr_1.2.1         maketools_1.3.2     cowplot_1.2.0      
#> [49] fastmap_1.2.0       grid_4.6.1          cli_3.6.6          
#> [52] magrittr_2.0.5      DiagrammeR_1.0.12   pbivnorm_0.6.0     
#> [55] scales_1.4.0        TTR_0.24.4          quantmod_0.4.29    
#> [58] otel_0.2.0          zoo_1.8-15          kableExtra_1.4.1   
#> [61] evaluate_1.0.5      knitr_1.51          viridisLite_0.4.3  
#> [64] rlang_1.3.0         Rcpp_1.1.2          xtable_1.8-8       
#> [67] glue_1.8.1          xml2_1.6.0          svglite_2.2.2      
#> [70] rstudioapi_0.19.0   jsonlite_2.0.0      R6_2.6.1           
#> [73] systemfonts_1.3.2