JiaruiLi/SpectraREML
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add arbitrary fixed-effect contrasts

Browse Source
This commit is contained in:
Jiarui Li
2026-08-13 13:39:33 +08:00
parent 8557daecac
commit 48ab77f9ab
8 changed files with 515 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
README.md
View File

@@ -27,6 +27,8 @@ The engine diagonalizes the GRM once, rotates the common design and every unique
- Variable numbers of extra fixed-effect covariates per task through a CSR index.
- Satterthwaite t/F fixed-effect inference by default, with optional
Kenward-Roger F tests and a joint test of task-specific covariates.
- Named arbitrary linear contrasts over the common fixed-effect design, using
the Satterthwaite or Kenward-Roger method selected for the run.
- Atomic block output and safe resume/force semantics.
- Generic Python CLI for manifest creation, validation, execution, provenance signatures, status, and result export.
@@ -114,11 +116,14 @@ python python/spectra_reml.py make-manifest \
--tasks data/tasks.tsv \
--extra-offsets data/extra_offsets.i64.bin \
--extra-indices data/extra_indices.i32.bin \
--contrast-matrix data/contrasts.f64.bin \
--contrast-metadata data/contrasts.tsv \
--output-dir work/blocks \
--n-samples 3523 \
--n-base-covariates 8 \
--n-phenotype-rows 10000 \
--n-extra-covariate-rows 2000
--n-extra-covariate-rows 2000 \
--n-contrasts 5
```
Validate, run, and export results:
@@ -141,6 +146,14 @@ python python/spectra_reml.py finalize \
`results.tsv.gz` retains the complete per-task summary and stores the fixed-effect vector and row-wise packed lower covariance as JSON arrays.
Optional named contrasts are supplied by `--contrast-matrix`,
`--contrast-metadata`, and `--n-contrasts`. Contrast rows act on `base_x`;
task-specific covariate columns receive zero weights automatically. Finalized
results contain contrast estimates, adjusted standard errors, statistics,
numerator and denominator degrees of freedom, and p-values as JSON arrays.
This release uses manifest v2, run-signature v2, block v3, and finalized-output
v2 contracts. It intentionally does not read older contracts.
Fixed-effect inference defaults to Satterthwaite. Select Kenward-Roger or turn
inference off with `--fixed-effect-test kenward-roger` or
`--fixed-effect-test none`. Coefficient-wise standard errors, statistics,

45
docs/FORMAT.md
View File

@@ -84,6 +84,22 @@ indices[offsets[i]:offsets[i+1]]
An index must be in `[0, n_extra_covariate_rows)`, and a task cannot reference the same row twice.
## Named shared fixed-effect contrasts
Optional arbitrary hypotheses use two files. `--contrast-matrix` is a
row-major `float64` matrix with shape `n_contrasts x n_base_covariates`; its
columns follow `base_x` exactly. Task-specific extra fixed-effect columns are
assigned zero weights automatically. `--contrast-metadata` is a UTF-8 TSV:
```text
contrast_index contrast_id rhs
0 group_a_b 0
```
Indices are consecutive and zero based, IDs are unique, and every matrix row
must contain a nonzero weight. Each row tests `L beta = rhs` with the method
selected by `--fixed-effect-test` (Satterthwaite or Kenward-Roger).
## Block output
For block number `KKKKKK`:
@@ -96,6 +112,12 @@ block_KKKKKK.fixed_se.f64.bin
block_KKKKKK.fixed_stat.f64.bin
block_KKKKKK.fixed_ddf.f64.bin
block_KKKKKK.fixed_p.f64.bin
block_KKKKKK.contrast_estimate.f64.bin
block_KKKKKK.contrast_se.f64.bin
block_KKKKKK.contrast_stat.f64.bin
block_KKKKKK.contrast_numdf.f64.bin
block_KKKKKK.contrast_ddf.f64.bin
block_KKKKKK.contrast_p.f64.bin
block_KKKKKK.complete
```
@@ -119,6 +141,7 @@ grad_inf
fixed_test_method
fixed_test_status
fixed_test_offset
contrast_test_offset
extra_joint_num_df
extra_joint_den_df
extra_joint_f
@@ -142,12 +165,13 @@ The covariance array uses the row-wise packed lower triangle:
The `.complete` marker is written last and contains tab-separated key/value rows:
```text
format spectra-reml-block-v2
format spectra-reml-block-v3
block 0
tasks 256
beta_elements 4096
cov_elements 34816
fixed_test_elements 4096
contrast_test_elements 1280
```
Consumers must ignore blocks without `.complete`.
@@ -180,6 +204,13 @@ one numerator degree of freedom. For every task with extra covariates, the
summary also contains an F test of the joint null that all task-specific fixed
effects are zero.
The six `contrast_*` arrays share `contrast_test_offset` and follow the
contrast metadata order. A one-row Satterthwaite contrast reports a signed t
statistic; Kenward-Roger reports an F statistic with one numerator degree of
freedom. Both methods retain the signed `L beta` estimate and their adjusted
standard error. Negative `contrast_test_offset` means that no valid contrast
tests were emitted for that task.
`fixed_test_status` is one of:
```text
@@ -218,9 +249,13 @@ fixed_effect_se_json
fixed_effect_statistic_json
fixed_effect_denominator_df_json
fixed_effect_p_value_json
contrast_estimate_json
contrast_standard_error_json
contrast_statistic_json
contrast_numerator_df_json
contrast_denominator_df_json
contrast_p_value_json
```
Project-specific software can attach coefficient names. The public C++
inference API also accepts general linear hypotheses `L beta = rhs`; the batch
format currently emits coefficient-wise tests and the joint extra-covariate
test.
Contrast JSON arrays follow `contrast_index`; IDs and right-hand sides are
stored in the manifest-referenced contrast metadata TSV.

3
include/spectra_reml/batch_io.hpp
View File

@@ -20,6 +20,8 @@ struct BatchInputPaths {
std::filesystem::path tasks_tsv;
std::filesystem::path extra_offsets_i64;
std::filesystem::path extra_indices_i32;
std::filesystem::path contrast_matrix_f64;
std::filesystem::path contrast_metadata_tsv;
std::filesystem::path output_directory;
};
@@ -28,6 +30,7 @@ struct BatchDimensions {
std::size_t base_covariate_count = 0;
std::size_t phenotype_row_count = 0;
std::size_t extra_covariate_row_count = 0;
std::size_t contrast_count = 0;
};
struct BatchOptions {

4
include/spectra_reml/types.hpp
View File

@@ -164,6 +164,10 @@ struct RemlResult {
std::vector<double> beta_covariance_packed_lower;
FixedEffectInferenceResult fixed_effect_inference;
FixedEffectTestResult extra_fixed_effect_joint_test;
// Named shared contrasts supplied by the batch interface, in manifest
// contrast order. Every contrast acts on the common/base fixed effects;
// task-specific fixed-effect columns are assigned zero weight.
std::vector<FixedEffectTestResult> contrast_tests;
std::string error;
[[nodiscard]] bool has_estimates() const noexcept {

147
python/spectra_reml.py
View File

@@ -31,10 +31,10 @@ except ImportError as exc: # pragma: no cover
raise SystemExit("NumPy is required: {}".format(exc))
MANIFEST_FORMAT = "spectra-reml-manifest-v1"
BLOCK_FORMAT = "spectra-reml-block-v2"
RUN_SIGNATURE_FORMAT = "spectra-reml-run-signature-v1"
FINALIZE_FORMAT = "spectra-reml-finalize-v1"
MANIFEST_FORMAT = "spectra-reml-manifest-v2"
BLOCK_FORMAT = "spectra-reml-block-v3"
RUN_SIGNATURE_FORMAT = "spectra-reml-run-signature-v2"
FINALIZE_FORMAT = "spectra-reml-finalize-v2"
TASK_HEADER = (
"task_index",
"task_id",
@@ -59,6 +59,7 @@ SUMMARY_HEADER = (
"fixed_test_method",
"fixed_test_status",
"fixed_test_offset",
"contrast_test_offset",
"extra_joint_num_df",
"extra_joint_den_df",
"extra_joint_f",
@@ -75,6 +76,12 @@ BLOCK_PATTERNS = (
"block_*.fixed_stat.f64.bin",
"block_*.fixed_ddf.f64.bin",
"block_*.fixed_p.f64.bin",
"block_*.contrast_estimate.f64.bin",
"block_*.contrast_se.f64.bin",
"block_*.contrast_stat.f64.bin",
"block_*.contrast_numdf.f64.bin",
"block_*.contrast_ddf.f64.bin",
"block_*.contrast_p.f64.bin",
"block_*.complete",
)
@@ -180,6 +187,8 @@ def manifest_paths(manifest: Mapping[str, Any]) -> dict[str, Path | None]:
"extra_offsets_i64",
"extra_indices_i32",
"output_directory",
"contrast_matrix_f64",
"contrast_metadata_tsv",
)
missing = [name for name in required if name not in raw]
if missing:
@@ -199,6 +208,30 @@ def read_tasks(path: Path) -> list[dict[str, str]]:
return list(reader)
def read_contrast_metadata(path: Path) -> list[dict[str, str]]:
with path.open("r", encoding="utf-8-sig", newline="") as handle:
reader = csv.DictReader(handle, delimiter="\t")
expected = ("contrast_index", "contrast_id", "rhs")
if tuple(reader.fieldnames or ()) != expected:
fail("Contrast metadata header must be exactly: {}".format("\t".join(expected)))
rows = list(reader)
seen: set[str] = set()
for expected_index, row in enumerate(rows):
try:
index = int(row["contrast_index"])
rhs = float(row["rhs"])
except ValueError as exc:
fail("Invalid contrast metadata at row {}: {}".format(expected_index, exc))
if index != expected_index:
fail("Contrast indices must be consecutive and zero based.")
if not row["contrast_id"] or row["contrast_id"] in seen:
fail("Contrast IDs must be nonempty and unique.")
if not np.isfinite(rhs):
fail("Contrast rhs must be finite.")
seen.add(row["contrast_id"])
return rows
def validate_manifest(path: Path) -> dict[str, Any]:
path = resolved(path)
manifest = read_json(path)
@@ -214,6 +247,7 @@ def validate_manifest(path: Path) -> dict[str, Any]:
"extra_covariate_row_count",
"task_count",
"extra_index_count",
"contrast_count",
)
parsed: dict[str, int] = {}
for name in names:
@@ -240,6 +274,28 @@ def validate_manifest(path: Path) -> dict[str, Any]:
with paths["grm_id"].open("r", encoding="utf-8-sig") as handle:
if sum(bool(line.strip()) for line in handle) != n:
fail("GRM ID row count differs from sample_count.")
contrast_paths = (paths["contrast_matrix_f64"], paths["contrast_metadata_tsv"])
if parsed["contrast_count"]:
if any(path is None for path in contrast_paths):
fail("Both contrast_matrix_f64 and contrast_metadata_tsv are required.")
require_size(
paths["contrast_matrix_f64"],
parsed["contrast_count"] * parsed["base_covariate_count"] * 8,
"contrast matrix",
)
contrast_rows = read_contrast_metadata(paths["contrast_metadata_tsv"])
if len(contrast_rows) != parsed["contrast_count"]:
fail("Contrast metadata row count differs from contrast_count.")
contrast_matrix = np.fromfile(paths["contrast_matrix_f64"], dtype="<f8")
contrast_matrix = contrast_matrix.reshape(
parsed["contrast_count"], parsed["base_covariate_count"]
)
if not np.all(np.isfinite(contrast_matrix)):
fail("Contrast matrix contains a non-finite value.")
if np.any(np.all(contrast_matrix == 0.0, axis=1)):
fail("Every contrast row must contain at least one nonzero weight.")
elif any(path is not None for path in contrast_paths):
fail("Contrast paths must be omitted when contrast_count is zero.")
sources = manifest.get("source_files")
if not isinstance(sources, dict) or not sources:
@@ -263,6 +319,9 @@ def validate_manifest(path: Path) -> dict[str, Any]:
"extra_offsets_i64": paths["extra_offsets_i64"],
"extra_indices_i32": paths["extra_indices_i32"],
}
if parsed["contrast_count"]:
controls["contrast_matrix_f64"] = paths["contrast_matrix_f64"]
controls["contrast_metadata_tsv"] = paths["contrast_metadata_tsv"]
hashes = manifest.get("control_sha256")
if not isinstance(hashes, dict) or set(hashes) != set(controls):
fail("control_sha256 does not cover all control files.")
@@ -314,6 +373,8 @@ def make_manifest(args: argparse.Namespace) -> Path:
"tasks_tsv": resolved(args.tasks),
"extra_offsets_i64": resolved(args.extra_offsets),
"extra_indices_i32": resolved(args.extra_indices),
"contrast_matrix_f64": resolved(args.contrast_matrix) if args.contrast_matrix else None,
"contrast_metadata_tsv": resolved(args.contrast_metadata) if args.contrast_metadata else None,
"output_directory": resolved(args.output_dir),
}
tasks = read_tasks(paths["tasks_tsv"]) # type: ignore[arg-type]
@@ -323,7 +384,13 @@ def make_manifest(args: argparse.Namespace) -> Path:
for name in ("grm_id", "extra_covariate_f32"):
if paths[name] is not None:
sources[name] = file_identity(paths[name]) # type: ignore[arg-type]
controls = ("base_x_f64", "tasks_tsv", "extra_offsets_i64", "extra_indices_i32")
controls = ["base_x_f64", "tasks_tsv", "extra_offsets_i64", "extra_indices_i32"]
if args.n_contrasts:
if paths["contrast_matrix_f64"] is None or paths["contrast_metadata_tsv"] is None:
fail("--contrast-matrix and --contrast-metadata are required when --n-contrasts is positive.")
controls.extend(("contrast_matrix_f64", "contrast_metadata_tsv"))
elif paths["contrast_matrix_f64"] is not None or paths["contrast_metadata_tsv"] is not None:
fail("Contrast paths require a positive --n-contrasts.")
manifest = {
"format": MANIFEST_FORMAT,
"created_utc": utc_now(),
@@ -334,6 +401,7 @@ def make_manifest(args: argparse.Namespace) -> Path:
"extra_covariate_row_count": args.n_extra_covariate_rows,
"task_count": len(tasks),
"extra_index_count": extra_index_count,
"contrast_count": args.n_contrasts,
},
"paths": {name: str(value) if value is not None else None for name, value in paths.items()},
"source_files": sources,
@@ -371,6 +439,7 @@ def build_engine_command(engine: Path, manifest: Mapping[str, Any], args: argpar
"--n-base-covariates", str(dims["base_covariate_count"]),
"--n-phenotype-rows", str(dims["phenotype_row_count"]),
"--n-extra-covariate-rows", str(dims["extra_covariate_row_count"]),
"--n-contrasts", str(dims["contrast_count"]),
"--block-size", str(args.block_size),
"--threads", str(args.threads),
"--fixed-effect-test", args.fixed_effect_test,
@@ -379,6 +448,9 @@ def build_engine_command(engine: Path, manifest: Mapping[str, Any], args: argpar
command.extend(("--grm-id", str(paths["grm_id"])))
if int(dims["extra_index_count"]):
command.extend(("--extra-covariates", str(paths["extra_covariate_f32"])))
if int(dims["contrast_count"]):
command.extend(("--contrast-matrix", str(paths["contrast_matrix_f64"])))
command.extend(("--contrast-metadata", str(paths["contrast_metadata_tsv"])))
option_map = (
("max_iterations", "--max-iterations"),
("line_search_max_evals", "--line-search-max-evals"),
@@ -571,16 +643,25 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
fixed_stat_path = output / (stem + ".fixed_stat.f64.bin")
fixed_ddf_path = output / (stem + ".fixed_ddf.f64.bin")
fixed_p_path = output / (stem + ".fixed_p.f64.bin")
contrast_estimate_path = output / (stem + ".contrast_estimate.f64.bin")
contrast_se_path = output / (stem + ".contrast_se.f64.bin")
contrast_stat_path = output / (stem + ".contrast_stat.f64.bin")
contrast_numdf_path = output / (stem + ".contrast_numdf.f64.bin")
contrast_ddf_path = output / (stem + ".contrast_ddf.f64.bin")
contrast_p_path = output / (stem + ".contrast_p.f64.bin")
marker_values = read_marker(marker)
block_format = marker_values.get("format")
try:
declared_tasks = int(marker_values["tasks"])
beta_elements = int(marker_values["beta_elements"])
cov_elements = int(marker_values["cov_elements"])
fixed_test_elements = int(marker_values["fixed_test_elements"])
contrast_test_elements = int(marker_values["contrast_test_elements"])
except (KeyError, ValueError) as exc:
fail("Malformed completion counts in {}: {}".format(marker, exc))
if marker_values.get("format") != BLOCK_FORMAT or min(
declared_tasks, beta_elements, cov_elements, fixed_test_elements
if block_format != BLOCK_FORMAT or min(
declared_tasks, beta_elements, cov_elements, fixed_test_elements,
contrast_test_elements,
) < 0:
fail("Invalid completion marker: {}".format(marker))
require_size(beta_path, beta_elements * 8, "block beta")
@@ -589,6 +670,15 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
require_size(fixed_stat_path, fixed_test_elements * 8, "fixed-effect statistic")
require_size(fixed_ddf_path, fixed_test_elements * 8, "fixed-effect denominator df")
require_size(fixed_p_path, fixed_test_elements * 8, "fixed-effect p-value")
for path, label in (
(contrast_estimate_path, "contrast estimate"),
(contrast_se_path, "contrast SE"),
(contrast_stat_path, "contrast statistic"),
(contrast_numdf_path, "contrast numerator df"),
(contrast_ddf_path, "contrast denominator df"),
(contrast_p_path, "contrast p-value"),
):
require_size(path, contrast_test_elements * 8, label)
with summary_path.open("r", encoding="utf-8-sig", newline="") as handle:
reader = csv.DictReader(handle, delimiter="\t")
if tuple(reader.fieldnames or ()) != SUMMARY_HEADER:
@@ -602,9 +692,17 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
fixed_stat = np.fromfile(fixed_stat_path, dtype="<f8")
fixed_ddf = np.fromfile(fixed_ddf_path, dtype="<f8")
fixed_p = np.fromfile(fixed_p_path, dtype="<f8")
contrast_estimate = np.fromfile(contrast_estimate_path, dtype="<f8")
contrast_se = np.fromfile(contrast_se_path, dtype="<f8")
contrast_stat = np.fromfile(contrast_stat_path, dtype="<f8")
contrast_numdf = np.fromfile(contrast_numdf_path, dtype="<f8")
contrast_ddf = np.fromfile(contrast_ddf_path, dtype="<f8")
contrast_p = np.fromfile(contrast_p_path, dtype="<f8")
next_beta_offset = 0
next_cov_offset = 0
next_fixed_test_offset = 0
next_contrast_test_offset = 0
q = int(manifest["dimensions"]["contrast_count"])
for row in rows:
index = int(row["task_index"])
if index in results or index < 0 or index >= task_count:
@@ -612,6 +710,7 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
p = int(row["n_fixed"])
beta_offset, cov_offset = int(row["beta_offset"]), int(row["cov_offset"])
fixed_test_offset = int(row["fixed_test_offset"])
contrast_test_offset = int(row["contrast_test_offset"])
if beta_offset < 0 or cov_offset < 0:
if beta_offset != cov_offset:
fail("Only one output offset is negative for task {}".format(index))
@@ -641,6 +740,25 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
fixed_ddf_values = fixed_ddf[fixed_test_offset:fixed_test_offset + p].tolist()
fixed_p_values = fixed_p[fixed_test_offset:fixed_test_offset + p].tolist()
next_fixed_test_offset += p
if contrast_test_offset < 0:
contrast_estimate_values: list[float] = []
contrast_se_values: list[float] = []
contrast_stat_values: list[float] = []
contrast_numdf_values: list[float] = []
contrast_ddf_values: list[float] = []
contrast_p_values: list[float] = []
else:
if (contrast_test_offset != next_contrast_test_offset or
contrast_test_offset + q > contrast_estimate.size):
fail("Invalid contrast test offset for task {}".format(index))
section = slice(contrast_test_offset, contrast_test_offset + q)
contrast_estimate_values = contrast_estimate[section].tolist()
contrast_se_values = contrast_se[section].tolist()
contrast_stat_values = contrast_stat[section].tolist()
contrast_numdf_values = contrast_numdf[section].tolist()
contrast_ddf_values = contrast_ddf[section].tolist()
contrast_p_values = contrast_p[section].tolist()
next_contrast_test_offset += q
result: dict[str, Any] = dict(row)
result["beta_json"] = json.dumps(beta_values, separators=(",", ":"))
result["covariance_packed_lower_json"] = json.dumps(cov_values, separators=(",", ":"))
@@ -648,9 +766,16 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]:
result["fixed_effect_statistic_json"] = json.dumps(fixed_stat_values, separators=(",", ":"))
result["fixed_effect_denominator_df_json"] = json.dumps(fixed_ddf_values, separators=(",", ":"))
result["fixed_effect_p_value_json"] = json.dumps(fixed_p_values, separators=(",", ":"))
result["contrast_estimate_json"] = json.dumps(contrast_estimate_values, separators=(",", ":"))
result["contrast_standard_error_json"] = json.dumps(contrast_se_values, separators=(",", ":"))
result["contrast_statistic_json"] = json.dumps(contrast_stat_values, separators=(",", ":"))
result["contrast_numerator_df_json"] = json.dumps(contrast_numdf_values, separators=(",", ":"))
result["contrast_denominator_df_json"] = json.dumps(contrast_ddf_values, separators=(",", ":"))
result["contrast_p_value_json"] = json.dumps(contrast_p_values, separators=(",", ":"))
results[index] = result
if (next_beta_offset != beta_elements or next_cov_offset != cov_elements or
next_fixed_test_offset != fixed_test_elements):
next_fixed_test_offset != fixed_test_elements or
next_contrast_test_offset != contrast_test_elements):
fail("Block binary arrays contain unused elements: {}".format(stem))
missing = sorted(set(range(task_count)).difference(results))
if missing:
@@ -669,6 +794,9 @@ def finalize(manifest: Mapping[str, Any], output: Path) -> Path:
"beta_json", "covariance_packed_lower_json", "fixed_effect_se_json",
"fixed_effect_statistic_json", "fixed_effect_denominator_df_json",
"fixed_effect_p_value_json",
"contrast_estimate_json", "contrast_standard_error_json",
"contrast_statistic_json", "contrast_numerator_df_json",
"contrast_denominator_df_json", "contrast_p_value_json",
]
try:
with opener(temporary, "wt", encoding="utf-8", newline="") as handle:
@@ -771,11 +899,14 @@ def build_parser() -> argparse.ArgumentParser:
make.add_argument("--tasks", type=Path, required=True)
make.add_argument("--extra-offsets", type=Path, required=True)
make.add_argument("--extra-indices", type=Path, required=True)
make.add_argument("--contrast-matrix", type=Path)
make.add_argument("--contrast-metadata", type=Path)
make.add_argument("--output-dir", type=Path, required=True)
make.add_argument("--n-samples", type=positive_integer, required=True)
make.add_argument("--n-base-covariates", type=positive_integer, required=True)
make.add_argument("--n-phenotype-rows", type=positive_integer, required=True)
make.add_argument("--n-extra-covariate-rows", type=int, default=0)
make.add_argument("--n-contrasts", type=int, default=0)
for action in ("validate", "status"):
command = subparsers.add_parser(action)
command.add_argument("--manifest", type=Path, required=True)

185
src/batch_io.cpp
View File

@@ -34,7 +34,96 @@
namespace spectra::reml {
namespace {
constexpr const char* kOutputFormatVersion = "spectra-reml-block-v2";
constexpr const char* kOutputFormatVersion = "spectra-reml-block-v3";
std::vector<std::string> split_tab(const std::string& line);
std::uint64_t parse_u64(const std::string& text, const char* field,
std::size_t line_number);
struct NamedContrast {
std::string id;
double rhs = 0.0;
std::vector<double> weights;
};
std::vector<NamedContrast> read_named_contrasts(
const BatchInputPaths& paths, const BatchDimensions& dimensions) {
if (dimensions.contrast_count == 0) {
return {};
}
const ColMajorMatrix matrix = read_row_major_f64_matrix(
paths.contrast_matrix_f64, dimensions.contrast_count,
dimensions.base_covariate_count);
std::ifstream input(paths.contrast_metadata_tsv);
if (!input) {
throw std::runtime_error("cannot open contrast metadata TSV " +
paths.contrast_metadata_tsv.string());
}
std::string line;
if (!std::getline(input, line)) {
throw std::runtime_error("contrast metadata TSV is empty");
}
if (!line.empty() && line.back() == '\r') line.pop_back();
if (split_tab(line) !=
std::vector<std::string>{"contrast_index", "contrast_id", "rhs"}) {
throw std::runtime_error(
"contrast metadata header must be contrast_index, contrast_id, rhs");
}
std::vector<NamedContrast> contrasts;
contrasts.reserve(dimensions.contrast_count);
std::unordered_set<std::string> ids;
std::size_t line_number = 1;
while (std::getline(input, line)) {
++line_number;
if (!line.empty() && line.back() == '\r') line.pop_back();
if (line.empty()) continue;
const auto fields = split_tab(line);
if (fields.size() != 3) {
throw std::runtime_error("invalid contrast metadata field count at line " +
std::to_string(line_number));
}
const auto index = parse_u64(fields[0], "contrast_index", line_number);
if (index != contrasts.size()) {
throw std::runtime_error(
"contrast indices must be consecutive and zero based");
}
if (fields[1].empty() || !ids.insert(fields[1]).second) {
throw std::runtime_error(
"contrast IDs must be nonempty and unique at line " +
std::to_string(line_number));
}
std::size_t consumed = 0;
double rhs = std::numeric_limits<double>::quiet_NaN();
try {
rhs = std::stod(fields[2], &consumed);
} catch (const std::exception&) {
throw std::runtime_error("invalid contrast rhs at line " +
std::to_string(line_number));
}
if (consumed != fields[2].size() || !std::isfinite(rhs)) {
throw std::runtime_error("invalid contrast rhs at line " +
std::to_string(line_number));
}
NamedContrast contrast;
contrast.id = fields[1];
contrast.rhs = rhs;
contrast.weights.resize(dimensions.base_covariate_count);
bool any_nonzero = false;
for (std::size_t col = 0; col < dimensions.base_covariate_count; ++col) {
contrast.weights[col] = matrix(index, col);
any_nonzero = any_nonzero || contrast.weights[col] != 0.0;
}
if (!any_nonzero) {
throw std::runtime_error("contrast " + contrast.id +
" has all-zero weights");
}
contrasts.push_back(std::move(contrast));
}
if (contrasts.size() != dimensions.contrast_count) {
throw std::runtime_error("contrast metadata row count mismatch");
}
return contrasts;
}
std::size_t checked_product(std::size_t left, std::size_t right,
const char* description) {
@@ -295,10 +384,24 @@ void write_block_atomic(const std::filesystem::path& output_directory,
output_directory / (stem + ".fixed_stat.f64.bin");
const auto fixed_ddf_path = output_directory / (stem + ".fixed_ddf.f64.bin");
const auto fixed_p_path = output_directory / (stem + ".fixed_p.f64.bin");
const auto contrast_estimate_path =
output_directory / (stem + ".contrast_estimate.f64.bin");
const auto contrast_se_path =
output_directory / (stem + ".contrast_se.f64.bin");
const auto contrast_stat_path =
output_directory / (stem + ".contrast_stat.f64.bin");
const auto contrast_numdf_path =
output_directory / (stem + ".contrast_numdf.f64.bin");
const auto contrast_ddf_path =
output_directory / (stem + ".contrast_ddf.f64.bin");
const auto contrast_p_path =
output_directory / (stem + ".contrast_p.f64.bin");
const auto complete_path = output_directory / (stem + ".complete");
const std::array<std::filesystem::path, 8> final_paths = {
const std::array<std::filesystem::path, 14> final_paths = {
summary_path, beta_path, covariance_path, fixed_se_path,
fixed_stat_path, fixed_ddf_path, fixed_p_path, complete_path};
fixed_stat_path, fixed_ddf_path, fixed_p_path, contrast_estimate_path,
contrast_se_path, contrast_stat_path, contrast_numdf_path,
contrast_ddf_path, contrast_p_path, complete_path};
for (const auto& path : final_paths) {
if (std::filesystem::exists(path)) {
if (!overwrite) {
@@ -314,9 +417,16 @@ void write_block_atomic(const std::filesystem::path& output_directory,
std::vector<double> fixed_stat_values;
std::vector<double> fixed_ddf_values;
std::vector<double> fixed_p_values;
std::vector<double> contrast_estimate_values;
std::vector<double> contrast_se_values;
std::vector<double> contrast_stat_values;
std::vector<double> contrast_numdf_values;
std::vector<double> contrast_ddf_values;
std::vector<double> contrast_p_values;
std::vector<std::int64_t> beta_offsets(results.size(), -1);
std::vector<std::int64_t> covariance_offsets(results.size(), -1);
std::vector<std::int64_t> fixed_test_offsets(results.size(), -1);
std::vector<std::int64_t> contrast_test_offsets(results.size(), -1);
for (std::size_t index = 0; index < results.size(); ++index) {
const auto& fit = results[index].fit;
if (!fit.has_estimates()) {
@@ -342,6 +452,19 @@ void write_block_atomic(const std::filesystem::path& output_directory,
fixed_ddf_values.push_back(test.denominator_df);
fixed_p_values.push_back(test.p_value);
}
if (!fit.contrast_tests.empty()) {
contrast_test_offsets[index] = static_cast<std::int64_t>(
contrast_estimate_values.size());
for (const auto& test : fit.contrast_tests) {
contrast_estimate_values.push_back(test.estimate);
contrast_se_values.push_back(test.standard_error);
contrast_stat_values.push_back(test.statistic);
contrast_numdf_values.push_back(
static_cast<double>(test.numerator_df));
contrast_ddf_values.push_back(test.denominator_df);
contrast_p_values.push_back(test.p_value);
}
}
}
}
@@ -352,10 +475,18 @@ void write_block_atomic(const std::filesystem::path& output_directory,
const auto fixed_stat_temp = temporary_path(fixed_stat_path);
const auto fixed_ddf_temp = temporary_path(fixed_ddf_path);
const auto fixed_p_temp = temporary_path(fixed_p_path);
const auto contrast_estimate_temp = temporary_path(contrast_estimate_path);
const auto contrast_se_temp = temporary_path(contrast_se_path);
const auto contrast_stat_temp = temporary_path(contrast_stat_path);
const auto contrast_numdf_temp = temporary_path(contrast_numdf_path);
const auto contrast_ddf_temp = temporary_path(contrast_ddf_path);
const auto contrast_p_temp = temporary_path(contrast_p_path);
const auto complete_temp = temporary_path(complete_path);
const std::array<std::filesystem::path, 8> temporary_paths = {
const std::array<std::filesystem::path, 14> temporary_paths = {
summary_temp, beta_temp, covariance_temp, fixed_se_temp,
fixed_stat_temp, fixed_ddf_temp, fixed_p_temp, complete_temp};
fixed_stat_temp, fixed_ddf_temp, fixed_p_temp, contrast_estimate_temp,
contrast_se_temp, contrast_stat_temp, contrast_numdf_temp,
contrast_ddf_temp, contrast_p_temp, complete_temp};
try {
write_binary_doubles(beta_temp, beta_values);
write_binary_doubles(covariance_temp, covariance_values);
@@ -363,6 +494,12 @@ void write_block_atomic(const std::filesystem::path& output_directory,
write_binary_doubles(fixed_stat_temp, fixed_stat_values);
write_binary_doubles(fixed_ddf_temp, fixed_ddf_values);
write_binary_doubles(fixed_p_temp, fixed_p_values);
write_binary_doubles(contrast_estimate_temp, contrast_estimate_values);
write_binary_doubles(contrast_se_temp, contrast_se_values);
write_binary_doubles(contrast_stat_temp, contrast_stat_values);
write_binary_doubles(contrast_numdf_temp, contrast_numdf_values);
write_binary_doubles(contrast_ddf_temp, contrast_ddf_values);
write_binary_doubles(contrast_p_temp, contrast_p_values);
{
std::ofstream summary(summary_temp, std::ios::trunc);
if (!summary) {
@@ -374,6 +511,7 @@ void write_block_atomic(const std::filesystem::path& output_directory,
"\tcov_offset\tsigma_g2\tsigma_e2\th2\tlogL\titerations"
"\tline_search_steps\tgrad_inf\tfixed_test_method"
"\tfixed_test_status\tfixed_test_offset"
"\tcontrast_test_offset"
"\textra_joint_num_df\textra_joint_den_df"
"\textra_joint_f\textra_joint_p\tfixed_test_error\terror\n";
summary << std::setprecision(17);
@@ -396,6 +534,7 @@ void write_block_atomic(const std::filesystem::path& output_directory,
<< to_string(fit.fixed_effect_inference.method) << '\t'
<< to_string(fit.fixed_effect_inference.status) << '\t'
<< fixed_test_offsets[index] << '\t'
<< contrast_test_offsets[index] << '\t'
<< fit.extra_fixed_effect_joint_test.numerator_df << '\t'
<< fit.extra_fixed_effect_joint_test.denominator_df << '\t'
<< fit.extra_fixed_effect_joint_test.statistic << '\t'
@@ -421,6 +560,9 @@ void write_block_atomic(const std::filesystem::path& output_directory,
<< "beta_elements\t" << beta_values.size() << '\n'
<< "cov_elements\t" << covariance_values.size() << '\n'
<< "fixed_test_elements\t" << fixed_se_values.size()
<< '\n'
<< "contrast_test_elements\t"
<< contrast_estimate_values.size()
<< '\n';
complete.flush();
if (!complete) {
@@ -434,6 +576,12 @@ void write_block_atomic(const std::filesystem::path& output_directory,
rename_checked(fixed_stat_temp, fixed_stat_path);
rename_checked(fixed_ddf_temp, fixed_ddf_path);
rename_checked(fixed_p_temp, fixed_p_path);
rename_checked(contrast_estimate_temp, contrast_estimate_path);
rename_checked(contrast_se_temp, contrast_se_path);
rename_checked(contrast_stat_temp, contrast_stat_path);
rename_checked(contrast_numdf_temp, contrast_numdf_path);
rename_checked(contrast_ddf_temp, contrast_ddf_path);
rename_checked(contrast_p_temp, contrast_p_path);
rename_checked(summary_temp, summary_path);
rename_checked(complete_temp, complete_path);
} catch (...) {
@@ -624,6 +772,7 @@ void run_task_batch(const BatchInputPaths& paths,
if (tasks.empty()) {
throw std::runtime_error("tasks TSV contains no tasks");
}
const auto named_contrasts = read_named_contrasts(paths, dimensions);
std::cerr << "Reading and diagonalizing the GRM for "
<< dimensions.sample_count << " samples...\n";
const SpectralGrm spectral = read_and_diagonalize_gcta_grm(
@@ -750,6 +899,7 @@ void run_task_batch(const BatchInputPaths& paths,
item.fit = fit_ai_reml_spectral(
phenotype, design, spectral.eigenvalues, options.reml);
std::vector<FixedEffectHypothesis> hypotheses;
const bool has_extra_joint = !task.extra_covariate_rows.empty();
if (!task.extra_covariate_rows.empty()) {
FixedEffectHypothesis extra_joint;
extra_joint.contrast = ColMajorMatrix(
@@ -761,14 +911,35 @@ void run_task_batch(const BatchInputPaths& paths,
}
hypotheses.push_back(std::move(extra_joint));
}
for (const auto& named : named_contrasts) {
FixedEffectHypothesis hypothesis;
hypothesis.contrast = ColMajorMatrix(1, fixed_count);
for (std::size_t col = 0;
col < dimensions.base_covariate_count; ++col) {
hypothesis.contrast(0, col) = named.weights[col];
}
hypothesis.rhs = {named.rhs};
hypotheses.push_back(std::move(hypothesis));
}
item.fit.fixed_effect_inference = infer_fixed_effects_spectral(
phenotype, design, spectral.eigenvalues, item.fit,
options.reml.fixed_effect_test, hypotheses,
options.reml.rank_tolerance_relative,
options.reml.covariance_floor_relative);
if (!item.fit.fixed_effect_inference.hypothesis_tests.empty()) {
const auto& tests =
item.fit.fixed_effect_inference.hypothesis_tests;
std::size_t contrast_begin = 0;
if (has_extra_joint && !tests.empty()) {
item.fit.extra_fixed_effect_joint_test =
item.fit.fixed_effect_inference.hypothesis_tests.front();
tests.front();
contrast_begin = 1;
}
if (tests.size() >= contrast_begin + named_contrasts.size()) {
item.fit.contrast_tests.assign(
tests.begin() + static_cast<std::ptrdiff_t>(contrast_begin),
tests.begin() + static_cast<std::ptrdiff_t>(
contrast_begin +
named_contrasts.size()));
}
} catch (const std::exception& exception) {
item.fit.status = FitStatus::numerical_error;

21
src/main.cpp
View File

@@ -29,12 +29,15 @@ void print_usage(std::ostream& output) {
<< " --tasks PATH TSV with task_index,task_id,phenotype_row,n_extra_covariates\n"
<< " --extra-offsets PATH int64 CSR offsets, m+1 entries\n"
<< " --extra-indices PATH int32 zero-based extra-covariate row indices\n"
<< " --contrast-matrix PATH optional row-major q x p0 float64 matrix\n"
<< " --contrast-metadata PATH optional TSV: contrast_index,contrast_id,rhs\n"
<< " --out-dir PATH task-set output directory\n\n"
<< "Required dimensions:\n"
<< " --n-samples N\n"
<< " --n-base-covariates P\n"
<< " --n-phenotype-rows R\n"
<< " --n-extra-covariate-rows R\n\n"
<< " --n-contrasts Q shared named contrast rows (default 0)\n\n"
<< "Conditionally required:\n"
<< " --extra-covariates PATH row-major rows x n float32 matrix; required\n"
<< " when any task uses extra covariates\n\n"
@@ -174,11 +177,14 @@ void reject_unknown_options(const Arguments& arguments) {
"--tasks",
"--extra-offsets",
"--extra-indices",
"--contrast-matrix",
"--contrast-metadata",
"--out-dir",
"--n-samples",
"--n-base-covariates",
"--n-phenotype-rows",
"--n-extra-covariate-rows",
"--n-contrasts",
"--block-size",
"--threads",
"--fixed-effect-test",
@@ -230,6 +236,8 @@ int main(int argc, char** argv) {
paths.tasks_tsv = require_value(arguments, "--tasks");
paths.extra_offsets_i64 = require_value(arguments, "--extra-offsets");
paths.extra_indices_i32 = require_value(arguments, "--extra-indices");
paths.contrast_matrix_f64 = optional_value(arguments, "--contrast-matrix");
paths.contrast_metadata_tsv = optional_value(arguments, "--contrast-metadata");
paths.output_directory = require_value(arguments, "--out-dir");
BatchDimensions dimensions;
@@ -244,6 +252,19 @@ int main(int argc, char** argv) {
dimensions.extra_covariate_row_count = parse_size(
require_value(arguments, "--n-extra-covariate-rows"),
"--n-extra-covariate-rows", true);
dimensions.contrast_count = parse_size(
optional_value(arguments, "--n-contrasts", "0"),
"--n-contrasts", true);
const bool any_contrast_path = !paths.contrast_matrix_f64.empty() ||
!paths.contrast_metadata_tsv.empty();
const bool all_contrast_paths = !paths.contrast_matrix_f64.empty() &&
!paths.contrast_metadata_tsv.empty();
if ((dimensions.contrast_count == 0 && any_contrast_path) ||
(dimensions.contrast_count > 0 && !all_contrast_paths)) {
throw std::invalid_argument(
"--contrast-matrix and --contrast-metadata are both required "
"when --n-contrasts is positive, and must be omitted otherwise");
}
BatchOptions options;
options.resume = arguments.flags.count("--resume") != 0U;

122
tests/python/test_spectra_reml_cli.py
View File

@@ -21,7 +21,7 @@ SPEC.loader.exec_module(cli)
class SpectraRemlCliTests(unittest.TestCase):
def fixture(self, root: Path) -> tuple[Path, dict]:
def fixture(self, root: Path, with_contrasts: bool = False) -> tuple[Path, dict]:
n, p0 = 3, 2
np.zeros(n * (n + 1) // 2, dtype="<f4").tofile(root / "grm.bin")
(root / "grm.id").write_text("F1 I1\nF2 I2\nF3 I3\n", encoding="utf-8")
@@ -36,9 +36,18 @@ class SpectraRemlCliTests(unittest.TestCase):
)
np.asarray([0, 0, 1], dtype="<i8").tofile(root / "offsets.bin")
np.asarray([0], dtype="<i4").tofile(root / "indices.bin")
if with_contrasts:
np.asarray([[0, 1], [1, -1]], dtype="<f8").tofile(
root / "contrasts.bin"
)
(root / "contrasts.tsv").write_text(
"contrast_index\tcontrast_id\trhs\n"
"0\tslope\t0\n"
"1\tintercept_minus_slope\t0\n",
encoding="utf-8",
)
manifest_path = root / "manifest.json"
args = cli.build_parser().parse_args(
[
command = [
"make-manifest",
"--manifest", str(manifest_path),
"--grm-bin", str(root / "grm.bin"),
@@ -55,19 +64,31 @@ class SpectraRemlCliTests(unittest.TestCase):
"--n-phenotype-rows", "2",
"--n-extra-covariate-rows", "1",
]
if with_contrasts:
command.extend(
[
"--contrast-matrix", str(root / "contrasts.bin"),
"--contrast-metadata", str(root / "contrasts.tsv"),
"--n-contrasts", "2",
]
)
args = cli.build_parser().parse_args(command)
cli.make_manifest(args)
return manifest_path, cli.validate_manifest(manifest_path)
@staticmethod
def write_blocks(root: Path) -> None:
def write_blocks(root: Path, contrast_count: int = 0) -> None:
root.mkdir(parents=True, exist_ok=True)
first_contrast_offset = 0 if contrast_count else -1
second_contrast_offset = contrast_count if contrast_count else -1
(root / "block_000000.summary.tsv").write_text(
"\t".join(cli.SUMMARY_HEADER) + "\n"
"0\ttrait_a\tconverged\t2\t0\t0\t0\t1\t1\t0.5\t-1\t4\t8\t1e-8"
"\tsatterthwaite\tok\t0\t0\tnan\tnan\tnan\t\t\n"
"\tsatterthwaite\tok\t0\t{}\t0\tnan\tnan\tnan\t\t\n"
"1\ttrait_b\tconverged_boundary\t3\t1\t2\t3\t0\t1\t0\t-2\t3\t6\t1e-9"
"\tsatterthwaite\tboundary_conditional\t2\t1\t8\t2\t0.2\t\t\n",
"\tsatterthwaite\tboundary_conditional\t2\t{}\t1\t8\t2\t0.2\t\t\n".format(
first_contrast_offset, second_contrast_offset
),
encoding="utf-8",
)
np.asarray([1, 2, 3, 4, 5], dtype="<f8").tofile(root / "block_000000.beta.f64.bin")
@@ -84,9 +105,22 @@ class SpectraRemlCliTests(unittest.TestCase):
np.asarray([0.5, 0.2, 0.1, 0.05, 0.01], dtype="<f8").tofile(
root / "block_000000.fixed_p.f64.bin"
)
contrast_elements = 2 * contrast_count
contrast_arrays = {
"contrast_estimate": np.arange(1, contrast_elements + 1, dtype="<f8"),
"contrast_se": np.full(contrast_elements, 0.25, dtype="<f8"),
"contrast_stat": np.arange(2, contrast_elements + 2, dtype="<f8"),
"contrast_numdf": np.ones(contrast_elements, dtype="<f8"),
"contrast_ddf": np.full(contrast_elements, 12.0, dtype="<f8"),
"contrast_p": np.full(contrast_elements, 0.05, dtype="<f8"),
}
for suffix, values in contrast_arrays.items():
values.tofile(root / "block_000000.{}.f64.bin".format(suffix))
(root / "block_000000.complete").write_text(
"format\t{}\nblock\t0\ntasks\t2\nbeta_elements\t5\ncov_elements\t9\n"
"fixed_test_elements\t5\n".format(cli.BLOCK_FORMAT),
"fixed_test_elements\t5\ncontrast_test_elements\t{}\n".format(
cli.BLOCK_FORMAT, contrast_elements
),
encoding="utf-8",
)
@@ -107,9 +141,83 @@ class SpectraRemlCliTests(unittest.TestCase):
"--n-extra-covariate-rows",
"--line-search-expansion",
"--zoom-safeguard",
"--n-contrasts",
):
self.assertEqual(command.count(option), 1)
def test_contrast_manifest_and_engine_command(self) -> None:
with tempfile.TemporaryDirectory() as temporary:
root = Path(temporary)
_, manifest = self.fixture(root, with_contrasts=True)
self.assertEqual(manifest["dimensions"]["contrast_count"], 2)
args = cli.build_parser().parse_args(
[
"run", "--manifest", str(root / "manifest.json"),
"--engine", "spectra_reml", "--fixed-effect-test",
"kenward-roger",
]
)
command = cli.build_engine_command(Path("spectra_reml"), manifest, args)
self.assertEqual(command.count("--contrast-matrix"), 1)
self.assertEqual(command.count("--contrast-metadata"), 1)
self.assertEqual(command[command.index("--n-contrasts") + 1], "2")
self.assertEqual(command[command.index("--fixed-effect-test") + 1], "kenward-roger")
def test_new_contract_rejects_missing_contrast_fields(self) -> None:
with tempfile.TemporaryDirectory() as temporary:
root = Path(temporary)
manifest_path, _ = self.fixture(root)
manifest = cli.read_json(manifest_path)
del manifest["dimensions"]["contrast_count"]
del manifest["paths"]["contrast_matrix_f64"]
del manifest["paths"]["contrast_metadata_tsv"]
cli.atomic_json(manifest_path, manifest)
with self.assertRaises(RuntimeError):
cli.validate_manifest(manifest_path)
def test_both_fixed_effect_methods_are_forwarded(self) -> None:
with tempfile.TemporaryDirectory() as temporary:
root = Path(temporary)
_, manifest = self.fixture(root, with_contrasts=True)
for method in ("satterthwaite", "kenward-roger"):
args = cli.build_parser().parse_args(
[
"run", "--manifest", str(root / "manifest.json"),
"--engine", "spectra_reml", "--fixed-effect-test", method,
]
)
command = cli.build_engine_command(Path("spectra_reml"), manifest, args)
self.assertEqual(
command[command.index("--fixed-effect-test") + 1], method
)
def test_contrast_finalize_contract(self) -> None:
with tempfile.TemporaryDirectory() as temporary:
root = Path(temporary)
manifest_path, manifest = self.fixture(root, with_contrasts=True)
engine = root / "spectra_reml"
engine.write_bytes(b"synthetic engine")
args = cli.build_parser().parse_args(
[
"run", "--manifest", str(manifest_path), "--engine",
str(engine), "--resume", "--fixed-effect-test",
"satterthwaite",
]
)
with mock.patch.object(
cli.subprocess, "run", return_value=SimpleNamespace(returncode=0)
):
cli.run_engine(manifest, args)
output = Path(manifest["paths"]["output_directory"])
self.write_blocks(output, contrast_count=2)
result = root / "contrast_results.tsv.gz"
cli.finalize(manifest, result)
with gzip.open(result, "rt", encoding="utf-8", newline="") as handle:
rows = list(csv.DictReader(handle, delimiter="\t"))
self.assertEqual(json.loads(rows[0]["contrast_estimate_json"]), [1.0, 2.0])
self.assertEqual(json.loads(rows[1]["contrast_estimate_json"]), [3.0, 4.0])
self.assertEqual(json.loads(rows[0]["contrast_numerator_df_json"]), [1.0, 1.0])
def test_signature_recovery_and_finalize(self) -> None:
with tempfile.TemporaryDirectory() as temporary:
root = Path(temporary)