Functional clustering methods for binary longitudinal data with temporal heterogeneity

In this paper, the DP is used as a prior distribution to cluster $\boldsymbol{\alpha}_{i}(t)$, and this clustering procedure is equivalent to clustering the set of $(\boldsymbol{\phi}_{i},\boldsymbol{\gamma}_{i})$, where $\boldsymbol{\phi}_{i}=(\boldsymbol{\phi}_{i1},\dots,\boldsymbol{\phi}_{ip})^{\top}$ and $\boldsymbol{\gamma}_{i}=(\boldsymbol{\gamma}_{i1},\dots,\boldsymbol{\gamma}_{ip})^{\top}$, as implied in (4). This process assigning the DP prior to the set of $(\boldsymbol{\phi}_{i},\boldsymbol{\gamma}_{i})$ is expressed as

where $\nu>0$ and $\mathcal{H}_{0}$ is a base distribution that randomly generates cluster-level parameters for $(\boldsymbol{\phi}_{i},\boldsymbol{\gamma}_{i})$. As another representation of the DP, it is worthwhile to look at the stick-breaking process (Sethuraman, 1994; Ishwaran and James, 2001) that allows the truncation of the summation in the DP after a large $K$ component, i.e.,

where $\boldsymbol{\phi}_{k}^{*}=(\boldsymbol{\phi}_{k1}^{*},\dots,\boldsymbol{\phi}_{kp}^{*})^{\top}$ and $\boldsymbol{\gamma}_{k}^{*}=(\boldsymbol{\gamma}_{k1}^{*},\dots,\boldsymbol{\gamma}_{kp}^{*})^{\top}$ are the parameters for cluster $k$, $\delta_{(\boldsymbol{\phi}_{k}^{*},\boldsymbol{\gamma}_{k}^{*})}(\cdot)$ is a Dirac measure at $(\boldsymbol{\phi}_{k}^{*},\boldsymbol{\gamma}_{k}^{*})$, $\pi_{k}$ is the probability mass at atom $(\boldsymbol{\phi}_{k}^{*},\boldsymbol{\gamma}_{k}^{*})$, and $K$ is a finite truncation for the maximum number of clusters. The equation in (7) implies that the model in (6) explores latent subpopulations by limiting the maximum number of subpopulations to $K$, not to infinity; see Ishwaran and James (2001) for theoretical arguments. Meanwhile, the set of cluster-level parameters, $\boldsymbol{\phi}_{k}^{*}$ and $\boldsymbol{\gamma}_{k}^{*}$, is drawn from the base distribution $\mathcal{H}_{0}$, and the random weight $\pi_{k}$ derives from a set of random variables that each follows a beta distribution, i.e.,

where $\mathbf{V}=\{V_{1},\dots,V_{K}\}$ and $V_{K}=1$, which guarantees the sum of all random weights is equal to one. Then, we can write $P(C_{i}=k|\mathbf{V})=\pi_{k}(\mathbf{V})$. The specification of $\nu$ may affect clustering performance. As $\nu$ goes to 0, the concentration toward the existing clusters gets stronger by decreasing the probability that a vector $(\boldsymbol{\phi}_{i},\boldsymbol{\gamma}_{i})$ forms a new cluster. In this work, we set $\nu=1$ by default so that every subject has the equal probability for shaping a new cluster.

This section discusses the specification of prior distributions of each model component. The set of indicator variables $\boldsymbol{\gamma}_{kl}^{*}$ has a beta-binomial prior distribution, i.e.,

for $k=1,\dots,K$ and $l=1,\dots,p$, where $B(\cdot,\cdot)$ denotes the beta function. If $a=b=1$, this prior distribution allocates equal probabilities for the number of active knots (Scott and Berger, 2010). This choice has been shown work successfully for function estimation with knot selection (Jeong and Park, 2016; Jeong et al., 2017).

For the basis coefficients $\boldsymbol{\phi}_{(\boldsymbol{\gamma}_{k}^{*})}^{\star}$ that have varying dimension in each iteration, we consider the following mixture of $g$-priors,

By characterizing the scale parameter with the total number of subjects, the prior in (9) corresponds to a Zellner-Siow prior, which is a multivariate Cauchy prior marginally for $\boldsymbol{\phi}_{(\boldsymbol{\gamma}_{k}^{*})}^{\star}$ (Liang et al., 2008). Therefore, the base measure $\mathcal{H}_{0}$ is constructed by combining a beta-binomial distribution and a multivariate Cauchy distribution.

The prior in (9) has several desirable properties. First, the prior distribution is invariant to linear transformations of the design matrix (Zellner, 1986). This means that the posterior distribution of the varying-coefficients is not affected by linear transformations of the basis functions in (3). More importantly, the prior in (9) utilizes the population-level covariance, which is determined by assuming that all individuals belong to the $k$th cluster. This specific structure enhances the convergence efficiency of MCMC because the cluster-level basis coefficients $\boldsymbol{\phi}_{(\boldsymbol{\gamma}_{k}^{*})}^{\star}$ of empty clusters are sampled by taking advantage of the summed information of all subjects. As a result, the prior can naturally begin a new cluster to which a few subjects may belong. Indeed, it is well known that choosing a reasonable scale for the base measure $\mathcal{H}_{0}$ is very important in using the DP prior (Gelman et al., 2015). In this regard, our prior construction in (9) has a clear advantage over the related studies, which use the DP prior but do not fully account the concern of scale (Ray and Mallick, 2006; Rodriguez and Dunson, 2014). Furthermore, having the right scale may be difficult with other penalty priors. For example, the Bayesian P-spline is a widely used Bayesian approach for nonparametric regression (Lang and Brezger, 2004). Since its covariance structure plays an important role in smoothing, however, the prior distribution cannot be simply modified to have a reasonable scale for the DP prior.

The remaining specification of priors for the fixed-dimensional parameter is standard. As for the fixed effects $\boldsymbol{\beta}$, we assign a multivariate normal distribution whose covariance matrix is positive definite, i.e.,

For practical purposes, $\mathbf{P}$ can be chosen as a diagonal matrix with large diagonal entries. For random effects, a multivariate normal distribution is used to generate the effects,

where $\boldsymbol{\Psi}$ is the covariance matrix of the random effects and has an inverse-Wishart prior,

In our study, $u$ and $\mathbf{D}$ are fixed in advance to make the prior distribution diffuse.

Given the prior distributions in Section 3.2, we propose a sampling algorithm used to simulate the target posterior distribution,

where $\boldsymbol{\gamma}^{*}=\{\boldsymbol{\gamma}^{*}_{1},\dots,\boldsymbol{\gamma}^{*}_{K}\}$, $\boldsymbol{\phi}_{(\boldsymbol{\gamma}^{*})}^{\star}=\{\boldsymbol{\phi}_{(\boldsymbol{\gamma}_{1}^{*})}^{\star},\dots,\boldsymbol{\phi}_{(\boldsymbol{\gamma}_{K}^{*})}^{\star}\}$, $\mathbf{b}=\{\mathbf{b}_{1},\dots,\mathbf{b}_{N}\}$, $\boldsymbol{\tau}=\{\tau_{1},\dots,\tau_{K}\}$, $\mathbf{C}=\{C_{1},\dots,C_{N}\}$, $\mathbf{L}=\{\mathbf{L}_{1},\dots,\mathbf{L}_{N}\}$, and $\mathbf{Y}=\{\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\}$ denoting $\mathbf{Y}_{i}$ as a set of binary responses for subject $i$. To simulate the target posterior distribution in (10), a standard Gibbs sampler based on its full conditional distributions cannot be implemented because the dimension $|\boldsymbol{\gamma}^{*}_{k}|\times 1$ of $\boldsymbol{\phi}_{(\boldsymbol{\gamma}^{*})}^{\star}$ depends on another model component $\boldsymbol{\gamma}^{*}$. In such a varying-dimensional problem, PCG sampling avoids the need of jumping between spaces of different dimensions through marginalization, permutation, and trimming, thereby making it implementable with the expectation of faster convergence; refer to Section 4 in Park and van Dyk (2009). In this study, we consider marginalizing the random effects $\mathbf{b}$ and the basis coefficients $\boldsymbol{\phi}_{(\boldsymbol{\gamma}^{*})}^{\star}$ in (10), thereby producing the following marginal distributions,

One iteration of the PCG sampler is shown in Algorithm 1. Steps 1 and 2 are marginalized by using (11), while Steps 3 and 4 are marginalized by using (12). To maintain the target stationary distribution of Algorithm 1, the sampling steps are permuted in a specific order. Trimming is used to remove the redundant samples of components. For more applications of the PCG sampling including other varying dimensional cases, refer to Park and van Dyk (2009); Jeong and Park (2016); Jeong et al. (2017); Park and Jeong (2018); Park et al. (2019). Because the target stationary distribution of Algorithm 1 is maintained in a specific order, the change of the order of sampling steps may not guarantee the stationarity of a Markov chain, and care must be taken not to change the sampling order; refer to van Dyk and Park (2008). The details of Algorithm 1 are given in A.

Throughout the simulation, we mainly consider three groups of varying coefficients, and the total number $N$ of subjects and the number of subjects in each cluster will be specified later based on simulation setups. The values of an underlying effect modifier $t$, $\{t_{ij}:i=1,\dots,N,j=1,\dots,n_{i}\}$, are randomly generated from a uniform distribution between 0 and 1, and the values of known covariates for subject $i$, i.e., $\mathbf{W}_{i}$, $\mathbf{X}_{i}$, and $\mathbf{Z}_{i}$, are independently generated from a standard normal distribution, except that the first column of both $\mathbf{W}_{i}$ and $\mathbf{Z}_{i}$ is set to a column vector of 1’s. Within the range of $t$ between 0 and 1, $\alpha_{kl}(t)$ denotes a varying-coefficient function of the $l$th covariate in cluster $k$. The varying coefficients for three clusters are constant, linear, or nonlinear, as described below:

The true values of the other model parameters are set to $\boldsymbol{\beta}=(\beta_{1},\beta_{2})^{\top}=(1,-1)^{\top}$ and

The latent response $L_{ij}$ of the probit varying-coefficient mixed model is drawn from

where $\boldsymbol{\alpha}_{C_{i}}(t_{ij})=(\alpha_{C_{i}1}(t_{ij}),\alpha_{C_{i}2}(t_{ij}))^{\top}$ and it is used to generate a series of binary data such that $Y_{ij}=I(L_{ij}>0)$ for observation $j$ on subject $i$.

In this section, we demonstrate the performance of the proposed method under various simulation setups. We consider three different scenarios: Scenario I, where each cluster has 400 subjects ($N=1200$), Scenario II, where each cluster has 200 subjects ($N=600$), and Scenario III, where the three clusters have 600, 400, and 200 subjects, respectively ($N=1200$). For each scenario, three different values of the concentration parameter are also considered to examine the robustness of the DP prior: $\nu\in\{0.1,1,10\}$. The proposed method is applied to each combination of the scenarios for the sample size and the concentration parameter with 300 replications of the datasets. We run 20,000 iterations of the PCG sampler, discarding the first half of the draws as burn-in and using the second half for our posterior analysis.

The side-by-side boxplots in Figure 1 illustrate clustering performance for the 300 replicated datasets in terms of precision, recall, and F1-score. The clustering labels of all subjects are chosen by the posterior modes. As expected, clusters with larger sample sizes perform better in clustering. The metrics show similar results across different values of the concentration parameter $\nu$, indicating that clustering performance is robust to the specification of the hyperparameter for the DP prior. We calculate coverage probabilities of the pointwise 95% credible bands with the 300 replicates, where the bands are specified by the 2.5% and 97.5% posterior quantiles. Figure 2 shows that the coverage probabilities are close to the nominal value of $0.95$, which validates the uncertainty quantification through the posterior distribution. The coverage probabilities are consistent for different values of $\nu$, further supporting the robustness of our proposed method against the hyperparameter specification. With 300 replicated datasets, Table 1 shows the root-mean-square errors (RMSE) of the posterior median and coverage probabilities of the 95% credible intervals for the fixed-dimensional parameters. The results indicate that the fixed-dimensional parameters are also insensitive to the hyperparameter specification.

To assess the efficiency of our proposed PCG sampler, we calculate the multivariate effective sample size (ESS) (Vats et al., 2019) for all simulation settings, as shown in Figure 3. The second half of the chain with 20,000 iterations is used to calculate the multivariate ESS, along with the running time of the sampler on a server equipped with CentOS7 and two Sky Lake CPUs @ 2.60GHz. The target parameters have a dimension of 185, including two varying-coefficient functions with 30 knots for each of three clusters and five fixed-dimensional parameters. As shown in Figure 3, the multivariate ESS is approximately 4,000–5,000 out of 10,000 iterations, implying that the proposed sampling algorithm exhibits reasonable convergence characteristics. Figure 3 also shows the multivariate ESS divided by time (seconds), providing the number of independent draws obtained per unit time. The results demonstrate that roughly 1.5 to 2.5 independent draws are obtained per second. Furthermore, the concentration hyperparameter has no significant influence on the multivariate ESS, indicating the robustness of our proposed method to the specification of the hyperparameter.

Figure 4 overlays the pointwise posterior medians of the varying coefficients of the 300 replicated datasets using the default concentration parameter value of $\nu=1$. Our results show that the estimated posterior medians become closer to the true functions as sample sizes increase. We also examined the posterior medians obtained with alternative values of $\nu$, specifically $\nu=0.1$ and $\nu=10$, but found that the results were similar and therefore omitted here.

To evaluate the necessity of subpopulation modeling, we compared our proposed method with two competing approaches that do not account for functional clustering: Jeong et al. (2017) and the mgcv package (Wood, 2017). Jeong et al. (2017) is a fully Bayesian method for estimating probit varying-coefficient mixed models with a homogeneous population using free-knot splines (Smith and Kohn, 1996). The mgcv package uses penalized quasi-likelihood to estimate the same model. To avoid redundancy, we compared the two competing methods with our proposed method using replicated datasets under Scenario I specified in Section 4.2

In Figure 5, we present the estimates of varying coefficients obtained by two competing methods that do not account for functional clustering. The estimated trends appear to be the average of the varying coefficients across the clusters in Scenario I, resulting in significant estimation bias for the cluster-specific effects. Furthermore, ignoring subpopulations leads to significant estimation bias for the fixed-dimensional parameters, as demonstrated in Table 2, despite the fact that they are common across all clusters. This suggests that accounting for subpopulation effects is crucial even for parameters that are shared among clusters.

We also investigate the performance of the proposed method in a homogeneous population with a single true cluster. To generate the simulation datasets, we set all subjects to have the same values of $\alpha_{11}(t)$ and $\alpha_{12}(t)$, as described in Section 4.1. Specifically, we set $\boldsymbol{\alpha}_{C_{i}}(t_{ij})=(\alpha_{11}(t_{ij}),\alpha_{12}(t_{ij}))^{\top}$ for all $i$, resulting in a single cluster. The number of subjects is set to $N=300$, and all other simulation settings are identical to those specified in Section 4.1.

Figure 6 displays the pointwise posterior medians of the varying coefficients using the proposed method and two competitors under the homogeneous population assumption. The proposed method shows a reasonably small estimation bias compared to Jeong et al. (2017) despite some deviation for the incorrectly clustered subjects, implying its applicability without knowing the population structure. In contrast, mgcv yields larger estimation bias than the other two methods. Figure 7 presents the coverage probabilities of the pointwise 95% credible or confidence bands of the varying coefficients, further demonstrating the worse performance of mgcv. Table 3 summarizes the results of the fixed-dimensional parameters, indicating that the RMSEs of the proposed method are slightly larger than those for Jeong et al. (2017). Considering the flexibility of the proposed method in accounting for potentially heterogeneous populations, however, it can be deemed more useful than Jeong et al. (2017).

In this section, we consider the German Socioeconomic Panel (GSOEP) data (Riphahn et al., 2003). The dataset consists of repeated observations from 7,293 subjects in Germany for the years 1984–1988, 1991, and 1994. The response variable of interest is working status (employed=1; otherwise=0) and covariates consist of $A_{ij}$ (age), $M_{ij}$ (marital status; married=1, otherwise=0), $K_{ij}$ (children under the age of 16 in the household; yes=1, no=0), $H_{ij}$ (degree of handicap; 0 to 100 in percent), and $S_{ij}$ (personal health satisfaction; 0 to 10). We confine our samples to 893 subjects under the age of 53 with Abitur degrees in order to examine the varying effect of having young children in the household on working status as a function of age for people with secondary education.

According to our preliminary analysis, which assumes varying effects for all covariates, we have decided to treat handicap, personal health satisfaction, and marital status as constant effects in subsequent analyses. As a result, we aim to model the varying effects of the intercept and the presence of children under the age of 16, while treating the remaining covariates as fixed effects. Our target model is then given by

The reduced model complexity based on the preliminary analysis leads to faster computation and greater stability.

Similar to Section 4.3, we compare the proposed method with the model in Jeong et al. (2017) that ignores heterogeneity among subjects; mgcv is not considered because Jeong et al. (2017) outperforms it (see Section 4.4). The corresponding simple model is given by

As shown in the next section, the simple model fails to account for heterogeneity among samples and can be obtained by pooling the results of the proposed target model.

We ran the proposed PCG sampler with three over-dispersed initial values. Figure 8 shows the convergence characteristics of the sampler by using the $R^{1/2}$ diagnostic for the fixed-dimensional parameters and the fixed points of the varying coefficients (Gelman and Rubin, 1992). Because all $R^{1/2}$ statistics are below 1.1, we combine the second halves of three chains each with 150,000 iterations through a label switching algorithm. Then, after thinning every 50th sample, our posterior inference is based on 4,500 mixed samples.

According to our posterior analysis, there are two main groups and a few minor groups. The group membership is determined by the posterior modes of the cluster labels. The interpretation of the analysis focuses on the two main groups. The first largest cluster, Group 1, accounts for 89.4% of subjects, the second largest, Group 2, accounts for 7.2%, and the remaining clusters account for 3.4%. Some characteristics of the two major groups are summarized in Table 4, demonstrating that these groups are made up of heterogeneous subjects. Specifically, Group 1 has a much lower proportion of females who tend to be more responsible for parenting than Group 2. In addition, Group 1 has a higher proportion of white-collar workers, civil servants, and university graduates with high job security than Group 2. Such difference in characteristics results in the different posterior estimates of varying-coefficient functions, as shown in Figure 9.

Figure 9 shows the posterior summaries of the varying-coefficient functions resulting from functional clustering. The first row of Figure 9 corresponds to the group-level varying-intercept functions. The intercept function of Group 1 is significantly positive and keeps increasing up to early 50s, implying that a posterior probability of being employed becomes higher as one tends to be old while holding all covariates constant. In contrast, Group 2 has a slightly positive but constant intercept function in all ages. The second row of Figure 9 shows the varying-coefficient function for having children below age 16 in the household. The 95% pointwise posterior intervals for Group 1 includes 0, which implies that having children below age 16 in the household does not significantly affect the probability of employment. Unlike Group 1, the existence of young children in Group 2’s household significantly decreases the probability of employment until his/her mid 40s. The probability is further decreased when the employee’s age tends to be younger. This is due in part to the fact that Group 2 has a higher proportion of females than Group 1, and females were more responsible for child care in the late twentieth century.

Table 5 shows the posterior summaries of fixed-dimensional parameters, where ${\rm Var}(b_{i})=\psi$ represents the variance of a random effect, and $\beta_{H}$, $\beta_{S}$, and $\beta_{M}$ represent the coefficients of fixed effects, $H_{ij}$, $S_{ij}$, and $M_{ij}$, respectively. Based on the fact that the 95% posterior intervals include zero, we decide that the fixed effects have no significant influence on the employment status when the varying effects of having young children in the household for heterogeneous subpopulations are accounted for in the model. When the heterogeneous subpopulation assumption is ignored, the results obtained by Jeong et al. (2017) appear to be similar, but with narrower 95% credible intervals. This is due to the fact that Jeong et al. (2017) does not fully account for the variability of heterogeneous subpopulations.

Our proposed method identifies two major subpopulations with different characteristics, as shown in Table 4, and these subpopulations show different age-varying effects of having young children in a household on working status, as illustrated in Figure 9. When the heterogeneous subpopulation assumption is ignored, however, the single population model proposed by Jeong et al. (2017) estimates the varying-coefficient functions applied to the entire population, as shown in the first column of Figure 10. In the presence of heterogeneous subpopulations, such an approach would fail to separate subpopulations with different characteristics, leading to erroneous conclusions. This is confirmed by producing the pooled varying-coefficient functions estimated by the proposed method, as shown in the second column of Figure 10. These findings demonstrate the proposed model’s validity and utility in accounting for a heterogeneous population.