Genomic SEM in umx (umxGSEM)

In umx: Structural Equation Modeling and Twin Modeling in R

knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = FALSE   # stub: set TRUE when verifying locally with current OpenMx/umx
)
library(umx)
umx_set_auto_plot(FALSE)

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)

1. DWLS weights: (W = \mathrm{diag}(1/\mathrm{diag}(V))) after reordering (V) into OpenMx residual order (all variances, then free covariances).
2. OpenMx data: observedStats = list(cov = S, useWeight = W, asymCov = V) — never legacy type = "acov".
3. Identification: genetic scale prefers unit factor variance (free first loading), not unit loading.
4. 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).
5. 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):

• Box mirror: https://utexas.app.box.com/s/sounavy84gwygj0j2askcyaoo2ostbgu
• Script: https://ibg.colorado.edu/cdrom2021/Day06-nivard/GenomicSEM_practical/GenomicSEM_Practical.R

Toy (subset) SNPs — for smoke tests only

[umxGSEM_dl_tutorial_files()] copies the IBG practical 1k-SNP files (and EUR LD scores / HapMap3) into a cache directory. 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). ldsc on this list is pedagogical, not a result.
• For multivariate GWAS they still pair SNP files with a full-SNP premade covstruc — in umx that is data(Psych_LDSC), not ldsc rebuilt from 1k SNPs.
• Full 1000G MAF file is not the .subset file.

tut = umxGSEM_dl_tutorial_files()
snps = umxGSEM_sumstats(
  files = tut$sumstats,
  ref = tut$ref,
  trait.names = c("SCZ", "BIP", "MDD"),
  se.logit = TRUE
)
gwas = umxGSEM_GWAS(covstruc = Psych_LDSC, SNPs = snps,
  traits = c("SCZ", "BIP", "MDD"), maxSNPs = 50)

Full-fat SNPs — production

1. 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.
2. Munge to LDSC/HapMap3 format (SNP, A1, A2, effect, P, …) with [umxGSEM_munge()]. HapMap3 list: Zenodo https://zenodo.org/records/7773502/files/w_hm3.snplist.gz (gunzip to w_hm3.snplist). Do not use [umxGSEM_dl_tutorial_files()] for this.
3. LD scores for European-ancestry LDSC if you run ldsc locally: https://alkesgroup.broadinstitute.org/LDSCORE/eur_w_ld_chr.tar.bz2 (~500 MB → ~1.2 GB) and weights_hm3_noMHC.tgz (~150 MB). Ancestry must match. Skip if you already have a covstruc.
4. Run multivariable LDSC → covstruc (S, V, I, N, m) — same shape as Psych_LDSC.
5. 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.)
6. Trait order in munge / ldsc / sumstats / SEM-GWAS must stay identical.
7. 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)

See the toy chunk above (umxGSEM_dl_tutorial_files → umxGSEM_sumstats → umxGSEM_GWAS). Pass traits= on the LDSC object; do not V[1:6, 1:6] (wrong pairs).

• umxGSEM_dl_tutorial_files() — cache the ~45 MB workshop pack (not production HapMap3 / Alkes / 1000G plink).
• 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()

umx documentation built on Oct. 11, 2026, 5:07 p.m.