---
title: "Getting started with hsquared"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with hsquared}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

> **Experimental 0.9.0 release — not production.** CRAN submission, acceptance,
> and public availability are separate states. Version tracks *covered*
> capability, not surface area. `public_covered_count` is **7** (multivariate
> + explicit 0.7 genomic GREML). Fitting needs a local Julia and an
> `HSquared.jl` checkout; it is not registered in Julia General. Do not use
> install-by-name (`Pkg.add("HSquared")`). Read
> [Can I fit and report this?](https://itchyshin.github.io/hsquared/articles/current-limits.html)
> before putting a number in a paper.

```{r f1-animal-path, echo = FALSE, fig.alt = "Flowchart of the default univariate animal model: pedigree and phenotype enter the hsquared formula, Julia REML returns variance components, heritability, and breeding values. Genomic and multivariate paths are not shown.", out.width = "100%"}
knitr::include_graphics("../man/figures/animal-model-path.svg")
```

<p class="hs-banner" role="note"><strong>Default univariate path shown.</strong> Genomic and QTL routes are not this arrow. Multivariate <code>cbind()</code> is a separate path — read <a href="https://itchyshin.github.io/hsquared/articles/current-limits.html">Can I fit and report this?</a> before reporting it.</p>

`hsquared` is the R interface for a Julia-backed animal model. The first
working path is one Gaussian trait, one `animal()` term, REML. This page uses
a four-animal pedigree so the first block runs without Julia. **n = 4 is a
syntax demo, not a number for a paper.**

```{r golden-data}
library(hsquared)

ped <- data.frame(
  id   = c("sire", "dam", "off1", "off2"),
  sire = c(NA, NA, "sire", "sire"),
  dam  = c(NA, NA, "dam", "dam")
)
dat <- data.frame(
  id     = c("sire", "dam", "off1", "off2"),
  sex    = c("m", "f", "m", "f"),
  weight = c(42, 38, 40, 37)
)
```

## Check the model (no Julia)

```{r golden-validate}
hsquared(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat,
  control = hs_control(engine = "validate")
)

model_spec(
  weight ~ sex + animal(1 | id, pedigree = ped),
  data = dat
)
```

`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](https://julialang.org/downloads/),
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:

```{r golden-fit-path, eval = FALSE}
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.

```{r golden-fit, eval = FALSE}
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?](https://itchyshin.github.io/hsquared/articles/current-limits.html)
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](https://itchyshin.github.io/hsquared/articles/gryphon-worked-example.html)
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](https://itchyshin.github.io/hsquared/articles/fitting-models.html).
Continue in the same order: [Choose a model](https://itchyshin.github.io/hsquared/articles/current-limits.html),
[Fit](https://itchyshin.github.io/hsquared/articles/fitting-models.html), [Diagnose](https://itchyshin.github.io/hsquared/articles/visualizing-models.html), then
[Report](https://itchyshin.github.io/hsquared/articles/current-limits.html). Read
[Progression & evidence](https://itchyshin.github.io/hsquared/articles/progression-evidence.html) when you need the
history behind a route label rather than a current reporting decision.

## Twin-package flow

```text
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](https://itchyshin.github.io/HSquared.jl/dev/quickstart.html).
Status before reporting:
[engine validation status](https://itchyshin.github.io/HSquared.jl/dev/validation-status.html).