# fastGMA documentation (full text)

This file collects the fastGMA documentation (https://chaoning.github.io/fastGMA/) and the output of `fastgma -h <n>` for every mode in one plain-text file, for AI assistants and offline reading. Each page is given in its Markdown source, starting with its front matter (title, description). The `fastgma -h <n>` output at the end is the authoritative list of options and defaults.

================================================================
Page: https://chaoning.github.io/fastGMA/
================================================================
---
layout: full
homepage: true
disable_anchors: true
description: Fast genomic mixed-model analysis for GWAS and GP of complex traits
---

## fastGMA

`fastGMA` is a C++ toolkit for genomic mixed-model analysis of complex traits. It provides efficient methods for:

- genome-wide association studies (GWAS),
- variance component and heritability estimation, and
- genomic prediction.

## Modules

| Mode | Purpose | Guide |
| --- | --- | --- |
| `--make-grm`, `--process-grm` | Build additive, dominance and epistatic GRMs; merge, convert, invert and decompose GRMs | [GRM](documentation/03_GRM/) |
| `--uvlmm` | Univariate linear mixed model with one or more GRMs | [uvLMM](documentation/04_UVLMM/) |
| `--remma` | Rapid linear mixed model for non-additive analysis | [REMMA](documentation/05_REMMA/) |
| `--glmm` | Generalized linear mixed model for binary traits | [GLMM](documentation/06_GLMM/) |
| `--mvlmm` | Multivariate linear mixed model for several correlated traits | [mvLMM](documentation/07_MVLMM/) |
| `--longwas` | Longitudinal data analysis with random regression models | [lonGWAS](documentation/08_LONGWAS/) |
| `--fast-reml` | REML estimation of genetic, GxE and heteroscedastic noise variances | [fastREML](documentation/09_FASTREML/) |

`fastgma -h` lists all modes, and `fastgma -h <n>` prints the full option list of mode `n`.

## Inputs and outputs

- **Genotypes:** PLINK binary files (`--bfile <prefix>` for `.bed`, `.bim`, `.fam`) or a feature matrix (`--ffile <prefix>` for `.fmat`, `.bim`, `.fam`), for every module except fastREML, which reads only GRMs.
- **Phenotypes and covariates:** a whitespace-separated table with a header row and the individual ID in the first column. Traits, covariates and class variables are selected by column name.
- **GRMs:** three files that share one prefix, given as `--grm <prefix>`:
  - `<prefix>.grm.bin`: the lower triangle of the GRM including the diagonal, stored row by row ($$A_{11}, A_{21}, A_{22}, A_{31}, \ldots$$) as little-endian float32.
  - `<prefix>.grm.id`: one line per individual (`FID IID`), in the row order of the matrix.
  - `<prefix>.grm.N.bin`: the number of SNPs used for each pair of individuals, in the same layout as `.grm.bin`.

  GRMs are built with `fastgma --make-grm` (see [GRM](documentation/03_GRM/)). This is the GCTA binary GRM format, so a GRM built by `fastgma` can be used directly in GCTA, and a GRM built by GCTA can be used directly in every `fastgma` module.

## Quick start

```bash
# 1. Additive GRM
fastgma --make-grm --bfile test --out test

# 2. Mixed-model GWAS with heritability estimation
fastgma --uvlmm --grm test --bfile test \
  --data pheno.txt --trait height --covar age --class sex \
  --out height
```

See [Installation](documentation/02_installation/) to download `fastgma`.

## Contact

If you have questions, please open an issue on [GitHub](https://github.com/chaoning/fastGMA/issues) or email ningchao91@gmail.com.


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/02_installation/
================================================================
---
layout: page
permalink: /documentation/02_installation/
title: Installation
description: ~
order: 2
---

`fastGMA` is a C++ command-line tool for genomic mixed-model analysis. It is distributed as a prebuilt 64-bit Linux executable.

## 1. Prebuilt Linux executable

A prebuilt 64-bit Linux executable is available here:
[**fastGMA Linux Executable**](https://github.com/chaoning/fastGMA/raw/refs/heads/main/app/linux/fastgma)

After downloading it, make it executable and check that it runs:

```bash
chmod +x fastgma
./fastgma -h
```

This is the quickest option if you only need to run the program. The executable is statically linked against MKL and the Intel OpenMP runtime, so it needs no other libraries on the 64-bit Linux systems it has been tested on. MKL's static libraries still call `dlopen`, so check the binary on your target system. It sets `KMP_AFFINITY=disabled` unless you set `KMP_AFFINITY` yourself, which avoids an affinity error of the Intel OpenMP runtime on some systems.

## 2. Quick validation

Confirm that the executable starts correctly:

```bash
./fastgma -h
```

For an end-to-end test, run the example commands on each module page with the data in `example/`.

## 3. Troubleshooting

- If the executable aborts with an Intel OpenMP affinity error, run it with `KMP_AFFINITY=disabled`.
- If it fails on your system, please open an issue with the error message and the output of `uname -a`.


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/03_GRM/
================================================================
---
layout: page
permalink: /documentation/03_GRM/
title: GRM
description: Build and process genomic relationship matrices
order: 3
---

The GRM module builds the genomic relationship matrices (GRMs) used as random-effect covariance matrices by `UVLMM`, `REMMA`, `fastREML` and the other mixed-model modules. It also post-processes existing GRMs. The files it writes use the GCTA format, so a GRM built by `fastGMA` can be used in GCTA, and a GRM built by GCTA can be used in `fastGMA`.

There are two entry points:

| Command | Short form | Purpose |
| --- | --- | --- |
| `fastgma --make-grm [options]` | `-m` | Compute a GRM from genotypes or a feature matrix |
| `fastgma --process-grm [options]` | `-d` | Merge, convert, invert, decompose or filter an existing GRM |

As elsewhere in `fastGMA`, a GRM is referred to by its prefix: `--grm output/test` means the files `output/test.grm.bin`, `output/test.grm.id` and so on.

## 1. Building a GRM (`--make-grm`)

### GRM types

Choose the matrix with `--grm-type`. Let $$x$$ be the genotype dosage (0/1/2 copies of the A1 allele), $$p$$ the A1 frequency and $$q = 1 - p$$.

| `--grm-type` | Matrix | Genotype coding for each SNP | GCTA equivalent |
| --- | --- | --- | --- |
| `add` (default) | Additive | $$x - 2p$$ | `--make-grm` |
| `dom` | Dominance, breeding-value ("statistical") coding | $$0 \to -2p^2$$, $$1 \to 2pq$$, $$2 \to -2q^2$$ | `--make-grm-d` |

Under Hardy–Weinberg equilibrium the dominance code is uncorrelated with the additive code, so fitting the `add` and `dom` GRMs together splits the genetic variance into the classical additive and dominance variances.

The scaling is set with `--code-type`:

| `--code-type` | Additive GRM | GCTA equivalent |
| --- | --- | --- |
| `2` (default) | each SNP standardized by its own SD, then averaged: $$A_{jk} = \frac{1}{N_{jk}} \sum_i z_{ij} z_{ik}$$, with $$z = (x - 2p)/\sqrt{2pq}$$ | `--make-grm` |
| `1` | VanRaden-type: $$A_{jk} = \sum_i (x_{ij} - 2p_i)(x_{ik} - 2p_i) \big/ \sum_i 2p_i q_i$$ | `--make-grm-alg 1` |

Both additive scalings are special cases of `--alpha`, which assumes that the variance of a SNP effect is proportional to $$(2pq)^\alpha$$. With $$v_i = 2p_i q_i$$, the additive GRM is

$$A_{jk} = \frac{\sum_i v_i^{\alpha} (x_{ij} - 2p_i)(x_{ik} - 2p_i)}{\sum_i v_i^{1+\alpha}},$$

so that the expected diagonal is 1 (for `--ffile` input see "Inputs" below). `--alpha -1` gives `--code-type 2` and `--alpha 0` gives `--code-type 1`, with identical output files. Under $$\alpha = -1$$ every SNP explains the same share of the genetic variance; values above $$-1$$ shift the variance toward common SNPs and values below $$-1$$ toward rare SNPs. Speed et al. (2017, *Nature Genetics* 49:986–992) recommended $$\alpha = -0.25$$ for human complex traits. To choose $$\alpha$$ for your data, build GRMs over a grid of values (for example −1, −0.75, −0.5, −0.25, 0) and compare the REML log-likelihoods. `--alpha` cannot be combined with `--code-type` and applies only to `add`. With missing genotypes, each element is rescaled by $$M/N_{jk}$$ as for the other scalings (see "Quality control and missing genotypes" below).

For the dominance GRM the two `--code-type` choices apply, using the variance of the dominance coding instead of $$2pq$$. An epistatic GRM (for example additive-by-additive) is built from existing GRMs with `--process-grm --make-epis` (see Section 3).

### Inputs

Provide exactly one of:

| Option | Input | Files read |
| --- | --- | --- |
| `--bfile <prefix>` | PLINK binary genotypes | `<prefix>.bed`, `<prefix>.bim`, `<prefix>.fam` |
| `--ffile <prefix>` | Any numeric feature-by-sample matrix (gene expression, methylation, dosages, ...) | `<prefix>.fmat`, `<prefix>.bim`, `<prefix>.fam` |

With `--ffile`:

- `<prefix>.fmat` has one row per feature and one whitespace-separated column per sample, in `.fam` order. It has no header and no row names.
- `<prefix>.bim` needs one line per feature; only the line count is used.
- Each feature is centred by its mean over the observed values. With the default `--alpha -1` (`--code-type 2`) it is also scaled by its standard deviation, so the GRM does not depend on the units of the features. Negative values are fine.
- `--code-type 1` and `--alpha` use the observed variance of feature $$i$$ as $$v_i$$ in place of $$2p_i q_i$$. Multiplying one feature by a constant $$c$$ multiplies its contribution to the GRM by $$c^{2(1+\alpha)}$$, so only $$\alpha = -1$$ is independent of the units. Use other values only when all features are on a common scale (for example log-expression from one assay). $$\alpha$$ then sets how strongly features with a large variance dominate the GRM; unlike for SNPs, it has no population-genetic interpretation. With complete data the mean diagonal is exactly 1 for every $$\alpha$$.
- Missing values are written as any of the `--missing-geno` tokens (default `NA NaN nan N/A ...`).
- `--maf` does not apply, and only `--grm-type add` is supported.

### Selecting individuals and SNPs

| Option | File or value | Effect |
| --- | --- | --- |
| `--keep <file>` | one individual per line: `IID`, or `FID IID` | use only these individuals |
| `--remove <file>` | same format | leave these individuals out |
| `--extract <file>` | SNP IDs in the first column | use only these SNPs |
| `--exclude <file>` | SNP IDs in the first column | leave these SNPs out |
| `--chr <list>` | chromosomes, e.g. `1-22`, `1,3,5` or `2 7 19` | use only these chromosomes |
| `--not-chr <list>` | same format | leave these chromosomes out |

A line with one column is matched to the IID; a line with two or more columns is matched to FID and IID, as in PLINK. A leading `chr` in chromosome names is ignored. The options can be combined: individuals are kept and then removed, and SNPs are extracted, excluded and then restricted by chromosome. The output keeps the `.fam` order of the individuals. Allele frequencies and missing rates are computed in the selected individuals, so the result is identical to running on a file subset with `plink --keep-allele-order --make-bed`. The SNP options require `--bfile`.

A leave-one-chromosome-out (LOCO) GRM for chromosome `c` is built with `--not-chr c`, and the GRM of one chromosome with `--chr c`.

### Allele frequencies

By default `p` is estimated from the selected individuals. `--read-freq <file>` uses given A1 frequencies instead, for example from a reference population, so that GRMs built on different subsets of individuals are on the same scale. The file needs a header line with the columns `SNP`, `A1` and `FREQ` (other columns are ignored), which is the format written by `--write-freq`:

```text
CHR	SNP	A1	A2	FREQ	NOBS
0	rs1155740	T	C	0.27383863	409
```

`FREQ` is the frequency of the allele in the `A1` column. If that allele is the `.bim` A2, the frequency is converted to `1 − FREQ`. SNPs that are absent from the file, or whose `A1` matches neither `.bim` allele, are excluded and counted in the log. The frequencies are used for centring, scaling and the `--maf` filter; the missing-rate filter still uses the analysed individuals. With the same frequencies, the GRM is identical to GCTA `--update-freq`.

`--write-freq` writes `<out>.freq` with one line per selected SNP: the A1 frequency used (estimated or read) and `NOBS`, the number of individuals with an observed genotype. With `--npart`, only part 1 writes the file. Both options require `--bfile`.

### Quality control and missing genotypes

A SNP (or feature) is excluded when any of the following holds:

- its minor allele frequency $$\min(p, 1 - p)$$ is below `--maf` (BED input only)
- its missing rate is above `--missing-rate`
- it has zero variance

Missing genotypes are handled as in GCTA. A missing call contributes nothing to the cross-product, and each pair of samples is scaled by `N_jk`, the number of SNPs observed in both samples. These counts are written to `.grm.N.bin`.

> **Note.** Because each element has its own denominator, a GRM built from data with missing genotypes is not guaranteed to be positive semi-definite (GCTA has the same property). With complete data it always is. If you need to invert such a GRM, use a larger `--val` with `--process-grm --make-inv`.

### Options

| Option | Default | Meaning |
| --- | --- | --- |
| `--bfile` / `--ffile` | — | Input prefix (exactly one is required) |
| `--out` | — | Output prefix (required) |
| `--grm-type` | `add` | `add` or `dom` |
| `--code-type` | `2` | `2`: per-SNP standardization; `1`: VanRaden-type scaling |
| `--alpha` | `-1` | Additive GRM with SNP-effect variance proportional to $$(2pq)^\alpha$$; excludes `--code-type` |
| `--maf` | `0.01` | Minimum minor allele frequency (BED input) |
| `--missing-rate` | `0.05` | Maximum missing rate per SNP/feature |
| `--keep` / `--remove` | — | Individuals to use / leave out (see above) |
| `--extract` / `--exclude` | — | SNPs to use / leave out (BED input) |
| `--chr` / `--not-chr` | — | Chromosomes to use / leave out (BED input) |
| `--read-freq` | — | A1 frequencies to use instead of estimating them (BED input) |
| `--write-freq` | off | Write `<out>.freq` (BED input) |
| `--val` | `0` | Value added to the diagonal of the output |
| `--npart <n> <i>` | `1 1` | Compute only part `i` of `n` of the lower triangle (see below) |
| `--npart-snp` | `20` | Number of chunks the SNPs are read in (exactly this many, unlike the association modules); raise it to lower memory use |
| `--missing-geno` | `NA NaN nan N/A ...` | Tokens treated as missing in `--ffile` input |
| `--threads` | `10` | Number of threads |

The full genotype matrix is never held in memory: the SNPs are read in `--npart-snp` chunks. Memory is dominated by the GRM block (about $$8n^2$$ bytes for $$n$$ individuals) plus one chunk (about $$8nm/S$$ bytes for $$m$$ SNPs and $$S$$ = `--npart-snp` chunks), so raise `--npart-snp` when $$m$$ is large. For very large samples, split the lower triangle into row blocks with `--npart` and merge them afterwards:

```bash
for i in 1 2 3 4; do
  fastgma --make-grm --bfile test --npart 4 $i --out output/test
done
fastgma --process-grm --merge --grm output/test --npart 4
```

Each job writes `output/test.4_<i>.grm.*` and its own log `output/test.4_<i>.log`. The merge step concatenates the `.grm.id`, `.grm.bin` and `.grm.N.bin` blocks into `output/test.grm.*` and copies `.grm.N` from the first block. The result is identical to a single full run.

## 2. Output files

| File | Content |
| --- | --- |
| `<out>.grm.bin` | Lower triangle including the diagonal, row by row ($$A_{11}, A_{21}, A_{22}, A_{31}, \ldots$$), little-endian float32 |
| `<out>.grm.id` | One line per sample: `FID<TAB>IID`, in matrix row order |
| `<out>.grm.N.bin` | Per-pair SNP counts `N_jk` in the same layout as `.grm.bin` (float32). GCTA needs this file for `--mgrm` merging, `--grm-cutoff`, `--make-bK-sparse` and `--adj-grm` |
| `<out>.grm.N` | Text file: number of SNPs used and the scale factor |
| `<out>.freq` | With `--write-freq`: `CHR SNP A1 A2 FREQ NOBS` for every selected SNP |

These are exactly the files GCTA writes (plus the text `.grm.N`).

Every module that reads a GRM accepts both value widths of `.grm.bin` (float32 as written now, or float64 from older fastGMA versions; the width is detected from the file size) and a `.grm.id` with one column (IID) or two columns (FID, IID). Samples are matched by IID (the second column, or the only column).

Each run also writes a log. It is appended to `<out>.log`, so several steps that share one `--out` prefix collect in one file. A partitioned run (`--npart n i`) logs to `<out>.<n>_<i>.log`. Set the environment variable `FASTGMA_LOG_LEVEL=debug` for more detail.

Reading a GRM in Python:

```python
import numpy as np

ids = [line.split()[-1] for line in open("output/test.grm.id") if line.strip()]
n = len(ids)
tri = np.fromfile("output/test.grm.bin", dtype=np.float32)
assert tri.size == n * (n + 1) // 2

A = np.zeros((n, n))
A[np.tril_indices(n)] = tri
A = A + A.T - np.diag(np.diag(A))
```

## 3. Processing an existing GRM (`--process-grm`)

Pass `--process-grm` together with exactly one operation:

| Operation | What it does | Main output |
| --- | --- | --- |
| `--merge` | Concatenate the `--npart` blocks of a partitioned run | `<grm>.grm.bin/.id/.N.bin/.N` |
| `--pca` | Leading eigenvectors of the GRM | `<out>.eigenval`, `<out>.eigenvec` |
| `--make-epis` | Elementwise product of two GRMs (e.g. additive × additive) | `<out>.grm.bin/.id` |
| `--reformat` | Convert to another storage format (`--out-fmt`) | see the format table below |
| `--make-inv` | Inverse of the GRM after adding `--val` to the diagonal (LAPACK Cholesky) | inverse in the `--out-fmt` format, plus `<out>.grm.id` |
| `--out-close` | Greedily select individuals to drop until no pair exceeds `--cut-value` | `<out>.grm.pruned_iid` |
| `--group` | Split samples into connected groups of relatives (pairs above `--cut-value`) | `<grm>.grm.group`, `<grm>.grm.group.size`, `<grm>.grm.sp.bin` |

| Option | Default | Meaning |
| --- | --- | --- |
| `--grm` | — | Input GRM prefix |
| `--out` | — | Output prefix |
| `--npart` | — | Number of partitions (`--merge`); must be greater than 1 |
| `--prefix1`, `--prefix2` | — | The two input GRMs (`--make-epis`); their `.grm.id` must have the same order |
| `--out-fmt` | `2` | Output format for `--reformat` / `--make-inv` (table below) |
| `--val` | `0.001` | Diagonal value added before `--make-inv` (this default is separate from the `--make-grm` default of `0`) |
| `--num-pc` | `20` | Number of leading principal components written by `--pca`; `0` writes all |
| `--cut-value` | `0.05` | Relatedness threshold for `--out-close` / `--group` |
| `--sparse` | off | The input is a sparse `.grm.sp.bin` (`--reformat`, `--pca`, `--out-close`; not `--make-epis`, `--make-inv` or `--group`) |
| `--threads` | `10` | Number of threads |

Output formats for `--out-fmt`:

| `--out-fmt` | File | Layout |
| --- | --- | --- |
| `0` | `<out>.grm.matrix` | Full $$n \times n$$ text matrix |
| `1` | `<out>.grm.index_triplet` | `row col value` per line, lower triangle, 1-based indices |
| `2` | `<out>.grm.iid_triplet` | `IID_row IID_col value` per line, lower triangle |
| `3` | `<out>.grm.bin` (dense) or `<out>.grm.sp.bin` (sparse) | Dense: GCTA float32 lower triangle. Sparse: records of `int64 row, int64 col, float64 value`, 0-based |

Notes:

- **`--merge`.** The blocks are read from `<grm>.<n>_<i>.grm.*` for `i = 1..n` and the merged GRM is written to `<grm>.grm.*`. `.grm.N.bin` is merged when the blocks have it. `--out` only names the log file (without it, the log goes to the console only).
- **`--pca`.** The files follow the GCTA layout. `<out>.eigenval` has one eigenvalue per line, largest first. `<out>.eigenvec` has one line per sample in `.grm.id` order with the columns `FID IID PC1 PC2 ...` (tab-separated, no header line). With a one-column `.grm.id`, the IID is written in both ID columns. `--num-pc 20` (default) computes only the 20 leading eigenpairs, which is faster than a full decomposition; `--num-pc 0` writes all of them.
- **`--make-epis`.** The two GRMs must have the same `.grm.id` order. The output is float32 `<out>.grm.bin` plus `<out>.grm.id`; `.grm.N` is the concatenation of the two input `.grm.N` files when both exist. No `.grm.N.bin` is written.
- **`--group`.** The output goes next to the input GRM, not to `--out`. Besides the group tables, it writes a sparse GRM (`.grm.sp.bin`) that keeps only within-group pairs.
- **`--make-inv`.** The inverse is computed with a LAPACK Cholesky factorization (`dpotrf`, `dpotri`). A GRM centred by sample allele frequencies has rank at most $$n - 1$$, so it cannot be inverted without a diagonal value. That is why `--val` defaults to `0.001` here, while `--make-grm` itself adds nothing. With missing genotypes the GRM can have slightly negative eigenvalues (Section 1). If the matrix is not positive definite, `--make-inv` stops, reports the smallest eigenvalue and asks to rerun with `--val` above a given value. With `--out-fmt 3` the inverse is stored as float32, like every `.grm.bin`.

## 4. Compatibility with GCTA

| fastGMA | GCTA |
| --- | --- |
| `--make-grm` | `--make-grm` |
| `--make-grm --code-type 1` | `--make-grm-alg 1` |
| `--make-grm --grm-type dom` | `--make-grm-d` |
| `--keep`, `--remove`, `--extract`, `--exclude`, `--chr` | same options |
| `--read-freq` | `--update-freq` |
| `--process-grm --merge` | concatenating `--make-grm-part` outputs |
| `--process-grm --pca` | `--pca` (same `.eigenval` / `.eigenvec` layout) |

With the same SNPs, the values agree with GCTA to float32 precision (largest difference about 1e-7), with or without missing genotypes.

- **GCTA GRMs in fastGMA.** Pass the GCTA prefix directly, for example `--grm gcta_out`.
- **fastGMA GRMs in GCTA.** `--reml`, `--mlma`, `--pca`, `--grm-cutoff`, `--make-bK-sparse` and `--mgrm` all read the files directly.
- **Files from older fastGMA versions** (float64 `.grm.bin`, one-column `.grm.id`) are still read by every module (see Section 2).

Differences to keep in mind when comparing the two tools:

| Item | fastGMA | GCTA |
| --- | --- | --- |
| Default `--maf` | `0.01` | `0` (only monomorphic SNPs removed) |
| Default missing-rate filter | `0.05` | none |
| Chromosomes used | all SNPs in the input | autosomes only by default |
| Partition file names | `<out>.<n>_<i>.grm.*` | `<out>.part_<n>_<i>.grm.*` |
| X chromosome GRM | not supported | `--make-grm-xchr` |
| Feature-matrix input (`--ffile`) | supported | not available |

To reproduce a GCTA GRM exactly, use `--maf 0` (or the same `--maf` in both tools), `--missing-rate 1`, and only autosomal SNPs (`--chr 1-22` for human data).

**LDAK.** With the same SNPs and complete genotypes, `--alpha` $$\alpha$$ gives the kinship matrix of LDAK `--calc-kins-direct --ignore-weights YES --power` $$\alpha$$ up to a constant factor: LDAK divides by the observed mean diagonal, so its matrix has mean diagonal exactly 1, whereas fastGMA (like GCTA) divides by the expected value $$\sum_i v_i^{1+\alpha}$$. After dividing the fastGMA matrix by its mean diagonal, the two agree to float32 precision. With missing genotypes they also differ in the treatment of missing calls (LDAK uses one denominator for all pairs instead of `N_jk`).

## 5. Examples

### Additive GRM for uvLMM

```bash
fastgma --make-grm --bfile test --out output/test
fastgma --uvlmm --grm output/test --bfile test \
  --data test_simu_pheno.txt --trait pheno --out output/test_uvlmm
```

### Leave-one-chromosome-out GRMs

```bash
for c in $(seq 1 22); do
  fastgma --make-grm --bfile test --not-chr $c --out output/test_loco$c
done
```

`output/test_loco<c>` is used when testing the SNPs of chromosome `c`.

### Choosing α

Build a GRM for each value on a grid and compare the REML log-likelihoods of the null model. Without `--bfile`, `--uvlmm` fits only the null model:

```bash
for a in -1 -0.75 -0.5 -0.25 0; do
  fastgma --make-grm --bfile test --alpha $a --out output/alpha_$a
  fastgma --uvlmm --grm output/alpha_$a --data test_simu_pheno.txt --trait pheno --out output/alpha_$a
done
grep "logL (REML" output/alpha_*.log
```

All models have the same fixed effects and the same number of variance components, so the log-likelihoods can be compared directly; the highest one gives the best-supported $$\alpha$$. For the example data:

| $$\alpha$$ | `G` | `residual` | Heritability | REML logL |
| --- | --- | --- | --- | --- |
| $$-1$$ | 0.2204 | 0.6551 | 0.2533 | −187.9204 |
| $$-0.75$$ | 0.2221 | 0.6534 | 0.2553 | −187.9097 |
| $$-0.5$$ | 0.2237 | 0.6518 | 0.2571 | −187.8994 |
| $$-0.25$$ | 0.2251 | 0.6504 | 0.2587 | −187.8895 |
| 0 | 0.2263 | 0.6491 | 0.2601 | −187.8802 |

With 427 individuals the log-likelihood changes by only 0.04 over the grid, so these data cannot distinguish the values of $$\alpha$$; estimating $$\alpha$$ needs large samples.

### Additive, dominance and epistatic GRMs for a multi-GRM model

```bash
fastgma --make-grm --grm-type add --bfile test --out output/add
fastgma --make-grm --grm-type dom --bfile test --out output/dom
fastgma --process-grm --make-epis --prefix1 output/add --prefix2 output/add --out output/axa

fastgma --uvlmm --model multi-grm --mgrm output/add output/dom output/axa \
  --bfile test --data test_simu_pheno.txt --trait pheno --out output/test_mgrm
```

### GRM from gene expression

```bash
# expr.fmat: genes x samples; expr.fam: samples; expr.bim: one line per gene
fastgma --make-grm --ffile expr --missing-rate 0.2 --out output/expr
```

### PCA, related-sample pruning and inverse

```bash
fastgma --process-grm --pca --grm output/test --out output/test             # top 20 PCs
fastgma --process-grm --pca --num-pc 0 --grm output/test --out output/test  # all PCs
fastgma --process-grm --out-close --grm output/test --cut-value 0.125 --out output/test
fastgma --process-grm --make-inv --grm output/test --out-fmt 3 --out output/test.inv
```

If the genotypes contain missing calls, `--make-inv` may stop with "not positive definite". The message gives the smallest eigenvalue and the `--val` needed; rerun with a larger value, for example `--val 0.02`.

### Use the GRM in GCTA

```bash
gcta64 --grm output/test --pheno test.phen --reml --out output/test_gcta
gcta64 --grm output/test --make-bK-sparse 0.05 --out output/test_sp
```


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/04_UVLMM/
================================================================
---
layout: page
permalink: /documentation/04_UVLMM/
title: uvLMM
description: Univariate linear mixed model for single trait
order: 4
---

`UVLMM` (univariate linear mixed model) is the single-trait GWAS engine of `fastGMA`. For one quantitative phenotype it first fits a **null model** to estimate the variance components (and heritability), then performs a **genome-wide single-SNP association scan** while controlling for the polygenic background and any additional random effects.

It is invoked through the top-level dispatcher with `--uvlmm` (short form `-u`):

```bash
fastgma --uvlmm [options]
```

Run `fastgma -h 3` for the full option list.

## 1. What it does

For each SNP, `UVLMM` tests the additive effect under a mixed model

$$
y = Xb + s\,g + u + e
$$

where $$Xb$$ are the fixed effects (intercept, covariates, class variables), $$s\,g$$ is the SNP being tested ($$g$$ its genotypes, $$s$$ its effect), $$u$$ is the polygenic random effect captured by one or more GRMs, and $$e$$ is the residual. The null model (no SNP) is fitted once to estimate the variance components. The genome-wide scan then uses a fast approximate test for every SNP and an exact test only for promising SNPs.

**Approximate test.** The exact test needs $$g'Pg$$ for every SNP. `UVLMM` approximates it by $$g'Pg \approx \gamma\, g'g$$. $$\gamma$$ is the mean of $$g'Pg / g'g$$ over `--num-random-snp` random SNPs that pass the MAF and missing-rate filters. The random SNPs are drawn from the whole genotype file with `--seed` (default `42`), so $$\gamma$$ is the same in reruns and in every `--parallel` / `--snp-set` job. SNPs with approximate p below `--p-approx-cut` are re-tested exactly. `--p-approx-cut 1` gives the exact test for every SNP; `--p-approx-cut 0` keeps the approximate test only.

## 2. Model and solver

The model is chosen with `--model` (default `lmm`). The fitting algorithm is chosen with `--solver`:

| `--model` | Model | `--solver` | Requires |
| --- | --- | --- | --- |
| `lmm` (default) | One GRM as random effect | `eigen` (default), `mme` or `general` | `--grm` |
| `lm` | Linear model, no random effect | not used | — |
| `multi-grm` | Several GRMs as random effects | `general` (default) or `mme` | `--mgrm` |
| `repeat` | GRMs plus a permanent (individual) effect for repeated records | `mme` (default) or `general` | `--mgrm` |

The three solvers:

| `--solver` | Method | Notes |
| --- | --- | --- |
| `eigen` | The GRM is eigendecomposed once and the data are rotated. REML and the per-SNP tests then cost $$O(n)$$ per evaluation. | `--model lmm` only. Needs one record per individual. `--random` is not allowed. Supports `--refit` / `--no-refit`. |
| `mme` | AI-REML on Henderson's mixed model equations (WEMAI as for `general`). The equations have one row per fixed effect and per random-effect level, so their size does not grow with the number of records; each GRM is inverted once. | Allows repeated records and `--random`. Needs positive-definite GRMs (see `--val` below). Same estimates and tests as `general`. Exact re-tests use the null variances. |
| `general` | Dense AI-REML on the $$n \times n$$ covariance matrix (WEMAI: the AI matrix is blended with the EM matrix when needed to keep variances positive). | Allows repeated records and `--random`. Exact re-tests use the null variances. |

**`mme` or `general`.** Both fit the same model; they differ in cost. Each `general` iteration factorizes the $$n \times n$$ covariance matrix of the $$n$$ records, while `mme` works with a matrix of size $$p + \sum_k q_k$$ ($$p$$ fixed effects, $$q_k$$ levels of random effect $$k$$). With many records per individual, as for repeated measures, `mme` is much faster and needs far less memory: on `example/mouse_long` (19,392 records of 1,212 mice, `--model repeat`, 8 threads) the null model and the scan of 11,833 SNPs take about 4 s and 0.2 GB in total. Each iteration costs about the cube of the matrix size in both, so `mme` is faster when $$p + \sum_k q_k$$ is smaller than the number of records. Each GRM adds a block for all individuals, so with three GRMs and two records per individual `mme` is slower than `general`. With one record per individual and one GRM, `eigen` is the fastest. When `--solver general` is used and the mixed model equations are smaller than the number of records, the log suggests `--solver mme`.

**GRM diagonal (`--val`).** `mme` needs the inverse of each GRM, and a centred GRM is singular (its rows sum to zero). `mme` therefore adds `--val` (default `0.001`) to the GRM diagonal when it reads the GRM, as `--longwas` does. The other solvers add nothing unless `--val` is given. To compare `mme` with `general` on the same model, give both the same `--val`. If a GRM is still not positive definite (for example, duplicated individuals), the run stops and suggests a larger `--val`.

Rules checked at start-up:

- `--solver eigen` or `--solver mme` with `--model lm` is an error. `--solver general` with `--model lm` is ignored with a warning.
- `--solver eigen` with `--model multi-grm` or `--model repeat` is an error.
- `--random` with `--solver eigen` is an error. Use `--model lmm --solver mme` (or `general`) to add non-genetic random effects.
- An individual with more than one record is an error with `--solver eigen`. Use `--solver mme` (or `general`), or `--model repeat` for a permanent effect.
- `--model repeat` cannot be combined with `--random`.
- `--refit` is accepted only with `--solver eigen`.

**Per-SNP refit (`--solver eigen`).** By default (`--refit`), each SNP that is re-tested exactly gets its own REML fit of the variance components with the SNP in the model. Score, likelihood-ratio and Wald p-values are then reported. With `--no-refit` the exact test uses the null variances, so the Wald p equals the score p. `--refit --p-approx-cut 1` refits every SNP.

## 3. Inputs

### Phenotype / covariate table (`--data`)

A whitespace- or tab-separated table with a header row. The first column must be the individual ID. Columns are selected by name. Only individuals present in both the data table and the GRM are retained.

| Option | Meaning |
| --- | --- |
| `--data` | Data file (**required**) |
| `--trait` | Column name of the trait — exactly one (**required**; UVLMM is univariate) |
| `--covar` | Column name(s) of continuous covariates |
| `--class` | Column name(s) of categorical covariates |
| `--random` | Column name(s) of additional non-genetic random effects (`--solver mme` or `general`) |
| `--missing-data` | Tokens treated as missing in the data table (default `NA NaN nan N/A ...`) |

### GRM input

| Option | Meaning |
| --- | --- |
| `--grm` | Prefix of a single GRM (`--model lmm`) |
| `--mgrm` | Prefixes of one or more GRMs (`--model multi-grm` / `--model repeat`) |
| `--val` | Value added to the GRM diagonal when it is read (default `0.001` for `--solver mme`, `0` otherwise; see Section 2) |

All `--mgrm` GRMs must list the same individuals in the same order. As elsewhere in `fastGMA`, GRMs are referenced by prefix (e.g. `--grm output/test`), not by a single file name. GRMs built by `fastgma --make-grm` and by GCTA are both accepted; see [GRM](../03_GRM/) for how to build them.

### Genotype input

| Option | Meaning |
| --- | --- |
| `--bfile` | PLINK binary genotype prefix (`.bed/.bim/.fam`) |
| `--ffile` | Feature-matrix prefix (`.fam/.bim/.fmat`). Each `.fmat` row is one feature and each column one sample — suitable for imputed dosages, expression, or any continuous per-sample matrix. |
| `--missing-geno` | Tokens treated as missing in the `.fmat` file (default `NA NaN nan N/A ...`) |

Give only one of `--bfile` / `--ffile`. A `--ffile` feature is centred by its mean, missing values get the mean, `--maf` is not applied and `af` is the feature mean. If no genotype file is supplied, the command fits the null model, writes `.var`, `.h2` and `.sln`, and stops before the SNP scan.

## 4. Outputs

All files share the `--out` prefix and are tab-separated with a header row:

| File | Contents |
| --- | --- |
| `<out>.var` | Estimated variance components with standard errors (inverse AI matrix at the estimates) |
| `<out>.h2` | Per-component proportion $$h^2_i = \sigma^2_i d_i / \sum_k \sigma^2_k d_k$$ ($$d_i$$ = mean GRM diagonal, 1 for non-GRM effects), SE by the delta method (not written for `--model lm`) |
| `<out>.sln` | Fixed-effect estimates; with `--predict`, also the predicted random effects (BLUP) and their reliabilities |
| `<out>.res` | With `--predict`: each record split into its fixed, random and residual parts |
| `<out>.assoc` | One row per SNP with the association results |
| `<out>.log` | Log of the run (appended; `<out>.<m>_<i>.log` with `--parallel`, `<out>.<start>_<end>.log` with `--snp-set`) |

`<out>.var` and `<out>.h2` are keyed by component name rather than a positional
index: `G` for the `--grm` effect, the GRM prefix basename for each
`--mgrm` entry, `permanent` for the `--model repeat` permanent-environment
effect, the column name for each `--random` effect, and `residual` last.
`--model lm` writes only the `residual` row.

```
component	variance	se
G	1.26776	0.206283
residual	2.39709	0.105419
```

`<out>.h2` adds a trailing `total` row holding the combined proportion of the
GRM components — the number usually quoted as the heritability. A
`permanent` or `--random` effect contributes variance but is not counted here:

```
component	h2	se
G	0.349279	0.0403661
residual	0.650721	0.0403661
total	0.349279	0.0403661
```

Each GRM variance is put on the phenotypic scale by the mean diagonal $$d_k$$ of its
kernel $$Z G_k Z'$$ over the analysed records (for one record per individual, the mean
diagonal of the GRM over the analysed individuals); a `permanent` / `--random` effect
and the residual have $$d = 1$$:

$$
h^2_k = \frac{\sigma^2_k d_k}{V_p}, \qquad V_p = \sum_k \sigma^2_k d_k \quad \bigl(\text{the mean of } \operatorname{diag}(V)\bigr)
$$

This is the definition used by `--fast-reml` (`variance_prop`) and the PVE of GEMMA,
and the same for every model and solver, so `--grm A` and `--mgrm A` give the same
`.h2`. It matters for GRMs whose mean diagonal is far from 1 (for example an $$A \circ A$$
epistatic GRM). For a standard additive GRM $$d$$ is $$1 + F$$ averaged over the sample
(close to 1); GCTA `--reml` reports the unscaled `V(G)/Vp`, which can be computed from
`<out>.var` as $$\sigma^2_G / (\sigma^2_G + \sigma^2_e)$$. The $$d_k$$ are written to the log
("mean diagonal over the records"). SEs use the delta method on the inverse AI matrix.

Values that are not finite (for example an SE that could not be derived from the
AI matrix) print as `NA`.

`<out>.sln` has the columns `type term levels effect se`, and with `--predict` also `rel`. `type` is `0` for the intercept (`mu`), `1` for a covariate and `2` for a class level. A fixed-effect column removed as linearly dependent (for example the reference level of a class variable) is listed with effect `0` and se `NA`. With `--predict`, type `3` rows hold the predicted genetic values (BLUP) and type `4` rows the predicted `--random` / `permanent` effects. For these rows `term` is the variance component as named in `<out>.var` (`G`, the GRM prefix basename, `permanent` or the `--random` column), `levels` is the individual ID or level, `se` is the square root of the prediction error variance (PEV), and `rel` is the reliability $$1 - \mathrm{PEV} / (\sigma^2_k K_{jj})$$ ($$K_{jj}$$ = GRM diagonal, `1` for a non-genetic effect); `rel` is `NA` on the fixed-effect rows. GRM rows cover every GRM individual, including those without records, whose values are predicted through their relationships. Numbers keep 10 significant digits. From the [repeated-records](#repeated-records) example with `--predict` (first rows of each kind):

```
type	term	levels	effect	se	rel
0	mu	1	19.7737604	0.07620444974	NA
2	Sex	1	-2.65796099	0.09581998468	NA
2	Sex	2	0	NA	NA
2	week	1	-13.23082676	0.05535989006	NA
...
3	mouse_long	1419	0.4080401921	0.3347456694	0.71612076
...
4	permanent	1000	-0.7212183419	0.4809213132	0.9077624809
...
```

`<out>.res` (with `--predict`) has one row per record: `id observed fixed random residual`, where `fixed` is $$x_i'\hat b$$, `random` is the sum of the predicted random effects of the record and `residual` is $$\hat e_i$$, so that `observed` = `fixed` + `random` + `residual`. `observed` - `fixed` is the phenotype corrected for the fixed effects. From the same run:

```
id	observed	fixed	random	residual
1419	4.9	3.884972647	0.4058137606	0.6092135924
1419	8.4	7.348924696	0.4058137606	0.6452615437
...
```

The `<out>.assoc` columns depend on the solver.

`--solver eigen` (`--model lmm` default):

```
chr snp pos a1 a2 af N cc_par cc_gra beta se p_score p_lrt p_wald
```

`cc_par` and `cc_gra` are the convergence diagnostics of the per-SNP refit (norm of the last parameter change and of the gradient). They are written only with `--refit` (the default) and are `nan` for SNPs that were not re-tested exactly. With `--no-refit` the header is `chr snp pos a1 a2 af N beta se p_score p_lrt p_wald`. For a SNP that was not re-tested exactly, `beta` and `se` are the approximate values and the three p-value columns all hold the approximate p.

`--model lm` (exact least-squares test for every SNP):

```
chr snp pos a1 a2 af N beta se p_score p_lrt p_wald
```

`--solver general` or `mme` (`multi-grm`, `repeat`, or `--model lmm` with either solver): only the Wald test is computed, so a single `p` column is written:

```
chr snp pos a1 a2 af N beta se p
```

Columns: chromosome `chr`, SNP ID `snp`, position `pos`, effect allele `a1` and other allele `a2`, allele frequency of `a1` (or feature mean) `af`, `N` (number of individuals with a non-missing genotype), the effect estimate `beta` with its standard error `se`, and the p-values. SNPs that fail the MAF/missing-rate filters report `nan`. With `--p-cut` below `1`, only SNPs whose final p (`p_wald`, or `p` for the general and mme solvers) is at or below the cut are written.

## 5. Key options

### Null model

| Option | Default | Meaning |
| --- | --- | --- |
| `--model` | `lmm` | `lmm`, `lm`, `multi-grm` or `repeat` |
| `--solver` | `eigen` for `lmm`; `general` for `multi-grm`; `mme` for `repeat` | `eigen`, `mme` or `general` (see Section 2) |
| `--val` | `0.001` for `mme`, `0` otherwise | Value added to the GRM diagonal (see Section 2) |
| `--maxiter0` | `100` | Max iterations for the null-model fit |
| `--continue` | off | Start from the variances in an existing `<out>.var` |
| `--predict` | off | Also write the predicted random effects and their reliabilities to `<out>.sln`, and `<out>.res` |
| `--threads` | `10` | Number of MKL/OpenMP threads |

### SNP scan

| Option | Default | Meaning |
| --- | --- | --- |
| `--maf` | `0.01` | Minimum minor allele frequency $$\min(p, 1 - p)$$ (PLINK input only) |
| `--missing-rate` | `0.05` | Keep a SNP if its missing rate is at or below this |
| `--p-approx-cut` | `1e-3` | SNPs with approximate p below this are re-tested exactly. `0` = approximate test only; `1` = exact test for every SNP |
| `--refit` / `--no-refit` | `--refit` | `--solver eigen`: re-estimate the variance components for every exactly re-tested SNP |
| `--maxiter` | `30` | Max iterations for the per-SNP refit |
| `--p-cut` | `1` | Write only SNPs whose final p-value is at or below this (`1` = all, including `nan` rows) |

### Restricting / parallelizing the scan

| Option | Meaning |
| --- | --- |
| `--snp-set <first> <last>` | Test only the SNPs from the first to the last SNP ID (inclusive). The output prefix becomes `<out>.<start>_<end>` (0-based start, exclusive end), for example `<out>.<start>_<end>.assoc` and the log `<out>.<start>_<end>.log`. |
| `--parallel <m> <i>` | Split the SNPs into `m` equal parts and test part `i` (1-based). The output prefix becomes `<out>.<m>_<i>`, for example `<out>.<m>_<i>.assoc` and the log `<out>.<m>_<i>.log`. |

`--snp-set` and `--parallel` cannot be used together; `--start-pos` / `--end-pos` apply only without them. Each job fits the same null model and writes all its files (`.var`, `.h2`, `.sln`, `.assoc`, `.log`) to its own prefix. Because `--seed` is fixed, all jobs use the same $$\gamma$$, and the merged `.assoc` files match a single full run.

### Advanced (rarely needed)

Grouped under `Advanced` in `--help`; defaults are almost always fine.

| Option | Default | Meaning |
| --- | --- | --- |
| `--cc-par0` / `--cc-gra0` | `1e-7` / `1e-6` | Null-model parameter / gradient convergence criteria |
| `--cc-logl0` | `5e-5` | Null-model convergence criterion for the change in log-likelihood |
| `--cc-par` / `--cc-gra` | `1e-7` / `1e-6` | Per-SNP refit parameter / gradient convergence criteria |
| `--start-pos` / `--end-pos` | `0` / `-1` | 0-based index of the first SNP and one past the last SNP to test (`-1` = up to the last SNP); output prefix unchanged |
| `--npart-snp` | `20` | Minimum number of chunks the tested SNPs are read in; chunks are kept to about 128 MB of genotype codes (memory only; does not change results) |
| `--num-random-snp` | `2000` | Random SNPs for the approximate-test calibration (capped at the number of SNPs) |
| `--seed` | `42` | Random seed for drawing the calibration SNPs |

Set the environment variable `FASTGMA_LOG_LEVEL=debug` to log per-block progress.

## 6. Examples

### Standard single-GRM LMM scan

```bash
fastgma --uvlmm --model lmm \
  --grm output/test \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --covar age \
  --class sex \
  --out output/test_uvlmm
```

`--model lmm` and `--solver eigen` are the defaults, so they can be omitted.

### Exact test with a per-SNP refit for every SNP

```bash
fastgma --uvlmm \
  --grm output/test \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --p-approx-cut 1 \
  --out output/test_uvlmm_exact
```

### Variance components only (no SNP scan)

```bash
fastgma --uvlmm --model lmm \
  --grm output/test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --out output/test_uvlmm
```

### Single GRM with a non-genetic random effect

```bash
fastgma --uvlmm --model lmm --solver general \
  --grm output/test \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --random litter \
  --out output/test_uvlmm_litter
```

### Repeated records

```bash
fastgma --make-grm --bfile example/mouse_long/plink --out output/mouse_long
fastgma --uvlmm --model repeat \
  --mgrm output/mouse_long \
  --bfile example/mouse_long/plink \
  --data example/mouse_long/phe.balance.txt \
  --trait trait \
  --class Sex week \
  --out output/mouse_long_repeat
```

16 weekly body weights of each of 1,212 mice: the GRM, a permanent effect per mouse and the residual. `--solver mme` is the default for `--model repeat`.

### Multiple GRMs

```bash
fastgma --uvlmm --model multi-grm \
  --mgrm output/grm_add output/grm_dom \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --out output/test_uvlmm_mgrm
```

### Split a scan into 10 parallel jobs

```bash
fastgma --uvlmm --model lmm \
  --grm output/test \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait pheno \
  --parallel 10 1 \
  --out output/test_uvlmm
```

Repeat with part indices `2` through `10`, then concatenate `output/test_uvlmm.10_*.assoc` (keeping one header).


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/05_REMMA/
================================================================
---
layout: page
permalink: /documentation/05_REMMA/
title: REMMA
description: Non-additive effects analysis with rapid linear mixed model
order: 5
---

`REMMA` (Random Effect Model for Multiple Associations) tests SNPs for additive, dominance and additive-by-additive epistatic effects. It first fits a **null model** with one or more GRMs as random effects (for example additive, dominance and epistatic GRMs). It then runs one test, per SNP or per SNP pair, with the variance components of that null model.

It is invoked through the top-level dispatcher with `--remma` (there is no short form):

```bash
fastgma --remma <test> [options]
```

For the full option list, run:

```bash
fastgma -h 4
```

## 1. What it does

The null model is

$$
\begin{aligned}
y &= Xb + Z u_1 + \cdots + Z u_K + e \\
u_k &\sim N(0, \sigma^2_k G_k), \qquad e \sim N(0, \sigma^2_e I) \\
V &= \textstyle\sum_k \sigma^2_k Z G_k Z' + \sigma^2_e I
\end{aligned}
$$

$$Xb$$ holds the fixed effects (intercept, covariates, class variables). Each $$u_k$$ is a random effect with covariance matrix $$G_k$$, one per GRM given with `--mgrm`. $$Z$$ maps records to individuals, so repeated records per individual are allowed and are used by every test. The variance components are estimated once by AI-REML (WEMAI: the AI matrix is blended with the EM matrix when needed to keep the variances positive). The start values are scale-aware: the phenotypic variance is split evenly over the $$K$$ GRMs and the residual. 

`--solver general` (default) fits on the dense $$n \times n$$ $$V$$ of the $$n$$ records; `--solver mme` fits the same model on the mixed model equations, whose size $$p + Kq$$ ($$p$$ fixed effects, $$q$$ individuals) does not grow with the number of records. `mme` is faster when $$p + Kq$$ is smaller than the number of records, for example with more than about three records per individual for three GRMs (see [uvLMM](../04_UVLMM/), Section 2). Both give the same estimates and tests. `mme` adds `--val` (default `0.001`) to each GRM diagonal so that the GRM can be inverted; give `general` the same `--val` to fit the identical model.

All tests below reuse $$V$$ from this fit, and the variances are not re-estimated for each SNP.

Let $$P = V^{-1} - V^{-1} X (X' V^{-1} X)^{-1} X' V^{-1}$$, $$S = Z'PZ$$ and $$s = Z'Py$$. Genotypes are coded per individual as the number of copies of the A1 allele (`0/1/2`). Missing genotypes are replaced by the SNP mean $$2p$$ for the additive code. With `--ffile` the code is the feature value, centred by its mean (missing values get the mean).

## 2. Tests

Choose exactly one test. The run stops with an error if none or more than one is given.

| Flag | Test | Model per SNP / pair |
| --- | --- | --- |
| `--add` | Additive | Exact GLS statistic with the null-model `V` |
| `--dom` | Dominance | Joint GLS fit: fixed effects + additive + dominance term |
| `--epis` | Additive-by-additive epistasis | Approximate scan of the interaction alone for every SNP pair; pairs below `--p-approx-cut` are refitted with both main effects |

### Additive test (`--add`)

Let $$x$$ be the centred SNP code (the dosage minus $$2p$$). With $$V$$ fixed, the GLS estimate of the SNP effect and its standard error are

$$
\hat\beta = \frac{x's}{x'Sx}, \qquad \operatorname{se}(\hat\beta) = \sqrt{\frac{1}{x'Sx}}, \qquad \frac{\hat\beta^2}{\operatorname{se}(\hat\beta)^2} \sim \chi^2_1
$$

This is the exact GLS statistic for $$y = Xb + Zxa + \cdots + e$$ at the null $$V$$. It needs only $$s$$ and $$S$$, so it is very fast. It is the same test as [uvLMM](../04_UVLMM/) `--model multi-grm --p-approx-cut 1` with the same GRMs, which gives the same `beta` and `se`; `REMMA` computes p from the chi-square distribution and uvLMM from the F distribution whose denominator degrees of freedom are the number of records minus the number of fixed effects minus 1, so p differs slightly in small samples. `--add` gives the additive results with the same null model as `--dom` and `--epis`.

### Dominance test (`--dom`)

For each SNP, `REMMA` fits

$$
y = Xb + Z x_a a + Z x_d d + Z u_1 + \cdots + Z u_K + e
$$

by generalized least squares with the fixed $$V$$ from the null model. $$x_a$$ is the centred additive code $$x - 2p$$. $$x_d$$ is the statistical dominance code $$0 \to -2p^2$$, $$1 \to 2pq$$, $$2 \to -2q^2$$, the same as in the dominance GRM (`--grm-type dom`) (see [GRM](../03_GRM/)). Here $$p$$ is the A1 frequency and $$q = 1 - p$$. A missing genotype gets the mean dominance code of the observed genotypes, so it carries no information (like the mean-imputed additive code). Wald tests (chi-square, 1 df) are reported for both $$a$$ and $$d$$ from the same fit.

With this coding, $$a$$ is the allele substitution effect, the quantity estimated by an ordinary additive GWAS, and $$d$$ is the dominance deviation. The genotypic additive effect (half the difference between the two homozygotes) is $$a - d(q - p)$$.

**Non-identifiable SNPs.** If one of the three genotype classes is absent among the analysed individuals, or the two codes are collinear after adjusting for the fixed effects, the dominance effect cannot be estimated. `REMMA` then reports `beta_dom`, `se_dom` and `p_dom` as `nan`, and the `_add` columns hold the additive-only test (the same statistic as `--add`). The log counts these SNPs.

**Small genotype classes (`--min-geno-count N`).** A dominance effect estimated from a very small genotype class (for example a few minor-allele homozygotes) can give an extreme `p_dom`. With `--min-geno-count N` (`N > 0`), the joint fit is only done when each genotype class `0`/`1`/`2` has at least `N` analysed individuals with a non-missing genotype; otherwise the SNP is treated exactly like a non-identifiable SNP (`_dom` columns `nan`, `_add` columns from the additive-only test). The default `0` turns the rule off (only an absent class triggers the fallback). The log reports the setting and how many tested SNPs it affected. The rule applies to `--dom` only; `--add` and `--epis` ignore the option with a warning.

### Epistasis test (`--epis`)

For a SNP pair $$i, j$$ the interaction code is $$z = x_i \circ x_j$$, the elementwise product of the two centred additive codes, centred again by its mean. The scan has two stages.

**Stage 1: approximate scan without main effects.** Every pair is tested with the interaction alone, $$y = Xb + Z z c + Z u_1 + \cdots + Z u_K + e$$:

$$
\hat c = \frac{z's}{z'Sz}, \qquad \operatorname{se}(\hat c) = \sqrt{\frac{1}{z'Sz}}
$$

Computing $$z'Sz$$ for every pair is slow, so `REMMA` approximates it by $$z'Sz \approx \gamma\, z'z$$:

1. It draws `--num-random-snp` random SNP pairs from all SNPs in the genotype file, with `--seed` (default `42`), and computes $$z'Sz / z'z$$ for each.
2. $$\gamma$$ is the mean of this ratio over the pairs that pass the MAF and missing-rate filters and give a finite, positive ratio. The ratio does not depend on the phenotype, so no p-value filter is applied (as in `UVLMM` and `MVLMM`). At least 10% of the random pairs must be usable, or the run stops.

With $$\gamma$$, $$z's$$ and $$z'z$$ of a block of pairs are matrix products of the SNP codes, so stage 1 costs about as much as reading the genotypes.

**Stage 2: full model below `--p-approx-cut`.** Pairs with approximate p below `--p-approx-cut` (default `1e-6`) are refitted by GLS with both main effects,

$$
y = Xb + Z x_i a_i + Z x_j a_j + Z z c + Z u_1 + \cdots + Z u_K + e
$$

and the Wald test of $$c$$ is reported. If the two SNPs are in complete LD, so their main effects are collinear, the redundant main effect $$x_j$$ is dropped and $$[X, x_i, z]$$ is fitted. For SNPs that are not in LD, $$z$$ is nearly uncorrelated with $$x_i$$ and $$x_j$$, so the two stages give similar p-values and stage 1 is a good filter.

With human genotypes (`example/human`, 427 individuals, three GRMs, the 5,137,150 pairs whose first SNP is one of the first 20), every pair with full-model p below a tenth of the cut reached stage 2; pairs just below the cut can be missed when the two SNPs are in LD. With the default `1e-6`, 25 pairs were refitted. All 4 pairs with full-model p < 1e-7 were among them, and 19 of the 25 pairs with full-model p < 1e-6.

Genome-wide significance for a scan of all pairs is far smaller (about 4e-13 for 500,000 SNPs). The cost of stage 2 grows with the cut: it refits about a fraction of all pairs equal to the cut, at $$O(q^2)$$ each for $$q$$ individuals, so `1e-3` refits 1.25e8 pairs for 500,000 SNPs. Raise the cut (for example `1e-4`) only when pairs with moderate p-values are of interest.

`--p-approx-cut 1` fits the full model for every pair (no $$\gamma$$ is estimated; `--num-random-snp` and `--seed` are not used). This costs $$O(q^2)$$ per pair for $$q$$ individuals. `--p-approx-cut 0` keeps the stage 1 results only. `--marginal`, the old option for the full model, is accepted as `--p-approx-cut 1` with a warning.

Because the random pairs come from the whole genotype file and the seed is fixed, every job of a split scan uses the same $$\gamma$$, and reruns give identical results.

**Which pairs are tested.** A run (or a `--snp-set` / `--parallel` / `--start-pos`/`--end-pos` job) tests every pair $$(i, j)$$ with $$j > i$$ and SNP $$i$$ in the tested range; $$j$$ runs up to the last SNP of the genotype file. Jobs that cover the SNPs once therefore test every pair exactly once. `--parallel m i` splits the SNPs so that each part holds about the same number of pairs (early parts hold fewer SNPs). The tested range is held in memory (about $$16qm$$ bytes for $$m$$ SNPs and $$q$$ individuals), and the later SNPs are read in blocks of the same size, but at least 1,024 SNPs.

## 3. Inputs

### Phenotype / covariate table (`--data`)

A whitespace- or tab-separated table with a header row. The first column must be the individual ID; it is matched to the IID of the `.grm.id` files and of the `.fam` file. Columns are selected by name. Only individuals present in both the table and the GRMs are used, and they must all be in the `.fam` file. An individual may have several records.

| Option | Meaning |
| --- | --- |
| `--data` | Phenotype/data file (**required**) |
| `--trait` | Column name of the trait — exactly one (**required**) |
| `--covar` | Column name(s) of continuous covariates |
| `--class` | Column name(s) of categorical covariates |
| `--missing-data` | Tokens treated as missing in the table (default `NA NaN nan N/A ...`) |

Records with a missing value in any used column are removed.

### GRM input

| Option | Meaning |
| --- | --- |
| `--mgrm` | One or more GRM prefixes (**required**). Each GRM is one random effect in the null model. |

All GRMs must list the same individuals in the same order. GRMs are referenced by prefix (`--mgrm output/add` reads `output/add.grm.bin` and `output/add.grm.id`). GRMs built by `fastgma --make-grm` and by GCTA are both accepted. A typical set is:

| Random effect | How to build it |
| --- | --- |
| Additive | `fastgma --make-grm --grm-type add` |
| Dominance | `fastgma --make-grm --grm-type dom` |
| Additive × additive | `fastgma --process-grm --make-epis --prefix1 <add> --prefix2 <add>` |

See [GRM](../03_GRM/) for details. A single additive GRM is also allowed.

### Genotype input

| Option | Meaning |
| --- | --- |
| `--bfile` | PLINK binary genotype prefix (`.bed/.bim/.fam`) |
| `--ffile` | Feature-matrix prefix (`.fam/.bim/.fmat`). Each `.fmat` row is one feature and each column one sample, for imputed dosages, expression or any other continuous per-sample matrix |
| `--missing-geno` | Tokens treated as missing in the `.fmat` file (default `NA NaN nan N/A ...`) |

Give exactly one of `--bfile` and `--ffile`. `--ffile` works with `--add` and `--epis`; `--dom` needs the genotype classes `0/1/2` and therefore `--bfile`. For a feature, `--maf` is not applied and `af` is the feature mean.

## 4. Outputs

All files share the `--out` prefix and are tab-separated with a header row. With `--snp-set` the prefix becomes `<out>.<start>_<end>` (0-based start, exclusive end); with `--parallel m i` it becomes `<out>.<m>_<i>`. The suffix applies to every file of that job.

| File | Written by | Contents |
| --- | --- | --- |
| `<out>.var` | all tests | Null-model variance components with standard errors (inverse AI matrix at the estimates) |
| `<out>.h2` | all tests | Per-component proportion of variance $$\sigma^2_k d_k / V_p$$ ($$d_k$$ = mean GRM diagonal), SE by the delta method |
| `<out>.sln` | all tests | Fixed-effect estimates of the null model |
| `<out>.assoc` | all tests | Association results |
| `<out>.random.assoc` | `--epis` with `--p-approx-cut` below 1 | Stage 1 results for the random pairs used to estimate $$\gamma$$ |
| `<out>.log` | all tests | Log of the run (appended; `<out>.<m>_<i>.log` with `--parallel`, `<out>.<start>_<end>.log` with `--snp-set`) |

### Null-model files

`<out>.var` and `<out>.h2` have one row per GRM, named by the basename of its `--mgrm` prefix, and `residual` last. `<out>.h2` adds a `total` row: the summed proportion of all GRM components.

The proportions put each GRM variance on the phenotypic scale by the mean diagonal $$d_k$$ of its kernel $$Z G_k Z'$$ over the analysed records (the mean of the diagonal elements $$G_k(j,j)$$ over the individuals with records, each counted once per record):

$$
h^2_k = \frac{\sigma^2_k d_k}{V_p}, \qquad V_p = \sum_k \sigma^2_k d_k + \sigma^2_e
$$

so GRMs of different scale (for example an $$A \circ A$$ epistatic GRM, whose mean diagonal can be far from 1) are comparable. This is the definition used by [fastREML](../09_FASTREML/) and [uvLMM](../04_UVLMM/). The $$d_k$$ are written to the log. SEs are from the delta method on the inverse AI matrix. `<out>.var` holds the unscaled variances, so $$\sigma^2_k / \sum \sigma^2$$ (the convention of GCTA `--reml`) can be computed from it.

The example output below is from the mouse data (`example/mouse/geno/plink`, `example/mouse/pheno/pheno_simu.txt`, `--trait trait --covar e1 --class c1`) with the GRMs `add` (`--grm-type add`), `dom` (`--grm-type dom`) and `axa` ($$\text{add} \circ \text{add}$$):

```
component	variance	se
add	1.15076	0.214824
dom	0.00432849	0.0642098
axa	0.0225108	0.186052
residual	3.20265	0.197495
```

```
component	h2	se
add	0.265346	0.0396146
dom	0.000985636	0.014621
axa	0.00531368	0.0439219
residual	0.728355	0.0492958
total	0.271645	0.0492958
```

(Mean diagonals in this example: `add` 1.0139, `dom` 1.00126, `axa` 1.03794.)

`<out>.sln` has the columns `type term levels effect se`. `type` is `0` for the intercept (`mu`), `1` for a covariate and `2` for a class level. A class level removed because it is linearly dependent (the reference level) is shown with effect `0` and se `NA`. Numbers keep 10 significant digits.

```
type	term	levels	effect	se
0	mu	1	-0.3537821355	0.09129408205
1	e1	1	-0.6772010854	0.05123134883
2	c1	0	0.6430544761	0.1281130157
2	c1	1	1.346277836	0.1284049036
2	c1	2	0	NA
```

### `--add`

```
chr snp pos a1 a2 af N beta se p
```

### `--dom`

```
chr snp pos a1 a2 af N beta_add se_add p_add beta_dom se_dom p_dom
```

Columns: `chr`, `snp`, `pos`, `a1` and `a2` come from the `.bim` file; `af` is the A1 allele frequency in the analysed individuals (the feature mean with `--ffile`) and `N` the number of analysed individuals with a non-missing genotype. `beta`/`se`/`p` are the additive results. `beta_add`/`se_add`/`p_add` and `beta_dom`/`se_dom`/`p_dom` are the additive and dominance terms from the joint fit (see Section 2 for non-identifiable SNPs).

With the default `--p-cut 1`, every SNP in the tested range gets one row. SNPs that fail the MAF or missing-rate filter, or give a non-finite statistic, report `nan` for the effect, SE and p-value.

### `--epis`

```
chr_1 snp_1 pos_1 a1_1 a2_1 af_1 N_1 chr_2 snp_2 pos_2 a1_2 a2_2 af_2 N_2 beta se p
```

The first seven columns describe the first SNP of the pair, the next seven the second SNP, and `beta`/`se`/`p` the interaction effect. `<out>.random.assoc` has the same columns and lists every random pair, with `nan` for pairs that fail the filters.

Only pairs in which both SNPs pass the filters are written. Pairs with stage 1 p below `--p-approx-cut` carry `beta`, `se` and `p` of the interaction in the full model with both main effects; the others carry the stage 1 values. With the default `--p-cut 1`, every tested pair is written, which is about $$M^2/2$$ rows for $$M$$ SNPs. Use `--p-cut` to keep the file small.

## 5. Options

### Main options

| Option | Default | Meaning |
| --- | --- | --- |
| `--data` | — | Phenotype/covariate table (**required**) |
| `--mgrm` | — | GRM prefix(es) of the null model (**required**) |
| `--bfile` / `--ffile` | — | PLINK binary prefix, or feature-matrix prefix (`--add`, `--epis`); exactly one is **required** |
| `--missing-geno` | `NA NaN nan N/A ...` | Missing-value tokens in the `.fmat` file |
| `--out` | — | Output prefix (**required**) |
| `--trait` | — | Trait column, exactly one (**required**) |
| `--covar`, `--class` | none | Covariate and class columns |
| `--missing-data` | `NA NaN nan N/A ...` | Missing-value tokens in the data file |
| `--add` / `--dom` / `--epis` | — | The test (exactly one) |
| `--maf` | `0.01` | Keep SNPs with $$\min(p, 1 - p)$$ at or above this, in `[0, 0.5]` (`--bfile` only) |
| `--missing-rate` | `0.05` | Keep SNPs with missing rate at or below this, in `[0, 1]` |
| `--min-geno-count` | `0` | `--dom`: minimum number of individuals in each genotype class for the dominance test; otherwise `_dom` is `nan` and the additive-only test is reported (`0` = off) |
| `--snp-set <first> <last>` | — | Test SNPs from the first to the last SNP ID (inclusive); with `--epis`, pairs whose first SNP is in this set |
| `--parallel <m> <i>` | — | Split the scan into `m` parts and run part `i` (1-based); with `--epis`, parts hold about equal numbers of pairs |
| `--p-cut` | `1` | Write SNPs/pairs whose final p is at or below this. For dominance a row is kept if `p_add` or `p_dom` passes. `1` writes all |
| `--p-approx-cut` | `1e-6` | `--epis`: refit pairs with both main effects when the stage 1 p is below this; `1` = full model for every pair, `0` = stage 1 only |
| `--solver` | `general` | Null-model algorithm: `general` (dense $$V$$) or `mme` (mixed model equations; faster when $$p + Kq$$ is below the number of records) |
| `--val` | `0.001` for `mme`, `0` for `general` | Value added to the GRM diagonal when it is read |
| `--threads` | `10` | Number of threads |

`--snp-set`, `--parallel` and `--start-pos`/`--end-pos` are mutually exclusive.

### Advanced

| Option | Default | Meaning |
| --- | --- | --- |
| `--npart-snp` | `20` | Minimum number of chunks the tested SNPs are read in; chunks are kept to about 128 MB of genotype codes (memory only; `--add` and `--dom`, not `--epis`) |
| `--start-pos` | `0` | First tested SNP (0-based index in the `.bim` file) |
| `--end-pos` | `-1` | End of the tested SNPs (0-based, exclusive). `-1`, or a value above the number of SNPs, means the last SNP |
| `--num-random-snp` | `10000` | `--epis`: random SNP pairs for the $$\gamma$$ calibration; capped at the number of SNPs |
| `--seed` | `42` | `--epis`: random seed for the calibration pairs |
| `--maxiter0` | `100` | Max iterations for the null-model fit |
| `--cc-par0` | `1e-7` | Null-model convergence: norm of the parameter change |
| `--cc-gra0` | `1e-6` | Null-model convergence: norm of the gradient |
| `--cc-logl0` | `5e-5` | Null-model convergence: change in REML log-likelihood |

Set the environment variable `FASTGMA_LOG_LEVEL=debug` to log per-chunk progress.

## 6. Examples

### Build the GRMs

```bash
fastgma --make-grm --grm-type add --bfile test --out output/add
fastgma --make-grm --grm-type dom --bfile test --out output/dom
fastgma --process-grm --make-epis --prefix1 output/add --prefix2 output/add --out output/axa
```

### Additive test with three GRMs

```bash
fastgma --remma --add \
  --bfile test \
  --mgrm output/add output/dom output/axa \
  --data test_simu_pheno.txt \
  --trait pheno \
  --covar age \
  --class sex \
  --threads 4 \
  --out output/test_remma_add
```

### Dominance test

```bash
fastgma --remma --dom \
  --bfile test \
  --mgrm output/add output/dom output/axa \
  --data test_simu_pheno.txt \
  --trait pheno \
  --covar age \
  --class sex \
  --out output/test_remma_dom
```

Add `--start-pos 0 --end-pos 200` to test only the first 200 SNPs.

### Epistasis scan, split into parallel jobs

```bash
for i in 1 2 3 4; do
  fastgma --remma --epis \
    --bfile test \
    --mgrm output/add output/dom output/axa \
    --data test_simu_pheno.txt \
    --trait pheno \
    --num-random-snp 1000 \
    --p-cut 1e-3 \
    --parallel 4 $i \
    --out output/test_remma_epis
done
```

Each job writes `output/test_remma_epis.4_<i>.*`. Each job tests the pairs whose first SNP is in its part, and the parts hold about equal numbers of pairs. Together they cover every pair once. Concatenate the `.assoc` files (keeping one header) to get the genome-wide result.

### Epistasis with the full model for every pair

```bash
fastgma --remma --epis --p-approx-cut 1 \
  --bfile test \
  --mgrm output/add output/dom output/axa \
  --data test_simu_pheno.txt \
  --trait pheno \
  --start-pos 0 --end-pos 10 \
  --out output/test_remma_epis_full
```

This fits the full model for every pair, so keep the range small.


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/06_GLMM/
================================================================
---
layout: page
permalink: /documentation/06_GLMM/
title: GLMM
description: Generalized linear mixed model for binary traits
order: 6
---

`GLMM` (generalized linear mixed model) is the GWAS engine of `fastGMA` for **binary (case/control) traits**. It fits a logistic **null model** with one or more GRMs as polygenic random effects, then tests every SNP for an **additive effect** with a score test. SNPs with a small p-value get a **saddlepoint approximation (SPA)**, which keeps the test calibrated when cases and controls are unbalanced.

It is invoked through the top-level dispatcher with `--glmm`:

```bash
fastgma --glmm [options]
```

Run `fastgma -h 5` to print the full option list.

## 1. What it does

### Null model

The trait $$y$$ must be coded `0`/`1`, and both cases and controls must be present. The null model is the logistic mixed model

$$
\operatorname{logit}(\mu) = \eta = Xb + Z g_1 + \cdots + Z g_K, \qquad g_k \sim N(0, \sigma^2_k G_k)
$$

where $$\mu = P(y = 1)$$, $$Xb$$ are the fixed effects (intercept, covariates, class variables), $$G_k$$ are the GRMs and $$Z$$ maps records to individuals. An individual may have several records.

> **Repeated records.** The model has no permanent environmental effect: the records of one individual are treated as conditionally independent given its genetic effects. Non-genetic causes shared by the records of an individual are absorbed into $$g$$ or ignored, which can overstate the evidence. The log warns when the data contain repeated records.

The fit starts, as in SAIGE and GMMAT, from the logistic regression without random effects (IRLS), with every $$\sigma^2_k = 1$$. It then uses penalized quasi-likelihood (PQL). Each outer iteration forms the working response and weights

$$
\begin{aligned}
y^* &= \eta + \frac{y - \mu}{\mu(1 - \mu)} \\
W &= \operatorname{diag}\{\mu(1 - \mu)\} \\
V &= W^{-1} + \textstyle\sum_k \sigma^2_k Z G_k Z'
\end{aligned}
$$

and estimates the $$\sigma^2_k$$ from the linear mixed model on $$y^*$$ with AI-REML. As in SAIGE and GMMAT, a variance that an AI step takes below `1e-5` is set to `0` (a boundary estimate, logged); a later step can take it back above `0`, which is also logged. Only if the AI step is not finite is the AI matrix blended with the EM matrix (WEMAI) to keep the variances positive. The fixed effects $$b$$ and the random effects $$g_k$$ are then updated for the next iteration. The outer loop stops when the changes in the variances and in $$b$$ are both below `--cc-par0`. There is no residual variance: the working residual covariance is $$W^{-1}$$.

The null model stores and inverts dense $$n \times n$$ matrices ($$n$$ = number of records), so memory grows with $$n^2$$ and time with $$n^3$$.

### SNP test

The test uses $$y^*$$, $$\mu$$ and $$V$$ from the final linearization of the null model. Let $$P = V^{-1} - V^{-1}X(X'V^{-1}X)^{-1}X'V^{-1}$$, $$S = Z'PZ$$ and $$s = Z'Py^*$$, and let $$x$$ be the centred SNP code of the individuals (the A1 count minus $$2p$$, with missing genotypes set to $$2p$$). The score of the SNP is $$T = x's$$ with variance $$x'Sx$$; at convergence $$Py^* = y - \mu$$, so $$T = x'Z'(y - \mu)$$ is the score test of GMMAT (`glmm.score`). `GLMM` reports

$$
\hat\beta = \frac{x's}{x'Sx}, \qquad \operatorname{se}(\hat\beta) = \frac{1}{\sqrt{x'Sx}}, \qquad p_\text{normal} = P\left(\chi^2_1 > \frac{T^2}{x'Sx}\right)
$$

$$\hat\beta$$ is the one-step (score) estimate of the effect on the logit (log-odds) scale. It is close to the maximum-likelihood estimate for small effects.

**Saddlepoint approximation.** With few cases (or few controls), the score of a rare variant is far from normal, and $$p_\text{normal}$$ can be much too small. For SNPs with $$p_\text{normal}$$ below `--spa-cut` (default `0.05`), `GLMM` computes $$p$$ by the saddlepoint approximation of SAIGE (Zhou et al. 2018). The score is written as $$T = \sum_r \tilde g_r (y_r - \mu_r)$$ over the records, with the record codes adjusted for the covariates, $$\tilde g = g - X(X'WX)^{-1}X'Wg$$, and the $$y_r$$ treated as independent Bernoulli($$\mu_r$$). $$T$$ is rescaled from its variance $$x'Sx$$ (which includes the GRMs) to the Bernoulli variance $$\tilde g' W \tilde g$$; SAIGE uses one variance ratio estimated from random markers, while `GLMM` uses the exact ratio of each SNP. For these SNPs, `se` is recomputed from the SPA p-value, $$\operatorname{se} = |\hat\beta| / \sqrt{\chi^2_{1, p}}$$, so that `beta`, `se` and `p` agree. For all other SNPs, $$p = p_\text{normal}$$. `--spa-cut 0` turns the approximation off and `--spa-cut 1` applies it to every SNP.

In a simulation with 3,000 individuals, 2% cases and 4,000 null rare variants (MAF 0.002–0.02), the normal approximation gave 21 SNPs with p < 1e-3 and 9 with p < 1e-4 (expected 4 and 0.4); SPA gave 5 and 0. With 50% cases the two agreed.

## 2. Inputs

### Phenotype / covariate table (`--data`)

A whitespace- or tab-separated table with a header row. The first column is the individual ID. It is matched against the IID column of the `.grm.id` files and of the `.fam` file. Columns are selected by name.

| Option | Meaning |
| --- | --- |
| `--data` | Phenotype / covariate table (**required**) |
| `--trait` | Column name of the binary trait (**required**). Give exactly one column. Values must be `0` or `1`. |
| `--covar` | Column name(s) of continuous covariates |
| `--class` | Column name(s) of categorical covariates |
| `--missing-data` | Tokens treated as missing in the table (default `NA NaN nan N/A ...`) |

Records are kept when the ID is in the GRMs and none of the used columns is missing. Every kept ID must also be in the `.fam` file. Linearly dependent fixed-effect columns are removed automatically, and the log lists them. A trait with only cases or only controls is rejected.

### GRM input (`--grm`, `--mgrm`)

Give one of:

| Option | Meaning |
| --- | --- |
| `--grm` | One GRM prefix, for example `--grm output/test` for `output/test.grm.bin` and `output/test.grm.id`. The component is named `G` |
| `--mgrm` | One or more GRM prefixes, one random effect each (for example additive and dominance GRMs). The components are named after the GRMs, as in [REMMA](../05_REMMA/) |

All GRMs must list the same individuals in the same order. Build them with `fastgma --make-grm`; GCTA GRMs are also accepted. See [GRM](../03_GRM/).

### Genotype input (`--bfile`, `--ffile`)

| Option | Meaning |
| --- | --- |
| `--bfile` | PLINK binary genotype prefix (`.bed/.bim/.fam`) |
| `--ffile` | Feature-matrix prefix (`.fam/.bim/.fmat`). Each `.fmat` row is one feature and each column one sample, for imputed dosages or any other continuous per-sample matrix |
| `--missing-geno` | Tokens treated as missing in the `.fmat` file (default `NA NaN nan N/A ...`) |

Give exactly one of `--bfile` and `--ffile`. A `--ffile` feature is centred by its mean, missing values get the mean, `--maf` is not applied, and `af`, `af_case` and `af_ctrl` are feature means.

## 3. Outputs

All files share the `--out` prefix and are tab-separated with a header row. With `--snp-set` the prefix becomes `<out>.<start>_<end>` (0-based start, exclusive end); with `--parallel m i` it becomes `<out>.<m>_<i>`.

| File | Contents |
| --- | --- |
| `<out>.var` | Polygenic variance(s) with standard errors (logit scale) |
| `<out>.sln` | Fixed-effect estimates of the null model (logit scale); with `--predict`, also the predicted genetic values (BLUP) and their reliabilities |
| `<out>.res` | With `--predict`: the linear predictor of each record split into its fixed and random parts, and the fitted probability |
| `<out>.assoc` | One row per SNP with the association results |
| `<out>.log` | Log of the run (appended; `<out>.<m>_<i>.log` with `--parallel`, `<out>.<start>_<end>.log` with `--snp-set`) |

### `<out>.var`

The header is `component variance se`, followed by one row per GRM component (`G` with `--grm`) with the variance and its SE on the logit scale. The SE comes from the inverse AI matrix of the working model at the final linearization; a variance on the boundary (`0`) has the SE `NA`. There is no `residual` row, because the working residual covariance is $$W^{-1}$$ (dispersion fixed at 1). No `.h2` file is written.

### `<out>.sln`

The columns are `type term levels effect_logit se`, and with `--predict` also `rel`, as in [uvLMM](../04_UVLMM/) but with the estimate column named for its scale: `type` is `0` for the intercept (`mu`), `1` for a covariate and `2` for a class level; a dependent column is listed with effect `0` and se `NA`. The estimates and their SEs come from the final linearization (the $$y^*$$, $$W$$ and $$V$$ of the SNP test): $$\hat b = (X' V^{-1} X)^{-1} X' V^{-1} y^*$$ with SEs from $$(X' V^{-1} X)^{-1}$$.

With `--predict`, type `3` rows hold the predicted genetic values of every GRM individual (including those without records, predicted through their relationships), one block per GRM. `term` is the component as named in `<out>.var` and `levels` the individual ID. As in uvLMM, they come from the working model of the final linearization: $$\hat g_k = \sigma^2_k G_k Z' P y^*$$ on the logit scale, `se` is the square root of the prediction error variance $$\mathrm{PEV}_k = \sigma^2_k \operatorname{diag}(G_k) - \sigma^4_k \operatorname{diag}(G_k Z'PZ G_k)$$, and `rel` is the reliability $$1 - \mathrm{PEV} / (\sigma^2_k G_{jj})$$ (`NA` on the fixed-effect rows and for a variance on the boundary). Numbers keep 10 significant digits.

### `<out>.res`

With `--predict`, one row per record: `id observed fixed random mu`. `observed` is the `0`/`1` trait, `fixed` is $$x_i'\hat b$$ and `random` the sum of the predicted genetic values of the record, both on the logit scale, and `mu` $$= 1 / (1 + e^{-(\text{fixed} + \text{random})})$$ is the fitted probability. At convergence the mean of `mu` equals the fraction of cases (the score equation of the intercept).

### `<out>.assoc`

```
chr	snp	pos	a1	a2	af	N	N_case	N_ctrl	af_case	af_ctrl	beta	se	p_normal	p
```

| Column | Meaning |
| --- | --- |
| `chr`, `snp`, `pos` | Chromosome, SNP ID and position from the `.bim` file |
| `a1`, `a2` | Counted (effect) allele and the other allele |
| `af` | Frequency of `a1` in the analysed individuals |
| `N` | Number of analysed individuals with a non-missing genotype |
| `N_case`, `N_ctrl` | Number of case and control records with a non-missing genotype |
| `af_case`, `af_ctrl` | Frequency of `a1` over those case and control records |
| `beta`, `se` | Additive effect (logit scale) and its standard error (from the SPA p-value when SPA was applied) |
| `p_normal` | p-value of the score test by the normal (chi-square) approximation |
| `p` | Final p-value: SPA when `p_normal` is below `--spa-cut`, otherwise `p_normal` |

The effect, SE and p-value columns are `nan` when the SNP fails the `--maf` or `--missing-rate` filter or gives a non-finite statistic. If the saddlepoint equation cannot be solved for a SNP, `p` keeps `p_normal`; the log counts these SNPs. With `--p-cut` below `1`, a row is written when `p` is at or below the cut.

### Convergence

The log reports each PQL iteration and whether the null model converged. Two warnings can appear:

- **Outer loop not converged.** The PQL loop reached `--maxiter0` iterations. The results use the last estimates and may be unreliable; increase `--maxiter0`.
- **Inner loop not converged.** The inner AI-REML loop also stops after `--maxiter0` iterations. The log reports in how many outer iterations this happened.

With few cases the polygenic variance is often estimated at the boundary: in a simulation with 3,000 individuals and 2% cases, `G` was set to `0` in the first iteration and the null model converged after 2 iterations. With `G = 0` the null model is the logistic regression, and the SNP test is its score test (with SPA).

## 4. Options

### Main options

| Option | Default | Meaning |
| --- | --- | --- |
| `--data` | — | Phenotype / covariate table (**required**) |
| `--trait` | — | Binary trait column, exactly one (**required**) |
| `--grm` / `--mgrm` | — | One GRM prefix, or one or more GRM prefixes; exactly one is **required** |
| `--bfile` / `--ffile` | — | PLINK binary prefix, or feature-matrix prefix; exactly one is **required** |
| `--out` | — | Output prefix (**required**) |
| `--covar`, `--class` | none | Covariate and class columns |
| `--missing-data` | `NA NaN nan N/A ...` | Missing-value tokens in the data file |
| `--missing-geno` | `NA NaN nan N/A ...` | Missing-value tokens in the `.fmat` file |
| `--maf` | `0.01` | Keep SNPs with $$\min(p, 1 - p)$$ at or above this, in `[0, 0.5]` (`--bfile` only) |
| `--missing-rate` | `0.05` | Keep SNPs with missing rate at or below this |
| `--spa-cut` | `0.05` | Saddlepoint approximation for SNPs with `p_normal` below this; `0` = off, `1` = every SNP |
| `--predict` | off | Also write the predicted genetic values and their reliabilities to `<out>.sln`, and `<out>.res` |
| `--snp-set <first> <last>` | — | Test SNPs from the first to the last SNP ID (inclusive) |
| `--parallel <m> <i>` | — | Split the SNPs into `m` parts and test part `i` (1-based) |
| `--p-cut` | `1` | Write SNPs whose `p` is at or below this; `1` writes all |
| `--threads` | `10` | Number of MKL/OpenMP threads |

`--snp-set`, `--parallel` and `--start-pos`/`--end-pos` are mutually exclusive.

### Advanced

| Option | Default | Meaning |
| --- | --- | --- |
| `--npart-snp` | `20` | Minimum number of chunks the tested SNPs are read in; chunks are kept to about 128 MB of genotype codes (memory only; does not change results) |
| `--start-pos` | `0` | First SNP to test (0-based `.bim` index) |
| `--end-pos` | `-1` | Stop before this SNP index; `-1` (or a value past the last SNP) means the end of the file |
| `--maxiter0` | `100` | Max iterations, used for both the outer PQL loop and the inner AI-REML loop |
| `--cc-par0` | `1e-7` | Convergence criterion for the change in variances and in fixed effects |
| `--cc-gra0` | `1e-6` | Convergence criterion for the REML score in the inner loop |

Set the environment variable `FASTGMA_LOG_LEVEL=debug` to log the inner AI-REML iterations.

## 5. Examples

### Build the GRM and run the scan

```bash
fastgma --make-grm --bfile test --out output/test

fastgma --glmm \
  --grm output/test \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait disease \
  --covar age \
  --class sex \
  --threads 4 \
  --out output/test_glmm
```

### Additive and dominance GRMs

```bash
fastgma --make-grm --grm-type dom --bfile test --out output/test_dom

fastgma --glmm \
  --mgrm output/test output/test_dom \
  --bfile test \
  --data test_simu_pheno.txt \
  --trait disease \
  --out output/test_glmm_mgrm
```

### Split a scan into parallel jobs

Every job refits the same null model, so the results match a single full run.

```bash
for i in 1 2 3 4; do
  fastgma --glmm --grm output/test --bfile test \
    --data test_simu_pheno.txt --trait disease \
    --parallel 4 $i --out output/test_glmm
done

(head -1 output/test_glmm.4_1.assoc
 tail -n +2 -q output/test_glmm.4_{1,2,3,4}.assoc) > output/test_glmm.assoc
```


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/07_MVLMM/
================================================================
---
layout: page
permalink: /documentation/07_MVLMM/
title: mvLMM
description: Multivariate linear mixed model for several correlated traits
order: 7
---

`MVLMM` (multivariate linear mixed model) is the multi-trait GWAS engine of `fastGMA`. For two or more quantitative traits it first fits a **multivariate null model** to estimate the genetic and residual covariance matrices between the traits (and from them the heritabilities and genetic correlations), then runs a **genome-wide scan** that tests each SNP jointly on all traits.

It is invoked through the top-level dispatcher with `--mvlmm`:

```bash
fastgma --mvlmm [options]
```

Run `fastgma -h 6` for the full option list.

## 1. What it does

With $$t$$ traits, the records of all traits are modelled jointly:

$$
\begin{aligned}
y &= Xb + Zu + e \\
\operatorname{Var}(u) &= \Sigma_g \otimes G \\
\operatorname{Var}(e) &= \Sigma_e \otimes I
\end{aligned}
$$

$$y$$ stacks the records of all traits, trait by trait. $$Xb$$ holds the fixed effects of each trait (an intercept plus that trait's covariates and class variables). $$u$$ is the polygenic effect with GRM $$G$$, and $$e$$ is the residual. $$\Sigma_g$$ and $$\Sigma_e$$ are $$t \times t$$ covariance matrices: their diagonals are the genetic and residual variances of the traits, and their off-diagonals the covariances between traits. Each extra GRM (`--mgrm`) or non-genetic random effect (`--random`) adds one more $$t \times t$$ covariance matrix.

Each trait uses its own non-missing records, so the traits do not need to be measured on the same individuals (with `--solver general`). Records of different traits from the same individual share the residual covariance.

**Null model.** All covariance matrices are estimated by REML with the average-information (AI) algorithm. When the AI step would leave the valid parameter space (a covariance matrix that is not positive definite), the update is chosen, by the REML log-likelihood, among an AI step blended with the EM step, the AI step shortened to 90% of the distance to the boundary, and an AI step in the Cholesky parameterisation. This lets the fit reach estimates on the boundary, such as a genetic correlation near ±1 or a variance near 0. If a covariance matrix ends close to singular, the log says so. Iteration stops when the parameter change is below `--cc-par0`, the gradient norm is below `--cc-gra0`, or the change in log-likelihood is below `--cc-logl0`, or after `--maxiter0` iterations.

**SNP test.** For each SNP, the centred count of `a1` is added as a fixed effect with one coefficient per trait, and the covariance matrices are kept at their null-model values. The $$t$$ effects are tested jointly ("no effect on any trait") with a score/Wald statistic that follows a chi-square distribution with $$t$$ degrees of freedom.

**Approximate test.** The exact test needs $$F = S'PS$$ (the $$t \times t$$ information matrix of the SNP effects) for every SNP. `MVLMM` approximates it by $$F \approx D\, \Gamma\, D$$, where $$D = \operatorname{diag}(\sqrt{x_a'x_a})$$, $$x_a$$ is the centred SNP over the records of trait $$a$$ (the same $$x$$ for every trait with complete data, so that $$F \approx (x'x)\, \Gamma$$), and $$\Gamma$$ is the mean of $$D^{-1} F D^{-1}$$ over `--num-random-snp` random SNPs that pass the MAF and missing-rate filters. $$F$$ does not depend on the phenotypes, so no SNP is excluded from the calibration for being associated. The random SNPs are drawn from the whole genotype file with `--seed` (default `42`), so $$\Gamma$$ is the same in reruns and in every `--parallel` / `--snp-set` job. SNPs with approximate p below `--p-approx-cut` are re-tested exactly. `--p-approx-cut 1` gives the exact test for every SNP; `--p-approx-cut 0` keeps the approximate test only.

> **Choosing `--p-approx-cut`.** The exact test is computed in batches and is cheap for samples below about 10,000 individuals: use `--p-approx-cut 1` there. The approximation works well in large samples of distantly related individuals and less well in small or structured samples. In the 1,062-individual mouse example, the approximate and exact `-log10 p` correlated at 0.94 and the only SNP with exact p < 1e-3 was not re-tested at the default cut (it was at `--p-approx-cut 0.01`). In a human sample of 427 individuals the correlation was 0.99996. If you use the approximate test on structured data, choose a larger cut such as `0.01`.

## 2. Solvers

The fitting algorithm is chosen with `--solver` (default `general`):

| `--solver` | Method | Notes |
| --- | --- | --- |
| `general` (default) | Dense AI-REML on the covariance matrix of all records (one row per record and trait) | Allows missing trait values, several GRMs (`--mgrm`), non-genetic random effects (`--random`), per-trait covariates and `--predict` |
| `eigen` | The GRM is eigendecomposed once ($$G = U \Lambda U'$$) and $$y$$ and $$X$$ are rotated by $$U'$$. The covariance matrix becomes block diagonal with one $$t \times t$$ block $$\lambda_i \Sigma_g + \Sigma_e$$ per individual, so REML costs $$O(n)$$ per evaluation | Needs complete (balanced) data and exactly one `--grm`. Supports `--refit` |

The two solvers maximise the same likelihood. On data that both can fit they give the same `.var`, `.sln`, `.h2`, `.rg` and `.assoc` files.

Rules checked at start-up. Each violation stops the run with an error; the eigen-solver errors suggest `--solver general`:

- `--solver eigen` needs every individual with records to be observed for every trait (with its covariates). Remove the incomplete individuals from the data file, or use `--solver general`.
- `--solver eigen` needs exactly one GRM given with `--grm` (no `--mgrm`).
- `--solver eigen` does not allow non-genetic random effects. `--random` must be omitted or name only the ID column.
- `--solver eigen` needs the same `--covar` and `--class` columns for all traits.
- `--predict` is not available with `--solver eigen`.
- `--refit` is accepted only with `--solver eigen`.

**Per-SNP refit (`--solver eigen --refit`).** Each SNP that is tested exactly (every SNP with `--p-approx-cut 1`) gets its own REML fit of $$\Sigma_g$$ and $$\Sigma_e$$ with the SNP in the model, starting from the null estimates. The `p` column is then the joint Wald test with the re-estimated covariance matrices. An approximate likelihood-ratio p-value `p_lrt` (ML log-likelihoods evaluated at the REML estimates of the two models) and the convergence diagnostics `cc_par` / `cc_gra` are added. Without `--refit` the covariance matrices stay at the null estimates for every SNP.

## 3. Inputs

### Data table (`--data`)

A whitespace- or tab-separated table with a header row. The first column must be the individual ID. It must match the IID (second column) of the `.grm.id` file and of the PLINK `.fam` file. Rows whose ID is not in the GRM(s) are dropped.

Columns are selected by **header name**. A token that is not a header name but a non-negative integer is read as a 0-based column index (column `0` is the ID).

| Option | Meaning |
| --- | --- |
| `--data` | Data file (**required**) |
| `--trait` | Column names of the traits (**required**; two or more, one is allowed) |
| `--covar` | Continuous covariates. Give once for all traits, or once per trait in `--trait` order |
| `--class` | Categorical covariates. Give once for all traits, or once per trait in `--trait` order |
| `--random` | ID column of each random effect (see below). Default: the ID column for every GRM |
| `--missing-data` | Tokens treated as missing in the data file (default `NA NaN nan N/A ...`; giving the option replaces the list) |

An intercept is always fitted for each trait. Per-trait covariates: `--trait y1 y2 --covar E1 E2 --covar E1` fits `E1` and `E2` for `y1` and `E1` for `y2`. `--covar E1` alone fits `E1` for both traits. `--covar` or `--class` given a number of times other than 1 or the number of traits is an error. Fixed-effect columns that are linearly dependent on earlier ones (for example one level of each class variable) are dropped and reported in the log.

A record of a trait is dropped when the trait, one of that trait's covariates or class variables, or one of the `--random` columns is missing. The other traits of the same individual are kept.

Each individual may have **at most one record per trait**. A repeated ID within a trait stops the run with an error; average or remove repeated records first.

### GRMs and random effects (`--grm`, `--mgrm`, `--random`)

| Option | Meaning |
| --- | --- |
| `--grm` | Prefix of a GRM. Its covariance block is named `G` |
| `--mgrm` | Prefixes of further GRMs (for example a dominance GRM). Each block is named after the prefix basename. `--solver general` only |
| `--random` | One ID column per GRM (in the order: `--mgrm` prefixes first, then `--grm`), then the columns of non-genetic random effects |

At least one of `--grm` / `--mgrm` is required. The first entries of `--random` link the GRMs to a data column (normally the ID column). Further entries are non-genetic random effects, such as a cage or batch column, with independent levels. A non-genetic block is named after its column.

| Model | Options |
| --- | --- |
| One GRM | `--grm output/test` |
| Dominance + additive GRM | `--mgrm output/test_dom --grm output/test` |
| GRM + non-genetic effect in column `cage` | `--grm output/test --random 0 cage` |

In the last row, `0` is the ID column (its header name works too). When `--random` is given it must list at least one column per GRM.

GRMs are referenced by prefix and must be in GCTA format (`.grm.bin`, `.grm.id`). GRMs built by `fastgma --make-grm` and by GCTA are both accepted; see [GRM](../03_GRM/).

### Genotype input

| Option | Meaning |
| --- | --- |
| `--bfile` | PLINK binary genotype prefix (`.bed/.bim/.fam`) |
| `--ffile` | Feature-matrix prefix (`.fam/.bim/.fmat`). Each `.fmat` row is one feature and each column one sample, for imputed dosages, expression or any other continuous per-sample matrix |
| `--missing-geno` | Tokens treated as missing in the `.fmat` file (default `NA NaN nan N/A ...`) |

Give exactly one of `--bfile` and `--ffile`. Every individual with records must be in the `.fam` file. A `--ffile` feature is centred by its mean, missing values get the mean, `--maf` is not applied and `af` is the feature mean. The SNP scan always runs; there is no variance-components-only mode.

## 4. Outputs

All files share the `--out` prefix and are tab-separated with a header row. Values that are not finite print as `nan` (`NA` in `.sln`).

| File | Contents |
| --- | --- |
| `<out>.var` | Estimated (co)variance components with standard errors |
| `<out>.var_cov` | Sampling covariance matrix of the components (inverse AI matrix) |
| `<out>.h2` | Heritability of each trait, SE by the delta method |
| `<out>.rg` | Genetic correlation of each pair of traits, SE by the delta method |
| `<out>.sln` | Fixed-effect estimates; with `--predict`, also the predicted random effects |
| `<out>.assoc` | One row per SNP with the joint association test |
| `<out>.log` | Log of the run (appended) |

### `<out>.var`

One row per element of the lower triangle of each covariance matrix:

```
component	row	col	variance	se
G	1	1	1.14819	0.16287
G	2	1	0.674167	0.113107
G	2	2	0.687829	0.11156
residual	1	1	0.961924	0.0498522
residual	2	1	0.212369	0.0385441
residual	2	2	1.15062	0.0563373
```

| Column | Meaning |
| --- | --- |
| `component` | `G` for the `--grm` effect, the prefix basename for each `--mgrm` GRM, the column name for each non-genetic `--random` effect, and `residual` last |
| `row`, `col` | Trait numbers (1-based, in `--trait` order). `row` = `col` = $$i$$ is the variance of trait $$i$$; `row` = $$i$$, `col` = $$j$$ is the covariance between traits $$i$$ and $$j$$ |
| `variance`, `se` | Estimate and standard error (square root of the diagonal of the inverse AI matrix at the estimates) |

The blocks appear in the order `--mgrm` GRMs, `--grm`, non-genetic effects, `residual`.

### `<out>.var_cov`

The full sampling covariance matrix of the `.var` estimates, in the same order. The header row is `parameter` followed by one label per parameter, and each row starts with its label. Labels are `<component>_<row>_<col>`, for example `G_1_1`, `G_2_1`, `residual_2_2`. Values are written with 10 significant digits, as in `.sln` (`.var`, `.h2`, `.rg` and `.assoc` use 6).

### `<out>.h2` and `<out>.rg`

```
trait	h2	se
y1	0.548342	0.0404541
y2	0.378117	0.0420733
```

```
trait1	trait2	rg	se
y2	y1	0.758612	0.0555244
```

Each GRM block $$m$$ is weighted by $$d_m$$, the mean diagonal element of its GRM over the analysed records (pooled over all traits; with complete data, the mean diagonal of the GRM over the individuals), as in `UVLMM`, `REMMA` and `--fast-reml`; non-genetic effects and the residual have $$d = 1$$. The `h2` of trait $$i$$ is $$\sum_{m \in \text{GRM}} d_m S_{m,ii} \big/ \sum_{m} d_m S_{m,ii}$$, where $$S_{m,ii}$$ is its variance in block $$m$$ and the denominator runs over all blocks. A non-genetic random effect such as `cage` therefore counts in the denominator only. The `rg` between traits $$i$$ and $$j$$ is computed from the weighted sum of the GRM blocks, $$S_g = \sum_{m \in \text{GRM}} d_m S_m$$: $$r_g = S_{g,ij} / \sqrt{S_{g,ii}\, S_{g,jj}}$$; with a single GRM, $$d$$ cancels and `rg` is unchanged. The $$d_m$$ are written to the log. Both SEs use the delta method with `<out>.var_cov`. `.rg` has one row per pair (the later trait in `--trait` order first) and only a header for a single trait.

### `<out>.sln`

```
trait	type	term	levels	effect	se
y1	0	mu	1	1.44219	0.093429
y1	1	cov1	1	0.414169	0.0498185
y1	2	sex	1	-0.347646	0.110603
y1	2	sex	2	0	NA
```

| Column | Meaning |
| --- | --- |
| `trait` | Trait name |
| `type` | `0` intercept, `1` covariate, `2` class level, `3` GRM random effect, `4` non-genetic random effect |
| `term` | Fixed effects: column name (or SNP ID). Random effects: the variance component as named in `<out>.var` |
| `levels` | Class level, or individual ID / level for random effects (`1` for intercept and covariates) |
| `effect`, `se` | Estimate and standard error (for random effects, the square root of the prediction error variance, PEV) |
| `rel` | With `--predict`: reliability $$1 - \mathrm{PEV} / \mathrm{prior\ variance}$$ of a random-effect level; `NA` on fixed-effect rows |

A fixed-effect column dropped as linearly dependent has `effect` `0` and `se` `NA` (here `sex` level `2`). Numbers keep 10 significant digits (the rows above are from an earlier version that printed 6). With `--predict`, one row per trait and level is added for each random effect (`type` 3 or 4). For GRM random effects, `--predict` also keeps the GRM individuals without records, so their genetic values are predicted too.

### `<out>.assoc`

One row per SNP. For the traits `y1` and `y2`:

```
chr	snp	pos	a1	a2	af	N	beta_y1	se_y1	beta_y2	se_y2	cov_y2_y1	p
```

With `--solver eigen --refit`:

```
chr	snp	pos	a1	a2	af	N	cc_par	cc_gra	beta_y1	se_y1	beta_y2	se_y2	cov_y2_y1	p_lrt	p
```

| Column | Meaning |
| --- | --- |
| `chr`, `snp`, `pos`, `a1`, `a2` | Chromosome, SNP ID, position, effect allele and other allele from the `.bim` file |
| `af` | Frequency of `a1` among the individuals with records (the feature mean with `--ffile`) |
| `N` | Number of individuals with records and a non-missing genotype |
| `cc_par`, `cc_gra` | `--refit` only: norm of the last parameter change and of the gradient of the per-SNP REML fit (`nan` for SNPs that were not refitted) |
| `beta_<trait>`, `se_<trait>` | Effect of one copy of `a1` on the trait, and its standard error |
| `cov_<ti>_<tj>` | Sampling covariance between the effect estimates of traits `ti` and `tj`, for every pair with `ti` after `tj` in `--trait` order (for three traits: `cov_y2_y1 cov_y3_y1 cov_y3_y2`) |
| `p_lrt` | `--refit` only: approximate likelihood-ratio p-value (chi-square, `t` df) |
| `p` | Joint p-value (chi-square, `t` df) |

For a SNP that was not re-tested exactly, `beta`, `se`, `cov` and `p` come from the approximate test. SNPs that fail `--maf` or `--missing-rate` and SNPs whose test is singular (for example monomorphic among the records of one trait) have `nan` in the effect and p-value columns.

With `--p-cut` below `1`, only SNPs with a finite `p` at or below the cut are written. The default `1` writes every SNP, including the `nan` rows.

## 5. Key options

### Null model

| Option | Default | Meaning |
| --- | --- | --- |
| `--solver` | `general` | `general` or `eigen` (see Section 2) |
| `--maxiter0` | `100` | Max iterations for the null-model fit |
| `--continue` | off | Start from the estimates in an existing `<out>.var` (ignored with a warning if they do not fit the model or are not positive definite) |
| `--predict` | off | Also write the predicted random effects to `<out>.sln` (`--solver general`) |
| `--threads` | `10` | Number of MKL/OpenMP threads |

### SNP scan

| Option | Default | Meaning |
| --- | --- | --- |
| `--maf` | `0.01` | Minimum minor allele frequency `min(af, 1 − af)` (`--bfile` only) |
| `--missing-rate` | `0.05` | Keep a SNP if its missing rate is at or below this |
| `--p-approx-cut` | `1e-3` | SNPs with approximate p below this are re-tested exactly. `0` = approximate test only; `1` = exact test for every SNP |
| `--refit` | off | `--solver eigen`: re-estimate the covariance matrices for every exactly tested SNP |
| `--p-cut` | `1` | In `[0, 1]`. `1` writes every SNP (including `nan` rows); below `1`, only SNPs with a finite `p` at or below the cut |

### Restricting / parallelizing the scan

| Option | Meaning |
| --- | --- |
| `--snp-set <first> <last>` | Test only the SNPs from the first to the last SNP ID (inclusive). The output prefix becomes `<out>.<start>_<end>` (0-based start, exclusive end) |
| `--parallel <m> <i>` | Split the SNPs into `m` equal parts and test part `i` (1-based). The output prefix becomes `<out>.<m>_<i>` |
| `--start-pos` / `--end-pos` | 0-based index of the first SNP and one past the last SNP (`-1` = up to the last SNP). Advanced; output prefix unchanged |

`--snp-set` and `--parallel` cannot be used together or with `--start-pos` / `--end-pos`. Every job fits the same null model and writes all output files (`.var`, `.var_cov`, `.h2`, `.rg`, `.sln`, `.assoc`) to its own prefix. The log is `<out>.<m>_<i>.log` for `--parallel`, `<out>.<start>_<end>.log` for `--snp-set` and `<out>.log` otherwise. Because `--seed` is fixed, all jobs use the same $$\Gamma$$, and the merged `.assoc` files match a single full run.

### Advanced (rarely needed)

Grouped under `Advanced` in `--help`.

| Option | Default | Meaning |
| --- | --- | --- |
| `--num-random-snp` | `2000` | Random SNPs for the approximate-test calibration (capped at the number of SNPs) |
| `--seed` | `42` | Random seed for drawing the calibration SNPs |
| `--npart-snp` | `20` | Minimum number of chunks the tested SNPs are read in; chunks are kept to about 128 MB of genotype codes (memory only; does not change results) |
| `--cc-par0` / `--cc-gra0` | `1e-7` / `1e-6` | Null model: parameter-change / gradient-norm convergence criteria |
| `--cc-logl0` | `5e-5` | Null model: log-likelihood-change convergence criterion |
| `--maxiter` | `30` | `--refit` only: max iterations per SNP |
| `--cc-par` / `--cc-gra` / `--cc-logl` | `1e-7` / `1e-6` / `5e-5` | `--refit` only: per-SNP convergence criteria |

Set the environment variable `FASTGMA_LOG_LEVEL=debug` to log per-block progress.

## 6. Examples

Build the GRM first (see [GRM](../03_GRM/)):

```bash
fastgma --make-grm --bfile test --out output/test
```

In the examples, `test_simu_pheno.txt` has the ID in its first column, the traits `y1`, `y2` and `y3`, and the covariates `E1` and `E2`.

### Two-trait GWAS with one GRM

```bash
fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --out output/test_mvlmm
```

### Exact test for every SNP, write only SNPs with p ≤ 1e-3

```bash
fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --p-approx-cut 1 \
  --p-cut 1e-3 \
  --out output/test_mvlmm_exact
```

### Eigen solver with a per-SNP refit (complete data)

```bash
fastgma --mvlmm --solver eigen --refit \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --p-approx-cut 1 \
  --out output/test_mvlmm_eigen
```

Every individual with a record must have both traits and `E1`.

### Per-trait covariates, three traits

`E1` and `E2` for `y1`, `E1` for `y2` and `E2` for `y3`:

```bash
fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 y3 \
  --covar E1 E2 --covar E1 --covar E2 \
  --out output/test_mvlmm_cov
```

`--class` is given the same way.

### GRM plus a non-genetic random effect

```bash
fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --random 0 cage \
  --out output/test_mvlmm_cage
```

`<out>.var` then has the blocks `G`, `cage` and `residual`.

### Dominance and additive GRMs, with predicted random effects

```bash
fastgma --make-grm --bfile test --grm-type dom --out output/test_dom

fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --mgrm output/test_dom \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --predict \
  --out output/test_mvlmm_mgrm
```

`<out>.var` then has the blocks `test_dom`, `G` and `residual`, and `.h2` / `.rg` combine `test_dom` and `G`.

### Split a scan into 10 parallel jobs

```bash
fastgma --mvlmm \
  --data test_simu_pheno.txt \
  --grm output/test \
  --bfile test \
  --trait y1 y2 \
  --covar E1 \
  --parallel 10 1 \
  --out output/test_mvlmm
```

Repeat with part indices `2` through `10`, then concatenate `output/test_mvlmm.10_*.assoc` (keeping one header).

## 7. Run time and memory

`--solver general` builds dense matrices with one row per record, $$N$$ rows in total (up to $$tn$$ for $$t$$ traits and $$n$$ individuals). It holds about two $$N \times N$$ matrices plus the GRMs themselves (8 bytes per element), and each REML iteration costs $$O(N^3)$$. It is suited to data sets with up to a few tens of thousands of records. `--solver eigen` costs one $$O(n^3)$$ eigendecomposition and then $$O(n)$$ per REML iteration; the exact SNP test costs $$O(n^2)$$ per SNP for the rotation.

Wall time with 4 threads on a shared node (exact test for every SNP unless noted):

| Data | `--solver general` | `--solver eigen` |
| --- | --- | --- |
| Mouse, 2 traits with missing values (2,358 records, 1,295 individuals, 1,407 SNPs) | 5.7 s | — |
| Mouse, 2 traits, complete cases (1,062 individuals, 1,407 SNPs) | 6.4 s | 0.97 s |
| Human, 2 traits (427 individuals, 256,868 SNPs) | 22.0 s | 12.3 s |
| Human, as above, default `--p-approx-cut 1e-3` | 6.9 s | 5.4 s |

For the univariate counterpart, see [uvLMM](../04_UVLMM/).


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/08_LONGWAS/
================================================================
---
layout: page
permalink: /documentation/08_LONGWAS/
title: lonGWAS
description: Longitudinal data analysis with random regression models
order: 8
---

`LonGWAS` is the longitudinal GWAS module of `fastGMA`. It analyses a trait measured repeatedly on the same individuals over time, such as body weight recorded every week. It first fits a **random regression null model** to estimate the genetic and permanent environmental covariance functions. It then tests every SNP for an effect on the **trajectory** of the trait, not only on its mean.

It is invoked through the top-level dispatcher with `--longwas`:

```bash
fastgma --longwas [options]
```

Run `fastgma -h 7` (or `fastgma --longwas --help`) for the full option list.

**Quick start.** With a GRM and a long-format data file (one row per record):

```bash
fastgma --make-grm --bfile example/mouse_long/plink --out output/mouse_long
fastgma --longwas --bfile example/mouse_long/plink --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt --trait trait --tpoint week --class Sex \
  --out output/mouse_longwas
```

This writes the variance components (`.var`), the heritability at each time point (`.h2`), the fixed effects (`.sln`) and one association test per SNP (`.assoc`); `--h2-cum` adds the heritability of the trait summed over a period (`.h2cum`). Section 7 has more examples.

## 1. Model

Let $$y_{ij}$$ be record $$j$$ of individual $$i$$, taken at time $$t_{ij}$$. Time is mapped to $$[-1, 1]$$ and expanded in normalized Legendre polynomials $$\phi_k$$. `LonGWAS` fits

$$
\begin{aligned}
y_{ij} = {} & \textstyle\sum_{k=0}^{f} b_k \phi_k(t_{ij}) && \text{fixed mean curve} \\
& + c && \text{time-independent fixed effects} \\
& + \textstyle\sum_{k=0}^{f} s_k x_i \phi_k(t_{ij}) && \text{SNP effect on the curve (tested)} \\
& + \textstyle\sum_{k=0}^{a} u_{ik} \phi_k(t_{ij}) && \text{additive genetic random regression} \\
& + \textstyle\sum_{k=0}^{p} \mathit{pe}_{ik} \phi_k(t_{ij}) && \text{permanent environmental random regression} \\
& + e_{ij}
\end{aligned}
$$

$$
u \sim N(0, K_a \otimes G), \qquad \mathit{pe} \sim N(0, K_p \otimes I), \qquad e \sim N(0, \sigma^2_e I)
$$

where $$x_i$$ is the centred genotype of individual $$i$$, $$G$$ is the additive GRM, and $$f$$, $$a$$ and $$p$$ are the orders set by `--forder`, `--aorder` and `--porder`. $$K_a$$ and $$K_p$$ are the $$(a+1) \times (a+1)$$ and $$(p+1) \times (p+1)$$ covariance matrices of the genetic and permanent environmental coefficients, and $$\sigma^2_e$$ is one residual variance for all records. `--time-covar` and `--time-class` add further fixed curves to the mean curve, and $$c$$ holds the time-independent fixed effects of `--covar` and `--class`.

The Legendre basis is

$$
\phi_k(x) = \sqrt{\frac{2k + 1}{2}}\, P_k(x), \qquad x = \frac{2(t - t_{\min})}{t_{\max} - t_{\min}} - 1
$$

with $$P_k$$ the Legendre polynomial of degree $$k$$, and $$t_{\min}$$, $$t_{\max}$$ the smallest and largest time point among the analysed records. Because time is rescaled to $$[-1, 1]$$, only the relative spacing of the time points matters: weeks 1–16 and days 7–112 give the same fit.

### Covariance functions and heritability

With $$\phi_a(t) = (\phi_0(t), \ldots, \phi_a(t))^\top$$, the additive genetic variance at time $$t$$ is $$\phi_a(t)^\top K_a \phi_a(t)$$, and the genetic covariance between two time points is $$\phi_a(t_1)^\top K_a \phi_a(t_2)$$ (likewise for the permanent environment with $$K_p$$). `LonGWAS` writes $$K_a$$, $$K_p$$ and $$\sigma^2_e$$ to `<out>.var`, and the heritability at each observed time point,

$$
h^2(t) = \frac{d\, \phi_a(t)^\top K_a \phi_a(t)}{d\, \phi_a(t)^\top K_a \phi_a(t) + \phi_p(t)^\top K_p \phi_p(t) + \sigma^2_e},
$$

to `<out>.h2`. Here $$d$$ is the mean diagonal element of the GRM as fitted (including `--val`) over the records, as in the other modules; it puts the additive genetic variance on the scale of the phenotypic variance of a record. It is written to the log. The genetic correlation between two time points follows from `<out>.var` as

$$
r_g(t_1, t_2) = \frac{\phi_a(t_1)^\top K_a \phi_a(t_2)}{\sqrt{\phi_a(t_1)^\top K_a \phi_a(t_1)\; \phi_a(t_2)^\top K_a \phi_a(t_2)}}.
$$

**Cumulative heritability (`--h2-cum`).** For a trait summed over a period, as 305-day milk yield is the sum of the daily yields of a test-day model, let $$S$$ be the $$m$$ observed time points in $$[t_1, t_2]$$, $$w_a = \sum_{t \in S} \phi_a(t)$$ and $$w_p = \sum_{t \in S} \phi_p(t)$$. The cumulative breeding value of an individual is $$w_a^\top u_i$$ (the sum of its genetic values at the time points of $$S$$), and

$$
h^2_{\text{cum}} = \frac{d\, w_a^\top K_a w_a}{d\, w_a^\top K_a w_a + w_p^\top K_p w_p + m\, \sigma^2_e},
$$

with independent residuals at the time points. The mean over the period gives the same value, and $$m = 1$$ gives $$h^2(t)$$. The standard errors of $$h^2(t)$$ and $$h^2_{\text{cum}}$$ are both from the delta method with the inverse AI matrix.

## 2. Fitting and testing

### Null model

The null model (no SNP) is fitted once by REML with a weighted AI/EM algorithm (WEMAI). At each iteration it uses the smallest EM weight that keeps $$K_a$$ and $$K_p$$ positive definite and $$\sigma^2_e$$ positive; if no blend does, the EM step is halved until the update is valid. The fit stops when one of the three convergence criteria (`--cc-par0`, `--cc-gra0`, `--cc-logl0`) is met; if it does not converge within `--maxiter0` iterations, the log warns and the last estimates are still written and used.

**Starting values.** By default, the diagonals of $$K_a$$ and $$K_p$$ and the residual variance start at the phenotypic variance divided by `aorder + porder + 3`, and the covariances at `0`. With `--continue`, the `variance` column of an existing `<out>.var` is used instead, provided the file has the expected rows for the chosen `--aorder`/`--porder` and gives positive definite matrices; otherwise the defaults are used with a warning. Without `--continue`, an existing `<out>.var` is not read and is overwritten.

### Solvers

The two solvers read different data layouts (Section 4):

| `--solver` | Data file | Method |
| --- | --- | --- |
| `general` (default) | **Long format**: one row per record (`--trait`, `--tpoint`). Individuals may have different numbers of records and different time points | REML on the mixed model equations, with the permanent-environment equations absorbed |
| `eigen` | **Wide format**: one row per individual and one column per time point (`--trait w1:w16`, `--tpoint-values`). Every individual needs a value in every time column | The GRM is eigendecomposed once and the data rotated, so the covariance matrix becomes block diagonal with $$T \times T$$ blocks ($$T$$ time points) |

On the same balanced data both solvers give the same estimates and p-values. `eigen` is faster and needs less memory: on the example data (1,212 mice, 16 weeks, 4 threads) it takes 2.8 s and 124 MB against 15.9 s and 585 MB for `general` on the long version of the same data. `--solver eigen` does not allow `--predict` or `--use-vmat`; it is the only solver that supports `--refit`.

**Data size.** `general` keeps two dense matrices whose size grows with the square of the number of individuals: the absorbed mixed-model equations, of dimension $$p_X + (a+1)N$$ ($$p_X$$ fixed-effect columns, $$N$$ individuals in the GRM), and $$F'PF$$ for the SNP test, of dimension $$(f+1)n$$ ($$n$$ individuals with records). Each takes $$8 d^2$$ bytes for dimension $$d$$; with the default orders this is about 7 GB each for 10,000 individuals. The number of records hardly matters. `eigen` needs $$O(n^2)$$ memory for the GRM and its eigenvectors only.

### SNP test

The SNP scan uses the null-model variance components (unless `--refit` is given). For each SNP the $$f + 1$$ regression coefficients $$s_0, \ldots, s_f$$ are estimated jointly, and $$H_0\colon s_0 = \cdots = s_f = 0$$ is tested with a Wald chi-square test with $$f + 1$$ degrees of freedom, so a SNP that changes the level, the slope or the shape of the curve can be detected.

The exact test needs $$F'PF$$ for the SNP design $$F$$; `LonGWAS` approximates it by $$\Gamma\, x'x$$, where $$x$$ is the centred genotype and the $$(f+1) \times (f+1)$$ matrix $$\Gamma$$ is calibrated on `--num-random-snp` random SNPs (drawn with `--seed`, default `42`) that pass the filters and give a well-defined test. $$F'PF / x'x$$ does not depend on the phenotype, so no p-value filter is applied (as in `UVLMM` and `MVLMM`). At least 10% of the random SNPs must qualify, or the run stops. SNPs with approximate p below `--p-approx-cut` are re-tested exactly; `--p-approx-cut 1` gives the exact test for every SNP. With `--solver eigen --refit`, each exactly tested SNP also gets its own REML fit of the variance components with the SNP in the model.

## 3. Model terms

| Term | Set by | Default | Description |
| --- | --- | --- | --- |
| Mean curve | `--forder` | `2` | Fixed regression of order $$f$$ on time (always included) |
| Time-varying fixed effects | `--time-covar`, `--time-class` | none | Each covariate, or each level of a class variable, gets its own fixed curve of order $$f$$ |
| Time-independent fixed effects | `--covar`, `--class` | none | One effect per covariate or class level, constant over time |
| SNP effect | `--forder` | `2` | Fixed regression of order $$f$$, tested with $$f + 1$$ degrees of freedom |
| Additive genetic effect | `--aorder` | `2` | Random regression of order $$a$$, covariance $$K_a \otimes G$$ |
| Permanent environmental effect | `--porder` | `2` | Random regression of order $$p$$, covariance $$K_p \otimes I$$ |
| Residual | — | — | One homogeneous residual variance |

Order `0` means a constant (intercept-only) term; order `1` adds a linear term, order `2` a quadratic term, and so on. The same `--forder` is used for the fixed curves and for the SNP effect. The highest order must be smaller than the number of distinct time points. The number of (co)variance parameters, $$(a+1)(a+2)/2 + (p+1)(p+2)/2 + 1$$, grows quadratically with `--aorder` and `--porder`.

### Choosing the orders

- **`--aorder`, `--porder`.** Fits with the same `--forder` (and the same fixed effects) can be compared by their REML log-likelihood, printed in the log as `logL = ... at the estimates`, for example with a likelihood-ratio test or AIC. On the example data, `--aorder 1 --porder 1` (7 parameters) gives logL −14,192.5 and the default `--aorder 2 --porder 2` (13 parameters) −13,768.6, so the quadratic random regressions fit far better. REML log-likelihoods with different `--forder` are not comparable, because the fixed effects differ.
- **High orders need much information.** With `--forder 4 --aorder 4 --porder 4` (31 parameters), neither solver converges in 100 iterations on the example data: the estimate of $$K_a$$ has eigenvalues 1.82, 0.078, 0.0099, 0.0019 and 0.0006, so the highest-order coefficients are hardly identified and every iteration needs the EM blend to stay positive definite. Raise the orders step by step and keep those the data support.
- **`--forder`** sets both the shape of the mean curve and the degrees of freedom of the SNP test; more coefficients detect more complex SNP effects but spread the test over more degrees of freedom.

## 4. Inputs

### Data file (`--data`)

The file is whitespace- or tab-separated and its first line is a header. **Columns are selected by name.** A 0-based column number is still accepted when no header column has that name, with a deprecation warning.

**Long format (`--solver general`).** One row per record, so an individual measured at 16 time points has 16 rows. `--trait` names the trait column and `--tpoint` the time column.

```
ID      Sex  week  trait
1419    1    1     4.9
1419    1    2     8.4
1419    1    3     10.5
...
1422    2    1     4.6
```

**Wide format (`--solver eigen`).** One row per individual and one column per time point. `--trait` lists the time columns in time order, either by name or as a range of adjacent columns (`w1:w16`, as in `fastREML`), and `--tpoint-values` gives a file with the time of each of these columns: one number per line, no header, in the same order.

```
ID    Sex  w1   w2   w3    ...  w16          weeks.txt:  1
1419  1    4.9  8.4  10.5  ...  ...                      2
1422  2    4.6  8.3  10.9  ...  ...                      ...
...                                                      16
```

The number of time values must equal the number of `--trait` columns, and the values must be distinct numbers. Every individual needs a value in every time column: a missing value stops the run, naming the individual and the column. Individuals whose time-independent covariates are missing are dropped before this check.

| Option | Meaning |
| --- | --- |
| `--data` | Data file (**required**) |
| `--id-col` | Column of the individual ID (default: the first column) |
| `--trait` | Long format: the trait column, exactly one (**required**). Wide format: the time columns in time order, e.g. `w1:w16` |
| `--tpoint` | Long format: the time point column (**required**); must be numeric |
| `--tpoint-values` | Wide format: file with the time of each `--trait` column, one value per line, no header (**required**) |
| `--covar` | Time-independent continuous covariates |
| `--class` | Time-independent class variables |
| `--time-covar` | Continuous covariates with a time-varying effect |
| `--time-class` | Class variables with a time-varying effect |
| `--missing-data` | Tokens treated as missing in the data file (default `NA NaN nan N/A ...`) |

Notes:

- Records whose ID is not in the GRM are dropped. IDs are matched to the IID of `.grm.id`.
- Long format: the records do not need to be sorted. A record with a missing value in any used column is dropped. A non-numeric time point that is not a missing token stops the run.
- Long format: time points are compared exactly, so values that differ only by rounding (for example `0.5` and `0.5000000001` from an export) are different time points. `LonGWAS` warns when two time points differ by less than $$10^{-6}$$ of the time range; round the time column in that case.
- Dependent fixed-effect columns (for example the last level of a class variable, which is confounded with the mean curve) are removed automatically and reported in the log.

### GRM input (`--grm`)

| Option | Meaning |
| --- | --- |
| `--grm` | Prefix of the additive GRM (**required**) |
| `--val` | Value added to the GRM diagonal to keep it positive definite (default `0.001`) |

As elsewhere in `fastGMA`, the GRM is referenced by prefix (`--grm output/test` reads `output/test.grm.bin` and `output/test.grm.id`). Build it with `fastgma --make-grm`; GCTA GRMs are also accepted. See [GRM](../03_GRM/). `--val` is added to the GRM diagonal when the GRM is read, for both solvers. If the GRM is still not positive definite (for example duplicated or near-identical individuals), the run stops and suggests a larger `--val`.

### Genotype input (`--bfile` / `--ffile`)

| Option | Meaning |
| --- | --- |
| `--bfile` | PLINK binary genotype prefix (`.bed/.bim/.fam`). The SNP effect is estimated per copy of the `a1` allele (column 5 of the `.bim` file) |
| `--ffile` | Feature-matrix prefix (`.fam/.bim/.fmat`). Each `.fmat` row is one feature and each column one sample, whitespace-separated, for example imputed dosages. Features are centred by their mean and tested like SNPs; the `--maf` filter does not apply, and `af` in `<out>.assoc` is the feature mean |
| `--missing-geno` | Tokens treated as missing in the `.fmat` file (default `NA NaN nan N/A ...`; PLINK BED encodes missing calls itself) |

Give at most one of `--bfile` and `--ffile`. Without either, only the null model is fitted and no `.assoc` is written.

## 5. Outputs

All files share the `--out` prefix and are tab-separated with a header row.

| File | Contents |
| --- | --- |
| `<out>.var` | Estimated (co)variance components of the null model with standard errors |
| `<out>.h2` | Heritability at each observed time point, with standard errors |
| `<out>.h2cum` | With `--h2-cum`: cumulative heritability over each requested period, with standard errors |
| `<out>.sln` | Fixed-effect estimates; with `--predict`, also the predicted random regression coefficients |
| `<out>.assoc` | One row per SNP with the association results (only with `--bfile` or `--ffile`) |
| `<out>.log` | Log of the run (appended; see below for `--snp-set` / `--parallel`) |

With `--snp-set` the output prefix becomes `<out>.<start>_<end>` (0-based start, exclusive end); with `--parallel m i` it becomes `<out>.<m>_<i>`. The suffix applies to every file of that job (`.var`, `.h2`, `.h2cum`, `.sln`, `.assoc` and `.log`).

### `<out>.var`

The lower triangle of $$K_a$$ (`G`), then the lower triangle of $$K_p$$ (`permanent`), then the residual variance. `row` and `col` are 1-based indices of the Legendre coefficients (`1` = $$\phi_0$$, `2` = $$\phi_1$$, ...). Rows with `row = col` are variances; the others are covariances. Standard errors come from the inverse AI matrix; values that are not finite print as `NA`. With the default orders, from the [standard scan](#standard-scan) example:

```
component	row	col	variance	se
G	1	1	1.50281	0.378199
G	2	1	0.519915	0.148035
G	2	2	0.241273	0.0679531
G	3	1	-0.202554	0.0586236
G	3	2	-0.0605725	0.022651
G	3	3	0.0397128	0.0117012
permanent	1	1	7.89243	0.365713
permanent	2	1	3.07146	0.150173
permanent	2	2	1.56857	0.0756402
permanent	3	1	-0.923908	0.0588248
permanent	3	2	-0.315397	0.0253611
permanent	3	3	0.155919	0.0147515
residual	1	1	1.46686	0.0190305
```

### `<out>.h2`

One row per distinct observed time point, in increasing order; `h2` is defined in Section 1, and the SE is from the delta method with the inverse AI matrix. From the same run:

```
time	h2	se
1	0.044509	0.0156255
2	0.0559842	0.0156134
3	0.0776433	0.0188936
...
16	0.121616	0.0295892
```

### `<out>.h2cum`

One row per `--h2-cum` interval: `from to n_times h2 se var_g var_pe var_e`, where `n_times` is the number of observed time points in the interval and `var_g` ($$d\, w_a^\top K_a w_a$$), `var_pe` and `var_e` ($$m \sigma^2_e$$) are the variance components of the summed trait. From the [cumulative heritability](#cumulative-heritability) example:

```
from	to	n_times	h2	se	var_g	var_pe	var_e
1	16	16	0.154694	0.0342144	182.8	975.416	23.4697
9	16	8	0.146791	0.033394	82.1053	465.494	11.7349
```

### `<out>.sln`

```
type	term	levels	effect	se	rel
```

`rel` is written with `--predict` only.

| `type` | `term` | `levels` |
| --- | --- | --- |
| `0` | Mean curve: `mu_t0`, `mu_t1`, ... (one per Legendre coefficient) | `1` |
| `1` | Covariate name; for `--time-covar`, `<name>_t0`, `<name>_t1`, ... | `1` |
| `2` | Class variable name; for `--time-class`, `<name>_t0`, `<name>_t1`, ... | Class level |
| `3` | Additive genetic coefficient `G_t0`, `G_t1`, ... (with `--predict`) | Individual ID |
| `4` | Permanent environmental coefficient `permanent_t0`, `permanent_t1`, ... (with `--predict`) | Individual ID |

A fixed-effect column removed as dependent is listed with effect `0` and se `NA`. With `--predict`, the `G_tk` rows cover every individual in the GRM, including those without records, and `rel` is the reliability of each coefficient, $$1 - \mathrm{PEV} / (K_a[k,k]\, G_{jj})$$ for `G_tk` and $$1 - \mathrm{PEV} / K_p[k,k]$$ for `permanent_tk`; it is `NA` on the fixed-effect rows. Numbers keep 10 significant digits. The genetic value of an individual at time $$t$$ is $$\sum_k \hat u_{ik} \phi_k(t)$$. Example (default orders, `--class Sex`, from the [prediction](#predict-the-random-regression-coefficients) example):

```
type	term	levels	effect	se	rel
0	mu_t0	1	24.22101506	0.09258007857	NA
0	mu_t1	1	5.592146947	0.03832815573	NA
0	mu_t2	1	-2.072129875	0.01685774373	NA
2	Sex	1	-0.3747445202	0.06230991615	NA
2	Sex	2	0	NA	NA
3	G_t0	1419	0.6258066228	0.5522730352	0.6848795603
...
```

### `<out>.assoc`

For the default `--forder 2`, from the [standard scan](#standard-scan) example:

```
chr	snp	pos	a1	a2	af	N	beta_t0	se_t0	beta_t1	se_t1	beta_t2	se_t2	cov_t1_t0	cov_t2_t0	cov_t2_t1	p
1	UNC010515443	3658709	G	W	0.495462	1212	0.378564	0.228033	0.224081	0.098128	-0.0403298	0.0408635	0.0187006	-0.00664458	-0.00205035	0.143166
...
```

With `--solver eigen --refit`, the columns `cc_par cc_gra` follow `N`. The number of `beta_t*`, `se_t*` and `cov_t*` columns depends on `--forder`.

| Column | Meaning |
| --- | --- |
| `chr`, `snp`, `pos`, `a1`, `a2` | Chromosome, SNP ID, position and alleles from the `.bim` file; effects are per copy of `a1` |
| `af` | Frequency of `a1` in the analysed individuals (with `--ffile`: the feature mean) |
| `N` | Number of analysed individuals with a non-missing genotype |
| `cc_par`, `cc_gra` | Convergence diagnostics of the per-SNP refit (norm of the last parameter change and of the gradient); `nan` for SNPs not tested exactly |
| `beta_tk`, `se_tk` | Estimated SNP effect on Legendre coefficient $$k$$ and its standard error |
| `cov_tk_tl` | Sampling covariance of `beta_tk` and `beta_tl` |
| `p` | Wald chi-square p-value with `forder + 1` degrees of freedom |

**The SNP effect at a time point.** The effect of one copy of `a1` at time $$t$$ is $$\beta(t) = \sum_k \beta_k \phi_k(t)$$, with $$\phi_k$$ the Legendre basis of Section 1 (using $$t_{\min}$$, $$t_{\max}$$ of the analysed records, which are written to the log). Its variance is $$\phi(t)^\top \Sigma_\beta\, \phi(t)$$, where $$\Sigma_\beta$$ has `se_tk`$$^2$$ on the diagonal and `cov_tk_tl` off the diagonal, so the SNP effect curve and a pointwise confidence band follow from one `.assoc` row.

SNPs that fail the `--maf` or `--missing-rate` filter are written with `nan` results. For SNPs that were not re-tested exactly, the effects, SEs, covariances and `p` are the approximate values. With the default `--p-cut 1` every SNP is written, including those with `nan` p; with `--p-cut` below `1`, only SNPs with a finite `p` at or below the cut are written.

## 6. Key options

### Model

| Option | Default | Meaning |
| --- | --- | --- |
| `--solver` | `general` | `general` (long format) or `eigen` (wide format) (Section 2) |
| `--forder` | `2` | Order of the fixed mean curve, the time-varying fixed effects and the SNP effect |
| `--aorder` | `2` | Order of the additive genetic random regression |
| `--porder` | `2` | Order of the permanent environmental random regression |

The orders must be `0` or larger (see [Choosing the orders](#choosing-the-orders)).

### Null model

| Option | Default | Meaning |
| --- | --- | --- |
| `--maxiter0` | `100` | Maximum number of iterations |
| `--continue` | off | Start from the estimates in an existing `<out>.var` |
| `--predict` | off | Also write the predicted random regression coefficients to `<out>.sln` (`--solver general` only) |
| `--h2-cum <from> <to>` | — | Heritability of the trait summed over the observed time points in `[from, to]`, with SE, to `<out>.h2cum`; may be given several times |
| `--threads` | `10` | Number of MKL/OpenMP threads |

### SNP scan

| Option | Default | Meaning |
| --- | --- | --- |
| `--maf` | `0.01` | Minimum minor allele frequency (`--bfile` only) |
| `--missing-rate` | `0.05` | Maximum SNP missing rate |
| `--p-approx-cut` | `1e-3` | SNPs with approximate p below this are re-tested exactly; `1` = exact test for every SNP |
| `--refit` | off | `--solver eigen` only: re-estimate the variance components for every exactly tested SNP |
| `--p-cut` | `1` | In `[0, 1]`. Write only SNPs with a finite `p` at or below this; `1` writes every SNP |
| `--seed` | `42` | Random seed for sampling the calibration SNPs |
| `--snp-set <first> <last>` | — | Test SNPs from the first to the last SNP ID (inclusive) |
| `--parallel <m> <i>` | — | Split the SNPs into `m` parts and test part `i` (1-based) |

`--snp-set`, `--parallel` and `--start-pos`/`--end-pos` are mutually exclusive.

### Advanced

| Option | Default | Meaning |
| --- | --- | --- |
| `--start-pos` | `0` | 0-based index of the first SNP to test |
| `--end-pos` | `-1` | 0-based index one past the last SNP to test; `-1` (or a value beyond the last SNP) tests to the end |
| `--npart-snp` | `20` | Minimum number of chunks the tested SNPs are read in; chunks are kept to about 128 MB of genotype codes (memory only; does not change results) |
| `--num-random-snp` | `2000` | Random SNPs for the approximate-test calibration; capped at the number of SNPs |
| `--use-vmat` | off | `--solver general`: build the full phenotypic covariance matrix `V` (records × records) for the test instead of using the mixed-model equations. Needs much more memory |
| `--cc-par0` | `1e-7` | Null model: convergence when the norm of the parameter change is below this |
| `--cc-gra0` | `1e-6` | Null model: convergence when the norm of the gradient is below this |
| `--cc-logl0` | `5e-5` | Null model: convergence when the change in log-likelihood is below this |
| `--maxiter` | `30` | Per-SNP refit (`--solver eigen --refit`): max iterations |
| `--cc-par` / `--cc-gra` / `--cc-logl` | `1e-7` / `1e-6` / `5e-5` | Per-SNP refit convergence criteria |

Set the environment variable `FASTGMA_LOG_LEVEL=debug` to log per-chunk progress.

## 7. Examples

The examples use the mouse body-weight data in `example/mouse_long/` (1,212 mice weighed in weeks 1–16; Gray et al. 2015, *Genetics* 201: 213–228) and are run from the repository root:

| File | Format | Contents |
| --- | --- | --- |
| `phe.balance.txt` | Long, `ID Sex week trait` | 19,392 records: every mouse has one record in each of the 16 weeks |
| `phe.unbalance.txt` | Long, `ID Sex week trait` | 15,499 records: a subset of `phe.balance.txt` with about 20% of the records removed, leaving 6–16 records per mouse |
| `phe.wide.txt` | Wide, `ID Sex w1 ... w16` | The data of `phe.balance.txt`, one row per mouse |
| `weeks.txt` | One value per line | The weeks of the columns `w1` ... `w16` (1 to 16), for `--tpoint-values` |

First build the GRM:

```bash
fastgma --make-grm --bfile example/mouse_long/plink --out output/mouse_long
```

### Standard scan

```bash
fastgma --longwas \
  --bfile example/mouse_long/plink \
  --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt \
  --trait trait \
  --tpoint week \
  --class Sex \
  --threads 4 \
  --out output/mouse_longwas
```

This fits quadratic curves (`--forder 2 --aorder 2 --porder 2`, the defaults) with a sex effect that is constant over time. The null model converges in 11 iterations (logL −13,768.6).

### Variance components only

```bash
fastgma --longwas \
  --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt \
  --trait trait \
  --tpoint week \
  --class Sex \
  --out output/mouse_longwas_vc
```

Without `--bfile` or `--ffile`, `LonGWAS` writes `.var`, `.h2` and `.sln` and skips the SNP scan.

### Balanced data with the eigen solver

```bash
fastgma --longwas --solver eigen \
  --bfile example/mouse_long/plink \
  --grm output/mouse_long \
  --data example/mouse_long/phe.wide.txt \
  --trait w1:w16 \
  --tpoint-values example/mouse_long/weeks.txt \
  --class Sex \
  --threads 4 \
  --out output/mouse_longwas_eigen
```

The wide file holds the same data as `phe.balance.txt`, and the [standard scan](#standard-scan) with `--data example/mouse_long/phe.balance.txt` gives the same variance components and p-values with `--solver general` (Section 2 compares time and memory). Add `--refit` to re-estimate the variance components for the exactly tested SNPs.

A missing value in a time column stops the run and names the first one, for example:

```
[error] The data are not balanced: 1 of the 19392 individual x time-point cells (1212 individuals x 16 time points) have no value; e.g. individual 1422 has none in column w5 (time 5). --solver eigen needs a value in every time column for every individual (after removing individuals with a missing covariate). Fill in or remove the missing values, or use the long format with --solver general.
```

### Linear curves with a sex-specific mean curve

```bash
fastgma --longwas \
  --bfile example/mouse_long/plink \
  --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt \
  --trait trait \
  --tpoint week \
  --time-class Sex \
  --forder 1 --aorder 1 --porder 1 \
  --maxiter0 400 \
  --threads 4 \
  --out output/mouse_longwas_linear
```

`--time-class Sex` gives each sex its own fixed linear curve, and the SNP test has 2 degrees of freedom. The fit converges slowly: it needs 313 iterations here, so `--maxiter0` is raised from the default `100`. Stopping at 100 iterations gives a warning, but the variance components are then already within 0.3% of the converged values.

### Cumulative heritability

```bash
fastgma --longwas \
  --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt \
  --trait trait \
  --tpoint week \
  --class Sex \
  --h2-cum 1 16 --h2-cum 9 16 \
  --out output/mouse_longwas_cum
```

This adds `<out>.h2cum` with the heritability of the body weight summed over weeks 1–16 and over weeks 9–16 (0.155 and 0.147, against 0.04–0.13 for the single weeks in `<out>.h2`): summing averages out the residual of the single records.

### Predict the random regression coefficients

```bash
fastgma --longwas \
  --bfile example/mouse_long/plink \
  --grm output/mouse_long \
  --data example/mouse_long/phe.unbalance.txt \
  --trait trait \
  --tpoint week \
  --class Sex \
  --predict \
  --continue \
  --threads 4 \
  --out output/mouse_longwas
```

`<out>.sln` then also holds `G_tk` and `permanent_tk` for every individual. `--continue` starts the fit from `output/mouse_longwas.var` of the first example.

### Split a scan into parallel jobs

```bash
for i in 1 2 3 4; do
  fastgma --longwas \
    --bfile example/mouse_long/plink \
    --grm output/mouse_long \
    --data example/mouse_long/phe.unbalance.txt \
    --trait trait \
    --tpoint week \
    --class Sex \
    --parallel 4 $i \
    --p-cut 0.05 \
    --out output/mouse_longwas
done
```

Each job writes `output/mouse_longwas.4_<i>.assoc` with only the SNPs with `p <= 0.05`. Because `--seed` is fixed, all jobs use the same calibration.


================================================================
Page: https://chaoning.github.io/fastGMA/documentation/09_FASTREML/
================================================================
---
layout: page
permalink: /documentation/09_FASTREML/
title: fastREML
description: REML variance components and heritability for G, GxE and noise-by-environment effects
order: 9
---

`fastREML` estimates genetic (G), gene-by-environment (GxE), noise-by-environment (NxE) and residual variance components, and the corresponding heritabilities, directly from a dense GRM. Small samples are fitted exactly; large samples (up to UK Biobank scale) are fitted without ever forming the GRM in memory, by streaming it from disk.

It is invoked through the top-level dispatcher with `--fast-reml`:

```bash
fastgma --fast-reml [options]
```

Run `fastgma -h 8` for the full option list.

## 1. What it does

For one quantitative trait and $$q$$ interacting environments $$W$$ (an $$n \times q$$ matrix), `fastREML` fits

$$
\begin{aligned}
y &= Xb + g + u_{\text{gxe}} + u_{\text{nxe}} + e \\
V = \operatorname{var}(y) &= \sigma^2_g K + \sigma^2_{\text{gxe}} K_{\text{gxe}} + \sigma^2_{\text{nxe}} N_E + \sigma^2_e I
\end{aligned}
$$

where

| Term | Covariance | Meaning |
| --- | --- | --- |
| $$g$$ | $$\sigma^2_g K$$ | additive genetic effect; $$K$$ is the GRM |
| $$u_{\text{gxe}}$$ | $$\sigma^2_{\text{gxe}} K_{\text{gxe}}$$, with $$K_{\text{gxe}} = K \circ (W W'/q)$$ | gene-by-environment interaction ($$\circ$$ = element-wise product) |
| $$u_{\text{nxe}}$$ | $$\sigma^2_{\text{nxe}} N_E$$, with $$N_E = \operatorname{diag}(W_i W_i'/q)$$ | noise-by-environment: residual variance that changes with the environment |
| $$e$$ | $$\sigma^2_e I$$ | residual |

$$X$$ contains the intercept, the `--covar` and `--class` columns and the environment main effects (`--env-int` columns). The environment columns are standardized (mean 0, variance 1) unless `--no-standardize-env` is given. `--no-noisebye` drops the NxE term, giving a three-component model.

The variance components are estimated by **restricted maximum likelihood (REML)** with the average-information (AI) algorithm. By default, the fixed effects stay in the model:

$$
-2 \log L = \log\lvert V\rvert + \log\lvert X' V^{-1} X\rvert + y' P y, \qquad P = V^{-1} - V^{-1} X (X' V^{-1} X)^{-1} X' V^{-1}
$$

With `--residualize`, the trait is first regressed on $$X$$, and the residuals are then fitted with $$P = V^{-1}$$. This is the fastGxE-compatible variant. Because it ignores the degrees of freedom used by $$X$$, it gives somewhat smaller variance components. On the mouse example below, with 41 fixed-effect columns, G is 0.873 instead of 0.913.

Heritabilities are the proportions of the phenotypic variance, $$V_p = \sigma^2_g d_g + \sigma^2_{\text{gxe}} d_{\text{gxe}} + \sigma^2_{\text{nxe}} d_{\text{nxe}} + \sigma^2_e$$, where $$d_g$$, $$d_{\text{gxe}}$$ and $$d_{\text{nxe}}$$ are the mean diagonal elements of the kernels. The SEs use the delta method with the inverse AI matrix.

## 2. Exact and implicit computation

| Path | Used when | How |
| --- | --- | --- |
| exact | $$n \le$$ `--exact-threshold` (default 8,000) and no `--force-implicit` | $$K$$ is read into memory; $$V$$ is factorized (Cholesky) and inverted every iteration; traces are exact. |
| implicit | otherwise | $$K$$ is never formed. Products $$KB$$ and $$K_{\text{gxe}} B$$ stream the lower triangle of the GRM from disk (or from RAM with `--in-memory`). $$V^{-1} B$$ is solved by block preconditioned conjugate gradients (PCG), traces are estimated with `--num-rand-vec` random probe vectors (Hutchinson estimator), and $$\log\lvert V\rvert$$ is estimated by stochastic Lanczos quadrature. |

The implicit path is stochastic: its estimates vary with `--seed` and `--num-rand-vec` within the probe noise. With a fixed seed and thread count it is exactly reproducible. With the default strict REML it also solves for $$V^{-1} X$$, which adds $$p$$ columns (the number of fixed-effect columns) to each PCG solve.

Starting values come from a randomized method-of-moments (Haseman-Elston) estimate using the same number of probes (`<out>.mom`).

**Boundary estimates.** When the REML optimum of a component is zero (for example NxE when there is no noise heterogeneity), the Newton step would make it negative. That component is then held at a tiny positive value ($$10^{-6}$$ times the sum of the components), and the other components take the Newton step. The log shows `boundary: NxE = ...` on those iterations.

## 3. Inputs

### Data table (`--data`)

A whitespace- or tab-separated table with a header row. The first column must be the individual ID. Only individuals in both the data table and the GRM, with no missing value in the used columns, are analysed, and each individual must appear once.

| Option | Meaning |
| --- | --- |
| `--trait` | Column name of the trait (exactly one) |
| `--covar` | Continuous covariates (column names; ranges such as `pc1:pc10` are allowed) |
| `--class` | Categorical covariates (one-hot coded) |
| `--env-int` | Interacting environment columns (**required**; ranges allowed, e.g. `e1:e40`). These are also fitted as fixed main effects. They must not also be given in `--covar` or `--class`. |
| `--missing-data` | Tokens treated as missing (default `NA NaN nan N/A ...`) |

### GRM (`--grm`)

A dense GRM prefix in GCTA format: `<grm>.grm.bin` (float32, or legacy float64; the file size tells which) and `<grm>.grm.id`. GRMs built by `fastgma --make-grm` and by GCTA are both accepted; see [GRM](../03_GRM/).

## 4. Outputs

All files use the `--out` prefix and are tab-separated with a header row. Values that are not defined are written as `nan`.

| File | Columns |
| --- | --- |
| `<out>.var` | `component variance se diag_mean variance_prop h2 h2_se`: one row each for `G`, `GxE`, `NxE` (absent with `--no-noisebye`) and `residual`, then a `total` row whose `h2`/`h2_se` hold the total genetic heritability (G + GxE). `diag_mean` is the mean kernel diagonal and `variance_prop` the proportion of $$V_p$$. `h2` is defined for the genetic components only. |
| `<out>.mom` | `component mom_variance mom_used diag_mean variance_prop`: the method-of-moments estimate and the positive starting values actually used, plus an `h2_total` row |
| `<out>.iter` | `iter sigma_g2 sigma_gxe2 sigma_nxe2 sigma_e2 neg2logL score_norm delta_norm ai_condition pcg_iters pcg_max_rel_resid`: one row per REML iteration (values at the start of the iteration; `pcg_iters` is 0 and `pcg_max_rel_resid` is `nan` on the exact path) |
| `<out>.diag` | `iid K_diag Kgxe_diag Nxe_diag V_diag`: per-individual kernel diagonals and the diagonal of `V` at the estimates |
| `<out>.samples` | `iid grm_index analysis_index`: the analysed individuals with their GRM row |
| `<out>.log` | the run log |

Example `<out>.var` (mouse example, default strict REML):

```
component	variance	se	diag_mean	variance_prop	h2	h2_se
G	0.912913	0.169888	1.0139	0.241628	0.241628	0.0365331
GxE	0.697459	0.20997	1.01353	0.184534	0.184534	0.0546111
NxE	3.80856e-06	0.619244	1	9.94221e-07	nan	nan
residual	2.1982	0.590571	1	0.573837	nan	nan
total	nan	nan	nan	nan	0.426162	0.0607458
```

## 5. Key options

### Model and likelihood

| Option | Default | Meaning |
| --- | --- | --- |
| `--residualize` | off | fastGxE-compatible: residualize the trait on `X`, then use `P = V^-1` |
| `--no-noisebye` | off | Drop the NxE component |
| `--no-standardize-env` | off | Use the environment columns as given |
| `--log-trait` | off | Log-transform the trait (must be positive) |
| `--int-trait` | off | Rank-based inverse-normal transform of the residualized trait |

### Convergence

| Option | Default | Meaning |
| --- | --- | --- |
| `--maxiter0` | 20 | Maximum REML iterations |
| `--cc-par0` | 1e-4 | Stop when the relative parameter change is below this |
| `--cc-gra0` | 1e-2 | … or the score norm is below this |
| `--cc-logl0` | 1e-2 | … or the change in `-2logL` is below this |

### Computation

| Option | Default | Meaning |
| --- | --- | --- |
| `--threads` | 10 | Number of threads |
| `--exact-threshold` | 8000 | Use the exact path when $$n$$ is at most this |
| `--force-implicit` | off | Always use the implicit path |
| `--num-rand-vec` | 32 | Random probe vectors for the trace and method-of-moments estimates |
| `--seed` | 42 | Seed of the probe vectors |
| `--in-memory` | off | Implicit path: load the kept GRM triangle into RAM once |
| `--grm-single` | off | With `--in-memory`: store the GRM in single precision (half the RAM) |
| `--pcg-tol` / `--pcg-maxiter` | 1e-5 / 500 | PCG relative residual tolerance / iteration limit. An iteration whose PCG solves stop above `--pcg-tol` logs a warning and cannot end REML as converged. |
| `--precond-rank` / `--precond-oversample` | 0 / 10 | Rank (0 = diagonal only) and oversampling of a randomized low-rank preconditioner for `K` |
| `--no-pcg-warmstart` | off | Do not start each iteration's PCG from the previous solution |
| `--logl-lanczos` | 15 | Lanczos steps for the `log|V|` estimate (0 = do not compute `-2logL` on the implicit path) |
| `--matvec-block` | 0 (auto) | Right-hand-side columns per GRM pass; caps implicit-path memory |
| `--dry-run` | off | Read the data, cache the kernel diagonals, time one GRM product, then exit |

## 6. Examples

All examples use the simulated mouse data (`example/mouse/pheno/pheno_gxe_simu.txt`: 1,304 individuals, 40 environments `e1`–`e40`, trait `trait`) and a GRM made with `fastgma --make-grm --bfile example/mouse/geno/plink --out mouse`.

### Default: strict REML, exact path

```bash
fastgma --fast-reml --grm mouse \
  --data example/mouse/pheno/pheno_gxe_simu.txt \
  --trait trait --env-int e1:e40 \
  --maxiter0 50 --threads 2 \
  --out mouse_gxe
```

It converges in 4 iterations (less than 1 s). NxE sits on the boundary; the heritability of G is 0.242 and the total genetic heritability (G + GxE) is 0.426.

### fastGxE-compatible residualized fit

```bash
fastgma --fast-reml --grm mouse \
  --data example/mouse/pheno/pheno_gxe_simu.txt \
  --trait trait --env-int e1:e40 \
  --maxiter0 50 --threads 2 --residualize \
  --out mouse_gxe_res
```

### Implicit (streaming) path, as used for large samples

```bash
fastgma --fast-reml --grm mouse \
  --data example/mouse/pheno/pheno_gxe_simu.txt \
  --trait trait --env-int e1:e40 \
  --maxiter0 50 --threads 2 \
  --force-implicit --num-rand-vec 100 \
  --out mouse_gxe_implicit
```

Add `--in-memory` (and `--grm-single`) when the kept GRM triangle fits in RAM. Otherwise every product re-streams the GRM file.


================================================================
fastgma -h <n>
================================================================
Output of `fastgma -h <n>` for every mode. fastGMA 3077860 (built 2026-10-27, oneMKL 2025.0.2)

==> fastgma -h 1 <==

Compute dense genomic relationship matrices (GRMs) in GCTA format:
lower-triangular
float32 `<out>.grm.bin`, "FID IID" `<out>.grm.id`, per-pair SNP counts
`<out>.grm.N.bin`, plus a text `<out>.grm.N` (SNPs used, scale factor).

Defaults:
--grm-type add
--code-type 2
--maf 0.01
--missing-rate 0.05
--val 0

Quick start:
Additive GRM from PLINK BED
fastgma --make-grm --bfile test --out test

Additive GRM from feature matrix
fastgma --make-grm --ffile test --out test

Dominance GRM (statistical coding)
fastgma --make-grm --grm-type dom --bfile test --out test

Variance-scaled alternative
fastgma --make-grm --bfile test --code-type 1 --out test

Additive GRM with SNP-effect variance proportional to (2pq)^alpha
fastgma --make-grm --bfile test --alpha -0.25 --out test

Leave-one-chromosome-out GRM for chromosome 6, with frequencies saved
fastgma --make-grm --bfile test --not-chr 6 --write-freq --out test_loco6

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
  -m,     --make-grm          Calculate GRM
          --threads INT       Number of threads (default: 10)
          --grm-type TEXT:{add,dom} [add]
                              GRM type (add, dom)
          --bfile TEXT Excludes: --ffile
                              Prefix for PLINK binary PED files
          --ffile TEXT Excludes: --bfile
                              Prefix for feature-by-sample matrix files
          --out TEXT REQUIRED Prefix for output file
          --code-type INT:{1,2} [2]  Excludes: --alpha
                              Genotype coding type (1 or 2)
          --alpha FLOAT Excludes: --code-type
                              Additive GRM: SNP-effect variance proportional to (2pq)^alpha
                              (default: -1)
          --npart INT:POSITIVE [[1,1]]  x 2
                              Partition GRM elements: <npart> <ipart>
          --npart-snp INT:POSITIVE [20]
                              Partition SNPs into n parts
          --maf FLOAT [0.01]  MAF cutoff value
          --missing-rate FLOAT [0.05]
                              Maximum allowed missing rate
          --keep TEXT         Keep only these individuals (one IID, or FID IID, per line)
          --remove TEXT       Remove these individuals (one IID, or FID IID, per line)
          --extract TEXT      Use only the SNPs listed in the first column (--bfile)
          --exclude TEXT      Leave out the SNPs listed in the first column (--bfile)
          --chr TEXT ...      Use only these chromosomes, e.g. 1-22 or 1,3,5 (--bfile)
          --not-chr TEXT ...  Leave out these chromosomes, e.g. 6 for a LOCO GRM (--bfile)
          --read-freq TEXT    Use the A1 frequencies in this file (header with SNP, A1, FREQ;
                              --bfile)
          --write-freq        Write the A1 frequency of every selected SNP to <out>.freq
                              (--bfile)
          --val FLOAT [0]     Value added to the GRM diagonal (default: 0)
          --missing-geno TEXT [[NA,Na,na,NAN,NaN,nan,-NAN,-NaN,-nan,<NA>,<na>,N/A,n/a]]  ...
                              Missing genotype values (space-separated)

==> fastgma -h 2 <==
Quick start:
Merge partitioned GRMs
fastgma --process-grm --merge --grm test --npart 5

PCA
fastgma --process-grm --pca --grm test --out test # test.eigenval, test.eigenvec
(top 20 PCs)
fastgma --process-grm --pca --num-pc 0 --grm test --out test # all PCs

Epistatic GRM from two dense GRMs
fastgma --process-grm --make-epis --prefix1 test1 --prefix2 test2 --out test

Reformat GRM
fastgma --process-grm --reformat --grm test --out-fmt 0 --out test #
test.grm.matrix
fastgma --process-grm --reformat --grm test --out-fmt 1 --out test #
test.grm.index_triplet
fastgma --process-grm --reformat --grm test --out-fmt 2 --out test #
test.grm.iid_triplet
fastgma --process-grm --reformat --grm test --out-fmt 3 --out test #
test.grm.bin
fastgma --process-grm --reformat --sparse --grm test --out-fmt 0 --out test #
test.grm.matrix
fastgma --process-grm --reformat --sparse --grm test --out-fmt 1 --out test #
test.grm.index_triplet
fastgma --process-grm --reformat --sparse --grm test --out-fmt 2 --out test #
test.grm.iid_triplet
fastgma --process-grm --reformat --sparse --grm test --out-fmt 3 --out test #
test.grm.sp.bin

GRM inverse
fastgma --process-grm --make-inv --grm test --out-fmt 0 --out test.inv #
test.inv.grm.matrix (diag + --val, default 0.001)
fastgma --process-grm --make-inv --grm test --out-fmt 1 --out test.inv #
test.inv.grm.index_triplet
fastgma --process-grm --make-inv --grm test --out-fmt 2 --out test.inv #
test.inv.grm.iid_triplet
fastgma --process-grm --make-inv --grm test --out-fmt 3 --out test.inv #
test.inv.grm.bin

Output one IID from each close pair
fastgma --process-grm --out-close --grm test --out test

Group related IIDs
fastgma --process-grm --group --grm test --cut-value 0.05

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
  -d,     --process-grm REQUIRED
                              Process an existing GRM
          --threads INT       Number of threads (default: 10)
          --merge             Merge GRM files
          --pca               Perform PCA
          --make-epis         Calculate epistatic GRM
          --reformat          Reformat GRM
          --make-inv          Calculate the inverse of GRM
          --out-close         Drop individuals until no pair exceeds --cut-value; writes
                              <out>.grm.pruned_iid
          --group             Group individuals by GRM relationships
          --grm TEXT          GRM prefix without the .grm suffix
          --out TEXT          Output prefix without the .grm suffix
          --prefix1 TEXT      First dense GRM prefix without the .grm suffix for --make-epis
          --prefix2 TEXT      Second dense GRM prefix without the .grm suffix for --make-epis
          --npart INT         Partition count (for merging)
          --out-fmt INT [2]   Output format: 0=matrix, 1=index_triplet, 2=iid_triplet, 3=binary
          --cut-value FLOAT [0.05]
                              Relatedness threshold for --out-close and --group
          --val FLOAT         Value added to the diagonal before --make-inv (default: 0.001)
          --num-pc INT        Number of leading PCs written by --pca (0 = all; default: 20)
          --sparse [0]        Indicate sparse GRM

==> fastgma -h 3 <==
Univariate linear mixed model GWAS (UVLMM)

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --threads INT [10]  Number of threads
          --model TEXT:{lmm,lm,multi-grm,repeat} [lmm]
                              Model to fit: lmm (single GRM, default), lm (no random effect),
                              multi-grm (multiple GRMs), repeat (repeated measures)
          --solver TEXT:{general,eigen,mme}
                              Fitting algorithm: eigen (eigendecomposition of the single --grm,
                              one record per individual; default for --model lmm), mme (AI-REML
                              on the mixed model equations, whose size is the number of fixed
                              effects plus random-effect levels, not records; default for
                              --model repeat and the choice for many repeated records; needs
                              positive-definite GRMs) or general (AI-REML on the dense n x n V;
                              default for --model multi-grm). Not used by --model lm.
          --data TEXT         Data file (required)
          --trait TEXT REQUIRED
                              Column name for the trait (single trait only)
          --covar TEXT ...    Column names for covariates
          --class TEXT ...    Column names for class variates
          --random TEXT ...   Column names for non-genetic random effects
          --missing-data TEXT ...
                              Missing value in data file
          --grm TEXT          Prefix for additive GRM file
          --mgrm TEXT ...     Prefixes for multiple GRM files
          --val FLOAT         Value added to the GRM diagonal when it is read (default 0.001
                              for --solver mme, which inverts the GRM and a centred GRM is
                              singular; 0 for the other solvers)
          --bfile TEXT        Prefix for PLINK binary files
          --ffile TEXT        Prefix for feature-matrix files (.fam, .bim, .fmat). Each .fmat
                              row is one feature and each column one sample,
                              whitespace-separated; suitable for imputed dosages, expression,
                              or any continuous per-sample matrix.
          --missing-geno TEXT ...
                              Missing value in genotype file
          --out TEXT REQUIRED Output prefix
          --predict           Output predicted values for random effects
          --continue          Continue from last saved variances
          --maxiter0 INT [100]
                              Max iterations for null model
          --snp-set TEXT x 2  Start and end SNP IDs to test (inclusive)
          --p-approx-cut FLOAT [0.001]
                              SNPs whose approximate p is below this are re-tested exactly (0 =
                              approximate test only; 1 = exact test for every SNP)
          --refit, --no-refit{false}
                              Eigen solver: re-estimate the variance components for every SNP
                              that is re-tested exactly (default on; --no-refit tests at the
                              null variances). --refit --p-approx-cut 1 refits every SNP.
          --p-cut FLOAT [1]   Write only SNPs whose final p-value (p_wald; p for the general
                              solver) is <= this (1 = all)
          --maf FLOAT [0.01]  Minimum minor allele frequency, min(p, 1-p) (PLINK input)
          --missing-rate FLOAT [0.05]
                              Keep a SNP if its missing rate is <= this
          --maxiter INT [30]  Max iterations for SNP testing
          --parallel INT x 2  Partition into m parts and test part i

Advanced:
          --cc-par0 FLOAT [1e-07]
                              Convergence criterion for parameters (null)
          --cc-gra0 FLOAT [1e-06]
                              Convergence criterion for gradient (null)
          --cc-logl0 FLOAT [5e-05]
                              Convergence criterion for the change in log-likelihood (null)
          --cc-par FLOAT [1e-07]
                              Convergence criterion for parameters
          --cc-gra FLOAT [1e-06]
                              Convergence criterion for gradient
          --npart-snp INT [20]
                              Minimum number of chunks the tested SNPs are read in; chunks are
                              kept to about 128 MB (memory only)
          --num-random-snp INT [2000]
                              Random SNPs for approximate-test calibration (auto-capped to
                              #SNPs)
          --seed UINT [42]    Random seed for drawing the calibration SNPs (--num-random-snp);
                              fixed so results are reproducible
          --start-pos INT [0]
                              0-based index of the first SNP to test
          --end-pos INT [-1]  0-based index one past the last SNP to test (-1 = all)

==> fastgma -h 4 <==
REMMA: Random Effect Model for Multiple Associations

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --data TEXT         Phenotype/covariate table (required; header line, first column =
                              individual ID)
          --mgrm TEXT ... REQUIRED
                              GRM file prefix(es) of the null model; one or more
          --bfile TEXT        PLINK BED file prefix (give one of --bfile and --ffile)
          --ffile TEXT        Prefix for feature-matrix files (.fam, .bim, .fmat): each .fmat
                              row is one feature and each column one sample,
                              whitespace-separated; centred by the feature mean, no MAF filter
                              (--add, --epis)
          --missing-geno TEXT ...
                              Missing value tokens in the --ffile matrix (PLINK BED encodes
                              missing calls itself)
          --out TEXT REQUIRED Output file prefix
          --trait TEXT REQUIRED
                              Trait column name (single trait)
          --covar TEXT ...    Covariate column names
          --class TEXT ...    Class (categorical) covariate column names
          --missing-data TEXT ...
                              Missing value tokens in the data file (default: NA NaN nan N/A
                              ...)
          --add               Additive test (exact GLS statistic with the null-model V)
          --dom               Dominance test, statistical coding: joint additive + dominance
                              GLS per SNP
          --epis              Additive x additive epistasis test of SNP pairs (i, j > i) with
                              SNP i in the tested range: approximate scan without main effects,
                              full model with both main effects below --p-approx-cut
          --maf FLOAT [0.01]  Keep SNPs with minor-allele frequency min(p, 1-p) >= maf
          --missing-rate FLOAT [0.05]
                              Keep SNPs with missing rate <= cutoff
          --min-geno-count INT [0]
                              Minimum number of individuals in each genotype class for the
                              dominance test (0 = off)
          --snp-set TEXT x 2  Test SNPs from the first to the second SNP ID (inclusive); epis:
                              pairs whose first SNP is in this set
          --parallel INT x 2  Split the scan into m parts and run part i (1-based): m i; epis
                              parts hold about equal numbers of pairs
          --p-cut FLOAT [1]   Write SNPs/pairs whose final p-value is <= p-cut (dominance:
                              p_add or p_dom); 1 = all
          --p-approx-cut FLOAT [1e-06]
                              epis: refit pairs with both main effects when the approximate p
                              is below this; 1 = full model for every pair
          --solver TEXT:{general,mme} [general]
                              Null-model algorithm: general (AI-REML on the dense n x n V) or
                              mme (AI-REML on the mixed model equations, whose size is the
                              number of fixed effects plus GRM individuals, not records; the
                              choice for many repeated records; needs positive-definite GRMs)
          --val FLOAT         Value added to the GRM diagonal when it is read (default 0.001
                              for --solver mme, which inverts the GRMs; 0 for general)
          --threads INT [10]  Number of threads

Advanced:
          --npart-snp INT [20]
                              Minimum number of chunks the tested SNPs are read in; chunks are
                              kept to about 128 MB (memory only; --add/--dom)
          --start-pos INT [0]
                              First tested SNP (0-based index)
          --end-pos INT [-1]  End of the tested SNPs (0-based, exclusive; -1 = last SNP)
          --num-random-snp INT [10000]
                              epis: random SNP pairs for the gamma calibration (capped at
                              #SNPs)
          --seed UINT [42]    epis: random seed for the gamma SNP-pair sample
          --maxiter0 INT [100]
                              Max iterations for the null model
          --cc-par0 FLOAT [1e-07]
                              Null-model convergence: norm of the parameter change
          --cc-gra0 FLOAT [1e-06]
                              Null-model convergence: norm of the gradient
          --cc-logl0 FLOAT [5e-05]
                              Null-model convergence: change in REML log-likelihood

==> fastgma -h 5 <==
GLMM: Generalized Linear Mixed Model GWAS (binary traits)

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --data TEXT         Phenotype/covariate table (required; header line, first column =
                              individual ID)
          --grm TEXT          GRM file prefix (one GRM; give --grm or --mgrm)
          --mgrm TEXT ...     GRM file prefixes, one random effect each (give --grm or --mgrm)
          --bfile TEXT        PLINK BED file prefix (give one of --bfile and --ffile)
          --ffile TEXT        Prefix for feature-matrix files (.fam, .bim, .fmat): each .fmat
                              row is one feature and each column one sample,
                              whitespace-separated; centred by the feature mean, no MAF filter
          --missing-geno TEXT ...
                              Missing value tokens in the --ffile matrix (PLINK BED encodes
                              missing calls itself)
          --out TEXT REQUIRED Output file prefix
          --predict           Also write the predicted genetic values (BLUP, logit scale) and
                              their reliabilities to <out>.sln, and <out>.res
          --trait TEXT REQUIRED
                              Binary (0/1) trait column name (single trait)
          --covar TEXT ...    Covariate column names
          --class TEXT ...    Class (categorical) covariate column names
          --missing-data TEXT ...
                              Missing value tokens in the data file (default: NA NaN nan N/A
                              ...)
          --maf FLOAT [0.01]  Keep SNPs with minor-allele frequency min(p, 1-p) >= maf
          --missing-rate FLOAT [0.05]
                              Keep SNPs with missing rate <= cutoff
          --snp-set TEXT x 2  Test SNPs from the first to the second SNP ID (inclusive)
          --parallel INT x 2  Split the SNPs into m parts and test part i (1-based): m i
          --p-cut FLOAT [1]   Write SNPs whose p is <= p-cut; 1 = all
          --spa-cut FLOAT [0.05]
                              Saddlepoint approximation for SNPs whose normal-approximation p
                              is below this (0 = off, 1 = every SNP)
          --threads INT [10]  Number of threads

Advanced:
          --npart-snp INT [20]
                              Minimum number of chunks the tested SNPs are read in; chunks are
                              kept to about 128 MB (memory only)
          --start-pos INT [0]
                              First tested SNP (0-based index)
          --end-pos INT [-1]  End of the tested SNPs (0-based, exclusive; -1 = last SNP)
          --maxiter0 INT [100]
                              Max outer (and inner) PQL iterations for the null model
          --cc-par0 FLOAT [1e-07]
                              Null-model convergence: change in variances and fixed effects
          --cc-gra0 FLOAT [1e-06]
                              Null-model convergence (inner loop): norm of the REML score

==> fastgma -h 6 <==
MVLMM: Multivariate Linear Mixed Model GWAS

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --data TEXT         Data file (required; whitespace-separated, header row, first
                              column = individual ID)
          --grm TEXT          Additive GRM file prefix
          --bfile TEXT        PLINK BED file prefix (give one of --bfile and --ffile)
          --ffile TEXT        Prefix for feature-matrix files (.fam, .bim, .fmat): each .fmat
                              row is one feature and each column one sample,
                              whitespace-separated (dosages, expression, ...); centred by the
                              feature mean, no MAF filter
          --out TEXT REQUIRED Output file prefix
          --mgrm TEXT ...     Further GRM file prefixes (general solver; fitted before --grm)
          --missing-data TEXT ...
                              Missing value tokens in the data file (default: NA, NaN, nan,
                              N/A, ...)
          --missing-geno TEXT ...
                              Missing value tokens in the --ffile matrix (PLINK BED encodes
                              missing calls itself)
          --trait TEXT ... REQUIRED
                              Trait column names (two or more traits; one is allowed)
          --covar TEXT ...    Covariate column names; repeat once per trait in --trait order,
                              or give once for all traits
          --class TEXT ...    Class column names; repeat once per trait in --trait order, or
                              give once for all traits
          --random TEXT ...   ID column per random effect: one per GRM (--mgrm first, then
                              --grm), then non-genetic effects (default: the ID column for
                              every GRM)
          --solver TEXT [general]
                              general: any data (missing trait values, several GRMs, other
                              random effects); eigen: eigendecomposition of the GRM, faster,
                              needs complete (balanced) data, one --grm and no other random
                              effects
          --refit             Eigen solver only: re-estimate the variance components for the
                              SNPs that are tested exactly (all SNPs with --p-approx-cut 1) and
                              add an approximate LRT p-value column p_lrt
          --predict           Predict random effects (general solver)
          --continue          Start the null model from the estimates in <out>.var
          --p-approx-cut FLOAT [0.001]
                              SNPs whose approximate p-value is below this are re-tested
                              exactly; 1 = exact test for every SNP, 0 = approximate test only
          --maf FLOAT [0.01]  Minimum minor allele frequency
          --missing-rate FLOAT [0.05]
                              Maximum genotype missing rate (SNP kept if <= cutoff)
          --p-cut FLOAT [1]   Write only SNPs with p <= this value (1 = all)
          --snp-set TEXT x 2  First and last SNP ID (inclusive) to test
          --parallel INT x 2  Split the scan into m parts and run part i (1-based): --parallel
                              m i
          --threads INT [10]  Number of threads

Advanced:
          --num-random-snp INT [2000]
                              Random SNPs for the approximate-test calibration (capped at
                              #SNPs)
          --seed UINT [42]    Random seed for sampling the --num-random-snp SNPs
          --npart-snp INT [20]
                              Minimum number of chunks the tested SNPs are read in; chunks are
                              kept to about 128 MB (memory only)
          --start-pos INT [0]
                              0-based index of the first SNP to test
          --end-pos INT [-1]  0-based index one past the last SNP to test (-1 = all)
          --maxiter0 INT [100]
                              Max iterations for the null model
          --cc-par0 FLOAT [1e-07]
                              Null model: convergence criterion for the parameter change
          --cc-gra0 FLOAT [1e-06]
                              Null model: convergence criterion for the gradient norm
          --cc-logl0 FLOAT [5e-05]
                              Null model: convergence criterion for the logL change
          --maxiter INT [30]  --refit: max iterations per SNP
          --cc-par FLOAT [1e-07]
                              --refit: convergence criterion for the parameter change
          --cc-gra FLOAT [1e-06]
                              --refit: convergence criterion for the gradient norm
          --cc-logl FLOAT [5e-05]
                              --refit: convergence criterion for the logL change

==> fastgma -h 7 <==
LonGWAS: Longitudinal GWAS using random regression

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --data TEXT         Data file (long format: one row per record, header row with
                              column names)
          --grm TEXT          Additive GRM prefix (<prefix>.grm.bin / .grm.id)
          --bfile TEXT        PLINK BED prefix; without --bfile or --ffile only the null model
                              is fitted
          --ffile TEXT        Prefix for feature-matrix files (.fam, .bim, .fmat): each .fmat
                              row is one feature (e.g. an imputed dosage) and each column one
                              sample, whitespace-separated; tested like a SNP
          --val FLOAT [0.001]
                              Value added to the GRM diagonal to keep it positive definite
          --out TEXT REQUIRED Output file prefix
          --solver TEXT [general]
                              general: MME-based REML, any (unbalanced) data; eigen:
                              eigendecomposition of the GRM, needs one record per individual at
                              each of the same time points
          --missing-data TEXT ...
                              Missing value tokens in the data file (default: NA, NaN, nan,
                              N/A, ...)
          --missing-geno TEXT ...
                              Missing value tokens in the --ffile matrix (PLINK BED encodes
                              missing calls itself)
          --id-col TEXT       Column with the individual IDs (default: first column)
          --trait TEXT ... REQUIRED
                              Trait column (long format); with --solver eigen (wide format) the
                              time columns in time order, e.g. w1:w16
          --tpoint TEXT       Time point column (long format, --solver general)
          --tpoint-values TEXT
                              Wide format (--solver eigen): file with the time of each --trait
                              column, one value per line, no header
          --covar TEXT ...    Time-independent covariate columns
          --class TEXT ...    Time-independent class columns
          --time-covar TEXT ...
                              Covariate columns with a time-varying (order --forder) effect
          --time-class TEXT ...
                              Class columns with a time-varying (order --forder) effect per
                              level
          --forder INT [2]    Order of the fixed regressions and of the SNP effect
          --aorder INT [2]    Additive genetic regression order
          --porder INT [2]    Permanent environmental regression order
          --predict           Write the predicted random regression coefficients (G keeps all
                              GRM individuals)
          --h2-cum FLOAT x 2  Heritability of the trait summed over the observed time points in
                              [from, to] (e.g. 305-day yield), with SE, to <out>.h2cum: from
                              to; may be repeated
          --continue          Start the null model from the estimates in <out>.var
          --snp-set TEXT x 2  First and last SNP ID to test (inclusive)
          --parallel INT x 2  Split the SNPs into m parts and test part i (1-based): m i
          --p-approx-cut FLOAT [0.001]
                              SNPs with an approximate p below this are re-tested exactly; 1 =
                              exact test for every SNP
          --refit             (--solver eigen) re-estimate the variance components for every
                              SNP that is tested exactly
          --maf FLOAT [0.01]  Minimum minor allele frequency
          --missing-rate FLOAT [0.05]
                              Maximum SNP missing rate
          --p-cut FLOAT [1]   Write only SNPs with p <= this value (1 = all)
          --seed UINT [42]    Random seed for sampling the calibration SNPs
          --maxiter0 INT [100]
                              Max iterations for null model
          --threads INT [10]  Number of threads

Advanced:
          --start-pos INT [0]
                              0-based index of the first SNP to test
          --end-pos INT [-1]  0-based index one past the last SNP to test (-1 = all)
          --npart-snp INT [20]
                              Minimum number of chunks the tested SNPs are read in; chunks are
                              kept to about 128 MB (memory only)
          --num-random-snp INT [2000]
                              Random SNPs for the approximate-test calibration (capped at the
                              number of SNPs)
          --use-vmat          (--solver general) build the full V matrix for the SNP test
          --cc-par0 FLOAT [1e-07]
                              Null model: convergence when the norm of the parameter change is
                              below this
          --cc-gra0 FLOAT [1e-06]
                              Null model: convergence when the norm of the gradient is below
                              this
          --cc-logl0 FLOAT [5e-05]
                              Null model: convergence when the change in logL is below this
          --maxiter INT [30]  Per-SNP refit (--solver eigen --refit): max iterations
          --cc-par FLOAT [1e-07]
                              Per-SNP refit (--solver eigen --refit): parameter convergence
          --cc-gra FLOAT [1e-06]
                              Per-SNP refit (--solver eigen --refit): gradient convergence
          --cc-logl FLOAT [5e-05]
                              Per-SNP refit (--solver eigen --refit): logL convergence

==> fastgma -h 8 <==
fastREML - dense-GRM variance components and heritability

fastgma [OPTIONS]

OPTIONS:
  -h,     --help              Print this help message and exit
          --fast-reml         Run the fastREML variance-component estimator
          --grm TEXT          Dense GRM prefix (expects <grm>.grm.bin / .grm.id) (required)
          --data TEXT         Phenotype/covariate data file with a header line (required)
          --out TEXT REQUIRED Output prefix
          --trait TEXT ... REQUIRED
                              Trait column name (exactly one)
          --covar TEXT ...    Numeric covariate column names (supports a:b ranges)
          --class TEXT ...    Categorical covariate column names (one-hot)
          --env-int TEXT ... REQUIRED
                              Interacting environment column names
          --threads INT [10]  Number of threads
          --maxiter0 INT [20]
                              Maximum REML iterations
          --num-rand-vec INT [32]
                              Stochastic trace probe vectors (implicit path)
          --seed UINT [42]    Random seed for MoM / trace probe vectors
          --cc-par0 FLOAT [0.0001]
                              Relative parameter-change convergence threshold
          --cc-gra0 FLOAT [0.01]
                              Gradient-norm convergence threshold
          --cc-logl0 FLOAT [0.01]
                              |Delta(-2logL)| convergence threshold (OR-ed with
                              --cc-par0/--cc-gra0)
          --logl-lanczos INT [15]
                              Lanczos steps for the SLQ log|V| / -2logL estimate (implicit
                              path; 0 = off)
          --pcg-tol FLOAT [1e-05]
                              PCG relative residual tolerance
          --pcg-maxiter INT [500]
                              Maximum PCG iterations
          --precond-rank INT [0]
                              Randomized low-rank preconditioner rank (0 = diagonal only)
          --precond-oversample INT [10]
                              Oversampling columns for the low-rank range finder
          --no-pcg-warmstart  Disable warm-starting each REML iteration's PCG from the previous
                              solution
          --exact-threshold INT [8000]
                              Materialize K and run exact REML when n <= this
          --matvec-block INT [0]
                              RHS columns per streaming GRM pass (0 = auto: one pass per PCG
                              round; caps implicit-path memory)
          --force-implicit    Always use the implicit (streaming) path
          --dry-run           Load data, cache kernel diagonals, time one streaming matvec,
                              then exit
          --in-memory         Load the whole GRM into RAM once (avoids re-streaming it per
                              matvec)
          --grm-single        With --in-memory, store the GRM as float (halves RAM, ample for
                              REML)
          --residualize       fastGxE-compatible: residualize the trait on the fixed effects
                              and environments, then fit with P = V^-1 (default: strict REML
                              keeping the fixed effects in P)
          --no-noisebye       Drop the noise-by-environment component (3-component model)
          --log-trait         Log-transform the trait before residualization (requires positive
                              values)
          --int-trait         Rank-based inverse-normal transform the residualized trait (Blom)
          --no-standardize-env
                              Do not standardize the environment columns
          --missing-data TEXT ...
                              Override missing-value labels