---
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: <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

From the practical (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).
- 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):

```{r gsem-sumstats-toy, eval = FALSE}
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

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`, …) — in umx: `umxGSEM_dl_RefList()` + `umxGSEM_munge()` (OpenMx C++ munge). HapMap3 list: Zenodo `w_hm3.snplist` (see `?umxGSEM_dl_RefList`).
3. **LD scores** for European-ancestry LDSC: `eur_w_ld_chr/` (Bulik-Sullivan / LDSC standard bundle; also linked from GenomicSEM install materials).
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)

```{r gwas-smoke, eval = FALSE}
# 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

```{r session, eval = TRUE}
sessionInfo()
```
