From conversation to a reproducible analysis

Work with an assistant to inspect data, specify a gsuite model, run it and prepare a report.

You bring the scientific question and study knowledge. An assistant can help inspect data, run R code, examine diagnostics and prepare a report using gsuite.

This illustrative conversation uses the single-trait pedigree example; it is not a completed analysis. To run it, the assistant needs access to your files, R and the existing development integration. That interface is not provided by the seven public source packages; greml/gsolve R packaging is deferred. Website rendering does not execute the code.

1. Start with the question and data

You might ask:

I have body-weight records, sex, cohort and herd, plus a pedigree. Help me estimate additive genetic variation and predict animal effects. First inspect the data and explain any issues that affect the analysis.

The assistant checks the column meanings and measurement units with you, then inspects missing values, duplicated records, factor levels and phenotype distributions. It checks that observation IDs match the pedigree, parent IDs are represented consistently, and pedigree relationships are valid. Repeated measurements require particular attention: several records per animal may need a different model from one observation per animal.

For a reproducible demonstration, use the Prepare ordinary data frames section and the data preparation at the start of Single trait in the pedigree tutorial. That supplies pedigree and data, including BW, sex, cohort, herd and id. For your study, use your own data frames instead.

str(data)
summary(data)
colSums(is.na(data))
anyDuplicated(pedigree$id)
setdiff(data$id, pedigree$id)
table(table(data$id))

These are initial checks, not a complete data-validation procedure. The assistant records issues and proposed changes in a data-check summary. Decisions such as excluding records or recoding missing parents belong in the analysis script, with an explanation you can review.

2. Translate the question into a model

You might ask:

Treat sex, cohort and herd as fixed effects. Use the pedigree for additive genetic animal effects. Explain the assumptions and show me the model code.

For the tutorial’s single-trait model:

\[ y = X\beta + Za + e,\qquad \operatorname{Var}(a)=\sigma_a^2 A,\qquad \operatorname{Var}(e)=\sigma_e^2 I. \]

Here \(A\) is the pedigree relationship matrix. The model assumes independent animal and residual effects and a common residual variance. Whether herd should be fixed, and whether additional environmental effects are needed, depends on the study design and your scientific question.

The qgtools-style separation of effects, relationships and covariance components makes each decision explicit:

Model decision	gsuite specification
Fixed effects and animal grouping	Formula
Observations and explanatory variables	data
Relationships between animals	Pedigree kernel
Genetic and residual variance components	gvc() descriptors
Estimate covariance parameters	task = "reml"

library(gsuite)

PED <- gs_pedigree(pedigree, id = "id", sire = "sire", dam = "dam")
vcs <- list(
  animal = gvc(variable = "id", traits = "BW", kernel = PED),
  residual = gvc(variable = "residual", traits = "BW")
)

Starting variances are omitted here: gsuite uses its existing deterministic initializers. The pedigree tutorial also shows the short form, which expresses the same model with fewer arguments. The assistant checks the supported model combinations before suggesting an extension; a general specification does not imply that every model is implemented.

3. Run and examine the fit

You might ask:

Run this model, retain the script and fit, and check convergence before interpreting the estimates. Explain any failures rather than silently changing the model.

fit <- gfit(
  formula = BW ~ sex + cohort + herd + (1 | id),
  data = data,
  vc = vcs,
  task = "reml"
)

summary(fit)
fit_diagnostics <- diagnostics(fit)
fit_diagnostics$status
fit_diagnostics$converged
fit_diagnostics$initialization
stopifnot(isTRUE(fit_diagnostics$status$ok),
  isTRUE(fit_diagnostics$converged))
coef(fit)
head(gs_ranef(fit)$animal)

This data-frame route defaults to CHOLMOD selected inversion and requires a build with CHOLMOD support. The check above stops the example before prediction plots and reporting if status or convergence is unsuccessful.

After checking status and convergence, examine estimates and residual patterns for unusual values or possible confounding.

plot(fitted(fit), residuals(fit),
     xlab = "Fitted value", ylab = "Residual")
abline(h = 0, lty = 2)

Record failures and unavailable diagnostics alongside the results.

4. Turn the analysis into a report

You might ask:

Create a Quarto report with the data checks, model and assumptions, convergence diagnostics, variance estimates, animal predictions and residual plots. Explain what the results support and what remains uncertain.

The report uses the actual saved results. It identifies the data used, exclusions, software versions and analysis settings, and links conclusions to tables and figures. It separates predicted animal effects from observed phenotypes and avoids presenting a failed fit as a completed analysis.

The working directory can remain small:

File	Purpose
analysis.R	Data preparation, checks and model fitting
results.rds	Extracted estimates, predictions and diagnostics for reporting
report.qmd	Methods, results, figures and interpretation
session-info.txt	R and package versions

Save the reporting values rather than relying on an in-memory fit to survive between R sessions:

results <- list(
  summary = summary(fit),
  diagnostics = diagnostics(fit),
  fixed_effects = coef(fit),
  animal_effects = gs_ranef(fit)$animal,
  record_mapping = fit$record_mapping,
  fitted = fitted(fit),
  residuals = residuals(fit)
)
saveRDS(results, "results.rds")
writeLines(capture.output(sessionInfo()), "session-info.txt")

The conversation helps develop the analysis. The scripts, extracted results and report preserve the work so you or a colleague can inspect and reproduce it without reconstructing the conversation.

Continue with more complex models

The same collaboration applies to multitrait and longitudinal analyses: discuss the question, inspect the record structure, specify effects and covariances, run a supported model and assess the results. gsuite remains an R package with explicit interfaces; no conversational service is required inside the package.