This vignette explains the formula-based model specification interface.
Model-fitting and prediction calls are illustrative and deliberately not
evaluated: no fitted results or convergence claims are reported.
Only data preparation and descriptive checks are executed.
We compare six basic model classes on the same DT_example observation
layout, changing one modeling decision at a time:
| Model | Predictor | Family | Random effects |
|---|---|---|---|
| Linear | Ordinary formula | Gaussian identity | None |
| Linear mixed | Ordinary formula | Gaussian identity | Additive |
| Generalized linear | Ordinary formula | Explicit response family/link | None |
| Generalized linear mixed | Ordinary formula | Explicit response family/link | Additive |
| Nonlinear | Expression with declared parameters | Usually Gaussian; others allowed | None |
| Nonlinear mixed | Expression with parameter models | Explicit response family/link | Additive and/or within parameters |
| There is an unavoidable distinction: **a model with a random genetic covariance | |||
| is a mixed model**. The three fixed-only examples cannot also have random | |||
| GxE effects with a relationship matrix. They provide the no-random-effect | |||
| baseline; their mixed counterparts introduce the same diagonal GxE covariance. | |||
Treating Name and Env:Name as fixed effects would be another baseline, but |
|||
| would not impose genetic covariance, shrinkage, or borrowing through relationships. |
DT_example contains 185 observations on 41 potato lines in three environments
in the installed example data. We retain Yield, Weight, Env, Name, and
the original observation layout.
For the linear-versus-nonlinear comparison, use centered Weight as a
descriptive covariate x. It is another measured trait, not an experimentally
assigned dose or a time axis. The association must not be interpreted as a
causal response curve. In particular, do not use it as a predictor in an
application if it is derived from the response or unavailable at prediction
time. Its strong association with Yield makes these examples an interface
illustration, not an independent validation of predictive performance.
library(sommer)
library(Matrix)
data(DT_example, package="enhancer")
DT <- DT_example
DT$Env <- droplevels(factor(DT$Env))
DT$Name <- droplevels(factor(DT$Name))
weight_center <- mean(DT$Weight)
DT$x <- DT$Weight - weight_center
yield_target <- 12
DT$HighYield <- as.integer(DT$Yield >= yield_target)
data.frame(
observations=nrow(DT),
individuals=nlevels(DT$Name),
environments=nlevels(DT$Env),
yield_weight_correlation=cor(DT$Yield, DT$Weight)
)
#> observations individuals environments yield_weight_correlation
#> 1 185 41 3 0.9747398
with(DT, table(Env, HighYield))
#> HighYield
#> Env 0 1
#> CA.2011 6 37
#> CA.2012 51 26
#> CA.2013 47 18
HighYield is an illustrative binary endpoint: whether the observed yield
reaches a prespecified target of 12. It is not a new independent measurement.
The binomial examples model the probability of reaching that target, not
the expected value of Yield. A real analysis would choose a scientifically
meaningful target before examining results.
No Poisson interpretation is inferred merely because Yield is integer-valued:
an integer-valued quality or yield measurement is not automatically an event
count. The design also has blocks; block effects are omitted here to isolate
the interface comparison, not because they are necessarily negligible.
Keep the existing sommer convention: A is the known relationship covariance
shape and Ai is its precision, passed through Gu.
A <- A_example
name_levels <- levels(DT$Name)
stopifnot(
!is.null(rownames(A)),
!is.null(colnames(A)),
all(name_levels %in% rownames(A)),
all(name_levels %in% colnames(A))
)
A <- A[name_levels, name_levels, drop=FALSE]
stopifnot(isSymmetric(A), all(is.finite(A)))
invisible(chol(A))
Ai <- solve(A)
Ai <- as(as(as(Ai, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(Ai, "inverse") <- TRUE
data.frame(
min_diagonal=min(diag(A)),
max_diagonal=max(diag(A)),
max_absolute_off_diagonal=max(abs(A[row(A) != col(A)]))
)
#> min_diagonal max_diagonal max_absolute_off_diagonal
#> 1 1 1 0
In the installed example, A_example is an identity matrix. It illustrates
the relationship-matrix interface, but does not demonstrate borrowing between
related individuals. To demonstrate that behavior, replace A with a genuine
positive-definite pedigree or genomic relationship matrix with matching Name
labels, then rebuild Ai. We do not invent relationships for these potato lines
or silently regularize a singular matrix.
In each mixed example, the covariance declaration is:
~ vsm(dsm(Env), ism(Name), Gu=Ai)
For an environment-major ordering, its statistical meaning is
$$ b \sim N(0,D\otimes A), \qquad D=\operatorname{diag}(\tau_1^2,\ldots,\tau_E^2). $$
Equivalently,
\begin{cases} \tau_e^2 A_{gg'}, & e=e',\\ 0, & e\ne e'. \end{cases} $$
Thus environments have different genetic variances and zero genetic
cross-environment covariance. Within an environment, relationships among
individuals follow A. sommer’s normalized dsm() shape and product-level
vsm() scale together represent these environment-specific variances;
dsm() alone is not a list of absolute variances.
This is the diagonal multienvironment genetic model. We do not also add a
shared Name main effect in the primary examples. Adding
vsm(ism(Name), Gu=Ai) would contribute covariance across environments, so the
total genetic covariance would no longer be diagonal across environments.
The interface separates these responsibilities:
| Component | Meaning |
|---|---|
| Main formula | The complete link-scale predictor, including nonlinear expressions |
parameters |
Regression and random-effect models for named nonlinear parameters |
curves |
Named parametric curves referenced by curve_term() |
family, random, rcov |
Observation distribution and covariance declarations |
method, control, start |
Statistical fitting method, numerical settings, initialization |
Unless stated otherwise, every main formula below describes
$$ g(\mu_i)=\eta_i,\qquad \mu_i=E(Y_i\mid\text{random effects}). $$
An identity link is still a link; it simply makes \(\eta_i=\mu_i\).
Parameter links are separate: parameter(link="log") means that the
parameter, not the response mean, is exponentiated from its working predictor.
The public mmes() arguments include parameters=list(),
curves=list(), start=NULL, and method=mmes_method().
Ordinary linear models need neither parameter nor curve declarations.
Gaussian identity-link linear models, with or without random effects, use
the established Gaussian engines and retain their numerical controls and
variance-component output. For example,
representation=list(fixedRank="drop") permits a reduced fixed-effect basis
when the ordinary fixed design is rank deficient.
The method defaults are
mmes_method(objective="reml", integration="auto", approximation="auto").
Gaussian linear mixed models therefore retain REML by default. For generalized
mixed models, PQL uses a REML working fit; this is not exact GLMM REML.
integration="auto" selects an eligible route rather than promising exact
integration for every model. An unavailable objective, approximation, or
integration combination is an error, not permission to change the requested
model silently.
For comparability, the examples below explicitly request ML. This does not change the REML/PQL defaults or imply that a generic REML likelihood exists for all generalized and nonlinear models. Fixed-only generalized and nonlinear likelihood calls require explicit ML; they reject an ineligible REML request rather than silently changing it.
library(sommer)
exact_ml <- mmes_method(objective="ml", integration="exact")
laplace_ml <- mmes_method(objective="ml", integration="laplace")
fit_control <- mmes_control(
optimizer=list(maxit=1000L),
mode=list(maxit=100L, tol=1e-8)
)
integration="exact" means no random-effect integration for fixed-only
models, and exact Gaussian integration where available for mixed models.
mode controls apply when a conditional mode solve is needed.
optimizer controls outer optimization; these iteration budgets are not
evidence of convergence. mmes_control() retains its existing numerical
modules, except for the removed task module, and adds the optimizer and
mode modules shown above. Statistical method selection belongs in method;
initial values belong in start.
The default optimizer$gradient="auto" uses native analytic derivatives
where available. The expression reference route uses numerical outer scores
with symbolic expression derivatives for conditional curvature
and reports that choice in its diagnostics; requesting an unavailable
gradient="analytic" is an error. Numerical differentiation and an exact
Gaussian integral are separate issues.
Conditional mode settings, including mode$maxit, mode$tol, and
mode$stepTol, govern the inner solve, not the outer optimizer.
mode$warmStart is the sole warm-start setting: its default "auto" permits
native GLMM mode reuse but selects deterministic cold starts ("off") for
nonlinear fits. The likelihood module contains only backend, trace, and
blocks; it does not own mode initialization. Coefficient covariance and
other post-fit inference are separate from optimization: fitting does not
silently request an additional coefficient-covariance calculation.
$$ Y_i\sim N(\mu_i,\sigma^2),\qquad \mu_i=\beta_0+\beta_{\operatorname{Env}(i)}+\beta_x x_i. $$
fit_lm <- mmes(
Yield ~ Env + x,
family=gaussian(link="identity"),
data=DT,
method=exact_ml,
control=fit_control
)
This estimates environment means and one common slope. There is no random genetic term, no relationship covariance in the fitted model, and no genetic shrinkage. The formula is linear in its unknown coefficients.
$$ Y_i\mid b\sim N(\mu_i,\sigma^2),\qquad \mu_i=\beta_0+\beta_{\operatorname{Env}(i)}+\beta_x x_i+b_{eg}, \qquad b\sim N(0,D_b\otimes A). $$
fit_lmm <- mmes(
Yield ~ Env + x,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
rcov=~ units,
family=gaussian(link="identity"),
data=DT,
method=exact_ml,
control=fit_control
)
The only statistical change from the linear model is the addition of the environment-specific genetic deviations. They shift the conditional mean vertically. The relationship matrix operates on those deviations.
$$ H_i\sim\operatorname{Bernoulli}(p_i),\qquad \operatorname{logit}(p_i) =\beta_0+\beta_{\operatorname{Env}(i)}+\beta_x x_i, $$
where \(H_i=1\) when Yield >= 12.
fit_glm <- mmes(
HighYield ~ Env + x,
family=binomial(link="logit"),
data=DT,
method=exact_ml,
control=fit_control
)
The probability curve is nonlinear in x after applying the inverse-logit,
but the predictor is linear in its unknown coefficients. This is a GLM,
not a nonlinear regression model. There is still no random covariance.
No artificial random effect is needed to specify this fixed-only model.
$$ H_i\mid b\sim\operatorname{Bernoulli}(p_i),\qquad \operatorname{logit}(p_i) =\beta_0+\beta_{\operatorname{Env}(i)}+\beta_x x_i+b_{eg}, \qquad b\sim N(0,D_b\otimes A). $$
fit_glmm <- mmes(
HighYield ~ Env + x,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
family=binomial(link="logit"),
data=DT,
method=laplace_ml,
control=fit_control
)
The covariance formula is unchanged from the Gaussian mixed model, but its variance is now on the log-odds scale, not the yield scale. It does not imply a diagonal covariance matrix for observed binary outcomes.
The binomial likelihood determines conditional observation variance
\(p_i(1-p_i)\); an independent Gaussian rcov=~units error is not added to it.
Laplace approximates integration over the Gaussian genetic effects.
Replace the linear contribution \(\beta_x x_i\) with a smooth saturating association \(a\tanh(kx_i)\):
$$ Y_i\sim N(\mu_i,\sigma^2),\qquad \mu_i=\beta_0+\beta_{\operatorname{Env}(i)} +a\tanh(kx_i),\qquad k>0. $$
Unlike a growth curve requiring nonnegative time, tanh() accepts the
negative and positive values of centered Weight. It is used here as a
simple shape comparison, not as a biological model of potato development.
fit_nlm <- mmes(
Yield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, random=NULL, link="identity"),
k=parameter(fixed=~ 1, random=NULL, link="log")
),
family=gaussian(link="identity"),
data=DT,
start=c(a=3, k=1),
method=exact_ml,
control=fit_control
)
Declared symbols a and k are nonlinear parameters; their enclosing
expression contributes to the predictor, while Env remains an ordinary
fixed effect. Unresolved or ambiguous symbols are rejected rather than
guessed. Use the supported expression grammar, not arbitrary R code, in
nonlinear formula contributions.
Here a is unconstrained and \(\log(k)=\alpha_k\). Restricting \(k>0\) avoids the
equivalent sign reversal \((a,k)\leftrightarrow(-a,-k)\). The start supplies
natural-scale parameter values, so k=1 means a working intercept of zero.
The same public natural-scale convention applies to curve starts, even when
their internal coordinates are logarithmic.
The model is nonlinear in the unknown \(k\), even though observations have an identity link. It has no random effects and needs no artificial random term.
For small \(kx\), \(a\tanh(kx)\approx akx\): nearly linear data can identify the product more strongly than its individual factors. A pleasant formula does not remove weak identification, local optima, or the need for sensible starts.
$$ \mu_i=\beta_0+\beta_{\operatorname{Env}(i)} +a\tanh(kx_i)+b_{eg},\qquad b\sim N(0,D_b\otimes A). $$
fit_nlmm_additive <- mmes(
Yield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(fixed=~ 1, link="log")
),
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
rcov=~ units,
family=gaussian(link="identity"),
data=DT,
start=c(a=3, k=1),
method=exact_ml,
control=fit_control
)
This is nonlinear in its fixed parameters, but Gaussian and linear in its
random effects. Integration over b remains exact Gaussian integration;
estimating a and k still needs nonlinear optimization.
$$ \begin{aligned} \mu_i&=\beta_0+\beta_{\operatorname{Env}(i)} +a\tanh(k_{eg}x_i),\\ \log(k_{eg})&=\alpha_k+u_{eg},\qquad u\sim N(0,D_k\otimes A). \end{aligned} $$
fit_nlmm_rate <- mmes(
Yield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(
fixed=~ 1,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
link="log"
)
),
rcov=~ units,
family=gaussian(link="identity"),
data=DT,
start=c(a=3, k=1),
method=laplace_ml,
control=fit_control
)
Moving the same covariance declaration into k changes its meaning:
There is no top-level random term here: the parameter random effect already
makes the model mixed.
fit_nlmm_both <- mmes(
Yield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(
fixed=~ 1,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
link="log"
)
),
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
rcov=~ units,
family=gaussian(link="identity"),
data=DT,
start=c(a=3, k=1),
method=laplace_ml,
control=fit_control
)
These declarations specify independent random blocks: the top-level additive effect and the log-rate effect have separate covariance parameters, even though they use the same formula and relationship matrix. Repeated formulas do not silently share effects or create cross-covariance. Repeating a formula does not declare a joint additive/parameter covariance.
These are syntax examples, not evidence that this small, unbalanced dataset supports estimating both genetic mechanisms. Shape variation requires suitable within-group covariate information; sparse replication may make this model weakly identified.
To ask whether environments change the common curvature, rather than only the vertical level, change the parameter model:
k=parameter(
fixed=~ Env,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
link="log"
)
This specifies \(\log(k_{eg})=\alpha_0+\alpha_e+u_{eg}\).
Keeping Env in the main formula also permits environment-level shifts;
removing it there imposes a different scientific restriction.
Natural parameter starts such as start=c(a=3, k=1) initialize parameter
intercepts and set their other regression coefficients to zero. Thus k=1
also initializes a constant rate when fixed=~ Env; it does not supply
different environment-specific rates.
For explicit regression starts, use start=list(coefficients=...) with a
complete named vector of mean coefficients on their working/link scale.
Names have the form fixed:(Intercept), fixed:Env2, a:(Intercept),
k:(Intercept), or k:Env2, following the actual model-matrix column names.
For example, k:(Intercept)=0 means the natural rate starts at one, and
k:Env2=log(2) doubles that rate in environment 2. Include every mean
coefficient, not just the changed contrast. These are coefficient starts,
not random-effect modes or covariance starts.
The six familiar labels are not six separate engines. Nonlinear predictors can also have non-Gaussian observations.
$$ \operatorname{logit}(p_i) =\beta_0+\beta_{\operatorname{Env}(i)}+a\tanh(kx_i). $$
fit_gnlm <- mmes(
HighYield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(fixed=~ 1, link="log")
),
family=binomial(link="logit"),
data=DT,
start=c(a=3, k=1),
method=exact_ml,
control=fit_control
)
This is not an ordinary GLM: the link predictor is nonlinear in unknown k.
There are no random effects to integrate; exact_ml does not mean that
nonlinear optimization has a closed-form solution or a guaranteed global optimum.
$$ \begin{aligned} \operatorname{logit}(p_i) &=\beta_0+\beta_{\operatorname{Env}(i)} +a\tanh(k_{eg}x_i),\\ \log(k_{eg})&=\alpha_k+u_{eg},\qquad u\sim N(0,D_k\otimes A). \end{aligned} $$
fit_gnlmm <- mmes(
HighYield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(
fixed=~ 1,
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
link="log"
)
),
family=binomial(link="logit"),
data=DT,
start=c(a=3, k=1),
method=laplace_ml,
control=fit_control
)
Two links now coexist without being confused:
The family has not changed its responsibility. The nonlinear expression and its parameter models have changed how the link predictor is constructed.
If a continuous positive-response model were scientifically justified,
Yield could instead be modeled with a Gamma-log likelihood:
fit_gamma_nlmm <- mmes(
Yield ~ Env + a * tanh(k * x),
parameters=list(
a=parameter(fixed=~ 1, link="identity"),
k=parameter(fixed=~ 1, link="log")
),
random=~ vsm(dsm(Env), ism(Name), Gu=Ai),
family=Gamma(link="log"),
data=DT,
start=c(a=3, k=1),
method=laplace_ml,
control=fit_control
)
This is a distributional what-if, not a recommendation for the recorded integer-valued quality measurements. It assumes
$$ \log(\mu_i)=\beta_0+\beta_{\operatorname{Env}(i)} +a\tanh(kx_i)+b_{eg}. $$
Consequently, \(\mu_i=\exp\{\beta_0+\beta_e+a\tanh(kx_i)+b_{eg}\}\).
The displayed nonlinear expression is not directly the response mean.
The same Gamma family could accompany the ordinary Yield ~ Env + x
predictor, with or without the additive random term.
Let \(B_i=\beta_0+\beta_{\operatorname{Env}(i)}\), \(b\sim N(0,D_b\otimes A)\), and \(u\sim N(0,D_k\otimes A)\).
| Example | Conditional predictor | Genetic covariance acts on | Integration |
|---|---|---|---|
| Linear | \(\mu_i=B_i+\beta_xx_i\) | None | None |
| Linear mixed | \(\mu_i=B_i+\beta_xx_i+b_{eg}\) | Additive yield deviation | Exact Gaussian |
| Generalized linear | \(\operatorname{logit}(p_i)=B_i+\beta_xx_i\) | None | None |
| Generalized linear mixed | \(\operatorname{logit}(p_i)=B_i+\beta_xx_i+b_{eg}\) | Additive log-odds deviation | Laplace |
| Nonlinear | \(\mu_i=B_i+a\tanh(kx_i)\) | None | None |
| Nonlinear mixed, additive | \(\mu_i=B_i+a\tanh(kx_i)+b_{eg}\) | Additive yield deviation | Exact Gaussian |
| Nonlinear mixed, shape | \(\mu_i=B_i+a\tanh(\exp(\alpha_k+u_{eg})x_i)\) | Log-rate deviation | Laplace |
| Generalized nonlinear mixed | \(\operatorname{logit}(p_i)=B_i+a\tanh(\exp(\alpha_k+u_{eg})x_i)\) | Log-rate deviation | Laplace |
The repeated dsm(Env) and relationship matrix preserve the covariance
pattern, not the units or interpretation of its estimated variance.
Also, diagonal latent GxE covariance need not translate to a simple
Kronecker covariance for nonlinear response means or observed responses.
Do not compare Gaussian Yield likelihoods or AIC values directly with
binomial HighYield likelihoods: they describe different responses.
Gaussian and continuous Gamma likelihoods should likewise not be used to
settle how an intrinsically discrete measurement was generated.
Within a common response and justified likelihood, model comparisons still
need convergence, identification, and appropriate validation checks.
For the expression reference route, independent residuals are specified by
rcov=NULL or rcov=~ units; structured residual covariance is not supported
there. The Gamma-log example estimates a dispersion parameter rather than
adding an independent Gaussian error to the Gamma response.
Name the prediction estimand explicitly:
predict(fit_gnlmm, newdata=DT, type="response", target="conditional")
predict(fit_gnlmm, newdata=DT, type="response", target="zero")
predict(fit_gnlmm, newdata=DT, type="response", target="population")
conditional: evaluate at fitted random-effect modes.zero: set all random deviations to zero, including deviations inside k.population: integrate over the fitted joint Gaussian random-effect prior.These are distinct requests, not aliases. A target that is unavailable for a particular fitting route must raise an explicit error; it must not return a conditional or zero-deviation prediction labeled as population-averaged. In particular, the target vocabulary does not promise population integration for every formula, likelihood, and covariance combination.
Native nonlinear and native Laplace fits support the explicit target
interface. The established Gaussian and PQL engines retain their
D/classify prediction interface and reject explicit target requests.
For expression-reference population prediction, joint predictor quadrature
is admitted by dimension, cost, accuracy and memory; there is no fixed
two-record cutoff. Unsupported integration requests error.
Choose prediction arguments for the fitted route rather than assuming that
the calls above apply unchanged to every mmes() result.
These targets generally differ for nonidentity links and random nonlinear
parameters. For example, averaging rates before evaluating tanh() is not
the same as integrating the resulting response means. Population prediction
is not an uncertainty interval or integration over parameter-estimation error.
For new observations, x must use the saved training weight_center, not a
center recomputed from the prediction batch.
| Contract | Native curve | Reference expressions |
|---|---|---|
| Random coordinates | Independent or explicit joint curve-parameter blocks; proper latent RR maps | Independent declared parameter records and additive records; RR rejected |
| Conditional Hessian | Curve derivatives and native family curvature | Symbolic expression chain rules and full native family curvature |
| Outer score | Analytic | Numerical differences |
| Integration | Exact Gaussian additive or conditional-mode Laplace | Exact Gaussian additive, no-integral fixed-only, or conditional-mode Laplace |
| Memory | Sparse native admission plus bounded dense post-fit work | Dense: 30 declarations, 60 mean coefficients, 512 modes for nonexact integration |
| Numerical controls | Native eligibility checks | No nondefault representation/solver/update controls, weights, custom contrasts or vcc |
familyParameters controls fixed/estimated dispersion, NB2 shape, ZIP
structural-zero predictors and ordered ordinal thresholds using the same
family registry on both routes. A family link and a parameter link act at
different stages and are each applied once.
Legacy scalar expression links outside that registry retain
their reference likelihood via a chain-rule adapter to the native density;
this is not mixed-family link support. Ordinary fixed-only
lm/glm links are unchanged.
Fitting never computes uncertainty. post_nonlinear(), vcov(), confint(),
predict(se.fit=TRUE) and summary(inference=TRUE) explicitly request local
observed-information/delta approximations. They require a reconstructed
likelihood matching the saved fit, stationarity, stable derivatives, positive
identified curvature and interior covariance/nuisance parameters. Conditional
modes are re-certified more tightly for post-fit numerical derivatives;
this does not change the saved fitting controls or statistical model. Ordinary
boundary-variance tests are not supported.
vcov(type="conditional") holds marginal parameters fixed;
vcov(type="random") adds mode sensitivity and is an approximate empirical
Bayes PEV, not an exact frequentist PEV.
Specification coefficients remain fixed-regression coefficients on their
declared link scales; coef(type="group") returns natural parameter values
per retained predictor row. coef(type="all") and vcov(type="all") include
covariance/family reporting coordinates. Confidence intervals are local normal
approximations. Prediction intervals add family observation variance and, for
population targets, integrated prior variance; categorical probabilities
support confidence intervals only. Grouped-binomial predictions are
proportions and require trials for new-data prediction intervals.
Binomial simulations also return proportions using the declared integer trials.
Seeded simulations restore the caller’s RNG.
profile_nonlinear() measures operations without an outer refit, and reports
the dense reference path distinctly from native sparse execution.
Ordinary Gaussian models retain representation=list(fixedRank="drop");
this control does not remove unidentified nonlinear coordinates.
This is a breaking replacement of the earlier specification interface, not a
compatibility layer. Put nonlinear scientific expressions in the main formula,
declare named parameter models in parameters, and register native curves
in curves. There is no separate nonlinear= argument. Put statistical
choices in mmes_method(), starts in start, and numerical settings in
mmes_control(); the former task control module is not part of this API.
Old calls must be rewritten rather than relying on compatibility mapping.
The syntax does not remove capability checks: unsupported expressions,
families, covariance structures, methods, or prediction targets are errors.
The general expression route is a dense numerical reference implementation,
currently limited to 60 mean coefficients and 512 random coefficients for
nonexact integration, with
an additional workspace check against execution$denseMemoryMB.
It does not promise the scalability of the specialized native engines;
see ?parameter for supported expressions and route-specific limits.
Use convergence diagnostics, scientifically defensible response models,
identifiability checks, and validation appropriate to the target population.
In particular, syntax alone cannot make the identity A_example demonstrate
relationship borrowing, make a measured Weight association causal, or make
a fixed-only fit contain random GxE.