Getting started with hsquared

Check the model (no Julia)

hsquared(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat,
  control = hs_control(engine = "validate")
)
#> Validated the v0.1 animal-model contract; no model was fitted (`engine = "validate"`). Julia fit target: `HSquared.fit_animal_model(y, X, Z, Ainv; method = :REML)`. The validated spec is returned invisibly as a named list; assign it to inspect the parsed contract (e.g. `spec$bridge$target`).

model_spec(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat
)
#> <hs_model_spec>
#>   response: weight
#>   family: gaussian(identity)
#>   method: REML
#>   observations: 4; fixed columns: 2; animal IDs: 4
#>   animal term: animal(1 | id, pedigree = ped); group: id
#>   random design: sparse Z with 4 nonzero entries
#>   Ainv: build_in_julia
#>   fitting: not run by `model_spec()`

hsquared(..., control = hs_control(engine = "validate")) confirms the formula and data. model_spec() shows the parsed design and the engine target that would be called. Nothing is fitted.

Fit (needs the Julia engine)

Installing hsquared (pak::pak("itchyshin/hsquared")) gives you the R interface only. Fitting needs a local Julia, the JuliaCall R package, and an HSquared.jl checkout (git clone https://github.com/itchyshin/HSquared.jl). Register the checkout with Sys.setenv(HSQUARED_JULIA_PROJECT = "/path/to/HSquared.jl"), or per call:

fit <- hsquared(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat,
  control = hs_control(
    engine_control = list(julia_project = "/path/to/HSquared.jl")
  )
)

Once the engine is registered, the same call as the validate block attempts a fit. n = 4 may not converge. If fit_diagnostics() says converged is FALSE, do not read heritability() as a result – a near-zero value is not h^2 = 0. heritability() warns when the fit failed. The successful first path on this page is the validate block above.

fit <- hsquared(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat
)

fit_diagnostics(fit)
heritability(fit)          # warns if the fit did not converge
variance_components(fit)
breeding_values(fit)

Without the engine, that default call errors with install guidance rather than silently degrading. Until then, control = hs_control(engine = "validate") is the successful path.

ML is not implemented. REML = FALSE is rejected; use REML = TRUE (the default). That is the same rule for the default fit and for control = hs_control(engine = "validate").

After this works

Read Can I fit and report this? before putting a number in a paper. Point estimates are reportable only on covered routes; intervals are experimental and not coverage-calibrated.

The gryphon worked example is a later teaching dataset (supplied relationship, compare to sommer), not this first example.

formula_status() lists which terms parse and which fit. permanent(), common_env(), and maternal_genetic() are parsed; they fit only with a named target =. The default hsquared() call names that target and stops. See Fitting models. Continue in the same order: Choose a model, Fit, Diagnose, then Report. Read Progression & evidence when you need the history behind a route label rather than a current reporting decision.

Twin-package flow

R formula + data
-> R validation and model specification
-> R bridge payload: y, X, sparse Z, method, IDs, pedigree metadata
-> HSquared.jl sparse engine
-> compact result
-> hsquared_fit object
-> R summaries, extractors, plots, and diagnostics

R owns the friendly interface. Julia owns the fit.

Engine-side walkthrough (pedigree utilities, bridge-shaped specs, experimental fitters): HSquared.jl Quick start. Status before reporting: engine validation status.