predict(), fitted(), and residuals() are part of the planned v0.1
fitted-object contract for univariate hsquared_fit objects. They are
univariate-only: the multivariate model (a cbind() response)
fits multiple traits jointly and is intentionally out of v0.1 response-scale
scope, so these methods stop with a scope message pointing to
breeding_values(), genetic_covariance(), and residual_covariance().
Out-of-sample newdata is not implemented; supplying it is an error rather
than a silent in-sample result.