Bayesian nonparametric multiway regression for clustered binomial data

TL;DR

This paper presents a Bayesian nonparametric multiway regression model for clustered binomial data, effectively capturing complex structures and dependencies in multi-dimensional datasets, demonstrated on periodontal disease data.

Contribution

It introduces a novel nonparametric regression approach with covariate-dependent clustering for multiway binomial data, utilizing low-rank assumptions for dimensionality reduction.

Findings

Model outperforms competitors in periodontal disease data analysis.

Flexible probit stick-breaking formulation captures covariate dependence.

Efficient Gibbs sampling enables practical posterior computation.

Abstract

We introduce a Bayesian nonparametric regression model for data with multiway (tensor) structure, motivated by an application to periodontal disease (PD) data. Our outcome is the number of diseased sites measured over four different tooth types for each subject, with subject-specific covariates available as predictors. The outcomes are not well-characterized by simple parametric models, so we use a nonparametric approach with a binomial likelihood wherein the latent probabilities are drawn from a mixture with an arbitrary number of components, analogous to a Dirichlet Process (DP). We use a flexible probit stick-breaking formulation for the component weights that allows for covariate dependence and clustering structure in the outcomes. The parameter space for this model is large and multiway: patients $\times$ tooth types $\times$ covariates $\times$ components. We reduce its effective dimensionality, and account for the multiway structure, via low-rank assumptions. We illustrate how this can improve performance, and simplify interpretation, while still providing sufficient flexibility. We describe a general and efficient Gibbs sampling algorithm for posterior computation. The resulting fit to the PD data outperforms competitors, and is interpretable and well-calibrated. An interactive visual of the predictive model is available at http://ericfrazerlock.com/toothdata/ToothDisplay.html , and the code is available at https://github.com/lockEF/NonparametricMultiway .

Figures

Tables

Table 1. Table 1: Mean LPPL for logistic simulation under different modeling approaches. For the logistic models, separate corresponds to independent coefficients and errors for each tooth type, shared β 𝛽 \beta corresponds to shared coefficients for each tooth type, and shared β , ϵ 𝛽 italic-ϵ \beta,\epsilon corresponds to shared errors and coefficients for each tooth type.

Stick-breaking	No $𝔼$	Rank($𝔼$) = 1	Rank($𝔼$) = 2
No $𝔹$	-3958	-3524	-3511
Rank($𝔹$)=1	-3477	-3135	-3176
Rank($𝔹$)=2	-3496	-3156	-3201
Full $𝔹$	-3625	-3324	-3393
Logistic	Separate	Shared $\beta$	Shared $\beta$, $ϵ$
	-3674	-3628	-2945

Table 2. Table 2: Mean LPPL for low-rank probit simulation under different modeling approaches. For the logistic models, separate corresponds to independent coefficients and errors for each tooth type, shared β 𝛽 \beta corresponds to shared coefficients for each tooth type, and shared β , ϵ 𝛽 italic-ϵ \beta,\epsilon corresponds to shared errors and coefficients for each tooth type.

Stick-breaking	No $𝔼$	Rank($𝔼$) = 1	Rank($𝔼$) = 2
No $𝔹$	-3965	-3981	-3979
Rank($𝔹$)=1	-3277	-3301	-3325
Rank($𝔹$)=2	-3303	-3390	-3376
Full $𝔹$	-3456	-3477	-3465
Logistic	Separate	Shared $\beta$	Shared $\beta$, $ϵ$
	-4016	-4214	-4405

Table 3. Table 3: Mean LPPL for full probit simulation under different modeling approaches. For the logistic models, separate corresponds to independent coefficients and errors for each tooth type, shared β 𝛽 \beta corresponds to shared coefficients for each tooth type, and shared β , ϵ 𝛽 italic-ϵ \beta,\epsilon corresponds to shared errors and coefficients for each tooth type.

Stick-breaking	no $𝔼$	Rank($Z$) = 1	Rank($Z$) = 2
No $𝔹$	-3961	-3975	-3981
Rank($𝔹$)=1	-3401	-3389	-3395
Rank($𝔹$)=2	-3123	-3147	-3167
Full $𝔹$	-3105	-3138	-3145
Logistic	Separate	Shared $\beta$	Shared $\beta$, $ϵ$
	-4005	-4341	-4583

Table 4. Table 4: Mean LPPL for low-rank probit DP models.

	No $𝔼$	Rank($𝔼$) = 1	Rank($𝔼$) = 2
Rank($𝔹$)=1	-3953	-3625	-3645
Rank($𝔹$)=2	-3946	-3647	-3644

Equations

$$Y_{ij}\sim \text{mbox}Binomial(n_{ij},p_{ij})$$

$$p_{ij}=\sum _{h=1}^{\infty }{\pi }_h{\delta }_{{\theta }_h},$$

$${\pi }_h=V_h\prod _{l<h}(1-V_l),$$

$$V_h=\Phi (Z_h),$$

$$V_{ijh}=\Phi \left(Z_{jh}+X[i,:]\mathbb{B}[:,h]\right).$$

$$V_{ijh}=\Phi \left(Z_{jh}+X[i,:]\mathbb{B}[:,j,h]\right),$$

$$\mathbb{B}[:,j,h]\stackrel{iid}{\sim }N(\bar{\mathbb{B}}[:,h],\Sigma )\text{ for }j=1,\dots ,J,$$

$$\begin{array}{rl}\mathbb{B}=& [[B_1,B_2,B_3]]\\ \to \mathbb{B}[d,j,h]=& \sum _{r=1}^R{B}_1[d,r]B_2[j,r]B_3[h,r]\end{array}$$

$$V_{ijh}=\Phi \left(Z_{jh}+X[i,:]\mathbb{B}[:,j,h]+{\mathbb{E}}_{ijh}\right).$$

$$\mathbb{E}[i,:,h]\sim N(0,{\Sigma }_h).$$

$$\mathbb{E}=$$

$$Z_{ijh}^{\ast }\sim \text{mbox}Normal\left(Z_{jh}+X[i,:]\mathbb{B}[:,j,h]+{\mathbb{E}}_{ijh},1\right).$$

$$\{C_{ij}:i,j\},\{Z_{ijh}^{\ast }:i,j,h\},\{{\theta }_h:h\},\{Z_{jh}:j,h\},B_1,B_2,B_3,E_1,E_2,E_3,\alpha ,\text{ and }{\sigma }^2.$$

$${\tilde{\mathbb{E}}}^{(t)}$$

$${\tilde{V}}_{jh}^{(t)}$$

$${\tilde{\pi }}_{jh}^{(t)}$$

$$p\left(Y_{i,:}\mid Y_{[i],:}\right)\approx \frac{1}{T}\sum _{t=1}^T\prod _{j=1}^J\left(\genfrac{}{}{0ex}{}{n_{ij}}{Y_{ij}}\right)({\tilde{p}}_j^{(t)}{)}^{Y_{ij}}(1-{\tilde{p}}_j^{(t)}{)}^{n_{ij}-Y_{ij}}.$$

$$\text{mbox}LPPL=\sum _{i=1}^n\text{mbox}logp\left(Y_{i,:}\mid Y_{[i],:}\right).$$

$$\text{mbox}log\left(\frac{p_{ij}}{1-p_{ij}}\right)=X_i\beta +{ϵ}_i,$$

$$\text{mbox}log\left(\frac{p_{ij}}{1-p_{ij}}\right)=X_i{\beta }_j+{ϵ}_{ij}.$$

$$P(Y_{ij}\mid X)=P(Y_{ij}\mid X_i\beta )\forall i,.$$

$$\text{mbox}log\left(\frac{p_{ij}}{1-p_{ij}}\right)=X_{ij}{\beta }_j,$$

$$\text{mbox}log\left(\frac{p_{ij}}{1-p_{ij}}\right)=X_{ij}{\beta }_j+{ϵ}_{ij},$$

$${\mathbf{Y}}_i^{(1)},\dots ,{\mathbf{Y}}_i^{(1000)}\stackrel{iid}{\sim }p\left({\mathbf{Y}}_i\mid {\mathbf{X}}_i,\{{\mathbf{Y}}_j,{\mathbf{X}}_j{\}}_{j\ne i}\right).$$

$$\hat{P}(\varphi <x)=\frac{1}{290\ast 4}\sum _{i=1}^{290}\sum _{j=1}^4\text{mathbbm}{1}_{\{{\varphi }_{ij}<x\}}.$$

$$P\left(C_{ij}=h\mid {\theta }_h,Y_{ij}\right)\propto {\pi }_{ijh}{\theta }_h^{Y_{ij}}(1-{\theta }_h{)}^{n_{ij}-Y_{ij}}.$$

$${\tilde{Y}}_h=\sum _{\{i,j:C_{ij}=h\}}{Y}_{ij},{\tilde{n}}_h=\sum _{\{i,j:C_{ij}=h\}}{n}_{ij}.$$

$$\text{mbox}Beta\left(a+{\tilde{Y}}_h,b+{\tilde{n}}_h-{\tilde{Y}}_h\right).$$

$$Z_{ijh}^{\ast }\sim \{\begin{array}{l}N(Z_{jh}+X[i,:]\mathbb{B}[:,j,h]+{\mathbb{E}}_{ijh},1)\text{mathbbm}{1}_{{\mathbb{R}}^{-}}\text{ if }h<C_{ij}\\ N(Z_{jh}+X_i{\beta }_{jh}+{\mathbb{E}}_{ijh},1)\text{mathbbm}{1}_{{\mathbb{R}}^{+}}\text{ if }h=C_{ij}\end{array}$$

$$E_{23}[:,r]=\text{mbox}vec\left(E_2[:,r]\otimes E_3[:,r]\right),$$

$$Z_E^{\ast \text{mat}}[i,:]=\text{mbox}vec({\mathbb{Z}}_E^{\ast }[i,:,:]).$$

Full text

Key words: Bayesian nonparametrics; binomial regression; Dirichlet processes; parafac/candecomp; tensor factorization

1 Introduction

Clustered, or longitudinal count data are ubiquitous in medicine and public health. Here, count responses, i.e., bio-markers for the corresponding disease, are collected in clinical studies/trials for each subject at multiple locations for the same subject (leading to clustered data), and/or at multiple occasions (time points), leading to longitudinal data. Statistical inference in this setup is typically considered under a generalized linear (or non-linear) mixed model specification, with appropriate distributional assumptions for the count responses (such as Poisson, negative binomial, etc), and the random (effect and error) terms. Our motivating data example is from a clinical study (Fernandes et al., 2009) of periodontal disease (PD) among Type-2 diabetic Gullah speaking African-Americans (henceforth, GAAD data), conducted at the Medical University of South Carolina. Here, the clinical attachment level (or, CAL), a clinical marker of PD, was recorded at each of the six sites of a subject’s tooth producing a clustered (multivariate) setup, and we are interested in assessing the covariate-response relationships on the ‘number of diseased, or missing tooth sites specific to a tooth-type’, such as incisors, canines, premolars, and molars. With the total count of tooth-sites $n_{ij}$ fixed for the $j$th tooth-type $(j=1,\dots,4)$ for a subject, and the simplifying assumption that the probability of a diseased, or missing tooth site is the same within each tooth-type, this reduces to the classical binomial regression problem for clustered data, where the probability of disease/missingness can be explained by the covariates via an appropriate link-function (say, probit, logit, etc). Related estimation and inference can be easily carried out under a variety of classical and Bayesian framework using popular software, such as SAS, or R.

However, a careful observation of the (raw) count histogram (see Figure 1) reveals the inadequacy of a standard maximum-likelihood (ML) Binomial fit (without any covariates). The bimodal behavior of the fitted Binomial model results from $n_{ij}$ taking two values (in particular, 48 for molars, premolars and incisors; 24 for canines), and a point estimate of the binomial probability $p_{ij}$ fails to explain the rugged histogram, including the bumps at the two extremities. It is likely that an estimated $p_{ij}$ distributed as a mixture of densities in $(0,1)$ may have explained the response histogram better. However, the (automatic) determination of the number of mixture components in a finite mixture model (FMM) is often challenging, mostly due to over- and under-fitting issues (Geoffrey and Peel, 2000). Our focus is Bayesian, and most Bayesian implementation relies on popular model choice criteria, such as the Bayesian information criterion (BIC), pseudolikelihood information criterion (PLIC), etc, rendering a computationally demanding framework for searching across the cardinality spectrum of the mixture components. To alleviate this, increasingly popular Bayesian nonparametric (BNP) methods are available, where the basic idea is to fit a single model that can adapt to the observed data complexity, instead of fitting and comparing different models (Lai et al., 2018). The most popular BNP model is the Dirichlet Process (DP) mixture model (Antoniak, 1974), which is an infinite mixture model with a base distribution, such as the normal (MacEachern and Müller, 1998), Student’s $t$ (Wei and Li, 2012), hidden Markov models (Wei and Li, 2011), etc. The DP can be considered a special case of the class of stick-breaking priors (Ishwaran and James, 2001), which define a more general family of infinite mixture models.

We consider a stick-breaking prior for the binomial probabilities $p_{ij}$, with a uniform base distribution. The DP for binomial probabilities has been considered previously (Gopalan and Berry, 1998; Lock and Dunson, 2017); in our context, it provides a much more flexible support for the underlying distribution, and a much better fit to the observed counts (see Figure 1). However, our primary interest is the relation between patient-level covariates (age, smoking status, etc.) and the proportion of diseased sites across different tooth types for each patient. There is a large body of work on covariate-dependent modeling using Bayesian nonparametrics (MacEachern, 2000; Griffin and Steel, 2006). In particular, the probit stick-breaking (PSB) formulation (Chung and Dunson, 2009; Rodriguez and Dunson, 2011) is attractive for our context, because it admits a general framework for conditional dependence on covariates and intra-individual dependence, where the conditional distribution for the outcome is still an infinite mixture. The covariates determine the mixture weights via a probit model for each stick-breaking component. The full PSB representation involves separate coefficients for each triplet of (a) stick-breaking mixture component, (b) tooth type, and (c) covariate. This richly parametrized model is very flexible, but is inefficient if the coefficients are treated independently. To explore a more general dependency structure among the coefficients indexed by (a)–(c) facilitated by the Bayesian proposal that automates seamless exchange of information, we cast this into a multiway regression framework (Lock, 2018) that treats the coefficients as a three-way array, and use the CANDECOMP/PARAFAC (CP) low-rank tensor factorization (Kolda and Bader, 2009) for estimation. We use a similar approach to model the intra-individual correlation across components and tooth types. The use of low-rank multiway (i.e., tensor) modeling within this novel BNP context is our primary methodological contribution, which can be easily generalized to other scenarios.

In Section 2 we describe our model development, beginning with the straightforward DP model for a binomial outcome. In Section 3, we formulate an efficient Gibbs sampling approach to posterior inference. In Section 4, we describe three simulation studies to validate our model and illustrate its motivation. In Section 5, we present the application to the GAAD data. Finally, in Section 6, we make some concluding remarks and discuss extensions of our approach. Details on the Gibbs sampling approach are relegated to the Appendix.

2 Statistical Model

2.1 Binomial likelihood

Let $Y_{ij}$ denote the proportion of diseased sites, and let $n_{ij}$ denote the total number of sites measured, for subject $i$ $(i=1,\ldots,I)$ in tooth type $j$ $(j=1,\ldots,J)$. Then,

[TABLE]

for some latent binomial probability $p_{ij}$. In what follows, we describe our approach to modeling the $p_{ij}$’s, beginning with a straightforward PSB model in Section 2.2 that we extend to allow for multiway structure in Section 2.3.

2.2 Probit stick-breaking formulation

We proceed with our BNP model development for the without- and with-covariate cases. Without covariates, we assume the $p_{ij}$’s have a DP prior with a beta base distribution. That is, $p_{ij}\sim P$ where $P\sim DP(\mbox{Beta}(a,b),\alpha)$. Then, each $p_{ij}$ is drawn from a theoretically infinite number of realizations $\theta_{h}$ from Beta$(a,b)$, with corresponding probability weights $\pi_{h}$,

[TABLE]

where $\delta_{\theta_{h}}$ is a point mass at $\theta_{h}$. The distribution for the weights $\pi_{h}$ can be represented as a stick-breaking process (Sethuraman, 1994), with the $\pi_{h}$ generated as

[TABLE]

where $V_{h}\overset{iid}{\sim}\mbox{Beta}(1,\alpha)$. In what follows we use a uniform base distribution ($a=b=1$), and set the concentration parameter $\alpha$ to $1$ unless otherwise specified.

Now, we extend this framework to include covariates. Note that when $\alpha=1$, the distribution of the “stick-breaks" $V_{h}$ (2) is equivalently given by

[TABLE]

where $Z_{h}\overset{iid}{\sim}N(0,1)$ and $\Phi$ is the standard normal CDF. Given individual-level covariates $X:I\times D$, we extend (3) to allow the stick-breaks to depend on $X$ via coefficients $\mathbb{B}$, where the $h$’th column of $\mathbb{B}$ gives the coefficients for component $h$ of the stick-breaking process. That is, if $V_{ijh}$ gives the stick-break for individual $i$, tooth type $j$ and component $h$,

[TABLE]

where, $Z_{jh}$ is an intercept term for tooth type $j$ and component $h$, $X[i,:]$ is the $i$’th row of $X$ and $\mathbb{B}[:,h]$ is the $h$’th column of $\mathbb{B}$. This formulation closely mimics the general approach to covariate-dependent DPs of Chung and Dunson (2009), and Rodriguez and Dunson (2011).

For our prior choices, we use independent $N(0,1)$ priors on the coefficients $\mathbb{B}$. We also assume $Z_{jh}\overset{iid}{\sim}N(\alpha,1)$, where $\alpha\sim N(0,1)$ is a concentration parameter, such that smaller $\alpha$ favors a larger number of components with non-trivial probability under the stick-breaking process.

2.3 Multiway extensions

Note that model (5) in Section 2.2 involves three simplifying assumptions:

The effect of the individual-level covariates (age, gender, etc.), $\mathbb{B}[:,h]$, are consistent across tooth types $j$ for each component $h$. 2. 2.

The coefficients $\mathbb{B}$ are independent across components: \mathbb{B}[:,1]\rotatebox[origin={c}]{90.0}{\models}\mathbb{B}[:,2]\rotatebox[origin={c}]{90.0}{\models}\mathbb{B}[:,3]\cdots. 3. 3.

The observations for each individual are independent across types $j$, given covariates $X[i,:]$.

Assumption 1. can be relaxed with the incorporation of type-specific coefficients. We extend model (4) as follows:

[TABLE]

where $\mathbb{B}[:,j,h]$ gives type-specific coefficients for each component $h$. One ad-hoc approach to (5) is to give the entries of $\mathbb{B}$ independent $N(0,1)$ priors, i.e., model each type entirely independently. Another approach is to use a standard hierarchical modeling framework across types for each component, e.g.,

[TABLE]

with covariance matrix $\Sigma$. This approach assumes that covariate effects are similar across types for a given component $h$. However, the same component may have different implications by type; for example, the same atom $\theta_{h}$ may give a relatively high value for one type and a relatively low value for another type. Moreover, the assumption of independent effects across the components (assumption 2. above) is almost certainly not satisfied. For example, within a type covariate effects may be similar for components whose atoms are close together ($\theta_{h}\approx\theta_{h^{\prime}})$.

With this motivation, we consider a more general model for the dependence structure of the covariate effects $\mathbb{B}$. It helps to consider the full set of coefficients as the three-way array $\mathbb{B}:D\times J\times H$: where $D$ is the number of individual-level covariates, $J$ is the number of types, and $H$ is the number of stick-breaking components (without truncation, $H=\infty$). Under this framework, the CP representation (Carroll and Chang, 1970) restricts $\mathbb{B}$ to be of low CP-rank: rank$(\mathbb{B})=R$. That is, $\mathbb{B}$ has the form of a CP factorization:

[TABLE]

where $B_{1}:D\times R$, $B_{2}:J\times R$, and $B_{3}:H\times R$ give $R$-dimensional representations for the covariates, tooth types and components, respectively. Note that the CP factorization simply extends the traditional low-rank matrix factorization (e.g., via a singular value decomposition) to higher-order arrays. This approach has the dual advantage of accounting for the potential multiway dependence of the coefficients, and reducing the dimensionality of the model from $JDH$ to $R(J+D+H)$.

There is a fast-growing literature on Bayesian approaches to CP factorization for multiway coefficient arrays in regression problems (Guhaniyogi et al., 2017; Zhou et al., 2015; Johndrow et al., 2017). These approaches may be extended to our context. In particular, we put a $N(0,1)$ prior on the marginal entries of $\mathbb{B}$, as before, but have their joint pdf supported only on the space of arrays with CP-rank $\leq$ R (see Section 8 of Lock (2018)). This prior facilitates Gibbs sampling, and can easily be incorporated into the sampling scheme of Section 3 with additional steps to sample the terms of the factorization for $\bm{\mathbb{B}}$.

Assumption 3. above, on the conditional independence of outcomes across tooth types, can be relaxed by incorporating correlated individual-level errors $\mathbb{E}_{ijh}$:

[TABLE]

These error terms can be represented as a three-way array of dimension $\mathbb{E}:I\times J\times H$. A straightforward ad-hoc approach is to model the $\mathbb{E}_{ijh}$’s separately for each component $h$, with a $J\times J$ covariance matrix $\Sigma_{h}$ for correlation between tooth types:

[TABLE]

However, this model is limiting, as it does not account for the likely correlations across components $h=1,2,\ldots$, and assumes that the correlation between tooth types is consistent across components. Instead, as for $\mathbb{B}$, the low-rank CP factorization yields a flexible and parsimonious model for multiway dependence in $\mathbb{E}$. That is,

[TABLE]

where $E_{1}:I\times R_{e}$, $E_{2}:J\times R_{e}$, and $E_{3}:H\times R_{e}$. Thus, $R_{e}$ latent random effects for each subject are weighted across the tooth types via $E_{2}$, and across the components via $E_{3}$. To facilitate sampling from the posterior predictive distribution for a new subject, we consider a separable prior on each of $E_{1}$, $E_{2}$ and $E_{3}$, where the entries of $E_{2}$ and $E_{3}$ are independent $N(0,1)$, and the entries of  the $r$’th column of $E_{1}$ (corresponding to the $r$’th latent effect) are independent $N(0,\sigma^{2}_{r})$. Because the scale of $E_{1}$, $E_{2}$ and $E_{3}$ are not independently identifiable, $\sigma^{2}_{r}$ controls the overall scale for the $r$’th rank-1 effect, allowing for appropriate shrinkage of the individual effects. We use diffuse inverse-gamma priors IG$(0.1,0.1)$ for each $\sigma^{2}_{r}$. A diagram of the full hierarchical model is shown in Figure 2.

3 Posterior Computation and Inference

3.1 Augmented Gibbs sampling algorithm

The model described in Section 2 and illustrated in Figure 2 allows for a straightforward and efficient augmented Gibbs sampling scheme for posterior computation. Let $C_{ij}$ give the index of the component $h$ that generated $Y_{ij}$: $C_{ij}\in\{1,\ldots,H\}$, $p_{ij}=\theta_{C_{ij}}$. Extending the augmentation approach in Albert and Chib (1993), we use latent Gaussian random variables to represent the probit link:

[TABLE]

where $C_{ij}=h$ if and only if $Z^{*}_{ijh}>0$ and $Z^{*}_{ijh}<0$ for all $k<h$. Gibbs sampling proceeds by sampling from the conditional distributions for

[TABLE]

Full details of the sampling algorithm are given in the Appendix A.

3.2 Posterior predictive simulation

Given a new subject with covariates $\tilde{x}:1\times p$, we sample from their joint posterior predictive distribution via the generative model at each Gibbs sampling iteration. That is, if $(\cdot)^{(t)}$ denotes the $t$’th draw for a given parameter in the Gibbs sampling chain, we draw from the posterior predictive distribution for the subject outcomes over the $J$ types $(\tilde{Y}_{1},\ldots,\tilde{Y}_{J})$ as follows:

Generate subject latent effects $\tilde{E}_{1}^{(t)}=[\tilde{E}_{11}^{(t)}\,\tilde{E}_{12}^{(t)}\,\ldots,\tilde{E}_{1R_{e}}^{(t)}]$ where $\tilde{E}_{1r}^{(t)}\overset{\text{indep}}{\sim}N(0,\sigma^{2}_{r})$. 2. 2.

Determine component weights $\tilde{\pi}_{jh}^{(t)}$ for $j=1,\ldots,J$, $h=1,\ldots,H$ via

[TABLE] 3. 3.

Generate $\tilde{p}_{j}^{(t)}$ from $\sum_{h=1}^{H}\tilde{\pi}_{jh}^{(t)}\delta_{\theta_{h}^{(t)}}$ for $j=1,\ldots,J$ 4. 4.

Generate $\tilde{Y}_{j}^{(t)}$ from $\mbox{Binomial}(\tilde{n}_{j},\tilde{p}_{j}^{(t)})$ for $j=1,\ldots,J$.

3.3 Rank selection and evaluation

To select the ranks of the model, $R$ and $R_{e}$, we consider the posterior predictive density of the subjects under a cross-validation scheme. That is, we consider $p\left(Y_{i,:}\mid Y_{[i],:}\right)$, where $Y_{i,:}=(Y_{i1},\ldots,Y_{iJ})$ and $Y_{[i],:}$ is the data for all subjects but subject $i$. We generate samples $\tilde{p}_{j}^{(t)}$ from the posterior predictive distribution for $Y_{i,:}$ given $Y_{[i],:}$ as in 3.2, and then approximate the predictive likelihood as follows:

[TABLE]

To evaluate performance and select the ranks, we favor a larger log posterior predictive likelihood (LPPL):

[TABLE]

4 Simulations

Here, we describe three simulation studies, that are intended to illustrate the low-rank PSB regression model and demonstrate its flexibility. For all studies, the generated data matched the observed GAAD data in terms of number of patient ($N=290$), number of tooth types ($M=4$), and number of patient-level covariates ($p=6$). The first three columns of the covariate matrix $X:I\times D$ are generated independently from a $N(0,1)$ distribution, while the next three columns are generated as categorical indicators independently from a Bernoulli$(0.5)$ distribution; following which, each column is scaled to have mean [math] and variance $1$. The outcomes are generated by a binomial distribution $Y_{ij}\sim\mbox{Binomial}(n_{ij},p_{ij})$, where $n_{ij}$ is the number of sites observed for subject $i$, and tooth type $j$ in the GAAD data, and $p_{ij}$ is generated from an assumed model. In Section 4.1, $p_{ij}$ are given by a simple logistic model, in Section 4.2 by a low-rank PSB regression model, while in Section 4.3 by a PSB regression model where the coefficients are not low-rank.

4.1 Logistic simulation

We generate latent probabilities $p_{ij}$ according to the simple logistic regression model

[TABLE]

where $\beta=(0,1,-1,0,1,-1)$ and $\epsilon_{i}\overset{iid}{\sim}N(0,1)$. Note that the covariate effects are consistent across tooth types, and the latent error term $\epsilon_{i}$ is at the subject level. Thus, the underlying probabilities for a subject are the same across types, $p_{i1}=p_{i2}=p_{i3}=p_{i4}$, though the number of affected sites $Y_{ij}$ may differ.

Table 1 shows the LPPL (8) under various PSB models and logistic models for $p_{ij}$. Our logistic model for inference matches the likelihood in (9), with diffuse independent $N(0,100)$ priors for the coefficients $\beta$, and $\epsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2})$ with $\sigma^{2}\sim\mbox{IG}(0.1,0,1)$. We also consider analogous logistic models with independent errors ($\epsilon_{ij}$) or covariates ($\beta_{j}$) across tooth type:

[TABLE]

As expected, simple logistic regression models generally tend to outperform nonparametric approaches; in particular, the logistic model with shared coefficients and errors across types matches the generative model (9) precisely, and performs the best. However, low-rank probit stick-breaking models are competitive. The best performing stick-breaking model has a rank-$1$ coefficient array $\mathbb{B}$ and a rank-$1$ error array $\mathbb{E}$.

Figure 3 displays elements of the fit for the best performing PSB model when rank$(\mathbb{B})=1$ and rank$(\mathbb{E})=1$, illustrating the motivation for a low-rank approach. The loadings for each covariate, $B_{1}$ in (6), are closely proportional to their values in the logistic model $\beta=(0,1,-1,0,1,-0)$. The loadings for each tooth type are roughly equivalent, which makes sense, given that the covariate effects are equal across tooth types. Under the true generative model $X\beta$ suffices to describe the conditional distribution of $Y$ given $X$,

[TABLE]

Thus, the better performance of a rank-$1$ restriction on $\mathbb{B}$ is explained by the appropriate reduction in dimensionality over having a separate set of coefficients for each component of the stick-breaking process and each tooth type. A rank-$1$ model on the error terms is also appropriate, because the unidimensional shared error term $\epsilon_{i}$ suffices to describe the dependence of $Y_{i}$ across tooth types.

4.2 Low-rank stick breaking

Here, we generate data under a rank-$1$ stick breaking model. That is, we generate data according to the stick breaking model defined by Equations (1) and (2), where the $\theta_{h}$ are generated from a Beta$(1,1)$ distribution. The stick breaks $V_{h}$ are generated as in the probit formulation (5), and coefficient array $\mathbb{B}$ has a rank-$1$ multiway structure (6) where the entries of the loading vectors $B_{1},B_{2},B_{3}$ are generated independently from a $N(0,1)$ distribution. Thus, under the generative model, a rank$-1$ representation for $\mathbb{B}$ is correct, and the tooth types are conditionally independent given the covariates $X$.

Table 3 shows the LPPL (8) under different models for $p_{ij}$. The model with Rank($\mathbb{B}$)=1 and no $\mathbb{E}$ closely matches the true generative model, and performs the best. Note that the logistic models perform much worse in this context because the underlying relationship between $p_{ij}$ and $X$ does not have a logistic form, and is in general not well described by any monotone link. Moreover, even though the generative model is rank-$1$, the effect of $X$ on $p_{ij}$ may differ by tooth type via different type loadings $B_{2}$. Hence, the logistic models with shared coefficients across tooth types perform even worse.

4.3 Full-rank stick breaking

Here, we generate data under a PSB model, wherein, the coefficients have no shared structure across types or components. That is, we generate data as in Section 4.2, except the coefficient array that describes the stick-breaking process, $\mathbb{B}$, is composed entirely of independent $N(0,1)$ entries. Thus, the underlying model has no low-rank structure, and besides the shared atoms $\theta_{h}$, there is no dependence among tooth types within a subject.

Table 3 shows the LPPL (8) under different models for $p_{ij}$. The model with full $\mathbb{B}$ and no $\mathbb{E}$ closely matches the true generative model, and performs the best. However, a rank-$2$ model has comparable performance, and this demonstrates how a model with restricted rank can reasonably approximate a more complex nonparametric model.

5 Application: GAAD Data

5.1 Methods comparison

We compare the approaches described in Sections 2.2-2.3, with alternative parametric approaches, for modeling the GAAD data. As a parametric approach, we consider the logistic regression model

[TABLE]

with diffuse $N(0,100)$ priors on the $\beta_{j}$. We also consider a model with an additional error term for unobserved factors:

[TABLE]

where $\epsilon_{ij}\overset{iid}{\sim}N(0,\sigma^{2})$ and $\sigma^{2}\sim\mbox{IG}(0.1,0.1)$. We consider the following approaches, and compute the LPPL for each as in (8):

Logistic model, as in Equation (10). [LPPL:-10084] 2. 2.

Logistic model with error, as in Equation (11). [LPPL:-4176] 3. 3.

Marginal DP, as described in Section 2.2. [LPPL: -4222] 4. 4.

PSB model with equal type effects, as described in Section 2.2. [LPPL: -4028] 5. 5.

PSB model with separate type effects, as in Equation (5). [LPPL: -4005] 6. 6.

Probit DP with multiway $\mathbb{B}$ and individual effects $\mathbb{E}$, as in Equation (7): See Table 4

Note that the best performing model in terms of LPPL is the probit DP model with rank$(\mathbb{B})=1$ and rank$(\mathbb{E})=1$. Thus, there are clear advantages behind considering a BNP approach, and also advantages to simplifying the dimensionality of the model via low-rank constraints. We expand on the results for this model in Sections 5.2 and 5.3.

5.2 Validation

We perform a 10-fold cross-validation study of the PSB model with rank$(\mathbb{B})=1$ and rank$(\mathbb{E})=1$ to assess the calibration and coverage of the predictive model. We randomly partition the data for $N=290$ individuals into $10$ equal subgroups of size $29$, $\{\mathbf{Y}_{1},\mathbf{X}_{1},\ldots,\mathbf{Y}_{10},\mathbf{X}_{10}\}$. For each subgroup, we draw $1000$ samples from the posterior predictive distribution of $\mathbf{Y}_{i}$, given their individual-level covariates and data for the remaining subgroups:

[TABLE]

Figure 4 shows a histogram of the proportion of diseased sites for all out-of-sample posterior predictive simulations and for the observed proportions. The marginal distribution of the simulated values closely matches that for the observed values. To assess the calibration of posterior predictive distributions and their inferential accuracy, we consider the quantile of each observed value under its out-of-sample posterior predictive distribution $\phi_{ij}=\frac{1}{1000}\sum_{t=1}^{1000}\mathbbm{1}_{\{Y_{ij}<Y_{ij}^{(t)}\}}$, and compute the empirical CDF of the posterior predictive quantiles, as follows:

[TABLE]

Figure 4 shows a plot of the empirical CDF of the posterior predictive quantiles. This is a straight line approximating the identity function, demonstrating that the predictive model is well-calibrated, and does not under-estimate or over-estimate uncertainty.

5.3 Interpretation

Figure 5 summarizes the posterior for the PSB model with rank$(\mathbb{B})=1$ and rank$(\mathbb{E})=1$, conditioning on the full dataset with $N=290$ individuals. Because rank$(\mathbb{B})=1$, the covariate effects can be decomposed into a single vector of weights for each covariate, for each tooth type, and for each component of the stick-breaking process. Figure 5 shows the weights for covariates and tooth types, with $95\%$ credible intervals. The covariate weights show a strong positive effect of age and a weak positive effect of smoking, a strong negative effect of sex (F) and a weak negative effect of BMI, and a negligible effect of A1C. The tooth type weights are all consistently positive, demonstrating that the effect of the covariates are similar across types, but may vary in magnitude. The conditional distribution of diseased sites across tooth types for an individual depends entirely on their univariate covariate score $XB_{1}$. Figure 5 illustrates the marginal distribution across tooth types for certain values on this continuum; for a higher covariate score, the resulting distribution favors a higher proportion of diseased sites. Thus, status is improved for females, non-smokers, younger individuals, and those with a higher BMI, and these effects are generally consistent across tooth types. The findings for sex, smoking and age are consistent with the literature on periodontal disease (Chambrone et al., 2010). The finding for BMI is somewhat counter-intuitive, and results are mixed in the literature (Kongstad et al., 2009; Chaffee and Weston, 2010); an ad-hoc correlation analysis of these data similarly yields a negative correlation between BMI and proportion of diseased sites for all tooth types, with p-value$<0.01$ for canines and incisors.

Note that the use of a low-rank tensor model simplifies interpretation substantially over (e.g.) a model that allows the covariate coefficients to vary freely over different components of the PSB process and different tooth types. An online interactive visualization showing the posterior predictive distribution for proportion of diseased sites for any given tooth type and combination of individual-level predictors is available at http://ericfrazerlock.com/toothdata/ToothDisplay.html.

6 Discussion

For our motivating application to the GAAD data, our low-rank PSB model provides sufficient flexibility without sacrificing parsimony, and inherently accommodates the multiway outcome. While we have focused on our motivating application with a binomial likelihood, the methodological ideas can be extended to other contexts. For example, the model easily extends to other likelihoods; the hierarchical prior for the weights $\pi$ (Figure 2) and steps 3-12 of the Gibbs sampling algorithm in the Appendix A remain the same. Moreover, while in our context the outcomes are a $2$-way array (patients $\times$ tooth types), an analogous approach may be used when the outcomes or predictors are of higher order. Even if the data do not take the form of a multiway array, one can still use a PSB formulation where the coefficients (components $\times$ covariates) have a low-rank matrix factorization. Such an approach may be particularly useful for variable selection problems, in which some covariates have no effect on any of the component weights (Chung and Dunson, 2009). For this, one could use a sparsity-inducing prior, such as a spike-and-slab (Ishwaran et al., 2005) or horseshoe (Carvalho et al., 2010) prior, directly on the covariate loadings in the low-rank factorization (e.g., $B_{1}$ in (6)). All the above are possible avenues for future research, and will be considered elsewhere.

Acknowledgments

The authors thank the Center for Oral Health Research at the Medical University of South Carolina for providing the motivating data. Lock and Bandyopadhyay’s research were supported by grants ULI RR033183/KL2 RR0333182 and R01DE024984, respectively, from the US National Institutes of Health.

Appendix A Gibbs sampling algorithm

The Gibbs sampling algorithm for the full model given in Figure 2, using the parameter augmentation described in Section 3, proceeds as follows:

Update $C_{ij}$’s, where

[TABLE] 2. 2.

Update the atoms $\theta_{h}$ for $h=1,\ldots,H$. Let $\tilde{Y}_{h}$ and $\tilde{n}_{h}$ be the total number of affected sites and total number of sites, respectively, that belong to observations allocated to DP component $j$:

[TABLE]

The full conditional distribution of $\theta_{h}$ is then

[TABLE] 3. 3.

Update the $Z^{*}_{ijh}$’s:

[TABLE]

where $N(\cdot,1)\mathbf{1}_{\mathbb{R}^{-}}$ and $N(\cdot,1)\mathbf{1}_{\mathbb{R}^{+}}$ define a normal distribution truncated below or above [math], respectively. 4. 4.

Update the coefficient tensor $\mathbb{B}=[\![B_{1},B_{2},B_{3}]\!]$. Define the array $\mathbb{Z}^{*}_{B}$ where $\mathbb{Z}^{*}_{B}[i,j,h]=Z^{*}_{ijh}-Z_{jh}-\mathbb{E}_{ijh}$. Then, $\mathbb{Z}^{*}_{B}=\langle X,\mathbb{B}\rangle_{1}$ where $\mathbb{Z}^{*}_{B}:I\times J\times H$, $X:I\times D$, and $\langle\cdot,\cdot\rangle_{1}$ defines the contracted tensor product over dimension $D$. We sample from the full conditional distributions of $B_{1}$, $B_{2}$, and $B_{3}$ as in Section 8 of Lock (2018). 5. 5.

Update the individual effect tensor $\mathbb{E}=[\![E_{1},E_{2},E_{3}]\!]$. Define the array $\mathbb{Z}^{*}_{E}$ where $\mathbb{Z}^{*}_{E}[i,j,h]=Z^{*}_{ijh}-Z_{jh}-X[i,:]\mathbb{B}[:,j,h]$. The entries of $\mathbb{Z}^{*}_{E}$ are independent and normally distributed with mean array $\mathbb{E}$ and variance $1$. The full conditional distributions for $E_{1}$, $E_{2}$ and $E_{3}$ can each be found by an application of well-known results for the Bayesian linear model (Lindley and Smith, 1972). For example, consider $E_{1}$. Let the columns of $E_{23}:JH\times R_{e}$ be given by the outer product of the respective columns in $E_{2}$ and $E_{3}$:

[TABLE]

where vec is the vectorization operator. Let $Z^{*\text{mat}}_{E}:I\times JH$ give the matricization of $\mathbb{Z}^{*}_{E}$, i.e.,

[TABLE]

Then, $E_{1}$ regresses $Z^{*\text{mat}}_{E}$ on $E_{23}$, giving

[TABLE]

where $\Sigma_{E_{1}}$ is the diagonal covariance matrix with $\Sigma_{E_{1}}[r,r]=\sigma^{2}_{r}$. The full conditional distributions for $E_{2}$ and $E_{3}$ are analogous to that for $E_{1}$, but with diagonal prior covariance $\Sigma_{E_{2}}[r,r]=\Sigma_{E_{3}}[r,r]=1$. 6. 6.

Update the $\sigma^{2}_{r}$’s: $\sigma_{r}^{2}\sim\mbox{Inverse-Gamma}\left(0.1+\frac{1}{2}I,0.1+{frac}{1}{2}\sum_{i=1}^{I}E_{1}[i,r]^{2}\right)$ 7. 7.

Update the $Z_{jh}$’s. Define the array $\mathbb{Z}^{*}_{Z}$ where $\mathbb{Z}^{*}_{Z}[i,j,h]=Z^{*}_{ijh}-X[i,:]\mathbb{B}[:,j,h]-\mathbb{E}_{ijh}$.

[TABLE] 8. 8.

Update $\alpha$: $\alpha\sim\mbox{Normal}\left(\frac{\sum_{j=1}^{J}\sum_{h=1}^{H}Z_{jh}}{JH+1},\frac{1}{JH+1}\right)$ for $r=1,\ldots,R_{e}$. 9. 9.

Update the stick-breaks $V_{ijh}$ for $h=1,\ldots,H$. For $h=1,\ldots,H-1$ is

[TABLE]

with $V_{ijH}=1$. 10. 10.

Update the weights $\pi_{ijh}=V_{ijh}\prod_{l<h}(1-V_{ijl})$ for $h=1,\ldots,H$. 11. 11.

Update $p_{ij}=\theta_{C_{ij}}$ for each pair $(i,j)$. 12. 12.

Determine component weights $\pi_{ijh}$ via

[TABLE]