# Prior Elicitation for Generalised Linear Models and Extensions

Geoffrey R. Hosack1

(1CSIRO–Data61
Hobart, Tasmania 7001, Australia
E-mail<EMAIL_ADDRESS>)

###### Abstract

A statistical method for the elicitation of priors in Bayesian generalised
linear models (GLMs) and extensions is proposed. Probabilistic predictions are
elicited from the expert to parametrise a multivariate t prior distribution
for the unknown linear coefficients of the GLM and an inverse gamma prior for
the dispersion parameter, if unknown. The elicited predictions condition on
defined elicitation scenarios. Dependencies among scenarios are then elicited
from the expert by additionally conditioning on hypothetical experiments.
Elicited conditional medians efficiently parametrise a canonical vine copula
model of dependence that may be truncated for efficiency. The statistical
elicitation method permits prior parametrisation of GLMs with alternative
choices of design matrices or observation models from the same elicitation
session. Extensions of the method apply to multivariate data, data with
bounded support, semi-continuous data with point mass at zero, and count data
with overdispersion or zero-inflation. A case study elicits a prior for an
extended GLM embedded in a statistical model of overdispersed counts described
by a binomial-simplex mixture distribution. The elicited canonical vine model
of dependence is found to incorporate substantial information into the prior.
The procedures of the statistical elicitation method are implemented in the R
package eglm.

##  1 Introduction

Generalised linear models (GLMs; Nelder and Wedderburn, 1972; McCullagh and
Nelder, 1989) provide foundational support for applied statistical inference
in science, engineering and medicine. To begin the process of Bayesian
inference for the unknown parameters of a GLM, a prior must first be
specified. This first inferential step can entail a statistical and scientific
investigation similar to an empirically focused study (Garthwaite et al.,
2005; Kuhnert, 2011). A general and structured elicitation procedure that
supports Bayesian GLMs and related models is therefore needed. The statistical
method for elicitation presented here is applicable to GLMs that allow for a
range of response data with either known or unknown dispersion. Moreover, the
elicitation method accommodates extensions to statistical models for
multivariate, bounded support, semi-continuous, overdispersed or zero-inflated
data.

GLMs combine random and systematic components to describe data (Nelder and
Wedderburn, 1972, Section 2.1). Elicitation methods for GLMs are often based
on techniques originally developed for the normal linear model (Section 2.2).
General approaches to elicit unknown dispersion parameters, however, are
unavailable for non-normal GLMs (Section 2.3). Comparatively more effort has
been placed on elicitation of the unknown linear coefficients of GLMs. The
present approach for elicitation generalises the concept of an independent
conditional means prior for the linear coefficients (ICMP; Bedrick et al.,
1996). An ICMP uses an assumption of conditional independence given the
covariates to ease the complexity of the elicitation process (Section 2.4).
Statistical methods of elicitation for GLMs have been proposed that relax the
assumption of conditional independence among elicited responses (Section 2.5).
These elicitation methods each emphasise most or all of the elicitation
principles summarised by Kadane and Wolfson (1998) for predictive elicitation:
1) expert opinion is worthwhile to elicit; 2) only observable quantities
should be assessed via an indirect rather than direct elicitation of unknown
model parameters; 3) experts should be asked for quantiles or probabilities
and never be asked to estimate moments greater than the mean; 4) feedback
should be provided to the expert; and 5) assessments should be provided
conditionally and unconditionally on hypothetical data.

The statistical elicitation method outlined in Section 3 also incorporates
these same principles, and additionally the following. By design the method is
_justifiable_. The expert builds a probability model that re-expresses the
elicited probability assessments without loss of information. All elicited
quantities are univariate because experts typically do not accurately assess
joint probabilities (Garthwaite et al., 2005). Graphical and numerical
feedback mechanisms (e.g., O’Hagan et al., 2006) are provided that support the
expert’s ability to learn, adjust and critique assumptions during the
elicitation session.

The method is also _sequential_. The elicitation first addresses the
dispersion parameter that may be known or unknown (Section 4). Given this
characterisation, marginal elicitations for the mean responses proceed
analogously to an ICMP (Section 5.1). Dependence in the elicited mean
responses is then modelled by a canonical vine copula (Bedford and Cooke,
2001; Joe, 2015, see Section 5.2) that efficiently captures expert opinion
(Section 5.3). This approach avoids asking the expert to re-evaluate central
credible intervals conditionally, which is a more difficult task compared to
the elicitation of marginal credible intervals (Al-Awadhi and Garthwaite,
1998, 2001; Elfadaly and Garthwaite, 2020) and requires additional elicited
quantiles, but rather elicits only conditional medians instead. These
elicitations condition on hypothetical data that are sequentially introduced
(Section 5.4). Once introduced, the hypothetical data are retained to ease the
cognitive burden on the expert (sensu Kadane et al., 1980). Moreover, the
sequential introduction of hypothetical data supports truncated canonical vine
dependence models that further lessen the elicitation burden (Section 5.5).
Truncated vines are helpful in a situation where adding further information
into the conditioning set of hypothetical response data has no meaningful
effect on the beliefs or predictions of the expert, or where time and resource
constraints limit the full exploration of dependencies.

An additional design principle is _generality_. The statistical elicitation
method not only applies to a broad class of models (sensu Kadane and Wolfson,
1998; Mikkola et al., 2021), but also the same elicitation session can be used
to parametrise alternative GLMs. At the time of the elicitation steps
described by Sections 4 and 5, it is sufficient to have specified only the
mean domain, link function and the elicitation scenarios for consideration by
the expert. Alternative model matrices or observation models may be specified
after the elicitation session concludes to induce alternative prior
specifications for GLMs (Section 6). The generality of the statistical method
for elicitation is further illustrated in Section 7 with extensions to
multivariate, semi-continuous, bounded support, zero-inflated and
overdispersed data.

An elicitation case study demonstrates an overdispersed binomial model in
Section 8, where the elicited quantity is modelled with an extended GLM. The
conclusion summarises future directions (Section 9), such as the development
of open and free software that is an important aspect of statistical methods
for elicitation (Kadane et al., 1980; Garthwaite et al., 2005; Mikkola et al.,
2021). The freely available R package eglm (Hosack, 2023) implements the
procedures of the statistical elicitation method.

##  2 Review of GLMs and Elicitation Methods

###  2.1 Generalised Linear Models

The Generalised Linear Model (GLM) extends regression analysis developed for
the normal linear model to the case where observables follow an exponential
dispersion model (Jørgensen, 1997b). An exponential dispersion model for a
random variable $\psi$ has the density function

$f(\psi\mid\theta,\lambda)=a(\psi,\phi)\exp\left\\{-\lambda[\psi\theta-b(\theta)]\right\\}$
(1)

for functions $a(\cdot)$ and $b(\cdot)$. The mean is
$\mathbb{E}\left[\psi\right]=\mu\in\Omega\subseteq\mathbb{R}$, where
$\mu=\tau(\theta)=db(\theta)/d\theta$. The variance is
$\mathbb{V}\left[\psi\right]=\phi v(\mu)$, where $\phi=1/\lambda>0$ is the
dispersion parameter and $v(\mu)=d^{2}b(\tau^{-1}(\mu))/d\mu^{2}$ is the
variance function. An exponential dispersion model with mean $\mu$ and
dispersion parameter $\phi$ is denoted by $\psi\sim ED(\mu,\phi)$.

If the dispersion parameter is known then Eq. (1) corresponds to a natural
exponential family model (Morris, 1982; Jørgensen, 1997b, Ch. 2). If the model
is closed with respect to scale transformation such that
$cED(\mu,\phi)=ED(c\mu,c^{2-p}\phi)$ for $c>0$ given power parameter
$p\in(-\infty,0]\cup[1,\infty)$ with variance function $v(\mu)=\mu^{p}$, then
Eq. (1) corresponds to a Tweedie exponential dispersion model (Jørgensen,
1987). The three standard continuous observation models used in GLMs (e.g.,
Lindsey, 1997) are special cases: normal ($p=0$), gamma ($p=2$) and inverse
Gaussian ($p=3$). The discrete Poisson distribution is also a Tweedie model
($p=1$), as are the semi-continuous compound Poisson distributions with point
mass at zero ($1<p<2$).

A GLM (e.g., McCullagh and Nelder, 1989) is defined by:

  1. 1.

The observations $\psi_{j(i)}\sim ED(\mu_{i},\phi_{i})$ for $i=1,\ldots,n$,
$j(i)=1,\ldots,J(i)$ and $\phi_{i}=\phi/w_{i}$ with known weight $w_{i}$.

  2. 2.

The linear predictor $\eta_{i}=x_{i}^{\top}\beta+o_{i}$, $i=1,\ldots,n\geq p$
with $p$-dimensional vector $\beta$ of unknown parameters. The $p$-dimensional
vector of covariates $x_{i}^{\top}$ forms the $i$th design point, which is the
$i$th row of the $n\times p$ model matrix $X$ of full column rank. The offset
$o_{i}$ for the linear predictor at the $i$th design point is assumed known.

  3. 3.

A continuous link function $g$ invertible over the mean domain $\Omega$ such
that $g(\mu)=\eta$.

The GLM consists of two basic components (Nelder and Wedderburn, 1972): The
“systematic component” defined by $\mu=\mathbb{E}\left[\psi\mid
X,\beta,\phi\right]=g^{-1}(X\beta)$ with $\beta$ unknown, and the “random
component” that describes the uncertainty in $\psi$ conditional on $\mu$ and
the dispersion parameter $\phi$.

###  2.2 Elicitation for the Normal Linear Model

Elicitation procedures developed for the normal linear model have had a major
impact on the development of probabilistic elicitation methods for GLMs
(Garthwaite et al., 2005). An elicitation approach for the special case of a
normal linear regression model with unknown variance was initially proposed by
Kadane et al. (1980). The procedure occurs in two stages. In the first stage,
marginal quantiles are elicited at each design point defined by the covariates
that correspond to a hypothetical observation. The mean is estimated by
ordinary least squares and the degrees of freedom parameter then estimated by
averaging over a function of the elicited quantiles evaluated at each design
point. In the second stage, conditional medians and quartiles are elicited at
each of the design points given hypothetical observations. These results are
then used to obtain an averaged estimate of the scale matrix for $\beta$ (see
Kadane et al., 1980; Kadane and Wolfson, 1998, for details).

Various alternatives to the method of Kadane et al. (1980) have been
suggested. Al-Awadhi and Garthwaite (1998) note that unconditional estimates
of central credible intervals (spreads) should be more accurate than
conditional estimates, and propose using the method of Kadane et al. (1980) to
generate preliminary estimates of correlations that are then rescaled by
marginally estimated spreads. Al-Awadhi and Garthwaite (1998) also investigate
an alternative method that avoids elicitation of conditional spreads
altogether, but which does not ensure that the elicited covariance matrix is
positive definite. Cooney and White (2023) note that the method of Kadane et
al. (1980) requires a large number of elicited quantiles to parametrise a
multivariate normal. Their suggested approach asks fewer questions of an
expert by eliciting concordance probabilities for pairwise jointly distributed
unknowns to parametrise partial correlations. The elicited concordance
probabilities are not guaranteed to conform with a valid correlation matrix in
high dimensions. Garthwaite and Dickey (1988) describe an alternative method
that depends on the choice of design points made by the expert. This approach
limits application to regression problems with independent design points, and
so prohibits for example polynomial regression or models that require factor
coding.

###  2.3 Unknown Dispersion Parameters in GLMs

Probabilistic elicitation of unknown dispersion parameters within GLMs are
uncommon outside of the normal linear model. For a GLM with a gamma
observation model, Elfadaly and Garthwaite (2015) elicit an unknown dispersion
parameter conditional on an assumed mean value for $\mu$. This elicitation
approach relies on a scaling property of the gamma distribution with respect
to its mean, which arises as a special case of the Tweedie class (Section
2.1), such that $cED(\mu,\phi)=ED(c\mu,c^{2-p}\phi)=ED(c\mu,\phi)$ for $p=2$
and $c>0$. A prior for the ratio of a hypothetical observation $y$ relative to
the mean $\mu$ is then parametrised conditional on a given value for $\mu$.
However, this scaling property available for the gamma distribution does not
generalise to other exponential dispersion models commonly used in GLMs.

###  2.4 Independent Conditional Means Priors for GLMs

Bedrick et al. (1996) describe an approach to elicit the coefficients of the
linear predictor for GLMs within the natural exponential family, or
equivalently, where the observation distribution is described by an
exponential dispersion model with known dispersion parameter $\phi$.
Subjective probability distributions are elicited for the mean of the response
$\mu_{i}$ conditional on a set of covariates $x_{i}^{\top}$ for
$i=1,\ldots,n$. In its most basic formulation, the $n\times p$ model matrix
$X$ for the elicitation session is full rank with $n=p$, and chosen so that
the elicitations for $\mu$ may be thought of as conditionally independent
given $X$. Bedrick et al. (1996) demonstrate how a conditional independence
assumption for the elicited distribution of the unknown means can be
reasonable, and also suggest guidelines for the choice of $X$ such that the
independence assumption is credible. Given a link function $g$ with linear
model $g(\mu)=X\beta$, the induced prior for $\beta$ is known as an
independent conditional means prior (ICMP). A normal approximation of the
prior for $\beta$ can be a plausible assumption for ICMPs (Bedrick et al.,
1996). Hosack et al. (2021) apply an ICMP approach to elicit multivariate
normal priors for the unknown coefficients $\beta$ in GLMs. Denham and
Mengersen (2007) describe an approach for logistic regression that is related
to ICMPs, where Monte Carlo approximation is used to derive a prior based on
samples drawn from elicited beta distribution that are assumed independent
conditional on the model matrix $X$.

###  2.5 GLMs with Dependent Conditional Means Priors

The assumption of independence among elicited responses within an ICMP
approach may be overly restrictive for GLM priors. Conditional on $\phi$ and
$X$, Chen and Ibrahim (2003) describe an elicitation approach that elicits a
mean response and an associated degree of belief to parametrise a prior for
$\beta$. Based on asymptotic results for this prior conditional on an
imaginary sample, Bové and Held (2011) derive a class of multivariate normal
$g$-priors and also review other related normal $g$-priors proposed for GLMs.
For both the prior of Chen and Ibrahim (2003) and a normal $g$-prior, the
modelled a priori dependence among design points (conditional means) is
entirely determined by the choice of model matrix. Robert (2007, Ch. 4) notes
that for this reason $g$-priors may usefully substitute for prior information
where such information is sparse. However, the above priors do not elicit
dependent conditional mean priors in GLMs and are not further evaluated here.

Recent approaches have elicited dependent conditional mean responses in GLMs.
An elicitation method for piecewise-linear and non-normal GLMs developed by
Garthwaite et al. (2013) allows for limited dependence among conditional means
given a normal prior specification for $\beta$. The model structure is
restricted to assume no interactions among covariates, and the elicited
covariance matrix of the multivariate normal prior for $\beta$ is assumed
block-diagonal. For multinomial GLMs, Elfadaly and Garthwaite (2020) adapt the
elicitation algorithm of Kadane et al. (1980) in Section 2.2 to parametrise a
normal prior for $\beta$. First marginal medians are elicited for category
odds, then conditional central credible intervals are elicited in a second
stage. Although not working with covariates, Wilson (2018) propose the use of
$D$-vine copulas (Joe, 2015, Ch. 3) to elicit priors for the multinomial
distribution. The $D$-vine copula is related to the canonical vine
representation of dependency that will be introduced in Section 5. Unlike the
$D$-vine copula, however, it will be seen that the canonical vine copula
sequentially adds new conditioning information into an elicitation session in
such a way so that the expert is never asked to forget information once
introduced into the elicitation session.

##  3 Statistical Elicitation Method for GLMs

###  3.1 Definitions for Elicitation

The expert has the domain knowledge elicited by the session facilitator
(Garthwaite et al., 2005) according to a model specified by the statistician.
Each of these roles may be fulfilled by the same individual, or by different
individuals or by some combination thereof. More than one individual may
contribute to each role. Garthwaite et al. (2005) and O’Hagan et al. (2006)
provide guidance on preparation for an elicitation session. In particular, the
facilitator provides the expert with probabilistic training before the
elicitation. This training incorporates well-established statistical concepts
that may not be fully appreciated by all domain experts. The definition of
quantiles (or percentiles or fractiles) and their use in probability
assessments is particularly important. Often medians are used to assess
location while inter-quartile or inter-tertile distances assess spread.
Another important aspect of the statistical training is education on the
sensitivity of statistical estimators to sample size. Without proper training,
experts can either overemphasise (Tversky and Kahneman, 1974) or
underemphasise (Garthwaite et al., 2005) the information contained within a
sample for a summary statistic such as the sample mean. Education for an
expert unfamiliar with the Central Limit Theorem (CLT), for example, may
include simulated examples (e.g., Dinov et al., 2008; Zhang et al., 2022) as
part of the expert training.

In the structured elicitation session, the expert is asked by the facilitator
to envision or predict hypothetical realisations of the elicitation target
$\psi$ that arises from a GLM (Section 2.1). The facilitator ensures that the
elicitation target is relevant and interpretable to the expert. The
elicitation session is supported by elicitation scenarios that are
conceptually defined by a $n\times p^{\prime}$ matrix $U$ for the expert. The
$i$th row of $U$ defines the $i$th scenario, $u_{i}^{\top}$, that is a
$p^{\prime}\times 1$ vector of measurable covariates. Elicitation scenarios
may be described numerically or by visualisations that might, for example,
relate to geographical covariates (e.g., Denham and Mengersen, 2007; Hosack et
al., 2017). The facilitator ensures that the scenarios are relevant and also
numerically or graphically interpretable to the expert.

The statistician describes the elicitation target by a set of exponential
dispersion models $\mathcal{F}$ with shared mean domain $\Omega$ that is
consistent with the choice of link function mapping
$g:\Omega\rightarrow\mathbb{R}$. The set of exponential dispersion models
$\mathcal{F}$ allows for distributions with different variance functions at
the time of elicitation, if desired. The systematic component of the GLM
depends on the model matrix $X$ that is linked to the elicitation scenarios of
$U$. The design point $x_{i}^{\top}$, which is the $i$th row of $X$, is a
function of the covariates that define the $i$th elicitation scenario
$u_{i}^{\top}$ for $i=1:n$. The expert is thus shielded from choices of factor
coding, covariate scaling or other modification to $X$ made by the
statistician that should not affect the elicitation. Many alternative choices
of $X$ with full column rank equal to $p\leq n$ may be consistent with $U$.
The chosen set of these alternative model matrices is denoted by
$\mathcal{X}$.

There is process uncertainty about how the systematic component $\mu$ relates
to the scenarios of $U$ through the covariates $X$ as determined by the
unknown $\beta$. The random component accounts for observation uncertainty,
induced by sampling variability or measurement error or experimental error,
that is unexplained by the covariates. The observation uncertainty is
described by an exponential dispersion model (Section 2.1), where the
dispersion parameter $\phi$ may be known or unknown. The observation
uncertainty of a sample mean $\bar{\psi}=(1/w)\sum_{j=1}^{w}\psi_{j}$, given
$\mu_{i}$ and $\phi$ at scenario $u_{i}^{\top}$, is reduced with an increasing
number of samples $w$ (see Section 4.2 for use of this property). Garthwaite
and Dickey (1988) describe how sequential elicitation of observation and
process uncertainty in the normal linear model eases the assessment tasks for
experts, where the systematic component is described as the “long run mean”.
For GLMs, an exponential dispersion model has expectation $\mu_{i}$ at
scenario $u_{i}^{\top}$, and so the sample mean $\bar{\psi}$ at scenario
$u_{i}^{\top}$ converges by the law of large numbers to the systematic
component $\mu_{i}$ as $w\rightarrow\infty$. Observation uncertainty
introduced by the random component is therefore absent from the systematic
component, and experts need not correct for observation uncertainty in their
assessments of the systematic component.

###  3.2 Model Structure for Elicitation

Let $\psi$ denote a quantity of interest that arises from a GLM (Section 2.1),
and $k:l=\\{k,k+1,\ldots,l\\}$ denote a sequence of integers. The elicitation
is structured by the prior for GLMs proposed by West (1985),

$\displaystyle\psi_{j(i)}\mid\mu_{i},\phi\sim f(\mu_{i},\phi)\,,\quad
i=1:n\,,\quad j(i)=1:J(i)\,,$ (2a) $\displaystyle g(\mu)=\eta=X\beta+o\,,$
(2b) $\displaystyle\beta\mid\lambda\sim N(\delta,\Sigma/\lambda)\,,$ (2c)
$\displaystyle 1/\phi=\lambda\sim G(s/2,r/2)\,,$ (2d)

where $N(m,S)$ is a normal distribution with mean $m$ and positive definite
covariance matrix $S$, and $G(x\mid a,b)=b^{a}x^{a-1}\exp[-bx]/\Gamma(a)$ is a
gamma distribution.

Eq. (2a) describes the observation uncertainty of the GLM, where $f$ is an
exponential dispersion model, and Eq. (2b) defines the linear model with
$n\times p$ covariate matrix $X$ of rank $p\leq n$ (Section 2.1). Eq. (2c)
describes the process uncertainty by a conditional normal prior for $\beta$.
GLM prior specifications often use a multivariate normal for $\beta$ given
$\phi$ (O’Hagan and Forster, 2004). The normal prior requires specification of
only the first two moments, and it is computationally convenient given a
normal approximation of the likelihood (Gelman et al., 2014). It is sometimes
also supported by asymptotic arguments, as in some $g$-prior settings (see
Section 2.5). Eq. (2d) is the prior for the index parameter $\lambda=1/\phi$
that is, in general, gamma distributed. The case of known dispersion is also
accommodated, as described below.

First, if the dispersion parameter is unknown then $s,r\in(0,\infty)$ in Eq.
(2d) and $\lambda$ is gamma distributed. The marginal prior distribution for
$\beta$ is then generalised multivariate t (Kotz and Nadarajah, 2004, Ch. 5)
with location $\delta$ and scale $\Sigma$,

$St_{n}(\beta\mid\delta,\Sigma,r,s)=\frac{\Gamma((s+n)/2)}{\Gamma(s/2)(\pi
r)^{n/2}}\left|\Sigma\right|^{-\frac{1}{2}}\left[1+\frac{1}{r}(\beta-\delta)^{\top}\Sigma^{-1}\left(\beta-\delta\right)\right]^{-\frac{s+n}{2}}.$
(3)

Although a normal prior for $\beta$ often adequately captures prior
information, West (1985) and O’Hagan and Forster (2004) also note that a
multivariate t prior may be preferred because it allows for extreme values of
$\beta$. Rescaling $\Sigma$ by $\mathbb{E}[\lambda]=s/r$ in Eq. (3) obtains
the equivalent multivariate t distribution denoted by $\beta\sim
t_{n}(\delta,(r/s)\Sigma,s)=St_{n}(\delta,(r/s)\Sigma,s,s)$, see Section S1
for details.

Second, if the dispersion parameter is known then let $s,r\rightarrow\infty$
subject to $s/r=1/\phi$. In this limit, the induced prior for $\beta$ of Eq.
(3) corresponds to Eq. (2c) with $\lambda=s/r$.

###  3.3 Outline of Statistical Elicitation Method

The structured elicitation method is summarised in Procedure 1. Given the
choices of $g$, $U$, and $\mathcal{F}$, the elicitation procedure begins by
selecting, without loss of generality, an exponential dispersion model
$f\in\mathcal{F}$ with variance function $v(\mu)$. The elicitation addresses
the observation uncertainty of the random component in Section 4. The process
uncertainty of the systematic component is addressed in Section 5. During the
elicitation session, the facilitator conditions the expert assessments for the
systematic component on $U$ rather than $X$. After the elicitation session,
the statistician can then explore alternative choices of $X\in\mathcal{X}$
consistent with $U$. To ensure this flexibility, during the elicitation
session the model matrix $X$ is specified as a square full rank matrix. This
non-singular transformation, given the invertible link function $g$, ensures
no loss of information between the elicited systematic component $\mu\mid U$
and the induced prior on $\beta$ (Section 6).

Without loss of generality, the elicitation session in Section 5 further
assumes the saturated model $X=I$ in Eq. (2b). In a slight abuse of notation,
a known offset, if present, is absorbed into the linear predictor $\eta$ such
that $\eta=X\beta$ and $o=0$ in Eq. (2b). The resulting parametrisation for
the linear predictor induced by Eqs. (2b) and (2c) is

$\eta\mid X=I,\lambda\sim N(m,V/\lambda)\,,\quad
V=\operatorname{diag}[V]^{1/2}R\operatorname{diag}[V]^{1/2}\,,$ (4)

where $\operatorname{diag}[V]$ is a diagonal matrix formed from the main
diagonal entries of positive definite matrix $V$ that has positive definite
correlation matrix $R$. The marginal elicitation of the systematic component
described by $m$ and $\operatorname{diag}[V]$ proceeds in Section 5.1.
Elicitation of the correlation matrix $R$ progresses by means of a canonical
vine copula (Section 5.2) as described in Sections 5.3, 5.4 and 5.5.

Data: $\mathcal{F}$, $g$, $U$, $\mathcal{X}$

Result: Estimates $\hat{\delta}$ and $\hat{\Sigma}$, Eq. (2c); estimates
$\hat{s}$ and $\hat{r}$, Eq. (2d)

1 for  _Random component given $f\in\mathcal{F}$_ do

2 Elicit $s$ and $r$ in Eq. (2d) (Section 4 and Procedure 2);

3

4 end for

5for  _Systematic component given $g$ and $U$_ do

6 Elicit $m$ and $\operatorname{diag}[V]$ in Eq. (4) (Section 5.1 and
Procedure 3);

7 Elicit $R$ in Eq. (4) (Sections 5.2–5.5 and Procedure 4);

8

9 end for

10Parametrise $\hat{s}$, $\hat{r}$, $\hat{\delta}$ and $\hat{\Sigma}$ given
$f^{\prime}\in\mathcal{F}$ and $X\in\mathcal{X}$ (Eq. (26) in Section 6);

Return $\hat{s}$, $\hat{r}$, $\hat{\delta}$, $\hat{\Sigma}$

Procedure 1 Overview of elicitation and parametrisation procedures.

After the elicitation session concludes, the priors of Eqs. (2c) and (2d) are
parametrised (Section 6). First, the choice of an exponential dispersion model
$f^{\prime}\in\mathcal{F}$ parametrises the prior for $\lambda$ in Eq. (2d).
Second, the choice of model matrix $X\in\mathcal{X}$ parametrises the prior
for $\beta\mid\lambda$ in Eq. (2c). Choices of $X\in\mathcal{X}$ with either
$p=n$ or $p<n$ are supported, although the latter choice necessarily induces
information loss in the parametrised prior. Extensions of the statistical
elicitation method are presented in Section 7. The method is illustrated by
application that uses an extended GLM in Section 8.

##  4 Random Component

The random component of a GLM captures sources of observation uncertainty not
modelled by the covariates (Section 3.1). The case of known dispersion is
treated in Section 4.1. Section 4.2 presents an elicitation procedure for an
unknown dispersion parameter $\phi$ in Eq. (2) that applies to any exponential
dispersion model. Procedure 2 summarises the structured elicitation for the
random component.

###  4.1 Known Dispersion

If the dispersion parameter $\phi$ in Eq. (2) is known, then let
$s,r\rightarrow\infty$ in Eq. (2) such that $s/r=1/{\phi}$. The marginal
distribution of $\beta$ is then multivariate normal by Eq. (3) and $\eta\sim
N(m,{\phi}V)$ in Eq. (4).

###  4.2 Unknown Dispersion

As noted in Section 1, experts should only be asked to elicit observable
quantities and should not be asked to predict the variance of a response
directly (e.g., Kadane and Wolfson, 1998; Garthwaite and O’Hagan, 2000). The
elicitation of a prior for an unknown dispersion parameter $\phi$ is based on
the predicted sample mean $\bar{\psi}$ of a hypothetical experiment
conditional on the systematic component. This conditioning isolates the
observation uncertainty modelled by the random component of the GLM.

Let the elicitation target follow an exponential dispersion model with
variance function $v(\mu)$. If $\psi_{j}\mid\mu,\phi\sim ED(\mu,\phi)$ for
$j=1:w$ then the sample mean is closed under convolution (Jørgensen, 1987),

$\bar{\psi}=\frac{1}{w}\sum_{j=1}^{w}\psi_{j}\sim
ED\left(\mu,\frac{\phi}{w}\right)\,,$ (5)

and the central limit theorem applies as the sample size $w\rightarrow\infty$
such that

$({\bar{\psi}-\mu})/{\sqrt{\phi/w}}\overset{d}{\rightarrow}N(0,v(\mu))\,.$ (6)

From Eqs. (2d) and (6), the sample mean conditional on a realised systematic
component $\hat{\mu}_{0}\in\Omega$ is for large $w$ approximately univariate
generalised t or equivalently t distributed,

$\bar{\psi}\mid\hat{\mu}_{0},w\sim
St_{1}(\hat{\mu}_{0},v(\hat{\mu}_{0})/w,r,s)=t_{1}(\hat{\mu}_{0},v_{\phi}(\hat{\mu}_{0},w),s)\,,\,v_{\phi}(\hat{\mu}_{0},w)=rv(\hat{\mu}_{0})/(ws)\,,$
(7)

with scale $v_{\phi}(\hat{\mu}_{0},w)$. In Eq. (7), the realised systematic
component $\hat{\mu}_{0}$ is the predicted response for an arbitrary scenario,
denoted say by $u_{0}^{\top}$, in the absence of observation uncertainty.
Theoretically, any value $\hat{\mu}_{0}\in\Omega$ may be chosen by either the
facilitator or the expert. Practically, the expert must be able to hypothesise
a scenario of covariates for which the choice of $\hat{\mu}_{0}$ is plausible.
The conditional elicitation of $\bar{\psi}$ given $\hat{\mu}_{0}$ includes all
and only sources of uncertainty for prediction that are unexplained by the
covariates.

Conditional on a value $\hat{\mu}_{0}$ consistent with a hypothetical
scenario, the elicitation of the sample mean $\bar{\psi}$ given $w$
independent realisations of the target $\psi$ parametrises Eq. (2d). The
expert is asked to provide two central credible intervals of probabilities for
the conditional sample mean $\bar{\psi}$ given $\hat{\mu}_{0}$ and $w$. The
two central credible intervals have probabilities $\alpha_{1}<\alpha_{2}$ with
corresponding lower bounds $d_{1}>d_{2}$ defined by

$P\left(\bar{\psi}\leq
d_{k}\mid\hat{\mu}_{0},w\right)=P\left(\bar{\psi}>2\hat{\mu}_{0}-d_{k}\mid\hat{\mu}_{0},w\right)=\frac{1-\alpha_{k}}{2}\,,\quad
k\in\\{1,2\\}\,,\quad\alpha_{1}<\alpha_{2}\,.$

From Eq. (7), the quantile function for the sample mean is

$F^{-1}(q\mid\hat{\mu}_{0},w,r,s)=\hat{\mu}_{0}+\sqrt{\frac{rv(\hat{\mu}_{0})}{ws}}T_{1}^{-1}(q\mid
s)\,,\quad q\in(0,1)\,,$ (8)

where $T_{1}^{-1}(\cdot\mid s)$ is the quantile function of the univariate
Student’s t distribution with $s$ degrees of freedom, location of zero and
scale of one. Given the choices of $d_{k}$ and $\alpha_{k}$ for
$k\in\\{1,2\\}$, the elicited parameter $s$ is obtained as the solution to the
relation

$\frac{d_{1}-\hat{\mu}_{0}}{d_{2}-\hat{\mu}_{0}}=\frac{F^{-1}\left(\frac{1-\alpha_{1}}{2}\mid\hat{\mu}_{0},w,r,s\right)-\hat{\mu}_{0}}{F^{-1}\left(\frac{1-\alpha_{2}}{2}\mid\hat{\mu}_{0},w,r,s\right)-\hat{\mu}_{0}}=\frac{T_{1}^{-1}\left(\frac{1-\alpha_{1}}{2}\mid
s\right)}{T_{1}^{-1}\left(\frac{1-\alpha_{2}}{2}\mid
s\right)}<\frac{\Phi^{-1}\left(\frac{1-\alpha_{1}}{2}\right)}{\Phi^{-1}\left(\frac{1-\alpha_{2}}{2}\right)}\,,$
(9)

where $\Phi^{-1}(\cdot)$ is the quantile function of the limiting standard
normal as $s,r\rightarrow\infty$ in Eq. (8). The elicited bounds $d_{1}$ and
$d_{2}$ are held consistent with the generalised t approximation by the upper
bound in Eq. (9), see also Kadane et al. (1980) for application of this
relation to the elicitation of the degrees of freedom parameter in the normal
linear model. Given the degrees of freedom $s$ and the elicited quantiles
$d_{1}$ and $d_{2}$, the rate parameter $r$ is obtained from Eqs. (7) and (8)
by

$v_{\phi}(\hat{\mu}_{0},w)=\left(\frac{d_{k}-\hat{\mu}_{0}}{T_{1}^{-1}([1-\alpha_{k}]/2\mid
s)}\right)^{2}\quad\implies\quad
r=\frac{v_{\phi}(\hat{\mu}_{0},w)ws}{v(\hat{\mu}_{0})}\,,\quad
k\in\\{1,2\\}\,.$ (10)

The choices of $\alpha_{1}$ and $\alpha_{2}$ should be sufficiently different
to elicit both the degrees of freedom $s$ from Eq. (9), which controls the
heaviness of the tails, and the scale $v_{\phi}(\hat{\mu}_{0},w)$ from Eq.
(10).

Several advantages of this elicitation approach for unknown $\phi$ are now
summarised:

  1. 1.

the approximate density function of the sample mean is unimodal and symmetric,
which agrees with general expectations often expressed by experts (Winkler,
1967);

  2. 2.

the approximating t distribution can be made arbitrarily close to the
distribution of the sample mean under $f$ in Eq. (2a) by choice of $w$ large
enough, and the quality of this approximation is assessable;

  3. 3.

the completed elicitation can be used to support prior parametrisations for
alternative observation models (Section 6); and

  4. 4.

the same elicitation strategy for unknown dispersion applies to any
exponential dispersion model, whether the support is discrete, continuous or
even semi-continuous (see Section 7.2 for an example of the latter).

Techniques for evaluating the quality of the t approximation are discussed
below.

A check on the expert’s understanding of the amount of information in a sample
(Section 3.1) can be assessed by the facilitator through trialling different
values of $w$, and confirming that the expert correspondingly adjusts the
predictions for the sample mean. The facilitator may also wish to guide
choices of finite $w$ and $\hat{\mu}_{0}$ that are robust to the CLT
approximation of Eq. (6). Based on a facilitator’s point estimate for $\phi$,
a bound on the approximation error of Eq. (6) is available through the Berry-
Esseen inequality,

$\sup_{z}\left|P\left(\frac{\sqrt{w}(\bar{\psi}-\hat{\mu}_{0})}{\sigma}\leq
z\right)-\Phi(z)\right|\leq\frac{K\breve{\mu}^{3}}{\sqrt{w}}\leq\frac{K\sqrt{\breve{\mu}^{4}}}{\sqrt{w}}\,,$
(11)

where $\sigma=\sqrt{\phi v(\hat{\mu}_{0})}$ and
$\breve{\mu}^{k}=\mathbb{E}\left[\left|\psi-\hat{\mu}_{0}\right|^{k}\right]/\sigma^{k}$
for $k\in\\{3,4\\}$. A recent estimate is $K=0.469$ (Shevtsova, 2017). The
upper bound on the Berry-Esseen inequality in Eq. (11) uses the more
accessible kurtosis $\breve{\mu}^{4}$ or excess kurtosis $\breve{\mu}^{4}-3$
(see Section S2.1 for details). This upper bound may assist the facilitator in
suggesting values for $\hat{\mu}_{0}$ and $w$ that improve the generalised t
approximation for the sample mean.

The discrepancy of the sample mean distribution with respect to the
approximating distribution may also be considered. Two options are provided
that use Monte Carlo samples of $\lambda$ from Eq. (2d) followed by
composition sampling for the sample mean from Eq. (5). The latter step is
conditional on $\lambda$ and so only the ability to sample from natural
exponential families is required. First, the distribution function of the
sample mean ${F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)$ is estimated with
accompanying confidence intervals for the Monte Carlo error (see Section
S2.2). The estimated distribution function of the sample mean,
$\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)$, and its generalised t
approximation are then compared. A numerical measure of discrepancy is
provided by the Kolmogorov distance,

$\sup_{x}\left|\hat{F}(x\mid\hat{\mu}_{0},w,s,r)-T_{1}\left(\frac{x-\hat{\mu}_{0}}{\sqrt{v_{\phi}(\hat{\mu}_{0},w)}}\mid
s\right)\right|\,.$ (12)

Second, the expected logarithm of the posterior odds against the generalised t
approximation is estimated (see Section S2.3). Assuming equal prior odds, this
measure is equivalent to the Kullback-Leibler divergence (Kullback and
Leibler, 1951), or expected log Bayes factor, for the joint distribution of
the sample mean and dispersion relative to the approximation from Eq. (6),

$D_{\bar{\psi},\lambda}=\mathbb{E}_{\bar{\psi},\lambda\mid\hat{\mu}_{0},w,s,r}\left[\log\frac{f(\bar{\psi}\mid\hat{\mu}_{0},\frac{1}{w\lambda})G(\lambda\mid\frac{s}{2},\frac{r}{2})}{N(\bar{\psi}\mid\hat{\mu}_{0},\frac{1}{w\lambda})G(\lambda\mid\frac{s}{2},\frac{r}{2})}\right]\,.$
(13)

Jeffreys (1961, Appendix B) and Kass and Raftery (1995) suggest thresholds for
the acceptance or rejection of hypotheses based on the Bayes factor. Eqs. (12)
and (13) can be assessed before or after an expert has contributed bounds
$d_{1}$ and $d_{2}$. In the former case, the facilitator may suggest $w$ based
on experimentation for a plausible choice of $s$ and $r$. In the latter case,
if the elicited discrepancies are unacceptably large then the elicitation can
be iteratively reassessed for the unknown dispersion (Procedure 2).

##  5 Systematic Component

The elicitation of the systematic component and associated process uncertainty
conditions on the scenarios documented by $U$ (Section 3.1). The $n\times 1$
systematic component $\mu$ is indexed by scenario, such that $\mu_{i}=\mu\mid
u_{i}^{\top}$, where $u_{i}^{\top}$ is the $i$th row of $U$ that defines the
$i$th scenario for $i=1:n$. Elicitation for the systematic component begins
with the marginal distributions of $\mu$ (Section 5.1). The dependencies among
scenarios are then modelled with a canonical vine copula (Section 5.2). The
scenarios additionally condition on hypothetical values of the systematic
component $\hat{\mu}_{1:l}$ realised at the first $l$ scenarios for $1\leq
l<n$ (Section 5.3). At each successive tree level $\mathcal{T}_{l}$ in the
canonical vine, the expert elicits the conditional distribution of the
systematic component $\mu_{k}\mid\hat{\mu}_{1:l}$ for $k=(l+1):n$. The
canonical vine copula model for dependence is progressively parametrised with
conditional median elicitations such that the expert is never asked to forget
hypothetical realised data once admitted into the elicitation session, which
accords with the design principles of the elicitation method (Section 1). The
choice of the conditioning values $\hat{\mu}_{1:l}$ for the systematic
component is considered in Section 5.4. To optionally reduce the number of
conditional elicitations, $t$-truncated canonical vine dependence models may
be used (Section 5.5). A $t$-truncated vine assumes no dependence on
hypothetically realised values of the systematic component above tree level
$\mathcal{T}_{t}$.

###  5.1 Elicitation of Marginal Distributions

Given the marginal prior of Eq. (3) and the specification $X=I$ in Eq. (4),
the linear predictor is distributed as

$\eta\sim St_{n}(m,V,r,s)$ (14)

with $n\times 1$ location vector $m$ and $n\times n$ positive definite scale
matrix $V$. The parameters $s$ and $r$ are known from the elicitation of the
random component (Section 4). Given also the link function $g$ in Eq. (2b),
the density function for the systematic component is

$f(\mu\mid m,V,r,s)=St_{n}\left(g(\mu)\mid
m,V,r,s\right)\,|dg(\mu)/d\mu|\,,\quad\mu\in\Omega\,.$ (15)

From Eq. (15), expressions for the univariate density, distribution and
quantile functions of the systematic component at each scenario in $U$ are
available to provide graphical and numerical feedback to the expert, see Eqs.
(S2), (S3) and (S4) in Section S3.1.

The central credible intervals $(a,b)$ for $\mu=g^{-1}(\eta)$ defined by

$P\left(\mu_{i}\leq
a_{i}\right)=P\left(\mu_{i}>b_{i}\right)=\frac{1-\alpha}{2}\,,\quad(a_{i},b_{i})\in\Omega\,,\quad
i=1:n\,,\quad 0<\alpha<1\,,$

are elicited, where $P\left(\mu_{i}\leq z\right)=F_{\mu}(z\mid
m_{i},V_{ii},r,s)$ is the marginal distribution function of the systematic
component at elicitation scenario $u_{i}^{\top}$. The location and scale are

$m_{i}=\frac{g(a_{i})+g(b_{i})}{2}\,,\quad
V_{ii}=\left(\frac{g(b_{i})-m_{i}}{T_{1}^{-1}((1+\alpha)/2\mid
s)}\right)^{2}\frac{s}{r}\,,$

because $g$ is invertible. The elicited central credible interval
$(a_{i},b_{i})$ for the systematic component at scenario $u_{i}^{\top}$
uniquely parametrises the required marginal distribution of the linear
predictor $\eta_{i}$ for $i=1:n$.

The expert iteratively updates each marginal prior by adjustment of the
central credible interval $(a_{i},b_{i})$ until the relevant probabilities and
quantiles of interest at scenario $u_{i}^{\top}$, along with any graphical
feedback, adequately approximate the expert’s beliefs given the model
assumptions. This elicitation protocol is repeated for each elicitation
scenario $u_{i}^{\top}$, $i=1:n$, to induce the $n$ marginal priors for the
underlying linear predictor $\eta$. The above elicitation steps are summarised
in Procedure 3.

###  5.2 Dependence and Canonical Vine Copula

The conditional distribution of the linear predictor
$\eta_{l:k}\mid\eta_{1:(l-1)}$ for $1\leq l\leq k\leq n$, where
$\eta_{1:0}=\emptyset$, has the generalised t distribution (Kotz and
Nadarajah, 2004, Ch. 5),

$\eta_{l:k}\mid\eta_{1:(l-1)}\sim St_{k-l+1}\left(m_{l:k\mid
1:(l-1)},V_{l:k,l:k\mid 1:(l-1)},r+\zeta(\eta_{1:(l-1)}),s+l-1\right)$ (16)

with conditional location vector

$m_{l:k\mid
1:(l-1)}=m_{l:k}+V_{l:k,1:(l-1)}V_{1:(l-1),1:(l-1)}^{-1}\left(\eta_{1:(l-1)}-m_{1:(l-1)}\right)$
(17)

and conditional scale matrix

$V_{l:k,l:k\mid
1:(l-1)}=V_{l:k,l:k}-V_{l:k,1:(l-1)}V_{1:(l-1),1:(l-1)}^{-1}V_{1:(l-1),l:k}\,,$
(18)

where
$\zeta(\eta_{1:(l-1)})=(\eta_{1:(l-1)}-m_{1:(l-1)})^{\top}V_{1:(l-1),1:(l-1)}^{-1}(\eta_{1:(l-1)}-m_{1:(l-1)})$.
Conditional on the index parameter $\lambda$, the linear predictor $\eta$ is
multivariate normal with mean $m$ and covariance $V/\lambda$ by Eq. (4). The
conditional correlation (Joe, 2015, Ch. 2) for this multivariate normal
distribution is

$\displaystyle\rho_{l,k|1:(l-1)}=$
$\displaystyle\left({V_{k,l}-V_{l,1:(l-1)}V_{1:(l-1),1:(l-1)}^{-1}V_{1:(l-1),k}}\right)$
$\displaystyle\times\left({V_{l,l}-V_{l,1:(l-1)}V_{1:(l-1),1:(l-1)}^{-1}V_{1:(l-1),l}}\right)^{-1/2}$
(19)
$\displaystyle\times\left({V_{k,k}-V_{k,1:(l-1)}V_{1:(l-1),1:(l-1)}^{-1}V_{1:(l-1),k}}\right)^{-1/2},\quad
1\leq l\leq k\leq n\,.$

For $k\neq l$, Eq. (19) is equivalent to the partial correlation of variates
$\eta_{k}$ and $\eta_{l}$ conditional on $\lambda$ and the realised variates
$\eta_{1:(l-1)}$. If $l=1$ then Eq. (19) reduces to the marginal correlation
between $\eta_{k}$ and $\eta_{l}$, given $\lambda$, that is defined by
$\rho_{l,k}=V_{k,l}/\sqrt{V_{l,l}V_{k,k}}$.

Bedford and Cooke (2002) use a specific selection of $\binom{n}{2}$ partial
correlations to develop a regular vine copula that defines a positive definite
correlation matrix $R$, such that each partial correlation has range $(-1,1)$.
There are various formulations of regular vines that each provide a different
parametrisation of an equivalent positive definite correlation matrix $R$
(Joe, 2015, Ch. 3). On example is the _canonical_ vine copula (Bedford and
Cooke, 2001) that can be represented by an upper triangular matrix of partial
correlations,

$\mathrm{P}=\begin{bmatrix}1&\rho_{1,2}&\rho_{1,3}&\ldots&\rho_{1,n-1}&\rho_{1,n}\\\
&1&\rho_{2,3|1}&\ldots&\rho_{2,n-1|1}&\rho_{2,n|1}\\\ &&&\ddots&&\vdots\\\
&&&&1&\rho_{n-1,n|1:(n-2)}\\\ 0&&&&&1\\\ \end{bmatrix}.$

Each row of the matrix $\mathrm{P}$ corresponds to a level of dependency,
denoted by $\mathcal{T}_{l}$ and referred to as the tree level (Bedford and
Cooke, 2001, 2002; Joe, 2015). The tree level $\mathcal{T}_{l}$ is the set of
conditional correlations between $\eta_{l}$ and each member of the set
$\eta_{(l+1):n}$ conditioned on the set $\eta_{1:(l-1)}$. The canonical vine
has the same conditioning variables within a tree level, and exactly one
variable is added to the conditioning set when moving up one level (Bedford
and Cooke, 2001). Joe (2015) provides an algorithm for the mapping
$\mathrm{P}\rightarrow R$ (see Algorithm S3.2 in Section S3.2).

The canonical vine is particularly advantageous for conditional elicitation.
The canonical vine ensures that the expert is never asked to forget
hypothetical realised data admitted into the elicitation session, in accord
with the design principles of the elicitation method (Section 1). Regular
vines other than the canonical vine adopt a different partial correlation
structure $\mathrm{P}$ and introduce more than one conditioning variable at
lower tree levels (Joe, 2015, Ch. 3). On the opposite extreme from the
canonical vine is the $D$-vine, for example, where the number of unique
conditioning variables at each tree level is maximised within the regular vine
constraints. For regular vines other than the canonical vine, the expert is
asked at some point to forget conditioning values from prior elicitations,
only to have these previously considered conditioning values reintroduced at
the same tree level or at higher tree levels. This complication is avoided by
means of the canonical vine copula.

###  5.3 Elicitation of Canonical Vine Copula

At tree level $\mathcal{T}_{l}$, let $c_{k|1:l}$ denote for $k>l$ the
conditional median of the systematic component at scenario $u_{k}^{\top}$
given hypothetical values $\hat{\mu}_{1:l}=g^{-1}(\hat{\eta}_{1:l})$ of the
systematic component realised at the first $l$ scenarios. The conditional
median is defined by

$P\left(\mu_{k}\leq c_{k\mid
1:l}\left.\mid\hat{\mu}_{1:l}\right.\right)=0.5\,,\quad 1\leq l<k\leq
n\,,\quad c_{k\mid 1:l}\in\Omega\,.$

The sequence of elicitations for conditional medians progressively elicits the
upper triangular entries of the array

$C=\begin{pmatrix}c_{1}&c_{2\mid 1}&c_{3\mid 1}&\ldots&c_{n-1\mid 1}&c_{n\mid
1}\\\ &c_{2}&c_{3\mid 1:2}&\ldots&c_{n-1\mid 1:2}&c_{n\mid 1:2}\\\
&&&\ddots&&\vdots\\\ &&&&c_{n-1}&c_{n\mid 1:(n-1)}\\\ 0&&&&&c_{n}\\\
\end{pmatrix}\,.$

The elicitation sequence begins in the top row for tree level
$\mathcal{T}_{1}$ and continues until the second to bottom row for tree level
$\mathcal{T}_{n-1}$. Within each row (tree level) the conditional elicitations
may be completed in any order, but would typically be conducted in sequence.
The main diagonal entries of $C$ are designated as the elicitations at tree
level $\mathcal{T}_{0}$, where the marginal medians of the systematic
component are $\operatorname{diag}[C]=g^{-1}(m)$ because $g$ is invertible.
The vector $m$ was previously elicited in Section 5.1. Moreover, the
conditional median of the systematic component is $m_{k|1:l}=g(c_{k\mid
1:l})$.

Given hypothetically realised values of the systematic component
$\hat{\mu}_{1:l}=g^{-1}(\hat{\eta}_{1:l})$, Eq. (17) with correspondingly
adjusted indices shows that the difference between the conditional and
marginal medians for $\eta_{k}$ is

$\displaystyle m_{k|1:l}-m_{k}$
$\displaystyle=V_{k,1:l}V_{1:l,1:l}^{-1}\left(\hat{\eta}_{1:l}-m_{1:l}\right)\,,\quad
1\leq l<k\leq n\,.$ (20)

Let $h_{1:l}=V_{1:l,1:l}^{-1}\left(\hat{\eta}_{1:l}-m_{1:l}\right)$. Solving
for $V_{k,l}$ in Eq. (20) with realisations for the linear predictor such that
$\hat{\eta}_{1:l}\neq m_{1:l}$ then obtains

$V_{k,l}=\begin{cases}\left(m_{k\mid 1}-m_{k}\right)/h_{1}&\textrm{ if
}l=1\,,\\\ \left[\left(m_{k\mid
1:l}-m_{k}\right)-V_{k,1:(l-1)}h_{1:(l-1)}\right]\big{/}h_{l}&\textrm{
otherwise.}\end{cases}$ (21)

At tree level $\mathcal{T}_{l}$, the entries of $m$, $V$ and $h$ that appear
on the right hand side of Eq. (21) have already been derived from elicitations
completed for the marginal elicitations (Section 5.1) and conditional medians
elicited in $C$ at lower tree levels $\mathcal{T}_{1:(l-1)}$. Apart from the
conditional median $m_{k|1:l}$ that is the current focus of elicitation, there
are no other unknowns on the right hand side of Eq. (21) at tree level
$\mathcal{T}_{l}$.

Given $V_{k,l}$, the scale of $\eta_{k}$ conditional on $\hat{\eta}_{1:l}$ is

$V_{k,k\mid 1:l}=V_{k,k}-V_{k,1:l}V_{1:l,1:l}^{-1}V_{1:l,k}\,,$ (22)

which follows from Eq. (18) with indices adjusted to account for conditioning
on $\hat{\eta}_{1:l}$. From the identity $V_{k,k\mid 1:l}=V_{k,k\mid
1:(l-1)}(1-\rho_{l,k\mid 1:(l-1)}^{2})$ (Whittaker, 1990, Ch. 5) and the
previously elicited conditional scale entry $V_{k,k\mid 1:(l-1)}$ at the
preceding tree level $\mathcal{T}_{l-1}$, the constraint $|\rho_{l,k\mid
1:(l-1)}|<1$ (Section 5.2) requires

$0<V_{k,k\mid 1:l}\leq V_{k,k\mid 1:(l-1)}\,.$ (23)

The elicited conditional scale for $\eta_{k\mid 1:l}$ at the current tree
level $\mathcal{T}_{l}$ must be positive and cannot exceed the conditional
scale for $\eta_{k\mid 1:(l-1)}$ previously determined at tree level
$\mathcal{T}_{l-1}$. The consistency check of Eq. (23) is available for
feedback to the expert at tree level $\mathcal{T}_{l}$. Feasible bounds on the
conditional median $c_{k|1:l}=g^{-1}(m_{k|1:l})$ may be graphically or
numerically presented to the expert at the start of each elicitation (Section
S3.2).

At each tree level $\mathcal{T}_{l}$, the completed elicitation
$c_{k|1:l}=g^{-1}(m_{k\mid 1:l})$ induces the scale entry $V_{k,l}$ from Eq.
(21), which having passed the consistency check of Eq. (23) then defines the
partial correlation $\mathrm{P}_{lk}=\rho_{l,k\mid 1:(l-1)}$ by Eq. (19). Thus
the successive elicitation of each row of entries along the upper triangular
part of $C$ completes the corresponding row in $\mathrm{P}$. The univariate
conditional density, distribution and quantile functions of the systematic
component are also available to provide graphical and numerical feedback to
the expert, see Eqs. (S5), (S6) and (S7) in Section S3.1. The conditional
elicitation steps are summarised in Procedure 4. The choice of conditioning
hypothetical realisations $\hat{\mu}_{1:l}$ is discussed in Section 5.4.

###  5.4 Conditioning Values for Systematic Component

In Section 5.3, the conditioning values of the systematic component by Eq.
(21) need only be specified such that $g(\hat{\mu}_{j})=\hat{\eta}_{j}\neq
m_{j}$ with $\hat{\mu}_{j}\in\Omega$ for $j=1:l$. From the expert’s
perspective, however, it is important to choose conditioning values for the
target that are not overly surprising or extreme so as to avoid severe or
uneven constraints on the conditional scale parameters elicited at different
tree levels by Eq. (23).

One choice for the conditioning values of the systematic component is as
follows. At the start of elicitations for tree level $\mathcal{T}_{l}$, set
$\hat{\mu}_{l}$ equal to one of the bounds of the conditional central credible
interval $(a_{l\mid 1:(l-1)},b_{l\mid 1:(l-1)})$ of probability $\alpha$
defined by the elicited distribution of the systematic component from the
previous tree level $\mathcal{T}_{l-1}$. Choosing the same probability level
$\alpha$ as in the marginal elicitations (Section 5.1) efficiently builds on
the expert’s earlier experience at previous tree levels. The conditioning
values reflect the non-increasing, and typically decreasing, uncertainty in
the systematic component with increasing tree level by Eq. (23). The selection
of the upper versus lower bound from the conditional central credible interval
for the hypothetical conditioning variable may at each tree level be chosen by
the facilitator or the expert, or selected systematically or randomly.

The above approach avoids outliers and is always applicable, and so it is a
reasonable default method for the choice of conditioning values. But the
common choice of $\alpha$ in the above approach is not the only possibility.
An alternative choice for the realised values is demonstrated in the
elicitation case study of Section 8 given unknown dispersion for a bounded
elicitation target. In that example, the conditioning values $\hat{\mu}_{1:l}$
are selected from the bounds of conditional credible intervals of probability
$\alpha$ induced by Eq. (4) with fixed dispersion of one. The magnitude of the
differences between the conditioning values and elicited medians are then
determined by the elicited scale $V$ for the linear predictor, which screens
out noise introduced by the unknown dispersion.

The elicitation for the canonical vine copulas sequentially builds on
conditioning values, and thereby avoids scenarios where an expert is asked at
some point to forget conditioning values from prior elicitations (Section 1).
This sequential aspect means that the ordering of scenarios can also be an
important consideration for an elicitation session. The ordering of
conditional elicitations is discussed further in the context of truncated vine
dependence models (Section 5.5).

###  5.5 Truncated Canonical Vine Dependence Model

The $t$-truncated canonical vine replaces dependencies at tree levels greater
than $t$ with conditional independence assumptions (Joe, 2015). Truncation
after tree level $\mathcal{T}_{t}$ sets conditional correlations to zero in
upper tree levels $\mathcal{T}_{(t+1):(n-1)}$. The truncated partial
correlation matrix $\mathrm{P}(t)$ has entries

$\mathrm{P}_{i,j}(t)=\begin{cases}0&\textrm{if }i>t\textrm{ and }j>i\\\
\mathrm{P}_{i,j}&\textrm{otherwise}\end{cases}\,,\quad i,j=1,\ldots n\,,\quad
t\in\left\\{0,\ldots,n-1\right\\}\,.$

Application of Algorithm S3.2 to $\mathrm{P}(t)$ returns the corresponding
truncated correlation matrix $R(t)$. The truncated elicited scale matrix of
Eq. (4) is then

$V(t)=\operatorname{diag}[V]^{1/2}R(t)\operatorname{diag}[V]^{1/2}\,.$

The diagonal entries of $V(t)$ were previously elicited by the marginal
elicitations of Section 5.1 such that
$\operatorname{diag}[V(t)]=\operatorname{diag}[V]$ for any choice of
truncation. Procedure 4 allows for truncation during the elicitation session.

Truncation after level $\mathcal{T}_{0}$ means that the elicitation for the
systematic component need not progress beyond Section 5.1. If the dispersion
parameter is known then the induced prior is an ICMP (Section 2.4). Truncation
after level $\mathcal{T}_{n-1}$ corresponds to no truncation, that is,
$\mathrm{P}=\mathrm{P}(n-1)$, $R=R(n-1)$ and $V=V(n-1)$. A truncated vine
ignores dependencies that would be elicited by considering additional
realisations of the systematic component at higher levels.

Kurowicka (2011) suggests placing high absolute (partial) correlations in low
tree levels so that smaller absolute correlations occur in the higher tree
levels that are truncated. The sequence of introduced conditioning values
during the elicitation can therefore be important. Those scenarios expected to
have high dependency with other scenarios should have conditioning values
introduced at lower tree levels early in the elicitation session to minimise
the loss of information induced by a truncated dependence model. Depending on
time and resources the elicitation session may be stopped before all tree
levels are completed, and the truncated dependence model then invoked.
Nevertheless, it should be noted that dependency information may yet be
present at higher tree levels, which would not be revealed without conducting
the corresponding conditional elicitations. Truncated vine dependence models
are evaluated in the elicitation case study (Section 8).

##  6 Induced Prior

At the conclusion of the elicitation session, the facilitator has elicited the
set of parameters $\theta=\left\\{s,r,m,V\right\\}$ (Section 3 and Procedure
1). The parameters $s$ and $r$ are elicited in Section 4. The location $m$ and
scale $V$ are elicited in Section 5 for the systematic component $\mu$
conditional on $s$, $r$, link function $g$ and elicitation scenarios $U$. The
elicited scale matrix $V$ may be truncated at level $t\in\\{0,\ldots,n-1\\}$
in Section 5.5; formally, set $V=V(t)$ in the set $\theta$. Let $H$ denote the
model described by Eqs. (2) and (4) with elicited parameters $\theta$.

Section 3 allows for alternative choices of $n\times p$ model matrix
$X\in\mathcal{X}$ of full column rank or exponential dispersion model
$f^{\prime}\in\mathcal{F}$. Let $H^{\prime}$ denote an alternative model
described by the choice of exponential dispersion model $f^{\prime}$ with
variance function $v^{\prime}(\mu)$ in Eq. (2a), $X\in\mathcal{X}$ in Eq.
(2b), $\delta^{\prime}$ and $\Sigma^{\prime}$ in Eq. (2c), and $s^{\prime}$
and $r^{\prime}$ in Eq. (2d). If the dispersion is known, with value
$\phi=r/s$ under $H$ in Section 4.1, let $\phi^{\prime}=r^{\prime}/s^{\prime}$
denote its value under $H^{\prime}$. If the dispersion is unknown, then
consistency with the elicited scale $v_{\phi}(\hat{\mu}_{0},w)$ of the sample
mean distribution conditional on $\hat{\mu}_{0}$ and $w$ in Section 4.2
requires

$v_{\phi}(\hat{\mu}_{0},w)=\frac{rv(\hat{\mu}_{0})}{ws}=\frac{r^{\prime}v^{\prime}(\hat{\mu}_{0})}{ws^{\prime}}\,.$
(24)

The parameter set of the alternative model $H^{\prime}$ is
$\theta^{\prime}=\\{s^{\prime},r^{\prime},\delta^{\prime},\Sigma^{\prime}\\}$.

The estimation procedure minimises the the Kullback-Leibler divergence or mean
information between $H$ and $H^{\prime}$ per observation from $H$ for $\mu$
and $\lambda$,

$\displaystyle D(H:H^{\prime})=$
$\displaystyle\int\log\frac{f\left(\mu\mid\lambda,\theta,A\right)G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)}{f\left(\mu\mid\lambda,\theta^{\prime},X\right)G\left(\lambda\mid\frac{s^{\prime}}{2},\frac{r^{\prime}}{2}\right)}f\left(\mu\mid\lambda,\theta,A\right)G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)d\mu
d\lambda$ $\displaystyle=$ $\displaystyle\int\log\frac{N\left(\eta\mid
Am,\frac{r}{s\lambda}AVA^{\top}\right)G\left(\lambda\mid\frac{s}{2},\frac{s}{2}\right)}{N\left(\eta\mid
X\delta^{\prime},\frac{r^{\prime}}{s^{\prime}\lambda}X\Sigma^{\prime}X^{\top}\right)G\left(\lambda\mid\frac{s^{\prime}}{2},\frac{s^{\prime}}{2}\right)}$
(25) $\displaystyle\times N\left(\eta\mid
Am,\frac{r}{s\lambda}AVA^{\top}\right)G\left(\lambda\mid\frac{s}{2},\frac{s}{2}\right)d\eta
d\lambda\,,$

where $A$ is a projection matrix from $H$ onto the column space
$\mathbf{C}(X)$ of $X$ such that $P\left(g(\mu)\in\mathbf{C}(X)\right)=1$. The
last equality follows from the invariance property of information (Kullback,
1959, Ch. 2) given the nonsingular transformation $g:\mu\rightarrow\eta$. The
minimisation of Eq. (25) is with respect to $\theta^{\prime}$ and the
projection matrix $A$. The divergence $D\left(H:H^{\prime}\right)$ is
minimised by the parameter set
$\left\\{\hat{\theta},\hat{A}\right\\}=\mathop{\mathrm{arg\,min}}_{\left\\{\theta^{\prime},A\right\\}}D\left(H:H^{\prime}\right)$
(see Section S4 for details), where
$\hat{\theta}=\left\\{\hat{s},\hat{r},\hat{\delta},\hat{\Sigma}\right\\}$ and

$\begin{gathered}\hat{\delta}=\left(X^{\top}V^{-1}X\right)^{-1}X^{\top}V^{-1}m\,,\quad\quad\hat{\Sigma}=q^{-1}\left(X^{\top}V^{-1}X\right)^{-1}\,,\\\
\hat{A}=X\left(X^{\top}V^{-1}X\right)^{-1}X^{\top}V^{-1}\,,\quad\hat{s}=s\,,\quad\hat{r}=rq\,,\\\
q=\begin{cases}{\phi^{\prime}}/{{\phi}}&\textrm{if dispersion is known}\\\
{v(\hat{\mu}_{0})}/{v^{\prime}(\hat{\mu}_{0})}&\textrm{if dispersion is
unknown.}\end{cases}\end{gathered}$ (26)

If the model matrix $X$ is non-singular then $\hat{A}=I$ such that no
information elicited from the expert is lost as the prior is constructed for
the unknown coefficients (see Section S4). An optimal projection $\hat{A}$ is
induced otherwise.

##  7 Extensions

The generality of statistical models supported by the elicitation method is
illustrated with examples for multivariate responses (Section 7.1), semi-
continuous GLMs with point mass at zero (Section 7.2), extended GLM on the
unit interval (Section 7.3), and overdispersed and zero-inflated count data
(Section 7.4).

###  7.1 Multivariate Modelling

Consider $J$ multivariate observations $y_{j}=[y_{j,1},\ldots,y_{j,d+1}]$ of
$d+1$ categories, where for $j=1:J$ the $k$th entry of $y_{j}$ take the value
one for category $k$ with probability $\mathbb{E}[y_{j,k}]=p_{j,k}>0$ and is
zero otherwise. The probabilities obey the constraints
$\sum_{k=1}^{d}p_{j,k}<1$ and $p_{j,d+1}=1-\sum_{k=1}^{d}p_{j,k}$. If the
$d+1$ categories are unordered then the observations follow a multinomial
distribution, which is an exponential dispersion model with known sample size
$\lambda$ (Jørgensen, 1987). Elicitation for multinomial responses is of
interest (Elfadaly and Garthwaite, 2017; Wilson, 2018). Elfadaly and
Garthwaite (2020), for example, adapt the technique of Kadane et al. (1980) to
elicit an additive logistic normal prior for the unknown $\beta$ of a
multinomial GLM. The procedure proposed by Elfadaly and Garthwaite (2020)
elicits the odds of success $p_{j,k}/p_{j,d+1}$ for each category
$k\in\\{1:d\\}$ against a reference category $d+1$. Specification of a
multivariate normal distribution on the log odds induces a logistic normal
distribution (Aitchison and Shen, 1980) on the $d$ dimensional positive
simplex of $p_{j,1:d}$.

The statistical elicitation method presented here applies to the logistic
normal prior, where the scenarios accommodate the multivariate GLM structure.
Identify $d$ rows of the scenario matrix $U$ and the corresponding design
points of $X$ with each of the $J$ observations. Let $u(j,k)_{i}^{\top}$ and
$x(j,k)_{i}^{\top}$ denote dependence of scenario $i\in\\{1:n\\}$ on category
$k\in\\{1:d\\}$ and the set of covariates associated with observation $y_{j}$.
The linear predictor of the multinomial GLM is then
$\log(p_{j,k}/p_{j,d+1})=x(j,k)_{i}^{\top}\beta=\eta_{i}$. Fahrmeir and Tutz
(1994, Ch. 3) describe possible design matrices for multinomial GLMs, see also
Elfadaly and Garthwaite (2020). Define the elicitation target
$\psi(j,k)=p_{j,k}/p_{j,d+1}$ for scenario $u(j,k)_{i}^{\top}$. Given
logarithmic link function $g$ and dispersion ${\phi}=1$, Sections 4.1 and 5
then provide the elicitation procedures for an additive logistic normal prior.

The elicitation method presented here avoids asking experts to provide
estimates of conditional central credible intervals (Section 5.3), as
specified in the design principles of Section 1. Also, the method allows
reduction of the number of elicited quantities via truncated canonical vine
dependence models (Section 5.5). A logistic normal prior for the multinomial
distribution without covariates is obtained as a special case, where the model
matrix is the $d\times d$ identity matrix $X=I$. Other choices for the
logistic transformation can be accommodated. For example, the choice
$\psi(j,k)=p_{j,k}/(1-\sum_{k^{\prime}=1}^{k}p_{j,k^{\prime}})$ with log link
$g$ corresponds to a multiplicative logistic normal distribution (Aitchison,
1982). See O’Hagan and Forster (2004, Ch. 12) for characterisation of
multivariate normal priors given log contrasts applied to the multinomial
distribution.

If the $d+1$ categories are ordered then the data are ordinal. The sequential
model or continuation ratio model for ordinal data (Fahrmeir and Tutz, 1994)
is defined by $P(y_{j}=k\mid y_{j}\geq k,x(j,k)_{i}^{\top})=h(\eta_{i})$,
where $h$ is a binary response function and $y_{j}=k$ if $y_{j,k}=1$. The
sequential model is equivalent to a univariate GLM with unconstrained $\beta$
and known dispersion, where the likelihood is formed from conditional binary
responses defined by the $d+1$ categories (e.g., Tutz, 2012, Ch. 9). The
procedures of Sections 4.1 and 5 therefore apply. For example, the sequential
logit model (Fahrmeir and Tutz, 1994) is defined by $\log[P(y_{j}=k\mid
x(j,k)_{i}^{\top})/P(y_{j}>k\mid x(j,k)_{i}^{\top})]=\eta_{i}$. The
elicitation method progresses for target
$\psi(j,k)=p_{j,k}/(1-\sum_{k^{\prime}=1}^{k}p_{j,k^{\prime}})$ with log link
$g$ and known dispersion.

###  7.2 Power Parameter of Compound Poisson

Given power parameter $p\in(1,2)$, the compound Poisson distribution is an
exponential dispersion model and Tweedie model with point mass at zero
(Jørgensen, 1997b, Ch. 4; Section 2.1). In applications, an unknown power
parameter $p$ can be strongly correlated with an unknown dispersion parameter
$\phi$ for the compound Poisson family (Peters et al., 2009). For application
within the GLM framework, one strategy is to set the value of the power
parameter a priori (Peel et al., 2013). Consider the case where the power
parameter is fixed a priori to a value that is elicited from the expert
together with the prior of Eq. (2d) for the unknown dispersion.

The elicitation begins with the unknown dispersion $\phi=1/\lambda$ for the
class $\mathcal{F}$ of compound Poisson distributions. Conditional on the
choice of systematic component $\hat{\mu}_{0}$ and sample size $w$, the
elicitation for the sample mean $\bar{\psi}$ (Section 4.2) parametrises the
shape $s$ and scale $v_{\phi}(\hat{\mu}_{0},w)$ of the approximating
univariate t distribution by Eqs. (2d), (9) and (10). However, in Eq. (10) the
parameter $r={v_{\phi}(\hat{\mu}_{0},w)ws}/{v(\hat{\mu}_{0})}$ is unknown
because it depends on the power parameter $p$ through the variance function
$v(\hat{\mu}_{0})=\hat{\mu}_{0}^{p}$ of the compound Poisson. Another
elicitation step is therefore required.

Conditional on $\lambda$ and $\hat{\mu}_{0}$, the probability that a single
response equals zero is
$P(\psi=0\mid\hat{\mu}_{0},\lambda,p)=\exp[-\lambda\hat{\mu}_{0}^{2-p}/(2-p)]$
(see Dunn and Smyth, 2008). Given $\hat{\mu}_{0}$, let
$q_{0}=P(\psi=0\mid\hat{\mu}_{0},p)$ denote the probability that the response
is zero by marginalising over the unknown $\lambda$. From Eqs. (2d) and (10)
with $\hat{\mu}_{0}$, $w$, $s$ and $v_{\phi}(\hat{\mu}_{0},w)$ obtained in
Section 4.2, the probability of zero response has the unit gamma density
function (Grassia, 1977; Ratnaparkhl and Mosimann, 1990),

$UG(q_{0}\mid s,r_{p})=\frac{r_{p}^{s/2}q_{0}^{r_{p}-1}(-\log
q_{0})^{s/2-1}}{\Gamma(s/2)}\,,\quad 0<q_{0}<1\,,$ (27)

where $r_{p}=(2-p)v_{\phi}(\hat{\mu}_{0},w)ws/(2\hat{\mu}_{0}^{2})>0$.
Conditional on $\hat{\mu}_{0}$ and the previously elicited $s$ and
$v_{\phi}(\hat{\mu}_{0},w)$, an elicited median $c_{0}$ for the probability of
zero response
$q_{0}=P(\psi=0\mid\hat{\mu}_{0},s,r_{p})=P(\psi=0\mid\hat{\mu}_{0},s,r,p)$
then parametrises $r_{p}$ such that

$p=2\left(1-\frac{r_{p}\hat{\mu}_{0}^{2}}{v_{\phi}(\hat{\mu}_{0},w)ws}\right)\,,\quad
r=\frac{2r_{p}\hat{\mu}_{0}^{2-p}}{2-p}\,,\quad
0<r_{p}<\frac{v_{\phi}(\hat{\mu}_{0},w)ws}{2\hat{\mu}_{0}^{2}}\,.$ (28)

The elicited $r_{p}$ is a solution to $F(c_{0}\mid s,r_{p})=1/2$, where
$F(x\mid s,r_{p})=\int_{0}^{x}UG(z\mid s,r_{p})dz$. This solution is feasible
if the upper bound on $r_{p}$ in Eq. (28) holds so that $1<p<2$, as required
for the compound Poisson. Graphical and numerical feedback for $q_{0}$ based
on the unit gamma distribution and density functions, given the elicited $s$
and $r_{p}$ of Eq. (27), is available to the expert. The elicitation steps for
the power parameter of the compound Poisson distribution are summarised in
Procedure 5. Conditional on the elicited $p\in(1,2)$, elicitation of the
systematic component of a compound Poisson GLM then proceeds as in Section 5.
Elicitation of the power parameter of the compound Poisson supports zero-
inflated discrete data models in Section 7.4.

###  7.3 Extended Generalised Linear Models

Extended generalised linear models (Sweeting, 1981; Jørgensen, 1983, 1997b)
have an observation model that generalises Eq. (1) to the density

$f(\psi\mid\mu,\phi)=a(\psi,\phi)\exp\left[-d(\psi,\mu)/(2\phi)\right]\,,\quad
a(\psi,\phi)>0\,,$

with respect to Lebesgue measure for a suitable function $d(\cdot,\cdot)$.
This is a dispersion model denoted by $\psi\sim DM(\mu,\phi)$. As
$\phi\rightarrow 0$ in a dispersion model then $DM(\mu,\phi)\approx N(\mu,\phi
v(\mu))$ so that $\mu$ is interpreted as a position parameter. Elicitation for
the systematic component given $\lambda$ proceeds according to Section 5. But
the convolution formula of Eq. (5) does not apply to dispersion models in
general, and $\mu$ is not necessarily equal to the mean. The procedure of
Section 4.2 that considers the approximate distribution of the sample mean
$\bar{\psi}$ for elicitation of an unknown dispersion parameter in exponential
dispersion models therefore does not in general apply to extended GLMs.

The elicitation procedure for an unknown dispersion in GLMs (Section 4.2),
however, does apply to special cases of extended GLMs. One example is the
standard simplex distribution on the unit interval (Jørgensen, 1997a, b),
which is a dispersion model denoted by $\psi\sim S(\mu,\lambda)$ with

$\mu=\mathbb{E}[\psi]\,,\quad
d(\psi,\mu)=\frac{(\psi-\mu)^{2}}{\psi(1-\psi)\mu^{2}(1-\mu)^{2}}\,,\quad
v(\mu)=\mu^{3}(1-\mu)^{3}\,,\quad\lambda\in\mathbb{R}^{+}.$

The simplex model is of interest because no exponential dispersion model has
bounded support with index set $\lambda\in\Lambda=\mathbb{R}^{+}$ (Jørgensen,
1997b). Also, for ${\mu}$ small or large, the variance is approximately
$\mathbb{V}[\psi\mid{\mu},\lambda]\approx v({\mu})/\lambda$ (see Section S5.1)
such that Eq. (6) applies. Given $f$ as a simplex model in Eq. (2a), prior
elicitation for the index parameter $\lambda$ proceeds as in Section 4.2 for
the sample mean $\bar{\psi}$ with sample size $w$ conditional on small
$\hat{\mu}_{0}$ or $1-\hat{\mu}_{0}$. This strategy is applied to an
elicitation target embedded within an overdispersed binomial model (Section
7.4) by the case study of Section 8.

###  7.4 Overdispersion and Zero-Inflation for Discrete Data

Overdispersed discrete data are accommodated, given $\mu$ and $\phi$ in Eq.
(2), as follows. The data $y$ conditional on the elicitation target $\psi$ are
assumed Poisson. The unknown dispersion parameter $\phi$ for the target sets
the overdispersion of the count data. Examples of overdispersed Poisson models
include (Winkelmann, 2008): negative binomial with $y\mid\psi\sim
Poisson(\psi)$ and $\psi\sim Gamma(\mu,\phi)$; Poisson-Inverse Gaussian with
$y\mid\psi\sim Poisson(\psi)$ and $\psi\sim InverseGaussian(\mu,\phi)$; and
Poisson-Lognormal with $y\mid\psi\sim Poisson(e^{\psi})$ and $\psi\sim
N(\mu,\phi)$. Elicitation for the target $\psi$ proceeds in the first two
cases according to Section 4 for the dispersion parameter $\phi=1/\lambda$ and
Section 5 for the systematic component. These are examples of Poisson mixture
models that use the Tweedie class as a mixing distribution with power
parameter $p=2$ for the negative binomial and $p=3$ for the Poisson-Inverse
Gaussian distribution (Jørgensen, 1997b, Ch. 4). For these first two cases,
the elicitation target $\psi$ is the conditional mean of a Poisson model such
that $\mathbb{E}[y|\psi]=\psi$, whereas the third case is instead
$\mathbb{E}[y|\psi]=\exp[\psi]$. The Poisson-lognormal is appropriate if the
base $B$ logarithm of the conditional mean
$\psi_{B}=\log_{B}\mathbb{E}[y|\psi]$ is interpretable as an elicitation
target. The procedures of Sections 4 and 5 then apply, see Section S5.2 for
details.

Zero-inflation is another important category for discrete count data
(Winkelmann, 2008). Jørgensen (1997b, Ch. 4) suggests the compound Poisson-
Poisson mixture that features both overdispersion and zero-inflation due to
the zero point mass in the mixing distribution. This model assumes
$y\mid\psi\sim Poisson(\psi)$ with conditional mean
$\mathbb{E}[y\mid\psi]=\psi$ and compound Poisson mixing distribution with
power parameter $1<p<2$ for $\psi$. The dispersion and power parameters for
the target $\psi$ that set the overdispersion and zero-inflation of $y$ are
elicited according to Sections 4.2 and 7.2.

The simplex-binomial mixture (Jørgensen, 1997b, Ch. 5) for overdispersed
counts assumes $y\sim\linebreak Binomial(k,\psi)$, where $y$ is the number of
successes out of $k$ independent trials and $\mathbb{E}[y\mid\psi]=k\psi$,
with $\psi\sim S(\mu,\lambda)$ a standard simplex prior on the success
probability. The elicitation target $\psi$ follows the standard simplex
distribution with dispersion $\phi=1/\lambda$ (Section 7.3). Overdispersion in
$y$ is determined by $\lambda$. The simplex-binomial model with unknown
overdispersion is applied in the case study of Section 8.

##  8 Elicitation Case Study

The seagrass species Zostera nigricaulis provides an important ecological
habitat in Port Philip Bay, which is a large coastal embayment located within
an urbanised region of southeastern Australia. Dissolved inorganic nitrogen
($DIN$, mg per L) and total suspended solids ($TSS$, mg per L) are two factors
thought to influence the productivity of Z. nigricaulis in Port Philip Bay
(Hirst et al., 2017; Nayar et al., 2018). At low $DIN$ levels Z. nigricaulis
may exhibit low growth due to lack of nutrients, whereas increased epiphytes
load at high $DIN$ levels may suppress seagrass growth by shading and
competition for nutrients. Total suspended solids may suppress seagrass by
shading.

In this case study, a regional model for overdispersed binomial counts $y_{i}$
of annual mean seagrass percent cover in Port Philip Bay was elicited. The
model was defined by

$\displaystyle p(y\mid\psi)=\binom{N}{y}\psi^{y}(1-\psi)^{N-y}\,,\quad
y\in\\{0:N\\}\,,$ $\displaystyle\psi\mid\mu,\lambda\sim
S(\mu,\lambda)\,,\quad\log\frac{\mu}{1-\mu}=\eta=X\beta\,,\quad\beta\mid\lambda\sim
N(\delta,\Sigma/\lambda)\,,\quad\lambda\sim G(s/2,r/2)\,,$

with prior structure for $\psi$ according to Eq. (2). The elicitation target
$\psi$ follows the standard simplex distribution (Section 7.3) such that
$f(\psi\mid\mu,1/\lambda)=S(\psi\mid\mu,\lambda)$ in Eq. (2a). The target
$\psi$ is the annual mean probability of intersection with seagrass for each
of $N=49$ points in a 0.5 m2 quadrat. The annual average percent cover of
seagrass is $100\times\psi$. At the observation level, the model is
overdispersed binomial conditional on $\mu$ (Section 7.4). The linear
predictor allows for interactions and quadratic responses with the base 10
logarithm of the annual average $DIN$ in the bottom 1 metre of the water
column, and the annual average $TSS$ in the water column (Section S6). The
elicitation scenarios were defined by $DIN$ and $TSS$ (Table S1) with
corresponding design points that were chosen in consultation with the expert.

Figure 1: Predicted annual mean probability of seagrass ($\psi$) given $TSS$
equal to 0.10 mg/L. The elicited marginal central credible intervals of
probability 1/3 from three scenarios are shown by the three blue vertical
lines. The induced marginal medians and central credible intervals of
probability 0.8 provided as feedback at these three scenarios are given by the
orange points and vertical lines. The predicted median and central credible
intervals of probabilities 1/3 and 0.8 across the range of $DIN$ are given by
the black solid, dashed and dotted lines, respectively.

Probabilistic terms and definitions were reviewed with the expert before the
elicitation session (Section 3). For the random component, the expert was
provided graphical feedback of the approximating Student’s t quantiles and
probability density function according to Section 4.2. Central credible
intervals of probability $1/3$ and $0.90$ were elicited for the sample mean
percent cover. Suggested values of the systematic component $\hat{\mu}_{0}$
emphasised low percent covers between 1% and 10% and sample sizes for $w$
between 10 and 800 quadrats (Sections S6.1). The discrepancies of the
approximation for the final two scenarios considered are plotted in Figure S1,
where the estimated Kolmogorov distances were $0.035$ and $0.054$ for
$\hat{\mu}_{0}$ equivalent to 1% and 10%, respectively. The lower discrepancy
for the lower value of $\hat{\mu}_{0}$ reflects the improved approximation of
the generalised t distribution for small $\hat{\mu}_{0}$ (see Sections 7.3 and
S5.1). The random component accounted for variability introduced by covariates
not associated with $DIN$ or $TSS$, such as independent factors related to
localised disturbance events and the patchy distribution of the target
seagrass species in the study region.

For the systematic component, the elicited central credible interval
probability level was set to $1/3$ such that the probabilities of the
systematic component being above, below or within the interval were equal
(Section S6.2). Graphical and numerical feedback for the elicitation target of
seagrass percent cover was provided for the $1/3$ and $0.8$ central credible
intervals, median, feasible ranges of conditional medians, and plotted density
function. The conditional medians were assessed in terms of relative change
with respect to previously elicited medians, given newly introduced realised
values of the systematic component at successive tree levels. The prior
parametrisation for the unknown $\beta$ induced a concave relationship between
$DIN$ and seagrass percent cover (Figure 1), which is consistent with low
percent cover of seagrass at the extremal levels of $DIN$ and high percent
cover at intermediate levels of $DIN$ given low levels of $TSS$.

Figure 2: Kullback-Leibler divergence between the full dependence model and
all possible truncated canonical vine dependence models. Values above the
horizontal dashed line indicate substantial evidence (sensu Jeffreys, 1961)
against truncation.

The elicitation of dependencies for the systematic component provided
substantial information. The amount of information lost by reporting a
truncated canonical vine dependence model $H^{\prime}$ is assessed in Figure 2
by comparison with the elicited model $H$ that has no truncation. The
Kullback-Leibler divergence between $H$ and $H^{\prime}$ for the joint
distribution of the unknown systematic component and dispersion (see Section
S6.3 for details) is interpreted as an expected log Bayes factor in favour of
$H$ and against $H^{\prime}$ given that $H$ is true. The categories for
evidence suggested by Jeffreys (1961, Appendix B) indicate “substantial”
evidence against a vine dependence model truncated after level
$\mathcal{T}_{0}$ (Figure 2) that omits the elicited canonical vine
altogether. Truncation after levels $\mathcal{T}_{1:3}$ results in values near
the threshold of $\log\sqrt{10}\approx 1.15$ for substantial evidence.

##  9 Conclusion

Expert assessments are made measurable and accessible to scientific
investigation in the Bayesian approach. To support this aspiration, the
statistical elicitation method proposed here provides a unified and flexible
framework for informative prior specification of GLMs and extensions that
apply to a wide array of data types and applications. The elicitation method
is general and accommodates various strategies designed to facilitate
probabilistic expert elicitation.

For example, the development of graphical feedback approaches may usefully
enhance an expert’s understanding of the information contained within a sample
mean via the central limit theorem (e.g., Dinov et al., 2008; Zhang et al.,
2022), and so assist the facilitated elicitation of an unknown dispersion
parameter in Section 4.2. Advances in the communication of statistical
concepts may also assist elaborations here. For instance, different choices of
exponential dispersion models by Eq. (7) induce alternative approximating t
distributions for the sample mean through the corresponding variance
functions. Sample mean predictions elicited across multiple realisations of
the systematic component could then be evaluated for consistency with each of
the models, and so suggest preferences for $f$ in Eq. (2a). The general
elicitation method also allows for the choice of feedback techniques that
present quantiles to the expert. The optimal choice often depends on the
domain, time availability and background of the expert (e.g., Kadane and
Wolfson, 1998; Garthwaite et al., 2005; O’Hagan et al., 2006). The statistical
elicitation method for GLMs and extensions therefore allows for the choice of
quantiles to present as feedback to the expert.

A consistent and recurring theme among recommendations for statistical
elicitation methods is the need for software to support application (Kadane et
al., 1980; Garthwaite et al., 2013; Elfadaly and Garthwaite, 2020; Mikkola et
al., 2021). The procedures that support the statistical elicitation method for
Bayesian GLMs and extensions described here are implemented in the R package
eglm (Hosack, 2023).

## Acknowledgments

The author thanks Keith Hayes, David Peel, two anonymous reviewers and
associate editor for their constructive comments on a previous version of the
manuscript. The author is grateful for the contribution of scientific
expertise provided by Professor Gregory Jenkins, School of BioSciences,
University of Melbourne, in the elicitation case study that was conducted with
ethics approval from the CSIRO Social Science Human Research Ethics Committee
(application number 103/20).

##  Appendix S1 t Distributions

The generalised multivariate t distribution (Kotz and Nadarajah, 2004, Ch. 5)
has the representation $z\sim St_{p}(\delta,\Sigma,r,s)$ if for $p\times 1$
location vector $\delta$ and $p\times p$ covariance matrix $\Sigma$,

$z=\delta+v^{-1/2}y\,,\quad y\sim N(0,\Sigma)\,,\quad v\sim G(s/2,r/2)\,,$

where $\mathbb{E}[v]=s/r$ so that $rv/s\sim G(s/2,s/2)$. An equivalent
parametrisation is

$z=\delta+v^{-1/2}y\,,\quad y\sim N(0,(r/s)\Sigma)\,,\quad v\sim
G(s/2,s/2)\,.$

The latter parametrisation is equivalent to the multivariate t distribution,
see Kotz and Nadarajah (2004, Ch. 1) or Joe (2015, Ch. 2), with density

$t_{p}(z\mid\delta,(r/s)\Sigma,s)=\frac{\Gamma((s+p)/2)}{\Gamma(s/2)(\pi
s)^{p/2}}\left|(r/s)\Sigma\right|^{-\frac{1}{2}}\left[1+\frac{1}{s}(z-\delta)^{\top}(s/r)\Sigma^{-1}\left(z-\delta\right)\right]^{-\frac{s+p}{2}}.$
(S1)

By Eqs. (3) and (S1), the density is equivalently expressed by $z\sim
t_{p}(\delta,(r/s)\Sigma,s)=St_{p}(\delta,(r/s)\Sigma,s,s)=St_{p}(\delta,\Sigma,r,s)$.
The median of $z_{i}$ is $\delta_{i}$, which is also the mean of $z_{i}$ if
$s>1$. As the parameter $s\rightarrow\infty$ then the multivariate t
approaches the multivariate normal with location $m$ and covariance matrix
$\Sigma$ (Kotz and Nadarajah, 2004, Ch. 1). The limiting multivariate normal
for the generalised t is obtained by setting $r=ks$ for $k>0$ then letting
$s\rightarrow\infty$ so that $z\sim
St_{p}(\delta,\Sigma,ks,s)=St_{p}(\delta,k\Sigma,s,s)=t_{p}(\delta,k\Sigma,s)\rightarrow
N(\delta,k\Sigma)$. The univariate t distribution with zero location and scale
of one is equivalent to the Student’s t distribution.

##  Appendix S2 Elicitation of Prior for Unknown Dispersion

###  S2.1 Berry-Esseen Bound and Log Convexity Inequality

Guidance for the choice of sample size $w$ may be obtained by consideration of
the Berry-Esseen bound. If the third absolute central moment of
$\psi\mid\hat{\mu}_{0},\lambda$ is finite, then a bound on the approximation
error of Eq. (6) is available by Eq. (11). The first inequality in Eq. (11)
uses the Berry-Esseen inequality for i.i.d. random variables. This inequality
conditions on $\hat{\mu}_{0}$, $\phi$ and $v(\cdot)$ in Eq. (11). For
application, the facilitator would be required to hypothesise how the expert’s
third absolute central moment $\breve{\mu}^{3}$, given the choice of
$\hat{\mu}_{0}$, might behave if the dispersion $\phi$ is known or fixed to a
best estimate. This task may pose a challenge for intuition and introduce the
need for additional calculations (e.g., Zhang et al., 2022).

Alternatively, the second inequality of Eq. (11) provides an upper bound in
terms of the kurtosis, $\breve{\mu}^{4}$. Analytical expressions of this
measure of concentration in a distribution’s tails, or instead the related
excess kurtosis defined by $\breve{\mu}^{4}-3$, are often presented in
summaries of distributional properties for exponential dispersion models with
fixed $\phi$. For a nonnegative random variable $z$, the log convexity
inequality of Lyapunov (see DasGupta, 2008) is

$\mathbb{E}\left[z^{\alpha_{2}}\right]\mathbb{E}\left[z^{\beta-\alpha_{2}}\right]\leq\mathbb{E}\left[z^{\alpha_{1}}\right]\mathbb{E}\left[z^{\beta-\alpha_{1}}\right]\,,\quad
0\leq\alpha_{1}<\alpha_{2}\leq\beta/2\,.$

For a real-valued random variable $\psi$ with finite variance
$\sigma^{2}=\mathbb{V}[\psi]<\infty$, let $\mu=\mathbb{E}[\psi]$,
$z=(|\psi-{\mu}|)/\sigma$, $\alpha_{1}=2$, $\alpha_{2}=3$ and $\beta=6$ in the
above inequality to obtain

$\left(\breve{\mu}^{3}\right)^{2}=\left(\frac{\mathbb{E}\left[|\psi-\mu|^{3}\right]}{\sigma^{3}}\right)^{2}\leq\frac{\mathbb{E}\left[|\psi-\mu|^{4}\right]}{\sigma^{4}}={\breve{\mu}^{4}}\,.$

The second bound of Eq. (11) follows for $\mu=\hat{\mu}_{0}$ and
$\psi\mid\hat{\mu}_{0},\phi\sim ED(\hat{\mu}_{0},\phi)$ with variance
$\mathbb{V}[\psi\mid\hat{\mu}_{0},\phi]=\phi v(\hat{\mu}_{0})$.

###  S2.2 Approximate Distribution of Sample Mean

Given the ability to draw Monte Carlo samples from the exponential dispersion
model $f=ED(\mu,\phi)$, which is a natural exponential family model
conditional on $\phi$, the following procedure is used to generate $n$ samples
from the joint distribution of the sample mean and an unknown index parameter
$\lambda$. Given the parameters $s$ and $r$ in Eq. (2d) and the choice of
$\hat{\mu}_{0}$ and $w$, for $i=1,\ldots,N$:

  1. 1.

Draw realisation $\lambda_{i}\mid s,r\sim G(s/2,r/2)$ from Eq. (2d).

  2. 2.

Draw realisation $\bar{\psi}_{i}\mid\hat{\mu}_{0},\lambda_{i},w\sim
ED(\hat{\mu}_{0},1/(w\lambda_{i}))$ by Eq. (5).

The above steps generate $N$ paired samples $(\bar{\psi}_{i},\lambda_{i})$
from the joint distribution of $\bar{\psi}$ and $\lambda$ described by $f$ and
$G$.

A Monte Carlo estimate of the cumulative distribution function of the sample
mean, $\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)$, is constructed based on
the empirical cdf of the $N$ samples for the sample mean $\bar{\psi}$. This
estimate for the distribution ${F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)$
includes Monte Carlo error. The Dvoretzky, Kiefer and Wolfowitz inequality
(see Kosorok, 2008, Theorem 11.6) gives

$P\left(\sup_{x\in\mathbb{R}}\sqrt{N}\left|\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)-{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)\right|>\epsilon\right)\leq
2\exp[-2\epsilon^{2}]\,.$

A simultaneous confidence interval of probability $1-\alpha$ is then

$\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)-\epsilon\leq{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)\leq\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)+\epsilon\,,\quad\epsilon=\sqrt{-\log(\alpha/2)/(2N)}\,.$

The estimated cdf $\hat{F}(\bar{\psi}\mid\hat{\mu}_{0},w,s,r)$ and its double-
sided confidence interval is then graphically compared with the approximating
distribution function produced by the elicitation,
$T_{1}([\bar{\psi}-\hat{\mu}_{0}]/\sqrt{v_{\phi}(\hat{\mu}_{0},w)}\mid s)$,
where $T_{1}(\cdot)$ is the distribution function of the Student’s t. The
approximating t distribution corresponds to replacing the exponential
dispersion model $f(\mu,\phi/w)=ED(\mu,\phi/w)$ with the normal distribution
$N(\mu,\phi/w)$. Both choices describe the index parameter $\lambda=1/\phi$
via the gamma distribution $G(\lambda\mid s/2,r/2)$ by Eq. (2d). A numerical
measure of discrepancy is provided by the Kolmogorov distance of Eq. (12),
which may be compared against the numerical error bound $\epsilon$. The latter
is reduced by increasing the Monte Carlo sample size $N$.

###  S2.3 Kullback-Leibler Divergence for Sample Mean

Let $F_{\bar{\psi}}$ and $G_{\lambda}$ denote the distribution functions of
the conditional sample mean ($\bar{\psi}\mid\hat{\mu}_{0},\lambda$) and the
index parameter ($\lambda$), respectively. The Monte Carlo samples from
Section S2.2 allow an estimate of the Kullback-Leibler divergence from Eq.
(13),

$\displaystyle D_{\bar{\psi},\lambda}$
$\displaystyle=\int\int\log\frac{f(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))G(\lambda\mid
s/2,r/2)}{N(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))G(\lambda\mid
s/2,r/2)}dF_{\bar{\psi}}(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))dG_{\lambda}(\lambda\mid
s/2,r/2)$
$\displaystyle=\int\int\log\frac{f(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))}{N(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))}dF_{\bar{\psi}}(\bar{\psi}\mid\hat{\mu}_{0},1/(w\lambda))dG_{\lambda}(\lambda\mid
s/2,r/2)$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\log\frac{f(\bar{\psi}_{i}\mid\hat{\mu}_{0},1/(w\lambda_{i}))}{N(\bar{\psi}_{i}\mid\hat{\mu}_{0},1/(w\lambda_{i}))}=\hat{D}_{\bar{\psi},\lambda}^{N}\,.$

Confidence bounds for the Monte Carlo error (Robert and Casella, 2005, Ch. 3)
are available based on the variance estimate

$\hat{v}_{N}=\frac{1}{N^{2}}\sum_{i=1}^{N}\left[\log\frac{f(\bar{\psi}_{i}\mid\hat{\mu}_{0},1/(w\lambda_{i}))}{N(\bar{\psi}_{i}\mid\hat{\mu}_{0},1/(w\lambda_{i}))}-\hat{D}_{\bar{\psi},\lambda}^{N}\right]^{2}$

by the approximation
$({\hat{D}_{\bar{\psi},\lambda}^{N}-D_{\bar{\psi},\lambda}})/{\hat{v}_{N}^{1/2}}\sim
N(0,1)$ for Monte Carlo sample size $N\rightarrow\infty$. The divergence
$D_{\bar{\psi},\lambda}$ is the discrepancy from the true distribution to the
approximating normal-gamma distribution for the joint distribution of the
sample mean and index parameter $\lambda$. Given equal prior odds for the
choice of the true model versus its approximation, the above directed
divergence is the expected logarithm of the posterior odds in favour of the
true model and against the approximation.

##  Appendix S3 Systematic Component

###  S3.1 Univariate Density and Distribution Functions

Eq. (15) gives the marginal distribution of the systematic component (Section
5.1). The marginal density of the systematic component for scenario
$u_{i}^{\top}$ is

$f\left(\mu_{i}\mid m_{i},V_{i,i},r,s\right)=St_{1}\left(g(\mu_{i})\mid
m_{i},V_{i,i},r,s\right)\,|dg(\mu_{i})/d\mu_{i}|\,.$ (S2)

The univariate marginal distribution function of the systematic component is

$F(z\mid
m_{i},V_{i,i},r,s)=T_{1}\left(\left.\frac{g(z)-m_{i}}{(r/s)^{1/2}V_{i,i}^{1/2}}\,\textrm{sgn}\frac{dg}{d\mu}\right.\mid
s\right)\,,\quad z\in\Omega\,,$ (S3)

where $T_{1}(z\mid s)$ is the cumulative distribution function of the
Student’s t and $\textrm{sgn}\,h$ is the signum function that equals $+1$ if
$h>0$ and $-1$ if $h<0$. The sign of $\textrm{sgn}\,{dg}/{d\mu}$ is constant
within the mean domain because the link function $g$ is continuous and
invertible (Section 2.1). The corresponding quantile function is

$F^{-1}(q\mid
m_{i},V_{i,i},r,s)=g^{-1}\left(m_{i}+(r/s)^{1/2}V_{i,i}^{1/2}T_{1}^{-1}(q\mid
s)\,\textrm{sgn}\,\frac{dg}{d\mu}\right)\,,\quad q\in(0,1)\,.$ (S4)

Eq. (16) in Section 5.2 provides the conditional density of the linear
predictor given realised values of the linear predictor. The univariate
conditional distributions are available and may be used to provide graphical
and numerical feedback to the expert at tree level $\mathcal{T}_{l}$ (Section
5.3). For link function $g$, the univariate conditional density of the
systematic component at scenario $u_{k}^{\top}$ given realised values of the
systematic component $\hat{\mu}_{1:l}$ at the first $l$ scenarios is

$f\left(\mu_{k}\mid m_{k\mid 1:l},V_{k,k\mid
1:l},r,s\right)=St_{1}(g(\mu_{k})\mid m_{k\mid 1:l},V_{k,k\mid
1:l,1:l},r+\zeta(\hat{\eta}_{1:l}),s+l)\,|dg(\mu_{k})/d\mu_{k}|\,,$ (S5)

where $k>l$, $\hat{\eta}_{1:l}=g(\hat{\mu}_{1:l})$ and
$\zeta(\hat{\eta}_{1:l})=(\hat{\eta}_{1:l}-m_{1:l})^{\top}V_{1:l,1:l}^{-1}(\hat{\eta}_{1:l}-m_{1:l})$.
The univariate conditional distribution function of the systematic component
at scenario $u_{k}^{\top}$ conditional on realised values of the systematic
component $\hat{\mu}_{1:l}$ at the first $l$ scenarios is

$\displaystyle F_{\mu}\left(z\mid m_{k\mid 1:l},V_{k,k\mid 1:l},r,s\right)$
$\displaystyle=T_{1}\left(\left.\frac{g(z)-m_{k\mid
1:l}}{\left[(r+\zeta(\hat{\eta}_{1:l}))/(s+l)\right]^{1/2}V_{k,k\mid
1:l}^{1/2}}\,\textrm{sgn}\frac{dg}{d\mu}\right.\mid s+l\right)\,,$ (S6)

where $z\in\Omega$. The conditional quantile function with $q\in(0,1)$ is

$\displaystyle F_{\mu}^{-1}\left(q\mid m_{k\mid 1:l},V_{k,k\mid
1:l},r,s\right)$
$\displaystyle=g^{-1}\left[m_{k|1:l}+\sqrt{\frac{[r+\zeta(\hat{\eta}_{1:l})]V_{k,k\mid
1:l}}{s+l}}T_{1}^{-1}\left(q\mid
s+l\right)\,\textrm{sgn}\,\frac{dg}{d\mu}\right]\,.$ (S7)

###  S3.2 Algorithm S1 for $P\rightarrow R$ Mapping and Bounds for Conditional
Medians

[ht!] Data: $n\times n$ partial correlation matrix
$\mathrm{P}=\left[p_{lk}\right]$ where $p_{lk}=\rho_{l,k|1:(l-1)}$ for $l<k$
Result: $n\times n$ correlation matrix $R=\left[\rho_{l,k}\right]$ initialise
$\rho_{1k}=\rho_{k1}\leftarrow p_{1k}$ for $k=2:n$ ;

2 if  _$n >2$_ then

3 for  _$l=2:(n-1)$_ do

4 for  _$k=(l+1):n$_ do

5 $r\leftarrow p_{lk}$ ;

6 for  _$\gamma=(l-1):1$_ do

7 $r\leftarrow p_{\gamma l}p_{\gamma k}+r(1-p_{\gamma
l}^{2})^{1/2}(1-p_{\gamma k}^{2})^{1/2}$ ;

8

9 end for

10 $\rho_{l,k}=\rho_{k,l}\leftarrow r$ ;

11

12 end for

13

14 end for

15

16 end if

return matrix $R=\left[\rho_{l,k}\right]$ with entries of main diagonal equal
to one. S1(). Mapping from partial correlation matrix $\mathrm{P}$ to
correlation matrix $R$ for canonical vine, see Algorithm 27 of Joe (2015).

Joe (2015) provides the mapping $\mathrm{P}\rightarrow R$ shown in Algorithm
S3.2. At the start of the elicitation for scenario $u_{k}^{\top}$ at tree
level $\mathcal{T}_{l}$ with $1\leq l<k$, let the partial correlation matrix
$\mathrm{P}$ with the zero entry $\mathrm{P}_{l,k}$ replaced by
$\rho\in\\{-1,+1\\}$ be denoted by $\mathrm{P}^{(\rho)}$. Let $R^{(\rho)}$
denote the corresponding correlation matrix induced by Algorithm S3.2 applied
to $\mathrm{P}^{(\rho)}$. The corresponding scale matrix for the linear
predictor is
$V^{(\rho)}=\operatorname{diag}[V]^{1/2}R^{(\rho)}\operatorname{diag}[V]^{1/2}$,
where the diagonal entries of the scale matrix were previously elicited in
Section 5.1 such that
$\operatorname{diag}[V^{(\rho)}]=\operatorname{diag}[V]$. From Eq. (20), the
feasible bounds for $c_{k|1:l}$ are given by the range of

$g^{-1}\left(m_{k}+V^{(\rho)}_{k,1:l}\left[V^{(\rho)}_{1:l,1:l}\right]^{-1}\left(\hat{\eta}_{1:l}-m_{1:l}\right)\right)\textrm{
for }\rho\in\\{\pm 1\\}\,.$

These bounds on the conditional median are available for graphical or
numerical presentation to the expert.

##  Appendix S4 Induced Prior

Section 6 presents the elicited parameters $\theta=\left\\{s,r,m,V\right\\}$
for the saturated model $H$ described by Eqs. (2) and (4). The choice of
observation model $f^{\prime}\in\mathcal{F}$ in Eq. (2a) and model matrix
$X\in\mathcal{X}$ of full column rank in Eq. (2c) induces a new alternative
model $H^{\prime}$ with parameters
$\theta^{\prime}=\\{s^{\prime},r^{\prime},\delta^{\prime},\Sigma^{\prime}\\}$.
The conditional property of information (Kullback, 1959, Ch. 2) applied to Eq.
(25) obtains

$D(H:H^{\prime})=\int
D\left(H:H^{\prime}\mid\lambda\right)G\left(\lambda\mid\frac{s}{2},\frac{s}{2}\right)d\lambda+\int\log\frac{G\left(\lambda\mid\frac{s}{2},\frac{s}{2}\right)}{G\left(\lambda\mid\frac{s^{\prime}}{2},\frac{s^{\prime}}{2}\right)}G\left(\lambda\mid\frac{s}{2},\frac{s}{2}\right)d\lambda\,$
(S8)

where

$\displaystyle D\left(H:H^{\prime}\mid\lambda\right)$
$\displaystyle=\int\log\frac{N\left(\eta\mid
Am,\frac{r}{s\lambda}AVA^{\top}\right)}{N\left(\eta\mid
X\delta^{\prime},\frac{r^{\prime}}{s^{\prime}\lambda}X\Sigma^{\prime}X^{\top}\right)}N\left(\eta\mid
Am,\frac{r}{s\lambda}AVA^{\top}\right)d\eta$
$\displaystyle=\frac{1}{2}\log\frac{\left|S^{\prime}\right|}{\left|S\right|}+\frac{1}{2}\mathop{\mathrm{Tr}}\left[S\left(S^{\prime-1}-S\right)\right]+\frac{\lambda}{2}\mathop{\mathrm{Tr}}\left[S^{\prime-1}(\gamma-\gamma^{\prime})(\gamma-\gamma^{\prime})^{\top}\right]\,,$
(S9)

is the Kullback-Leibler divergence of a multivariate normal (Kullback, 1959,
Ch. 9) with

$\gamma=Am\,,\quad\gamma^{\prime}=X\delta^{\prime}\,,\quad
S=(r/s)AVA^{\top}\,,\quad
S^{\prime}=(r^{\prime}/s^{\prime})X\Sigma^{\prime}X^{\top}\,.$

With respect to $\gamma$ and $\gamma^{\prime}$, Eq. (S9) is minimised for the
substitutions $A=\hat{A}$ and $\delta^{\prime}=\hat{\delta}$ from Eq. (26),
whereby $\gamma=\gamma^{\prime}$. Then the last addend in Eq. (S9) is zero and
$D(H:H^{\prime}\mid\lambda)$ does not depend on $\lambda$. Thus, with respect
to $s^{\prime}$, $D(H:H^{\prime})$ is minimised by the substitution of
$s^{\prime}={s}$ so that the second integral of Eq. (S8) is zero. Next, note
that the substitutions $\Sigma^{\prime}=\hat{\Sigma}$, $s^{\prime}=\hat{s}$
and $r^{\prime}=\hat{r}$ into $S^{\prime}$ and $A=\hat{A}$ into $S$ obtain the
result $S=S^{\prime}=(r/s)X\left(X^{\top}V^{-1}X\right)^{-1}X^{\top}$. Eqs.
(S8) and (S9) are therefore zero given the substitutions from Eq. (26). The
result $D(H:H^{\prime})=0$ is the minimum by the convexity property of the
Kullback-Leibler divergence (Kullback, 1959, Ch. 2).

If $p<n$ and, contrary to Eq. (26), the identity map $A=I$ is assumed then
$D(H:H^{\prime})=\infty$ because $X\beta$ in model $H^{\prime}$ lies in a
subspace of the saturated model $H$. The model $H^{\prime}$ would then
contradict the expert assessments documented by $H$. If instead the assumption
that $A$ is a projection matrix onto $\mathbf{C}(X)$ is accepted, then the
optimal projection is $\hat{A}$ and the divergence is minimised such that
$D\left(H:H^{\prime}\right)=0$.

If $p=n$ then the optimal projection is the identity $\hat{A}=I$ with
$\hat{\delta}=X^{-1}m$ and $\hat{\Sigma}=X^{-1}VX^{-\top}$ by Eq. (26).
Specifying a square model matrix $X$ of full rank ensures a non-singular
transformation between $\beta$ and the elicited systematic component $\mu$
with no loss of information.

For the case of conditional means priors (CMPs) where $p<n$, Bedrick et al.
(1996) discuss how an expert is unlikely to provide responses within
$\textbf{C}(X)$. For a CMP with known dispersion in the normal linear model,
Bedrick et al. (1996) derive $\hat{\delta}$ and $\hat{\Sigma}$ in Eq. (26) by
first imagining $n-p$ predictors with associated coefficients $\beta^{\ast}$,
then taking the conditional distribution of $\beta\mid\beta^{\ast}=0$ as the
prior. The above interpretation in terms of $D\left(H:H^{\prime}\right)$ does
not depend on imagined coefficients $\beta^{\ast}$ and applies to the more
general Eq. (2).

##  Appendix S5 Extensions

###  S5.1 Variance Approximation for Standard Simplex Distribution

The standard simplex distribution with density

$S(y\mid\mu,\lambda)=\left(\frac{\lambda}{2\pi\left[y(1-y)\right]^{3}}\right)^{\frac{1}{2}}\exp\left[-\frac{\lambda}{2}d(y,\mu)\right]\,,\quad
y,\mu\in(0,1)\,,$

where $d(y,\mu)=\frac{(y-\mu)^{2}}{y(1-y)\mu^{2}(1-\mu)^{2}}$, is an example
of a dispersion model available for use by an extended GLM (Section 7.3). A
related distribution (Jørgensen, 1997a, b, Ch. 5) is

$S^{+}(y\mid\mu,\lambda)=\frac{1}{\Gamma(1/2,\lambda/2)\sqrt{\pi
e^{\lambda}}}\left[y(1-y)\right]^{-\frac{1}{2}}\exp\left[-\frac{\lambda}{2}d(y,\mu)\right]\,,\quad
y,\mu\in(0,1)\,,$ (S10)

where $d(y,\mu)=(y-\mu)^{2}/(y(1-y))$ and
$\Gamma(a,b)=\int_{b}^{\infty}t^{a-1}e^{-t}dt$ is the (upper) incomplete gamma
function. Jørgensen (1997a, b, Ch. 5) suggest calculation of the variance for
$y\sim S(\mu,\lambda)$ by use of the mixed moments
$\mathbb{E}[y^{b_{1}}(1-y)^{b_{2}}]$ (see Eq. 7.2 of Jørgensen (1997a) or Eq.
5.57 of Jørgensen (1997b)), which can be expressed as a ratio of normalising
constants of the standard simplex distribution and Eq. (S10) (see Table 1 of
Jørgensen (1997a) or Table 5.1 of Jørgensen (1997b)). Applying this
derivation, Jørgensen (1997a) notes that the variance of the standard simplex
is approximately $\mathbb{V}[y]=\mu^{3}(1-\mu)^{3}/\lambda$ for large
$\lambda$ or small $\mu$ or $1-\mu$.

This approximation to the variance of the standard simplex distribution can be
shown as follows. Since $\mathbb{E}[y]=\mu$ for the standard simplex
distribution, applying the approach described above for the variance of the
standard simplex distribution using mixed moments with $b_{1}=b_{2}=1$ obtains
the expression (cf. Jørgensen, 1997a, b, Ch. 5),

$\mathbb{V}[y]=\mu(1-\mu)-\left(\frac{\lambda}{2}\right)^{1/2}\exp\left[\frac{\lambda}{2\mu^{2}(1-\mu)^{2}}\right]\Gamma\left(\frac{1}{2},\frac{\lambda}{2\mu^{2}(1-\mu)^{2}}\right)\,.$
(S11)

The equivalence $\Gamma(1/2,z)=\sqrt{\pi}\textrm{erfc}(\sqrt{z})$ (DLMF, 2020,
Eq. 8.4.6), where $\textrm{erfc}(z)$ is the complementary error function,
suggests for large $\lambda$ or small $\mu$ or $1-\mu$ the approximation

$\displaystyle\Gamma\left(\frac{1}{2},\frac{\lambda}{2\mu^{2}(1-\mu)^{2}}\right)$
$\displaystyle=\sqrt{\pi}\textrm{erfc}\left(\sqrt{\frac{\lambda}{2\mu^{2}(1-\mu)^{2}}}\right)$
$\displaystyle\approx\left(\frac{2}{\lambda}\right)^{1/2}\exp\left[-\frac{\lambda}{2\mu^{2}(1-\mu)^{2}}\right]\left\\{\mu(1-\mu)-\left(\frac{\mu^{3}(1-\mu)^{3}}{\lambda}\right)\right\\}\,.$

The above approximation uses the first two terms of an asymptotic expansion of
the complementary error function $\textrm{erfc}(z)$ that holds as
$z\rightarrow\infty$ (DLMF, 2020, Eq. 7.12.1). Substituting the above
approximation into Eq. (S11) yields the approximation $\mathbb{V}[y]\approx
v(\mu)/\lambda$, where $v(\mu)=\mu^{3}(1-\mu)^{3}$ is the variance function of
the standard simplex distribution. This approximation holds for large
$\lambda$ or small $\mu$ or $1-\mu$.

###  S5.2 Lognormal with Unknown Scale

In Section 7.4, the mean of the Poisson distributed random variable $y$ given
$\psi$ is $z=e^{\psi}$. Since $\psi\sim N(\mu,\phi)$, the random variable $z$
follows a lognormal distribution, which is not an exponential dispersion
model. But $k\log z\sim N(k\mu,k^{2}\phi)$, given $\mu$ and $\phi$ in Eq. (2),
where the known constant $k=1/\log B$ accounts for the desired base $B$ of the
logarithm. A lognormal model is appropriate if the logarithm of the
elicitation target with base $B$ is interpretable, denoted by
$\psi_{B}=\log_{B}z=k\log z$. If the elicitation target $\psi_{B}$ is
described by Eq. (2) for $f$ normal and $g$ the identity link, then the
marginal distribution is $\psi_{B}\sim
St_{n}(X\delta_{B},X\Sigma_{B}X^{\top},r,s)$. Procedure 1 applied to the
transformed target $\psi_{B}$ then elicits $\delta_{B}$, $\Sigma_{B}$, $s$ and
$r$. The induced location and conditional scale for $\psi$ in Eq. (2c) are
then $\delta=\delta_{B}/k$ and $\Sigma=\Sigma_{B}/k^{2}$, and $\lambda\sim
G(s/2,r/2)$ in Eq. (2d).

##  Appendix S6 Elicitation Case Study

In the elicitation case study of Section 8, the potential response of annual
mean Zostera nigricaulis percent cover to changes in nitrogen levels and
suspended sediments is examined for a seagrass ecosystem within Port Phillip
Bay, Australia. The linear predictor

$\displaystyle\eta_{i}=$
$\displaystyle\beta_{0}+\beta_{1}\log_{10}DIN_{i}+\beta_{2}TSS_{i}+\beta_{3}\log_{10}DIN_{i}\times
TSS_{i}+\beta_{4}(\log_{10}DIN_{i})^{2}$
$\displaystyle+\beta_{5}TSS_{i}^{2}+\beta_{6}(\log_{10}DIN_{i})^{2}\times
TSS_{i}^{2}$ (S12)

allowed for interactions and quadratic responses with the base 10 logarithm of
the annual average $DIN$ in the bottom 1 metre of the water column, and the
annual average $TSS$ in the water column. A $n\times n$ model matrix $X$ of
full column rank with the polynomial basis expansion for Eq. (S12) was derived
from the $n=7$ elicitation scenarios in the elicitation scenario matrix $U$,
see Table S1. The $i$th row of $X$ is a design point $x_{i}^{\top}$ that
corresponds to the $i$th scenario in $U$.

Table S1: The elicitation scenarios of $U$ defined by the annual averages of
dissolved inorganic nitrogen in the bottom 1 metre of the water column ($DIN$
in mg per litre) and total suspended solids in the entire water column ($TSS$
in mg per litre).
Scenario | $DIN$ | $TSS$
---|---|---
1 | 0.0001 | 0.1000
2 | 0.0500 | 0.1000
3 | 0.5000 | 0.1000
4 | 0.0001 | 12.2500
5 | 0.0001 | 50.0000
6 | 0.0500 | 50.0000
7 | 0.5000 | 50.0000

The elicitation session occurred by video conference in November 2020. At the
start of the elicitation session, the contributing expert with scientific
expertise of seagrass ecosystems in coastal embayments of southeastern
Australia was educated in subjective probability (sensu Lindley, 2014, Ch. 3),
cognitive biases (e.g., Tversky and Kahneman, 1974) and engaged in practice
elicitation exercises (O’Hagan et al., 2006). The definition for the quantity
of interest of the elicitation (Section 8) was reviewed along with the choice
of the covariates that together defined the elicitation scenarios in $U$
(Table S1).

###  S6.1 Random Component

The elicitation of the unknown dispersion parameter (Section 4) determines the
amount of overdispersion in the binomial response for the seagrass study.
Central credible intervals of probabilities $1/3$ and $0.90$ were elicited for
the sample mean $\bar{\psi}$ conditioned on values of the systematic component
$\hat{\mu}_{0}$ emphasising small values from $0.01$ to $0.10$ with possible
sample sizes explored for $w$ between 10 and 800 quadrats. The elicited values
of $s$ and $r$ were compared across the different choices of $\hat{\mu}_{0}$
and $w$ and iteratively modified. The final elicited parameters for the gamma
prior of the index parameter $\lambda$ were $s=14.3$ and $r=118$. The final
scenario evaluated the sample mean prediction for a percent cover of 1%
(corresponding to $\hat{\mu}_{0}=0.01$) with a relatively small sample size
$w=10$. The elicited random component reflected observation uncertainty
introduced by localised disturbance events and the patchy distribution of the
seagrass.

The discrepancies between the estimated distribution functions of the sample
mean against the approximate generalised t approximation are shown for the
final two choices of $w$ and $\hat{\mu}_{0}$. In each case, Monte Carlo
realisations for $\psi$ were drawn by the inverse cumulative distribution
function method for the exact Studentised simplex density (Jørgensen, 1997a,
b, Ch. 5) derived from Eq. (2d) and the simplex density (Section 7.3),

$\displaystyle f(\psi\mid\mu,s,r)=$ $\displaystyle\int
S(\psi\mid\mu,\lambda)G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)d\lambda$
$\displaystyle=$
$\displaystyle\frac{\sqrt{a}}{\beta\left(\frac{b-a}{2a},\frac{1}{2}\right)}\left[\psi(1-\psi)\right]^{-3/2}\left[1+a\frac{(\psi-\mu)^{2}}{\psi(1-\psi)\mu^{2}(1-\mu)^{2}}\right]^{-\frac{b}{2a}}\,,$

where $a=1/r$ and $b=(s+1)/r$. The Monte Carlo estimate and its 0.95
confidence interval were calculated according to Section S2.2 for the
penultimate scenario of $\hat{\mu}_{0}=0.10$ and the final scenario of
$\hat{\mu}_{0}=0.01$ for $w=10$ (Figure S1). The estimated Kolmogorov
distances were $0.054$ and $0.035$, respectively, which reflects the better
approximation of the generalised t for the lower value of $\hat{\mu}_{0}$ in
the simplex distribution (Section 7.3).

(a) Penultimate scenario: $\hat{\mu}_{0}=0.10$.

(b) Final scenario: $\hat{\mu}_{0}=0.01$.

Figure S1: Comparisons of estimated cumulative distribution function for the
sample mean against the generalised t approximation for the final two
scenarios considered for the unknown dispersion with $w=10$. The estimated
distribution function and its 0.95 confidence interval are based on $N=2000$
Monte Carlo realisations for the sample mean $\bar{\psi}$.

###  S6.2 Systematic Component

For marginal elicitation of the systematic component (Section 5.1), the
numerical feedback included the scenario information encoded in the scenario
matrix $U$ (Table S1), which was defined by levels of dissolved inorganic
nitrogen and total suspended solids. The marginal centres and scales of the
transformed generalised t distribution were elicited according to Section 5.1.
Graphical and numerical feedback for the elicitation target of seagrass
percent cover were provided for the $1/3$ and $0.80$ central credible
intervals, along with the median and plotted probability density function.

These same quantities were also presented for the elicited conditional
distributions at tree levels $\mathcal{T}_{1:(n-1)}$ (Section 5.3).
Conditional medians were elicited where, in addition to the information in
$U$, hypothetical realised values for the systematic component were introduced
at tree levels greater than zero (Sections 5.2 and 5.3). The realised values
for the systematic component $\hat{\mu}_{1:l}$ were sequentially introduced
with increasing tree level $\mathcal{T}_{l}$, $l=1,\ldots,n-1$. Thus, a new
realised $\hat{\mu}_{l}$ was introduced for the $l$th scenario at tree level
$\mathcal{T}_{l}$. The conditioning values $\hat{\mu}_{1:l}$ were generated as
in Section 5.4, where their tail probabilities of $(1-\alpha)/2$ were derived
from Eq. (4) given known dispersion $\lambda=1$. This latter choice induced
smaller relative changes in the conditioning values relative to their
conditional medians than would have resulted from the alternative approach
presented in Section 5.4 because of the substantial overdispersion elicited in
this case study. Given this choice, the conditioning values were determined by
the elicited values of $m$ and $V$ in Eq. (2c). The magnitude of the random
component was in this sense screened out while selecting conditioning values
for the systematic component.

The conditioning values $\hat{\mu}_{1:l}$ were interpreted in terms of
relative change with respect to previously elicited medians. During the
elicitation session, the realised values of the systematic component were
interpreted as proportional changes with respect to the predicted median
responses at lower tree levels. The expert judged that negligible information
was provided at tree level $\mathcal{T}_{3}$ for scenarios 4 to 7, and the
conditional medians were unadjusted at that level. The rationale, to
paraphrase, was that scenarios 1 and 3 were predicted to have similar levels
of seagrass percent cover (see Figure 1), though for different reasons: Low
$DIN$ in scenario 1 suggested possible nutrient constraints on seagrass
growth, whereas high $DIN$ in scenario 3 suggested enhanced growth of
epiphytes and possible light suppression on seagrass. The additional
realisation of scenario 3 at tree level $\mathcal{T}_{3}$ did not add much
additional information. Scenario 4, however, had a higher amount of $TSS$ than
scenarios 1 and 3 (Table S1), and the conditional medians of scenarios 5 to 7
were updated at the next tree level $\mathcal{T}_{4}$.

###  S6.3 Kullback-Leibler Divergence for Truncated Vine Models

The directed divergences from the full dependence model to each of the
truncated dependence models are evaluated as follows. Define the truncated
model $H(t)$ with scale $V(t)$ given truncated dependence after level
$\mathcal{T}_{t}$ (Section 5.5) and corresponding alternative fitted models
$H^{\prime}(t)$ for $t\in\\{0,\ldots,n-1\\}$. From Eq. (26) and the
nonsingular $X$, the divergence $D(H(t):H^{\prime}(t))$ is minimised for the
choice of parameters

$\begin{gathered}\hat{\delta}(t)=X^{-1}m\,,\quad\quad\hat{\Sigma}(t)=\left(X^{\top}V(t)^{-1}X\right)^{-1}\,,\\\
\hat{A}(t)=I\,,\quad\hat{s}=s\,,\quad\hat{r}=r\,.\end{gathered}$ (S13)

Let $\hat{H}(t)$ denote the truncated model with parameters given by Eq. (S13)
and model matrix $X$. For the joint distribution of the systematic component
and dispersion, the Kullback-Leibler divergence between the full dependence
model $H$ described by Eqs. (2) and (4) the alternative models $\hat{H}(t)$
truncated at level $t\in\\{0,\ldots,n-1\\}$ is

$\displaystyle D\left(H:\hat{H}(t)\right)=$ $\displaystyle\int
D\left(H:\hat{H}(t)\mid\lambda\right)G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)d\lambda+\int\log\frac{G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)}{G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)}G\left(\lambda\mid\frac{s}{2},\frac{r}{2}\right)d\lambda$
(S14) $\displaystyle=$
$\displaystyle\frac{1}{2}\log\frac{\left|S^{\prime}(t)\right|}{\left|V\right|}+\frac{1}{2}\mathop{\mathrm{Tr}}\left[V\left(S^{\prime}(t)^{-1}-V^{-1}\right)\right]$
$\displaystyle+\frac{s}{2r}\mathop{\mathrm{Tr}}\left[S^{\prime}(t)^{-1}(m-\gamma^{\prime}(t))(m-\gamma^{\prime}(t))^{\top}\right]\,$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log\frac{\left|V(t)\right|}{\left|V\right|}+\frac{1}{2}\mathop{\mathrm{Tr}}\left[V\left(V(t)^{-1}-V^{-1}\right)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log\frac{\left|R(t)\right|}{\left|R\right|}+\frac{1}{2}\mathop{\mathrm{Tr}}\left[RR(t)^{-1}\right]-\frac{n}{2}\,,$

where $\gamma^{\prime}(t)=X\hat{\delta}(t)=m$ and
$S^{\prime}(t)=X\hat{\Sigma}(t)X^{\top}=V(t)$ is the scale matrix with
truncated dependence and corresponding correlation matrix $R(t)$. The above
uses the conditional property of information as in Eq. (S8). Results are shown
in Figure 2.

##  Appendix S7 Elicitation Procedures

Data: Set of exponential dispersion models $\mathcal{F}$.

Result: $s$ and $r$ in Eq. (2d); $v_{\phi}(\hat{\mu}_{0},w)$ in Eq. (10) if
dispersion is unknown.

1 Choose $f\in\mathcal{F}$ (Section 3.1);

2 if  _$\phi$ is known_ then

3 Let $s,r\rightarrow\infty$ such that $\lambda\triangleq 1/{\phi}=s/r$ for
known dispersion ${\phi}$ (Section 4.1);

4 return $s$, $r$;

5

6 else

7 Choose $\hat{\mu}_{0}$ and sample size $w$ for sample mean
$\bar{\psi}\mid\hat{\mu}_{0},w$ (Section 4.2);

8 repeat

9 if  _Candidate $\hat{\mu}_{0}$ unsatisfactory to expert_ then

10 Adjust $\hat{\mu}_{0}$;

11

12 if  _Approximation to sample mean distribution unsatisfactory_ then

13 Increase $w$;

14

15 Elicit lower bound $d_{1}$ of central credible interval with probability
$\alpha_{1}$ for $\bar{\psi}\mid\hat{\mu}_{0},w$;

16 Elicit lower bound $d_{2}<d_{1}$ of central credible interval with
probability $\alpha_{2}>\alpha_{1}$ for $\bar{\psi}\mid\hat{\mu}_{0},w$;

17 Parametrise $s$ by Eq. (9) and $v_{\phi}(\hat{\mu}_{0},w)$ with induced $r$
by Eq. (10);

18 Provide graphical and numerical feedback of approximate distribution via
Eq. (7);

19 Assess approximation accuracy (Section 4.2);

20

21 until  _Approximation to sample mean is accepted_ ;

22 return $s$, $r$, $v_{\phi}(\hat{\mu}_{0},w)$;

23

24 end if

Procedure 2 Elicitation of dispersion parameter for GLMs; see Section 7 for
alternative strategies applicable to extensions of GLMs.

Data: Elicitation scenario matrix $U$, link function $g$; $s$ and $r$ from
Procedure 2.

Result: Location vector $m$ and main diagonal of scale matrix
$\operatorname{diag}[V]$ for the linear predictor $\eta$.

1 for  _$i=1:n$_ do

2 repeat

3 Elicit central credible interval $(a_{i},b_{i})$ of probability $\alpha$ for
$\mu_{i}\mid u_{i}^{\top}$ (Section 5.1);

4 Parametrise $m_{i}$ and $V_{ii}$ conditional on $g$, $s$, $r$,
$(a_{i},b_{i})$ (Section 5.1);

5 Provide graphical and numerical feedback via Eqs. (S2), (S3) and (S4) ;

6

7 until  _expert accepts $f\left(\mu_{i}\mid m_{i},V_{i,i},r,s\right)$_;

8 end for

return $m$, $\operatorname{diag}[V]$.

Procedure 3 Elicitation of location vector $m$ and scale parameters
$\operatorname{diag}[V]$ for linear predictor $\eta$.

Data: Elicitation scenario matrix $U$ and link function $g$; $s$ and $r$ from
Procedure 2; location vector $m$ and diagonal scale $\operatorname{diag}[V]$
from Procedure 3; choice of uppermost tree level $\mathcal{T}_{t}$ for
elicitation with $t\in\left\\{0,\ldots,n-1\right\\}$.

Result: $n\times n$ positive definite correlation matrix $R$

1 Initialise $\mathrm{P}\leftarrow I$ with $I$ the $n\times n$ identity matrix
;

2 if  _$t >0$_ then

3 Select conditioning value $\hat{\mu}_{1}\in\left\\{a_{1},b_{1}\right\\}$
(Section 5.4) ;

4 $\hat{\eta}_{1}\leftarrow g(\hat{\mu}_{1})$;

5 for  _$k=2:n$_ do

6 repeat

7 Elicit feasible $c_{k|1}$ for ${\mu}_{k}\mid u_{k}^{\top},\hat{\mu}_{1}$ and
set $m_{k|1}\leftarrow g(c_{k|1})$ (Section 5.3) ;

8 Parametrise $V_{k,1}$ and $V_{k,k|1}$ from Eqs. (21) and (22) ;

9 Provide graphical and numerical feedback via Eqs. (S5), (S6) and (S7) ;

10

11 until  _expert accepts $f\left(\mu_{k}\mid m_{k\mid 1},V_{k,k\mid
1},r,s\right)$_;

12 $\mathrm{P}_{1k}=\rho_{1,k}\leftarrow V_{k,1}/\sqrt{V_{1,1}V_{k,k}}$;

13

14 end for

15 if  _$t >1$_ then

16 for  _$l=2:t$_ do

17 evaluate $(a_{l|1:(l-1)},b_{l|1:(l-1)})$ from $m_{l|1:(l-1)}$,
$V_{l,l|1:(l-1)}$ (Section 5.4) ;

18 select conditioning value
$\hat{\mu}_{l}\in\left\\{a_{l|1:(l-1)},b_{l|1:(l-1)}\right\\}$ ;

19 $\hat{\eta}_{l}\leftarrow g\left(\hat{\mu_{l}}\right)$;

20 for  _$k=(l+1):n$_ do

21 repeat

22 Elicit feasible $c_{k|1:l}$ for ${\mu}_{k}\mid
u_{k}^{\top},\hat{\mu}_{1:l}$ and set $m_{k|1:l}\leftarrow g(c_{k|1:l})$
(Section 5.3) ;

23 Parametrise $V_{k,l}$ and $V_{k,k|1:l}$ from Eqs. (21) and (22) ;

24 Provide graphical and numerical feedback via Eqs. (S5), (S6) and (S7) ;

25

26 until  _expert accepts $f\left(\mu_{k}\mid m_{k\mid 1:l},V_{k,k\mid
1:l},r,s\right)$_;

27 Evaluate $\mathrm{P}_{lk}=\rho_{l,k|1:(l-1)}$ by Eq. (19) with substitution
$\gamma=l$ ;

28

29 end for

30

31 end for

32

33 end if

34

35 end if

36$R\leftarrow P$ by Algorithm S3.2 ;

return $R$.

Procedure 4 Elicitation of positive definite correlation matrix $R$.

Data: Set of compound Poisson distributions $\mathcal{F}$ with $p\in(1,2)$ and
unknown dispersion; $s$ and $v_{\phi}(\hat{\mu}_{0},w)$ from Procedure 2;

Result: $r$ in Eq. (2d), $p\in(1,2)$;

1 Elicit $s$ and $v_{\phi}(\hat{\mu}_{0},w)$ by Section 4.2 and Procedure 2
for sample mean $\bar{\psi}\mid\hat{\mu}_{0},w$ ;

2 if  _unknown $p\in(1,2)$_ then

3 repeat

4 Elicit median for $q_{0}=P(\psi=0\mid\hat{\mu}_{0})$ to parametrise $r_{p}$
(Section 7.2) ;

5 Parametrise $r$ and $p$ by Eq. (28) ;

6 Provide graphical and numerical feedback ;

7

8 until  _expert accepts $q_{0}\sim UG(s,r_{p})$, see Eq. (27)_;

9 else

10 Given known $p\in(1,2)$, set $r\leftarrow
v_{\phi}(\hat{\mu}_{0},w)ws/\hat{\mu}_{0}^{p}$ ;

11

12 end if

return $s$, $r$, $p$.

Procedure 5 Elicitation of power parameter $p$ for compound Poisson model.

## References

  * Aitchison (1982) Aitchison, J. (1982).  The statistical analysis of compositional data.  Journal of the Royal Statistical Society: Series B (Methodological) 44(2), 139–160.
  * Aitchison and Shen (1980) Aitchison, J. and S. M. Shen (1980).  Logistic-normal distributions: Some properties and uses.  Biometrika 67, 261–272.
  * Al-Awadhi and Garthwaite (1998) Al-Awadhi, S. A. and P. H. Garthwaite (1998).  An elicitation method for multivariate normal distributions.  Communications in Statistics – Theory and Methods 27, 1123–1142.
  * Al-Awadhi and Garthwaite (2001) Al-Awadhi, S. A. and P. H. Garthwaite (2001).  Prior distribution assessment for a multivariate normal distribution: An experimental study.  Journal of Applied Statistics 28, 5–23.
  * Bedford and Cooke (2001) Bedford, T. and R. M. Cooke (2001).  Probability density decomposition for conditionally dependent random variables modeled by vines.  Annals of Mathematics and Artificial Intelligence 32, 245–268.
  * Bedford and Cooke (2002) Bedford, T. and R. M. Cooke (2002).  A new graphical model for dependent random variables.  The Annals of Statistics 30, 1031–1068.
  * Bedrick et al. (1996) Bedrick, E. J., R. Christensen, and W. Johnson (1996).  A new perspective on priors for generalized linear models.  Journal of the American Statistical Association 91, 1450–1460.
  * Bové and Held (2011) Bové, D. S. and L. Held (2011).  Hyper-_g_ priors for Generalized Linear Models.  Bayesian Analysis 6, 387–410.
  * Chen and Ibrahim (2003) Chen, M.-H. and J. G. Ibrahim (2003).  Conjugate priors for generalized linear models.  Statistica Sinica 13, 461–476.
  * Cooney and White (2023) Cooney, P. and A. White (2023).  Incorporating expert opinion on observable quantities into statistical models - a general framework.  arXiv:2302.06391v1.
  * DasGupta (2008) DasGupta, A. (2008).  Asymptotic Theory of Statistics and Probability.  New York: Springer.
  * Denham and Mengersen (2007) Denham, R. and K. Mengersen (2007).  Geographically assisted elicitation of expert opinion for regression models.  Bayesian Analysis 2, 99–136.
  * Dinov et al. (2008) Dinov, I. D., N. Christou, and J. Sanchez (2008).  Central limit theorem: New SOCR applet and demonstration activity.  Journal of Statistics Education 16(2).
  * DLMF (2020) DLMF (2020).  NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.1.0 of 2020-12-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Dunn and Smyth (2008) Dunn, P. K. and G. K. Smyth (2008).  Evaluation of Tweedie exponential dispersion models using Fourier inversion.  Statistics and Computing 18, 73–86.
  * Elfadaly and Garthwaite (2015) Elfadaly, F. G. and P. H. Garthwaite (2015).  Eliciting prior distributions for extra parameters in some generalized linear models.  Statistical Modelling 15, 345–365.
  * Elfadaly and Garthwaite (2017) Elfadaly, F. G. and P. H. Garthwaite (2017).  Eliciting Dirichlet and Gaussian copula prior distributions for multinomial models.  Statistics and Computing 27, 449–467.
  * Elfadaly and Garthwaite (2020) Elfadaly, F. G. and P. H. Garthwaite (2020).  On quantifying expert opinion about multinomial models that contain covariates.  Journal of the Royal Statistical Society. Series A 183, 959–981.
  * Fahrmeir and Tutz (1994) Fahrmeir, L. and G. Tutz (1994).  Multivariate Statistical Modelling Based on Generalized Linear Models.  Springer-Verlag New York Inc.
  * Garthwaite et al. (2013) Garthwaite, P. H., S. A. Al-Awadhi, F. G. Elfadaly, and D. J. Jenkinson (2013).  Prior distribution elicitation for generalized linear and piecewise-linear models.  Journal of Applied Statistics 40, 59–75.
  * Garthwaite and Dickey (1988) Garthwaite, P. H. and J. M. Dickey (1988).  Quantifying expert opinion in linear regression problems.  Journal of the Royal Statistical Society. Series B (Methodological) 50, 462–474.
  * Garthwaite et al. (2005) Garthwaite, P. H., J. B. Kadane, and A. O’Hagan (2005).  Statistical methods for eliciting probability distributions.  Journal of the American Statistical Association 100(470), 680–701.
  * Garthwaite and O’Hagan (2000) Garthwaite, P. H. and A. O’Hagan (2000).  Quantifying expert opinion in the UK water industry: an experimental study.  The Statistican 49, 455–477.
  * Gelman et al. (2014) Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2014).  Bayesian Data Analysis (Third ed.).  Boca Raton: CRC Press.
  * Grassia (1977) Grassia, A. (1977).  On a family of distributions with argument between 0 and 1 obtained by transformation of the gamma and derived compound distributions.  Australian Journal of Statistics 19, 108–114.
  * Hirst et al. (2017) Hirst, A. J., K. Giri, D. Ball, and R. S. Lee (2017).  Determination of the physical drivers of Zostera seagrass distribution using a spatial autoregressive lag model.  Marine and Freshwater Research 68, 1752–1763.
  * Hosack (2023) Hosack, G. R. (2023).  eglm: Elicitation for Generalised Linear Models (GLMs) and Extensions.  R package version 0.3.0, https://doi.org/10.5281/zenodo.7772711.
  * Hosack et al. (2017) Hosack, G. R., K. R. Hayes, and S. C. Barry (2017).  Prior elicitation for Bayesian generalised linear models with application to risk control option assessment.  Reliability Engineering & System Safety 167(Supplement C), 351 – 361.
  * Hosack et al. (2021) Hosack, G. R., A. Ickowicz, and K. R. Hayes (2021).  Quantifying the risk of vector-borne disease transmission attributable to genetically modified vectors.  Royal Society Open Science 8, 201525.
  * Jeffreys (1961) Jeffreys, H. (1961).  The Theory of Probability (Third ed.).  Oxford: Oxford University Press.
  * Joe (2015) Joe, H. (2015).  Dependence Modeling with Copulas.  Number 134 in Monographs on Statistics and Applied Probability. Boca Raton: CRC Press.
  * Jørgensen (1983) Jørgensen, B. (1983).  Maximum likelihood estimation and large-sample inference for generalized linear and nonlinear regression models.  Biometrika 70, 19–28.
  * Jørgensen (1987) Jørgensen, B. (1987).  Exponential dispersion models.  Journal of the Royal Statistical Society. Series B (Methodological) 49, 127–162.
  * Jørgensen (1997a) Jørgensen, B. (1997a).  Proper dispersion models.  Brazilian Journal of Probability and Statistics 11, 89–128.
  * Jørgensen (1997b) Jørgensen, B. (1997b).  The Theory of Dispersion Models.  Chapman & Hall, London.
  * Kadane et al. (1980) Kadane, J. B., J. M. Dickey, R. L. Winkler, W. S. Smith, and S. C. Peters (1980).  Interactive elicitation of opinion for a normal linear model.  Journal of the American Statistical Association 75(372), 845–854.
  * Kadane and Wolfson (1998) Kadane, J. B. and L. J. Wolfson (1998).  Experiences in elicitation.  The Statistician 47, 1–20.
  * Kass and Raftery (1995) Kass, R. W. and A. E. Raftery (1995).  Bayes factors.  Journal of the American Statistical Association 90, 773–795.
  * Kosorok (2008) Kosorok, M. R. (2008).  Introduction to Empirical Processes and Semiparametric Inference.  New York: Springer.
  * Kotz and Nadarajah (2004) Kotz, S. and S. Nadarajah (2004).  Multivariate t Distributions and their Applications.  Cambridge University Press.
  * Kuhnert (2011) Kuhnert, P. M. (2011).  Four case studies in using expert opinion to inform priors.  Environmetrics 22, 622–674.
  * Kullback (1959) Kullback, S. (1959).  Information Theory and Statistics (Dover republication, 1997 ed.).  New York: John Wiley and Sons.
  * Kullback and Leibler (1951) Kullback, S. and R. A. Leibler (1951).  On information and sufficiency.  The Annals of Mathematical Statistics 22, 79–86.
  * Kurowicka (2011) Kurowicka, D. (2011).  Optimal truncation of vines.  In D. Kurowicka and H. Joe (Eds.), Dependence Modeling: Vine Copula Handbook, Chapter 11, pp. 233–247. Singapore: World Scientific.
  * Lindley (2014) Lindley, D. V. (2014).  Understanding Uncertainty (Revised ed.).  Wiley.
  * Lindsey (1997) Lindsey, J. K. (1997).  Applying Generalized Linear Models.  Springer-Verlag New York, Inc.
  * McCullagh and Nelder (1989) McCullagh, P. and J. A. Nelder (1989).  Generalized Linear Models (2nd ed.).  Boca Raton, Florida USA: Chapman & Hall/CRC.
  * Mikkola et al. (2021) Mikkola, P., O. A. Martin, S. Chandramouli, M. Hartmann, O. A. Pla, O. Thomas, H. Pesonen, J. Corander, A. Vehtari, S. Kaski, P.-C. Burkner, and A. Klami (2021).  Prior knowledge elicitation: The past, present, and future.  arXiv:2112.01380v1.
  * Morris (1982) Morris, C. N. (1982).  Natural exponential families with quadratic variance functions.  The Annals of Statistics 10, 65–80.
  * Nayar et al. (2018) Nayar, S., M. G. K. Loo, J. E. Tanner, A. R. Longmore, and G. P. Jenkins (2018).  Nitrogen acquisition and resource allocation strategies in temperate seagrass Zostera nigricaulis: Uptake, assimilation and translocation processes.  Scientific Reports 8, 17151.