Skip to contents

Why 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] TRUE

That 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

Population parameters

lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(0, 10)
}))
#>   name lower  est upper   fix label backTransform        prior
#> 1  tka  -Inf 0.45   Inf FALSE  <NA>          <NA> dnorm(0, 10)

A single between subject variability term

m <- lotri({
  eta.ka ~ 0.3
  prior(eta.ka) ~ dgamma(2, 1)
})

attr(m, "lotriPriors")
#> [1] "dgamma(2, 1)"

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)" NA

The 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] NA

Because both the degrees of freedom and the size of the block are known, an improper prior is caught:

lotri({
  e1 + e2 ~ c(1,
              0.1, 1)
  prior(e1, e2) ~ invWishart(1)
})
#> Error:
#> ! prior 'invWishart' on a 2x2 block needs degrees of freedom greater than 1, but 1 was given

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.1

Correlated 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:

lotri({
  eta.cl + eta.v ~ c(0.3,
                     0.01, 0.1)
  eta.ka ~ 0.5
  prior(eta.cl, eta.v) ~ invWishart(4)
  om.eta.ka ~ 0.01
})
#> Error:
#> ! a model cannot have both degrees of freedom (ie 'invWishart()') and a normal prior (ie 'om.eta ~ 0.1') on its omegas; these are alternatives, not additions

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:

attr(m, "lotriPriors")
#> NULL
diag(m)
#> eta.cl 
#>    0.3

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:

lotriEst(lotri({
  tcl <- 1
  eta.cl ~ 0.3
  eta.v ~ 0.1
  tcl + om.eta.cl + om.eta.v ~ c(0.01,
                                 0.002, 0.005,
                                 0.001, 0.0005, 0.004)
}))$prior
#> [1] "multiNormal(c(1, 0.3, 0.1), lotri(tcl + om.eta.cl + om.eta.v ~ c(0.01, 0.002, 0.005, 0.001, 5e-04, 0.004)))"

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 prior

The 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] TRUE

The per-row line form builds up the block, exactly as it does for etas:

m2 <- lotri({
  tka <- 1
  tcl <- 3
  tv  <- 4
  tcl ~ 1
  tv ~ c(0.01, 1)
})

identical(m, m2)
#> [1] TRUE

When the parameters are uncorrelated the result is simply independent normal priors, since that is what an MVN with a diagonal covariance is:

lotriEst(lotri({
  tcl <- 3
  tv  <- 4
  tcl + tv ~ c(1,
               0, 1)
}))$prior
#> [1] "dnorm(3, 1)" "dnorm(4, 1)"

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:

  1. the R name, when R parameterizes the distribution the same way ‘Stan’ does – dnorm(), dlnorm(), dgamma(), dbeta()
  2. the camelCase name, which is the ‘Stan’ name written the way the rest of this package is written – invWishart(), lkjCorr(), studentT()
  3. 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:

lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(sd=10, mean=0)
}))$prior
#> [1] "dnorm(0, 10)"

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 above

Use studentT() with its own parameterization instead:

lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ studentT(3, 0, 10)
}))$prior
#> [1] "studentT(3, 0, 10)"

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 univariate

The 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_S

Bounds 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:

lotri({
  a <- c(-10, 1, 10)
  prior(a) ~ dlnorm(0, 1)
})
#> Error:
#> ! prior 'dlnorm' has positive support but 'a' has a lower bound of -10

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 above

So 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 above

A 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:

Priors 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)" NA

The 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:

  1. the prior itself, from the prior column of as.data.frame() (which is the $iniDf of a nlmixr2 model) together with that row’s lower and upper for any truncation, and
  2. 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 the multiNormal() family, which is what the shorthand above produces when the parameters are correlated. A covariance matrix prior such as lkjCorr() 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.