---
title: "Prior distributions in lotri"
author: "Matthew Fidler"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Prior distributions in lotri}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(lotri)
```

# 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(...)`:

```{r}
m <- lotri({
  tka <- 0.45
  tcl <- c(0, 1, 10)

  prior(tka) ~ dnorm(0, 10)
  prior(tcl) ~ dlnorm(1, 0.5)
})

lotriEst(m)
```

Because the statement names the parameter it applies to, prior lines are
**order independent**.  These two blocks are the same:

```{r}
m1 <- lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(0, 10)
})

m2 <- lotri({
  prior(tka) ~ dnorm(0, 10)
  tka <- 0.45
})

identical(m1, m2)
```

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](#the-normal-prior-shorthand).

# What a prior can be attached to

## Population parameters

```{r}
lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(0, 10)
}))
```

## A single between subject variability term

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

attr(m, "lotriPriors")
```

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

```{r}
m <- lotri({
  eta.cl + eta.v ~ c(0.1,
                     0.01, 0.2)
  prior(eta.cl, eta.v) ~ lkjCorr(2)
})

attr(m, "lotriPriors")
```

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:

```{r}
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")]
```

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:

```{r}
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")]
```

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:

```{r error=TRUE}
lotri({
  eta.cl + eta.v ~ c(0.3,
                     0.01, 0.1)
  ~invWishart(4)
  prior(eta.cl, eta.v) ~ invWishart(5)
})
```

Give the scale matrix explicitly only when it differs from the estimates:

```{r}
lotriEst(lotri({
  e1 + e2 ~ c(1,
              0.1, 1)
  prior(e1, e2) ~ invWishart(4, lotri(e1 + e2 ~ c(2,
                                                  0.5, 2)))
}))
attr(lotri({
  e1 + e2 ~ c(1,
              0.1, 1)
  prior(e1, e2) ~ invWishart(4, lotri(e1 + e2 ~ c(2,
                                                  0.5, 2)))
}), "lotriPriors")
```

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

```{r error=TRUE}
lotri({
  e1 + e2 ~ c(1,
              0.1, 1)
  prior(e1, e2) ~ invWishart(1)
})
```

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

```{r}
m <- lotri({
  eta.cl ~ 0.3
  eta.v ~ 0.1
  om.eta.cl ~ 0.01
  om.eta.v ~ 0.04
})

attr(m, "lotriPriors")
```

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

```{r}
diag(m)
```

Correlated omega priors work the same way, including the per-row line form:

```{r}
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]
```

An `om.` name has to match a real between subject variability; it never
quietly creates one:

```{r error=TRUE}
lotri({
  eta.ka ~ 0.3
  om.eta.nope ~ 0.1
})
```

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:

```{r error=TRUE}
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
})
```

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

```{r}
m <- lotri({
  tcl <- 1
  eta.cl ~ 0.3
  tcl + om.eta.cl ~ c(0.01,
                      0.002, 0.005)
})

lotriEst(m)$prior
```

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:

```{r}
attr(m, "lotriPriors")
diag(m)
```

Naming the block with `prior()` means the same thing:

```{r}
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
```

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:

```{r}
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
```

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

```{r error=TRUE}
lotri({
  eta.a ~ 1
  eta.b ~ 1
  prior(eta.a, eta.b) ~ lkjCorr(2)
})
```

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

```{r}
lotriEst(lotri({
  tka <- 1
  tka ~ 4
}))
```

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:

```{r}
m <- lotri({
  tka <- 1
  tcl <- 3
  tv  <- 4
  tcl + tv ~ c(1,
               0.01, 1)
})

lotriEst(m)$prior
```

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:

```{r}
.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)))
```

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

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

identical(m, m2)
```

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

```{r}
lotriEst(lotri({
  tcl <- 3
  tv  <- 4
  tcl + tv ~ c(1,
               0, 1)
}))$prior
```

## Transformations

`sd()`, `var()`, `cor()`, `cov()` and `chol()` all work here too, so the
prior can be written whichever way is most natural:

```{r}
lotriEst(lotri({
  tcl <- 3
  tv  <- 4
  tcl + tv ~ sd(2,
                0.5, 3)
}))$prior
```

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

```{r error=TRUE}
lotri({
  tka <- 1
  tka ~ 0
})
```

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

```{r}
a <- lotri({ tka <- 0.45; prior(tka) ~ dnorm(0, 10) })
b <- lotri({ tka <- 0.45; prior(tka) ~ normal(0, 10) })

identical(a, b)
lotriEst(a)$prior
```

```{r}
.camel <- lotri({ e1 ~ 1; prior(e1) ~ invWishart(2) })
.stan  <- lotri({ e1 ~ 1; prior(e1) ~ inv_wishart(2) })

identical(.camel, .stan)
attr(.camel, "lotriPriors")
```

Arguments may be positional or named, in any order:

```{r}
lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(sd=10, mean=0)
}))$prior
```

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

```{r error=TRUE}
lotri({
  tka <- 0.45
  prior(tka) ~ dt(3)
})
```

Use `studentT()` with its own parameterization instead:

```{r}
lotriEst(lotri({
  tka <- 0.45
  prior(tka) ~ studentT(3, 0, 10)
}))$prior
```

## The supported distributions

`lotriPriorDists()` returns the whole table, including the 'Stan' name for
each distribution, which is what a package generating 'Stan' code needs:

```{r}
d <- lotriPriorDists()
nrow(d)
head(d, 10)
```

The `kind` column says what a distribution may be attached to --
`univariate` for a single parameter, `matrix` and `multivariate` for a
covariance block:

```{r}
subset(d, kind == "matrix", select=c(name, stanName, parNames))
```

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

```{r}
m <- lotri({
  propSd <- c(0, 0.1)
  prior(propSd) ~ dcauchy(0, 5)
})

as.data.frame(m)[, c("name", "lower", "est", "upper", "prior")]
```

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:

```{r error=TRUE}
lotri({
  a <- c(-10, 1, 10)
  prior(a) ~ dlnorm(0, 1)
})
```

# What else is checked

Unknown distributions are rejected, with a suggestion when there is an
obvious near match:

```{r error=TRUE}
lotri({ a <- 1; prior(a) ~ dnorml(0, 1) })
```

So are the wrong number of arguments, and argument names that do not belong
to the distribution:

```{r error=TRUE}
lotri({ a <- 1; prior(a) ~ dnorm(0) })
```

```{r error=TRUE}
lotri({ a <- 1; prior(a) ~ dnorm(mu=0, sd=1) })
```

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:

```{r}
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")]
```

```{r}
as.expression(m)
```

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:

```{r}
identical(as.data.frame(eval(as.expression(m))), as.data.frame(m))
```

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

```{r}
m <- lotri({
  a ~ 1
  b ~ c(0, 1)
  c ~ c(0.5, 0, 1)
  prior(a) ~ dgamma(1, 1)
}, rcm=TRUE)

dimnames(m)[[1]]
attr(m, "lotriPriors")
```

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()`.

```{r}
m <- lotri({
  tka <- 0.45
  prior(tka) ~ dnorm(0, 10)
})

.df <- as.data.frame(m)
.p <- .df$prior[!is.na(.df$prior)]
.p

## map the canonical name to the Stan one
.fn <- as.character(str2lang(.p)[[1]])
lotriPriorDists()$stanName[lotriPriorDists()$name == .fn]
```

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.