same(), which is NONMEM's
$OMEGA BLOCK(n) SAME:lotri({
iov.cl1 + iov.v1 ~ c(0.1,
0.01, 0.2)
iov.cl2 + iov.v2 ~ same()
iov.cl3 + iov.v3 ~ same()
})
This is one estimated 2x2 covariance shared by three blocks, which is how inter-occasion variability is parameterized when every occasion draws its own random effects from one shared covariance -- and it is what lets those random effects be correlated.
same() repeats the immediately preceding block under new names and
takes no arguments. A further same() repeats that same original
block rather than the copy, the way NONMEM chains SAME. It works
with a condition (iov.cl2 + iov.v2 ~ same() | occ), it inherits the
fixed flags of the block it repeats, and it composes with the
cnd(same = n) nesting property that lotriSep() uses.
The condition column can now carry a :same: suffix. Rather
than adding a column to the data frame from as.data.frame(), a
repeated block records the element it mirrors in the existing
condition column, naming it:
<baseCondition>:same:<masterEta> # diagonal row
<baseCondition>:same:<masterEta1>:<masterEta2> # covariance row
The master is named rather than indexed because neta1/neta2 are
renumbered whenever parameters are added, dropped or reordered. A
copy keeps its master's est and fix, so a consumer that does not
understand the suffix still sees a numerically correct matrix -- just
with more free parameters than the model really has.
Four helpers are exported for reading that column, and code that
compares condition directly (condition == "id") should move to
them, since such a test misclassifies a repeated block:
lotriBaseCondition() strips the suffix,lotriIsSame() says whether a row mirrors another,lotriSameMap() maps each eta to the eta it mirrors,lotriSameBreak() drops the linkage for a block that has been
structurally changed.same() round trips through as.data.frame(), as.lotri() and
as.expression(), and the linkage is carried across lotriMat().
A repetition that no longer describes a real mirror -- because the
blocks it was built from were reordered or dropped, or because it
carries something same() cannot express -- is written out with its
explicit values instead. That loses the annotation, never the
numbers.
A prior cannot be put on a repeated block -- it is not a parameter of its own, so the prior belongs on the block it mirrors. Without this a model could carry two different priors on what is one estimated parameter.
same() cannot be combined with a cov function, nor with
rcm=TRUE on a matrix the permutation would actually reorder; both
would move a repeated block away from the block it repeats, and are
refused rather than silently producing a matrix whose repetition is
no longer true. A same() matrix that is already block diagonal is
left alone, since there rcm is a no-op, which keeps rxode2's
ini({}) -- which always passes rcm=TRUE -- working.
Prior distributions can now be specified in a lotri({}) (and
therefore ini({})) block with prior(name) ~ dist(...), ie:
lotri({
tka <- 0.45
tcl <- c(0, 1, 10)
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
prior(tka) ~ dnorm(0, 10)
prior(tcl) ~ dlnorm(1, 0.5)
prior(eta.cl, eta.v) ~ lkjCorr(2)
})
Because the statement names its target, prior lines may be given
anywhere in the block. Priors may be put on population estimates,
on individual etas, and on whole covariance blocks (with the
matrix-valued distributions like lkjCorr() and invWishart()).
Every distribution has three accepted spellings: the R name where R
parameterizes it the same way 'Stan' does (dnorm(), dlnorm(),
dgamma()), the camelCase name (invWishart(), lkjCorr(),
studentT()), and the 'Stan' name itself (inv_wishart(),
lkj_corr(), student_t()). The canonical spelling, which is what
is stored and printed back, is the R one where there is a faithful
one and the camelCase one otherwise. Both positional and named
arguments work. lotri validates the distribution name, its arity,
and its support against the parameter's bounds; it does not generate
any 'Stan' code.
$OMEGAP/$OMEGAPD pair of a NONMEM
NWPRI model:lotri({
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
prior(eta.cl, eta.v) ~ invWishart(4)
})
It works on a 1x1 block too (an inverse Wishart of dimension one is
an inverse gamma), and an improper nu <= p - 1 is an error.
TNPRI model needs. om. prepended to a between
subject variability names its omega element, so the normal prior
shorthand can be used on it:lotri({
eta.cl ~ 0.3
eta.v ~ 0.1
om.eta.cl ~ 0.01
om.eta.v ~ 0.04
})
The omega itself is untouched; only the prior is added. Correlated
omega priors work the same way, including the per row line form. An
om. name has to match a real between subject variability and never
creates one. prior(om.eta.cl) and prior(eta.cl) mean the same
thing; the om. spelling exists so the shorthand has a name to put
on the left of a ~, since eta.cl ~ ... already means the omega
value.
Degrees of freedom and a normal prior are alternative ways of putting a prior on an omega, so a model that gives both is an error.
~ gives it a normal prior with a zero mean and the given
variance.lotri({
tka <- 1
tcl <- 3
tv <- 4
tka ~ 4 # tka ~ N(0, sd=2)
tcl + tv ~ c(1, # (tcl, tv) ~ MVN(0, Sigma)
0.01, 1)
})
All of the matrix spellings work here, including the per row line form
(tcl ~ 1; tv ~ c(0.01, 1)) and the sd(), var(), cor(), cov()
and chol() transformations. The <- value remains the initial
estimate; it is not the prior mean. An uncorrelated block simply
becomes independent normal priors. A zero variance is an error.
Note this changes behavior: lotri({b <- 3; b ~ 0.4}) used to be a
"duplicated parameter" error and is now a normal prior on b. A name
that is not an estimate still specifies an eta as before.
Added lotriPriorDists() which returns the table of supported
distributions (including the 'Stan' name for each), so that
downstream packages can generate the corresponding 'Stan' code.
The lotriEst data frame and as.data.frame() output gained a
prior column.
prior() also takes the normal prior shorthand, so
prior(tka) ~ 0.1 means what tka ~ 0.1 means:
lotri({
tka <- 0.45
prior(tka) ~ 0.1 # dnorm(0.45, sqrt(0.1))
prior(tcl, tv) ~ c(1, # multivariate, centered on the estimates
0.01, 1)
})
Every matrix spelling works here as it does bare, including the per row line form, so a covariance can be built up a line at a time:
lotri({
tcl <- 3
tv <- 4
prior(tcl) ~ 1
prior(tv) ~ c(0.001, 1) # same 2x2 as `tcl ~ 1; tv ~ c(0.001, 1)`
})
An uncorrelated group becomes independent normal priors, and the mean
is what the model already says. Note the line form is the one place a
prior line is not order independent, since a row leans on the line
before it. The point is that a bare ~ cannot be used
everywhere -- piping onto a model reads tka ~ 0.1 as changing the
estimate -- so the prior() flag gives the shorthand a spelling that
works when piping too.
A namespaced distribution such as stats::dnorm(0, 1) is now an
error rather than being evaluated as a variance.
dnorm(), dcauchy(), ...) can now be put on a
single omega COVARIANCE (off-diagonal) element, not just a diagonal one --
a marginal, independent prior on that one cell, distinct from a whole-block
distribution like invWishart()/multiNormal():lotri({
eta.cl + eta.v ~ c(0.1,
0.01, 0.2)
prior(eta.cl, eta.v) ~ dnorm(0, 0.1)
})
The two names must covary with each other but need not be the model's
entire connected block -- prior(eta.cl, eta.v) works even when a third,
correlated eta is also in the block, unlike a whole-block prior which
still requires naming every member. prior(om.eta.cl, om.eta.v) ~ ...
is accepted identically. A block cannot carry both a whole-block
invWishart()/multiNormal() prior and a marginal prior on one of its
own cells -- these remain alternatives, not additions.
A condition written on a continuation no longer relocates unrelated parameters. The rows of one block share a level of variability because they covary, so a condition on a later row carries that block over -- but it used to carry the whole default level with it, so parameters declared long before, and unrelated to the block, silently changed level.
A line that cannot be parsed at the default level is no longer folded into the most recently seen condition. It is folded only when the block immediately before it is at that level, which is what the fallback was for; otherwise a parameter with no condition could land at a level of variability it was never given.
Labels are no longer lost, or attached to the wrong parameter, when a
block has more than one level of variability. Which level a trailing
label() belonged to was decided by whether the DEFAULT level had any
labels yet, so once it did, a conditioned line's label landed on the
default level and overwrote the one already there; and
as.data.frame() only ever read the labels of a single matrix, so a
conditioned object lost them entirely.
as.lotri() on a data frame whose only level of variability is not
id now keeps that level's name. It previously returned a bare
matrix, so an occasion-only model came back looking like an id level
one.
lotriEst(x, drop=TRUE) kept the lotriFix class only when there
were fixed elements, so a matrix carrying labels, priors or a
same() repetition came back unclassed with the attribute orphaned
and every consumer dispatching to the default method.
A plus-form block, or a block whose right hand side is arithmetic
(d ~ 0.1*2), written after a line-form block overwrote the rows
the line form had just written. A line-form block leaves the row
counter pointing at its FIRST row with the rest held separately, and
neither branch settled that, so
lotri({a ~ 1; b ~ c(0.1, 2); c1 + d1 ~ c(3, 0.2, 4)}) failed with
"length of 'dimnames' [1] not equal to array extent". The same
applied to a 1x1 block under a condition.
A conditioned block written in the line form lost every row after the
first. a ~ 1 | occ did not set the lastN counter the line form
accumulates through, so a following b ~ c(0.1, 2) | occ could not
see the row before it and the block was silently truncated to 1x1.
This showed up as eval(as.expression(x)) returning a smaller matrix
than x for any conditioned block of more than one parameter.
lotri() read its own condition-property list with a partially-matching
attr(). Because "lotri" is a prefix of "lotriLabels",
"lotriFix" and friends, combining a conditioned matrix with a labelled
one (as in lotri(lotri(b ~ 2) | occ, lotri({c ~ 1; label("z")})))
picked up the neighbouring attribute, attached it as the property list,
and gave the result a spurious lotri class and a bogus Properties:
line. The internal reads are now exact.
lotri now requires 'armadillo4r' 15.4.2 or newer. The 'Armadillo'
headers shipped with older 'armadillo4r' releases discard the return
value of a std::uniform_int_distribution draw, which the 'libc++'
shipped with 'clang' 23 now marks [[nodiscard]]. That made R CMD check report a significant compilation warning -- and so a check
WARNING -- on CRAN's clang-trunk flavour, even though none of
lotri's own sources were involved (#57).
src/nearPD.cpp and src/rcm.cpp now include <algorithm>,
<cmath> and <cstddef> for the std::fill/std::fill_n,
std::max, std::abs and std::size_t they use directly, instead
of relying on the 'Armadillo' headers to drag them in. 'libc++' has
been steadily dropping transitive includes, and clang 23's is where
that started breaking packages (#57).
The trace=TRUE iteration message in lotriNearPD() printed an
arma::uword with %lld; it now uses %llu with an explicit cast,
so the format matches the argument on every 'Armadillo' word-size
configuration (#57).
The test suite no longer calls structure() with the deprecated
special names .Dim, .Dimnames and .Names; it uses dim,
dimnames and names instead. This clears the "Found calls to
structure() using deprecated special names" NOTE on the r-devel
check flavours (#57).
A prior() given for a fixed (fix()ed) parameter is now an error
at resolution time, naming the parameter, instead of parsing and
building cleanly and only surfacing as a problem in whatever
consumes the priors. A fixed parameter is a constant, so it cannot
carry a prior; this applies to a univariate prior row, a member of a
multivariate prior, and an omega-block member on a fixed omega
element alike, and also to a block prior (lkjCorr(),
invWishart()) whose block has a fixed covariance even if every
variance in it is free. The implicit ~invWishart(4) shorthand,
which applies to every omega block in the model, quietly skips a
block that is entirely fixed instead (#52).
lotri() can now re-parse the named c() form that
lotriAsExpression()/as.expression() writes for a matrix above
nameEst's size threshold when a later row's names reach back into an
earlier, already-closed block -- as happens for a combined theta+omega
covariance matrix, where every omega row names all the preceding theta
rows too. Previously an om.-prefixed row in that position was always
misread as the om. normal prior shorthand instead of the next row of
the matrix itself (#53).
Labels now follow the matrix when rcm=TRUE re-orders it.
Previously the labels stayed in the order they were parsed in while
the matrix was permuted, so they were applied to the wrong
parameters.
lotriLabels are no longer dropped when matrices are combined (ie
lotri(mat1, mat2) or lotriMat()); they are now concatenated in C
along with the matrix itself.
as.expression() now works on a lotri object that has only
population estimates and no matrix; it used to fail with "second
argument must be a list".
lotriNearPD(x, only.values=TRUE) no longer fails with "unknown
c++ error" when x is already positive definite. The eigenvalues
were only returned from inside the branch that clamps a negative
eigenvalue, so an input that needed no clamping fell through and
tried to assign the full matrix into a length-n result.
Fixed rchk issues and small bugs found while linting
Fix unsigned integer underflow in rcm.cpp: when called with a
0-row matrix, n-1 would wrap to UWORD_MAX causing an infinite
loop; added an early-return guard for n == 0.
Fix type narrowing and underflow in nearPD.cpp: changed unsigned int n to arma::uword to prevent silent 32-bit truncation; added
n == 0 early-return guard to prevent out-of-bounds access on the
eigenvalue vector.
Fix variable shadowing in lotriProp.c: local double val renamed
to inVal to eliminate -Wshadow warning against the function
parameter.
Fix hard-coded buffer size in matlist.h: snprintf buffer
increased from 100 to 256 bytes and size passed via sizeof(out).
Use R_xlen_t instead of int for variables storing Rf_length()
and Rf_xlength() returns across matlist.h, lotriProp.c,
lotriBounds.c, and lotriLstToMat.c.
cpp4r and armadillo4r (issue #41)Before you could specify matrices as:
m <- lotri({
a + b ~ c(1,
0.5, 1)
})
Now you can specify per row as:
m <- lotri({
a ~ 1
b ~ c(0.5, 1)
})
This form is now the default when converting from a matrix to a lotri expression. In addition if the matrix is large enough (by default a 5x5 matrix), these would be named when changing them to an expression:
m <- lotri({
a ~ c(a=1)
b ~ c(a=0.5, b=1)
c ~ c(a=0.5, b=0.5, c=1)
d ~ c(a=0.5, b=0.5, c=0.5, d=1)
e ~ c(a=0.5, b=0.5, c=0.5, d=1,
e=1)
})
This way changing to an R parsed expression will be rendered in a more human readable format.
You can change the deparsing options that are used by default with
lotri with options(lotri.plusNames=TRUE) which prefers the a+b+c
syntax when deparsing. Otherwise, the line format is used by default.
The dimension number before naming the values in the line-format can
be controlled with options(lotri.nameEst=2) or some other dimension.
New option of cov added which check for matrix suitability for
covariance matrix. When cov=TRUE, off-diagonal elements in
covariance matrices may no longer be nonzero if the diagonal value
is zero (rxode2#481). This will also check tht the matrix is
non-positive definite on the non-diagonal terms. cov can also be a
function to allow correction of the matrix to a positive definite
matrix automatically.
New option of rcm; When enabled, and lotriIsBlockMat() is not
true, lotri() will permute the matrix to try to get a banded
matrix using the Reverse Cuthill McKee algorithm.
Change internals for lotri so that new rxode2 is no longer
required to be binary linked to lotri.
Add new function rcm() which permutes the matrix to get a band
matrix (if possible). This uses the Reverse Reverse Cuthill McKee
(RCM) algorithm.
Moved nmNearPD() to this package and renamed to lotriNearPD().
In addition to moving, this function will now retain the dimension
names.
New exported function lotriAsExpression() which has more fine
control than as.expression() and will work without converting the
matrix to a lotri form. This by default uses the new line form, but
can be changed back to the option(lotri.plusNames=TRUE). Also be
default it will name each element in a matrix when the dimension is
above 5x5. You can change that number by
option(lotri.nameEst=10) to increase it to be named above
10x10. If you do not like the naming you can also disable it with
option(lotri.nameEst=TRUE), or if you always want it on you can use `option(lotri.nameEst=FALSE)
lotri(n1+n2~omega)Can convert lotri objects to data-frames similar to the internal
data frame used in nlmixr(); These can then be converted back with
as.lotri()
Added the ability to add estimates to a lotri object. These
estimates are an attached data.frame to the original lotri matrix.
You can extract them or drop them with the function lotriEst()
Allow specifying fixed components in lotri() matrices.
Add cov, cor, sd, var, and chol options for matrix
specification. The final matrix will always be the covariance matrix
Add function lotriMatInv() which takes a symmetric block matrix
and converts it into a list of matrices. An sort of inverse
operation of lotriMat()
Add error for lotri(~c(40))
Added ability to flag fixed and unfixed components in a
matrix. Currently lotri only supports one type.
For the lotriFix objects, ie those created with population types
of estimates and covariance estimates, allow them to be converted to
an equivalent expression with as.expression() for the lotri object
and lotriDataFrameToLotriExpression() for the data.frame
Added lotriIsBlockMat() to check to see if a matrix is in a block
diagonal matrix form.
call.=FALSE or equivalent.rchklotriMat to mix named and unnamed matrices; When mixed, an
unnamed matrix will be returned.$lower and $upper gives default values even if it
wasn't specified.Matrix import and added lotriMat to create banded
matrices (faster than Matrix for now included repeated matrices
with list(matrix, rep)).upper and lower bounds
are specified|; Returns a list of matrices
where each condition is the name of the matrix returned.| id(lower=c(eta1=3))NEWS.md file to track changes to the package.