Public API#

The top-level evo_lmm namespace contains the supported public classes and functions. The signatures below are generated directly from the installed package.

Evolutionary random-effects LMMs backed by raw GRGL dosage operators.

class evo_lmm.AssociationResult(chrom, local_idx, score, beta, se, chisq, pvalue, chisq_linreg=None, pvalue_linreg=None, model_mask=None, frequencies=None, inverse_scale=nan, calibration_factor=nan)#

Bases: object

Compact BOLT-compatible association output for one chromosome block.

beta and se are in raw diploid-dosage effect units, matching the evolutionary model’s sigma_b2 scale. score is the calibrated inverse-variance score x_j^T V_loco^-1 y on the BOLT-normalised test column. Entries excluded by model_mask (monomorphic or covariate-collinear columns) carry nan statistics and pvalue = 1.

Parameters:

  • chrom (Any)

  • local_idx (ndarray)

  • score (ndarray)

  • beta (ndarray)

  • se (ndarray)

  • chisq (ndarray)

  • pvalue (ndarray)

  • chisq_linreg (ndarray | None)

  • pvalue_linreg (ndarray | None)

  • model_mask (ndarray | None)

  • frequencies (ndarray | None)

  • inverse_scale (float)

  • calibration_factor (float)

good()#

Return the boolean mask of variants with reportable statistics.

Return type:

ndarray

class evo_lmm.CalibrationResult(factor, std, ratio_of_medians, median_of_ratios, selected, tried, prospective, retrospective, inverse_scale, residuals=<factory>, quadratic_form=<factory>, screen_threshold=5.0, seed=0)#

Bases: object

Prospective/retrospective calibration of the evolutionary LMM statistic.

Parameters:

  • factor (float)

  • std (float)

  • ratio_of_medians (float)

  • median_of_ratios (float)

  • selected (tuple[tuple[Any, int], ...])

  • tried (int)

  • prospective (ndarray)

  • retrospective (ndarray)

  • inverse_scale (dict[Any, float])

  • residuals (dict[Any, ndarray])

  • quadratic_form (dict[Any, float])

  • screen_threshold (float)

  • seed (int)

factor#

The applied calibration factor: the ratio of prospective to retrospective statistic sums, or the ratio of medians when the jackknife standard error exceeds CALIBRATION_STD_LIMIT.

Type:

float

inverse_scale#

Per-chromosome VinvScaleFactor in raw-effect units. A tested chromosome’s calibrated inverse-variance score is (x_j^T H_loco^-1 P y) / sigma_b2 / inverse_scale[chrom].

Type:

dict[Any, float]

residuals#

H_loco^-1 P_C y for each left-out chromosome, reused by evo_lmm.association() so the solves are not repeated.

Type:

dict[Any, numpy.ndarray]

class evo_lmm.ConvergenceReport(status='unknown', converged=False, iterations=0, step_se_norm=nan, step_se_tol=nan, newton_decrement=nan, score_norm=nan, accepted_step=nan, ai_condition=nan, ai_damping=nan)#

Bases: object

How a fit ended, in sample-size-independent terms.

Both fitters – single-component evo_lmm.fit_reml() and evo_lmm.fit_multicomponent_reml() – judge convergence by the same rule and report it through this object, so the two are comparable.

step_se_norm is the criterion: max_i |step_i| / SE_i with step = AI^-1 score and SE = sqrt(diag(AI^-1)) – the largest move the next undamped Newton step would make in any coordinate, in units of that coordinate’s own standard error. Convergence is declared when it is at or below step_se_tol. newton_decrement is sqrt(score' AI^-1 score), logged alongside it. Both are dimensionless, invariant to a reparameterization of the coordinates, and order one at a fixed distance from the optimum regardless of sample size.

score_norm is a diagnostic, not a criterion: at a fixed statistical distance from the optimum it grows roughly like sqrt(n), so no absolute bound on it means the same thing at two sample sizes. It is reported because it is cheap and occasionally informative, never to be thresholded.

status names the exit, and is the field to branch on:

converged

step_se_norm <= step_se_tol at an accepted iterate.

converged_after_dense_finish

The criterion was met only after the exact dense finishing optimizer ran (single-component fitter).

stalled_near_tolerance

Every step halving was rejected, but the criterion was within ten times its tolerance, so the iterate is accepted as converged.

line_search_stalled

The iteration budget ran out and the final iteration’s step was rejected.

iteration_cap

The iteration budget ran out with steps still being accepted.

optimizer_stalled

A delegated optimizer (dense L-BFGS-B) returned without meeting the criterion, whatever its own verdict was.

unidentified

The criterion is met on every direction the information matrix locates, but at least one coordinate has a standard error wider than its own coordinate box, so the data do not place it anywhere. The fit is done; that coordinate’s value is not an estimate and must not be reported. A pure-null panel lands here: with no genetic variance there is no information about the frequency-shape parameter.

not_started

max_iter=0; the reported state is the initial point.

oracle

Not produced by an optimizer at all – an explicitly constructed covariance, used by oracle tests.

converged is the summary status in CONVERGED_STATUSES; it never distinguishes these cases on its own, which is why status exists.

Parameters:

  • status (str)

  • converged (bool)

  • iterations (int)

  • step_se_norm (float)

  • step_se_tol (float)

  • newton_decrement (float)

  • score_norm (float)

  • accepted_step (float)

  • ai_condition (float)

  • ai_damping (float)

class evo_lmm.DenseREMLOracle(genotypes, frequencies, covariates=None, *, model='simplified')#

Bases: object

Convenient exact REML oracle for tests and small dense simulations.

Parameters:

  • genotypes (np.ndarray)

  • frequencies (np.ndarray)

  • covariates (np.ndarray | None)

  • model (str)

evo_lmm.convergence_statistics(score, ai, *, se_limit=60.0)#

Return the scale-free convergence statistics (step_se_norm, decrement).

step_se_norm = max_i |step_i| / SE_i with step = AI^-1 score and SE_i = sqrt((AI^-1)_ii): the largest move the next undamped Newton step would make in any coordinate, measured in that coordinate’s own standard error. decrement = sqrt(score' AI^-1 score) is the Newton decrement, which approximates twice the objective decrease the step would buy.

Both are dimensionless and invariant to a reparameterization of the coordinates. They also measure distance from the optimum in standard-error units, so they stay order one at a fixed statistical distance regardless of n, whereas ||score||_inf at the same point grows roughly like sqrt(n). A fixed absolute score tolerance therefore demands getting sqrt(n) times closer as data are added, which is why it is not the gate.

The statistic is taken over the coordinates the information matrix locates; unidentified_coordinates() masks the rest, and a fit carrying any of them reports status="unidentified" so the value is not read as an estimate. Excluding them is not leniency: where the standard error is 1e10 log units, score-over-information is numerical noise divided by numerical noise, and thresholding it reports the floating-point floor rather than the fit.

When no coordinate is identified the result is inf, not zero: an iteration must never stop because its information matrix was momentarily unusable. That case is named by the unidentified status once the fit is over, not treated as convergence mid-loop. Non-finite input also returns inf.

Parameters:

  • score (ndarray)

  • ai (ndarray)

  • se_limit (float)

Return type:

tuple[float, float]

class evo_lmm.EvolutionaryLmmOps(chromosomes, frequencies=None, covariates=None, *, sample_filter=None, model='simplified', mutation_filter=None)#

Bases: object

Matrix-free projected raw-dosage operators for evolutionary LMMs.

Parameters:

  • chromosomes (Any) – An (N, M) dense dosage matrix, a mapping of chromosome labels to matrices/GRGs, or a sequence of matrices/(label, source) pairs.

  • frequencies (Any) – Sample allele frequencies, supplied as one vector, one vector per chromosome, or a mapping. GRG inputs may omit this and extract them through GRAPP’s frequency traversal.

  • covariates (np.ndarray | None) – Fixed-effect design. An intercept is added when absent and the basis is QR-orthonormalised once.

  • model (str) – Coordinate interpretation for transformed phi arrays.

  • sample_filter (Sequence[int] | None)

  • mutation_filter (Mapping[Any, Sequence[int]] | None)

classmethod from_dense(genotypes, frequencies, covariates=None, *, chrom_labels=None, model='simplified')#

Convenience constructor for dense tests and small simulations.

Parameters:

  • genotypes (ndarray | Sequence[ndarray])

  • frequencies (ndarray | Sequence[ndarray])

  • covariates (ndarray | None)

  • chrom_labels (Sequence[Any] | None)

  • model (str)

Return type:

EvolutionaryLmmOps

test_stats(chrom)#

Return prior-independent tested-genotype statistics for a chromosome.

These are the quantities the BOLT-style calibration and association formulas need: the projected raw-dosage norms, the BOLT normalisation, and the eligibility mask. They depend only on the genotypes and the covariate basis, never on the fitted evolutionary prior.

Parameters:

chrom (Any)

Return type:

TestVariantStats

apply_model_x(weights, theta=None, exclude_chrom=None)#

Apply P_C X diag(weights) to variant coefficients.

weights is normally a vector of model coefficients. When theta is a prior, it applies the model operator P_C X diag(sqrt(w(theta))) to those coefficients. A raw coefficient vector or chromosome mapping may also be supplied without theta. This method deliberately contains no sample-standardisation or 1/M normalization.

Parameters:

  • weights (ndarray | Mapping[Any, ndarray] | EvolutionaryPrior)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

model_scores(vector, theta=None, exclude_chrom=None)#

Return model scores, optionally for B_theta^T = diag(sqrt(w))X^T P_C.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_k(vector, theta, exclude_chrom=None)#

Apply K = P_C X diag(w(theta)) X^T P_C.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dk(vector, theta, parameter, exclude_chrom=None)#

Apply a first derivative of K in a transformed shape coordinate.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_h(vector, phi, exclude_chrom=None)#

Apply the projected shape matrix H = K + delta I.

Parameters:

  • vector (ndarray)

  • phi (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dh(vector, phi, parameter, exclude_chrom=None)#

Apply a first derivative of H in log_delta, log_tau, or logit_r.

Parameters:

  • vector (ndarray)

  • phi (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dh_matmat(values, phi, parameter, exclude_chrom=None)#

Apply a first derivative of H to several columns at once.

This is the matrix-RHS counterpart of apply_dh(). Keeping the columns batched is important for GRG inputs: one matmat traversal replaces one traversal per stochastic trace probe.

Parameters:

  • values (ndarray)

  • phi (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_k_matmat(values, theta, exclude_chrom=None)#

Apply K to batched right-hand sides without per-column traversals.

Parameters:

  • values (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

solve_ph(rhs_columns, phi, exclude_chrom=None, *, tol=0.0005, max_iter=None, initial=None, stats=None)#

Solve projected H z = P_C rhs for one or many right-hand sides.

initial is an optional warm start. Each column is validated against a zero start independently; a poor or invalid cached column is reset without affecting the other right-hand sides. If stats is passed, it is populated with aggregate iteration and warm-start diagnostics.

Parameters:

  • rhs_columns (ndarray)

  • phi (Any)

  • exclude_chrom (Any)

  • tol (float)

  • max_iter (int | None)

  • initial (ndarray | None)

  • stats (dict[str, Any] | None)

Return type:

ndarray

test_scores(chrom, vector)#

Return BOLT-normalised test-genotype scores for one chromosome.

Parameters:

  • chrom (Any)

  • vector (ndarray)

Return type:

ndarray

test_column(chrom, local_idx)#

Return one projected, BOLT-normalised test-genotype column.

Parameters:

  • chrom (Any)

  • local_idx (int)

Return type:

ndarray

chromosome_frequencies(chrom)#

Return sample allele frequencies for one chromosome in operator order.

Parameters:

chrom (Any)

Return type:

ndarray

local_indices(chrom)#

Return mutation identifiers for a chromosome in operator order.

Parameters:

chrom (Any)

Return type:

ndarray

kernel_trace(theta, exclude_chrom=None)#

Return tr(P_C K P_C) without constructing a dense kernel.

Parameters:

  • theta (Any)

  • exclude_chrom (Any)

Return type:

float

dense_kernel(theta, exclude_chrom=None)#

Materialise a kernel for a small dense input or a test oracle.

Parameters:

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

class evo_lmm.EvolutionaryPrior#

Bases: object

Common interface implemented by the two evolutionary priors.

coordinates(delta)#

Return unconstrained (log_delta, log_tau[, logit_r]) coordinates.

Parameters:

delta (float)

Return type:

ndarray

class evo_lmm.FitDiagnostics(convergence, trace_estimator, trace_probes, objective=nan, initialization='default', trace_operator_queries=0, trace_standard_errors=<factory>, cg_iterations=<factory>, cg_warm_start_hits=0, cg_warm_start_rejections=0, cg_initial_residual_norms=<factory>, cg_final_residual_norms=<factory>, random_seed=None, boundary_hits=(), warnings=())#

Bases: object

Numerical and identifiability diagnostics from a fit.

Convergence lives in ConvergenceReport under convergence; the flat accessors below delegate to it so both fitters expose the same names.

objective is the profiled REML objective or ``nan``. It is never a stand-in for something else: the stochastic paths do not evaluate a log-determinant, so they report nan rather than a surrogate, and convergence is judged by convergence.step_se_norm in every path.

Parameters:

  • convergence (ConvergenceReport)

  • trace_estimator (str)

  • trace_probes (int)

  • objective (float)

  • initialization (str)

  • trace_operator_queries (int)

  • trace_standard_errors (dict[str, float])

  • cg_iterations (list[int])

  • cg_warm_start_hits (int)

  • cg_warm_start_rejections (int)

  • cg_initial_residual_norms (list[float])

  • cg_final_residual_norms (list[float])

  • random_seed (int | None)

  • boundary_hits (tuple[str, ...])

  • warnings (tuple[str, ...])

class evo_lmm.FitResult(prior, sigma_b2, sigma_e2, delta, h2, log_likelihood, fixed_effects, projected_phenotype, ph_y, diagnostics, model, ops=None)#

Bases: object

Scientific-scale fit result for an evolutionary LMM.

Parameters:

  • prior (EvolutionaryPrior)

  • sigma_b2 (float)

  • sigma_e2 (float)

  • delta (float)

  • h2 (float)

  • log_likelihood (float)

  • fixed_effects (ndarray)

  • projected_phenotype (ndarray)

  • ph_y (ndarray)

  • diagnostics (FitDiagnostics)

  • model (str)

  • ops (Any)

property sigma_g2: float#

Compatibility alias; unlike GRAPP this is raw-effect sigma_b2.

blup()#

Return the projected genetic-value BLUP sigma_b2 K P_V y.

Return type:

ndarray

class evo_lmm.FullPrior(sigma_b2, tau, rho)#

Bases: EvolutionaryPrior

The full evolutionary prior with identifiable coupling rho^2.

Only rho^2 appears in the conditional variance, so the sign of rho is not identifiable. The public parameter accepts -1 <= rho <= 1 and the optimizer works directly with r = rho^2.

Parameters:

  • sigma_b2 (float)

  • tau (float)

  • rho (float)

class evo_lmm.SimplifiedPrior(sigma_b2, tau)#

Bases: EvolutionaryPrior

The exact rho^2 = 1 evolutionary prior.

Parameters:

  • sigma_b2 (float) – Per-locus focal-trait effect variance on raw dosage units.

  • tau (float) – Non-negative sigma_a^2 / W_S aggregate. tau=0 is the frequency-independent boundary and is useful for tests and diagnostics.

class evo_lmm.GrgSimulation(tree_sequence, grg, frequencies, effects, genetic_value, phenotype, prior, residual_variance, seed)#

Bases: object

A simulated GRG phenotype with effects sampled from an evolutionary prior.

Parameters:

  • tree_sequence (Any)

  • grg (Any)

  • frequencies (ndarray)

  • effects (ndarray)

  • genetic_value (ndarray)

  • phenotype (ndarray)

  • prior (EvolutionaryPrior)

  • residual_variance (float)

  • seed (int)

class evo_lmm.VariantData(frequencies, local_idx, raw_norm2=None, centered_norm2=None)#

Bases: object

Cached per-variant information used by an evolutionary operator.

frequencies are sample allele frequencies and local_idx are the mutation identifiers in the source chromosome. raw_norm2 and centered_norm2 are optional because GRGL operators can obtain them via graph traversals when needed.

Parameters:

  • frequencies (ndarray)

  • local_idx (ndarray)

  • raw_norm2 (ndarray | None)

  • centered_norm2 (ndarray | None)

evo_lmm.allele_frequency_q(frequencies)#

Return q = x_hat * (1 - x_hat) for sample allele frequencies.

Parameters:

frequencies (ndarray)

Return type:

ndarray

evo_lmm.association(result, y=None, *, use_loco=True, calibrate=True, calibration=None, calibration_variants=30, seed=0, screen_threshold=5.0, cg_tol=0.0005)#

Compute calibrated BOLT-style association statistics per chromosome.

The mixed-model statistic uses the fitted evolutionary LOCO covariance V_loco = sigma_b2 * (K_evo,loco + delta I) while the tested columns keep GRAPP’s independent BOLT normalisation. With calibrate=True the prospective/retrospective moment matching of evo_lmm.calibrate_association() supplies the per-chromosome inverse scale; otherwise the uncalibrated (factor = 1) scale is used, which is only appropriate for diagnostics.

beta and se are returned in raw diploid-dosage units. A single-variant linear-regression chi-square is reported alongside the mixed model statistic so inflation can be compared directly.

Parameters:

  • result (FitResult)

  • y (ndarray | None)

  • use_loco (bool)

  • calibrate (bool)

  • calibration (CalibrationResult | None)

  • calibration_variants (int)

  • seed (int)

  • screen_threshold (float)

  • cg_tol (float)

Return type:

list[AssociationResult]

evo_lmm.association_summary(results)#

Return mean chi-square and lambda_GC for the mixed model and linreg.

lambda_gc divides the median chi-square by the median of a one-degree-of-freedom chi-square, so a calibrated null gives values near one for both statistics.

Parameters:

results (Sequence[AssociationResult])

Return type:

dict[str, dict[str, float]]

evo_lmm.calibrate_association(result, y=None, *, count=30, seed=0, screen_threshold=5.0, cg_tol=0.0005, stats=None)#

Calibrate the LMM statistic against the fitted evolutionary covariance.

For every selected variant j on its own chromosome c the retrospective and prospective statistics are

retro = (N - C) * (x_j^T u_c)^2 / (||u_c||^2 * ||x_j||^2)

pro   = (N - C) * (x_j^T u_c)^2 / (x_j^T H_c^-1 x_j) / (y^T u_c)

with u_c = H_c^-1 P_C y and H_c the LOCO evolutionary shape matrix. Both are invariant to sigma_b2, so the factor is a pure shape correction; sigma_b2 enters only through inverse_scale.

Parameters:

  • result (FitResult)

  • y (ndarray | None)

  • count (int)

  • seed (int)

  • screen_threshold (float)

  • cg_tol (float)

  • stats (dict[Any, TestVariantStats] | None)

Return type:

CalibrationResult

evo_lmm.select_calibration_variants(result, *, count=30, seed=0, screen_threshold=5.0, stats=None)#

Select calibration variants with BOLT’s blocked GRAMMAR pre-screen.

Eligible variants are split into count equal blocks in chromosome order. One variant is drawn uniformly from each block and accepted when its all-chromosome GRAMMAR retrospective statistic falls below screen_threshold, which keeps strongly associated variants out of the calibration set. Returns the selected (chrom, local_idx) pairs and the number of candidates examined.

Parameters:
Return type:

tuple[list[tuple[Any, int]], int]

class evo_lmm.TestVariantStats(chrom, local_idx, centered_norm2, projected_norm2, test_scale, norm_scale, std_projected_norm2, model_mask)#

Bases: object

Per-variant tested-genotype statistics for one chromosome block.

All quantities are computed after the covariate projection P_C. The test_scale is the BOLT normalisation sqrt(centered_norm2 / (n - 1)) applied to raw diploid dosage columns, so norm_scale = 1 / test_scale converts a raw-dosage quantity into BOLT-normalised units. Nothing here depends on the fitted evolutionary prior.

Parameters:

  • chrom (Any)

  • local_idx (ndarray)

  • centered_norm2 (ndarray)

  • projected_norm2 (ndarray)

  • test_scale (ndarray)

  • norm_scale (ndarray)

  • std_projected_norm2 (ndarray)

  • model_mask (ndarray)

projected_norm2#

||P_C x_j||^2 in raw diploid dosage units.

Type:

numpy.ndarray

std_projected_norm2#

||P_C x_j||^2 in BOLT-normalised units, i.e. projected_norm2 * norm_scale**2.

Type:

numpy.ndarray

model_mask#

Variants eligible for association output; monomorphic and covariate-collinear columns are excluded exactly as in GRAPP.

Type:

numpy.ndarray

evo_lmm.dense_variant_data(genotypes, frequencies=None, *, local_idx=None)#

Build variant metadata from an (N, M) raw dosage matrix.

If frequencies are omitted, they are calculated as the mean dosage divided by ploidy two. This default is deliberately explicit in the docstring and is only intended for diploid dense test/simulation inputs.

Parameters:

  • genotypes (ndarray)

  • frequencies (ndarray | None)

  • local_idx (Sequence[int] | None)

Return type:

VariantData

evo_lmm.exact_average_information(y, H, derivatives, covariates=None)#

Return the direct quadratic-form average-information matrix.

Parameters:

  • y (ndarray)

  • H (ndarray)

  • derivatives (Sequence[ndarray])

  • covariates (ndarray | None)

Return type:

ndarray

evo_lmm.exact_reml_score(y, H, derivatives, covariates=None, *, scale=1.0)#

Return exact shape scores for V = scale * H.

The scale is normally the profiled sigma_b2. Omitting it evaluates the equivalent unit-scale score.

Parameters:

  • y (ndarray)

  • H (ndarray)

  • derivatives (Sequence[ndarray])

  • covariates (ndarray | None)

  • scale (float)

Return type:

ndarray

evo_lmm.fit_dense_reml(genotypes, frequencies, y, *, covariates=None, model='simplified', **kwargs)#

Fit the exact dense oracle in one call.

Parameters:

  • genotypes (ndarray)

  • frequencies (ndarray)

  • y (ndarray)

  • covariates (ndarray | None)

  • model (str)

  • kwargs (Any)

Return type:

FitResult

evo_lmm.fit_evolutionary_bolt_lmm(chrom_grgs, y, *, frequencies=None, covariates=None, model='simplified', initial=None, sample_filter=None, trace_probes=12, seed=0, max_iter=50, cg_tol=0.0005, trace_method='hutchinson', warm_start=True, initialization='default')#

Fit the CPU GRGL-backed evolutionary model across chromosomes.

Parameters:

  • chrom_grgs (Sequence[tuple[Any, Any]])

  • y (ndarray)

  • frequencies (Mapping[Any, ndarray] | Sequence[ndarray] | None)

  • covariates (ndarray | None)

  • model (str)

  • initial (Any)

  • sample_filter (Sequence[int] | None)

  • trace_probes (int)

  • seed (int)

  • max_iter (int)

  • cg_tol (float)

  • trace_method (str)

  • warm_start (bool)

  • initialization (str)

Return type:

FitResult

evo_lmm.fit_evolutionary_bolt(chrom_grgs, y, *, frequencies=None, covariates=None, model='simplified', initial=None, sample_filter=None, trace_probes=12, seed=0, max_iter=50, cg_tol=0.0005, trace_method='hutchinson', warm_start=True, initialization='default')#

Fit the CPU GRGL-backed evolutionary model across chromosomes.

Parameters:

  • chrom_grgs (Sequence[tuple[Any, Any]])

  • y (ndarray)

  • frequencies (Mapping[Any, ndarray] | Sequence[ndarray] | None)

  • covariates (ndarray | None)

  • model (str)

  • initial (Any)

  • sample_filter (Sequence[int] | None)

  • trace_probes (int)

  • seed (int)

  • max_iter (int)

  • cg_tol (float)

  • trace_method (str)

  • warm_start (bool)

  • initialization (str)

Return type:

FitResult

evo_lmm.fit_evolutionary_reml(ops, y, *, model=None, initial=None, delta=1.0, trace_probes=12, seed=0, max_iter=50, tol=1e-06, step_se_tol=0.01, cg_tol=0.0005, max_step=2.0, exact=None, trace_method='hutchinson', warm_start=True, initialization='default')#

Fit evolutionary shape parameters by profiled average-information REML.

sigma_b2 is profiled as y'P_H y / (N-rank(C)) and sigma_e2 is derived as delta*sigma_b2. Dense operators use exact traces by default; matrix-free operators use fixed Rademacher probes (Hutchinson) and warm-started projected CG solves.

The stochastic defaults are deliberately cheap and match GRAPP’s solver budget: trace_probes=12 and cg_tol=5e-4. They are appropriate for exploratory fits and for benchmark parity, not for final reported estimates. If a larger stochastic trace budget is desired, raise trace_probes while retaining the sketch solver tolerance cg_tol=5e-4; inspect FitDiagnostics.trace_standard_errors as a diagnostic rather than treating it as a separate convergence rule.

Convergence is declared when step_se_tol bounds max_i |step_i| / SE_i – see convergence_statistics(). The default 1e-2 means the next undamped Newton step would move no coordinate by more than one percent of its own standard error. tol keeps its two remaining roles: the line search accepts a trial displacement below it, and it is passed as gtol to the exact dense finishing optimizer. It is no longer the convergence gate, because ||score||_inf grows with n and a fixed absolute bound on it is a different requirement at every sample size. FitDiagnostics.status names how the fit ended; converged alone does not distinguish the criterion being met from the finishing optimizer’s loose back-stop.

Parameters:

  • ops (EvolutionaryLmmOps)

  • y (ndarray)

  • model (str | None)

  • initial (Any)

  • delta (float)

  • trace_probes (int)

  • seed (int)

  • max_iter (int)

  • tol (float)

  • step_se_tol (float)

  • cg_tol (float)

  • max_step (float)

  • exact (bool | None)

  • trace_method (str)

  • warm_start (bool)

  • initialization (str)

Return type:

FitResult

evo_lmm.fit_evolutionary_lmm(genotypes, y, frequencies=None, *, covariates=None, model='simplified', initial=None, sample_filter=None, trace_probes=12, seed=0, max_iter=50, cg_tol=0.0005, exact=None, trace_method='hutchinson', warm_start=True, initialization='default')#

Fit a simplified or full evolutionary LMM from dense matrices or GRGs.

Missing phenotypes are removed before constructing the operator. For GRGs, the corresponding individual sample filter is passed to GRAPP and allele frequencies are recomputed on the retained sample.

Parameters:

  • genotypes (Any)

  • y (ndarray)

  • frequencies (Any)

  • covariates (ndarray | None)

  • model (str)

  • initial (Any)

  • sample_filter (Sequence[int] | None)

  • trace_probes (int)

  • seed (int)

  • max_iter (int)

  • cg_tol (float)

  • exact (bool | None)

  • trace_method (str)

  • warm_start (bool)

  • initialization (str)

Return type:

FitResult

evo_lmm.fit_reml(ops, y, *, model=None, initial=None, delta=1.0, trace_probes=12, seed=0, max_iter=50, tol=1e-06, step_se_tol=0.01, cg_tol=0.0005, max_step=2.0, exact=None, trace_method='hutchinson', warm_start=True, initialization='default')#

Fit evolutionary shape parameters by profiled average-information REML.

sigma_b2 is profiled as y'P_H y / (N-rank(C)) and sigma_e2 is derived as delta*sigma_b2. Dense operators use exact traces by default; matrix-free operators use fixed Rademacher probes (Hutchinson) and warm-started projected CG solves.

The stochastic defaults are deliberately cheap and match GRAPP’s solver budget: trace_probes=12 and cg_tol=5e-4. They are appropriate for exploratory fits and for benchmark parity, not for final reported estimates. If a larger stochastic trace budget is desired, raise trace_probes while retaining the sketch solver tolerance cg_tol=5e-4; inspect FitDiagnostics.trace_standard_errors as a diagnostic rather than treating it as a separate convergence rule.

Convergence is declared when step_se_tol bounds max_i |step_i| / SE_i – see convergence_statistics(). The default 1e-2 means the next undamped Newton step would move no coordinate by more than one percent of its own standard error. tol keeps its two remaining roles: the line search accepts a trial displacement below it, and it is passed as gtol to the exact dense finishing optimizer. It is no longer the convergence gate, because ||score||_inf grows with n and a fixed absolute bound on it is a different requirement at every sample size. FitDiagnostics.status names how the fit ended; converged alone does not distinguish the criterion being met from the finishing optimizer’s loose back-stop.

Parameters:

  • ops (EvolutionaryLmmOps)

  • y (ndarray)

  • model (str | None)

  • initial (Any)

  • delta (float)

  • trace_probes (int)

  • seed (int)

  • max_iter (int)

  • tol (float)

  • step_se_tol (float)

  • cg_tol (float)

  • max_step (float)

  • exact (bool | None)

  • trace_method (str)

  • warm_start (bool)

  • initialization (str)

Return type:

FitResult

evo_lmm.exact_reml_loglikelihood(y, covariance, covariates=None, *, scale=1.0, include_constant=False)#

Evaluate the exact restricted log likelihood for a dense covariance.

covariance is the full V matrix. The profiled form used by the fitter is obtained by passing scale=1 and profiling that scale outside this function. A fixed intercept is used when covariates are omitted.

Parameters:

  • y (ndarray)

  • covariance (ndarray)

  • covariates (ndarray | None)

  • scale (float)

  • include_constant (bool)

Return type:

float

evo_lmm.profiled_average_information(y, H, derivatives, covariates=None, *, scale=1.0)#

Return the shape AI after eliminating the profiled scale coordinate.

The Schur complement is formed from the full coordinate vector (log_sigma_b2, shape...) with H_scale = H. scale is the current profiled sigma_b2.

Parameters:

  • y (ndarray)

  • H (ndarray)

  • derivatives (Sequence[ndarray])

  • covariates (ndarray | None)

  • scale (float)

Return type:

ndarray

evo_lmm.haseman_elston_initialization(ops, y, prior, *, probes=12, seed=0)#

Estimate (sigma_b2, sigma_e2, delta) by projected HE moments.

The moment solution is intentionally only an initializer. Exact dense traces are used for dense operators; spherical XTrace estimates are used for GRG operators.

Parameters:
Return type:

tuple[float, float, float]

evo_lmm.grg_variant_data(grg, *, frequencies=None, sample_filter=None, mutation_filter=None)#

Build cached metadata for one GRG chromosome.

GRGL retains mutation identifiers even when a mutation filter is supplied; the returned local_idx preserves those identifiers for test_column.

Parameters:

  • grg (Any)

  • frequencies (ndarray | None)

  • sample_filter (Sequence[int] | None)

  • mutation_filter (Sequence[int] | None)

Return type:

VariantData

evo_lmm.loco_solve(result, y, exclude_chrom)#

Solve the fitted projected LOCO shape system for one chromosome.

Parameters:

  • result (FitResult)

  • y (ndarray)

  • exclude_chrom (Any)

Return type:

ndarray

evo_lmm.loco_solutions(result, y=None)#

Return one projected P_H y solution per chromosome exclusion.

Parameters:
Return type:

dict[Any, ndarray]

evo_lmm.predict_blup(result)#

Return projected genetic-value BLUPs.

Parameters:

result (FitResult)

Return type:

ndarray

evo_lmm.prior_from_coordinates(model, coordinates, *, sigma_b2=1.0)#

Construct a unit/profiled prior and delta from optimizer coordinates.

The returned prior has the supplied sigma_b2 and coordinates are (log_delta, log_tau[, logit_r]). Exact boundary values are supported by allowing infinite transformed coordinates.

Parameters:

  • model (str)

  • coordinates (ndarray)

  • sigma_b2 (float)

Return type:

tuple[EvolutionaryPrior, float]

evo_lmm.prior_from_parameters(model, *, sigma_b2, tau, rho=1.0)#

Create a prior from scientific parameters with an explicit model name.

Parameters:

  • model (str)

  • sigma_b2 (float)

  • tau (float)

  • rho (float)

Return type:

EvolutionaryPrior

evo_lmm.restricted_log_likelihood(y, covariance, covariates=None, *, scale=1.0, include_constant=False)#

Evaluate the exact restricted log likelihood for a dense covariance.

covariance is the full V matrix. The profiled form used by the fitter is obtained by passing scale=1 and profiling that scale outside this function. A fixed intercept is used when covariates are omitted.

Parameters:

  • y (ndarray)

  • covariance (ndarray)

  • covariates (ndarray | None)

  • scale (float)

  • include_constant (bool)

Return type:

float

evo_lmm.sample_allele_frequencies(grg, *, sample_filter=None, adjust_missing=True)#

Extract sample allele frequencies from a GRGL-backed chromosome.

sample_filter is expressed in individual indices, matching GRAPP’s raw genotype operators. Missing alleles are excluded from the denominator by default, so the result remains a sample frequency rather than an imputed dosage frequency.

Parameters:

  • grg (Any)

  • sample_filter (Sequence[int] | None)

  • adjust_missing (bool)

Return type:

ndarray

evo_lmm.simulate_grg_lmm(prior, *, n_individuals=24, sequence_length=100000.0, population_size=1000.0, recombination_rate=1e-08, mutation_rate=2e-07, residual_variance=0.4, seed=7)#

Simulate a diploid tree sequence, convert it to a GRG, and draw a trait.

The returned effects are raw-dosage SNP effects sampled independently as Normal(0, prior.effect_variances(frequencies)). The genetic value is computed by the GRG-backed raw operator and projected off the intercept, matching the fitting boundary used by fit_evolutionary_bolt_lmm().

The temporary .trees file is removed after GRGL conversion. The GRG object and original in-memory tree sequence remain available in the result.

Parameters:

  • prior (EvolutionaryPrior)

  • n_individuals (int)

  • sequence_length (float)

  • population_size (float)

  • recombination_rate (float)

  • mutation_rate (float)

  • residual_variance (float)

  • seed (int)

Return type:

GrgSimulation

class evo_lmm.MultiComponentFit(prior, sigma_e2, h2, convergence, ops, objective=nan, ai_covariance=None, standard_errors=None, trace_method='hutchinson', trace_probes=0, cg_tol=nan, trace_standard_error=nan, phenotype=None, warnings=(), initialization='default', mom_raw_component_scales=None, mom_truncated=None)#

Bases: object

Profiled REML result for the partitioned simplified model.

Convergence is reported through evo_lmm.ConvergenceReport, the same object the single-component fitter uses, so the two fitters are judged by one criterion and read the same way. The flat accessors below delegate to it.

objective is the profiled REML objective or ``nan``. The AI path never evaluates a log-determinant, so it reports nan rather than a surrogate; an earlier revision put 0.5*||score||^2 in this slot, which was neither a likelihood nor the convergence criterion. Only the exact dense method returns a real objective here.

warnings reports each tau_c in a regime where it is unidentified. It is a report only: a boundary hit does not change status or converged.

Parameters:

  • prior (MultiComponentPrior)

  • sigma_e2 (float)

  • h2 (float)

  • convergence (ConvergenceReport)

  • ops (MultiComponentOps)

  • objective (float)

  • ai_covariance (ndarray | None)

  • standard_errors (dict[str, float] | None)

  • trace_method (str)

  • trace_probes (int)

  • cg_tol (float)

  • trace_standard_error (float)

  • phenotype (ndarray | None)

  • warnings (tuple[str, ...])

  • initialization (str)

  • mom_raw_component_scales (ndarray | None)

  • mom_truncated (ndarray | None)

class evo_lmm.MultiComponentOps(components)#

Bases: object

Sum of independently partitioned projected raw-dosage component operators.

Parameters:

components (Mapping[Any, EvolutionaryLmmOps])

project(values)#

Project vectors or matrices off the shared covariate basis.

Parameters:

values (ndarray)

Return type:

ndarray

classmethod from_operators(components)#

Construct a partition from already-adapted dense or GRGL operators.

Parameters:

components (Mapping[Any, EvolutionaryLmmOps])

Return type:

MultiComponentOps

kernel_trace(prior)#

Return tr(P_C K P_C) without materialising a dense kernel.

Parameters:

prior (MultiComponentPrior)

Return type:

float

apply_k(values, prior)#

Apply the summed partitioned kernel to one or many vectors.

Parameters:
Return type:

ndarray

derivative_kernels(prior)#

Return analytic derivatives in the same coordinate order as the prior.

Parameters:

prior (MultiComponentPrior)

Return type:

dict[str, ndarray]

apply_dh_matmat(values, prior, parameter)#

Apply a component derivative to batched vectors (dense/GRGL path).

Parameters:
Return type:

ndarray

apply_component_derivatives_matmat(values, prior)#

Apply all component derivatives to shared batched right-hand sides.

The projection and right-hand-side batch are shared across the returned derivatives; GRGL-backed component operators retain their matmat traversal rather than falling back to one traversal per probe.

Parameters:
Return type:

dict[str, ndarray]

apply_shape_matmat(values, prior)#

Apply H = I + sum_c K_c to batched projected vectors.

Parameters:
Return type:

ndarray

class evo_lmm.MultiComponentPrior(labels, components)#

Bases: object

One simplified prior per annotation category.

coordinates are ordered (log_sigma_b2_c, log_tau_c) pairs. A tau_c of zero is represented by -inf in transformed coordinates.

Parameters:
classmethod flat(labels, scales=None)#

Return the exact M0 boundary with all tau_c = 0.

Parameters:

  • labels (Sequence[Any])

  • scales (Sequence[float] | None)

Return type:

MultiComponentPrior

with_shared_tau(tau)#

Return M1’s shared-tau identifiability-crutch specification.

Parameters:

tau (float)

Return type:

MultiComponentPrior

evo_lmm.fit_multicomponent_reml(ops, y, *, initial=None, max_iter=200, method='ai', trace_method='hutchinson', trace_probes=12, seed=0, cg_tol=0.0005, tol=1e-06, step_se_tol=0.01, max_step=2.0, fit_tau=True, initialization='default')#

Fit by profiled REML.

The residual scale sigma_e2 is profiled; component scales are searched as ratios to that scale in log coordinates. The objective is REML, and the reported component sigma_b2 values are returned on the scientific scale after profiling. tol is the score-norm convergence tolerance; the default matches the single-component fitter. It is no longer the convergence gate: convergence is declared when step_se_tol bounds max_i |step_i| / SE_i, the same scale-free statistic the single-component fitter uses, and tol only accepts a vanishing line-search displacement.

fit_tau=False holds every tau_c at its value in initial and searches only the |c| scale coordinates. That is the pooled-shape mode used for per-gene fitting, where the shapes are estimated once across genes and must not be re-estimated per gene.

initialization="he" applies the joint projected Haseman–Elston (|c|+1) moment system to the component scale coordinates. The current category-specific tau_c values are retained: HE is linear in the covariance components conditional on those evolutionary weights and cannot identify the nonlinear shape coordinates by itself. Raw moment scales and negative-scale flags are retained on the returned fit for the required no-truncation audit; invalid or non-positive values fall back to the requested/default scale for optimisation.

Parameters:

  • ops (MultiComponentOps)

  • y (ndarray)

  • initial (MultiComponentPrior | ndarray | None)

  • max_iter (int)

  • method (str)

  • trace_method (str)

  • trace_probes (int)

  • seed (int)

  • cg_tol (float)

  • tol (float)

  • step_se_tol (float)

  • max_step (float)

  • fit_tau (bool)

  • initialization (str)

Return type:

MultiComponentFit

evo_lmm.profiled_reml_objective(ops, y, prior)#

Evaluate the exact dense profiled-REML objective for a fixed prior.

Returns (objective, sigma_e2) for V = sigma_e2 * (I + K). This is the small-dense reference used to verify the M0/M1/M2 nesting ladder.

Parameters:
Return type:

tuple[float, float]

class evo_lmm.CollapsedVariants(genotypes, frequencies, source_indices)#

Bases: object

Result of optional MAC-threshold burden construction.

Parameters:

  • genotypes (dict[Any, ndarray])

  • frequencies (dict[Any, ndarray])

  • source_indices (dict[Any, tuple[tuple[int, ...], ...]])

class evo_lmm.MoMResult(component_scales, residual_variance, raw_component_scales, truncated, system, trace_standard_errors=None)#

Bases: object

Joint method-of-moments estimates and the negative-estimate audit.

Parameters:

  • component_scales (ndarray)

  • residual_variance (float)

  • raw_component_scales (ndarray)

  • truncated (ndarray)

  • system (ndarray)

  • trace_standard_errors (dict[str, float] | None)

evo_lmm.collapse_mac(genotypes, *, ploidy=2, mac_threshold=10.0)#

Collapse variants with MAC below mac_threshold per category.

Filtering and collapsing occur before frequency recomputation. The collapsed column is a dosage burden and its frequency is its sample mean divided by ploidy * n.

Parameters:

  • genotypes (Mapping[Any, ndarray])

  • ploidy (int)

  • mac_threshold (float)

Return type:

CollapsedVariants

evo_lmm.flat_prior(labels, scales=None)#

Return named M0: a flat per-category prior with tau_c = 0.

Parameters:

  • labels (tuple[Any, ...])

  • scales (Mapping[Any, float] | None)

Return type:

MultiComponentPrior

evo_lmm.joint_mom_initialization(ops, y, prior, *, trace_method='exact', trace_probes=12, seed=0)#

Solve the (|c|+1) projected Haseman–Elston moment system.

The returned raw_component_scales are never truncated. The component_scales field applies the RareEffect boundary rule for an explicitly requested baseline comparison.

Parameters:
Return type:

MoMResult

evo_lmm.fit_rare_effect_baseline(ops, y, *, max_iter=100)#

Fit the named flat baseline marginally and apply the MoM-ratio rule.

Each category is fitted independently by exact restricted (REML) profiling at the tau=0 boundary. The joint adjustment is then computed from the partitioned moment system; no evolutionary weighting is introduced.

Parameters:
Return type:

RareEffectBaselineResult

class evo_lmm.RareEffectBaselineResult(marginal_scales, marginal_mom_scales, joint_mom_scales, adjusted_scales, negative_mom_fallback)#

Bases: object

Marginal per-category estimates plus the MoM-ratio adjustment.

marginal_scales are restricted (REML) marginal estimates: the objective in fit_rare_effect_baseline() carries the slogdet(B' V^-1 B) term. Whether RareEffect’s published pipeline uses ML or REML here is not verified in this repository, so do not describe this field as reproducing an ML convention.

Parameters:

  • marginal_scales (ndarray)

  • marginal_mom_scales (ndarray)

  • joint_mom_scales (ndarray)

  • adjusted_scales (ndarray)

  • negative_mom_fallback (ndarray)

evo_lmm.rare_effect_mom_ratio(marginal_scales, marginal_mom_scales, joint_mom_scales)#

Apply RareEffect’s marginal-ML × joint-MoM/marginal-MoM adjustment.

A non-positive marginal or joint MoM estimate triggers the published unadjusted-marginal fallback rather than silently truncating the result.

Parameters:

  • marginal_scales (ndarray)

  • marginal_mom_scales (ndarray)

  • joint_mom_scales (ndarray)

Return type:

RareEffectBaselineResult

class evo_lmm.HeritabilityEstimates(rare_effect: 'float', evolutionary: 'float', rare_genetic_variance: 'float', evolutionary_genetic_variance: 'float')#

Bases: object

Parameters:

  • rare_effect (float)

  • evolutionary (float)

  • rare_genetic_variance (float)

  • evolutionary_genetic_variance (float)

class evo_lmm.GeneComponentReport(gene, pooled_tau, sigma_b2_by_category)#

Bases: object

Gene-level empirical-Bayes report with pooled category shapes.

Parameters:

  • gene (Any)

  • pooled_tau (dict[Any, float])

  • sigma_b2_by_category (dict[Any, float])

class evo_lmm.FitReport(heritability, heritability_se, component_standard_errors, maf_decomposition, tau_profiles=None)#

Bases: object

Integrated estimand report with covariance-derived uncertainty.

Parameters:

  • heritability (HeritabilityEstimates)

  • heritability_se (float)

  • component_standard_errors (dict[str, float])

  • maf_decomposition (dict[Any, ndarray] | None)

  • tau_profiles (dict[Any, ProfileLikelihood] | None)

class evo_lmm.ProfileLikelihood(tau, objective, lower, upper)#

Bases: object

One-dimensional profile likelihood on a scientific parameter scale.

Parameters:

  • tau (ndarray)

  • objective (ndarray)

  • lower (float)

  • upper (float)

evo_lmm.boundary_lrt_pvalue(statistic, added_boundaries=1)#

Mixture-null p-value for independent non-negative boundary coordinates.

Parameters:

  • statistic (float)

  • added_boundaries (int)

Return type:

float

evo_lmm.delta_method_se(function, estimate, covariance, step=1e-05)#

Finite-difference delta-method standard error from an AI covariance.

Parameters:

  • function (Callable[[ndarray], float])

  • estimate (ndarray)

  • covariance (ndarray)

  • step (float)

Return type:

float

evo_lmm.genic_variance_by_maf(ops, prior, bins)#

Decompose projected genic variance by MAF bin and category.

Parameters:
Return type:

dict[Any, ndarray]

evo_lmm.heritability_conventions(ops, prior, sigma_e2)#

Return RareEffect’s uncentered-n and evo-lmm’s projected-d h².

Parameters:
Return type:

HeritabilityEstimates

evo_lmm.gene_component_report(gene, pooled_tau, sigma_b2_by_category)#

Construct the pooled-shape/per-gene-scale reporting unit.

Parameters:

  • gene (Any)

  • pooled_tau (dict[Any, float])

  • sigma_b2_by_category (dict[Any, float])

Return type:

GeneComponentReport

evo_lmm.fit_genes(genes, phenotype, pooled_tau, *, max_iter=100, trace_method='hutchinson', trace_probes=12)#

Fit per-gene scales with the category shapes pooled across genes.

pooled_tau is held fixed for every gene (fit_tau=False); it is a pooled estimate, not a per-gene starting point. Only the |c| scale coordinates are searched, so the reported pooled_tau is the shape the per-gene scales were actually conditioned on.

Parameters:

  • genes (dict[Any, MultiComponentOps])

  • phenotype (ndarray)

  • pooled_tau (dict[Any, float])

  • max_iter (int)

  • trace_method (str)

  • trace_probes (int)

Return type:

dict[Any, GeneComponentReport]

evo_lmm.fit_report(fit, *, maf_bins=None)#

Convert a fit into both estimands plus delta-method h² uncertainty.

Parameters:
Return type:

FitReport

evo_lmm.fit_tau_profiles(fit, tau_grids)#

Profile each category’s tau_c with other fitted parameters fixed.

Parameters:
Return type:

dict[Any, ProfileLikelihood]

evo_lmm.fit_parameter_profiles(fit, parameter_grids)#

Profile scientific-scale sigma_b2[label] and tau[label] grids.

Parameters:
Return type:

dict[str, ProfileLikelihood]

evo_lmm.profile_tau(tau_values, objective, *, confidence=0.95)#

Evaluate a tau profile and use the likelihood-ratio cutoff for bounds.

Parameters:

  • tau_values (Sequence[float])

  • objective (Callable[[float], float])

  • confidence (float)

Return type:

ProfileLikelihood

Priors#

class evo_lmm.SimplifiedPrior(sigma_b2, tau)

Bases: EvolutionaryPrior

The exact rho^2 = 1 evolutionary prior.

Parameters:

  • sigma_b2 (float) – Per-locus focal-trait effect variance on raw dosage units.

  • tau (float) – Non-negative sigma_a^2 / W_S aggregate. tau=0 is the frequency-independent boundary and is useful for tests and diagnostics.

class evo_lmm.FullPrior(sigma_b2, tau, rho)

Bases: EvolutionaryPrior

The full evolutionary prior with identifiable coupling rho^2.

Only rho^2 appears in the conditional variance, so the sign of rho is not identifiable. The public parameter accepts -1 <= rho <= 1 and the optimizer works directly with r = rho^2.

Parameters:

  • sigma_b2 (float)

  • tau (float)

  • rho (float)

Operators and data boundaries#

class evo_lmm.EvolutionaryLmmOps(chromosomes, frequencies=None, covariates=None, *, sample_filter=None, model='simplified', mutation_filter=None)

Bases: object

Matrix-free projected raw-dosage operators for evolutionary LMMs.

Parameters:

  • chromosomes (Any) – An (N, M) dense dosage matrix, a mapping of chromosome labels to matrices/GRGs, or a sequence of matrices/(label, source) pairs.

  • frequencies (Any) – Sample allele frequencies, supplied as one vector, one vector per chromosome, or a mapping. GRG inputs may omit this and extract them through GRAPP’s frequency traversal.

  • covariates (np.ndarray | None) – Fixed-effect design. An intercept is added when absent and the basis is QR-orthonormalised once.

  • model (str) – Coordinate interpretation for transformed phi arrays.

  • sample_filter (Sequence[int] | None)

  • mutation_filter (Mapping[Any, Sequence[int]] | None)

classmethod from_dense(genotypes, frequencies, covariates=None, *, chrom_labels=None, model='simplified')

Convenience constructor for dense tests and small simulations.

Parameters:

  • genotypes (ndarray | Sequence[ndarray])

  • frequencies (ndarray | Sequence[ndarray])

  • covariates (ndarray | None)

  • chrom_labels (Sequence[Any] | None)

  • model (str)

Return type:

EvolutionaryLmmOps

test_stats(chrom)

Return prior-independent tested-genotype statistics for a chromosome.

These are the quantities the BOLT-style calibration and association formulas need: the projected raw-dosage norms, the BOLT normalisation, and the eligibility mask. They depend only on the genotypes and the covariate basis, never on the fitted evolutionary prior.

Parameters:

chrom (Any)

Return type:

TestVariantStats

apply_model_x(weights, theta=None, exclude_chrom=None)

Apply P_C X diag(weights) to variant coefficients.

weights is normally a vector of model coefficients. When theta is a prior, it applies the model operator P_C X diag(sqrt(w(theta))) to those coefficients. A raw coefficient vector or chromosome mapping may also be supplied without theta. This method deliberately contains no sample-standardisation or 1/M normalization.

Parameters:

  • weights (ndarray | Mapping[Any, ndarray] | EvolutionaryPrior)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

model_scores(vector, theta=None, exclude_chrom=None)

Return model scores, optionally for B_theta^T = diag(sqrt(w))X^T P_C.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_k(vector, theta, exclude_chrom=None)

Apply K = P_C X diag(w(theta)) X^T P_C.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dk(vector, theta, parameter, exclude_chrom=None)

Apply a first derivative of K in a transformed shape coordinate.

Parameters:

  • vector (ndarray)

  • theta (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_h(vector, phi, exclude_chrom=None)

Apply the projected shape matrix H = K + delta I.

Parameters:

  • vector (ndarray)

  • phi (Any)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dh(vector, phi, parameter, exclude_chrom=None)

Apply a first derivative of H in log_delta, log_tau, or logit_r.

Parameters:

  • vector (ndarray)

  • phi (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_dh_matmat(values, phi, parameter, exclude_chrom=None)

Apply a first derivative of H to several columns at once.

This is the matrix-RHS counterpart of apply_dh(). Keeping the columns batched is important for GRG inputs: one matmat traversal replaces one traversal per stochastic trace probe.

Parameters:

  • values (ndarray)

  • phi (Any)

  • parameter (str)

  • exclude_chrom (Any)

Return type:

ndarray

apply_k_matmat(values, theta, exclude_chrom=None)

Apply K to batched right-hand sides without per-column traversals.

Parameters:

  • values (ndarray)

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

solve_ph(rhs_columns, phi, exclude_chrom=None, *, tol=0.0005, max_iter=None, initial=None, stats=None)

Solve projected H z = P_C rhs for one or many right-hand sides.

initial is an optional warm start. Each column is validated against a zero start independently; a poor or invalid cached column is reset without affecting the other right-hand sides. If stats is passed, it is populated with aggregate iteration and warm-start diagnostics.

Parameters:

  • rhs_columns (ndarray)

  • phi (Any)

  • exclude_chrom (Any)

  • tol (float)

  • max_iter (int | None)

  • initial (ndarray | None)

  • stats (dict[str, Any] | None)

Return type:

ndarray

test_scores(chrom, vector)

Return BOLT-normalised test-genotype scores for one chromosome.

Parameters:

  • chrom (Any)

  • vector (ndarray)

Return type:

ndarray

test_column(chrom, local_idx)

Return one projected, BOLT-normalised test-genotype column.

Parameters:

  • chrom (Any)

  • local_idx (int)

Return type:

ndarray

chromosome_frequencies(chrom)

Return sample allele frequencies for one chromosome in operator order.

Parameters:

chrom (Any)

Return type:

ndarray

local_indices(chrom)

Return mutation identifiers for a chromosome in operator order.

Parameters:

chrom (Any)

Return type:

ndarray

kernel_trace(theta, exclude_chrom=None)

Return tr(P_C K P_C) without constructing a dense kernel.

Parameters:

  • theta (Any)

  • exclude_chrom (Any)

Return type:

float

dense_kernel(theta, exclude_chrom=None)

Materialise a kernel for a small dense input or a test oracle.

Parameters:

  • theta (Any)

  • exclude_chrom (Any)

Return type:

ndarray

class evo_lmm.VariantData(frequencies, local_idx, raw_norm2=None, centered_norm2=None)

Bases: object

Cached per-variant information used by an evolutionary operator.

frequencies are sample allele frequencies and local_idx are the mutation identifiers in the source chromosome. raw_norm2 and centered_norm2 are optional because GRGL operators can obtain them via graph traversals when needed.

Parameters:

  • frequencies (ndarray)

  • local_idx (ndarray)

  • raw_norm2 (ndarray | None)

  • centered_norm2 (ndarray | None)

Fitting results#

class evo_lmm.FitResult(prior, sigma_b2, sigma_e2, delta, h2, log_likelihood, fixed_effects, projected_phenotype, ph_y, diagnostics, model, ops=None)

Bases: object

Scientific-scale fit result for an evolutionary LMM.

Parameters:

  • prior (EvolutionaryPrior)

  • sigma_b2 (float)

  • sigma_e2 (float)

  • delta (float)

  • h2 (float)

  • log_likelihood (float)

  • fixed_effects (ndarray)

  • projected_phenotype (ndarray)

  • ph_y (ndarray)

  • diagnostics (FitDiagnostics)

  • model (str)

  • ops (Any)

property sigma_g2: float

Compatibility alias; unlike GRAPP this is raw-effect sigma_b2.

blup()

Return the projected genetic-value BLUP sigma_b2 K P_V y.

Return type:

ndarray

class evo_lmm.FitDiagnostics(convergence, trace_estimator, trace_probes, objective=nan, initialization='default', trace_operator_queries=0, trace_standard_errors=<factory>, cg_iterations=<factory>, cg_warm_start_hits=0, cg_warm_start_rejections=0, cg_initial_residual_norms=<factory>, cg_final_residual_norms=<factory>, random_seed=None, boundary_hits=(), warnings=())

Bases: object

Numerical and identifiability diagnostics from a fit.

Convergence lives in ConvergenceReport under convergence; the flat accessors below delegate to it so both fitters expose the same names.

objective is the profiled REML objective or ``nan``. It is never a stand-in for something else: the stochastic paths do not evaluate a log-determinant, so they report nan rather than a surrogate, and convergence is judged by convergence.step_se_norm in every path.

Parameters:

  • convergence (ConvergenceReport)

  • trace_estimator (str)

  • trace_probes (int)

  • objective (float)

  • initialization (str)

  • trace_operator_queries (int)

  • trace_standard_errors (dict[str, float])

  • cg_iterations (list[int])

  • cg_warm_start_hits (int)

  • cg_warm_start_rejections (int)

  • cg_initial_residual_norms (list[float])

  • cg_final_residual_norms (list[float])

  • random_seed (int | None)

  • boundary_hits (tuple[str, ...])

  • warnings (tuple[str, ...])

class evo_lmm.AssociationResult(chrom, local_idx, score, beta, se, chisq, pvalue, chisq_linreg=None, pvalue_linreg=None, model_mask=None, frequencies=None, inverse_scale=nan, calibration_factor=nan)

Bases: object

Compact BOLT-compatible association output for one chromosome block.

beta and se are in raw diploid-dosage effect units, matching the evolutionary model’s sigma_b2 scale. score is the calibrated inverse-variance score x_j^T V_loco^-1 y on the BOLT-normalised test column. Entries excluded by model_mask (monomorphic or covariate-collinear columns) carry nan statistics and pvalue = 1.

Parameters:

  • chrom (Any)

  • local_idx (ndarray)

  • score (ndarray)

  • beta (ndarray)

  • se (ndarray)

  • chisq (ndarray)

  • pvalue (ndarray)

  • chisq_linreg (ndarray | None)

  • pvalue_linreg (ndarray | None)

  • model_mask (ndarray | None)

  • frequencies (ndarray | None)

  • inverse_scale (float)

  • calibration_factor (float)

good()

Return the boolean mask of variants with reportable statistics.

Return type:

ndarray

Association and calibration#

class evo_lmm.CalibrationResult(factor, std, ratio_of_medians, median_of_ratios, selected, tried, prospective, retrospective, inverse_scale, residuals=<factory>, quadratic_form=<factory>, screen_threshold=5.0, seed=0)

Bases: object

Prospective/retrospective calibration of the evolutionary LMM statistic.

Parameters:

  • factor (float)

  • std (float)

  • ratio_of_medians (float)

  • median_of_ratios (float)

  • selected (tuple[tuple[Any, int], ...])

  • tried (int)

  • prospective (ndarray)

  • retrospective (ndarray)

  • inverse_scale (dict[Any, float])

  • residuals (dict[Any, ndarray])

  • quadratic_form (dict[Any, float])

  • screen_threshold (float)

  • seed (int)

factor

The applied calibration factor: the ratio of prospective to retrospective statistic sums, or the ratio of medians when the jackknife standard error exceeds CALIBRATION_STD_LIMIT.

Type:

float

inverse_scale

Per-chromosome VinvScaleFactor in raw-effect units. A tested chromosome’s calibrated inverse-variance score is (x_j^T H_loco^-1 P y) / sigma_b2 / inverse_scale[chrom].

Type:

dict[Any, float]

residuals

H_loco^-1 P_C y for each left-out chromosome, reused by evo_lmm.association() so the solves are not repeated.

Type:

dict[Any, numpy.ndarray]

class evo_lmm.TestVariantStats(chrom, local_idx, centered_norm2, projected_norm2, test_scale, norm_scale, std_projected_norm2, model_mask)

Bases: object

Per-variant tested-genotype statistics for one chromosome block.

All quantities are computed after the covariate projection P_C. The test_scale is the BOLT normalisation sqrt(centered_norm2 / (n - 1)) applied to raw diploid dosage columns, so norm_scale = 1 / test_scale converts a raw-dosage quantity into BOLT-normalised units. Nothing here depends on the fitted evolutionary prior.

Parameters:

  • chrom (Any)

  • local_idx (ndarray)

  • centered_norm2 (ndarray)

  • projected_norm2 (ndarray)

  • test_scale (ndarray)

  • norm_scale (ndarray)

  • std_projected_norm2 (ndarray)

  • model_mask (ndarray)

projected_norm2

||P_C x_j||^2 in raw diploid dosage units.

Type:

numpy.ndarray

std_projected_norm2

||P_C x_j||^2 in BOLT-normalised units, i.e. projected_norm2 * norm_scale**2.

Type:

numpy.ndarray

model_mask

Variants eligible for association output; monomorphic and covariate-collinear columns are excluded exactly as in GRAPP.

Type:

numpy.ndarray

evo_lmm.calibrate_association(result, y=None, *, count=30, seed=0, screen_threshold=5.0, cg_tol=0.0005, stats=None)

Calibrate the LMM statistic against the fitted evolutionary covariance.

For every selected variant j on its own chromosome c the retrospective and prospective statistics are

retro = (N - C) * (x_j^T u_c)^2 / (||u_c||^2 * ||x_j||^2)

pro   = (N - C) * (x_j^T u_c)^2 / (x_j^T H_c^-1 x_j) / (y^T u_c)

with u_c = H_c^-1 P_C y and H_c the LOCO evolutionary shape matrix. Both are invariant to sigma_b2, so the factor is a pure shape correction; sigma_b2 enters only through inverse_scale.

Parameters:

  • result (FitResult)

  • y (ndarray | None)

  • count (int)

  • seed (int)

  • screen_threshold (float)

  • cg_tol (float)

  • stats (dict[Any, TestVariantStats] | None)

Return type:

CalibrationResult

evo_lmm.select_calibration_variants(result, *, count=30, seed=0, screen_threshold=5.0, stats=None)

Select calibration variants with BOLT’s blocked GRAMMAR pre-screen.

Eligible variants are split into count equal blocks in chromosome order. One variant is drawn uniformly from each block and accepted when its all-chromosome GRAMMAR retrospective statistic falls below screen_threshold, which keeps strongly associated variants out of the calibration set. Returns the selected (chrom, local_idx) pairs and the number of candidates examined.

Parameters:
Return type:

tuple[list[tuple[Any, int]], int]

evo_lmm.association(result, y=None, *, use_loco=True, calibrate=True, calibration=None, calibration_variants=30, seed=0, screen_threshold=5.0, cg_tol=0.0005)

Compute calibrated BOLT-style association statistics per chromosome.

The mixed-model statistic uses the fitted evolutionary LOCO covariance V_loco = sigma_b2 * (K_evo,loco + delta I) while the tested columns keep GRAPP’s independent BOLT normalisation. With calibrate=True the prospective/retrospective moment matching of evo_lmm.calibrate_association() supplies the per-chromosome inverse scale; otherwise the uncalibrated (factor = 1) scale is used, which is only appropriate for diagnostics.

beta and se are returned in raw diploid-dosage units. A single-variant linear-regression chi-square is reported alongside the mixed model statistic so inflation can be compared directly.

Parameters:

  • result (FitResult)

  • y (ndarray | None)

  • use_loco (bool)

  • calibrate (bool)

  • calibration (CalibrationResult | None)

  • calibration_variants (int)

  • seed (int)

  • screen_threshold (float)

  • cg_tol (float)

Return type:

list[AssociationResult]

evo_lmm.association_summary(results)

Return mean chi-square and lambda_GC for the mixed model and linreg.

lambda_gc divides the median chi-square by the median of a one-degree-of-freedom chi-square, so a calibrated null gives values near one for both statistics.

Parameters:

results (Sequence[AssociationResult])

Return type:

dict[str, dict[str, float]]

On this page
Show Source