--- title: "Genomic SEM in umx (umxGSEM)" subtitle: "Psych_LDSC common-factor check against GenomicSEM" author: "Timothy C. Bates" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{Genomic SEM in umx (umxGSEM)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} 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](https://doi.org/10.1038/s41562-019-0566-x). --- ## Data: `Psych_LDSC` Five traits from multivariate LDSC: **SCZ**, **BIP**, **MDD**, **EA**, **INSOM**. ```{r load-data} 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 `=~`. ```{r fit-umxGSEM} m1 = umxGSEM( model = "g ~= SCZ + BIP + MDD + EA + INSOM", covstruc = Psych_LDSC, estimation = "DWLS" ) umxSummary(m1) ``` Optional explicit lavaan syntax (same model): ```{r fit-lavaan-syntax} 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`): ```{r tmx-data} 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: - Script: ### 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: ```text 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. ```{r gsem-sumstats-toy, eval = FALSE} 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): (`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 ```{r session, eval = TRUE} sessionInfo() ```