Modeling for Dynamic Ordinal Regression Relationships: An Application to Estimating Maturity of Rockfish in California

We first describe our approach to regression in the context of a single distribution, that is, without any aspect of time. Let $\{(y_{i},\boldsymbol{x}_{i}):i=1,\dots,n\}$ denote the data, where each observation consists of an ordinal response $y_{i}$ along with a vector of covariates $\boldsymbol{x}_{i}=(x_{i1},\dots,x_{ip})$. The methodology is developed in DeYoreo and Kottas (2014a) for multivariate ordinal responses, however we work with a univariate response for notational simplicity and because this is the relevant setting for the application of interest. Our model assumes the ordinal responses arise as discretized versions of latent continuous responses, which is natural for many settings and particularly relevant for the fish maturity application, as maturation is a continuous variable recorded on a discrete scale. With $C$ categories, introduce latent continuous responses $(Z_{1},\dots,Z_{n})$ such that $Y_{i}=j$ if and only if $Z_{i}\in(\gamma_{j-1},\gamma_{j}]$, for $j=1,\dots,C$, and cut-offs $-\infty=\gamma_{0}<\gamma_{1}<\dots<\gamma_{C}=\infty$.

We focus on settings in which the covariates may be treated as random, which is appropriate, indeed necessary, for many environmental applications. In the fish maturity example, the body characteristics are interrelated and arise in the form of a data vector, and we are interested in various relationships, including but not limited to the way in which maturity varies with age and length. This motivates our focus on building a flexible model for the joint density $f(z,\boldsymbol{x})$, for which we apply a DP mixture of multivariate normals: $(z_{i},\boldsymbol{x}_{i})\mid G\stackrel{{\scriptstyle iid}}{{\sim}}\int\mathrm{N}(\cdot;\boldsymbol{\mu},\boldsymbol{\Sigma})dG(\boldsymbol{\mu},\boldsymbol{\Sigma})$, with $G\mid\alpha,G_{0}\sim\mathrm{DP}(\alpha,G_{0})$.

By the constructive definition of the DP (Sethuraman, 1994), a realization $G$ from a $\mathrm{DP}(\alpha,G_{0})$ is almost surely of the form $G=\sum_{l=1}^{\infty}p_{l}\delta_{\boldsymbol{\theta}_{l}}$. The locations $\boldsymbol{\theta}_{l}=$ $(\boldsymbol{\mu}_{l},\boldsymbol{\Sigma}_{l})$ are independent realizations from the centering distribution $G_{0}$, and the weights are determined through stick-breaking from beta distributed random variables. In particular, let $v_{l}\stackrel{{\scriptstyle iid}}{{\sim}}\text{beta}(1,\alpha)$, $l=1,2,\dots$, independently of $\{\boldsymbol{\theta}_{l}\}$, and define $p_{1}=v_{1}$, and for $l=2,3,\dots$, $p_{l}=v_{l}\prod_{r=1}^{l-1}(1-v_{r})$. Therefore, the model for $f(z,\boldsymbol{x})$ has an almost sure representation as a countable mixture of multivariate normals, and implies the following induced model for the regression function:

$$\mathrm{Pr}(Y=j\mid\boldsymbol{x};G)=\sum_{r=1}^{\infty}w_{r}(\boldsymbol{x})\int_{\gamma_{j-1}}^{\gamma_{j}}\text{N}(z;m_{r}(\boldsymbol{x}),s_{r})\,dz$$

(1)

with covariate-dependent weights $w_{r}(\boldsymbol{x})\propto p_{r}\text{N}(\boldsymbol{x};\boldsymbol{\mu}_{r}^{x},\boldsymbol{\Sigma}_{r}^{xx})$, covariate-dependent means $m_{r}(\boldsymbol{x})=\mu_{r}^{z}+\boldsymbol{\Sigma}_{r}^{zx}(\boldsymbol{\Sigma}_{r}^{xx})^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{r}^{x})$, and variances $s_{r}=\Sigma_{r}^{zz}-\boldsymbol{\Sigma}_{r}^{zx}(\boldsymbol{\Sigma}_{r}^{xx})^{-1}\boldsymbol{\Sigma}_{r}^{xz}$. Here, $\boldsymbol{\mu}_{r}$ is partitioned into ${\mu}_{r}^{z}$ and $\boldsymbol{\mu}_{r}^{x}$ according to $Z$ and $\boldsymbol{X}$, and $({\Sigma}_{r}^{zz},\boldsymbol{\Sigma}_{r}^{xx},\boldsymbol{\Sigma}_{r}^{zx},\boldsymbol{\Sigma}_{r}^{xz})$ are the components of the corresponding partition of covariance matrix $\boldsymbol{\Sigma}_{r}$.

This modeling strategy allows for non-linear, non-standard relationships to be captured, and overcomes many limitations of standard parametric models. In addition, the cut-offs may be fixed to arbitrary increasing values (which we recommend to be equally spaced and centered at zero) without sacrificing the ability of the model to approximate any distribution for mixed ordinal-continuous data. In particular, it can be shown that the induced prior model on the space of mixed ordinal-continuous distributions assigns positive probability to all Kullback-Leibler neighborhoods of any distribution in this space. This represents a key computational advantage over parametric models. We refer to DeYoreo and Kottas (2014a) for more details on model properties, a review of existing approaches to ordinal regression, and illustration of the benefits afforded by the nonparametric joint model over standard methods. This discussion refers to ordinal responses with three or more categories. For the case of binary regression, i.e., when $C=2$, additional restrictions are needed on the covariance matrix $\boldsymbol{\Sigma}$ to facilitate identifiability (DeYoreo and Kottas, 2014b).

In developing a model for a collection of distributions indexed in discrete time, we seek to build on previous knowledge, retaining the powerful and well-studied DP mixture model marginally at each time $t\in\mathcal{T}$, with $\mathcal{T}=\{1,2,\dots\}$. We thus seek to extend the DP prior to model $G_{\mathcal{T}}=\{G_{t}:t\in\mathcal{T}\}$, a set of dependent distributions such that each $G_{t}$ follows a DP marginally. The dynamic DP extension can be developed by introducing temporal dependence in the weights and/or atoms of the constructive definition, $G=$ $\sum_{l=1}^{\infty}p_{l}\delta_{\boldsymbol{\theta}_{l}}$.

The general formulation of the DDP introduced by MacEachern (1999, 2000) expresses the atoms $\boldsymbol{\theta}_{l\mathcal{S}}=\{\boldsymbol{\theta}_{l,t}:t\in S\}$, $l=1,2,\dots$ as independently and identically distributed (i.i.d.) sample paths from a stochastic process over $S$, and the latent beta random variables which drive the weights, $\boldsymbol{v}_{l\mathcal{S}}=\{v_{l,t}:t\in S\}$, $l=1,2,\dots$, as i.i.d. realizations from a stochastic process with beta$(1,\alpha_{t})$ marginal distributions. The distributions could be indexed in time, space, or by covariates, and $S$ represents the corresponding index set, being $\mathbb{Z}^{+}$ in our case. The DDP model for distributions indexed in discrete time expresses $G_{t}$ as $\sum_{l=1}^{\infty}p_{l,t}\delta_{{\boldsymbol{\theta}}_{l,t}}$, for $t\in\mathcal{T}$. The locations $\boldsymbol{\theta}_{l\mathcal{T}}=\{\boldsymbol{\theta}_{l,t}:t\in\mathcal{T}\}$ are i.i.d. for $l=1,2,\dots$, from a time series model for the kernel parameters. The stick-breaking weights $\boldsymbol{p}_{l\mathcal{T}}=\{p_{l,t}:t\in\mathcal{T}\}$, $l=1,2,\dots$, arise through a latent time series with beta$(1,\alpha_{t})$ marginal distributions, independently of $\boldsymbol{\theta}_{l\mathcal{T}}$.

The general DDP can be simplified by introducing dependence only in the weights, such that the atoms are not time dependent, or alternatively, the atoms can be time dependent while the weights remain independent of time. We refer to these as common atoms and common weights models, respectively. While the majority of DDP applications fall into the common weights category, we believe the most natural of the two simplifications for time series is to assume that the locations are constant over time, and introduce dependence in the weights. Although we will end up working with the more general version of the DDP with dependent atoms and weights, for the application and related settings, it seems plausible that there is a fixed set of factors that determine the region in which the joint density of body characteristics is supported, but dynamics are caused by changes in the relative importance of the factors. The construction of dependent weights requires dependent beta random variables, so that $p_{1,t}=v_{1,t}$, $p_{l,t}=v_{l,t}\prod_{r=1}^{l-1}(1-v_{r,t})$, for $l=2,3,\dots$, with each $\{v_{l,t}:t\in\mathcal{T}\}$ a realization from a time series model with beta($1,\alpha$) marginals. Equivalently, we can write $p_{1,t}=1-\beta_{1,t}$, $p_{l,t}=(1-\beta_{l,t})\prod_{r=1}^{l-1}\beta_{r,t}$, for $l=2,3,\dots$, with each $\{\beta_{l,t}:t\in\mathcal{T}\}$ a realization from a time series model with beta($\alpha,1$) marginals.

There have been many variations of the DDP model proposed in the literature. The common weights version was originally discussed by MacEachern (2000), in which a Gaussian process was used to generate dependent locations, with the autocorrelation function controlling the degree to which distributions which are “close” are similar, and how quickly this similarity decays. De Iorio et al. (2004) consider also a common weights model, in which the index of dependence is a covariate, a key application of DDP models. In the order-based DDP of Griffin and Steel (2006), covariates are used to sort the weights. Covariate dependence is incorporated in the weights in the kernel and probit stick-breaking models of Dunson and Park (2008) and Rodriguez and Dunson (2011), respectively, however these prior models do not retain the DP marginally. Gelfand et al. (2005) developed a DP mixture model for spatial data, using a spatial Gaussian process to induce dependence in DDP locations indexed by space. For data indexed in discrete time, Rodriguez and ter Horst (2008) apply a common weights model, with atoms arising from a dynamic linear model. Di Lucca et al. (2013) develop a model for a single time series of continuous or binary responses through a DDP in which the atoms are dependent on lagged terms. Xiao et al. (2015) construct a dynamic model for Poisson process intensities built from a DDP mixture with common weights and different types of autoregressive processes for the atoms. Taddy (2010) assumes the alternative simplification of the DDP with common atoms, and models the stick-breaking proportions $\{v_{l,t}:t\in\mathcal{T}\}$ using the positive correlated autoregressive beta process from McKenzie (1985). Nieto-Barajas et al. (2012) also use the common atoms simplification of the DDP, modeling a time series of random distributions by linking the beta random variables through latent binomially distributed random variables.

To generate a correlated series $\{\beta_{l,t}:t\in\mathcal{T}\}$ such that each $\beta_{l,t}\sim\mathrm{beta}(\alpha,1)$ marginally, we define a stochastic process

$$\mathcal{B}=\left\{\beta_{t}=\exp\left(-\frac{\zeta^{2}+\eta^{2}_{t}}{2\alpha}\right):\,t\in\mathcal{T}\right\},$$

(2)

which is built from a standard normal random variable $\zeta$ and an independent stochastic process ${\eta}_{\mathcal{T}}=\{\eta_{t}:t\in\mathcal{T}\}$ with standard normal marginal distributions. This transformation leads to marginal distributions $\beta_{t}\sim\text{beta}(\alpha,1)$ for any $t$. To see this, take two independent standard normal random variables $Y_{1}$ and $Y_{2}$, such that $W=$ $(Y_{1}^{2}+Y_{2}^{2})/2$ follows an exponential distribution with mean 1, and thus $B=\exp(-W/\alpha)$ $\sim\mathrm{beta}(\alpha,1)$. To our knowledge, this is a novel construction for a common atoms DDP prior model. The practical utility of the transformation in (2) is that it facilitates building the temporal dependence through Gaussian time-series models, while maintaining the DP structure marginally.

Because we work with distributions indexed in discrete time, we assume ${\eta}_{\mathcal{T}}$ to be a first-order autoregressive (AR) process, however alternatives such as higher order processes or Gaussian processes for spatially indexed data are possible. The requirement of standard normal marginal distributions on ${\eta}_{\mathcal{T}}$ leads to a restriction on the variance of the AR(1) model, such that $\eta_{l,t}\sim\text{N}(\phi\eta_{l,t-1},1-\phi^{2})$, $t=2,\dots,T$. Thus $|\phi|<1$, which implies stationarity for the stochastic process ${\eta}_{\mathcal{T}}$. Since $\mathcal{B}$ is a transformation of a strongly stationary stochastic process, it is also strongly stationary. Note that the correlation in $(\beta_{l,t},\beta_{l,t+k})$ is driven by the autocorrelation present in ${\eta}_{\mathcal{T}}$, and this induces dependence in the weights $(p_{l,t},p_{l,t+k})$, which leads to dependent distributions $(G_{t},G_{t+k})$.

We now explore this dependence, discussing some of the correlations implied by this prior model. First, consider the correlation of the beta random variables used to define the dynamic stick-breaking weights. Let $\rho_{k}=\text{corr}(\eta_{t},\eta_{t+k})$, which is equal to $\phi^{k}$ under the assumption of an AR(1) process for $\eta_{\mathcal{T}}$. The autocorrelation function associated with $\mathcal{B}$ is

$$\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)=\frac{\alpha^{1/2}(1-\rho^{2}_{k})^{1/2}(\alpha+1)^{2}(\alpha+2)^{1/2}}{\left\{(1-\rho^{2}_{k}+\alpha)^{2}-\alpha^{2}\rho^{2}_{k}\right\}^{1/2}}-\alpha(\alpha+2)$$

(3)

as described in Appendix A. Smaller values for $\alpha$ lead to smaller correlations for any fixed $\phi$ at a particular lag, and $\phi$ controls the strength of correlation, with large $\phi$ producing large correlations which decay slowly. The parameters $\phi$ and $\alpha$ in combination can lead to a wide range of correlations, however $\alpha\geq 1$ implies a lower bound near $0.5$ on the correlation for any lag $k$. In the limit, as $\alpha\rightarrow 0^{+}$, $\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)\rightarrow 0$, and as $\alpha\rightarrow\infty$, $\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)$ tends towards $0.5$ as $\rho_{k}\rightarrow 0^{+}$, and $1$ as $\rho_{k}\rightarrow 1^{-}$. Assuming $\rho_{k}=\phi^{k}$, gives $\lim_{\phi\to 1^{-}}\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)=1$ and $\lim_{\phi\to 0^{+}}\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)=\alpha^{1/2}(\alpha+1)(\alpha+2)^{1/2}-\alpha(\alpha+2).$ This tends upwards to $0.5$ quickly as $\alpha\rightarrow\infty$.

Note that $\text{corr}(\beta_{t},\beta_{t+k}\mid\alpha,\phi)$ is a function of $\rho^{2}_{k}$ but not $\rho_{k}$, which is to be expected since $\eta_{t}$ enters the expression for $\beta_{t}$ only through $\eta_{t}^{2}$, and thus $-\rho_{k}$ and $\rho_{k}$ have the same effect in the correlation. The same is true of the correlation in the DP weights, which is given below. We believe that $\phi\in(0,1)$ is a natural restriction, since we are building a stochastic process for distributions correlated in time through a transformation of an AR process, which intuition suggests should be positively correlated. However, all that is strictly required to preserve the DP marginals is $|\phi|<1$.

Assuming $\boldsymbol{\beta}_{l,\mathcal{T}}=\{\beta_{l,t}:t\in\mathcal{T}\}$ is generated by $\mathcal{B}$, from an underlying AR(1) process for ${\eta}_{\mathcal{T}}$ with coefficient $\phi$, we study the dependence induced in the resulting DDP weights at consecutive time points, $(p_{l,t},p_{l,t+1})$. The covariance is given by

$$\text{cov}(p_{l,t},p_{l,t+1}\mid\alpha,\phi)=\left\{\frac{\alpha^{3/2}(1-\phi^{2})^{1/2}}{(2+\alpha)^{1/2}\left\{(1-\phi^{2}+\alpha)^{2}-\alpha^{2}\phi^{2}\right\}^{1/2}}\right\}^{l-1}\\ \left\{1-\frac{2\alpha}{\alpha+1}+\frac{\alpha^{3/2}(1-\phi^{2})^{1/2}}{(2+\alpha)^{1/2}\left\{(1-\phi^{2}+\alpha)^{2}-\alpha^{2}\phi^{2}\right\}^{1/2}}\right\}-\frac{\alpha^{2l-2}}{(1+\alpha)^{2l}},$$

(4)

which can be divided by $\mathrm{var}(p_{l,t}\mid\alpha)=$ $\{2\alpha^{l-1}/((1+\alpha)(2+\alpha)^{l})\}-\{\alpha^{2l-2}/(1+\alpha)^{2l}\}$ to yield $\text{corr}(p_{l,t},p_{l,t+1}\mid\alpha,\phi)$; note that $\text{E}(p_{l,t}\mid\alpha)=$ $\text{E}(p_{l,t+1}\mid\alpha)$ and $\text{var}(p_{l,t}\mid\alpha)=$ $\text{var}(p_{l,t+1}\mid\alpha)$. Derivations are given in Appendix A. This expression is decreasing in index $l$, and larger values of $\phi$ lead to larger correlations in the weights at any particular $l$. Moreover, the decay in correlations with weight index is faster for small $\alpha$ and small $\phi$. As $\alpha\rightarrow 0^{+}$, $\text{corr}(p_{1,t},p_{1,t+1}\mid\alpha,\phi)\rightarrow 1$ for any value of $\phi$, and as $\alpha\rightarrow\infty$, $\text{corr}(p_{1,t},p_{1,t+1}\mid\alpha,\phi)$ is contained in $(0.5,1)$, with values closer to $1$ for larger $\phi$. Note that $\text{corr}(p_{l,t},p_{l,t+k}\mid\alpha,\phi)$ has the same expression as $\text{corr}(p_{l,t},p_{l,t+1}\mid\alpha,\phi)$, but with $\phi$ replaced by $\phi^{k}$; it is thus decreasing with the lag $k$, with the speed of decay controlled by $\phi$.

Finally, assume $G_{t}$ is a random distribution on $\mathbb{R}^{M}$, such that $G_{t}=$ $\sum_{l=1}^{\infty}p_{l,t}\delta_{\boldsymbol{\theta}_{l}}$, where the $\{p_{l,t}:t\in\mathcal{T}\}$ are defined through the dynamic stick-breaking weights $\{\beta_{l,t}:t\in\mathcal{T}\}$, and the $\boldsymbol{\theta}_{l}$ are i.i.d. from a distribution $G_{0}$ on $\mathbb{R}^{M}$. Consider two consecutive distributions $(G_{t},G_{t+1})$, and a measurable subset $A\subset\mathbb{R}^{M}$. The correlation of consecutive distributions, $\text{corr}(G_{t}(A),G_{t+1}(A)\mid\phi,\alpha,G_{0})$, is discussed in Appendix A.

We now utilize the method for incorporating dependence into the weights of the DP and the approach to ordinal regression involving DP mixtures of normals for the latent response-covariate distribution to further develop the DDP mixture modeling framework for dynamic ordinal regression.

In the original rockfish data source, maturity is recorded on an ordinal scale from $1$ to $6$, representing immature (1), early and late vitellogenesis $(2,3)$, eyed larvae $(4)$, and post-spawning $(5,6)$. Because scientists are not necessarily interested in differentiating between every one of these maturity levels, and to make the model output simple and more interpretable, we collapse maturity into three ordinal levels, representing immature $(1)$, pre-spawning mature $(2,3,4)$, and post-spawning mature $(5,6)$.

Many observations have age missing or maturity recorded as unknown. Exploratory analysis suggests there to be no systematic pattern in missingness. Further discussion with fisheries research scientists having expertise in aging of rockfish and data collection revealed that the reason for missing age in a sample is that otoliths (ear stones used in aging) were not collected or have not yet been aged. Maturity may be recorded as unknown because it can be difficult to distinguish between stages, and samplers are told to record unknown unless they are reasonably sure of the stage. Therefore, there is no systematic reason that age or maturity is not present, and it is reasonable to assume that the data are missing at random, or that the probability an observation is missing does not depend on the missing values, allowing us to ignore the missing data mechanism, and base inferences only on the complete data (e.g., Rubin, 1976; Gelman et al., 2004).

Age can not be treated as a continuous covariate, as there are approximately $25$ distinct values of age in over $2,200$ observations. Age is in fact an ordinal random variable, such that a recorded age $j$ implies the fish was between $j$ and $j+1$ years of age. This relationship between discrete recorded age and continuous age is obtained by the following reasoning. Chilipepper rockfish are winter spawning, and the young are assumed to be born in early January. The annuli (rings) of the otiliths are counted in order to determine age, and these also form sometime around January. Thus, for each ring, there has been one year of growth.

We therefore treat age much in the same way as maturity, using a latent continuous age variable. Let $U$ represent observed ordinal age, let $U^{*}$ represent underlying continuous age, and assume, for $j=1,2,\dots,$ that $U=j$ iff $U^{*}\in(j,j+1]$. Equivalently, $U=j$ iff $\log(U^{*})\in(\log(j),\log(j+1)]$, for $j=1,2,\dots$, and $U=0$ iff $\log(U^{*})\in(-\infty,0]$, so that the support of the latent continuous random variable corresponding to age is $\mathbb{R}$. Letting $W$ be the latent continuous random variable which determines $U$ through discretization, we assume $u_{t,i}=j$ iff $w_{t,i}\in(\log(j),\log(j+1)]$, for $j=0,1,\dots$, so that $W$ is interpretable as log-age on a continuous scale.

Considering year of sampling as the index of dependence, observations occur in years $1993$ through $2007$, indexed by $t=1,\dots,T=15$, with no observations in $2003$, $2005$, or $2006$. Let the missing years be given by $\boldsymbol{s}$, here $\boldsymbol{s}=\{11,13,14\}$, and let $\boldsymbol{s}^{c}=\{1,\dots,T\}\setminus\boldsymbol{s}$ represent all other years in $\{1,\dots,T\}$. This situation involving time points in which data is completely missing is not uncommon in these types of problems, and can be handled with our model for equally spaced time points. We retain the ability to provide inference and estimation in years for which no data was recorded, which we will see are reasonable and exhibit more uncertainty than in other years.

The common atoms DDP model presented in Section 2.3 was tested extensively on simulated data, and performed very well in capturing trends in the underlying distributions when there were no missing years or forecasting was not the focus. In these settings, the common atoms model is sufficiently powerful from an inferential perspective such that the need to turn to a more complex model is diminished. However, as a consequence of the need to force the same set of components to be present at each time point, density estimates at missing time points tend to resemble an average across all time points, which is not desirable when a trend or change in support is suggested. A more general model is required for this application, since it is important to make inferences in years for which data was not collected. We therefore consider a simple extension, adding dependence through a vector autoregressive model in the mean component of the DP atoms, such that $\boldsymbol{\theta}_{l}=(\boldsymbol{\mu}_{l},\boldsymbol{\Sigma}_{l})$ becomes $\boldsymbol{\theta}_{l,t}=(\boldsymbol{\mu}_{l,t},\boldsymbol{\Sigma}_{l})$.

Assuming $Z$ represents maturity $Y$ on a continuous scale, $W$ is interpretable as log-age, and $X$ represents length, a dependent nonparametric mixture model is applied to estimate the time-dependent distributions of the trivariate continuous random vectors $\boldsymbol{Y}^{*}_{ti}=(Z_{ti},W_{ti},X_{ti})$, for $t=1,\dots,T$, and $i=1,\dots,n_{t}$. In our notation, $T$ is the number of years or the final year containing data (and it is possible that not all years in $\{1,\dots,T\}$ contain observations), and $n_{t}$ is the sample size in year $t$.

We utilize the computationally efficient approach to inference which involves truncating the countable representation for each $G_{t}$ to a finite level $N$ (Ishwaran and James, 2001), such that the dependent stick-breaking weights are given by $p_{1,t}=1-\beta_{1,t}$, $p_{l,t}=(1-\beta_{l,t})\prod_{r=1}^{l-1}\beta_{r,t}$, for $r=1,\dots,N-1$, and $p_{N,t}=\prod_{l=1}^{N-1}\beta_{l,t}$, ensuring $\sum_{l=1}^{N}p_{l,t}=1$. Since $\alpha$ is not a function of $t$, the same truncation level is applied for all mixing distributions. In choosing the truncation level, we use the expression relating $N$ to the expectation of the sum of the first $N$ weights $w_{1},\dots,w_{N}$ generated from stick-breaking of beta$(1,\alpha)$ random variables. The expression is $\mathrm{E}(\sum_{j=1}^{N}w_{j}\mid\alpha)=$ $1-(\alpha/(\alpha+1))^{N}$, which can be further averaged over the prior for $\alpha$ to obtain $\mathrm{E}(\sum_{j=1}^{N}w_{j})$, with $N$ chosen such that this expectation is close to 1 up to the desired level of tolerance for the approximation. The hierarchical model can be expressed as follows:

	$\displaystyle y_{t,i}=j\leftrightarrow\gamma_{j-1}<z_{t,i}\leq\gamma_{j},\,t\in\boldsymbol{s}^{c},\,i=1,\dots,n_{t}$
	$\displaystyle u_{t,i}=j\leftrightarrow\log(j)<w_{t,i}\leq\log(j+1),\,t\in\boldsymbol{s}^{c},\,i=1,\dots,n_{t}$
	$\displaystyle\{\boldsymbol{y}^{*}_{t,i}\}\mid\{\boldsymbol{\mu}_{l,t}\},\{\boldsymbol{\Sigma}_{l}\},\{L_{t,i}\}\sim\prod_{t\in\boldsymbol{s}^{c}}\prod_{i=1}^{n_{t}}\text{N}(\boldsymbol{\mu}_{L_{t,i},t},\boldsymbol{\Sigma}_{L_{t,i}})$
	$\displaystyle\{L_{t,i}\}\mid\{\eta_{l,t}\},\{\zeta_{l}\}\sim\prod_{t\in\boldsymbol{s}^{c}}\prod_{i=1}^{n_{t}}\sum_{l=1}^{N}p_{l,t}\delta_{l}(L_{t,i})$
	$\displaystyle\zeta_{l}\stackrel{{\scriptstyle iid}}{{\sim}}\text{N}(0,1),\quad l=1,\dots,N-1$
	$\displaystyle\eta_{l,1}\stackrel{{\scriptstyle iid}}{{\sim}}\text{N}(0,1),\quad l=1,\dots,N-1$
	$\displaystyle\eta_{l,t}\mid\eta_{l,t-1},\phi\sim\text{N}(\phi\eta_{l,t-1},1-\phi^{2}),\quad l=1,\dots,N-1,\,t=2,\dots,T$
	$\displaystyle\boldsymbol{\mu}_{l,1}\mid\boldsymbol{m_{0}},\boldsymbol{V_{0}}\sim\text{N}(\boldsymbol{m_{0}},\boldsymbol{V_{0}}),\quad l=1,\dots,N$
	$\displaystyle\boldsymbol{\mu}_{l,t}\mid\boldsymbol{\mu}_{l,t-1},\Theta,\boldsymbol{m},\boldsymbol{V}\sim\text{N}(\boldsymbol{m}+\Theta\boldsymbol{\mu}_{l,t-1},\boldsymbol{V}),\quad l=1,\dots,N,\,t=2,\dots,T$
	$\displaystyle\boldsymbol{\Sigma}_{l}\mid\nu,\boldsymbol{D}\stackrel{{\scriptstyle iid}}{{\sim}}\text{IW}(\boldsymbol{\Sigma}_{l};\nu,\boldsymbol{D}),\quad l=1,\dots,N$		(5)

with priors on $\alpha$, $\boldsymbol{\psi}=(\boldsymbol{m},\boldsymbol{V},\boldsymbol{D})$, $\phi$, and $\Theta$. Recall that $\beta_{l,t}$ is defined through $(\eta_{l,t},\zeta_{l},\alpha)$, and the $\{\beta_{l,t}\}$ determine the $\{p_{l,t}\}$ through stick-breaking.

The parameters $\{\zeta_{l}\}$ and $\{\eta_{l,t}\}$ can be updated individually with slice samplers, which involves alternating simulation from uniform random variables and truncated normal random variables, implying draws for $\{\beta_{l,t}\}$, and hence $\{p_{l,t}\}$. The configuration variables $L_{t,i}$ are drawn from discrete distributions on $\{1,\dots,N\}$, with probabilities proportional to $p_{l,t}\text{N}(\boldsymbol{y}^{*}_{t,i};\boldsymbol{\mu}_{l,t},\boldsymbol{\Sigma}_{l})$ for $l=1,\dots,N$. The update for $\Sigma_{l}$ is $\text{IW}(\nu+M_{l},\boldsymbol{D}+\sum_{\{(t,i):L_{ti}=l\}}(\boldsymbol{y}^{*}_{t,i}-\boldsymbol{\mu}_{l,t})(\boldsymbol{y}^{*}_{t,i}-\boldsymbol{\mu}_{l,t})^{T})$, where $M_{l}=|\{t,i\}:L_{t,i}=l|$, and each $\boldsymbol{\mu}_{l,t}$ is updated from a normal distribution. The parameters $\alpha$ and $\phi$, given priors $\text{IG}(a_{\alpha},b_{\alpha})$, and uniform on $(0,1)$ or $(-1,1)$, respectively, can be sampled using Metropolis-Hastings steps. We assume that $\Theta$ is diagonal, with elements $(\theta_{1},\theta_{2},\theta_{3})$, however we advocate for a full covariance matrix $\boldsymbol{V}$. This implies that each element of $\boldsymbol{\mu}_{l,t}$ has a mean which depends only on the corresponding element of $\boldsymbol{\mu}_{l,t-1}$, however there exists dependence in the elements of $\boldsymbol{\mu}_{l,t}$, which seems reasonable for most applications. We assume uniform priors on $(0,1)$ or $(-1,1)$ for each element of $\Theta$, and update them with a Metropolis-Hastings step. Finally, the parameters $\boldsymbol{\psi}$ have closed-form full conditional distributions, given priors $\boldsymbol{m}\sim\text{N}(\boldsymbol{a_{m}},\boldsymbol{B_{m}})$, $\boldsymbol{V}\sim\text{IW}(a_{\boldsymbol{V}},\boldsymbol{B_{V}})$, $\boldsymbol{D}\sim\text{W}(a_{\boldsymbol{D}},\boldsymbol{B_{D}})$. The full conditionals and posterior simulation details are further described in Appendix B.

To implement the model, we must specify the parameters of the hyperpriors on $\boldsymbol{\psi}$. A default specification strategy is developed by considering the limiting case of the model as $\alpha\rightarrow 0^{+}$ and $\Theta\rightarrow\boldsymbol{0}$, which results in a single normal distribution for $\boldsymbol{Y}^{*}_{t}$. In the limit, with $\boldsymbol{Y}^{*}_{t}\mid\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}\sim\text{N}(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma})$ and $\boldsymbol{\mu}_{t}\mid\boldsymbol{m},\boldsymbol{V}\sim\text{N}(\boldsymbol{m},\boldsymbol{V})$, we find $\text{E}(\boldsymbol{Y}^{*}_{t})=\boldsymbol{a_{m}}$ and $\text{Cov}(\boldsymbol{Y}^{*}_{t})=\boldsymbol{B_{m}}+\boldsymbol{B_{V}}(a_{\boldsymbol{V}}-d-1)^{-1}+a_{\boldsymbol{D}}\boldsymbol{B_{D}}(\nu-d-1)^{-1}$, where $d$ is the response-covariate dimension, here $d=3$. The only covariate information we require is an approximate center (such as the midpoint of the data) and range, denoted by $c^{x}$ and $r^{x}$ for $X$, and analogously for $U$. We use $c^{x}$ and $r^{x}/4$ as proxies for the marginal mean and standard deviation of $X$. We also seek to scale the latent variables appropriately. The centers and ranges $c^{u}$ and $r^{u}$ provide approximate centers $c^{w}$ and ranges $r^{w}$ of latent log-age $W$. Since $Y$ is supported on $\{1,\dots,C\}$, latent continuous $Z$ must be supported on values slightly below $\gamma_{1}$ up to slightly above $\gamma_{C-1}$, so that $r^{z}/4$ is a proxy for the standard deviation of $Z$, where $r^{z}=(\gamma_{C-1}-\gamma_{1})$. Using these mean and variance proxies, we fix $a_{\boldsymbol{m}}=(0,c^{u},c^{x})$. Each of the three terms in $\text{Cov}(\boldsymbol{Y}^{*}_{t})$ can be assigned an equal part of the total covariance, for instance being set to $3^{-1}\text{diag}\{(r^{z}/4)^{2},(r^{w}/4)^{2},(r^{x}/4)^{2}\}$. For dispersed but proper priors, $\nu$, $a_{\boldsymbol{V}}$ and $a_{\boldsymbol{D}}$ can be fixed to small values such as $d+2$, and $\boldsymbol{B_{m}}$, $\boldsymbol{B_{V}}$, and $\boldsymbol{B_{D}}$ determined accordingly.

It remains to specify $\boldsymbol{m}_{0}$ and $\boldsymbol{V}_{0}$, the mean and covariance for the initial distributions $\boldsymbol{\mu}_{l,1}$. We propose a fairly conservative specification, noting that in the limit, $\text{E}(\boldsymbol{Y}^{*}_{1})=\boldsymbol{m}_{0}$, and $\text{Cov}(\boldsymbol{Y}^{*}_{1})=a_{\boldsymbol{D}}\boldsymbol{B_{D}}(\nu-d-1)^{-1}+\boldsymbol{V}_{0}$. Therefore, $\boldsymbol{m}_{0}$ can be specified in the same way as $\boldsymbol{a_{m}}$ but using only the subset of data at $t=1$, and $\boldsymbol{V}_{0}$ can be set to $\text{diag}\{(r^{z}_{1}/4)^{2},(r^{w}_{1}/4)^{2},(r_{1}^{x}/4)^{2}\}-a_{\boldsymbol{D}}\boldsymbol{B_{D}}(\nu-d-1)^{-1}$, where the subscript $1$ indicates the subset of data at $t=1$.

In simulation studies and the rockfish application, we observed a moderate to large amount of learning for all hyperparameters. For instance, for the rockfish data, the posterior distribution for $\phi$ was concentrated on values close to 1, indicating the DDP weights are strongly correlated across time. There was also moderate learning for $\alpha$ as its posterior distribution was concentrated around $0.5$, with small variance, shifted down relative to the prior which had expectation $2$. The posterior distribution for each element of $\boldsymbol{m}$ was reduced in variance and concentrated on values not far from those indicated by the prior mean. The posterior samples for the covariance matrices $\boldsymbol{V}$ and $\boldsymbol{D}$ supported smaller variance components than suggested by the prior.

Various simulation settings were developed to study both the common atoms version of the model and the more general version (DeYoreo, 2014, chapter 4). While we focus only on the fish maturity data application in this paper, our extensive simulation studies have revealed the inferential power of the model under different scenarios for the true latent response distribution and ordinal regression relationships.

We first discuss inference results for quantities not involving maturity. As illustrated with results for six years in Figure 3, the estimates for the density of length display a range of shapes. The interval estimates reflect the different sample sizes in these years; for instance, in $1994$ there are $271$ observations, whereas in $2000$ and $2004$ there are only $64$ and $37$ observations, respectively. A feature of our modeling approach is that inference for the density of age can be obtained over a continuous scale. The corresponding estimates are shown in Figure 4 for the same years as length.

The posterior mean surfaces for the bivariate density of age and length are shown in Figure 5 for all years, including the ones (years $2003$, $2005$ and $2006$) for which data is not available. The model yields more smooth shapes for the density estimates in these years. An ellipse with a slight “banana” shape appears at each year, though some nonstandard features and differences across years are present. In particular, the density in year $2002$ extends down farther to smaller ages and lengths; this year is unique in that it contains a very large proportion of the young fish which are present in the data. One can envision a curve going through the center of these densities, representing $\text{E}(X\mid U^{*}=u^{*};G_{t})$, for which we show posterior mean and $95\%$ interval bands for three years in Figure 6. The estimates from our model are compared with the von Bertalanffy growth curves for length-at-age, which are based on a particular function of age and three parameters (estimated here using nonlinear least squares). It is noteworthy that the nonparametric mixture model for the joint distribution of length and age yields estimated growth curves which are overall similar to the von Bertalanffy parametric model, with some local differences especially in year 2002. The uncertainty quantification in the growth curves afforded by the nonparametric model is important, since the attainment of unique growth curves by group (i.e. by location or cohort) is often used to suggest that the groups differ in some way, and this type of analysis should clearly take into account the uncertainty in the estimated curves.

The last year $2007$, in addition to containing few observations, is peculiar. There are no fish that are younger than age $6$ in this year, and most of the age $6$ and $7$ fish are recorded as immature, even though in all years combined, less than $10\%$ of age $6$ as well as age $7$ fish are immature. This year appears to be an anomaly. As there are no observations in $2005$ or $2006$, and a small number of observations in $2007$ which seem to contradict the other years of data, hereinafter, we report inferences only up to $2004$.

Inference for the maturation probability curves is shown over length and age in Figures 7 and 8. The probability that a fish is immature is generally decreasing over length, reaching a value near $0$ at around $350$ mm in most years. There is a large change in this probability over length in $2002$ and $2004$ as compared to other years, as these years suggest a probability close to $1$ for very small fish near $200$ to $250$ mm. Turning to age, the probability of immaturity is also decreasing with age, also showing differences in $2002$ and $2004$ in comparison to other years. There is no clear indication of a general trend in the probabilities associated with levels $2$ or $3$. Years $1995$-$1997$ and $1999$ display similar behavior, with a peak in probability of post-spawning mature for moderate length values near $350$ mm, and ages $6$-$7$, favoring pre-spawning mature fish at other lengths and ages. The last four years $2001$-$2004$ suggest the probability of pre-spawning mature to be increasing with length up to a point and then leveling off, while post-spawning is favored most for large fish. Post-spawning appears to have a lower probability than pre-spawning mature for any age at all years, with the exception of $1998$, for which the probability associated with post-spawning is very high for older fish.

The Pacific States Marine Fisheries Commission states that all Chilipepper rockfish are mature at around $4$-$5$ years, and at size $304$ to $330$ mm. A stock assessment produced by the Pacific Fishery Management Council (Field, 2009) fitted a logistic regression to model maturity over length, from which it appears that $90\%$ of fish are mature around $300$-$350$ mm. As our model does not enforce monotonicity on the probability of maturity across age, we obtain posterior distributions for the first age greater than $2$ (since biologically all fish under $2$ should be immature) at which the probability of maturity exceeds $90\%$, given that it exceeds $90\%$ at some point. That is, for each posterior sample we evaluate $\mathrm{Pr}(Y>1\mid u^{*};G_{t})$ over a grid in $u^{*}$ beginning at $2$, and find the smallest value of $u^{*}$ at which this probability exceeds $0.9$. Note that there were very few posterior samples for which this probability did not exceed $0.9$ for any age (only $4$ samples in $1993$ and $8$ in $2003$). The estimates for age at $90\%$ maturity are shown in the left panel of Figure 9. The model uncovers a (weak) U-shaped trend across years. Also noteworthy are the very narrow interval bands in $2002$. Recall that this year contains an abnormally large number of young fish. In this year, over half of fish age $2$ (that is, of age $2$-$3$) are immature, and over $90\%$ of age $3$ (that is, of age $3$-$4$) fish are mature, so we would expect the age at $90\%$ maturity to be above $3$ but less than $4$, which our estimate confirms. A similar analysis is performed for length (right panel of Figure 9) suggesting a trend over time which is consistent with the age analysis.

Due to the monotonicity in the maturity probability curve in standard approaches, and the fact that age and length are treated as fixed, the point at which maturity exceeds a certain probability is a reasonable quantity to obtain in order to study the age or length at which most fish are mature. However, since we are modeling age and length, we can obtain their entire distribution at a given maturity level. These are inverse inferences, in which we study, for instance, $f(x\mid Y=1;G_{t})$ as opposed to $\mathrm{Pr}(Y=1\mid x;G_{t})$. It is most informative to look at age and length for immature fish, as this makes it clear at which age or length there is essentially no probability assigned to the immature category. The posterior mean estimates for $f(u^{*},x\mid Y=1;G_{t})$ are shown in Figure 10.

Here, we discuss results from posterior predictive model checking. In particular, we generated replicate data sets from the posterior predictive distribution, and compared to the real data using specific test quantities (Gelman et al., 2004). Using the MCMC output, we simulate replicate data sets $\{(y_{ti},w_{ti},x_{ti})^{rep}\}$ of the same size as the original data. We then choose some test quantity, $T(\{y_{ti},w_{ti},x_{ti}\})$, and for each replicate data set at time $t$, determine the value of the test quantity and compare the distribution of test quantities with the value computed from the actual data. To obtain Figure 11 we computed, for each replicate sample, the proportion of age $6$ fish that were of maturity levels $1$ and $2$. Boxplots of these proportions are shown, with the actual proportions from the data indicated as blue points. The width of each box is proportional to the number of age $6$ fish in that year. Figure 12 refers to fish of at least age $7$ with length larger than $400$ mm. Finally, Figure 13 plots the sample correlation of length and age for fish of maturity level $2$. The results reported in Figures 11, 12 and 13 suggest that the model is predicting data which is very similar to the observed data in terms of practically important inferences.

Although results are not shown here, we also studied residuals with cross-validation, randomly selecting $20\%$ of the observations in each year and refitting the model, leaving out these observations. We obtained residuals $\tilde{y}_{ti}-\text{E}({Y}\mid W=\tilde{w}_{ti},X=\tilde{x}_{ti};G_{t})$ for each observation $(\tilde{y}_{ti},\tilde{w}_{ti},\tilde{x}_{ti})$ which was left out. There was no apparent trend in the residuals across covariate values, that is, no indication that we are systematically under or overestimating fish maturity of a particular length and/or age.