text:
\section{Posterior means of the patient fixed effects in the Australian GP data} \label{app:estimation result patient}

\begin{table}[ht!]
\centering
\begin{tabular}{|c|rrrrrrrr|}
\hline
\multirow{2}{*}{Variable} & \multicolumn{8}{c|}{Product} \\ \cline{2-9}
 & \multicolumn{1}{c}{1}   & \multicolumn{1}{c}{2}   & \multicolumn{1}{c}{3}   & \multicolumn{1}{c}{4}   & \multicolumn{1}{c}{5}   & \multicolumn{1}{c}{6}   & \multicolumn{1}{c}{7}   & \multicolumn{1}{c|}{8}   \\ \hline
Intercept & \cellcolor{gray!25}1.4026  & \cellcolor{gray!25}-1.2448 & -0.3887 & \cellcolor{gray!25}1.1099  & \cellcolor{gray!25}-2.4041 & -0.1309 & \cellcolor{gray!25}-1.7113 & 0.6777  \\
dagegp1   & \cellcolor{gray!25}0.1974  & -0.1334 & 0.0110  & 0.0752  & \cellcolor{gray!25}-0.5025 & -0.0202 & \cellcolor{gray!25}-0.2892 & 0.0644  \\
dagegp3   & -0.1309 & 0.0614  & -0.0620 & -0.0003 & \cellcolor{gray!25}0.3161  & -0.0032 & 0.0857  & 0.0109  \\
dagegp4   & \cellcolor{gray!25}-0.3881 & \cellcolor{gray!25}0.1814  & \cellcolor{gray!25}-0.2384 & -0.1042 & \cellcolor{gray!25}0.8045  & -0.0268 & \cellcolor{gray!25}0.3775  & 0.0017  \\
drfe2     & -0.0430 & 0.0008  & -0.0372 & -0.0140 & 0.0434  & -0.0189 & -0.0471 & 0.0070  \\
drfe3     & \cellcolor{gray!25}-0.2479 & -0.0529 & 0.0285  & 0.0762  & 0.0914  & 0.1058  & -0.0275 & \cellcolor{gray!25}0.1353  \\
drfe4     & -0.0199 & 0.1012  & -0.0094 & 0.0513  & 0.0666  & 0.0723  & -0.0719 & 0.0054  \\
dbleed1   & 0.0500  & \cellcolor{gray!25}-0.1327 & 0.0604  & -0.0887 & \cellcolor{gray!25}0.3982  & -0.0249 & \cellcolor{gray!25}-0.5218 & \cellcolor{gray!25}-0.2285 \\
dbleed2   & 0.0163  & -0.0748 & 0.0217  & -0.0227 & 0.0068  & 0.0397  & -0.0821 & -0.0248 \\
dbp1      & -0.0591 & 0.0000  & -0.0303 & 0.0288  & 0.0034  & 0.0306  & -0.0235 & \cellcolor{gray!25}-0.1427 \\
dbp3      & \cellcolor{gray!25}-0.9902 & \cellcolor{gray!25}0.2436  & 0.0060  & 0.0128  & \cellcolor{gray!25}0.2350  & \cellcolor{gray!25}-0.2954 & \cellcolor{gray!25}0.2541  & 0.0355  \\
drel1     & 0.0426  & -0.0107 & -0.0933 & -0.0022 & \cellcolor{gray!25}0.1565  & 0.0293  & 0.0290  & \cellcolor{gray!25}-0.3945 \\
drel3     & -0.0144 & 0.0272  & -0.0206 & 0.0004  & -0.0261 & 0.0086  & -0.0185 & 0.0190  \\
drel4     & -0.0916 & 0.0869  & 0.0691  & -0.0006 & -0.0094 & 0.0285  & 0.0033  & \cellcolor{gray!25}-0.2036 \\
dchild1   & \cellcolor{gray!25}-1.7332 & \cellcolor{gray!25}1.3001  & -0.0078 & -0.0903 & \cellcolor{gray!25}0.9163  & \cellcolor{gray!25}-0.9867 & \cellcolor{gray!25}0.5286  & -0.0371 \\
dchild2   & -0.0462 & 0.0350  & -0.0626 & -0.0410 & \cellcolor{gray!25}0.9787  & -0.0504 & \cellcolor{gray!25}0.5930  & -0.0539 \\
dfut1     & \cellcolor{gray!25}-0.3189 & -0.0043 & \cellcolor{gray!25}0.1956  & 0.0249  & \cellcolor{gray!25}0.6284  & -0.0775 & \cellcolor{gray!25}0.2036  & -0.1144 \\
dfut2     & \cellcolor{gray!25}-0.2846 & \cellcolor{gray!25}0.1912  & \cellcolor{gray!25}-0.2151 & \cellcolor{gray!25}-0.1990 & -0.0067 & 0.0366  & -0.1458 & 0.0115  \\
dfut4     & \cellcolor{gray!25}-0.3581 & 0.0496  & 0.0475  & 0.0106  & \cellcolor{gray!25}0.2882  & 0.0059  & 0.0128  & 0.0312  \\
dpil1     & \cellcolor{gray!25}0.4713  & \cellcolor{gray!25}0.3618  & -0.0934 & \cellcolor{gray!25}-0.2596 & -0.0112 & -0.0316 & -0.0441 & -0.0278 \\
dpil3     & \cellcolor{gray!25}-0.1822 & \cellcolor{gray!25}-0.2448 & 0.0300  & 0.0612  & 0.0552  & 0.0351  & 0.0469  & 0.0816  \\
dwt1      & 0.0815  & 0.0368  & \cellcolor{gray!25}-0.2582 & -0.0616 & 0.0318  & 0.0640  & -0.0127 & 0.0833  \\
dcomp2    & \cellcolor{gray!25}-0.3389 & \cellcolor{gray!25}-0.1969 & \cellcolor{gray!25}0.2133  & 0.0628  & \cellcolor{gray!25}0.2300  & -0.0034 & \cellcolor{gray!25}0.3100  & -0.0164 \\
dpay2     & -0.0245 & -0.0540 & -0.0202 & -0.0029 & 0.0080  & 0.0603  & 0.0072  & 0.0075  \\
dpay3     & 0.0321  & -0.0631 & -0.0683 & -0.0167 & -0.0365 & \cellcolor{gray!25}0.2909  & -0.0173 & -0.0046 \\
dsmk2     & \cellcolor{gray!25}-0.2666 & -0.0118 & -0.0266 & -0.0129 & -0.0038 & 0.0422  & 0.0867  & 0.0308  \\
dsmk3     & \cellcolor{gray!25}-0.5183 & -0.0143 & 0.0125  & 0.0248  & 0.0144  & -0.0531 & 0.0461  & 0.0332  \\ \hline
\end{tabular}
\caption{Regression coefficient posterior mean estimates for the attributes of a female patient for various products in the contraceptive discussion data.  Parameters whose 90\% credible interval does not include 0 are shown in grey.}
\label{table:estimation result for women attributes}
\end{table}

\section{Posterior means of the GP fixed effects in the Australian GP data} \label{app:estimation result GP}

\begin{table}[ht!]
\centering
\begin{tabular}{|c|rrrrrrrr|}
\hline
\multirow{2}{*}{Variable} & \multicolumn{8}{c|}{Product} \\ \cline{2-9}
                                                                      & \multicolumn{1}{c}{1}   & \multicolumn{1}{c}{2}   & \multicolumn{1}{c}{3}   & \multicolumn{1}{c}{4}   & \multicolumn{1}{c}{5}   & \multicolumn{1}{c}{6}   & \multicolumn{1}{c}{7}   & \multicolumn{1}{c|}{8}   \\ \hline
Female                                                                         & -0.0643                  & 0.0254                   & \cellcolor{gray!25}-0.4392                  & 0.0716                   & 0.0375                   & \cellcolor{gray!25}0.6089                   & \cellcolor{gray!25}-0.4569                  & -0.0248                  \\ \hline
Fellow                                                                         & -0.0182                  & -0.0990                  & 0.0682                   & 0.0375                   & 0.2107                   & 0.1086                   & -0.1452                  &  -0.0091                  \\  \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Family\\ planning\end{tabular}}     & \multirow{2}{*}{-0.0004} & \multirow{2}{*}{-0.0160} & \multirow{2}{*}{-0.1220} & \multirow{2}{*}{0.2292}  & \multirow{2}{*}{0.0446}  & \multirow{2}{*}{0.0360}  & \multirow{2}{*}{-0.0321} & \multirow{2}{*}{-0.0133} \\
                                                                               &                          &                          &                          &                          &                          &                          &                          &                          \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Bulk-\\ bill\end{tabular}}          & \multirow{2}{*}{-0.0199} & \multirow{2}{*}{-0.0372} & \multirow{2}{*}{0.0425}  & \multirow{2}{*}{-0.0389} & \multirow{2}{*}{-0.0594} & \multirow{2}{*}{ 0.0035} & \multirow{2}{*}{0.0503}  & \multirow{2}{*}{0.0060}  \\
                                                                               &                          &                          &                          &                          &                          &                          &                          &                          \\ \hline
Age                                                                            & 0.0087                   & 0.0079                   & \cellcolor{gray!25}0.0205                   & -0.0064                  & \cellcolor{gray!25}0.0176                   & -0.0043                  & 0.0090                  & -0.0101                  \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Australian\\ graduate\end{tabular}} & \multirow{2}{*}{0.0835}  & \multirow{2}{*}{0.0569}  & \multirow{2}{*}{-0.0080}  & \multirow{2}{*}{0.3457}  & \multirow{2}{*}{0.0989}  & \multirow{2}{*}{-0.2348} & \multirow{2}{*}{-0.0933} &  \cellcolor{gray!25} \\
                                                                               &                          &                          &                          &                          &                          &                          &                          &  \multirow{-2}{*}{\cellcolor{gray!25}0.5586}                        \\ \hline
Urban                                                                          & -0.0947                  & 0.0072                   & 0.0693                   & -0.0073                  & 0.0112                  & 0.0040                   & -0.0285                   & 0.1814                   \\ \hline
\end{tabular}
\caption{Regression coefficient posterior mean estimates for the GP fixed effects for various products in the contraceptive discussion data.  Parameters whose 90\% credible interval does not include 0 are shown in grey.}
\label{table:estimation result for doctors fixed effects}
\end{table}

\newpage
\section{Posterior means of the parameters in \texorpdfstring{$\bm{R}_{\boldsymbol \epsilon}$}{bReps}} \label{app:correlation}

\begin{table}[ht!]
\centering
\begin{tabular}{|c|c|c|}
\hline
Parameter & Posterior Mean & 95\% Credible Interval \\
\hline
$r_{21}$ & -0.1157   & (-0.2058, -0.0251)     \\
$r_{31}$ & -0.0502   & (-0.1446, { }0.0451)   \\
$r_{32}$ & { }0.1617 & ({ }0.0707, { }0.2506) \\
$r_{41}$ & -0.0481   & (-0.1512, { }0.0532)   \\
$r_{42}$ & { }0.0454 & (-0.0578, { }0.1512)   \\
$r_{43}$ & { }0.5891 & ({ }0.5151, { }0.6577) \\
$r_{51}$ & -0.2354   & (-0.3223, -0.1467)     \\
$r_{52}$ & -0.0255   & (-0.1203, { }0.0695)   \\
$r_{53}$ & { }0.1759 & ({ }0.0822, { }0.2657) \\
$r_{54}$ & { }0.2424 & ({ }0.1431, { }0.3398) \\
$r_{61}$ & { }0.4721 & ({ }0.3930, { }0.5472) \\
$r_{62}$ & -0.2702   & (-0.3572, -0.1810)     \\
$r_{63}$ & { }0.0144 & (-0.0753, { }0.1038)   \\
$r_{64}$ & { }0.0358 & (-0.0646, { }0.1367)   \\
$r_{65}$ & -0.0695   & (-0.1629, { }0.0266)   \\
$r_{71}$ & -0.2054   & (-0.3068, -0.1028)     \\
$r_{72}$ & -0.0540   & (-0.1559, { }0.0499)   \\
$r_{73}$ & { }0.1849 & ({ }0.0810, { }0.2877) \\
$r_{74}$ & { }0.1877 & ({ }0.0747, { }0.2973) \\
$r_{75}$ & { }0.5185 & ({ }0.4288, { }0.6016) \\
$r_{76}$ & -0.0051   & (-0.1138, { }0.1014)   \\
$r_{81}$ & -0.0196   & (-0.1147, { }0.0749)   \\
$r_{82}$ & -0.0491   & (-0.1412, { }0.0440)   \\
$r_{83}$ & { }0.0205 & (-0.0704, { }0.1091)   \\
$r_{84}$ & { }0.1044 & ({ }0.0059, { }0.2009) \\
$r_{85}$ & -0.0788   & (-0.1701, { }0.0131)   \\
$r_{86}$ & { }0.1841 & ({ }0.0930, { }0.2753) \\
$r_{87}$ & { }0.2046 & ({ }0.1036, { }0.3040) \\ \hline
\end{tabular}
\caption{Posterior mean estimates of the elements in the correlation matrix $\bm{R}_{\boldsymbol \epsilon}$ of the errors terms with their 95\% credible intervals.}
\label{table:correlation}
\end{table}

\newpage
\section{Posterior means of the parameters in \texorpdfstring{$\boldsymbol \Sigma_{\boldsymbol \alpha}$}{Sigba}} \label{app:covariance}

\begin{table}[ht!]
\centering
\begin{tabular}{|c|c|c|}
\hline
Parameter & Posterior Mean & 95\% Credible Interval \\
\hline
$\sigma^2_{\alpha_1}$ & { }0.5147  & ({ }0.3630, 0.7095) \\
$\sigma^2_{\alpha_2}$ & { }0.6492  & ({ }0.4745, 0.8690) \\
$\sigma^2_{\alpha_3}$ & { }1.3021  & ({ }0.9554, 1.7634) \\
$\sigma^2_{\alpha_4}$ & { }1.5329  & ({ }1.0905, 2.1402) \\
$\sigma^2_{\alpha_5}$ & { }1.0586  & ({ }0.7785, 1.4112) \\
$\sigma^2_{\alpha_6}$ & { }1.4533  & ({ }1.0695, 1.9627) \\
$\sigma^2_{\alpha_7}$ & { }1.9449  & ({ }1.3741, 2.7053) \\
$\sigma^2_{\alpha_8}$ & { }1.2432  & ({ }0.8962, 1.6945) \\
$\sigma_{(\alpha_2, \alpha_1)}$ & { }0.3009  & ({ }0.1813, 0.4417) \\
$\sigma_{(\alpha_3, \alpha_1)}$ & { }0.2733  & ({ }0.1148, 0.4542) \\
$\sigma_{(\alpha_3, \alpha_2)}$ & { }0.3026  & ({ }0.1314, 0.5004) \\
$\sigma_{(\alpha_4, \alpha_1)}$ & { }0.2481  & ({ }0.0731, 0.4489) \\
$\sigma_{(\alpha_4, \alpha_2)}$ & { }0.2659  & ({ }0.0681, 0.4870) \\
$\sigma_{(\alpha_4, \alpha_3)}$ & { }0.2525  & (-0.0200, 0.5434)   \\
$\sigma_{(\alpha_5, \alpha_1)}$ & { }0.0764  & (-0.0677, 0.2272)   \\
$\sigma_{(\alpha_5, \alpha_2)}$ & { }0.1866  & ({ }0.0322, 0.3566) \\
$\sigma_{(\alpha_5, \alpha_3)}$ & -0.0177    & (-0.2376, 0.2022)   \\
$\sigma_{(\alpha_5, \alpha_4)}$ & { }0.5141  & ({ }0.2705, 0.8058) \\
$\sigma_{(\alpha_6, \alpha_1)}$ & { }0.2012  & ({ }0.0305, 0.3969) \\
$\sigma_{(\alpha_6, \alpha_2)}$ & { }0.2195  & ({ }0.0360, 0.4232) \\
$\sigma_{(\alpha_6, \alpha_3)}$ & { }0.1101  & (-0.1427, 0.3672)   \\
$\sigma_{(\alpha_6, \alpha_4)}$ & { }0.2469  & (-0.0342, 0.5528)   \\
$\sigma_{(\alpha_6, \alpha_5)}$ & { }0.2856  & ({ }0.0539, 0.5430) \\
$\sigma_{(\alpha_7, \alpha_1)}$ & { }0.0611  & (-0.1420, 0.2722)   \\
$\sigma_{(\alpha_7, \alpha_2)}$ & { }0.2400  & ({ }0.0277, 0.4768) \\
$\sigma_{(\alpha_7, \alpha_3)}$ & { }0.2573  & (-0.0378, 0.5871)   \\
$\sigma_{(\alpha_7, \alpha_4)}$ & -0.2164    & (-0.5979, 0.1113)   \\
$\sigma_{(\alpha_7, \alpha_5)}$ & { }0.0148  & (-0.2596, 0.2956)   \\
$\sigma_{(\alpha_7, \alpha_6)}$ & { }0.2915  & (-0.0194, 0.6297)   \\
$\sigma_{(\alpha_8, \alpha_1)}$ & { }0.2570  & ({ }0.0982, 0.4404) \\
$\sigma_{(\alpha_8, \alpha_2)}$ & { }0.3344  & ({ }0.1562, 0.5423) \\
$\sigma_{(\alpha_8, \alpha_3)}$ & { }0.0750  & (-0.1594, 0.3157)   \\
$\sigma_{(\alpha_8, \alpha_4)}$ & { }0.2315  & (-0.0290, 0.5223)   \\
$\sigma_{(\alpha_8, \alpha_5)}$ & { }0.2537  & ({ }0.0449, 0.4889) \\
$\sigma_{(\alpha_8, \alpha_6)}$ & { }0.4904  & ({ }0.2432, 0.7802) \\
$\sigma_{(\alpha_8, \alpha_7)}$ & { }0.4139  & ({ }0.1235, 0.7563) \\ \hline
\end{tabular}
\caption{Posterior mean estimates of the elements in the covariance matrix $\boldsymbol \Sigma_{\boldsymbol \alpha}$ of the random effects with their 95\% credible intervals.  The variance parameter in $\boldsymbol \Sigma_{\boldsymbol \alpha}$ is denoted by $\sigma^2_{\alpha_i}$ while the $(i,j)$-th covariance parameter is given by $\sigma_{(\alpha_i, \alpha_j)}$.}
\label{table:covariance}
\end{table}
\section{Introduction} \label{sec:introduction}

Bayesian inference for the multivariate probit (MVP) model is usually performed using the data augmentation representation of \cite{chib1998analysis}, whereby the latent variables indicating the observed outcomes are normally distributed.  For unique identification of the regression parameters, the covariance matrix of these latent normal random variates is assumed to be a correlation matrix $\bm{R}_{\boldsymbol \epsilon}$.  However, Monte Carlo sampling for $\bm{R}_{\boldsymbol \epsilon}$ in a Bayesian context is difficult due to the restrictions on the diagonal entries and the requirement that the matrix $\bm{R}_{\boldsymbol \epsilon}$ must be positive definite.

This article presents three contributions, two  methodological and the third a subject matter one.  The first methodological contribution provides an improved method for sampling the potentially high dimensional correlation matrix $\bm{R}_{\boldsymbol \epsilon}$ within a Markov chain Monte Carlo (MCMC) algorithm.  In order to circumvent the positive definiteness restriction imposed on a correlation matrix, we adopt the reparameterisation strategy of \cite{smith2013bayesian} which re-expresses $\bm{R}_{\boldsymbol \epsilon}$ as an unconstrained Cholesky factor $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$.  This maps the manifold space of a correlation matrix to a Euclidean space, which  improves posterior simulation while keeping the number of unknown parameters the same.  A prior distribution is then specified on $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ such that the implied marginal densities of the correlation coefficients are uniform on $(-1,1)$.  We employ the Hamiltonian Monte Carlo (HMC) algorithm \citep{neal2011mcmc} to sample the high dimensional $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ efficiently, thereby avoiding the slow exploration of parameter space by random walk updates as in \cite{smith2013bayesian}.

The second methodological contribution is to introduce antithetic sampling, based on the work of \cite{hammersley1956new}, into the Metropolis-Hastings (MH) literature.  In order to implement this idea, we specify the proposal distribution of parameter update as a deterministic function.  Here, the generated samples will be super-efficient in terms of the reduction in variance of the Monte Carlo estimates compared to the same estimates constructed from uncorrelated samples.  Although the chain update proposal is deterministic, the convergence properties are not compromised when this is embedded within a larger system of MCMC sampling.  Our proposed methodology is motivated by the over-relaxation algorithm \citep{adler1981over, barone1990improving}, and is similar to the idea built within the framework of HMC in \cite{pakman2014exact}.  However, our proposed sampler is different from these methods in two main aspects.  First, there is no randomness in the proposal distribution for parameter updates in our method, whereas theirs still retain a certain degree of stochasticity.  Second, we introduce perfect negative correlation between successive MCMC samples via the deterministic proposal, while they suggest partial or zero dependence between the samples.  Results based on our real data application document a significant improvement of up to a 16 times performance gain in the mixing behaviour of the Markov chain, thereby lowering the autocorrelation between the iterates.  The computing time of the algorithm is also marginally reduced due to the deterministic sampling.

Our methodological development is motivated by the staged stated preference panel data collection described in \cite{fiebig2017consideration}, which is used to study the decision-making of Australian general practitioners (GPs) about female contraceptive products.  Here, the authors used the data from the third and final stage, whereas we explore outcomes from the second stage.  This second stage relates to the question of which particular contraceptive products GPs would discuss with a female patient, defined by a vignette that is part of the experimental design.  Separate univariate analyses on each product would ignore possible complex dependence structures that are useful in exploring which particular bundles of products are discussed with patients.  This is important here because in any correlated choice problem there may be multiple close substitutes, which makes joint rather than marginal probabilities more relevant.  Therefore, we model the GPs' choices by an MVP model.  Inspection of the resulting graphical model describing this interaction between products lends support to the suitability of a multivariate approach.  By using the MVP model, we are able to compute the joint probability of specific product bundles being discussed with a patient.  Posterior estimation of this probability, based on a patient with certain socio-economic and clinical characteristics, reveals differing views among the GPs in the sample on the suitability of long acting contraceptive choices.  This variability is known as medical practice variation in the health industry, whereby the decision making of GPs is influenced by both their personal characteristics such as gender, age and qualifications, as well as other unobservables that we model as random effects.

The rest of the paper is organised as follows. Section~\ref{sec:multivariate} describes the MVP model with random effects and reviews previous research associated with sampling $\bm{R}_{\boldsymbol \epsilon}$.  Section~\ref{sec:efficient} presents our proposed methodology of sampling $\bm{R}_{\boldsymbol \epsilon}$, and Section~\ref{sec:deterministic} outlines the antithetic sampling technique whose efficiency is illustrated via simulation studies in Section~\ref{sec:simulation}.  Section~\ref{sec:application} provides our analysis of the discussion preference data of contraceptive products by Australian GPs, and Section~\ref{sec:conclusion} concludes.  Appendices~\ref{app:sampling scheme}--\ref{app:covariance} provide further details on the contraceptive product data analysis.

\section{Multivariate probit model with random effects} \label{sec:multivariate}

The MVP model has been used extensively to model correlated binary data \citep{gibbons1998health, buchmueller2013preference}.  Let $\bm{y}_{it}=(y_{1,it}, \dotsc, y_{D,it})^\top$ be a vector of $D$ correlated binary outcomes for individual $i=1,\dotsc,P$ at time period $t$, for $t=1, \dotsc, T$.  The latent variable representation of the MVP model, using the data augmentation approach of \cite{albert1993bayesian}, is given by
\begin{gather}
\bm{y}^\ast_{it} = \boldsymbol \alpha_i + \bm{B} \bm{x}_{it} + \boldsymbol \epsilon_{it},
\label{eqn:MVP} \\
\boldsymbol \alpha_i = (\alpha_{1,i}, \dotsc, \alpha_{D,i})^\top \stackrel{iid}{\sim} \mathcal{N}({\bf 0}, \boldsymbol \Sigma_{\boldsymbol \alpha}), \\
\boldsymbol \epsilon_{it} = (\epsilon_{1,it}, \dotsc, \epsilon_{D,it})^\top \stackrel{iid}{\sim} \mathcal{N}({\bf 0}, \bm{R}_{\boldsymbol \epsilon}),
\end{gather}
for $i=1,\dotsc,P, t=1,\dotsc,T,$ where $\bm{y}^\ast_{it}=(y^\ast_{1,it}, \dotsc, y^\ast_{D,it})^\top$ is a continuous latent variable, $\boldsymbol \alpha_i$ is a $D$-vector of outcome-specific random effects for individual $i$ allowing for heterogeneity between individuals, $\bm{x}_{it}=(1, x_{1,it}, \dotsc, x_{K-1,it})^\top$ is an exogenous variable, $\bm{B}$ is a $D \times K$ matrix of regression coefficients and $\boldsymbol \epsilon_{it}$ is a $D$-vector correlated error term which models the dependence structure between outcomes.  The variable $\bm{x}_{it}$ is assumed to be uncorrelated with both $\boldsymbol \alpha_i$ and $\boldsymbol \epsilon_{it}$.  This is entirely appropriate in the stated preference case that is our motivating analysis but relaxing the assumption of exogenous $\bm{x}_{it}$ represents a useful extension.  In order for $\bm{B}$ to be uniquely identified \citep{chib1998analysis}, $\bm{R}_{\boldsymbol \epsilon}$ is set to be a correlation matrix.  The observed outcome $\bm{y}_{it}$ is defined to be dependent on the latent variable $\bm{y}^\ast_{it}$ via the relationship
\begin{equation}
y_{d,it} = \mathds{1}(y^\ast_{d,it}>0), \quad d=1, \dotsc, D,
\label{eqn:latent}
\end{equation}
where $\mathds{1}(\bm E)$ is an indicator function which takes value 1 if the event $\bm E$ occurs and 0 otherwise.  Let $\bm{y}=\{\bm{y}_{it}; i=1,\dotsc,P, t=1,\dotsc,T\}$ be the set of observed discrete outcomes.  The density of the latent variables $\bm{y}^\ast$ conditional on the random effects $\boldsymbol \alpha_{1:P}=(\boldsymbol \alpha_1,\dotsc,\boldsymbol \alpha_P)$ is given by
\begin{equation}
p(\bm{y}^\ast|\boldsymbol \alpha_{1:P}, \boldsymbol \theta) = \prod_{i=1}^P \prod_{t=1}^T \phi(\bm{y}^\ast_{it}; \boldsymbol \mu_{it}, \bm{R}_{\boldsymbol \epsilon}),
\label{eqn:likelihood}
\end{equation}
where $\boldsymbol \theta:=(\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \Sigma_{\boldsymbol \alpha})$ denotes the vector of model parameters, $\boldsymbol \mu_{it}=\boldsymbol \alpha_i + \bm{B} \bm{x}_{it}$ and $\phi$ is the multivariate normal density function.

Following the specification of the MVP model in (\ref{eqn:MVP})--(\ref{eqn:latent}), the posterior density is
\begin{equation}
\pi(\bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta|\bm{y}) = \frac{p(\bm{y}|\bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta)p(\bm{y}^\ast|\boldsymbol \alpha_{1:P},\boldsymbol \theta)p(\boldsymbol \alpha_{1:P}|\boldsymbol \theta)p(\boldsymbol \theta)}{p(\bm{y})},
\label{eqn:posterior}
\end{equation}
where $p(\bm{y})$ is the marginal likelihood, $p(\boldsymbol \theta)$ is the prior on the model parameters $\boldsymbol \theta$ and
\begin{equation}
p(\bm{y}|\bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta) = \prod_{i=1}^P \prod_{t=1}^T \prod_{d=1}^D \bigg(\mathds{1}(y_{d,it}=0)\mathds{1}(y^\ast_{d,it} \leq 0) + \mathds{1}(y_{d,it}=1)\mathds{1}(y^\ast_{d,it} > 0)\bigg).
\label{eqn:augmented posterior}
\end{equation}
Useful conjugate priors are available for $\bm{B}$ (or $\boldsymbol \beta=\text{vec}(\bm{B})$) and $\boldsymbol \Sigma_{\boldsymbol \alpha}$ which simplifies MCMC sampling, but it is difficult to posit a suitable prior for $\bm{R}_{\boldsymbol \epsilon}$.

\subsection{Prior choice for the correlation matrix \texorpdfstring{$\bm{R}_{\boldsymbol \epsilon}$}{bReps}} \label{subsec:prior}

\cite{barnard2000modeling} decompose a covariance matrix $\Sig_{\beps}$ as $\bm{S}\bm{R}_{\boldsymbol \epsilon}\bm{S},$ where $\bm{S}$ is a diagonal matrix of standard deviations and $\bm{R}_{\boldsymbol \epsilon}$ is a correlation matrix.  They show that if $\Sig_{\beps} \sim \mathcal{IW}(\nu, \bm{I})$, i.e. an inverse-Wishart distribution with degrees of freedom $\nu$ and scale matrix $\bm{I}$, then the density of $\bm{R}_{\boldsymbol \epsilon}$ is
\begin{equation}
p(\bm{R}_{\boldsymbol \epsilon}) \propto |\bm{R}_{\boldsymbol \epsilon}|^{\frac{1}{2}(\nu-1)(D-1)-1}\Bigg(\prod_{i=1}^D|\bm{R}_{\boldsymbol \epsilon}(-i;-i)|\Bigg)^{-\frac{\nu}{2}},
\label{eqn:marginally uniform prior}
\end{equation}
where $\bm{R}_{\boldsymbol \epsilon}(-i;-i)$ denotes the $i$-th principal submatrix of $\bm{R}_{\boldsymbol \epsilon}$, that is $\bm{R}_{\boldsymbol \epsilon}$ with its $i$-th row and column removed.  We follow \cite{barnard2000modeling} and take (\ref{eqn:marginally uniform prior}) as the prior for $\bm{R}_{\boldsymbol \epsilon}$, which induces a modified Beta distribution on each off-diagonal element $r_{ij}$ of $\bm{R}_{\boldsymbol \epsilon}, i \neq j$.  In particular, the marginal densities of the $r_{ij}$ are uniform on $(-1,1)$ when $\nu=D+1$, which means that posterior inference is invariant to the ordering of the binary outcomes $\bm{y}$.  Furthermore, recent results in \cite{wang2018on} establish that for such a choice of $\nu$, the corresponding matrix of partial correlations $\rho_{kl}$ has the LKJ distribution of \cite{lewandowski2009generating} with unit shape parameter.  This means that all $\rho_{kl}$ are marginally distributed according to a $\textrm{Beta}(\frac{D}{2}, \frac{D}{2})$ distribution over $(-1,1)$ with both shape parameters $\frac{D}{2}$, which is informative in high dimensions because the Beta density increasingly concentrates around zero.  The informativity of $\rho_{kl}$ is useful in practical applications, where more often than not a sparse structure on the partial correlation matrix is desirable to suggest conditional independence.

The dependence structures imposed by the marginally uniform prior are less studied in the literature.  Since analytical results for these properties are limited \citep{tokuda2011visualizing}, we briefly illustrate these graphically instead.  The results obtained are based on correlation matrices of dimension $D=4$ but they can be generalised to higher dimensions.  We generate $10^7$ samples from (\ref{eqn:marginally uniform prior}) with $\nu=D+1$ by normalising the covariance matrices drawn from an $\mathcal{IW}(D+1, \bm{I})$ distribution.  Figure~\ref{fig:corr_vs_par_corr} illustrates the pairwise dependence structures among the correlations $r_{ij}$ and the partial correlations $\rho_{kl}$ when the pairs share (top panels) or do not share (bottom panels) common indices. When there is a shared index, the density on $(r_{12}, r_{13})$ tends to support similar values in absolute terms (the visible cross pattern), which is less apparent when there is no common index in $(r_{12}, r_{34})$. However, both distributions have most of their density on the vertices corresponding to $|r_{ij}|\approx1$.  This means that inference for all pairs of $r_{ij}$ is skewed towards jointly extreme values a priori (the univariate margin for each $r_{ij}$ is still uniform on $(-1,1)$), although this effect diminishes with an increase in the number of observations.  In contrast, pairs of partial correlations $\rho_{kl}$ exhibit no dependence structure regardless of whether or not there is a common index.  Independence is also observed between $r_{ij}$ and $\rho_{kl}$, except when both parameters have the same indices $(r_{12},\rho_{12})$ in which case they are strongly positively correlated.
\begin{figure}[t!]
    \centering
    \includegraphics[trim= 1cm 3.5cm 1.25cm 5cm,clip,width=\textwidth]{corr_vs_par_corr}
    \caption{Bivariate density plots showing the dependence structures associated with the marginally uniform prior (\ref{eqn:marginally uniform prior}) on $\bm{R}_{\boldsymbol \epsilon}$ with $\nu=D+1$, for pairs of parameters sharing common indices (top panels) and without a common index (bottom panels).}
    \label{fig:corr_vs_par_corr}
\end{figure}

We now discuss related work on priors for $\bm{R}_{\boldsymbol \epsilon}$.  Let $\mathcal{R}^D$ be the space of all valid correlation matrices.  \cite{barnard2000modeling} also suggest a uniform prior over all correlation matrices in $\mathcal{R}^D$, which is equivalent to the LKJ prior with unit shape, as suggested by the \cite{stan2017stan}.  Note that the induced prior on the partial correlation matrix is the marginally uniform prior in (\ref{eqn:marginally uniform prior}) with $\nu=D+1$ (Figure \ref{fig:corr_vs_par_corr}).  This might not be a suitable prior for $\rho_{kl}$ since, as discussed above, this joint distribution for $\rho_{kl}$ exhibits dependence and has large mass on extreme values.  \cite{chib1998analysis} propose using a multivariate normal prior on the $r_{ij}$, with the support of the prior restricted to values of $r_{ij}$ which give a correlation matrix in $\mathcal{R}^D$, while \cite{liechty2004bayesian} introduce a mixture of normal distributions prior on $r_{ij}$ to express a priori knowledge of blocked structure in $\bm{R}_{\boldsymbol \epsilon}$.  However, these choices of normal priors do not imply that all marginal densities of the $r_{ij}$ are the same due to the constraints imposed on the $r_{ij}$ for the resulting $\bm{R}_{\boldsymbol \epsilon}$ to be in $\mathcal{R}^D$.

\subsection{Posterior sampling of \texorpdfstring{$\bm{R}_{\boldsymbol \epsilon}$}{bReps}} \label{subsec:posterior}

Posterior simulation for $\bm{R}_{\boldsymbol \epsilon}$ is challenging for two reasons: (i) the diagonal elements of $\bm{R}_{\boldsymbol \epsilon}$ must be 1 and, (ii) $\bm{R}_{\boldsymbol \epsilon}$ must be positive definite.  \cite{chib1998analysis} suggest sampling the $r_{ij}$ elements of $\bm{R}_{\boldsymbol \epsilon}$ in blocks using a random walk Metropolis-Hastings (RWMH) algorithm with a multivariate $t$ proposal density.  However, the resulting matrix obtained after each proposal is not guaranteed to be a valid correlation matrix in addition to the RWMH algorithm being notorious for its slow exploration of the parameter space.  Tuning the parameters of this proposal distribution also requires finding an approximate mode of the log posterior distribution and the observed Fisher information for every iteration, resulting in high computational overheads.  In the setting of hierarchical regression models, \cite{barnard2000modeling} adopt the Griddy-Gibbs sampler of \cite{ritter1992facilitating} to sample $\bm{R}_{\boldsymbol \epsilon}$.  Here, prior to the Gibbs step, one needs to solve a quadratic equation to determine the support for a single $r_{ij}$ (while keeping the rest fixed) which results in a valid correlation matrix.  The authors document the clear inefficiency in this sampling scheme when the prior in (\ref{eqn:marginally uniform prior}) is used due to its tendency to place more weight on the edges of $\mathcal{R}^D$ space.  Moreover, the design of drawing one $r_{ij}$ at a time becomes computationally prohibitive when $D$ is large.

\section{Efficient sampling for \texorpdfstring{$\bm{R}_{\boldsymbol \epsilon}$}{bReps} when using a marginally uniform prior} \label{sec:efficient}

This section describes an efficient way of sampling $\bm{R}_{\boldsymbol \epsilon}$ by utilising Hamiltonian dynamics \citep{duane1987hybrid}.  This involves reparameterising $\bm{R}_{\boldsymbol \epsilon}$ to enable sampling of parameters in an unconstrained space.  Due to the attractive properties of the marginally uniform prior in (\ref{eqn:marginally uniform prior}) with $\nu=D+1$ discussed in Section~\ref{subsec:prior}, we will use this prior hereafter.  Inference for the posterior distribution in (\ref{eqn:posterior}) can be performed using a Gibbs sampler (see Chapter 10 of \cite{greenberg2012introduction} for details). Our focus here is on the following non-standard conditional posterior distribution
\begin{equation}
\pi(\bm{R}_{\boldsymbol \epsilon}|\bm{y}, \bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta_{-\bm{R}_{\boldsymbol \epsilon}}) \propto \prod_{i=1}^P \prod_{t=1}^T \phi(\bm{y}^\ast_{it}; \boldsymbol \mu_{it}, \bm{R}_{\boldsymbol \epsilon}) \cdot p(\bm{R}_{\boldsymbol \epsilon}),
\label{eqn:conditional posterior}
\end{equation}
where $\boldsymbol \theta_{-\bm{\mathcal{S}}}$ is defined as $\boldsymbol \theta$, but excluding the parameters $\bm{\mathcal{S}}$.

\subsection{An unconstrained parameterisation} \label{subsec:unconstrained}

Because of the restrictions on sampling correlation coefficients on a confined space, we adopt the reparameterisation strategy in \cite{smith2013bayesian} which re-expresses $\bm{R}_{\boldsymbol \epsilon}$ via a positive definite matrix $\Sig_{\beps}$ as
\begin{equation}
\bm{R}_{\boldsymbol \epsilon} = \boldsymbol \Lambda_{\boldsymbol \epsilon}^{-1/2} \Sig_{\beps} \boldsymbol \Lambda_{\boldsymbol \epsilon}^{-1/2},
\label{eqn:reparameterisation}
\end{equation}
where $\boldsymbol \Lambda_{\boldsymbol \epsilon}=\textrm{diag}(\Sig_{\beps})$.  The covariance matrix $\Sig_{\beps}$ can then be written in terms of its Cholesky factorisation $\Sig_{\beps}=\bm{\mathcal{L}}_{\boldsymbol \epsilon}\Leps^\top,$ where $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ is a lower triangular matrix.  The diagonal elements of $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ are set to 1 so that the transformation of $\bm{R}_{\boldsymbol \epsilon}$ to $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ is one-to-one.  We define an operator $\textrm{vechL}$ which vectorises the strict lower triangle of a matrix by row.  The unknown parameter $\textrm{vechL}(\bm{\mathcal{L}}_{\boldsymbol \epsilon})=\{L_{ij};i=2, \dotsc, D, j<i\}$ lies in $\mathbb{R}^{D(D-1)/2}$ and is therefore unconstrained.  \cite{lindstrom1988newton} also implement the Cholesky factorisation on a covariance matrix to optimise the log-likelihood function of a linear mixed effects model.  Other possible reparameterisation methods for $\bm{R}_{\boldsymbol \epsilon}$ include using polar coordinates \citep{rapisarda2007parameterizing} and partial autocorrelations \citep{daniels2009modeling}, but we adopt the representation in (\ref{eqn:reparameterisation}) due to its computational tractability.

By using a change of variables, we can rewrite the density function in (\ref{eqn:conditional posterior}) in terms of $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ as
\begin{equation}
\pi(\bm{\mathcal{L}}_{\boldsymbol \epsilon}|\bm{y}, \bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta_{-\bm{\mathcal{L}}_{\boldsymbol \epsilon}}) \propto \pi(\bm{R}_{\boldsymbol \epsilon}|\bm{y}, \bm{y}^\ast, \boldsymbol \alpha_{1:P}, \boldsymbol \theta_{-\bm{R}_{\boldsymbol \epsilon}}) \cdot |{\bf J}|,
\label{eqn:unconstrained conditional posterior}
\end{equation}
where $|{\bf J}|=|\partial \textrm{vechL}(\bm{R}_{\boldsymbol \epsilon}) / \partial \textrm{vechL}(\bm{\mathcal{L}}_{\boldsymbol \epsilon})^\top|$ is the determinant of the Jacobian for the transformation.  We now note that for the transformation from $\bm{R}_{\boldsymbol \epsilon}$ to $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$, the prior on lower triangular Cholesky factor $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ whose diagonal entries are all fixed as ones, given by
\begin{equation}
p(\bm{\mathcal{L}}_{\boldsymbol \epsilon}) \propto p(\bm{R}_{\boldsymbol \epsilon}) \cdot |{\bf J}|,
\label{eqn:prior cholesky}
\end{equation}
induces a marginally uniform prior on all $r_{ij}$ for $\nu=D+1$.

\subsection{Sampling the Cholesky factor using HMC} \label{subsec:sampling}

HMC, popularised by \cite{neal2011mcmc}, has enjoyed considerable recent interest within the statistical literature due to its ability to generate credible but distant candidate parameters for the MH algorithm, thereby reducing autocorrelation in the posterior samples.  It does so by exploiting gradient information of the log posterior density to simulate a trajectory according to physical dynamics.

Given a target distribution of interest $\pi(\boldsymbol \vartheta)$, which in our case is the density in (\ref{eqn:unconstrained conditional posterior}), HMC introduces a fictitious momentum variable $\bm{u}$ into the physical system, which is assumed to follow a $\mathcal{N}({\bf 0}, \bm{M})$ pseudo-prior and targets the augmented distribution
\begin{equation}
\pi(\boldsymbol \vartheta, \bm{u}) \propto \exp(-\mathcal{H}(\boldsymbol \vartheta, \bm{u})),
\label{eqn:HMC}
\end{equation}
where $\mathcal{H}(\boldsymbol \vartheta, \bm{u}) = -\log \pi(\boldsymbol \vartheta) + \frac{1}{2} \bm{u}^\top \bm{M} ^{-1} \bm{u}$ is termed the Hamiltonian which is made up of potential energy and kinetic energy components.  The potential energy is derived from minus the log density of $\boldsymbol \vartheta$ under the target distribution while the kinetic energy is due to the movement of the momentum variable $\bm{u}$.  The Hamiltonian system is used to describe the evolution of $\boldsymbol \vartheta$ and $\bm{u}$ over time $t$ via the differential equations
\begin{equation}
\frac{d \boldsymbol \vartheta}{dt} = \frac{\partial \mathcal{H}}{\partial \bm{u}} \quad \textrm{and} \quad \frac{d \bm{u}}{dt} = -\frac{\partial \mathcal{H}}{\partial \boldsymbol \vartheta}.
\label{eqn:hamiltonian dynamics}
\end{equation}
The dynamics in (\ref{eqn:hamiltonian dynamics}) can be implemented in practice using the leapfrog method \citep{neal2011mcmc} and discretising continuous time by a stepsize $\varepsilon$ so that
\begin{equation}
\begin{aligned}
\bm{u}(t + \varepsilon/2) &= \bm{u}(t) - (\varepsilon/2) \frac{\partial \mathcal{H}}{\partial \boldsymbol \vartheta}(\boldsymbol \vartheta(t)) \\
\boldsymbol \vartheta(t+\varepsilon) &= \boldsymbol \vartheta(t) + \varepsilon \frac{\partial \mathcal{H}}{\partial \bm{u}}(\bm{u}(t+\varepsilon/2)) \\
\bm{u}(t+\varepsilon) &= \bm{u}(t+\varepsilon/2) - (\varepsilon/2) \frac{\partial \mathcal{H}}{\partial \boldsymbol \vartheta}(\boldsymbol \vartheta(t+\varepsilon)).
\end{aligned}
\label{eqn:leapfrog}
\end{equation}
\cite{neal2011mcmc} shows that properties of the Hamiltonian such as reversibility and volume preservation are maintained under the symplectic integrator in (\ref{eqn:leapfrog}).  Proposed values $\boldsymbol \vartheta'$ and $\bm{u}'$ obtained after a trajectory length of $\mathcal{T}=n\varepsilon$ by iterating procedures in (\ref{eqn:leapfrog}) $n$ times are then accepted with probability $\min \{1, \exp(\mathcal{H}(\boldsymbol \vartheta, \bm{u})-\mathcal{H}(\boldsymbol \vartheta', \bm{u}')) \}$. The invariant distribution of the Markov chain generated from the HMC algorithm is $\pi(\boldsymbol \vartheta, \bm{u})$ and samples from $\pi(\boldsymbol \vartheta)$ can be obtained by marginalising out the momentum $\bm{u}$.

In order to implement the HMC algorithm as described above, computation of the derivatives of (\ref{eqn:unconstrained conditional posterior}) with respect to the $L_{ij}$ is required for the leapfrog update.  Lemma~\ref{lemma:derivatives for gradient} derives the expressions for these gradients.
\begin{lemma}
Let $\bm{E}_k$ denote the matrix obtained by removing column $k$ from an identity matrix $\bm{I}$.  For the parameterisation of $\bm{R}_{\boldsymbol \epsilon}$ in (\ref{eqn:reparameterisation}),
\begin{enumerate}[font=\normalfont, label=(\roman*)]
\item $\displaystyle \frac{\partial \bm{R}_{\boldsymbol \epsilon}^{-1}}{\partial L_{ij}}=-\boldsymbol \Lambda_{\boldsymbol \epsilon}^{1/2} \bigg( \Sig_{\beps}^{-1} \frac{\partial \Sig_{\beps}}{\partial L_{ij}} \Sig_{\beps}^{-1} + \frac{\partial \boldsymbol \Lambda_{\boldsymbol \epsilon}^{-1/2}}{\partial L_{ij}} \boldsymbol \Lambda_{\boldsymbol \epsilon}^{1/2} \Sig_{\beps}^{-1} + \Sig_{\beps}^{-1} \boldsymbol \Lambda_{\boldsymbol \epsilon}^{1/2} \frac{\partial \boldsymbol \Lambda_{\boldsymbol \epsilon}^{-1/2}}{\partial L_{ij}} \bigg) \boldsymbol \Lambda_{\boldsymbol \epsilon}^{1/2}$.
\item $\displaystyle \frac{\partial \log|\bm{R}_{\boldsymbol \epsilon}(-k;-k)|}{\partial L_{ij}}=\text{tr}\bigg(\bm{R}_{\boldsymbol \epsilon}^{-1}(-k;-k) \bm{E}_k^\top \frac{\partial \bm{R}_{\boldsymbol \epsilon}}{\partial L_{ij}} \bm{E}_k\bigg)$.
\item $\displaystyle \frac{\partial \log |\bm{R}_{\boldsymbol \epsilon}|}{\partial L_{ij}}= - \frac{2L_{ij}}{\sum_{k=1}^i L^2_{ik}}$.
\end{enumerate}
\label{lemma:derivatives for gradient}
\end{lemma}
\begin{proof}
Lemma 1(i) and (ii) are respectively obtained using Theorems 1 and 2 in Chapter 8 of \cite{magnus1999matrix}, by expressing $\frac{\partial \Sig_{\beps}^{-1}}{\partial L_{ij}}$ in terms of $\frac{\partial \bm{R}_{\boldsymbol \epsilon}^{-1}}{\partial L_{ij}}$ using the chain rule, and writing $\bm{R}_{\boldsymbol \epsilon}(-k;-k)$ as $\bm{E}_k^\top \bm{R}_{\boldsymbol \epsilon} \bm{E}_k$.  Lemma 1(iii) is straightforward by noting that $|\bm{R}_{\boldsymbol \epsilon}|=|\boldsymbol \Lambda_{\boldsymbol \epsilon}|^{-1}$ since $|\boldsymbol \Sigma_\epsilon|=1$ from its Cholesky decomposition.
\end{proof}

\section{A deterministic proposal distribution} \label{sec:deterministic}

Various strategies have been proposed to reduce the variability in the Monte Carlo estimate of the expectation $\mathbb{E}[f(\boldsymbol \vartheta)]$ of a scalar function $f$ of parameter $\boldsymbol \vartheta$ with respect to some posterior distribution $\pi(\boldsymbol \vartheta)$, including the Rao-Blackwellisation \citep{robert2004monte} and the control variates \citep{dellaportas2012control, oates2017control}.  These techniques produce an efficient estimator of $\mathbb{E}[f(\boldsymbol \vartheta)]$ based on sampled $\boldsymbol \vartheta$ generated from an MCMC sampler.

Here, we focus on a particular class of methods which integrate variance reduction techniques dynamically within an MCMC sampling algorithm.  Let $\boldsymbol \vartheta=(\vartheta_1, \dotsc, \vartheta_n)^\top$ be a parameter vector with normal full conditional distributions $\vartheta_i|\boldsymbol \vartheta_{-i} \sim \mathcal{N}(\mu_i,\sigma_i^2)$, where the conditional mean $\mu_i$ and the conditional variance $\sigma^2_i$ may depend on $\boldsymbol \vartheta_{-i} = \{\vartheta_j:j=1, \dotsc,n, j \not = i\}$.  \cite{adler1981over} and \cite{barone1990improving} introduce an over-relaxation method where the update on $\boldsymbol \vartheta$ is performed using Gibbs sampling, and where the new value $\vartheta'_i$ for each margin of $\boldsymbol \vartheta$ is generated as
\begin{equation}
    \vartheta'_i = (1+\kappa)\mu_i - \kappa \vartheta_i + u \sigma_i \sqrt{1-\kappa^2}, \quad i=1,\dotsc,n,
    \label{eqn:over-relaxation}
\end{equation}
with $u\sim \mathcal{N}(0,1)$ being a standard normal random variable. Equation (\ref{eqn:over-relaxation}) allows for the introduction of dependence between successive samples via the constant antithetic parameter $\kappa$, which is required to be in the open interval $(-1,1)$ so that the Markov chain is ergodic and produces $\pi(\boldsymbol \vartheta)$ as its stationary distribution.  This scheme is exactly the conventional Gibbs sampler when $\kappa=0$. Variance reduction in estimating $\mathbb{E}[f(\boldsymbol \vartheta)]$ is achieved through the antithetic variable method \citep{hammersley1956new} by setting $\kappa>0$ so that the estimation bias in the previous sample is corrected in the opposite direction.   The rate of convergence for the over-relaxation method in (\ref{eqn:over-relaxation}) is studied in \cite{barone1990improving}, while \cite{green1992metropolis} establish that the asymptotic variance of the estimator for $\mathbb{E}[f(\boldsymbol \vartheta)]$ using this strategy for linear $f$ is proportional to $\frac{1-\kappa}{1+\kappa}$.

The inefficiency of an MCMC sampler in estimating $\mathbb{E}[f(\boldsymbol \vartheta)]$ is usually measured by the integrated autocorrelation time \citep{roberts2009examples}, which is defined as
\begin{equation*}
\textrm{IACT}_f = 1 + \sum_{j=1}^\infty \rho_{j,f},
\end{equation*}
where $\rho_{j,f}$ is the lag $j$ autocorrelation function of the MCMC iterates of $f(\boldsymbol \vartheta)$ after convergence.  Alternatively, one can measure the efficiency of the sampler by computing the effective sample size per MCMC iteration, which by definition is the reciprocal of the IACT.  A small value of the IACT is desirable in practice as it indicates that the Markov chain mixes well.  Motivated by the over-relaxation sampler and noting that the IACT can be less than 1 if some of the autocorrelations are negative, in which case a Monte Carlo estimator constructed is super-efficient, we introduce into the MH literature a deterministic design of the proposal distribution for $\boldsymbol \vartheta$
\begin{equation}
q(\boldsymbol \vartheta'|\boldsymbol \vartheta) = \delta_{\psi(\boldsymbol \vartheta)}(\boldsymbol \vartheta'),
\label{eqn:deterministic proposal}
\end{equation}
where $\psi$ is a mapping function which introduces negative correlation between samples and $\delta_{\psi(\boldsymbol \vartheta)}$ is the Dirac delta function at $\psi(\boldsymbol \vartheta)$.  In this case, the MH acceptance probability involves the ratio of $\pi(\boldsymbol \vartheta)$ evaluated at $\boldsymbol \vartheta'$ and $\boldsymbol \vartheta$.

When $\pi(\boldsymbol \vartheta)$ is a normal distribution, we propose setting
\begin{equation}
\psi(\boldsymbol \vartheta) = 2 \boldsymbol \mu_{\boldsymbol \vartheta} - \boldsymbol \vartheta,
\label{eqn:antithetic}
\end{equation}
where $\boldsymbol \mu_{\boldsymbol \vartheta}$ is the mean of $\pi(\boldsymbol \vartheta)$.  It is clear that (\ref{eqn:antithetic}) represents an example of the antithetic variable with perfect negative correlation, and also an instance of the over-relaxation method in (\ref{eqn:over-relaxation}) with $\kappa=1$, which is outside the range of values for which the Markov chain is ergodic.  Symmetry of the normal density gives $\pi(\boldsymbol \vartheta')=\pi(\boldsymbol \vartheta)$, which in turn translates to an acceptance probability of one.  Clearly, our proposed antithetic sampling will only yield an ergodic Markov chain when it is coupled with stochastic simulation of additional parameters that affect the value of the deterministic proposal $\psi(\boldsymbol \vartheta)$, in particular $\boldsymbol \mu_{\boldsymbol \vartheta}$.  Under this condition, the value of $\boldsymbol \mu_{\boldsymbol \vartheta}$ changes in every iteration of the update and this drives the exploration of $\boldsymbol \vartheta$ in the parameter space.  Furthermore, the dependence between $\boldsymbol \vartheta$ and other model parameters prevents exact periodicity from occurring, and thus the Markov chain is aperiodic.

The conditional posterior distribution of the random effects $\boldsymbol \alpha_{1:P}$ in our MVP model is normal and likewise for the regression parameters $\boldsymbol \beta$ when using a conjugate prior.  Therefore, we can employ the antithetic sampling method in (\ref{eqn:antithetic}) to improve the IACTs of $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$.  In fact, antithetic sampling of normal random variables can also be understood in terms of a HMC update.  Suppose that $\boldsymbol \vartheta \sim \mathcal{N}(\boldsymbol \mu_{\boldsymbol \vartheta}, \boldsymbol \Sigma_{\boldsymbol \vartheta})$, and the prior on the momentum variable $\bm{u}$ is chosen as $\mathcal{N}({\bf 0}, \boldsymbol \Sigma^{-1}_{\boldsymbol \vartheta})$.  \cite{pakman2014exact} show that the resulting Hamiltonian system can be solved analytically, with solution given by
\begin{equation}
\boldsymbol \vartheta(t) = \boldsymbol \mu_{\boldsymbol \vartheta} + \boldsymbol \Sigma_{\boldsymbol \vartheta} \bm{u}(0) \sin(t) + (\boldsymbol \vartheta(0)-\boldsymbol \mu_{\boldsymbol \vartheta}) \cos(t),
\label{eqn:solution for Hamilton}
\end{equation}
which is a linear combination of $\boldsymbol \mu_{\boldsymbol \vartheta}$, the initial value $\boldsymbol \vartheta(0)$ of $\boldsymbol \vartheta$ and the initial momentum $\bm{u}(0)$.  Note that (\ref{eqn:solution for Hamilton}) is a multivariate generalisation of (\ref{eqn:over-relaxation}) with $t=\textrm{cos}^{-1}(-\kappa)$.  Equation~(\ref{eqn:solution for Hamilton}) is thus equivalent to the antithetic sampler in (\ref{eqn:antithetic}) when setting $t=\pi$ radians.  Since there is no approximation error in the Hamiltonian dynamics for a normal distribution, an MH accept-reject step is not required in the HMC sampler, and the proposed value of $\boldsymbol \vartheta$ will always be accepted.  This equivalence relation was first observed by \cite{pakman2014exact}, but was not particularly useful in their framework of sampling from a truncated multivariate normal distribution.  Our proposal for antithetic sampling is different from theirs in the sense that it is entirely deterministic, and we choose $t=\pi$ radians to induce a perfect negative proposal correlation.  \cite{pakman2014exact}, on the other hand,  suggest setting $t=\frac{\pi}{2}$ radians, which is equivalent to drawing a fresh sample from a random number generator when it is applied to the setting of a normal distribution.  We refer to this approach as the independent sampler hereafter.

So far, our discussion has mainly focused on normal $\pi(\boldsymbol \vartheta)$.  This is because an analytic solution to the Hamiltonian system is only available for a normal distribution.
It is possible to extend the proposed antithetic sampler to more general distributions by obtaining an approximation of $\boldsymbol \mu_{\boldsymbol \vartheta}$ in order to propose a new value of $\boldsymbol \vartheta$, and then accept or reject the proposal in an MH algorithm to target the true $\pi(\boldsymbol \vartheta)$, as suggested in \cite{green1992metropolis}.  However, the application of this generalisation and its variants (e.g.~\cite{creutz1987overrelaxation}) is somewhat limited due to high rejection rates in the accept-reject step \citep{neal1998suppressing}. In this case, the HMC algorithm provides a way to overcome this shortcoming.

\section{Simulation studies} \label{sec:simulation}

We now study the efficiency of the antithetic variable technique described in Section~\ref{sec:deterministic}.  Two examples are presented.  The first examines the antithetic sampler in a more general setting, while the second is specific to the application in Section~\ref{sec:application}.  Reported IACT values of the parameters are computed using the \verb|coda| package \citep{plummer2006coda} in \verb|R|.
\vspace{\baselineskip}

\noindent
{\bf Example 1. }The stationary distribution $\pi(\boldsymbol \theta)$ is specified as a bivariate normal distribution with high correlation (0.99) between the variables.  We investigate the performance of three sampling schemes - the independent sampler, the over-relaxation algorithm with $\kappa=0.9$, and a coupling of the over-relaxation algorithm (on the first margin) with the antithetic sampler (on the second margin).  Note that this coupling strategy introduces stochasticity into the antithetic sampler, which is essential to produce an ergodic Markov chain.  The samplers are each run for 10\,000 iterations from the same initialised value $(2,2)$, and the update on each margin is performed conditional on the other.  Figure~\ref{fig:correlated normal} illustrates the trajectories of the first 50 samples generated.  Exploration of the target space is reduced to a random walk under the independent sampler.  In contrast, the other two samplers move between different contours of the density and explore the full support of the distribution in an elliptical manner, thereby reducing the IACT significantly.  The IACT decreases further when the over-relaxation algorithm on the second margin is replaced by antithetic sampling.  In this analysis, the mixing of both margins is improved by a factor of 1.75.
\begin{figure}[h!]
    \centering
    \includegraphics[trim= 0.5cm 1cm 6.5cm 9.5cm,clip,width=0.6\textwidth]{normal}
    \caption{Trajectories of the first 50 samples generated from the independent sampler (left), the over-relaxation algorithm with $\kappa=0.9$ (middle), and the over-relaxation algorithm coupled with the antithetic sampler (right).  The blue solid lines represent the 95\% confidence region of the bivariate normal distribution.}
    \label{fig:correlated normal}
\end{figure}
\vspace{\baselineskip}

\noindent
{\bf Example 2. }A simulated dataset is generated following the MVP model given in (\ref{eqn:MVP})--(\ref{eqn:latent}), with $D=8, P=162, T=16$ and values of the parameters $\boldsymbol \theta=(\boldsymbol \beta, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \Sigma_{\boldsymbol \alpha})$ set to be the posterior mean estimates of the parameters in Model 1 of the female contraceptive product analysis of Section~\ref{sec:application}.  To avoid hand-tuning the stepsize $\varepsilon$ and the trajectory length $\mathcal{T}$ for the HMC update of $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$, we utilise the No-U-Turn Sampler (NUTS) with the dual averaging scheme of \cite{hoffman2014no}.  We use the following non-informative prior distributions: $\boldsymbol \beta \sim \mathcal{N}({\bf 0}, 100\bm{I})$, $\boldsymbol \Sigma_{\boldsymbol \alpha} \sim \mathcal{IW}(9, \bm{I})$ and the prior distribution on the lower triangular Cholesky factor $\bm{\mathcal{L}}_{\boldsymbol \epsilon}$ given in (\ref{eqn:prior cholesky}).  The sampling scheme is run for 30\,000 iterations, with the first 5\,000 samples discarded as burn-in.  Appendix~\ref{app:sampling scheme} details the Gibbs sampling scheme.

Figure~\ref{fig:autocorrelation} compares graphically the marginal posterior densities and sample autocorrelations of randomly sampled random effects $\boldsymbol \alpha_{1:P}$ and the regression parameter $\boldsymbol \beta$ between independent and antithetic sampling.  Despite the absence of a stochastic component in the updates of $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$, the kernel density plots of these parameters indicate that the coupling of a stochastic MCMC scheme for the remaining parameters with the antithetic variable technique gives the same posterior distributions as those under independent sampling.  The autocorrelation plots show that the samples generated from antithetic sampling have positive dependence with a higher rate of decay over the number of lags, thereby demonstrating the superior mixing of the Markov chain.  The IACT values of the randomly sampled parameters are significantly lower, with improvement factors of 3.72 and 2.10 observed for $\alpha_{3,80}$ and $\beta_{182}$ respectively.  The box plot showing the distribution of the IACT values of $\bm{\alpha}_{1:P}$ also indicates that some of these parameters are super-efficient.  Furthermore, the log IACT ratios of the independent sampler compared to the antithetic sampler are well above 0, suggesting that all $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$ parameters experience efficiency gains.  Although perfect negative correlation is induced between successive samples by the deterministic proposal, this does not necessarily translate to an equivalent autocorrelation in the posterior samples.  Rather, the negative relationship is used to reduce the magnitude of positive autocorrelation present in the MCMC samples.  Note that convergence to the posterior distribution might be slow for poorly initialised values under antithetic sampling so we suggest using independent sampling during the burn-in period and later switching to the deterministic proposal.
\begin{figure}[t!]
\onehalfspacing
\centering
\includegraphics[trim= 0.2cm 0.75cm 1.25cm 0.75cm,clip,width=0.7\textwidth]{IS_vs_AS_ratio}
\caption{Marginal posterior densities of a randomly selected random effects term (top panel) and regression coefficient (bottom panel), and their sample autocorrelation plots under independent sampling (IS) and antithetic sampling (AS).  Rightmost column gives the distributions of the IACT values and the element-wise IACT ratios of IS to AS for $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$ on the log scale.}
\label{fig:autocorrelation}
\end{figure}

The remaining simulation experiments investigate the performance of the MVP model in the context of recovering the true parameters of the data generating process under different specifications of prior distribution on $\boldsymbol \theta$.  We use the posterior root-mean-square error (RMSE) defined by
\begin{equation}
\textrm{RMSE}(\boldsymbol \theta) = \sqrt{\frac{1}{N}\sum_{j=1}^N (\boldsymbol \theta^{[j]}-\boldsymbol \theta_{\textrm{true}})^2},
\label{eqn:RMSE}
\end{equation}
as the performance measure, where $\boldsymbol \theta^{[j]}$ is the $j$-th iterate from the $N$ posterior samples and $\boldsymbol \theta_{\textrm{true}}$ is the true value of $\boldsymbol \theta$.
The measure in (\ref{eqn:RMSE}) is defined for univariate $\boldsymbol \theta$.  For a multivariate $\boldsymbol \theta$, the posterior RMSE is calculated for each margin of $\boldsymbol \theta$.  All the results shown are based on 1\,000 different replicate sets of simulated data with the same true parameter values.

We first consider the conditionally conjugate hierarchical inverse-Wishart $\mathcal{HIW}(\lambda, \bm{A})$ prior of \cite{huang2013simple} with degrees of freedom $\lambda$ and positive scale parameter $\bm{A}=(A_1, \dotsc, A_D)^\top$ as an alternative to the inverse-Wishart prior on the $D \times D$ covariance matrix $\boldsymbol \Sigma_{\boldsymbol \alpha}$,
\begin{equation*}
\begin{gathered}
\boldsymbol \Sigma_{\boldsymbol \alpha}|a_1,\dotsc,a_D \sim \mathcal{IW}\Bigg(\lambda + D - 1, 2 \lambda\textrm{diag}\bigg(\frac{1}{a_1}, \dotsc, \frac{1}{a_D}\bigg)\Bigg), \\
a_i \stackrel{iid}{\sim} \mathcal{IG}(0.5, A^{-2}_i), \quad i=1,\dotsc, D,
\end{gathered}
\end{equation*}
where $\mathcal{IG}(a, b)$ is an inverse-Gamma distribution with shape $a$ and scale $b$.  The marginal prior of the standard deviation in $\boldsymbol \Sigma_{\boldsymbol \alpha}$ is a half-$t(\lambda, A_i)$ distribution, as suggested in \cite{gelman2006prior}.  In the simulation, we select $\lambda=2$ and choose a weakly informative scale parameter whereby $A_1 = A_2 = 0.23$ and $A_3 = \cdots = A_8 = 0.46$ so that approximately 95\% of the half-$t$ density is below 1 and 2 respectively.  This specification is relevant to the real data application in Section~\ref{sec:application}, where our prior belief is that the variability in the tendency of GPs to discuss pill contraceptives is lower compared to non-pill alternatives.  In contrast, the inverse-Wishart prior assumes the same variability for all variance parameters $\sigma^2_{\alpha_i}$ in $\boldsymbol \Sigma_{\boldsymbol \alpha}$.  Figure~\ref{fig:IW vs HIW} shows the distribution of the average RMSE ratio of each type of parameter in $\boldsymbol \Sigma_{\boldsymbol \alpha}$, based on 1\,000 replicate simulations, for the hierarchical inverse-Wishart prior versus the inverse-Wishart prior. Although the hierarchical inverse-Wishart prior is flexible enough to specify different strengths of prior on each $\sigma^2_{\alpha_i}$, Figure~\ref{fig:IW vs HIW} shows that in this case its performance is similar to the more restrictive inverse-Wishart prior.  This result is somewhat unsurprising considering that the estimated $\sigma^2_{\alpha_i}$ in the application example are more or less similar across the different contraceptive products (see Appendix~\ref{app:covariance}).  The distributions for the posterior RMSE ratio of the correlation coefficients and the partial correlations are concentrated around 1 since both the hierarchical inverse-Wishart prior with $\lambda=2$ and the inverse-Wishart prior with $D+1$ degrees of freedom and scale matrix $\bm{I}$ induce the same marginally uniform prior, i.e.~(\ref{eqn:marginally uniform prior}) with $\nu = D+1$, on the resulting correlation matrix $\bm{R}_{\boldsymbol \alpha}$, which in turn gives the same implied LKJ distribution on the partial correlations.
\begin{figure}[t!]
\centering
\begin{subfigure}{.475\textwidth}
  \centering
  \includegraphics[trim= 2.25cm 1.75cm 12.25cm 11.5cm,clip,width=0.8\textwidth]{hiw_vs_iw}
  \caption{Hierarchical inverse-Wishart prior versus inverse-Wishart prior on $\boldsymbol \Sigma_{\boldsymbol \alpha}$.}
\label{fig:IW vs HIW}
\end{subfigure}
\hfill
\begin{subfigure}{.475\textwidth}
  \centering
  \includegraphics[trim= 2.25cm 1.75cm 12.25cm 11.5cm,clip,width=0.8\textwidth]{hs_vs_n}
  \caption{Horseshoe shrinkage prior versus normal prior on $\boldsymbol \beta$.}
\label{fig:horseshoe vs normal}
\end{subfigure}
\caption{Distributions of the average posterior RMSE ratio of all parameters in (a) $\boldsymbol \Sigma_{\boldsymbol \alpha}$ or (b) $\boldsymbol \beta$, based on 1\,000 replicate analyses, under different prior choices. (a) Standard deviations, correlations and partial correlations for parameters in $\boldsymbol \Sigma_{\boldsymbol \alpha}$ for the hierarchical inverse-Wishart prior versus the inverse-Wishart prior on $\boldsymbol \Sigma_{\boldsymbol \alpha}$. (b) Sparse regression coefficients $\beta_i=0$ and non-sparse coefficients $\beta_i\neq 0$ for the horseshoe prior versus the $\mathcal{N}(\boldsymbol{0},100\boldsymbol{I})$ prior on $\boldsymbol \beta$.
}
\label{fig:prior beta and sigma}
\end{figure}

To identify sparse signals (coefficients which are significant) in the regression parameter $\boldsymbol \beta$, we employ the horseshoe shrinkage prior \citep{carvalho2010horseshoe} given by
\begin{equation*}
\beta_i|\lambda_i, \tau \sim \mathcal{N}(0, \tau^2 \lambda_i^2), \quad \lambda_i \sim \mathcal{C}^+(0,1), \quad \tau \sim \mathcal{C}^+(0,1),
\end{equation*}
where $\mathcal{C}^+(0,1)$ is a half-Cauchy distribution with location 0 and scale 1 restricted to positive support.  The simulation is carried out by setting 75\% of the smallest non-intercept regression coefficients (in absolute value) in $\boldsymbol \beta$ to 0, from which we generate the simulated datasets.  We model the prior on each intercept separately by a flat $\mathcal{N}(0, 100)$ distribution to avoid heavily penalising these parameters.  Gibbs sampling from the posterior distribution of $\boldsymbol \beta$ is implemented by adopting the latent variable formulation in \cite{makalic2016simple}.  Figure~\ref{fig:horseshoe vs normal} displays the results of comparing this prior specification for $\boldsymbol \beta$ to a $\mathcal{N}({\bf 0}, 100\bm{I})$ prior, again in terms of the average RMSE ratio over all regression parameters.  The horseshoe prior performs as well as the $\mathcal{N}({\bf 0}, 100\bm{I})$ prior on non-zero entries of $\boldsymbol \beta$, although the variability in the RMSE ratio is large.  On the other hand, the horseshoe prior outperforms the normal prior for those parameters whose true values are zero, reducing the RMSE by half.  This occurs as the horseshoe prior places a greater density around zero, which results in a more concentrated posterior distribution for parameters which are truly zero.  Therefore, it is an attractive default option when we expect sparsity in the regression parameters, as is the case for our analysis of the characteristics affecting the decision-making behaviour of GPs in the next section.

\section{Discussion of female contraceptive products by Australian GPs} \label{sec:application}
\subsection{Background and aims of study} \label{subsec:background}

In order to study the decision-making behaviour of Australian GPs, we obtain data from \cite{fiebig2017consideration} who design a stated preference experiment in which GPs are asked to select the contraceptive products that they would consider discussing with hypothetical female patients.  The GPs evaluate a sequence of vignettes where patients are defined in terms of socio-economic and clinical characteristics that are varied as part of the experimental design.  Table~\ref{table:attributes} in Appendix~\ref{app:attributes} contains the attributes of the patients with a description for each level of the categorical variables.  The GPs choose from a set of 9 products that they would discuss with the patient before deciding upon their most preferred product to be subsequently prescribed to the patient.  A sample of 162 GPs participated in the experiment where each subject makes choices for 16 different patients, resulting in 2\,592 observations.  The following covariate information is collected on the GPs themselves: age, gender, whether they are registered as a Fellow of the Royal Australian College of GPs, whether they have a certificate in family planning, whether they are an Australian medical graduate, whether their location of practice is in an urban area and whether they bulk-bill patients.  Analysis of this panel data is based on the set of binary outcomes as to whether or not to discuss each of the contraceptive products.  Due to low occurrences for the prescription of the hormonal patch which was yet to be released in the Australian market, we removed this product from the dataset leaving observations on the 8 remaining products.

The experiment is designed to mimic the choice problem faced by GPs in a consultation where they need to match a product with a particular patient.  In characterising such a decision problem, \cite{frank2007custom} distinguish between ``custom-made'' and ``ready-to-wear'' (or norm-based) choices.  A custom-made choice involves the GP undertaking a careful evaluation of the patient and then matching her to an appropriate product.  However, as new products are introduced, GPs face considerable costs in the process of gaining the knowledge and expertise required to discuss and prescribe these products.  This is particularly the case when more familiar products are available even though they may be somewhat inferior to the new products; an especially salient situation in the market for contraceptive products.  In such cases, some GPs will tend to adopt norms (here particular products) that work well for a broad class of patients and to place less weight on certain patient attributes that would indicate a different product that is potentially a better match.

Particular interest is in the dependence between the products.  That is, which products tend to be discussed together and which tend to form distinct clusters.  If GPs pursue custom-made strategies, then a considerable portion of the dependence between products will be explained by the attributes of the patient.  Conditional on the observable features of the patient and characteristics of the GPs, remaining dependencies will reflect the relationship between unobservables related to evaluations of the suitability of certain products for a particular patient, and how individual GP's product effects are correlated across products.  The proposed model is designed to capture these forms of heterogeneity and will permit a detailed analysis of the choices.

The prevalence of ready-to-wear choices is one possible explanation for the relatively low uptake of long acting reversible contraceptive (LARC) methods in Australia \citep{black2013australian}.  LARC methods are contraceptives that are administered less frequently than monthly and include hormonal implants, intrauterine contraception (IUC), both hormonal and copper-bearing, and contraceptive injections.  There is increasing support for the greater use of these more effective methods to reduce unintended pregnancies and abortion rates.  In our analysis below, we will use the model to explore a case where there is no clinical reason why at least one of these LARC methods should not be considered for discussion by GPs.  For ease of presentation, we will use the subscripts in Table~\ref{table:products} to denote the products.
\begin{table}[ht!]
\centering
\onehalfspacing
\begin{tabular}{|c|l|}
\hline
Subscript & \multicolumn{1}{c|}{Product} \\ \hline
1 & Combined pill \\
2 & Mini-pill \\
\cellcolor{gray!25}3 & \cellcolor{gray!25}Hormonal injection \\
\cellcolor{gray!25}4 & \cellcolor{gray!25}Hormonal implant \\
\cellcolor{gray!25}5 & \cellcolor{gray!25}Hormonal IUD \\
6 & Vaginal ring \\
\cellcolor{gray!25}7 & \cellcolor{gray!25}Copper IUD \\
8 & Condom \\ \hline
\end{tabular}
\caption{Correspondence of parameter subscripts to each female contraceptive product.  Long acting reversible contraceptive methods are shown in grey.}
\label{table:products}
\end{table}

\subsection{Analysis and results} \label{subsec:discussion}

We consider two different models for the data:
\begin{align}
\text{Model 1: }& \bm{y}^\ast_{it} = \boldsymbol \alpha_i + \bm{B} \bm{x}_{it} + \boldsymbol \epsilon_{it}, \label{eqn:model 1}\\
\text{Model 2: }& \bm{y}^\ast_{it} = \boldsymbol \alpha_i + \bm{B} \bm{x}_{it} + \bm{C}  \bm{z}_i + \boldsymbol \epsilon_{it}, \label{eqn:model 2}
\end{align}
for $i=1,\dotsc,P=162$ GPs and $t=1,\dotsc,T=16$ patients.  Here $\boldsymbol \alpha_i$ and $\bm{C}  \bm{z}_i$ respectively represent GP-specific random and fixed effects with $\bm{z}_i$ being a vector of GP characteristics, and $\bm{B} \bm{x}_{it}$ represents fixed effects of the patient.  We select a horseshoe prior on $\boldsymbol \beta=\textrm{vec}(\bm{B})$ and model the covariance matrix $\boldsymbol \Sigma_{\boldsymbol \alpha}$ of the random effects by the $\mathcal{HIW}(2,\bm{A})$ prior in Section~\ref{sec:simulation} where $\bm{A}=(0.23, 0.23, 0.46, \dotsc, 0.46)^\top$.  The scale is chosen to express the prior information that the variances of the random effects are expected to be small, with those for the pill products being less variable compared to the non-pill alternatives.  The difference between these two models is the presence of the GP-specific fixed effects in Model 2, which explain some of the relationships in the random effects of Model 1.  Let $\bm{X}=(X_1, \dotsc, X_D)^\top$ be a vector of normal random variables with covariance matrix given by $\boldsymbol \Sigma_{\bm{X}}$.  Recall that $X_i$ and $X_j$ are conditionally independent given the other random variables if the $(i,j)$-th entry of the precision matrix $\boldsymbol \Sigma^{-1}_{\bm{X}}$ is zero.

Figures~\ref{fig:correlation graph} and \ref{fig:random effect graph} give graphical summaries of the posterior distribution of the dependence structures of the latent variable $\bm{y}^\ast_{it}$ conditional on $\boldsymbol \alpha_i$ and $\bm{x}_{it}$ (as well as $\bm{z}_i$ for Model 2), and the random effects $\boldsymbol \alpha_i$ respectively.  All graphs are obtained by computing the 95\% credible interval of the posterior distribution for each entry of $\bm{R}_{\boldsymbol \epsilon}^{-1}$ and $\boldsymbol \Sigma^{-1}_{\boldsymbol \alpha}$, where an edge is formed between two nodes if the credible interval does not include 0. The absence of an edge between any two nodes indicates a potential conditional independence between the two variables given the rest.  The dependence structures associated with the latent variables are the same for both models.  This supports the use of the MVP model in order to capture the complex dependencies between different products that would otherwise be ignored in separate univariate analyses on each product.
\begin{figure}[ht!]
\centering
\begin{tikzpicture}[>=stealth',shorten >=1pt,auto,thick,main node/.style={circle,draw}]
	\draw (67.5:2cm)   node[main node, scale=0.9] (1) {$y^\ast_1$};
    \draw (22.5:2cm)   node[main node, scale=0.9] (2) {$y^\ast_2$};
    \draw (337.5:2cm)  node[main node, scale=0.9] (3) {$y^\ast_3$};
    \draw (292.5:2cm)  node[main node, scale=0.9] (4) {$y^\ast_4$};
    \draw (247.5:2cm)  node[main node, scale=0.9] (5) {$y^\ast_5$};
    \draw (202.5:2cm)  node[main node, scale=0.9] (6) {$y^\ast_6$};
    \draw (157.5:2cm)  node[main node, scale=0.9] (7) {$y^\ast_7$};
    \draw (112.5:2cm)  node[main node, scale=0.9] (8) {$y^\ast_8$};
    \path[every node/.style={font=\sffamily\small}]
    (1) edge[red, line width=0.08cm]  node [left] {} (5)
        edge[blue, line width=0.16cm] node [left] {} (6)
    (2) edge[blue, line width=0.055cm] node [left] {} (3)
        edge[red, line width=0.091cm]  node [left] {} (6)
    (3) edge[blue, line width=0.055cm] node [left] {} (2)
    	edge[blue, line width=0.2cm] node [left] {} (4)
    (4) edge[blue, line width=0.2cm] node [left] {} (3)
    	edge[blue, line width=0.082cm] node [left] {} (5)
    (5) edge[red, line width=0.08cm]  node [left] {} (1)
    	edge[blue, line width=0.082cm] node [left] {} (4)
		edge[blue, line width=0.176cm] node [left] {} (7)
        edge[red, line width=0.026cm]  node [left] {} (8)
    (6) edge[blue, line width=0.16cm] node [left] {} (1)
    	edge[red, line width=0.091cm]  node [left] {} (2)
    	edge[blue, line width=0.062cm] node [left] {} (8)
    (7) edge[blue, line width=0.176cm] node [left] {} (5)
    	edge[blue, line width=0.07cm] node [left] {} (8)
    (8) edge[red, line width=0.026cm]  node [left] {} (5)
    	edge[blue, line width=0.062cm] node [left] {} (6)
        edge[blue, line width=0.07cm] node [left] {} (7);
\end{tikzpicture}
\caption{Graphical model illustrating substantial dependence structure of the latent variables $\bm{y}^\ast$ conditional on the random effects and the covariates in both Model 1 and 2.  Edges between $y^\ast_i$ and $y^\ast_j$ are included if the 95\% credible interval of the marginal posterior distribution of the $(i,j)$-th entry of $\bm{R}_{\boldsymbol \epsilon}^{-1}$ does not contain 0.  Blue edges represent positive dependence while red edges represent negative dependence.  The thickness of the edges is proportional to the strength of the dependence.}
\label{fig:correlation graph}
\vspace{10pt}
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[>=stealth',shorten >=1pt,auto,thick,main node/.style={circle,draw}]
	\draw (67.5:2cm)  node[main node, scale=0.9] (1) {$\alpha_1$};
    \draw (22.5:2cm)  node[main node, scale=0.9] (2) {$\alpha_2$};
    \draw (337.5:2cm) node[main node, scale=0.9] (3) {$\alpha_3$};
    \draw (292.5:2cm) node[main node, scale=0.9] (4) {$\alpha_4$};
    \draw (247.5:2cm) node[main node, scale=0.9] (5) {$\alpha_5$};
    \draw (202.5:2cm) node[main node, scale=0.9] (6) {$\alpha_6$};
    \draw (157.5:2cm) node[main node, scale=0.9] (7) {$\alpha_7$};
    \draw (112.5:2cm) node[main node, scale=0.9] (8) {$\alpha_8$};
    \path[every node/.style={font=\sffamily\small}]
    (1) edge[blue, line width=0.17cm] node [left] {} (2)
    	edge[blue, line width=0.114cm] node [left] {} (3)
    (2) edge[blue, line width=0.17cm] node [left] {} (1)
    (3) edge[blue, line width=0.114cm] node [left] {} (1)
    	edge[blue, line width=0.074cm] node [left] {} (7)
    (4) edge[blue, line width=0.141cm] node [left] {} (5)
    	edge[red, line width=0.063cm]  node [left] {} (7)
    (5) edge[blue, line width=0.141cm] node [left] {} (4)
    (6) edge[blue, line width=0.098cm] node [left] {} (8)
    (7) edge[blue, line width=0.074cm] node [left] {} (3)
    	edge[red, line width=0.063cm]  node [left] {} (4)
    	edge[blue, line width=0.068cm] node [left] {} (8)
    (8) edge[blue, line width=0.098cm] node [left] {} (6)
    	edge[blue, line width=0.068cm] node [left] {} (7);
\end{tikzpicture}
\centering
\caption{Model 1.}
\label{fig:model 1 graph}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[>=stealth',shorten >=1pt,auto,thick,main node/.style={circle,draw}]
	\draw (67.5:2cm)  node[main node, scale=0.9] (1) {$\alpha_1$};
    \draw (22.5:2cm)  node[main node, scale=0.9] (2) {$\alpha_2$};
    \draw (337.5:2cm) node[main node, scale=0.9] (3) {$\alpha_3$};
    \draw (292.5:2cm) node[main node, scale=0.9] (4) {$\alpha_4$};
    \draw (247.5:2cm) node[main node, scale=0.9] (5) {$\alpha_5$};
    \draw (202.5:2cm) node[main node, scale=0.9] (6) {$\alpha_6$};
    \draw (157.5:2cm) node[main node, scale=0.9] (7) {$\alpha_7$};
    \draw (112.5:2cm) node[main node, scale=0.9] (8) {$\alpha_8$};
    \path[every node/.style={font=\sffamily\small}]
    (1) edge[blue, line width=0.165cm] node [left] {} (2)
    (2) edge[blue, line width=0.165cm] node [left] {} (1)
    (4) edge[blue, line width=0.134cm] node [left] {} (5)
    	edge[red, line width=0.042cm]  node [left] {} (7)
    (5) edge[blue, line width=0.134cm] node [left] {} (4)
    (6) edge[blue, line width=0.119cm] node [left] {} (8)
    (7) edge[red, line width=0.042cm]  node [left] {} (4)
    	edge[blue, line width=0.086cm] node [left] {} (8)
    (8) edge[blue, line width=0.119cm] node [left] {} (6)
    	edge[blue, line width=0.086cm] node [left] {} (7);
\end{tikzpicture}
\centering
\caption{Model 2.}
\label{fig:model 2 graph}
\end{subfigure}
\caption{Graphical models illustrating substantial dependence structure of the GP-specific random effects $\boldsymbol \alpha$ in each model.  Edges between $\alpha_i$ and $\alpha_j$ are included if the 95\% credible interval of the marginal posterior distribution of the $(i,j)$-th entry of $\boldsymbol \Sigma_{\boldsymbol \alpha}^{-1}$ does not contain 0.  Blue edges represent positive dependence while red edges represent negative dependence.  The thickness of the edges is proportional to the strength of the dependence.}
\label{fig:random effect graph}
\end{figure}

Figure~\ref{fig:correlation graph} is also instrumental in explaining the suitability of the contraceptive products for a patient in terms of substitute goods, which in consumer theory is defined as products with similar functions that can be used in place of each other.  For conciseness, we only focus on some important relationships illustrated in the graphical model.  The propensity to discuss pill products $(y^\ast_1, y^\ast_2)$ are independent of each other given the hormonal IUD and the vaginal ring $(y^\ast_5, y^\ast_6)$ by the Markov property since all paths from $y^\ast_1$ to $y^\ast_2$ pass through $(y^\ast_5, y^\ast_6)$, reflecting the use of these non-pill contraceptives as pill alternatives dictated by particular clinical conditions.  The clique formed between $(y^\ast_5, y^\ast_7, y^\ast_8)$ suggests dependence in the propensity to discuss the hormonal IUD, the copper IUD and the condoms.  In fact, the posterior correlation between the propensity scores for both the IUD methods $(y^\ast_5, y^\ast_7)$ is around 0.52 on average (see Appendix~\ref{app:correlation}), suggesting a high tendency of these products to be discussed together.  This also reflects the fact that these IUD methods are substitutes.  Noticeably, the propensity to discuss the hormonal injection and the hormonal implant $(y^\ast_4, y^\ast_5)$ exhibit the highest level of association as indicated by our model, with a mean posterior correlation of 0.59.  This indicates the likelihood of these two prominent LARC products being included together in discussions, and it is consistent with them being close substitutes for each other for many patients.

Figure~\ref{fig:random effect graph} can be interpreted in the same way as Figure~\ref{fig:correlation graph}, regarding the substitutability of different products but in the context of ready-to-wear choices.  This is because the random effects in (\ref{eqn:model 1}) characterise the persistence of GPs in discussing a particular product after observing the patient's attributes.  There are clear differences in the graphical structure when comparing Figures~\ref{fig:model 1 graph} and \ref{fig:model 2 graph}.  The changes in the dependence structure of the GP random effects arise because some of the persistence in product choices can be explained by GP characteristics.  For example, the tendency of GPs to include both the hormonal injection and the copper IUD $(\alpha_3, \alpha_7)$ as ready-to-wear choices is due to their age (see significance of GP characteristics in Appendix~\ref{app:estimation result patient}).  The posterior structure also provides some confidence that the random effects specification is useful in capturing important GP characteristics that are not directly observed.  Three clusters of products with substantial dependence in ready-to-wear choices are identified from the model after accounting for the observed GP characteristics.  Particularly relevant is the dependence between the hormonal IUD and the implant $(\alpha_4, \alpha_5)$.  There is positive correlation between these two LARCs, indicating the tendency for GP attitudes (either positive or negative) to be aligned.  A second cluster includes both of the pills $(\alpha_1, \alpha_2)$ which is consistent with these products being used as a ready-to-wear default.  GPs who are more likely to discuss the combined pill after conditioning on the patient's attributes behave similarly when considering the mini-pill.  Contraceptives that are not pill- or hormone-based form the final bundle.

Our models allow us to examine posterior predictions for a range of patients.  Since we are interested in the uptake of LARC products, we specify a particular female patient where there is no clinical reason why a LARC should not be considered for discussion.  Table~\ref{table:attributes} of Appendix~\ref{app:attributes} gives the attributes of this base-case patient.  Figure~\ref{fig:marginal probability} summarises the estimate of the predictive probability of a GP discussing a particular product, where the range of predictions shown is generated for all GPs in the sample based on Model 2.  For this particular base-case patient, there is considerable agreement amongst all GPs in the sample that the combined pill (product 1) is one of the most suitable products to be discussed, but they have much more variable views on the other products.  Amongst the LARCs (products 3, 4, 5 and 7), the hormonal injection (product 3) and the implant (product 4) are the products which are the most likely to be discussed, with the variability across GPs perhaps simply reflecting a view that they are good substitutes to each other, which is in fact what we find in Figure~\ref{fig:correlation graph}.  GPs could indeed have consistent views about the need to discuss LARCs, as they do with the combined pill, but they are divided on which of the LARC products to discuss.  To explore this possibility, the final column in Figure~\ref{fig:marginal probability} shows the predicted probability of the GPs discussing at least one of these two products, that is $\mathbb{P}(y_3 + y_4 \geq 1)$.  The results suggest that the GPs will discuss either product 3 or 4 (or both) with similar probability to the combined pill.  While this joint probability does indicate a median that is similar to that of discussing the combined pill, the variability across GPs remains much larger than that associated with the combined pill.  This evidence is consistent with the hypothesised resistance amongst some GPs to even discuss LARCs, let alone recommend them.
\begin{figure}[t!]
\onehalfspacing
\centering
\includegraphics[trim= 2cm 1.5cm 2cm 3.5cm,clip,width=0.6\textwidth]{marginal_probability}
\caption{Predicted probability of a GP discussing each product for a base-case patient for each of the 162 Australian GPs.}
\label{fig:marginal probability}
\end{figure}

\subsection{Comparing sampling schemes}

In order to investigate the performance of the antithetic sampler, Figure~\ref{fig:posterior application} illustrates marginal posterior distributions of those Model 2 parameters whose densities demonstrate the greatest visual differences between independent and antithetic sampling of the random effects $\boldsymbol \alpha_{1:P}$ and regression parameters $\boldsymbol \beta$.  The marginal posterior distributions of $\alpha_{7,110}$ and $\beta_{236}$ are effectively the same under both updating approaches.  This occurs because the mean of the conditional posterior distribution, which is a key ingredient in the deterministic antithetic sampler proposal, changes between iterations; a change largely driven by the stochastic update of the latent variable $\bm{y}^\ast$.  This outcome suggests that the posterior distribution of the other parameters remains adequately explored by the antithetic sampler.
\begin{figure}[t!]
    \centering
    \includegraphics[trim= 0.5cm 2.25cm 0.5cm 8cm,clip,width=0.8\textwidth]{posterior_application}
    \caption{Marginal posterior density estimates of those Model 2 parameters with the greatest visual differences  between using independent sampling (IS) and antithetic sampling (AS) for $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$.}
    \label{fig:posterior application}
\end{figure}

Table~\ref{table:IACT model 2} compares the performance between independent and antithetic sampling schemes when estimating Model 2.  The antithetic variable method generates samples marginally faster than independent sampling because it is deterministic.  Based on the results shown, we observe an improvement of 4.86 and 3.31 times performance gain on average in the mixing of $\boldsymbol \alpha_{1:P}$ and $\boldsymbol \beta$ respectively.  As a result of this, the mean IACT of $\bm{y}^\ast$ is also improved.
\begin{table}[t!]
\centering
\doublespacing
\begin{tabular}{|c|r|r|r|r|r|}
\hline
\multirow{2}{*}{Parameter} & \multicolumn{2}{c|}{Mean IACT} & \multicolumn{3}{c|}{IACT Ratio} \\ \cline{2-6}
 & \multicolumn{1}{c|}{IS} & \multicolumn{1}{c|}{AS} & \multicolumn{1}{c|}{Min} & \multicolumn{1}{c|}{Max} & \multicolumn{1}{c|}{Mean} \\ \hline
$\bm{y}^\ast$             & 3.6387  & 2.6686  & 0.8242 & 3.1419  & 1.2127 \\
$\boldsymbol \alpha_{1:P}$            & 16.8872 & 4.6456  & 1.4857 & 13.3424 & 4.8632 \\
$\boldsymbol \beta$                  & 15.0446 & 4.0105  & 1.4566 & 16.0173 & 3.3111 \\
$\textrm{vechL}(\bm{\mathcal{L}}_{\boldsymbol \epsilon})$        & 14.8292 & 14.5422 & 0.9338 & 1.1737  & 1.0191 \\
$\textrm{vechL}(\bm{R}_{\boldsymbol \epsilon})$          & 12.7311 & 12.5170 & 0.9147 & 1.1509  & 1.0180 \\
$\textrm{diag}(\boldsymbol \Sigma_{\boldsymbol \alpha})$    & 24.8056 & 14.6929 & 1.3130 & 2.0651  & 1.7222 \\
$\textrm{vechL}(\bm{R}_{\boldsymbol \alpha})$ & 9.5025  & 5.1716  & 1.4599 & 2.3336  & 1.8424 \\ \hline
Time per iteration & 0.0243 & 0.0239 & \multicolumn{1}{c|}{-} & \multicolumn{1}{c|}{-} & \multicolumn{1}{c|}{-} \\ \hline
\end{tabular}
\caption{Comparison of the performance between independent sampling (IS) and antithetic sampling (AS) in the contraceptive products preference data in terms of the speed (seconds per iteration), the mean IACT and the IACT ratio for each block of parameter.}
\label{table:IACT model 2}
\end{table}

\section{Conclusion} \label{sec:conclusion}

Many methods exist for fitting a multinomial logit model with random effects, such as simulated maximum likelihood \citep{gong2004mobility}, quadrature \citep{hartzel2001multinomial,hedeker2003mixed}, multinomial-Poisson transformation \citep{lee2017poisson}, and moment-based estimation \citep{perry2017fast}, among others.  Computational strategies for the MVP model, on the other hand, are less well studied.  In this article, we introduce a HMC sampling approach to generate the posterior samples of $\bm{R}_{\boldsymbol \epsilon}$.  This method requires reparameterising $\bm{R}_{\boldsymbol \epsilon}$ into an unconstrained Cholesky factor in order to circumvent the restrictive properties of a correlation matrix having diagonal entries of 1 and being positive definite.  Furthermore, we propose a novel antithetic variable technique to accelerate the mixing of the random effects and the regression parameters, where significant gains in efficiency are observed in our application.  Although our antithetic sampling deterministically specifies the proposal distribution within the Metropolis-Hastings update, the ergodicity of the Markov chain is unaffected when it is embedded within a larger system of stochastic updates.

Our application considers the discussion of female contraceptive products by Australian GPs based on outcomes from the second stage of the stated preference data from \cite{fiebig2017consideration}.  An examination of the correlation matrix underlying the choices reveals a complex dependence structure between the products, hence indicating the plausibility of our formulation to model these choices in a multivariate setting.  Our empirical study also suggests evidence of medical practice variation among the GPs, especially with regard to the inclusion of LARCs in the discussion with patients.  The combined pill was the most popular contraceptive choice among the patients, and it represented a likely ready-to-wear default for many GPs.  Without GPs even discussing LARCs, their uptake was likely to remain relatively constrained in such a context.

\section*{Acknowledgements}

David Gunawan, Denzil Fiebig and Robert Kohn were partially supported by the Australian Research Council Discovery Project scheme DP150104630 and Scott Sisson was partially supported by the Discovery Project grant DP160102544.  Vincent Chin, David Gunawan, Robert Kohn and Scott Sisson were also partially supported by the Australian Centre of Excellence for Mathematical \& Statistical Frontiers (ACEMS) grant CE140100049.

\bibliographystyle{chicago}
\interlinepenalty=10000

\section{Parameter expansion for the MVP model} \label{app:parameter}

We first describe the parameter expansion (PX) strategy in \cite{talhouk2012efficient} for the MVP
model based on the marginally uniform prior on $\bm{R}_{\boldsymbol \epsilon}$ of \cite{barnard2000modeling}.  We then illustrate
the extension of this method to the panel data setting where the resulting sampler fails to
converge to the correct posterior distribution.

Consider the MVP model as in Section~\ref{sec:multivariate}, but without the random effects,
\begin{gather}
\bm{y}^*_{i} = \bm{B} \bm{x}_{i} + \boldsymbol \epsilon_{i}, \\
\boldsymbol \epsilon_{i} = (\epsilon_{1,i}, \dotsc, \epsilon_{D,i})^\top \stackrel{iid}{\sim} \mathcal{N}({\bf 0}, \bm{R}_{\boldsymbol \epsilon}),
\end{gather}
for $i=1,\dotsc,n$.  The posterior density of interest given the set of observed discrete outcomes
$\bm{y}=(\bm{y}_1,\dotsc,\bm{y}_n)$ and prior distribution $p(\bm{B}, \bm{R}_{\boldsymbol \epsilon})$ is
\begin{equation}
\pi(\bm{y}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}|\bm{y}) \propto p(\bm{y}|\bm{y}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon})p(\bm{y}^*|\bm{B}, \bm{R}_{\boldsymbol \epsilon})p(\bm{B}, \bm{R}_{\boldsymbol \epsilon}).
\label{eqn:MVP posterior}
\end{equation}
By taking advantage of the non-identifiability of the regression parameter $\bm{B}$ when $\bm{R}_{\boldsymbol \epsilon}$ is
not set to be a correlation matrix, \cite{talhouk2012efficient} introduce an expansion parameter
$\bm{\mathcal{D}}=\textrm{diag}(\delta_1, \dotsc, \delta_D)$ and define a new set of latent variables $\bm{z}^*_i=\bm{\mathcal{D}}\bm{y}^*_i, i=1, \dotsc, n$ such that
the posterior density under this transformation is
\begin{equation}
\pi(\bm{z}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{\mathcal{D}}|\bm{y}) \propto p(\bm{y}|\bm{z}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{\mathcal{D}})p(\bm{z}^*|\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{\mathcal{D}})p(\bm{B}, \bm{R}_{\boldsymbol \epsilon},\bm{\mathcal{D}}),
\label{eqn:MVP PX posterior}
\end{equation}
where $\bm{z}^*=(\bm{z}^*_1, \dotsc, \bm{z}^*_n)$.  The marginal posterior distribution of $(\bm{B}, \bm{R}_{\boldsymbol \epsilon})$ in (\ref{eqn:MVP PX posterior}) is the same as
(\ref{eqn:MVP posterior}) provided that the marginal prior of $p(\bm{B}, \bm{R}_{\boldsymbol \epsilon},\bm{\mathcal{D}})$ is $p(\bm{B}, \bm{R}_{\boldsymbol \epsilon})$.  Using the prior distribution
on $\bm{R}_{\boldsymbol \epsilon}$ in (\ref{eqn:marginally uniform prior}) with $\nu=D+1$ and choosing
\begin{equation}
\delta_i^2 \stackrel{iid}{\sim} \mathcal{IG}((D+1)/2, \bm{R}_{\boldsymbol \epsilon}^{-1}(i,i)/2),
\label{eqn:inverse gamma}
\end{equation}
where $\bm{R}_{\boldsymbol \epsilon}^{-1}(i,i)$ is the $i$-th diagonal entry of $\bm{R}_{\boldsymbol \epsilon}$, the covariance matrix $\boldsymbol \Sigma_{\boldsymbol \epsilon}=\bm{\mathcal{D}} \bm{R}_{\boldsymbol \epsilon} \bm{\mathcal{D}}$ is shown
to have an $\mathcal{IW}(D+1, \bm{I})$ distribution.  With this reparameterisation of the constrained $\bm{R}_{\boldsymbol \epsilon}$
into an unconstrained $\boldsymbol \Sigma_{\boldsymbol \epsilon}$, posterior sampling from the conditional posterior distribution
$\pi(\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{\mathcal{D}}|\bm{z}^*, \bm{y})$ can be circumvented by a change of variables technique.  That is, by letting
$\boldsymbol \Sigma_{\boldsymbol \epsilon}=\bm{\mathcal{D}} \bm{R}_{\boldsymbol \epsilon} \bm{\mathcal{D}}$ and $\bm{\gamma}=\bm{\mathcal{D}} \bm{B}$, iterates from the posterior distribution in (\ref{eqn:MVP posterior}) can be obtained
by the following steps:
\begin{enumerate}
\item Sample $\bm{y}^*$ from a truncated multivariate normal distribution (Step 1 in Appendix~\ref{app:sampling scheme}).
\item Sample $\bm{\mathcal{D}}$ following (\ref{eqn:inverse gamma}) and set $\bm{z}_i^* = \bm{\mathcal{D}} \bm{y}_i^*, i=1,\dotsc,n$.
\item Sample $(\bm{\gamma}, \boldsymbol \Sigma_{\boldsymbol \epsilon}) \sim \pi(\bm{\gamma}, \boldsymbol \Sigma_{\boldsymbol \epsilon}|\bm{z}^*, \bm{y})$.
\item Set $\bm{B}=\bm{\mathcal{D}}^* \bm{\gamma}, \bm{R}_{\boldsymbol \epsilon}=\bm{\mathcal{D}}^* \boldsymbol \Sigma_{\boldsymbol \epsilon} \bm{\mathcal{D}}^*$ and $\bm{y}_i^*=\bm{\mathcal{D}}^* \bm{z}_i^*, i=1, \dotsc, n$ where $\bm{\mathcal{D}}^* = \textrm{diag}(\boldsymbol \Sigma_{\boldsymbol \epsilon})^{-1/2}$.
\end{enumerate}

The posterior density of the MVP model with random effects in (\ref{eqn:MVP})--(\ref{eqn:latent}), on the other
hand, is
\begin{equation}
\begin{aligned}
\pi(\bm{y}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \alpha_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha}|\bm{y}) \propto & p(\bm{y}|\bm{y}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \alpha_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha}) \times \\
& p(\bm{y}^*|\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \alpha_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha})p(\boldsymbol \alpha_{1:P}|\boldsymbol \Sigma_{\boldsymbol \alpha})p(\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \Sigma_{\boldsymbol \alpha}).
\end{aligned}
\label{eqn:MVP with RE posterior}
\end{equation}
Following \cite{talhouk2012efficient}, we define new sets of latents variables $\bm{z}^*_{it} = \bm{\mathcal{D}} \bm{y}^*_{it}$ and
$\bm{w}_i = \bm{\mathcal{D}} \boldsymbol \alpha_i$ for $i=1,\dotsc,P$ and $t=1,\dotsc,T$, which gives the new posterior density as
\begin{equation}
\begin{aligned}
\pi(\bm{z}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{w}_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha}, \bm{\mathcal{D}}|\bm{y}) \propto & p(\bm{y}|\bm{z}^*, \bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{w}_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha}, \bm{\mathcal{D}}) p(\bm{z}^*|\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \bm{w}_{1:P}, \boldsymbol \Sigma_{\boldsymbol \alpha}, \bm{\mathcal{D}}) \times \\ &p(\bm{w}_{1:P}|\boldsymbol \Sigma_{\boldsymbol \alpha}, \bm{\mathcal{D}})p(\bm{B}, \bm{R}_{\boldsymbol \epsilon}, \boldsymbol \Sigma_{\boldsymbol \alpha}, \bm{\mathcal{D}}),
\end{aligned}
\label{eqn:MVP with RE PX posterior}
\end{equation}
where $\bm{z}^*=(\bm{z}_{it};i=1, \dotsc,P, t=1, \dotsc, T)$ and $\bm{w}_{1:P}=(\bm{w}_1, \dotsc, \bm{w}_P)$.  Assuming a conjugate
prior on $\boldsymbol \Sigma_{\boldsymbol \alpha}$, we then generate the samples from the posterior distribution in (\ref{eqn:MVP with RE PX posterior}) by a change
of variables ($\boldsymbol \Sigma_{\boldsymbol \epsilon}=\bm{\mathcal{D}} \bm{R}_{\boldsymbol \epsilon} \bm{\mathcal{D}}$ and $\bm{\gamma}=\bm{\mathcal{D}} \bm{B}$) according to
\begin{enumerate}
\item Sample $\bm{y}^*$ from a truncated multivariate normal distribution (Step 1 in Appendix~\ref{app:sampling scheme}).
\item Sample $\boldsymbol \alpha_{1:P}$ from a multivariate normal distribution (Step 4 in Appendix~\ref{app:sampling scheme}).
\item Sample $\boldsymbol \Sigma_{\boldsymbol \alpha}$ from an inverse-Wishart distribution (Step 6 in Appendix~\ref{app:sampling scheme}).
\item Sample $\bm{\mathcal{D}}$ following (\ref{eqn:inverse gamma}) and set $\bm{z}_{it}^* = \bm{\mathcal{D}} \bm{y}_{it}^*$ and $\bm{w}_i = \bm{\mathcal{D}} \boldsymbol \alpha_i$ for $i=1,\dotsc,P$ and
$t=1,\dotsc,T$.
\item Sample $(\bm{\gamma}, \boldsymbol \Sigma_{\boldsymbol \epsilon}) \sim \pi(\bm{\gamma}, \boldsymbol \Sigma_{\boldsymbol \epsilon}|\bm{z}^*, \bm{y})$.
\item Set $\bm{B}=\bm{\mathcal{D}}^* \bm{\gamma}, \bm{R}_{\boldsymbol \epsilon}=\bm{\mathcal{D}}^* \boldsymbol \Sigma_{\boldsymbol \epsilon} \bm{\mathcal{D}}^*$ and $\bm{y}_{it}^*=\bm{\mathcal{D}}^* \bm{z}_{it}^*, i=1, \dotsc, P, t=1,\dotsc,T$ where $\bm{\mathcal{D}}^* = \textrm{diag}(\boldsymbol \Sigma_{\boldsymbol \epsilon})^{-1/2}$.
\end{enumerate}
We choose the conjugate matrix normal prior on $\bm{B}$, the marginally uniform prior in (\ref{eqn:marginally uniform prior}) on
$\bm{R}_{\boldsymbol \epsilon}$ and the hierarchical inverse-Wishart prior on $\boldsymbol \Sigma_{\boldsymbol \alpha}$.  Kernel density plots comparing the
posterior distribution of the variance parameter in $\boldsymbol \Sigma_{\boldsymbol \alpha}$ between the HMC algorithm and the
parameter expansion technique are given in Figure~\ref{fig:hmc vs px}.  Using the posterior distribution from
the HMC algorithm as a benchmark, we notice that the parameter expansion method tends to
overestimate the variance of the random effects and the posterior densities are flatter.  The
observed differences are due to the PX technique, not the choice of kernel density bandwidth.
\begin{figure}[h!]
\onehalfspacing
\centering
\includegraphics[trim= 0.5cm 1cm 1.5cm 1.75cm,clip,width=0.75\textwidth]{hmc_vs_px_hiw}
\caption{Kernel density plots comparing the posterior density of the variance parameter $\sigma^2_{\alpha_i}$ in $\boldsymbol \Sigma_{\boldsymbol \alpha}$ between the HMC algorithm (blue solid) and the parameter expansion technique (red dashed).}
\label{fig:hmc vs px}
\end{figure}

redpajama_set_name: RedPajamaArXiv


Free BBQ food: Burgers, dogs, and more!
-Kicker is a Punk band from Oakland. Featuring current and past members of Neurosis, and Dystopia, along with veteran road dog Pete The Roadie on vocals. We're here to drink your beer, eat your chips, and bludgeon your ears with Punk.
-High and Tight is an American rock band formed and residing in NorCal's Bay Area. Influences are drawn from their collective musical backgrounds including; garage punk, california skate rock, street punk, 60s texas garage, surf, Boston Rock, 70s/80s metal, and Little Richard (of course). Since its inception in 2009, High and Tight has played over 100 shows up and down the west coast and been an integrel part of the San Francisco/Oakland music scene.


text:
Find Out Why Jennifer Lawrence Cursed Bradley Cooper Backstage At The SAG Awards
So far it's been a banner awards season for the cast of American Hustle, taking home Golden Globes as well as the SAG Award for Outstanding Performance By a Cast in a Motion Picture. OK! News: Soph...
OK! Magazine
okmagazine.com

redpajama_set_name: RedPajamaCommonCrawl