Association analysis with glma

Run quantitative, conditional, binary, LOCO and marker-set analyses on a small simulated genotype panel.

This tutorial uses the artificial PLINK panel distributed with gsuite: 96 individuals, 48 markers and six chromosome-local regions. gsim creates two continuous traits from declared marker effects. The small example is intended to explain the interface and results; it is not a power, calibration or performance study.

Prepare the genotype panel

library(gbase)
library(glma)
library(gsim)

example_dir <- system.file("extdata", "summary-statistic", package = "gsuite")
bed <- file.path(example_dir, "reference.bed")
Glist <- gprep(bedfiles = bed)

length(Glist$ids)
sum(lengths(Glist$rsids))

Glist fixes the selected people, markers, alleles and source-file order used by every individual-level analysis below.

Simulate two quantitative traits

Declare a few nonzero effects, then let gsim generate phenotypes from the existing genotype panel.

markers <- unlist(Glist$rsids, use.names = FALSE)
traits <- c("trait_A", "trait_B")
effects <- matrix(0, length(markers), length(traits),
  dimnames = list(markers, traits))
effects[c(3, 11, 27, 35), "trait_A"] <- c(.45, -.30, .35, .25)
effects[c(3, 19, 27, 43), "trait_B"] <- c(.30, .40, -.25, .35)

simulation <- gsim(
  Glist = Glist,
  nt = 2,
  architecture = "fixed",
  beta = effects,
  h2 = c(.45, .35),
  vg = c(.45, .35),
  re = .15,
  standardize_W = FALSE,
  scale_effects = TRUE,
  seed = 20260929,
  compute_sumstats = FALSE,
  chunk_size = 16
)

set.seed(20260930)
covariates <- matrix(rnorm(length(Glist$ids)), ncol = 1,
  dimnames = list(Glist$ids, "batch"))
phenotypes <- simulation$Y + .25 * covariates[, 1]
rownames(phenotypes) <- Glist$ids
colnames(phenotypes) <- traits

Named phenotype rows and covariate rows are aligned to the selected Glist individuals. An intercept is added by glma; do not include one in covariates.

Scan one or several quantitative traits

linear <- glma_linear(
  phenotype = phenotypes,
  Glist = Glist,
  covariates = covariates,
  block_size = 16
)

linear

The phenotype matrix produces one glma_fit per named trait while BED reading and genotype decoding are shared:

<glma_multi_fit> 2 quantitative traits: trait_A, trait_B

The principal result is an ordinary data frame:

association_A <- linear$trait_A$associations
head(association_A[c("marker", "chromosome", "a1", "a2",
  "beta", "se", "statistic", "p", "status")])

successful_A <- subset(association_A, status == "ok")
successful_A[order(successful_A$p), ][1:5, ]

For a linear scan, beta is the allele-1 dosage effect, se its standard error, and statistic the t statistic. Every selected marker remains in the table. A marker that cannot be tested has a non-ok status and unavailable estimates rather than being silently removed.

Condition on selected markers

condition_on uses 1-based positions in selected Glist marker order. The specified genotypes enter the fixed-effect design together.

conditional <- glma_linear(
  phenotype = phenotypes[, "trait_A"],
  Glist = Glist,
  covariates = covariates,
  condition_on = 3L,
  block_size = 16
)

conditional$method
conditional$diagnostics
head(conditional$associations)

This answers whether each tested marker contributes after adjustment for the conditioning marker. It is not a multi-marker omnibus test.

Analyse a binary trait

Create a reproducible binary outcome from the first simulated genetic value and the same covariate.

probability <- plogis(
  -.4 + as.numeric(scale(simulation$G[, 1])) + .25 * covariates[, 1]
)
binary <- rbinom(length(probability), 1, probability)
names(binary) <- Glist$ids

logistic <- glma_logistic(binary, Glist, covariates, block_size = 16)
logistic

head(logistic$associations[c("marker", "case_n", "control_n",
  "beta", "se", "statistic", "p", "status")])

Here beta is an allele-1 log-odds effect and statistic is a Wald z statistic. Separation, sparse outcomes and convergence failures are reported in status. This route is ordinary logistic regression, not binary mixed-model association.

Account for a genomic background with LOCO

Both methods use the joint full-panel engine. The infinitesimal method fixes one Gaussian/ridge prior with equal variance on empirically standardized genotypes; it does not compare elastic or mixture priors. The standalone interface has retired the previous chromosome-specific variance-ratio search. Its chunk_ridge, observation_pcg and marker_coordinates algorithms and variance-search controls are no longer accepted. Use the joint product backends and controls shown here.

infinitesimal <- glma_infinitesimal_loco(
  phenotypes[, "trait_A"],
  Glist,
  covariates,
  algorithm = "decoded_full",
  block_size = 16,
  controls = list(scaling_factor = 1)
)

infinitesimal
infinitesimal$fit[c("genetic_variance", "residual_variance")]
infinitesimal$fit$loco
infinitesimal$diagnostics

The explicit scaling factor keeps this small teaching example deterministic. For a real dataset, omit it to use the native cross-chromosome structure check and ridge-reference calibration when required. Streamed storage supports the decoded backends; packed backends require prepared storage. Threads are local to each operation.

The variance procedure is explicit:

1. Project the phenotype and genotypes on the fixed covariates. Blocked RHE supplies the starting heritability, bounded to [.01, .99]. Strong-marker adjustment affects this initialization only; the final fit retains the original phenotype and all requested markers. Gaussian/ridge-only fitting fixes marker power at -1, giving equal prior weights on empirically standardized markers.
2. Fit a shared three-point heritability grid at half, three-quarters and the starting estimate; at the lower boundary use .01, .03 and .05. Each grid point fits the full phenotype and Gaussian Monte Carlo probes: ten through 40,000 samples, three above that. The grid and the held-out Gaussian fit share genotype scans. Stabilization requires at least four sweeps and changes no greater than .005 in both heritability and normalized validation MSE. The held-out fit supplies this stopping diagnostic; it does not select among priors for the Gaussian method.
3. Linearly interpolate a bracketed zero of the Monte Carlo equation. If the grid does not bracket a zero, retain the initial estimate and report heritability_grid_used_initial = TRUE. Interpolate the phenotype/residual inner product at the accepted heritability h, then set ve = crossprod(y_residual, r(h)) / (N - p) and vg = h * ve / (1 - h), where p is the fixed-design rank including the intercept. This is an approximate stochastic variance estimator.
4. Hold that global vg/ve pair fixed for the global and chromosome-excluded Gaussian refits. Each retained marker has prior variance vg / m_active on the standardized scale. Reuse projected chunk Grams and full-response grid effects for warm starts. fit$loco reports the shared variance pair, active marker count and fitting sweeps for each chromosome.

The association scan tests each marker against its chromosome-excluded residual, after fixed-covariate projection, and returns effects, standard errors, chi-square statistics and p-values. These are scaled residual-offset test estimates; exact mixed-model GLS effects are not asserted. The estimator differs from the retired per-chromosome variance search, so old statistical and timing records do not establish equivalence. Broad mixed-model null-tail qualification remains open.

Diagnostics identify fitting_engine = "joint_full_panel_ridge", variance_estimator = "blocked_rhe_mc_grid", variance_shared_across_chromosomes = TRUE, the marker power, candidate count, heritability and fitting/scan phases.

The alternative prior-selected background chooses among native prediction priors and returns marker effects after the selected global and LOCO refits:

prior_selected <- glma_prior_selected_loco(
  phenotypes[, "trait_A"],
  Glist,
  covariates,
  algorithm = "decoded_full",
  block_size = 16,
  controls = list(
    scaling_factor = 1,
    validation_fraction = .2,
    split_seed = 20260932,
    maximum_sweeps = 200L
  )
)

prior_selected$fit
prior_selected$diagnostics
head(prior_selected$associations)

The explicit scaling_factor = 1 keeps this tiny example deterministic and bypasses automatic structure calibration. It is not a general recommendation for real datasets.

Fit several genetic background components

sets defines disjoint named lists of marker IDs for either LOCO method. Unassigned fitted markers form an automatic background group. This is one joint background model and one association scan; it is distinct from aggregating association evidence over marker sets with VEGAS. The first implementation requires prepared genotypes and supports at most 32 components, including background. Unknown markers, overlaps, empty sets and sets without any retained background marker are rejected.

component_sets <- list(group_A = markers[seq_len(16)])
component_fit <- glma(phenotypes[, "trait_A"], Glist,
  covariates = covariates, sets = component_sets,
  method = "infinitesimal_loco", block_size = 16,
  controls = list(scaling_factor = 1))
component_fit$fit$components
component_fit$fit$loco_components

# Same groups; one ridge/elastic family selected across all groups.
component_prior <- glma(phenotypes[, "trait_A"], Glist,
  covariates = covariates, sets = component_sets,
  method = "prior_selected_loco", block_size = 16,
  controls = list(scaling_factor = 1))

For group c, K_c = Z_c Z_c'/m_c, where Z_c contains empirically standardized marker columns and m_c is its global fitted marker count. The genomic covariance is the sum of vg_c * K_c, plus residual variance times the identity. greml’s existing full-moment multi-component RHE estimates the genetic variance profile and residual component jointly, including cross-component moments. Rank-deficient or ill-conditioned moment systems fail explicitly. Raw RHE estimates and adjusted initial values are retained. A positive ridge floor of 1e-8 times the projected phenotype mean square approximates boundary components and is reported; this is not an exact zero-component estimate.

The shared scalar Monte Carlo grid refines the total genetic scale while keeping the RHE proportions fixed. Both grouped methods fix marker power at -1. Prior-selected LOCO chooses one family across groups with distinct group scales; it does not select a separate family for each set. This approximate procedure is not a multidimensional REML or Monte Carlo variance optimization.

LOCO removes the tested chromosome from every group while retaining global kernel normalizations and variance parameters. It does not redistribute removed prior variance among retained markers. A group wholly on that chromosome has zero active contribution. fit$components gives global marker counts and variance estimates; fit$loco_components gives active counts, global normalization counts, fixed variances and active prior variance sums. For grouped fits, fit$loco$genetic_variance is their retained prior variance sum. Diagnostics name component_rhe_profile_mc_grid and report the ridge floor, RHE moment condition estimate and fixed variance-profile policy. Existing no-set calls retain their previous estimator and normalization. The earlier one-component timings do not measure this extension.

Aggregate marker evidence with VEGAS

VEGAS uses checked summary statistics and sparse LD from the same Glist. Prepare the LD resource once and derive aligned effect/standard-error summaries from the linear scan.

ld_prefix <- tempfile("glma-tutorial-ld-")
LD_Glist <- gprep(
  Glist = Glist,
  task = "sparseld",
  out_prefix = ld_prefix,
  max_distance_bp = 0,
  max_distance_variants = 8,
  r2 = 0,
  block_size = 16,
  nthreads = 1
)

statistics <- data.frame(
  rsids = association_A$marker,
  ea = association_A$a1,
  nea = association_A$a2,
  b = association_A$beta,
  seb = association_A$se
)

sets <- split(markers, rep(paste0("region", 1:6), each = 8))

evidence <- vegas(
  statistics,
  LD_Glist,
  sets = sets,
  threads = 4,
  verbose = TRUE
)

evidence
evidence$stat$trait[c("set", "marker_count", "statistic",
  "score", "p_value", "p_available")]

Each marker set defines one region. The one-call form opens one sparse-LD reader and prepares eigenvalue spectra in bounded batches across marker sets. For repeated analyses, prepare the spectra once and reuse the ordinary serializable R object for any number of checked traits:

prepared <- prepare_vegas(LD_Glist, sets, threads = 4)
saveRDS(prepared, "vegas-reference.rds")

evidence_A <- vegas(statistics, prepared = prepared)
evidence_B <- vegas(statistics_B, prepared = prepared)

Evaluation does not reopen the LD resource or repeat the eigendecompositions. Memory remains bounded by the preparation thread count rather than the total number of sets. gsea::prepare_vegas() and gsea::vegas() re-export the same interfaces for pathway workflows.

For full within-gene LD computed directly from BED, prepare once through gbase and reuse the stored spectra in glma:

eigen_Glist <- gbase::gprep(Glist, task = "eigen", sets = sets, nthreads = 4)
prepared_BED <- prepare_vegas(eigen_Glist)
# Align checked statistics to prepared_BED$markers, including both alleles.
saveRDS(prepared_BED, "vegas-bed-reference.rds")

Each trait supplied to vegas(stat, prepared = prepared_BED) must contain every row of prepared_BED$markers in its exact marker and allele order. This metadata is the union of the requested gene markers. The saved spectra retain full gene membership; removing a marker requires preparation again. Constant positive weights within each gene can rescale the stored spectrum; nonuniform weights require sparse-LD preparation. Sparse-LD cutoffs and windows do not apply to the BED spectra. Tail defaults remain glma’s existing approximations; this route does not assert identical P-values to qgg’s different tail evaluator.

Remove the temporary sparse-LD files when the example is finished:

unlink(unlist(LD_Glist$sparseLD$paths, use.names = FALSE))

Save and report results

All results are ordinary R lists and data frames:

saveRDS(list(linear = linear, conditional = conditional,
  logistic = logistic, infinitesimal = infinitesimal,
  prior_selected = prior_selected, vegas = evidence),
  "glma-analysis.rds")

write.csv(association_A, "trait-A-associations.csv", row.names = FALSE)
write.csv(evidence$stat$trait, "trait-A-marker-sets.csv", row.names = FALSE)

For real analyses, inspect the complete status column, method-specific fit information and diagnostics before filtering by p-value. Multiple-testing correction, Manhattan/QQ plotting and scientific interpretation remain explicit analysis decisions rather than automatic transformations of the native output.

The larger connected genomic analysis shows linear GWAS and VEGAS results for 10,000 simulated individuals and 50,000 markers. The glma validation page records the focused numerical comparisons and their limits.