2022

[1]Mark J. Meyer

[1]Department of Mathematics and Statistics, Georgetown University,
Washington, DC, 20057, USA

2]Duke University School of Medicine, Durham, NC, 27710, USA

# Bayesian Analysis of Multivariate Matched Proportions with Sparse Response

<EMAIL_ADDRESS>Haobo Cheng  Katherine Hobbs Knutson  *  [

###### Abstract

Multivariate matched proportions (MMP) data appears in a variety of contexts
including post-market surveillance of adverse events in pharmaceuticals,
disease classification, and agreement between care providers. It consists of
multiple sets of paired binary measurements taken on the same subject. While
recent work proposes non-Bayesian methods to address the complexities of MMP
data, the issue of sparse response, where no or very few “yes” responses are
recorded for one or more sets, is unaddressed. The presence of sparse response
sets results in underestimates of variance, loss of coverage, and lowered
power in existing methods. Bayesian methods have not previously been
considered for MMP data but provide a useful framework when sparse responses
are present. In particular, the Bayesian probit model provides an elegant
solution to the problem of variance underestimation. We examine three
approaches built on that model: a naïve analysis with flat priors, a penalized
analysis using half-Cauchy priors on the mean model variances, and a
multivariate analysis with a Bayesian functional principal component analysis
(FPCA) to model the latent covariance. We show that the multivariate analysis
performs well on MMP data with sparse responses and outperforms existing non-
Bayesian methods. In a re-analysis of data from a study of the system of care
(SOC) framework for children with mental and behavioral disorders, we are able
to provide a more complete picture of the relationships in the data. Our
analysis provides additional insights into the functioning on the SOC that a
previous univariate analysis missed.

###### keywords:

Multivariate Probit Regression, Bayesian Inference, Bayesian FPCA, Penalized
Bayesian Regression, Systems of Care, Pediatric Mental and Behavioral
Disorders

##  1 Introduction

Multivariate matched proportions (MMP) data arises when multiple sets of
paired binary measurements are taken on the same subject. An example of such
data comes from a study of care coordination within a System of Care (SOC)
framework (Knutson2018, ). An SOC informs the treatment of children with
severe mental and behavioral disorders by coordinating care between six
different components: mental health, primary care, the education system, child
welfare, juvenile justice, and developmental disability services. However,
primary care frequently provides the sole treatment for children with mental
health disorders making it unclear if those children are receiving the
potential benefits of the SOC. In particular, primary care may not be
initiating contact with SOC components at the same rate as psychiatric care.
To begin to investigate primary care’s tendency to initiate coordination
within an SOC, Knutson et al. Knutson2018  present the results of a
retrospective chart study of youth referred by pediatricians to a child
psychiatrist at an urban community health center. The authors assess
coordination by examining the documented contacts in the medical record
between primary care and each of the remaining components of the SOC. The
study then compares those contacts to the documented contacts between the
psychiatrist, or specialty care, and the SOC components made after the initial
psychiatric evaluation.

We summarize the contact data for all of the SOC components in Table 1 where
the column headers denote the assessor, primary care or specialty care, and
the column sub-headers denote the components. The original analysis of this
data did not account for the multivariate nature of the data nor did it
account for the sparse responses in contacts made with two components. Sparse
response is most noticeable for contact made with developmental disabilities
by specialty care where no contacts were made (see Table 1). The sparse
response of juvenile justice requires noting that not all combinations of
zeros and ones appear in the data, only $0,0$ and $0,1$. Alternatively, we can
represent the data in a series of population averaged tables, as in Table 2,
and diagnose sparse response by noting the zero cell counts in both
developmental disabilities and juvenile justice tables. Since each component
has a corresponding set of paired binary measurements, two of the five sets in
this data exhibit sparse response.

Table 1: Multivariate matched proportions data from the SOC study for care-specific contact with developmental disabilities (DD), mental health (MH), juvenile justice (JJ), child welfare (CW), and education system (ED). Table values are 1 if contact was made with the agency in the column sub-header by the assessor in the column header, and 0 if no contact was made. The Pattern Count column denotes the number of times the pattern appears in the data. This table excludes patterns that do not appear in the data. Primary Care |  | Specialty Care |  | Pattern
---|---|---|---|---
DD | MH | JJ | CW | ED |  | DD | MH | JJ | CW | ED |  | Count
0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 1 |  | 17
0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 0 |  | 16
0 | 1 | 0 | 0 | 0 |  | 0 | 1 | 0 | 0 | 0 |  | 5
0 | 0 | 0 | 0 | 1 |  | 0 | 0 | 0 | 0 | 1 |  | 5
0 | 1 | 0 | 0 | 0 |  | 0 | 1 | 0 | 0 | 1 |  | 3
0 | 0 | 0 | 0 | 0 |  | 0 | 1 | 0 | 0 | 0 |  | 3
0 | 0 | 0 | 0 | 0 |  | 0 | 1 | 0 | 0 | 1 |  | 3
0 | 1 | 0 | 1 | 0 |  | 0 | 0 | 0 | 0 | 0 |  | 2
0 | 0 | 0 | 0 | 1 |  | 0 | 0 | 0 | 0 | 0 |  | 2
0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 1 | 1 |  | 2
1 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 0 |  | 1
0 | 1 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 0 |  | 1
0 | 1 | 1 | 0 | 0 |  | 0 | 1 | 1 | 1 | 0 |  | 1
0 | 1 | 0 | 0 | 1 |  | 0 | 1 | 0 | 1 | 1 |  | 1
0 | 1 | 0 | 1 | 0 |  | 0 | 1 | 0 | 1 | 1 |  | 1
0 | 0 | 0 | 0 | 1 |  | 0 | 1 | 0 | 0 | 1 |  | 1
0 | 1 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 1 |  | 1
0 | 0 | 0 | 0 | 1 |  | 0 | 1 | 0 | 1 | 1 |  | 1
0 | 0 | 0 | 0 | 1 |  | 0 | 0 | 1 | 1 | 1 |  | 1
1 | 0 | 0 | 0 | 1 |  | 0 | 0 | 0 | 0 | 1 |  | 1
0 | 0 | 0 | 1 | 1 |  | 0 | 1 | 0 | 0 | 1 |  | 1
0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 1 | 0 |  | 1
1 | 0 | 0 | 0 | 0 |  | 0 | 1 | 0 | 0 | 1 |  | 1
0 | 0 | 0 | 1 | 0 |  | 0 | 1 | 0 | 0 | 1 |  | 1
0 | 0 | 0 | 1 | 0 |  | 0 | 1 | 0 | 1 | 1 |  | 1
0 | 1 | 0 | 1 | 0 |  | 0 | 1 | 1 | 1 | 1 |  | 1
Table 2:  Population averaged tables for all SOC components. ‘Yes’ denotes
contact was made with the component, ‘No’ indicates contact was not made.
Within each table, the diagonal counts within each table represent the number
of concordant pairs while the off-diagonal counts represent the number of
discordant pairs.

Development Disabilities |  | Mental Health
---|---|---
Specialty | Primary Care |  | Specialty | Primary Care
Care | Yes | No |  | Care | Yes | No
Yes | $n_{111}=0$ | $n_{121}=0$ |  | Yes | $n_{112}=12$ | $n_{122}=12$
No | $n_{211}=3$ | $n_{221}=71$ |  | No | $n_{212}=4$ | $n_{222}=46$

Child Welfare |  | Education
---|---|---
Specialty | Primary Care |  | Specialty | Primary Care
Care | Yes | No |  | Care | Yes | No
Yes | $n_{114}=3$ | $n_{124}=7$ |  | Yes | $n_{115}=11$ | $n_{125}=31$
No | $n_{214}=4$ | $n_{224}=60$ |  | No | $n_{215}=2$ | $n_{225}=30$

Juvenile Justice
---
Specialty | Primary Care
Care | Yes | No
Yes | $n_{113}=1$ | $n_{123}=2$
No | $n_{213}=0$ | $n_{223}=71$

Several authors consider non-Bayesian approaches for analyzing MMP data.
Klingenberg and Agresti Kling2006  use marginal probability models to test for
simultaneous marginal heterogeneity. They employ a generalized estimating
equation (GEE) with an identity link and working independence to estimate
effects and construct a multivariate version of McNemar’s test (McNemar1947,
). Some examine methods for only pairs of correlated proportions Consonni2008
; Saeki2017 . Others propose multiple testing approaches to account for
correlation between sets of matched proportions (Westfall2010, ; Xu2013, ).
For example, Westfall, Troendle, and Pennello Westfall2010  obtain loose
control of the family wise error rate using Bonferroni-Holm on bootstrapped
$p$-values. Lui and Chang Lui2013  show that mixed-effects exponential risk
models can be used to estimate the risk ratios for each set of matched
proportions based on consistent estimators of the risk. In later work, the
same authors consider a mixed-effects logistic regression model with separate
random effects of pair and set Lui2016 . The authors show the estimated effect
of treatment for each outcome using this approach reduces to the $\log$ of the
well-known univariate Cochran-Mantel-Haenszel (CMH) estimator for stratified
$2\times 2$ tables (Cochran1950, ; Mantel1959, ), after integrating out the
random effects. Jiang and Xu Jiang2017  explore power and sample size
calculations for this and the GEE-based marginal probability model. As is
common in the analysis of matched proportions, each of these methods exclude
the concordant pairs, relying only on information from the discordant pairs
(Agresti2013, ).

The literature on Bayesian methods for matched proportions is limited to the
univariate case. For example, Altham Altham1971  examines closed form results
for a multinomial model with Dirichlet priors while Broemeling and Gregurich
Broemeling1996  discuss the Gibbs sampling approach to this model. Ghosh et
al. Ghosh2000  present an item response model with logit, probit, and
complementary $\log$-$\log$ links. This formulation allows for a hierarchical
Bayesian modeling of the data. The authors show that, for such models, the
concordant pairs also contribute to estimation and inference. While these
approaches are useful for univariate matched proportions, to our knowledge, no
previous work considers Bayesian methods for the analysis of MMP data.

More importantly, two of the existing methods for MMP data cannot directly
handle sparse responses, as we will now illustrate. The GEE-based approach
estimates marginal probabilities using the identity link under the assumption
of working independence (Kling2006, ). Using a model with one parameter for
each marginal probability, the authors show that the estimates equal the
sample probabilities of the columns of Table 1. But when a column is sparse,
as is the case for the developmental disabilities component of the SOC data,
the estimated probability is zero. The estimated variance will also be zero,
as will the associated covariances. The variance of the difference in marginal
probabilities between care-type depends on both this variance and the
corresponding covariance of the sparse column, resulting in underestimation of
the variance of the difference. The bootstrap-based method pre-differences the
paired outcomes, i.e. takes the difference between paired columns of Table 1,
and estimates the differences in the marginal probabilities using sample
averages (Westfall2010, ). In general, these differences can take on the
values $-1,0,$ and 1. But when the response is sparse, as in the developmental
disabilities component, the differences can only be 0 or 1 meaning the lower
bound of a bootstrapped interval is bounded below by 0. Similarly, for the
juvenile justice component, only $-1$ and 0 are possible realizations from the
data resulting in a bootstrapped interval that is bounded above by 0. Once
again, the uncertainty of the difference in marginal probabilities for sparse
outcomes is underestimated.

Lui and Chang do both provide adjustments for sparse data in both their
exponential risk model (ERM) and their mixed effects regression approach
(Lui2013, ; Lui2016, ). The ERM approach, for example, bases its estimation
off of the $K$-specific tables in Table 2. Specifically, the estimate is the
ratio of the counts of discordant pairs in each table:
$n_{21k}\big{/}n_{12k}$. When one of these elements is zero, the authors
suggest adding $1/2$ to the numerator and the denominator making their
sparsity adjusted estimator
$\left(n_{21k}+\frac{1}{2}\right)\big{/}\left(n_{12k}+\frac{1}{2}\right)$. Lui
and Chang do not evaluate the performance of their estimators in the presence
of sparse response and instead evaluate the performance of their testing
procedure in non-sparse simulated data settings (Lui2013, ; Lui2016, ). Such a
correction is philosophically Bayesian in its thinking as it recognizes, a
priori, that the true population parameter is likely not zero (or infinite)
and adjusts the estimate accordingly. In fact, this estimator is an empirical
Bayes estimate of the marginal risk for a single matched set under Jeffreys’
prior when using the Dirichlet model proposed by Altham1971 .

Adjustments of this kind are not generalizable either: it is not applicable to
the GEE or Bootstrap approaches as it adjusts the counts from Table 2 whereas
both the GEE and Bootstrap rely on the individual-level data as in Table 1.
Other common adjustments like data augmentation, e.g. adding one success and
one failure in the case of binary data, do not easily generalize to the
multivariate case either. Augmenting sparse cells in Table 2 separately
ignores the multivariate nature of the data, breaking potential underlying
correlations. Adding covariate patterns to Table 1 must take into account the
impact on the multivariate structure and can lead to severe sample size
inflation. At the extreme, augmentation would involve adding one of every
possible type of covariate pattern to Table 1. For even a small number of
sets, such an addition would greatly inflate the sample size far beyond what
augmentation strategies aim for, potentially biasing the results. Thus, there
is a need for an MMP method that can accommodate sparse responses more
adequately in addition to the need for a Bayesian approach.

Using the Albert and Chib AlbertChib1993 ; AlbertChib1995  probit model as a
basis, we consider several approaches for the Bayesian modeling of
multivariate matched proportions. The latent variable representation of the
probit model solves the issue of variance estimation in the existing methods:
by convention, the latent-space variance is fixed at one. This alone may not
be sufficient in the presence of sparse response where estimates can easily
diverge, so we also consider shrinkage priors on model components. Finally,
noting that the first two approaches assume independence in the latent space,
we examine a multivariate analysis. We refer to the three Bayesian models we
examine as the naïve analysis which uses a simple Bayesian probit model, the
penalized analysis which uses half-Cauchy priors on the variance components of
a global ridge-like prior Gelman2006 , and the multivariate analysis. Our
multivariate analysis makes use of a Bayesian functional principal components
analysis or FPCA VanDerLinde2008  to model a general latent covariance matrix.
The Bayesian FPCA can accommodate a wide range of covariance structures, not
just those arising from functional data, making it an appealing approach for
dealing MMP data.

In this manuscript, we demonstrate that all Bayesian models have better
coverage and power in the presence of sparse responses than the existing
methods for MMP data while maintaining similar bias and interval widths.
Comparing between the Bayesian models, we find that the multivariate analysis
tends to perform best in terms of coverage, although the naïve and penalized
analyses do perform similarly well. Finally, we reanalyze the SOC data using
the multivariate analysis and include the two components with sparse responses
which were excluded from the original analysis. Our reanalysis confirms some
of the original results but also identifies additional findings the univariate
analysis missed.

##  2 Methods

We begin by defining an overarching notation for the models we consider. At
the pair-level, let $x_{ijk}$ denote the response for pair $i$, at observation
$j$, for the $k$th set where $i=1,\ldots,n$, $j=1,2$, and $k=1,\ldots,K$. In
context of the SOC data, $i$ denotes the subjects and $j$ denotes care-type,
primary care or specialty care. The index $k$ denotes the SOC component with
which contact is made, $x_{ijk}=1$, or not made, $x_{ijk}=0$. Combining across
$k$, we let ${\bf
x}_{ij}=\left[\begin{array}[]{ccc}x_{ij1}&\cdots&x_{ijK}\end{array}\right]^{\prime}$
be the vector of responses at observation $j$ for pair $i$. Subject $i$’s
stacked vector of responses is then ${\bf x}_{i}=\left[\begin{array}[]{cc}{\bf
x}_{i1}&{\bf x}_{i2}\end{array}\right]^{\prime}$. Realizations of this vector
make up the rows of Table 1.

The target of inference is a vector consisting of the differences in the
$k$-specific marginal probabilities of contact, $P(x_{ijk}=1)=\theta_{jk}$,
which both the GEE-based Kling2006  and bootstrap-based Westfall2010  methods
estimate. First let
$\mbox{\boldmath$\theta$}=\left[\begin{array}[]{cccccc}\theta_{11}&\cdots&\theta_{1K}&\theta_{21}&\cdots&\theta_{2K}\end{array}\right]$
be the vector of marginal probabilities. Then let $\rho_{k}$ be the $k$th
difference in the marginal probabilities, thus
$\rho_{k}=\theta_{1k}-\theta_{2k}$. The vector of $\rho_{k}$’s is $\rho$,
$\mbox{\boldmath$\rho$}=\left[\begin{array}[]{ccc}\rho_{1}&\cdots&\rho_{K}\end{array}\right]$.
Alternatively, define ${\bf L}$ to be a block matrix of the form ${\bf
L}=\left({\bf I}_{K}|-{\bf I}_{K}\right)$ where ${\bf I}_{K}$ is a $K\times K$
identity matrix. The vector representing the differences in marginal
probabilities is then $\mbox{\boldmath$\rho$}={\bf
L}\mbox{\boldmath$\theta$}$.

###  2.1 Naïve Analysis

We first consider a naïve analysis using probit regression as first described
by Albert and Chib AlbertChib1993 ; AlbertChib1995 . The regression framework
makes use of the latent variable representation of the probit model. In the
Bayesian context, this approach can be considered a data augmentation step
where the latent variables are the augmented data Gelman2013 . We define the
following mapping

$\displaystyle x_{ijk}=\left\\{\begin{array}[]{cc}0&\text{ if }z_{ijk}<0\\\
1&\text{ if }z_{ijk}\geq 0\end{array}\right.,$ (3)

where the latent space model for MMP data is
$z_{ijk}=\beta_{jk}+\epsilon_{ijk}$. In general, $\epsilon_{ijk}\sim
N(0,\sigma_{z}^{2})$ and, by convention, $\sigma_{z}^{2}$ is fixed at 1.
Because the model is saturated, the relationship between the latent space
coefficients, $\beta_{jk}$, and the marginal probabilities, $\theta_{jk}$, is
$\beta_{jk}=\Phi^{-1}(\theta_{jk})$ where $\Phi(\cdot)$ is the standard normal
CDF. This relationship, and indeed the probit model itself, is induced by
assuming the error terms are normal. Letting
$\mbox{\boldmath$\beta$}=\left[\begin{array}[]{cccccc}\beta_{11}&\cdots&\beta_{1K}&\beta_{21}&\cdots&\beta_{2K}\end{array}\right]$,
we place the a flat, non-informative prior on $\beta$, specifically
$\pi(\mbox{\boldmath$\beta$})\propto 1$. The resulting posterior conditionals
are:

$\displaystyle z_{ijk}$ $\displaystyle|x_{ijk}=0,\beta_{jk}\sim
N(\beta_{jk},1)_{\\{z_{ijk}<0\\}}$ $\displaystyle z_{ijk}$
$\displaystyle|x_{ijk}=1,\beta_{jk}\sim N(\beta_{jk},1)_{\\{z_{ijk}\geq 0\\}}$
$\beta$ $\displaystyle|{\bf Z},{\bf W}\sim N\left[({\bf W}^{\prime}{\bf
W})^{-1}{\bf W}^{\prime}{\bf Z},({\bf W}^{\prime}{\bf W})^{-1}\right]$

where ${\bf W}$ is a design matrix built by stacking $2K\times 2K$ identity
matrices on top of each other, one for each subject, and ${\bf Z}$ is the
vector of all $z_{ijk}$’s. The subscript notation $\\{\\}$ on the normal
densities indicates truncation to the region defined in the brackets.
Estimates of the difference in marginal probabilities, $\rho$, are obtained
from the transformed posterior samples of $\beta$:
$\mbox{\boldmath$\rho$}={\bf L}\Phi(\mbox{\boldmath$\beta$})$ where the
function $\Phi(\cdot)$ is applied to the vector $\beta$ in component-wise
fashion.

In the above conditionals, $\sigma_{z}^{2}$ is fixed at one. Since the latent
variable representation is scale-invariant, $\sigma_{z}^{2}$ is not
necessarily identifiable. Fixing $\sigma_{z}^{2}$ at one solves the
identifiability issue in the latent space. But this convention provides an
additional benefit in the presence of sparse response by stabilizing the model
variance. Without this constraint, the posterior chains for $\beta_{jk}$ and,
subsequently, $\rho_{k}$, would not converge in the presence of sparse
responses. Even with this constraint, the posterior chains can be unstable
because in the latent space, the support of each $\beta_{jk}$ is $\mathbb{R}$.
To prevent the model components from diverging, we penalize the $\beta_{jk}$
and shrink them toward zero.

###  2.2 Penalized Analysis

There are many choices for shrinkage options in the Bayesian framework, but we
select normal priors with a half-Cauchy hyper-prior on the scale. Gelman
Gelman2006  proposed the use of half-$t$ priors for variance components, of
which the half-Cauchy is a special case. Polson and Scott Polson2012  show
that the half-Cauchy performs well as global-shrinkage prior and recommend
their routine use in normal hierarchical models. The model specification for
the penalized analysis begins in the same fashion with the mapping in Equation
(3) and with the same latent model: $z_{ijk}=\beta_{jk}+\epsilon_{ijk}$ for
$\epsilon_{ijk}\sim N(0,1)$. Instead of the flat prior, we place the following
global hierarchical prior on $\beta_{jk}$:

$\displaystyle\beta_{jk}$ $\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\lambda)$ $\displaystyle\lambda$ $\displaystyle\sim
IG\left(\frac{1}{2},\frac{1}{\mu}\right)$ (4) $\displaystyle\mu$
$\displaystyle\sim IG\left(\frac{1}{2},\frac{1}{A^{2}}\right)$

where $IG$ denotes the inverse gamma distribution and $A$ is a fixed hyper-
parameter. These hierarchical inverse gamma priors are the mixture
representation of the half-Cauchy prior Wand2011 . That is, if $\lambda\sim
IG\left(1/2,1/\mu\right)$ and $\mu\sim IG\left(1/2,1/A^{2}\right)$, then
$\sqrt{\lambda}\sim HC(A)$ where $HC$ denotes the half-Cauchy.

The mixture-model representation of the half-Cauchy induces a model with fully
identifiable conditional posterior distributions:

$\displaystyle z_{ijk}$ $\displaystyle|x_{ijk}=0,\beta_{jk}\sim
N(\beta_{jk},1)_{\\{z_{ijk}<0\\}}$ $\displaystyle z_{ijk}$
$\displaystyle|x_{ijk}=1,\beta_{jk}\sim N(\beta_{jk},1)_{\\{z_{ijk}\geq 0\\}}$
$\beta$ $\displaystyle|Z,{\bf W}\sim N\left[({\bf W}^{\prime}{\bf
W}+\lambda^{-1}{\bf I}_{2K})^{-1}{\bf W}^{\prime}{\bf Z},({\bf W}^{\prime}{\bf
W}+\lambda^{-1}{\bf I}_{2K})^{-1}\right]$ $\displaystyle\lambda$
$\displaystyle|\mbox{\boldmath$\beta$},\mu\sim
IG\left(K+\frac{1}{2},\frac{1}{\mu}+\frac{1}{2}\mbox{\boldmath$\beta$}^{\prime}\mbox{\boldmath$\beta$}\right)$
$\displaystyle\mu$ $\displaystyle|\lambda\sim
IG\left(1,\frac{1}{A^{2}}+\frac{1}{\lambda}\right)$

where ${\bf I}_{2K}$ is the $2K\times 2K$ identity matrix, ${\bf W}$ is as
defined in Section 2.1. This model provides global shrinkage of the latent-
space coefficients, $\beta_{jk}$, across all $j$ and $k$ and stabilizes the
chain when the response is sparse. Estimates of $\rho$ are once again obtained
from the transformed samples of $\beta$ using $\mbox{\boldmath$\rho$}={\bf
L}\Phi(\mbox{\boldmath$\beta$})$.

The choice of the hyper-parameter $A$ can be regarding as placing a weak upper
bound on the variance of $\beta$ Gelman2006 . The parameter can be interpreted
in context as the standard deviation of the corresponding model component
Gelman2006 . Since we are modeling in a latent space, our choice of the value
of $A$ is less clear as the scale is not directly dependent on a covariate.
Thus, we select a relatively large value, $A=10$, to give a diffuse prior.

###  2.3 Multivariate Analysis

Both of the preceding methods treat the data as independent univariate
samples, however each subject’s ${\bf x}_{i}$ may have a multivariate
distribution. Chib and Greenberg ChibGreenberg1998  discuss a multivariate
version of the Albert and Chib AlbertChib1993 ; AlbertChib1995  model which,
in the MMP context, begins with the same mapping as in Equation (3). The
latent model is then ${\bf
z}_{i}=\mbox{\boldmath$\beta$}+\mbox{\boldmath$\epsilon$}_{i}$ where ${\bf
z}_{i}=\left[\begin{array}[]{cccccc}z_{i11}&\cdots&z_{i1K}&z_{i21}&\cdots&z_{i2K}\end{array}\right]^{\prime}$
and
$\mbox{\boldmath$\epsilon$}_{i}=\left[\begin{array}[]{cccccc}\epsilon_{i11}&\cdots&\epsilon_{i1K}&\epsilon_{i21}&\cdots&\epsilon_{i2K}\end{array}\right]^{\prime}$.
The latent error vector, $\mbox{\boldmath$\epsilon$}_{i}$, is then assumed to
be multivariate normal, $\mbox{\boldmath$\epsilon$}_{i}\sim N({\bf
0},\mbox{\boldmath$\Sigma$})$. The matrix $\Sigma$ is a correlation matrix.

Drawing posterior samples for $\Sigma$ can be difficult. Chib and Greenberg
suggest a Metropolis-Hastings step to obtain draws from the posterior of
$\Sigma$ and others authors have further discussed the same issue Liu2001 ;
Zhang2006 ; Webb2008 . In all of these approaches, the MCMC sampler also
requires obtaining draws from a truncated multivariate normal posterior, the
posterior conditional distribution of ${\bf z}_{i}$, for each subject. Using
existing approaches for random drawing from truncated multivariate normal
distributions is time intensive. In our preliminary model exploration, a model
of this kind took over 11 hours to generate 20,000 posterior samples when
$K=2$ (run on a computer with a 3.7 GHz 6-Core Intel Core i5 processor and 32
GB of memory).

To accommodate a multivariate structure without inducing such a computational
burden, we turn to recent developments in the functional data literature. In
particular, Goldsmith and Kitago Goldsmith2016  and Meyer et al.
MeyerMorris2022  both show that the Bayesian FPCA VanDerLinde2008 , can be
used to model a variety of covariance structures for both Gaussian
Goldsmith2016  and ordinal MeyerMorris2022  functional data including
independent, exponential, and compound symmetric (also known as exchangeable)
covariance structures. For MMP data, we may expect exchangeability within
observation level $j$, what might be called a block compound symmetric
structure, or independence within $j$ but dependence between $k$. A more
general form or unstructured form may result as well. An additional challenge
is that it is hard to know what the structure of the latent correlation matrix
might be. Thus, for the multivariate analysis, we propose using an FPCA in the
latent space to model possible multivariate associations in MMP data.

Our proposed model performs the FPCA on the latent errors setting
$\mbox{\boldmath$\epsilon$}_{i}=\mbox{\boldmath$\Omega$}\mbox{\boldmath$\psi$}{\bf
c}_{i}^{\prime}+\mbox{\boldmath$\varepsilon$}_{i}$ for the $L\times 1$ vector
of subject scores ${\bf c}_{i}$, the $2K\times L$ matrix of basis coefficients
$\Psi$, and the $2K\times 2K$ matrix of basis-functions $\Omega$. The model
for each subject’s latent vector is then

$\displaystyle{\bf
z}_{i}=\mbox{\boldmath$\beta$}+\mbox{\boldmath$\Omega$}\mbox{\boldmath$\Psi$}{\bf
c}_{i}^{\prime}+\mbox{\boldmath$\varepsilon$}_{i},$

where $\mbox{\boldmath$\varepsilon$}_{i}$ is a $2K\times 1$ vector of normal,
independent errors with variance equal to $\sigma_{\varepsilon}^{2}{\bf
I}_{2K}$. Consistent with previous work using Bayesian FPCA, we use $L=2$
scores for the expansion Goldsmith2016 ; MeyerMorris2022 . Letting $\Omega$ be
a matrix of cubic B-spline basis functions, the prior on $\Psi$ is then
$\mbox{\boldmath$\Psi$}\sim
N[0,(\mbox{\boldmath$\Lambda$}_{\psi}\otimes\mbox{\boldmath$\mathcal{P}$})^{-1}]$
where $\mbox{\boldmath$\Lambda$}_{\psi}$ is an $L\times L$ diagonal matrix of
tuning parameters $\lambda_{\ell}^{-1}$ and $\mathcal{P}$ is a $2K\times 2K$
penalty matrix of the form
$\mbox{\boldmath$\mathcal{P}$}=\xi\mbox{\boldmath$\mathcal{P}$}_{0}+(1-\xi)\mbox{\boldmath$\mathcal{P}$}_{2}$.
The matrices $\mbox{\boldmath$\mathcal{P}$}_{0}$ and
$\mbox{\boldmath$\mathcal{P}$}_{2}$ are the zeroth and second derivative
penalty matrices, respectively. The value $\xi$ is a control parameter that
balances smoothness and shrinkage with values near zero favoring shrinkage—we
take $\xi=0.01$. Finally, the subject scores, ${\bf c}_{i}$, have mean zero
normal priors with variance equal to ${\bf I}_{L}$, an $L\times L$ identity
matrix.

We use the same prior on $\beta$ as in the penalized analysis, found in
Equation (4). Combining with the FPCA specification, we can obtain the full
conditional posterior distributions for each model component. Let $\Psi_{ijk}$
denote the $j,k$th sub-components of the vector
$\mbox{\boldmath$\Omega$}\mbox{\boldmath$\Psi$}{\bf c}_{i}^{\prime}$. The full
conditionals are

$\displaystyle z_{ijk}$ $\displaystyle|x_{ijk}=0,\beta_{jk},\Psi_{ijk}\sim
N(\beta_{jk}+\Psi_{ijk},1)_{\\{z_{ijk}<0\\}}$ $\displaystyle z_{ijk}$
$\displaystyle|x_{ijk}=1,\beta_{jk},\Psi_{ijk}\sim
N(\beta_{jk}+\Psi_{ijk},1)_{\\{z_{ijk}\geq 0\\}}$ $\beta$
$\displaystyle|Z,{\bf W}\sim
N\left[\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\beta$}}{\bf W}^{\prime}({\bf
Z}-\mbox{\boldmath$\psi$}),\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\beta$}}\right]$
$\displaystyle\lambda$ $\displaystyle|\mbox{\boldmath$\beta$},\mu\sim
IG\left(K+\frac{1}{2},\frac{1}{\mu}+\frac{1}{2}\mbox{\boldmath$\beta$}^{\prime}\mbox{\boldmath$\beta$}\right)$
$\displaystyle\mu$ $\displaystyle|\lambda\sim
IG\left(1,\frac{1}{A^{2}}+\frac{1}{\lambda}\right)$ $\Psi$ $\displaystyle|\sim
N\left[\frac{1}{\sigma_{\varepsilon}^{2}}\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\Psi$}}({\bf
C}\otimes\mbox{\boldmath$\Omega$})^{\prime}({\bf Z}-{\bf
W}\mbox{\boldmath$\beta$}),\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\Psi$}}\right]$
$\displaystyle{\bf c}_{i}$ $\displaystyle|\sim
N\left[\frac{1}{\sigma_{\varepsilon}^{2}}\mbox{\boldmath$\Sigma$}_{{\bf
c}}(\mbox{\boldmath$\Omega$}\mbox{\boldmath$\Psi$})^{\prime}({\bf
z}_{i}-\mbox{\boldmath$\beta$}),\mbox{\boldmath$\Sigma$}_{{\bf c}}\right]$
$\displaystyle\sigma_{\varepsilon}^{2}$ $\displaystyle|\sim
IG\left[1+nK,1+\frac{1}{2}({\bf Z}-{\bf
W}\mbox{\boldmath$\beta$}-\mbox{\boldmath$\psi$})^{\prime}({\bf Z}-{\bf
W}\mbox{\boldmath$\beta$}-\mbox{\boldmath$\psi$})\right]$
$\displaystyle\lambda_{\ell}$ $\displaystyle|\sim
IG\left(2K,K+\frac{1}{2}\mbox{\boldmath$\Psi$}_{\ell}^{\prime}\mbox{\boldmath$\mathcal{P}$}\mbox{\boldmath$\Psi$}_{\ell}\right)$

for
$\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\beta$}}=\left(\frac{1}{\sigma_{\varepsilon}^{2}}{\bf
W}^{\prime}{\bf W}+\lambda^{-1}{\bf I}_{2K}\right)^{-1}$,
$\mbox{\boldmath$\Sigma$}_{\mbox{\boldmath$\Psi$}}=\left(\frac{1}{\sigma_{\varepsilon}^{2}}{\bf
C}^{\prime}{\bf
C}\otimes\mbox{\boldmath$\Omega$}^{\prime}\mbox{\boldmath$\Omega$}+\mbox{\boldmath$\Lambda$}_{\psi}\otimes\mbox{\boldmath$\mathcal{P}$}\right)^{-1}$,
$\mbox{\boldmath$\Sigma$}_{{\bf
c}}=\left[\frac{1}{\sigma_{\varepsilon}^{2}}(\mbox{\boldmath$\Omega$}\mbox{\boldmath$\Psi$})^{\prime}(\mbox{\boldmath$\Omega$}\mbox{\boldmath$\Psi$})+{\bf
I}_{L}\right]^{-1}$, ${\bf C}$ equal to the matrix of ${\bf c}_{i}$ vectors
stacked on top of each other, $\mbox{\boldmath$\Psi$}_{\ell}$ equal to the
$\ell$th column of $\Psi$, and $\psi$ the vector containing all $\Psi_{ijk}$s.
Note that we draw samples from $\sigma_{\varepsilon}^{2}$ to improve the
estimation of other model components, but fix $\sigma_{\varepsilon}^{2}$ at
one when sampling the latent variable.

This expansion allows us to treat the latent variables as we would in the
univariate analysis and does not require sampling from a multivariate
truncated normal nor does it necessitate a Metropolis-Hastings step to draw
latent correlations. As with the naïve and penalized analyses, estimates of
$\rho$ are easily obtained from the posterior of $\beta$ with
$\mbox{\boldmath$\rho$}={\bf L}\Phi(\mbox{\boldmath$\beta$})$. This model is
much more computationally efficient than other multivariate probit models
taking between 60 and 70 seconds, depending on $K$, to generate 20,000
posterior samples (on the same computer, one with a 3.7 GHz 6-Core Intel Core
i5 processor and 32 GB of memory).

##  3 Simulation Study

We base our simulation study off of the sample size of and effect sizes
observed in the SOC data. For each dataset, we generate data for $n=75$
subjects where $K=2,3,4,$ and 5 sets of binary measurements are taken on each
subject. We aim to induce sparse responses within sets and achieve this by
first generating probabilities for $\theta$ from a multivariate normal
distribution. In the event a draw from this distribution is negative, it is
set to zero. Values of ${\bf x}_{i}$ are then drawn from a Bernoulli
distribution. The “true” values of $\theta$ vary with $K$:

$\displaystyle K=2$
$\displaystyle:\mbox{\boldmath$\theta$}=\left[\begin{array}[]{cccc}0.05&\theta_{12}&0.25&0.005\end{array}\right],$
$\displaystyle K=3$
$\displaystyle:\mbox{\boldmath$\theta$}=\left[\begin{array}[]{cccccc}0.05&\theta_{12}&0.1&0.25&0.005&0.15\end{array}\right],$
$\displaystyle K=4$
$\displaystyle:\mbox{\boldmath$\theta$}=\left[\begin{array}[]{cccccccc}0.05&\theta_{12}&0.005&0.1&0.25&0.005&0.05&0.15\end{array}\right],\text{
and}$ $\displaystyle K=5$
$\displaystyle:\mbox{\boldmath$\theta$}=\left[\begin{array}[]{cccccccccc}0.05&\theta_{12}&0.005&0.1&0.25&0.25&0.005&0.05&0.15&0.35\end{array}\right].$

The $k=2$ component is set to exhibit sparse response since
$\theta_{22}=0.005$. We add an additional sparse responses for $K=4$ and 5,
although sparse responses may occur for any component in any given simulated
draw. To evaluate the effect of sparse response as the degree of difference
between observations changes, we vary $\theta_{12}$ from 0.05 to 0.2 by 0.1.
All of these values for $\theta$ encompass the range of values we see when
estimating $\rho$ for the SOC data.

To evaluate the MMP methods, we determine the coverage, power, bias, and
interval width for the $k=2$ component, $\rho_{2}=\theta_{21}-\theta_{22}$,
which is simulated to exhibit sparse response. We compare all three Bayesian
models from Section 2 to the GEE-based approach (Kling2006, ), the Bootstrap-
based method (Westfall2010, ), and a modified ERM approach (Lui2013, ). The
ERM estimates risk ratios Lui2013 , however, for comparison with the GEE,
Bootstrap, and the Bayesian approaches, we estimate the risk difference using
their model. All Bayesian estimates are based on 20,000 total samples,
discarding the first 10,000 while all Bootstrap estimates are based on 10,000
resamples. Additional details on the GEE, ERM, and Bootstrap approaches are in
the Section 1 of the Supplementary Materials.

For each combination of $K$ and $\theta_{12}$, we generate 1000 simulated
datasets, retaining only those datasets that exhibit sparsity for the second
set—as determined by the sparsity exhibited in either Table 1 or 2. The
parameters of the simulation result in between 294 and 422 datasets exhibiting
sparse responses (out of 1000), depending on the value of both $\theta_{12}$
and $K$. We then randomly sample 200 of these datasets for evaluation. All
Bayesian models perform similarly in terms of bias. The multivariate model has
higher coverage for most values of $\theta_{12}$ and the power advantage of
the simpler methods decreases as $K$ increases (see Section 2 of the
Supplementary Materials). For these reasons and to ensure we capture the
multivariate nature of the data, we suggest using the multivariate analysis.

Figure 1: Coverage of $\rho_{2}$ for the multivariate probit analysis (MVP),
exponential risk model (ERM), GEE-based model (GEE), and the Bootstrap-based
model (Boot) by $\theta_{12}$ and $K$. All intervals are 95% intervals.

Figure 1 contains the coverage comparison between the four models under
comparison for all values of $\theta_{12}$ and $K$. Coverage curves for the
multivariate probit analysis (MVP) are in solid gold with circles, the
exponential risk model (ERM) results are in dashed red with triangles, the
GEE-based model (GEE) results are in dotted green with plus signs, and the
Bootstrap-based model (Boot) results are in dotted-dashed blue with x’s.
Consistent with other evaluations of intervals for binary
data,AgrestiCoull1998 ; AgrestiCaffo2000  we find that the coverage for all
methods varies by effect size. The multivariate analysis has the highest
coverage for almost every combination of $\theta_{12}$ and $K$, is typically
at or above nominal (95%), and never drops below 92.5%. The comparator
methods, however, have coverage that is frequently below nominal and, in some
case, at or below 90%. Coverage is more inconsistent for the ERM and GEE,
particularly when $\theta_{12}$ is smaller. As $\theta_{12}$ increases, all
methods have coverage near nominal.

Figure 2: Power to detect $\rho_{2}$ for the multivariate probit analysis
(MVP), exponential risk model (ERM), GEE-based model (GEE), and the Bootstrap-
based model (Boot) by $\theta_{12}$ and $K$. All tests were conducted at the
5% level.

A gain in coverage can come at a loss in power, so we also investigate the
power of each method to test the null $H_{0}:\rho_{2}=0$. The results are in
Figure 2 where significance is determined by exclusion from the 95% interval
for the multivariate analysis and Bootstrap-based model and by the resulting
$p$-values being less than 0.05 for the GEE and ERM. Power curves for the
models follow the same color, line, and point schemes as in the depiction of
the coverage curves. The multivariate analysis has the highest power at each
$\theta_{12}$, regardless of the value of $K$. The GEE and Bootstrap are
second with nearly identical power while the ERM has the lowest power. All
models have power that goes to one as $\theta_{12}$ increases.

Figure 3: Bias in estimating $\rho_{2}$ for the multivariate probit analysis
(MVP), exponential risk model (ERM), GEE-based model (GEE), and the Bootstrap-
based model (Boot) by $\theta_{12}$ and $K$.

Figure 4: Width of intervals on $\rho_{2}$ for the multivariate probit
analysis (MVP), exponential risk model (ERM), GEE-based model (GEE), and the
Bootstrap-based model (Boot) by $\theta_{12}$ and $K$. All intervals are 95%
intervals.

Finally, we assess the bias of each method, in Figure 3, and interval width,
in Figure 4. The same color, line, and point schemes used in the coverage and
power curves is used for the bias and width curves Bias is virtually
identically between the four methods regardless of the value of $\theta_{12}$
or $K$. In general, bias tends to increase as $\theta_{12}$ and increases
slightly with $K$ as well, though more so with changes in the effect size.
Interval width is very similar between the four methods, although the
multivariate analysis tends to have the widest intervals followed by the ERM,
the GEE-based model, and then the Bootstrap-based method. For all methods,
width increases as $\theta_{12}$ increases but only exhibits minimal changes
as $K$ increases. Histograms of bias collapsing across $\theta_{12}$ are in
Section 2 of the Supplementary Materials.

##  4 Re-analysis of the SOC Data

The goal of the SOC data is to determine if primary care physicians are
initiating contact with the same agencies within the SOC framework, and with
similar frequency, to psychiatric, or specialty, care. Knutson et al. argues
that it is important to identify services for improved coordination to help
the overburdened mental health system and their original findings suggest
there is room for improvement among primary care givers Knutson2018 . They
separately analyze only the non-sparse sets using standard univariate methods
including logistic regression Knutson2018 . Because of the sparsity in the
juvenile justice and developmental disabilities components, the authors do not
examine these outcomes. Of the three outcomes they analyze, only contact with
the education system, is deemed significantly different Knutson2018 . We now
re-analyze all outcomes, simultaneously, using the multivariate analysis.

In total, we have 74 patients available for analysis. The most common primary
mental health diagnosis was mood disorder followed by attention deficit
hyperactivity disorder, adjustment disorder, and anxiety disorder Knutson2018
. Additional primary diagnoses include post-traumatic stress disorder,
substance abuse, and learning disorder Knutson2018 . Each patient was
initially seen by primary care, who made recommendations for contact (or not)
with components of the SOC. Separately, specialty care later examined each
patient and independently made recommendations for contact within the SOC. We
only seek to estimate and analyze $\rho$, the differences in the marginal
probability of contact between primary care and specialty care for each SOC
component. Table 3 presents the results of the multivariate analysis including
posterior median values for each $\rho_{k}$, 95% credible intervals (CrI), and
posterior probabilities that $\rho_{k}$ is greater than zero, $P(\rho_{k}>0)$.
In Table 3, we abbreviate each component as follows: developmental
disabilities (DD), mental health (MH), juvenile justice (JJ), child welfare
(CW), and education system (ED). In total, we generated 20,000 draws from the
posterior, retaining the last 10,000 for estimation and inference. All model
parameters were judged to have converged using trace plots (see Section 3 of
the Supplementary Materials) and the Gelman-Rubin diagnostic GelmanRubin1992
(see Table 3 where values of $\hat{R}$ and its upper 95% CrI near one suggest
convergence).

We begin by noting what is consistent in our analysis with the univariate
analysis from Knutson et al.Knutson2018 . The authors found a significant
association for the education component of the SOC whereby primary care
contacted the education system less than specialty care. We confirm this
result noting that the difference in the marginal probability of contact with
the education system between primary care and specialty care is
$\rho_{ED}=-0.394$ (95% CrI $[-0.533,-0.240]$, $P[\rho_{ED}>0]=0.000$).
Knutson et al.Knutson2018  also found no significant associations for the
child welfare and mental health components. Our analysis confirms this result
for contact with child welfare ($\rho_{CW}=-0.040$, 95% CrI $[-0.150,0.063]$,
$P[\rho_{CW}>0]=0.224$) but by analyzing the data using a multivariate
analysis, we have a more nuanced view of the relationship between care type
and contact with mental health. Specifically, the difference in contact with
mental health between primary care and specialty care is $\rho_{MH}=-0.110$
and its 95% CrI includes zero $(-0.259,0.035)$. But only 6.8% of the posterior
distribution of $\rho_{MH}$ sits above zero which indicates that a large
portion of the posterior distribution is negative.

Table 3: Results from the re-analysis of the SOC Data using the multivariate analysis. Gelman-Rubin $\hat{R}$ values and upper 95% credible bounds near one suggest convergence. $P(\rho>0)$ is the posterior probability that $\rho$ is greater than zero. Component abbreviations: developmental disabilities (DD), mental health (MH), juvenile justice (JJ), child welfare (CW), and education system (ED).  |  | Cred. Int. |  | Gelman-Rubin
---|---|---|---|---
Component | $\rho$ | 2.5% | 97.5% | $P(\rho>0)$ | $\hat{R}$ | Upper 95%
DD | 0.035 |  $-$0.002 | 0.096 | 0.970 | 1.003 | 1.008
MH |  $-$0.110 |  $-$0.259 | 0.035 | 0.068 | 1.000 | 1.001
JJ |  $-$0.024 |  $-$0.089 | 0.031 | 0.170 | 1.002 | 1.006
CW |  $-$0.040 |  $-$0.150 | 0.063 | 0.224 | 1.001 | 1.002
ED |  $-$0.394 |  $-$0.533 |  $-$0.240 | 0.000 | 1.000 | 1.001

Due to sparse response, contact with the developmental disabilities and
juvenile justice components was not assessed in the original univariate
analysis Knutson2018 . But using the proposed multivariate analysis, we can
conduct inference on both. We do not find a significant association in contact
with the juvenile justice component ($\rho_{JJ}=-0.024$, 95% CrI
$[-0.089,0.031]$, $P[\rho_{JJ}>0]=0.170$). The relationship between care type
and contact with the developmental disabilities component is more nuanced. We
find a positive effect suggesting that the marginal probability of contact is
$\rho_{DD}=0.035$ more for primary care than for specialty care. The 95% CrI
does exclude zero, $-0.002,0.096$, although it is on the edge and the
posterior probability that $\rho_{DD}$ is greater than zero is 97.0%. Thus, a
sizable portion of the posterior distribution of $\rho_{DD}$ is positive. Both
the naïve and penalized analyses provide similar inferential results to the
multivariate analysis (see Section 3 of the Supplementary Materials).

Overall, our analysis suggests there are inefficiencies in the system of care
framework, a consistent finding with that of Knutson et al. Knutson2018 . But
our results are more nuanced and suggest multiple areas for targeted
improvement. In the study, primary care is recommending contact with both the
education and mental health systems at lower rates than they should be;
potentially leaving important actors out of the care for children with mental
health disorders. Conversely, primary care is contacting disability services
more than specialty care; potentially overburdening those services with
children who do not need them. Many children only receive mental health care
through their primary care physicians,Costello2014  so it is important these
rates of contact to align to ensure the children benefit from the SOC
frameworkKnutson2018 .

##  5 Discussion

Methods for MMP data are limited to non-Bayesian approaches that do not
perform well in the presence of sparse response. In this manuscript, we
propose several Bayesian analyses for MMP data including a multivariate
analysis that models the latent covariance using a Bayesian FPCA. The Bayesian
FPCA gives a low-dimensional expansion using well-behaved basis functions, in
this case B-splines, that do not overly burden the multivariate analysis with
parameters, whereas simulating from a general latent correlation would add a
potentially large number of parameters to the model. Prior work demonstrates
the ability of the Bayesian FPCA to adequately model a variety of complex
covariance structures Goldsmith2016 ; MeyerMorris2022 . These structures are
not unique to functional data and the technique can be used for multivariate
data in general. The Bayesian FPCA is also useful for modeling latent
covariances in multivariate probit models since those structures may be
difficult to observe and therefore hard to pre-simplify with an assumption,
e.g. a compound symmetric covariance instead of an unstructured one.
Additionally, since the multivariate analysis is much more computationally
efficient than other multivariate probit models, it is a useful tool for
analyzing MMP data.

As we demonstrate in simulation, the multivariate analysis outperforms the
existing methods in the presence of sparse response regardless of the effect
size, as controlled by $\theta_{12}$, or the number of outcomes, $K$. It
routinely has coverage in excess of nominal and does not dip below 92% when
constructing 95% intervals. Conversely, the existing methods for MMP data have
inconsistent coverage that in some cases can be as low as 84% (again when
constructing 95% intervals). At the same time, the multivariate analysis has
the highest power without sacrificing bias or coverage and with only minimally
larger interval widths. Consequently, we recommend the use of the multivariate
analysis for MMP data with sparse response. This ensures that the multivariate
structure of the data is modeled, but we note that both the naïve and
penalized models perform similarly to the multivariate analysis. All Bayesian
methods we consider out perform the existing methods in the presence of sparse
responses.

The re-analysis of the SOC data confirms earlier findings, specifically that
contact with the education system is lower among primary care than specialty
care. But our multivariate analysis reveals several additional results that
warrant further investigation. Among the components that were previously
analyzed, contact with the mental health component may also be lower for
primary care than for specialty care. The main advantage our multivariate
analysis provides is in the ability to also assess the sparse components which
were not examined in the original univariate analysis. What we see is that
there is some evidence of increased contact with disability services among
primary care compared to specialty care. Taken all together, the SOC study and
our re-analysis should provide the foundation for further examination of the
coherence of SOC frameworks. This data can be thought of as a pilot study that
can provide a good source for developing future, larger studies to further
examine the differences in contact within the SOC between primary care and
specialty care. Our Bayesian multivariate analysis can seamlessly incorporate
prior information from this study by replacing the mean of the prior on
$\beta$ with values derived from the analysis presented in Table 3.

Supplementary Information

The Supplementary Materials referenced in Sections 3 and 4 can be found
alongside this manuscript.

Acknowledgments

Partial funding for this work was provided by internal Georgetown University
Summer Academic Research Grants.

Compliance with Ethical Standards

The authors have no potential or actualy conflicts of interest to declare. No
additional data was collected for this research. The original study was
retrospective in nature and was approved by the Boston University
Institutional Review Board Knutson2018 .

## References

  *   * (1) Knutson, K.H., Meyer, M.J., Thakrar, N., Stein, B.D.: Care coordination for youth with mental health disorders in primary care. Clinical Pediatrics 57, 5–10 (2018). https://doi.org/10.1177/0009922817733740
  * (2) Klingenberg, B., Agresti, A.: Multivariate extension of McNemar’s test. Biometrics 62, 921–928 (2006). https://doi.org/10.1111/j.1541-0420.2006.00525.x
  * (3) McNemar, Q.: Note on the sampling error of the difference between correlated proportions or percentages. Pyschometrika 12, 153–157 (1947). https://doi.org/10.1007/BF02295996
  * (4) Consonni, G., La Rocca, L.: Tests based on intrinsic priors for the equality of two correlated proportions. Journal of the American Statistical Association 103, 1260–1269 (2008). https://doi.org/10.1198/01621450800000043
  * (5) Saeki, H., Tango, T., Wang, J.: Statistical inference for noninferiority of difference in proportions of clustered matched-pair data from multiple raters. Journal of Biopharmaceutical Statistics 27, 70–83 (2017). https://doi.org/10.1080/10543406.2016.1148709
  * (6) Westfall, P.H., Troendle, J.F., Pennello, G.: Multiple McNemar tests. Biometrics 66, 1185–1191 (2010). https://doi.org/10.1111/j.1541-0420.2010.01408.x
  * (7) Xu, J., Yu, M.: Sample size determination and re-estimation for matched pair designs with multiple binary endpoints. Biometrical Journal 55, 430–443 (2013). https://doi.org/10.1002/bimj.201100231
  * (8) Lui, K.-J., Chang, K.-C.: Testing and estimation of proportion (or risk) ratio under the matched-pair design with multiple binary endpoints. Biometrical Journal 55, 603–616 (2013). https://doi.org/10.1002/bimj.201200224
  * (9) Lui, K.-J., Chang, K.-C.: Notes on testing noninferiority in multivariate binary data under the matched-pair design. Statistical Methods in Medical Research 25, 1272–1289 (2016). https://doi.org/10.1177/0962280213477022
  * (10) Cochran, W.G.: The comparison of percentages in matched samples. Biometrika 37, 256–266 (1950). https://doi.org/10.2307/2332378
  * (11) Mantel, N., Haenszel, W.: Statistical aspects of the analysis of data from retrospective studies of disease. Journal of the National Cancer Institute 22, 719–748 (1959). https://doi.org/10.1093/jnci/22.4.719
  * (12) Jiang, Y., Xu, J.: A comparative study of matched pair designs with two binary endpoints. Statistical Methods in Medical Research 26, 2526–2542 (2017). https://doi.org/10.1177/0962280215601136
  * (13) Agresti, A.: Categorical Data Analysis, 3${}^{\text{rd}}$ edn. John Wiley & Sons, Hoboken, NJ (2013)
  * (14) Altham, P.M.E.: The analysis of matched proportions. Biometrika 58, 561–576 (1971). https://doi.org/10.2307/2334391
  * (15) Broemeling, L.D., Gregurich, M.A.: A Bayesian alternative to the analysis of matched categorical responses. Communications in Statistics - Theory and Methods 25, 1429–1445 (1996). https://doi.org/10.1080/03610929608831777
  * (16) Ghosh, M., Chen, M.-H., Ghosh, A., Agresti, A.: Hierarchical Bayesian analysis of binary matched pairs data. Statistica Sinica 10, 647–657 (2000)
  * (17) Albert, J., Chib, S.: Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88, 669–679 (1993). https://doi.org/10.1080/01621459.1993.10476321
  * (18) Albert, J., Chib, S.: Bayesian residual analysis for binary response models. Biometrika 82, 747–759 (1995). https://doi.org/10.1093/biomet/82.4.747
  * (19) Gelman, A.: Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 1, 513–533 (2006). https://doi.org/10.1214/06-BA117A
  * (20) der Linde, A.V.: Variational Bayesian functional PCA. Computational Statistics & Data Analysis 53, 517–533 (2008). https://doi.org/10.1016/j.csda.2008.09.015
  * (21) Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B.: Bayesian Data Analysis, 3${}^{\text{rd}}$ edn. Chapman and Hall-CRC, Boca Raton, FL (2013)
  * (22) Polson, N.G., Scott, J.G.: On the half-Cauchy prior for a global scale parameter. Bayesian Analysis 7, 887–902 (2012). https://doi.org/10.1214/12-BA730
  * (23) Wand, M.P., Ormerod, J.T., Padoan, S..A., Frühwirth, R.: Mean field variational bayes for elaborate distributions. Bayesian Analysis 6, 847–900 (2011). https://doi.org/10.1214/11-BA631
  * (24) Chib, S., Greenberg, E.: Analysis of multivariate probit models. Biometrika 85, 347–361 (1998). https://doi.org/10.1093/biomet/85.2.347
  * (25) Liu, C.: Discussion. Journal of Computational and Graphical Statistics 10, 75–81 (2001). https://doi.org/10.1198/10618600152418746
  * (26) Zhang, X., Boscardin, W.J., Belin, T.R.: Sampling correlation matrices in Bayesian models with correlated latent variables. Journal of Computational and Graphical Statistics 15, 880–896 (2006). https://doi.org/10.1198/106186006X160050
  * (27) Webb, E.L., Forster, J.J.: Bayesian model determination for multivariate ordinal and binary data. Computational Statistics & Data Analysis 52, 2632–2649 (2008). https://doi.org/10.1016/j.csda.2007.09.008
  * (28) Goldsmith, J., Kitago, T.: Assessing systematic effects of stroke on motor control using hierarchical function-on-scalar regression. Journal of the Royal Statistical Society, Series C 65, 215–236 (2016). https://doi.org/%****␣bmmp.bbl␣Line␣450␣****10.1111/rssc.12115
  * (29) Meyer, M.J., Morris, J.S., Gazes, R.P., Coull, B.A.: Ordinal probit functional outcome regression with application to computer-use behavior in rhesus monkeys. Annals of Applied Statistics 16, 537–550 (2022). https://doi.org/10.1214/21-AOAS1513
  * (30) Agresti, A., Coull, B.A.: Approximate is better than ’exact’ for interval estimation of binomial proportions. The American Statistician 52, 119–126 (1998). https://doi.org/10.1080/00031305.1998.10480550
  * (31) Agresti, A., Caffo, B.: Simple and effective confidence intervals for proportions and differences of proportions result from adding two successes and two failures. The American Statistician 54, 280–288 (2000). https://doi.org/10.1080/00031305.2000.10474560
  * (32) Gelman, A., Rubin, D.B.: Inference from iterative simulation using multiple sequences. Statistical Science 7, 457–511 (1992). https://doi.org/10.1214/ss/1177011136
  * (33) Costello, E.J., He, J.-P., Sampson, N.A., Kessler, R.C., Merikangas, K.R.: Services for adolescents with psychiatric disorders: 12-month data from the National Comorbidity Survey-Adolescent. Psychiatric Services 65, 359–366 (2014). https://doi.org/10.1176/appi.ps.201100518