
Prior distributions in lotri
Matthew Fidler
2026-08-15
Source:vignettes/lotri-priors.Rmd
lotri-priors.RmdWhy priors are in lotri
A lotri({}) block is the ini({}) block of a
nlmixr2 model. Until now it could only say what a parameter
is – an estimate, optional bounds, a fix flag, a
label – but not what is believed about it beforehand.
Bayesian estimation methods need that missing piece. This adds a way to attach a prior distribution to any parameter in the block.
lotri deliberately stops at parsing, validating and
storing the prior. It does not evaluate the density and it does not
generate any ‘Stan’ code; that is the job of the package doing the
estimation. What lotri guarantees is that by the time the
estimation code sees a prior, the distribution exists, has the right
number of arguments, and does not contradict the parameter’s bounds.
The syntax
A prior is given with prior(name) ~ dist(...):
m <- lotri({
tka <- 0.45
tcl <- c(0, 1, 10)
prior(tka) ~ dnorm(0, 10)
prior(tcl) ~ dlnorm(1, 0.5)
})
lotriEst(m)
#> name lower est upper fix label backTransform prior
#> 1 tka -Inf 0.45 Inf FALSE <NA> <NA> dnorm(0, 10)
#> 2 tcl 0 1.00 10 FALSE <NA> <NA> dlnorm(1, 0.5)Because the statement names the parameter it applies to, prior lines are order independent. These two blocks are the same:
m1 <- lotri({
tka <- 0.45
prior(tka) ~ dnorm(0, 10)
})
m2 <- lotri({
prior(tka) ~ dnorm(0, 10)
tka <- 0.45
})
identical(m1, m2)
#> [1] TRUEThat matters in practice: you can keep all the priors together at the bottom of a long block instead of interleaving them with the estimates.
prior() is the general form and works for any
distribution. Normal priors are common enough to have a shorthand as
well, which is covered below in the normal prior shorthand.
What a prior can be attached to
A whole covariance block
Correlation and covariance matrices have their own distributions, and those apply to the block rather than to any one element. Give every name in the block:
m <- lotri({
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
prior(eta.cl, eta.v) ~ lkjCorr(2)
})
attr(m, "lotriPriors")
#> [1] "lkjCorr(2)" NAThe prior is stored on the first diagonal element of the block it belongs to, which is why only the first entry is filled in above.
Degrees of freedom for an omega block
The scale matrix of the Wishart family is optional, because the block it is put on already is that matrix. So an inverse Wishart prior on an omega block is written by giving only its degrees of freedom:
m <- lotri({
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
eta.ka ~ 0.3
prior(eta.cl, eta.v) ~ invWishart(4)
prior(eta.ka) ~ invWishart(2)
})
as.data.frame(m)[, c("name", "est", "prior")]
#> name est prior
#> 1 eta.cl 0.10 invWishart(4)
#> 2 (eta.cl,eta.v) 0.01 <NA>
#> 3 eta.v 0.20 <NA>
#> 4 eta.ka 0.30 invWishart(2)This is the pair that a NONMEM NWPRI model writes as
$OMEGAP and $OMEGAPD: the omega values you
already gave are the prior scale matrix, and the number here is the
degrees of freedom.
It works on a 1x1 block as well, since an inverse Wishart of dimension one is an inverse gamma.
When every block shares the same degrees of freedom, a one-sided
~ sets them all at once instead of naming each block:
m <- lotri({
eta.cl + eta.v ~ c(0.3,
0.01, 0.1)
eta.ka ~ 0.5
~invWishart(4)
})
as.data.frame(m)[, c("name", "est", "prior")]
#> name est prior
#> 1 eta.cl 0.30 invWishart(4)
#> 2 (eta.cl,eta.v) 0.01 <NA>
#> 3 eta.v 0.10 <NA>
#> 4 eta.ka 0.50 invWishart(4)Both blocks got it – the 2x2 and the 1x1. Each is still checked individually, so a degrees of freedom that is proper for one block and not another is caught. Naming a block as well is a duplicate rather than an override:
lotri({
eta.cl + eta.v ~ c(0.3,
0.01, 0.1)
~invWishart(4)
prior(eta.cl, eta.v) ~ invWishart(5)
})
#> Error:
#> ! more than one prior given for 'eta.cl, eta.v'Give the scale matrix explicitly only when it differs from the estimates:
lotriEst(lotri({
e1 + e2 ~ c(1,
0.1, 1)
prior(e1, e2) ~ invWishart(4, lotri(e1 + e2 ~ c(2,
0.5, 2)))
}))
#> NULL
attr(lotri({
e1 + e2 ~ c(1,
0.1, 1)
prior(e1, e2) ~ invWishart(4, lotri(e1 + e2 ~ c(2,
0.5, 2)))
}), "lotriPriors")
#> [1] "invWishart(4, lotri(e1 + e2 ~ c(2, 0.5, 2)))"
#> [2] NABecause both the degrees of freedom and the size of the block are known, an improper prior is caught:
A prior on the omega values themselves
The two NONMEM prior flavours want different things from an omega. An
NWPRI model gives it degrees of freedom, which is the
invWishart() above. A TNPRI model instead puts
a normal prior on the omega elements, jointly with the thetas,
so those elements need names of their own.
Prepending om. to a between subject variability names
its omega element:
m <- lotri({
eta.cl ~ 0.3
eta.v ~ 0.1
om.eta.cl ~ 0.01
om.eta.v ~ 0.04
})
attr(m, "lotriPriors")
#> [1] "dnorm(0.3, 0.1)" "dnorm(0.1, 0.2)"om.eta.cl ~ 0.01 reads exactly like the shorthand for a
population estimate: a normal prior with a variance of 0.01, centered on
the omega value the model already gives. The omega itself is untouched –
only the prior was added:
diag(m)
#> eta.cl eta.v
#> 0.3 0.1Correlated omega priors work the same way, including the per-row line form:
m2 <- lotri({
eta.cl + eta.v ~ c(0.3,
0.01, 0.1)
om.eta.cl ~ 0.01
om.eta.v ~ c(0.001, 0.02)
})
attr(m2, "lotriPriors")[1]
#> [1] "multiNormal(c(0.3, 0.1), lotri(om.eta.cl + om.eta.v ~ c(0.01, 0.001, 0.02)))"An om. name has to match a real between subject
variability; it never quietly creates one:
lotri({
eta.ka ~ 0.3
om.eta.nope ~ 0.1
})
#> Error:
#> ! prior given for unknown parameter(s): 'om.eta.nope'Naming the eta directly means the same thing, so
prior(om.eta.cl) and prior(eta.cl) are
interchangeable. The om. spelling exists so the shorthand
has a name to put on the left of the ~, since
eta.cl ~ ... already means the omega value itself.
The two are alternatives, not additions, so a model that gives an omega both degrees of freedom and a normal prior is rejected:
One joint block over thetas and omega elements
A TNPRI variance matrix does not stop at the omegas – it
covers the thetas and the omega elements together, with
covariances between them. A block may therefore name both:
m <- lotri({
tcl <- 1
eta.cl ~ 0.3
tcl + om.eta.cl ~ c(0.01,
0.002, 0.005)
})
lotriEst(m)$prior
#> [1] "multiNormal(c(1, 0.3), lotri(tcl + om.eta.cl ~ c(0.01, 0.002, 0.005)))"The mean is what the model already says: the estimate for
tcl, the omega value for om.eta.cl. The prior
is stored once, on the first name of the block, because the block spans
two places – the estimates and the omega – and there is no one row that
owns it. The covariance keeps every name, so the members are recovered
from it rather than from where it is stored:
Naming the block with prior() means the same thing:
lotriEst(lotri({
tcl <- 1
eta.cl ~ 0.3
prior(tcl, om.eta.cl) ~ multiNormal(c(1, 0.3),
lotri(tcl + om.eta.cl ~ c(0.01,
0.002, 0.005)))
}))$prior
#> [1] "multiNormal(c(1, 0.3), lotri(tcl + om.eta.cl ~ c(0.01, 0.002, 0.005)))"Unlike a prior on the omega itself, a joint block does not have to be
one covariance block – a TNPRI matrix covers whichever
elements it likes:
One block at a time
The names have to be exactly one block. Two unrelated diagonal elements are two 1x1 blocks, not one 2x2 block, so this is an error:
lotri({
eta.a ~ 1
eta.b ~ 1
prior(eta.a, eta.b) ~ lkjCorr(2)
})
#> Error:
#> ! 'eta.a, eta.b' is not a single covariance block, so it cannot share a priorThe normal prior shorthand
Normal priors are by far the most common, so they have a shorthand
that reuses the matrix syntax you already know. Putting a
population estimate on the left of a ~
gives it a normal prior:
lotriEst(lotri({
tka <- 1
tka ~ 4
}))
#> name lower est upper fix label backTransform prior
#> 1 tka -Inf 1 Inf FALSE <NA> <NA> dnorm(1, 2)The number on the right is a variance, so
tka ~ 4 is a normal prior with a standard deviation of 2.
The mean is the <- estimate, so the prior is centered on
what the model already says the parameter is – the pair a NONMEM
NWPRI model writes as $THETAP and
$THETAPV. Write prior(tka) ~ dnorm(mu, sd)
when the prior is not centered there.
This is unambiguous because a name cannot be both an estimate and an eta – that combination used to be an error.
More than one parameter
The full matrix syntax works, so a covariance between the priors is written the same way a covariance between etas is:
m <- lotri({
tka <- 1
tcl <- 3
tv <- 4
tcl + tv ~ c(1,
0.01, 1)
})
lotriEst(m)$prior
#> [1] NA
#> [2] "multiNormal(c(3, 4), lotri(tcl + tv ~ c(1, 0.01, 1)))"
#> [3] "multiNormal(c(3, 4), lotri(tcl + tv ~ c(1, 0.01, 1)))"That is a multivariate normal whose mean vector is the estimates. The
covariance is kept as the lotri expression that built it,
which is valid R and round trips exactly.
The matrix means exactly what it means for etas: the off-diagonal is a covariance, not a correlation. A prior block and an eta block written the same way give the same matrix:
.eta <- lotri({ a + b ~ c(1,
0.5, 2) })
.prior <- lotri({
a <- 1
b <- 2
a + b ~ c(1,
0.5, 2)
})
## pull the covariance back out of the stored prior
identical(unname(as.matrix(eval(str2lang(lotriEst(.prior)$prior[1])[[3]]))),
unname(as.matrix(.eta)))
#> [1] TRUEThe per-row line form builds up the block, exactly as it does for etas:
When the parameters are uncorrelated the result is simply independent normal priors, since that is what an MVN with a diagonal covariance is:
Transformations
sd(), var(), cor(),
cov() and chol() all work here too, so the
prior can be written whichever way is most natural:
lotriEst(lotri({
tcl <- 3
tv <- 4
tcl + tv ~ sd(2,
0.5, 3)
}))$prior
#> [1] "multiNormal(c(3, 4), lotri(tcl + tv ~ c(4, 0.5, 9)))"
#> [2] "multiNormal(c(3, 4), lotri(tcl + tv ~ c(4, 0.5, 9)))"sd(2, ..., 3) gives variances of 4 and 9, as the stored
prior shows.
Zero variance
A zero variance is a point mass rather than a prior, so it is rejected rather than quietly accepted:
lotri({
tka <- 1
tka ~ 0
})
#> Error:
#> ! a normal prior on 'tka' cannot have zero variance; did you mean 'fix()'?Distribution names
There are three spellings of every distribution, and all of them are accepted on input:
- the R name, when R parameterizes the distribution
the same way ‘Stan’ does –
dnorm(),dlnorm(),dgamma(),dbeta() - the camelCase name, which is the ‘Stan’ name
written the way the rest of this package is written –
invWishart(),lkjCorr(),studentT() - the ‘Stan’ name itself –
inv_wishart(),lkj_corr(),student_t()
The canonical one – what gets stored and printed back – is the R name when there is a faithful one, and the camelCase name otherwise. Whichever you write, you get the same thing:
a <- lotri({ tka <- 0.45; prior(tka) ~ dnorm(0, 10) })
b <- lotri({ tka <- 0.45; prior(tka) ~ normal(0, 10) })
identical(a, b)
#> [1] TRUE
lotriEst(a)$prior
#> [1] "dnorm(0, 10)"
.camel <- lotri({ e1 ~ 1; prior(e1) ~ invWishart(2) })
.stan <- lotri({ e1 ~ 1; prior(e1) ~ inv_wishart(2) })
identical(.camel, .stan)
#> [1] TRUE
attr(.camel, "lotriPriors")
#> [1] "invWishart(2)"Arguments may be positional or named, in any order:
Why dt() is not studentT()
The one place where an R name is refused is
dt(). R’s dt(x, df, ncp) is the standardized
(or noncentral) t, while studentT(nu, mu, sigma) (the
‘Stan’ student_t) is a location-scale t. They are different
distributions, so aliasing them would silently change the model:
lotri({
tka <- 0.45
prior(tka) ~ dt(3)
})
#> lotri syntax error:
#> =================================================================================
#> :001: tka <- 0.45
#> lotri error:
#> bad matrix expression: 'tka ~ dt(3)'
#> matrix expression should be 'name ~ c(lower-tri)'
#> :002: prior(tka) ~ dt(3)
#> =================================================================================
#> Error:
#> ! lotri syntax errors aboveUse studentT() with its own parameterization
instead:
The supported distributions
lotriPriorDists() returns the whole table, including the
‘Stan’ name for each distribution, which is what a package generating
‘Stan’ code needs:
d <- lotriPriorDists()
nrow(d)
#> [1] 47
head(d, 10)
#> rName stanName camelName name
#> 1 dnorm normal normal dnorm
#> 2 <NA> std_normal stdNormal stdNormal
#> 3 <NA> exp_mod_normal expModNormal expModNormal
#> 4 <NA> skew_normal skewNormal skewNormal
#> 5 <NA> student_t studentT studentT
#> 6 dcauchy cauchy cauchy dcauchy
#> 7 <NA> double_exponential doubleExponential doubleExponential
#> 8 dlogis logistic logistic dlogis
#> 9 <NA> gumbel gumbel gumbel
#> 10 <NA> skew_double_exponential skewDoubleExponential skewDoubleExponential
#> parNames nPar nReq support kind
#> 1 mean,sd 2 2 real univariate
#> 2 0 0 real univariate
#> 3 mu,sigma,lambda 3 3 real univariate
#> 4 xi,omega,alpha 3 3 real univariate
#> 5 nu,mu,sigma 3 3 real univariate
#> 6 location,scale 2 2 real univariate
#> 7 mu,sigma 2 2 real univariate
#> 8 location,scale 2 2 real univariate
#> 9 mu,beta 2 2 real univariate
#> 10 mu,sigma,tau 3 3 real univariateThe kind column says what a distribution may be attached
to – univariate for a single parameter, matrix
and multivariate for a covariance block:
subset(d, kind == "matrix", select=c(name, stanName, parNames))
#> name stanName parNames
#> 29 lkjCorr lkj_corr eta
#> 30 lkjCorrCholesky lkj_corr_cholesky eta
#> 31 wishart wishart nu,Sigma
#> 32 invWishart inv_wishart nu,Sigma
#> 33 wishartCholesky wishart_cholesky nu,L_S
#> 34 invWishartCholesky inv_wishart_cholesky nu,L_SBounds and truncated priors
The bounds are not repeated in the prior. They already live on the parameter, so a symmetric distribution on a parameter bounded below by zero is a half distribution:
m <- lotri({
propSd <- c(0, 0.1)
prior(propSd) ~ dcauchy(0, 5)
})
as.data.frame(m)[, c("name", "lower", "est", "upper", "prior")]
#> name lower est upper prior
#> 1 propSd 0 0.1 Inf dcauchy(0, 5)Here propSd has lower = 0 and a Cauchy
prior, so it is a half-Cauchy. A package generating ‘Stan’ code has
everything it needs to emit the T[0, ] truncation.
Because the bounds are known, lotri can also check that
the prior does not contradict them. A distribution with positive support
on a parameter that allows negative values is an error:
What else is checked
Unknown distributions are rejected, with a suggestion when there is an obvious near match:
lotri({ a <- 1; prior(a) ~ dnorml(0, 1) })
#> lotri syntax error:
#> =================================================================================
#> :001: a <- 1
#> lotri error:
#> bad matrix expression: 'a ~ dnorml(0, 1)'
#> matrix expression should be 'name ~ c(lower-tri)'
#> :002: prior(a) ~ dnorml(0, 1)
#> =================================================================================
#> Error:
#> ! lotri syntax errors aboveSo are the wrong number of arguments, and argument names that do not belong to the distribution:
lotri({ a <- 1; prior(a) ~ dnorm(0) })
#> lotri syntax error:
#> =================================================================================
#> :001: a <- 1
#> lotri error:
#> 'dnorm' is missing argument(s): sd
#> :002: prior(a) ~ dnorm(0)
#> =================================================================================
#> Error:
#> ! lotri syntax errors above
lotri({ a <- 1; prior(a) ~ dnorm(mu=0, sd=1) })
#> lotri syntax error:
#> =================================================================================
#> :001: a <- 1
#> lotri error:
#> 'dnorm' has no argument 'mu'; valid argument(s): mean, sd
#> :002: prior(a) ~ dnorm(mu = 0, sd = 1)
#> =================================================================================
#> Error:
#> ! lotri syntax errors aboveA matrix-valued distribution on a single parameter, a univariate one on a block, a prior on a parameter that does not exist, and two priors on the same parameter are all errors as well.
Getting the priors back out
Priors round trip. They appear in the prior column of
the estimate data frame and of as.data.frame(), and they
are printed back as valid lotri code:
m <- lotri({
tka <- 0.45
label("Ka")
tcl <- c(0, 1, 10)
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
eta.ka ~ 0.3
prior(tka) ~ dnorm(0, 10)
prior(tcl) ~ dlnorm(1, 0.5)
prior(eta.ka) ~ dgamma(2, 1)
prior(eta.cl, eta.v) ~ lkjCorr(2)
})
as.data.frame(m)[, c("name", "est", "condition", "prior")]
#> name est condition prior
#> 1 tka 0.45 <NA> dnorm(0, 10)
#> 2 tcl 1.00 <NA> dlnorm(1, 0.5)
#> 3 eta.cl 0.10 id lkjCorr(2)
#> 4 (eta.cl,eta.v) 0.01 id <NA>
#> 5 eta.v 0.20 id <NA>
#> 6 eta.ka 0.30 id dgamma(2, 1)
as.expression(m)
#> lotri({
#> tka <- 0.45
#> label("Ka")
#> tcl <- c(0, 1, 10)
#> eta.cl ~ 0.1
#> eta.v ~ c(0.01, 0.2)
#> eta.ka ~ 0.3
#> prior(tka) ~ dnorm(0, 10)
#> prior(tcl) ~ dlnorm(1, 0.5)
#> prior(eta.cl, eta.v) ~ lkjCorr(2)
#> prior(eta.ka) ~ dgamma(2, 1)
#> })Note that the block prior comes back as
prior(eta.cl, eta.v), with the whole block recovered, not
just the element it was stored on. Since the deparsed form is valid
input, it can be fed straight back in:
identical(as.data.frame(eval(as.expression(m))), as.data.frame(m))
#> [1] TRUEPriors are matched to parameters by name, never by
position, so they are unaffected by the matrix re-ordering that
lotri does when rcm=TRUE (which is what
nlmixr2 uses):
m <- lotri({
a ~ 1
b ~ c(0, 1)
c ~ c(0.5, 0, 1)
prior(a) ~ dgamma(1, 1)
}, rcm=TRUE)
dimnames(m)[[1]]
#> [1] "c" "a" "b"
attr(m, "lotriPriors")
#> [1] NA "dgamma(1, 1)" NAThe matrix has been re-ordered, and the prior has followed
a rather than staying on the first row.
Using this downstream
A package that generates ‘Stan’ code needs two things from
lotri, and both are available:
- the prior itself, from the
priorcolumn ofas.data.frame()(which is the$iniDfof anlmixr2model) together with that row’slowerandupperfor any truncation, and - the ‘Stan’ spelling of the distribution, from
lotriPriorDists().
m <- lotri({
tka <- 0.45
prior(tka) ~ dnorm(0, 10)
})
.df <- as.data.frame(m)
.p <- .df$prior[!is.na(.df$prior)]
.p
#> [1] "dnorm(0, 10)"
## map the canonical name to the Stan one
.fn <- as.character(str2lang(.p)[[1]])
lotriPriorDists()$stanName[lotriPriorDists()$name == .fn]
#> [1] "normal"which is enough to write
target += normal_lpdf(tka | 0, 10);.
A prior must never be silently ignored
The one thing a consumer of this column must not do is skip it. If a method cannot use a prior, quietly ignoring it means the fit does something other than what the model says, with nothing to tell the user. A prior that is present should either be used or be an error.
rxode2 provides the assertions for this, so an
estimation method can declare what it supports in one line:
-
assertRxUiNoPriors()for a method that cannot use priors at all -
assertRxUiNormalPriors()for a method that supports priors, but only normal ones –dnorm(),stdNormal()and themultiNormal()family, which is what the shorthand above produces when the parameters are correlated. A covariance matrix prior such aslkjCorr()is rejected.
Note that these live in rxode2 rather than here, because
what counts as supported is a property of the estimation method, not of
the specification.