An application of Dirichlet process in clustering subjects via variance shift models: A course-evaluation study

Abstract

In this article, the Dirichlet process (DP) is applied to cluster subjects with longitudinal observations. The basis of clustering is the ability of subjects to adapt themselves to new circumstances. Indeed, the basis of clustering depends on the time of changing response variability. This is done by providing a random change-point time in the variance structure of mixed-effects models. The DP is assumed as a prior for the distribution of the random change point. The discrete nature of the DP is utilized to cluster subjects according to the time of adaption. The proposed model is useful to identify groups of subjects with distinctive time-based progressions or declines. Transition mixed-effects models are also used to account for the serial correlation among observations over time. A joint modelling approach is utilized to handle the bias created in these models. The Gibbs sampling technique is adopted to achieve parameter estimates. Performance of the proposed method is evaluated via conducting a simulation study. The usefulness of the proposed model is assessed on a course-evaluation dataset.

Keywords

course-evaluation data Dirichlet process prior Markov chain Monte Carlo Model-based clustering transition mixed-effects models variance shift models

1 Introduction

Cluster analysis is commonly employed to determine the intrinsic grouping in a set of unlabelled data. In this way, objects within a cluster are more similar to each other, in comparison to those in the other clusters. The notation of similarity can differ depending on the purpose of study. Clustering techniques have been applied for different purposes in many fields such as marketing, biology, medicine, econometrics and social sciences. Many clustering algorithms have been adopted to address different concepts of similarity measures. There are good reviews of clustering methods (see, e.g., Xu and Wunsch, 2005).

The model-based or distribution-based clustering approach, which utilizes a mixture of distributions, is one of the most adaptable clustering approaches in facing with datasets involving complex structures (see, e.g., McLachlan and Peel, 2000; McNicholas, 2016). In this approach, clusters consist of subjects belonging most likely to the same distribution. This clustering approach can be flexibly applied to cluster data with different sources of dependency. A very important and widely used type of data having a dependent structure is longitudinal data, which involves repeatedly measuring subjects recorded through time. There are two major sources of dependency in longitudinal data: intra-class correlation and serial correlation (Diggle et al., 2002).

In recent years, clustering of longitudinal data, sometimes called clustering of trajectories, has been the focus of attention. Overviews of model-based clustering approaches in longitudinal researches are provided by Vermunt (2010), Frühwirth-Schnatter (2011) and Bock (2014). Because of the rich structure of longitudinal data, different clustering purposes can be assumed to achieve homogeneous clusters in different studies. For example, Juarez and Steel (2010) assume a transition mixed-effects model and consider clustering based on constructing homogeneous groups in which the common behaviour may arise either from the dynamics (i.e., the coefficient of the lagged response), the effects of the covariates or from the equilibrium level of the time series. Suarez and Ghosal (2016) and Vogt and Linton (2016) among many others consider functional data clustering where individual observations are viewed as realizations of a random function. In this setting, finding representative curve patterns corresponding to different shapes and variations is the purpose.

A very prominent feature of the longitudinal data is the ability to monitor the evolution of underlying variables over time. Clustering subjects through time is among the most empirically important topics in econometrics, medicine, and social and behavioural sciences (see, e.g., Nagin, 1999; Weisburd et al., 2004; Juarez and Steel, 2010; Nie et al., 2010; Airila et al., 2014; Theodore et al., 2015; Carbonneau et al., 2016; Collins et al., 2016; Graziane et al., 2016). On the other hand, the behaviour change of subjects through time has always been the centre of attention (see, e.g., Nigg, 2001; Tobia and Inauen, 2010; Fleig et al., 2015). However, to the best of my knowledge, none of these researches has pointed out the problem of clustering subjects based on the time of changing subject's behaviour, as discussed in this article. Indeed, the idea of this article is motivated by the problem of clustering longitudinal data including the evaluation scores of 1 154 professors teaching at the Isfahan University of Technology (IUT) in Iran, during the years 1998–2014. It is seen that the variability of evaluation scores for each professor is large at first and then is reduced after the passage of time and gaining more experience. Since the way each professor teaches depends on many latent and unobserved factors, reductions in the variability for different individuals arise at different time points. Therefore, it can be important to cluster professors based on the time of reaching these change points. Accessing the stability of teacher's behaviour and effectiveness has been the focus of much attention over the last decade (e.g., Goldhaber and Hansen, 2013; Morgan et al., 2014).

There are many other situations in which the subject shows a very irregular and uncontrolled behaviour at the start of a process. For example, changes in monetary policy may affect firms’ stock volatility (see, e.g., Xu et al., 2016). Finding groups of firms with common change points in variances can be of great importance in practice to explore the effective features of firms that can contribute to the formation of groups. In these situations, there is usually a large variability for the response variable (the variable of interest) at the beginning of the process, whereas, after the passage of time and adaption to the new circumstance, the behaviour of subjects will be more stable. Usually, subjects who have a similar performance in adapting with the conditions show similar manners. Thus, identifying groups of subjects with a similar ability in attaining a stable time point in which their response variability is controlled may be essential. The reverse problem may also happen, as usually found in the literature of degradation. The reliability prediction based on degradation modelling is a useful method to estimate reliability for some highly reliable components or systems when there are rare failures (Wang and Coit, 2007). It is usual that the variability in any given degradation measure increases with the usage time (Wang and Coit, 2007; Rathod et al., 2011). Clustering of components of a large and complex system (or population) based on the time of starting or speeding up the degradation can be important for grading products, understanding the effective factors on degradation of parts and achieving more accurate reliability predictions. Also, the same process can occur regarding some medical or behavioural features of individuals, such as reaction time, memory and fluid intelligence measures, which have been shown to have more variability by aging, as mentioned in Morse (1993) and references therein.

Categories obtained in clustering based on the time of changing variability can then be used for further statistical inference in clusters with more homogeneous subjects. This is particularly important when doing some specific statistical inference is the key point in the data analysis process, and because of a very non-homogeneous population, the resulting inference cannot efficiently be used for most elements of the population. Moreover, by clustering subjects, it may be feasible to specify common factors effective in the creation of clusters. It may also be possible to test if subjects categorized in the same cluster show similar responses to different stimulants. Various methods have been designed to this purpose; as a good reference, see Nagin (2005). Although regression models with shifts in variance components have been investigated vastly as variance shift models (e.g., Babadi et al., 2014; Gumedze and Chatora, 2014; Li et al., 2015; Xu et al., 2016), none of these researches has yet pointed out the problem of clustering subjects based on the time of changing response variability. These current researches have been mainly designed to detect observations with large variations and to make models robust.

To formulate the above clustering idea, I assume that the variation of response measurements will change after certain time lag. This assumption corresponds to a variance shift model, in which the variance of residual term will change. In practice, changes usually happen at different time points for different subjects. Therefore, it is more realistic to assume a random change point in the structure of a mixed-effects model. In some applications, it is usually important that different subjects can be categorized based on behaviour of their responses over time. In these studies, it usually takes time for a subject to gain experience and be adapted to the conditions; thus the variability of response will usually be reduced after a change-point. Similarly, the reverse problem may happen such that the variance of response measured on parts may be increased after a change-point time due to the degradation. Being able to cluster subjects based on this similarity measure, that is, the time of changing response variability, the Dirichlet process (DP; Ferguson, 1973) is considered as a prior for the distribution of the random change point. According to the DP, the distribution $G$ , as a random measure, is distributed on the space of all possible distribution functions. The DP consists of two components: the base distribution $G_{0}$ , which represents our primary guess about the true shape of $G$ , and a scale parameter $M$ , which shows our uncertainty about $G_{0}$ . Then, the property of DPs in producing discrete distributions is applied to cluster subjects.

A prominent feature of assuming the DP as a prior for the distribution of the random change point is that the number of clusters can be inferred based on the data structure and not in an a priori manner. Applications of the DP to cluster analysis of data with special structures have appeared in many disciplines (e.g., Pennell and Dunson, 2007; Li et al., 2010; Heinzl et al., 2012; Heinzl and Tutz, 2013; Wang and Wang, 2013; Rikhtehgaran and Kazemi, 2016).

To study the association of responses over time, transition mixed-effects models are used by assuming lagged responses in the mean structure of the model (Crouchley and Davies, 2001; Diggle et al., 2002). The bias created in these models is handled by joint modelling of the initial responses and subsequent responses (Kazemi and Crouchley, 2006).

Estimates of model parameters are obtained by the use of the Gibbs sampler algorithm (Geman and Geman, 1984), in a Bayesian perspective. The performance of the proposed model is evaluated by conducting a simulation study as well as applying the model on a real dataset.

The remainder of this article is organized as follows. Sections 2 introduces the course-evaluation dataset. Section 3 specifies the transition mixed-effects model. In Section 4, the semi-parametric approach of clustering subjects with longitudinal observations is introduced. Section 5 presents the proposed transition variance shift model. Section 6 includes a simulation study. In Section 7, the usefulness of the proposed model is shown in clustering of a course-evaluation dataset. The last section includes conclusions.

2 Motivating data-set: Course-evaluation data

The data includes 12 354 records from 1 154 official or contractual professors teaching at the IUT in Iran during the years 1998–2014. The course evaluation is recorded via an electronic questionnaire, which consists of a series of questions in order to evaluate the instruction of a given course. The dataset is unbalanced such that the number of teaching semesters varies between 1 and 32. The data includes all course-evaluation scores for 14 departments of the university, including engineering and the basic science. Among these records, $72.73 %$ belongs to the undergraduate courses and the rest is for the graduate courses. Also, the $98.83 %$ of the records relate to the pure courses versus applied ones. The evaluation scores varies in the range 9.5–20. Among 1 154 professors, 255 individuals have only one record and thus are ignored in data analysis. Also, 37 professors have the history of more than 32 terms teaching at the IUT and thus are ignored in the current study. Therefore, 862 professors are considered for further study. In Figure 1, the left panel depicts the profile of evaluation scores against the number of teaching semesters. As is seen, after the passage of time and gaining more experience about class management and teaching, the variability of evaluation scores has reduced. Also, as is seen in the right panel of Figure 1, the average of evaluation scores is gradually increased by increasing the experience. Also, a jump is revealed from the first to the second semester, which can show the adaption of professors with the circumstances after one term teaching.

It is assumed that the variability of evaluation scores for each professor is large at first and is reduced after the passage of time and gaining more experience. Since the way of teaching each professor depends on many latent and unobserved factors, they get to the stability at different time points. It can be important to cluster professors based on the time of reaching these stabilization points. For example, it is possible to examine the relationship between the grade-point-average (GPA) of classes and their corresponding evaluation scores concerning the professors assigned to the same cluster instead of assessing this relationship among the population of whole professors. Therefore, the relationship between these two factors can be evaluated more precisely. Inference about this relationship may be an important issue in judgement about course-evaluation systems. In the next sections, I develop statistical models that aim to flexibly model the data structure and cluster subjects based on the time of reaching stabilization points.

3 Specification of transition mixed-effects models

Let $Y_{it}$ denotes the $t$ -th measurement taken on the $i$ -th individual, $t = 1, \dots, T_{i}$ , $i = 1, \dots, n$ . The transition mixed-effects model (Diggle et al., 2002; Hsiao, 2002) for longitudinal data is given by

Y_{it} = x_{it}^{'} β + γ y_{i, t - 1} + α_{i} + ε_{it},

(3.1)

where the

p

-dimensional vactor

x_{it}

includes known covariates and the vector

β

includes regression coefficients. In this model, the

α_{i}

’s denote random intercepts which represent between-subject variability due to the unobserved time-invariant covariates and the

ε_{it}

’s are residual terms which represents within-subject variability due to unobserved time-varying covariates and measurement error. The covariate

y_{i, t - 1}

is the lagged-response variable for the

i

-th subject which represents state dependence. The coefficient

γ

measures the effects of previous outcomes on their current values. It is usually assumed that random effects and

ε_{it}

’s are independent for all

i

and

t

, and are independent of covariates. Also, the usual assumptions are that the random effects

α_{i} \overset{iid}{\sim} N (0, σ_{α}^{2})

and the residuals

ε_{it} \overset{iid}{\sim} N (0, σ_{ε}^{2})

In fitting transition models, the presence of random effects which are added to capture the between-subject variability may cause the so-called initial conditions problem. This problem happens because of the correlation between the random effects and the initial state $Y_{i 0}$ . Ignoring this correlation, in the case of small to moderate $T_{i}$ ’s and large $n$ , produces seriously biased parameter estimates (Kazemi and Davies, 2002). A pragmatic solution is to consider the following equation for the $Y_{i 0}$ and then jointly model it with the subsequent responses (Kazemi and Crouchley, 2006). More specifically, I consider

Y_{i 0} = x_{i 0}^{'} β_{0} + u_{i 0}, ​ ​ ​ i = 1, \dots, n,

(3.2)

where

p_{0}

-dimensional vectors

x_{i 0}

’s are initial period covariates and the

β_{0}

includes corresponding regression coefficients. I further assume that the

u_{i 0} = ϕ α_{i} + ε_{i 0}

, where

ε_{i 0} \overset{iid}{\sim} N (0, σ_{ε_{1}}^{2})

. Combining Equations (3.1) and (3.2), the full model is given by

Y_{i} = {\tilde{X}}_{i} θ + α_{i} w_{i} + ε_{i}, ​ ​ i = 1, \dots, n,

(3.3)

where

Y_{i} = {(Y_{i 0}, Y_{i 1}, \dots, Y_{{iT}_{i}})}^{'}

θ = {(β_{0}^{'}, γ, β^{'})}^{'}

w_{i} = (ϕ, 1_{T_{i}}^{'})^{'}

, the model matrix

{\tilde{X}}_{i} = {({\tilde{X}}_{1_{i}}^{'}, {\tilde{X}}_{2_{i}}^{'})}^{'}

, where

{\tilde{X}}_{1_{i}} = (x_{i 0}^{'}, 0, 0_{p}^{'})

and

{\tilde{X}}_{2_{i}} = (0_{T_{i} \times p_{0}}, y_{i, - 1}, X_{i})

, where

y_{i, - 1}

corresponds to

y_{i}

lagged by one-time period and

X_{i} = (x_{i 1}, \dots, x_{{iT}_{i}})^{'}

. Also, it is usually assumed that

ε_{i} \overset{ind}{\sim} N_{T_{i} + 1} (0, σ_{ε}^{2} I_{T_{i} + 1})

Maximum likelihood estimates of model parameters can easily be obtained (see, e.g., Kazemi and Crouchley, 2006).

4 Clustering subjects using Dirichlet processes

In many longitudinal studies, clustering of subjects with longitudinal observations is the purpose. In the model-based clustering approach, a mixture distribution with normal components for random effects is usually applied to cluster subjects. In these cases, the determination of the number of mixture components, or say subjects clusters is usually the main concern. To deal with this issue, the Bayesian semi-parametric approach of the DP is applied in the literature (e.g., Li et al., 2010; Heinzl et al., 2012; Heinzl and Tutz, 2013; Wang and Wang, 2013). Ferguson (1973) introduced the DP as a random probability measure defined on the space of all possible distribution functions. Indeed, the random probability measure $G$ follows a DP with the base distribution $G_{0}$ defined on a measurable space $(ϒ, A)$ and a scale parameter $M$ ( $> 0$ ), if for any finite measurable partition $A_{1}, \dots, A_{k}$ of $A$ , we have

(G (A_{1}), \dots, G (A_{n})) \sim Dir ({MG}_{0} (A_{1}), \dots, {MG}_{0} (A_{n})),

(4.1)

where

Dir

denotes the Dirichlet distribution. Let denote the distribution of

G

DP (M, G_{0})

. The base distribution

G_{0}

is the expectation of the random process and states our belief about the true shape of

G

. The scale parameter

M

, which appears in the variance of the random measure

G

, expresses the amount of uncertainty about

G_{0}

. Distributions drawn from a DP are discrete with probability one. This feature of the DP is the basis of using DP in the model-based clustering approach. Sethuraman (1994) developed an alternative representation for the DP called stick-breaking which utilizes the discreteness property of the DP to introduce the probability mass function for this random probability measure. In this representation,

\begin{matrix} G = \sum_{j = 1}^{\infty} π_{j} δ_{ξ_{j}}, \\ ξ_{j} \overset{iid}{\sim} G_{0}, \\ π_{j} = γ_{j} \prod_{i = 1}^{j - 1} (1 - γ_{i}), \\ γ_{j} \overset{iid}{\sim} Beta (1, M), \end{matrix}

(4.2)

where

δ_{ξ_{j}}

denotes a degenerated distribution of unit mass centred at

ξ_{j}

. I denote the process of generating

π_{j}

’s by

π_{j} \sim SB (M)

. In practice, the above infinite summation has to be truncated up to a finite integer

C

As was mentioned, the base distribution $G_{0}$ is the expected value of the process. Thus, the discrete distributions drawn from the DP are around the base distribution. However, the variation of $G$ around $G_{0}$ is controlled by the dispersion parameter $M$ in a stochastic way such that, if $M$ gets large then $G$ becomes close to $G_{0}$ .

The discrete nature of the DP makes it useless in those data modelling where a continuous distribution is required. However, this restriction is relaxed by the introduction of Dirichlet process mixture (DPM) models (Escobar, 1994; MacEachern, 1994) by adding a hierarchy level to the model. Indeed, it is assumed that

\begin{matrix} Y_{i} \overset{ind}{\sim} f (y_{i} | ξ_{i}), ​ i = 1, \dots, n, \\ ξ_{i} | G \overset{iid}{\sim} G, ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ \\ G \sim DP (M, G_{0}) . ​ ​ ​ ​ ​ ​ ​ ​ \end{matrix}

(4.3)

Thus, the DPM model has an interpretation as an infinite mixture model by assuming the stick-breaking representation of the DP.

Nevertheless, the discreteness property of draws from the DP can be used for clustering purposes. Thus, in the DPM model, several $ξ_{i}$ ’s are prone to take similar values, such that the number of distinct values of $ξ_{i}$ , denoted by $K$ , is less than or equal to $n$ . These distinct values, which I denote by $ξ_{j}^{*}$ , $j = 1, \dots, K$ , present clusters. The number of clusters, $K$ , can indirectly be controlled by the dispersion parameter $M$ .

4.1 The identifiability issue

When Markov chain Monte Carlo (MCMC) simulation methods such as the Gibbs sampler and the Metropolis–Hastings (MH) algorithm (Hastings, 1970) are applied to achieve Bayes estimates of parameters for mixture models, the so-called label switch problem takes place (Frühwirth-Schnatter, 2006). Indeed, because of the invariance of the likelihood function with respect to the permutation of the component labels in the mixture model, the marginal posterior distribution of the parameters for all components is identical. Therefore, the Bayesian inference is not appropriate. To handle the label switch problem, the constraint $ξ_{1}^{*} < \dots < ξ_{K}^{*}$ is usually placed on the model, where $ξ_{j}^{*} \overset{iid}{\sim} G_{0}$ , $j = 1, \dots, K$ .

5 The proposed transition variance shift model

For clustering subjects based on the time of changing response variability, I first consider the transition model given in Equation (3.1) and then allowing for subject-specific change points $c_{i}$ , $i = 1, \dots, n$ , in the variance structure of the residuals. Indeed, I let $ε_{it} | c_{i} \overset{ind}{\sim} N (0, σ_{ε_{it}}^{2})$ , where the scale parameter $σ_{ε_{it}}^{2}$ is set to be $σ_{ε_{1}}^{2}$ for $t \leq c_{i}$ , and equals $σ_{ε_{2}}^{2}$ for $t > c_{i}$ . To attain a flexible modelling structure together with being able to cluster subjects, the DP is considered as a prior for unknown distribution of $c_{i}$ ’s. More specifically, it is assumed that $c_{i} | G \overset{iid}{\sim} G$ where $G \sim DP (M, G_{0})$ and $G_{0}$ is a discrete uniform distribution on ${1, \dots, T^{*}}$ , where $T^{*}$ is the maximum of $T_{i}$ ’s, $i = 1, \dots, n$ . I assume $α_{i} | c_{i} \overset{ind}{\sim} N (κ c_{i}, σ_{α}^{2})$ . If $κ$ is estimated near zero, it means that there is no significant correlation between these two random effects.

To handle the initial conditions problem, I consider Equation (3.2) and assume that the $Y_{i 0}$ is dependent on the $c_{i}$ through $α_{i}$ . Then, full model is given by Equation (3.3) and letting $ε_{i} \overset{ind}{\sim} N_{T_{i} + 1} (0, Σ_{c_{i}})$ where $Σ_{c_{i}}$ is a positive-definite blocked-diagonal matrix of the form $Σ_{c_{i}} = bldiag (σ_{ε_{1}}^{2} I_{c_{i} + 1}, σ_{ε_{2}}^{2} I_{T_{i} - c_{i}})$ .

The above random intercept model is inspired by the longitudinal study on course-evaluation data whereby individual profiles show subject-specific shift in variances over time. Note, however, that based on the purpose of study, different clustering strategies can be considered by assuming the DP as the prior of the random effects of models with different structures.

5.1 Bayesian estimation

The parameter estimates can be obtained from the freely available software Openbugs (Lunn et al., 2009). Being applicable in this software, I use an equivalent structure introduced by Ishwaran and James (2001) which introduces labels $λ_{i}$ ’s, $i = 1, \dots, n$ , and considers the following model

\begin{matrix} Y_{i} | c_{j}^{*}, ϕ, α_{i}, θ, λ_{i}, σ_{ε_{1}}^{2}, σ_{ε_{2}}^{2} ​ \overset{ind}{\sim} N_{T_{i} + 1} ({\tilde{X}}_{i} θ + α_{i} w_{i}, Σ_{c_{λ_{i}}^{*}}), \\ λ_{i} \overset{ind}{\sim} G (λ_{i}) = \sum_{j = 1}^{C} π_{j} Δ_{j} (λ_{i}), ​ ​ where ​ ​ π_{j} \sim SB (M), \\ c_{j}^{*} \overset{iid}{\sim} DU {1, \dots, T^{*}}, ​ j = 1, \dots, C, \end{matrix}

(5.1)

where

α_{i} | c_{λ_{i}}^{*} \overset{ind}{\sim} N (κ c_{λ_{i}}^{*}, σ_{α}^{2})

and

DU

denotes a discrete uniform distribution.

Fitting complex models, such as the proposed transition variance shift model, requires the use of MCMC methods. To obtain Bayesian parameter estimates using the MCMC simulation methods, I apply conditionally conjugate priors with hyper parameters being chosen such that the corresponding priors be vague. The Gibbs sampler simulates iteratively from the complete conditional posterior distribution of each unknown parameter, given the current values of all other model parameters and observations. Then, means of MCMC samples are used to achieve Bayes estimates.

In the Appendix, I provide the chosen priors for the model parameters and all full conditionals to implement the Gibbs sampling.

6 A simulation study

The data generating process is organized to the mixed model

Y_{it} = β_{0} + β_{1} x_{it} + γ y_{i, t - 1} + α_{i} + ε_{it}, ​ ​ i = 1, \dots, 300, ​ t = 1, \dots, 30 .

(6.1)

I set

β_{0} = 0

β_{1} = 1

and

γ = 0.5

. Also, I assume that the

x_{it}

’s, which produce part of response variability, are independent and identically distributed (i.i.d.) random variables drawn from

N (0, σ_{x}^{2})

. Also,

α_{i}

’s and

ε_{it}

’s, which, respectively, represent between- and within-subject variabilities, are i.i.d. random variables each drawn from

N (0, σ_{α}^{2})

and

N (0, σ_{ε_{it}}^{2})

. Parameters

σ_{x}^{2}

and

σ_{α}^{2}

are both set to one. For each subject,

σ_{ε_{it}}^{2}

is set to

25

for

t \leq c_{i}

and is set to

1

otherwise. Indeed, the standard deviation of residuals becomes five times larger after the change points to make the change points distinguishable. For each subject,

c_{i}

is randomly generated from the five equally distanced values between 0 and 30, that is, the set

{5, 10, 15, 20, 25}

. Therefore, each subject is randomly assigned to a cluster. In this way, 61, 52, 59, 70 and 58 subjects are assigned to the first through fifth clusters, respectively. Model given in Equation (3.2) is considered to generate

y_{i 0}

’s with

β_{00} = 0

and

β_{01} = 1

. I set

ϕ = 3.98

, leading to

corr (Y_{i 0}, α_{i} | c_{i}) = 0.62

, which shows a moderate correlation between

Y_{i 0}

and

α_{i}

when conditioning on the change point

c_{i}

, that is, for all individuals in the same cluster. Figure 2 shows generated data points for 30 randomly selected subjects.

I fit the proposed model of Section 5, where $G_{0}$ is let to be the discrete uniform distribution on ${1, \dots, 30}$ . To implement the Gibbs sampler, the following independent priors are adopted: $N (0, 100)$ for each regression coefficient, $ϕ$ and $κ$ , $U (0, 1)$ for $γ$ , and $IG (0.01, 0.01)$ for variance components. To be sure that in the fitted model, $σ_{ε_{1}}^{2}$ is larger than $σ_{ε_{2}}^{2}$ , I assume that $σ_{ε_{1}}^{2} = ε \times σ_{ε_{2}}^{2}$ , where $ε \sim G (0.01, 0.01) I (1, \infty)$ . I set the scale parameter $M = 1$ and the truncation parameter $C = 10$ . Then, by means of the Openbugs software, 150 000 MCMC samples are generated after considering 1 000 burn-in. Results, after convergence is achieved, are reported in Table 1. Since $α_{i}$ and $c_{i}$ , for $i = 1, \dots, n$ , are assumed to be independent in the process of generating the set dataset, the $κ$ is estimated near zero. In each Gibbs sampler run, sizes of 10 primary assumed clusters are computed. Then, Bayes estimates of change points together with the Bayes median (BMed), the Bayes 2.5 $%$ quantile (Q1) and the Bayes 97.5 $%$ quantile (Q3) of the sizes of 10 primary assumed clusters are reported in Table 2. Based on these results, it is seen that among 10 assumed clusters, there are only 5 active (non-empty) clusters. Other clusters have been inactive at least for $97.5 %$ of the Gibbs sampler runs. Thus, the number of clusters is estimated to be five. Then, the likelihood of each subject based on the Bayes estimates of parameters obtained from the Openbugs software is computed in the Matlab software (The Mathworks Inc., 2012) for clusters 1 to 5. Then, each subject is assigned to the cluster with the largest individual likelihood. After that, the size of clusters is computed, as reported in Table 2. Next, the estimated cluster of each subject is compared with its true cluster. The ratio of true clustering is obtained to be $94 %$ . Figure 3 depicts profiles of subjects belonging to each cluster. The estimated change points are also depicted with vertical lines.

Table 1:

Bayesian estimation results for the simulated dataset

$β_{0}$	$β_{1}$	$γ$	$σ_{ε_{1}}^{2}$	$σ_{ε_{2}}^{2}$	$σ_{α}^{2}$	$β_{00}$	$β_{01}$	$ϕ$	$κ$
$-$ 0.1727	0.9948	0.4999	24.9400	0.9936	0.9289	$-$ 0.4601	0.9093	3.9107	0.0110
(0.1379)	(0.0150)	(0.006)	(0.5283)	(0.0224)	(0.0880)	(0.3304)	(0.2815)	(0.3202)	(0.0082)

Note: Bayesian standard deviations are given in parentheses.

Table 2:

Bayesian estimates of change points, the BMed and (Q1,Q3) for the size of each cluster, obtained from the Openbugs software, are reported for the simulated dataset. Also, the size of each cluster computed by the use of Matlab software are reported

	A	B	C	D	E	F	G	H	I	J
$c_{j}^{*}$	5.418	6.138	6.858	10.580	15.440	16.150	16.870	20.590	25.480	27.260
BMed	62	0	0	53	59	0	0	67	59	0
(Q1, Q3)	(60, 63)	(0, 0)	(0, 0)	(50, 56)	(56, 62)	(0, 0)	(0, 0)	(64, 71)	(57, 61)	(0, 0)
Size of cluster	66	-	-	51	64	-	-	63	56	-

Comparing the estimated change points $c_{1}^{*}, c_{4}^{*}, c_{5}^{*}, c_{8}^{*}$ and $c_{9}^{*}$ with the true values and also noting to the good estimates of model parameters reported in Table 1, together with achieving 94 $%$ true ratio for clustering, all confirm the usefulness of the proposed model.

7 The analysis of course-evaluation data

In order to classify professors based on the time of changing evaluation-scores variability, I fitted model given in Equation (5.1) by assuming no covariate in the structure of the model, that is, the pure values of evaluation scores are considered for the clustering purpose and no confounding effects of covariates are extracted from the evaluation scores. I also let $c_{i} | G \overset{iid}{\sim} DP (M, G_{0})$ , where $G_{0}$ is assumed to be a discrete uniform distribution on the set ${1, \dots, 32}$ . I set $M = 1$ for the dispersion parameter of the DP to prevent over-clustering.

Table 3:

Bayes estimates for the course-evaluation data

$β_{0}$	$β_{00}$	$ϕ$	$γ$	$σ_{ε_{1}}^{2}$	$σ_{ε_{2}}^{2}$	$σ_{α}^{2}$
16.610	-	8.009	0.051	1.638	0.409	0.753
(0.054)	-	(9.764)	(0.003)	(0.057)	(0.012)	(0.045)

Note Bayesian standard deviations are given in parentheses.

To implement the Gibbs sampler, let priors be the same as those adopted in the simulation study. To improve the rate of convergence to the target posterior distributions, I set $κ = 0$ . I use the Openbugs software by setting 120 000 burn-in and then running 80 000 samples. Results, after the convergence is achieved, are reported in Table 3.

The stick-breaking representation of the DP is approximated by $C = 10$ components. Then, the number of active clusters, $K$ , is considered as the number of clusters. To show that $C = 10$ components suffices to approximate the infinite mixture, the estimated posterior probability density function of parameter $K$ and the history of produced values of $K$ are shown in panels of Figure 4. Lack of tendency towards larger number of clusters indicates that good approximations of the DP are achieved by the use of the truncated stick-breaking representation.

Table 4:

Bayes estimates of change points and the BMed obtained from the Openbugs software are reported for the course-evaluation data. Also, the size of each cluster computed by the use of Matlab software and the ratio of censored subjects (RCS) in each cluster are reported

	A	B	C	D	E	F	G	H	I	J
$c_{j}^{*}$	1.337	4.700	7.897	11.610	13.740	16.390	19.180	22.050	24.960	30.510
BMed	3466	66	45	94	0	0	0	0	0	79
Size of cluster	150	215	85	94	-	-	-	-	-	318
RCS	0	0.7302	0.6706	0.5426	-	-	-	-	-	0.8019

Table 4 shows the Bayes estimates of the stabilization time points and the BMed of size of primary assumed clusters from the MCMC outputs. It is seen that clusters A, B, C, D and J have been active in the most Gibbs sampler runs, while other clusters have been inactive at least for the 50 $%$ of Gibbs sampler runs. Afterwards, the parameter estimates of model are used in the Matlab software to classify professors. Then, each professor is assigned to the cluster with the largest individual likelihood. The size of each cluster is reported in Table 4. As is seen, about $16.69 %$ of professors reach their stabilization points during the first two semesters. However, about $35.37 %$ of them have not reached their stability even after 30 semesters. It should be noted that because of the unbalanced data, there may be professors in each cluster who only have been recorders till the stabilization time of the corresponding cluster. For these professors, it is inferred that their stabilization time points are at least the stabilization point of their corresponding clusters. Indeed, their records are censored in their underlying clusters. The ratio of these professors in each cluster are reported as the RCS in Table 4. Figure 5 shows profiles of evaluation scores in each of the estimated clusters. In these graphs, professors with censored records are ignored. The estimated stabilization time points are also depicted with vertical lines. Now, these clusters can be used for further statistical inference.

8 Conclusion

In this article, I extended the application of the DP prior for clustering subjects with longitudinal observations. This was specifically done by incorporating the DP as the prior for the unknown distribution of the random change point in the structure of the variance components of a transition mixed-effects model. Indeed, I assumed that the measure of similarity for model-based clustering of subjects is the time of changing response variability. In some applications, mainly in medicine, industry, and social and behavioural sciences, it usually takes time for a subject to gain experience and be adapted to the conditions and thus the variability of response will usually be reduced after a change point. Similarly, in the current literature, for example, the new parts usually start with small variations, while after the passage of time and the degradation, their variability will increase. Therefore, it is usually important that different subjects can be categorized based on the behaviour of their responses over time. Identifying groups of subjects with a similar ability in adaption to the circumstances is important for post-processing issues. Some examples are exploring common factors that can contribute to the formation of clusters, testing if subjects of the same cluster show the same responses to the stimulants and assessing some relationships among variables for the more homogeneous subjects in clusters.

It was shown by a simulation study that the ratio of true clustered subjects was high in the proposed model. In a real data application, the proposed model was used to cluster professors based on the stabilization term in which, the variability of the course-evaluation scores was reduced. Clustering of professors based on this similarity concept could be important for investigating different hypotheses in each cluster with more homogeneous individuals instead of the non-homogeneous population of whole individuals. One of the most challenging issues about the course-evaluation system is to assess the relationship between the course-evaluation scores and the GPA of classes.

The proposed model can also be applied to the field of medicine, economy or the industry (degradation). As another example than the course-evaluation study, I can point out the work of Xu et al. (2016) in detecting and identifying a common variance change point among different countries.

In their work, the impact of the gross domestic product (GDP) growth rate on the variation of unemployment rate among countries was investigated. This impact was addressed by regressing unemployment rate on GDP growth rate and assuming a common change point in the variance of residuals. Nevertheless, fitting a mixed-effects model and assuming a country-specific change point in the variance structure may produce more flexibility to the model and improve the efficiency of the model parameters’ estimates. On the other hand, applying the proposed model could cluster countries based on the time when the impact of GDP growth rate on the variation of unemployment rate becomes smaller. Then, further issues such as finding common features in countries of the same cluster may be informative.

Further, the methodology can be improved by clustering subjects based on the time of changing other features of models. The model can also be extended to the variance shift models with several jumping points for each subject.

Acknowledgments

The author is grateful to the academic affairs of the Isfahan University of Technology for making the course-evaluation data available.

Appendix

The following distributions are adopted as priors for the proposed model:

Inverse-Gamma prior, $IG (τ_{1_{j}}, τ_{2_{j}})$ , for $σ_{ε_{j}}^{2}$ for $j = 1, 2$ ,

Inverse-Gamma prior, $IG (δ_{1}, δ_{2})$ , for $σ_{α}^{2}$ ,

Multivariate normal prior, $N_{p_{1}} (θ_{0}, Λ)$ , for $θ$ , where $p_{1} = p + p_{0} + 1$ ,

Normal prior, $N (μ_{κ}, σ_{κ}^{2})$ , for $κ$ ,

Normal prior, $N (ϕ_{0}, σ_{ϕ}^{2})$ , for $ϕ$ .

Full conditionals required for running the Gibbs sampler for the proposed model are derived as follows:

$θ | σ_{ε_{1}}^{2}, σ_{ε_{2}}^{2}, ϕ, α, σ^{2}, y \sim N_{p_{1}} (μ_{θ}, Σ_{θ}),$ with the mean vector

μ_{θ} = Σ_{θ} (θ_{0}^{'} Λ^{- 1} + \frac{1}{σ_{ε 1}^{2}} \sum_{i = 1}^{n} {\tilde{r}}_{i, c_{λ_{i}}^{*}}^{'} {\tilde{X}}_{i, c_{λ_{i}}^{*}} + \frac{1}{σ_{ε 2}^{2}} \sum_{i = 1}^{n} ​ {\tilde{r}}_{i, T_{i} - c_{λ_{i}}^{*}}^{'} {\tilde{X}}_{i, T_{i} - c_{λ_{i}}^{*}}), ​ ​ ​

and the covariance matrix

Σ_{θ} = {(Λ^{- 1} + \frac{1}{σ_{ε_{1}}^{2}} \sum_{i = 1}^{n} {\tilde{X}}_{i, c_{λ_{i}}^{*}}^{'} {\tilde{X}}_{i, c_{λ_{i}}^{*}} + \frac{1}{σ_{ε_{2}}^{2}} \sum_{i = 1}^{n} {\tilde{X}}_{i, T_{i} - c_{λ_{i}}^{*}}^{'} {\tilde{X}}_{i, T_{i} - c_{λ_{i}}^{*}})}^{- 1},

where each matrix with index

c_{λ_{i}}^{*}

is consisted of the first to the

(c_{λ_{i}}^{*} + 1)

-th rows and columns of the corresponding matrix. Similarly, matrices with index

T_{i} - c_{λ_{i}}^{*}

are consisted of the

(c_{λ_{i}}^{*} + 2)

-th to the

(T_{i} + 1)

-th rows and columns of the corresponding matrices. Also,

{\tilde{r}}_{i, c_{λ_{i}}^{*}} = y_{i, c_{λ_{i}}^{*}} - α_{i} w_{c_{λ_{i}}^{*}}

and

{\tilde{r}}_{i, T_{i} - c_{λ_{i}}^{*}} = y_{i, T_{i} - c_{λ_{i}}^{*}} - α_{i} w_{T_{i} - c_{λ_{i}}^{*}}

$σ_{ε_{1}}^{2} | θ, α, φ, y \sim IG (τ_{1_{1}}^{*}, τ_{2_{1}}^{*}),$ where

τ_{1_{1}}^{*} = τ_{1_{1}} + (n + \sum_{i = 1}^{n} c_{λ_{i}}^{*}) / 2,

and

τ_{2_{1}}^{*} = τ_{2_{1}} + \frac{1}{2} \sum_{i = 1}^{n} {(r_{i, c_{λ_{i}}^{*}} - α_{i} w_{c_{λ_{i}}^{*}})}^{'} (r_{i, c_{λ_{i}}^{*}} - α_{i} w_{c_{λ_{i}}^{*}}),

where

r_{i, c_{λ_{i}}^{*}} = y_{i, c_{λ_{i}}^{*}} - {\tilde{X}}_{i, c_{λ_{i}}^{*}} θ

$σ_{ε_{2}}^{2} | θ, α, φ, y \sim IG (τ_{1_{2}}^{*}, τ_{2_{2}}^{*}),$ where

τ_{1_{2}}^{*} = τ_{1_{2}} + \sum_{i = 1}^{n} (T_{i} - c_{λ_{i}}^{*}) / 2,

and

τ_{2_{2}}^{*} = τ_{2_{2}} + \frac{1}{2} {\sum_{i = 1}^{n} (r_{i, T_{i} - c_{λ_{i}}^{*}} - α_{i} w_{T_{i} - c_{λ_{i}}^{*}})}^{'} (r_{i, T_{i} - c_{λ_{i}}^{*}} - α_{i} w_{T_{i} - c_{λ_{i}}^{*}})

$φ | σ_{ε}^{2}, θ, α, y \sim N (μ_{φ}^{*}, σ_{φ}^{2 *}),$ where $μ_{φ}^{*} = σ_{φ}^{2 *} (\frac{1}{σ_{φ}^{2}} φ_{0} + \frac{1}{σ_{ε 1}^{2}} \sum_{i = 1}^{n} r_{i 0} α_{i})$ , $σ_{φ}^{2 *} = {(\frac{1}{σ_{φ}^{2}} + \frac{1}{σ_{ε 1}^{2}} {SS}_{α})}^{- 1}$ , the $r_{i 0}$ ’s are fitted residuals for $t = 0$ and ${SS}_{α} = \sum_{i = 1}^{n} α_{i}^{2}$ .

$α_{i} | σ_{ε_{1}}^{2}, σ_{ε 2}^{2}, σ_{α}^{2}, φ, y \sim N (a_{i}^{- 1} b_{i}, a_{i}^{- 1}),$ where $b_{i} = \frac{1}{σ_{ε 1}^{2}} w_{c_{λ_{i}}^{*}}^{'} r_{i, c_{λ_{i}}^{*}} + \frac{1}{σ_{ε 2}^{2}} w_{T_{i} - c_{λ_{i}}^{*}}^{'} r_{i, T_{i} - c_{i}} + \frac{1}{σ_{α}^{2}} κ c_{λ_{i}}^{*}$ and $a_{i} = \frac{1}{σ_{ε 1}^{2}} w_{c_{λ_{i}}^{*}}^{'} w_{c_{λ_{i}}^{*}} + \frac{1}{σ_{ε 2}^{2}} w_{T_{i} - c_{λ_{i}}^{*}}^{'} w_{T_{i} - c_{i}} + \frac{1}{σ_{α}^{2}}$ .

$σ_{α}^{2} | α, y \sim IG (δ_{1}^{*}, δ_{2}^{*}),$ where $δ_{1}^{*} = δ_{1} + \frac{n}{2}$ and $δ_{2}^{*} = δ_{2} + \frac{1}{2} {SS}_{α}$ .

$κ | α, c, y \sim N (μ_{κ}^{*}, σ_{κ}^{2 *}),$ where

μ_{κ}^{*} = σ_{κ}^{2 *} (\frac{1}{σ_{κ}^{2}} μ_{κ} + \frac{1}{σ_{α}^{2}} \sum_{i = 1}^{n} α_{i} c_{λ_{i}}^{*}), ​ ​ ​

and

σ_{κ}^{2 *} = {(\frac{1}{σ_{κ}^{2}} + \frac{1}{σ_{α}^{2}} \sum_{i = 1}^{n} c_{λ_{i}}^{2 *})}^{- 1} .

By the conjugacy of the multinomial and the Dirichlet distributions, I have

π_{j} = γ_{j}^{*} \prod_{r = 1}^{j - 1} (1 - γ_{r}^{*}), ​ j = 1, \dots, C - 1,

where

π_{1} = γ_{1}^{*}

and

γ_{j}^{*} | ​ λ, y \overset{iid}{\sim} Beta (1 + M_{j}, M + \sum_{r = j + 1}^{C} M_{r})

$λ_{i} | θ, σ_{ε}^{2}, φ, α_{i}, c_{λ_{i}}^{*}, π, y \sim \sum_{j = 1}^{C} π_{j} φ_{T_{i} + 1} (y_{i} | {\tilde{X}}_{i} θ + α_{i} w_{i}, Σ_{c_{λ_{i}}^{*}}) Δ_{j} (λ_{i}) .$

For each random effect, let $λ_{r}^{*}$ ’s, $r = 1, \dots, m$ be $m$ unique values of $λ_{i}$ ’s. Then,

\begin{matrix} f (c_{λ_{r}^{*}}^{*} ​ | ​ α, σ_{ε_{1}}^{2}, σ_{ε_{2}}^{2}, φ, θ, y) \propto {σ_{ε_{1}}^{2}}^{- \frac{1}{2} (n + M_{r} c_{λ_{r}^{*}}^{*})} \times {σ_{ε_{2}}^{2}}^{- \frac{1}{2} \sum_{{i : λ_{i} = λ_{r}^{*}}} (T_{i} - c_{λ_{r}^{*}}^{*})} \\ \times \exp {- \frac{1}{2} \sum_{{i : λ_{i} = λ_{r}^{*}}} \frac{1}{σ_{ε_{1}}^{2}} (r_{i, c_{λ_{r}^{*}}^{*}} - α_{i} w_{c_{λ_{i}^{*}}^{*}})^{'} (r_{i, c_{λ_{r}^{*}}^{*}} - α_{i} w_{c_{λ_{i}^{*}}^{*}}) ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ \\ + \frac{1}{σ_{ε_{2}}^{2}} (r_{i, T_{i} - c_{λ_{r}^{*}}^{*}} - α_{i} w_{T_{i} - c_{λ_{i}^{*}}^{*}})^{'} (r_{i, T_{i} - c_{λ_{r}^{*}}^{*}} - α_{i} w_{T_{i} - c_{λ_{i}^{*}}^{*}})}, ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ \end{matrix}

for

c_{λ_{r}^{*}}^{*} \in {1, \dots, T^{*}}

, where

M_{r}

is the number of

λ_{i}

’s that are equal to

λ_{r}^{*}

The MH algorithm can be used to simulate from this complete conditional posterior distribution.

References

Airila

, Hakanen

, Luukkonen

, Lusa

, Punakallio

, Leino-Arjas

(2014) Developmental trajectories of multisite musculoskeletal pain and depressive symptoms: The effects of job demands and resources and individual factors. Psychology and Health , 29, 1421–41.

Babadi

, Rasekh

, Rasekhi

, Zare

, Zadkarami

(2014) A variance shift model for detection of outliers in the linear measurement error model. Abstract and Applied Analysis , 2014, 9. DOI: http://dx.doi.org/10.1155/2014/396875.

Bock

H-H

(2014) Model-based clustering methods for time series, In Gaul W, Geyer- Schulz A, Baba Y and Okada A eds. German–Japanese Interchange of Data Analysis Results , pages 3–15. Cham: Springer.

Carbonneau

, Boivin

, Brendgen

, Nagin

, Tremblay

(2016) Comorbid development of disruptive behaviours from age 1/12 to 5 years in a population birth-cohort and association with school adjustment in first grade. Journal of Abnormal Child Psychology , 44, 677–90.

Collins

, O'Connor

, Supplee

(2016) Behaviour problems in elementary school among low-income boys: The role of teacher–child relationships. The Journal of Educational Research , DOI: 10.1080/00220671.2015.1039113

Crouchley

, Davies

(2001) A comparison of GEE and random effects models for distinguishing heterogeneity, nonstationarity and state dependence in a collection of short binary event series. Statistical Modelling: An International Journal , 1, 271–85.

Diggle

, Heagerty

, Liang

, Zeger

(2002) Analysis of Longitudinal Data , 2nd edition New York, NY: Oxford University Press.

Escobar

(1994) Estimating normal means with a Dirichlet process prior. Journal of the American Statistical Association , 89, 268–77.

Ferguson

(1973) A Bayesian analysis of some nonparametric problems. Annals of Statistics , 1, 209–30.

10.

Fleig

, Kuper

, Lippke

, Schwarzer

, Wiedemann

(2015) Cross-behaviour associations and multiple health behaviour change: A longitudinal study on physical activity and fruit and vegetable intake. Journal of Health Psychology , 20, 525–34.

11.

Frühwirth-Schnatter

(2006) Finite Mixture and Markov Switching Models . New York, NY: Springer-Verlag.

12.

Frühwirth-Schnatter

(2011) Panel data analysis: A survey on model-based clustering of time series. Advances in Data Analysis and Classification , 5, 251–80.

13.

Geman

, Geman

(1984) Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence , 6, 721–41.

14.

Goldhaber

, Hansen

(2013) Is it just a bad class? Assessing the long-term stability of estimated teacher performance. Economica , 80, 589–612.

15.

Graziane

, Beer

, Snitz

, Chang

C-C H

, Ganguli

(2016) Dual trajectories of depression and cognition: A longitudinal population-based study. The American Jou- rnal of Geriatric Psychiatry , 24, 364–73.

16.

Gumedze

, Chatora

(2014) Detection of outliers in longitudinal count data via overdispersion. Computational Statistics and Data Analysis , 79, 192–202.

17.

Hastings

(1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika , 57, 97–109.

18.

Heinzl

, Tutz

(2013) Clustering in linear mixed models with approximate Dirichlet process mixtures using EM algorithm. Statistical Modelling , 13, 41–67.

19.

Heinzl

, Fahrmeir

, Kneib

(2012) Additive mixed models with Dirichlet process mixture and P-spline priors. Adva- nces in Statistical Analysis , 96, 47–68.

20.

Hsiao

(2002) Analysis of Panel Data , 2nd edition. Cambridge, UK: Cambridge University Press.

21.

Ishwaran

, James

(2001) Gibbs sampling methods for stick-breaking priors. Journal of American Statistical Association , 96, 161–73.

22.

Juarez

, Steel

MFJ

(2010) Model-based clustering of non-Gaussian panel data based on skew-t distributions. Journal of Business and Economic Statistics , 28, 52–66.

23.

Kazemi

, Crouchley

(2006) Modelling the initial conditions in dynamic regression models of panel data with random effects. In Baltagi, BH, ed. Panel Data Econometrics, Theoretical Contributions and Empirical Applications . Amsterdam: Elsevier. pages 91–117.

24.

Kazemi

, Davies

(2002) The asymptotic bias of MLEs for dynamic panel data models. In Stasinopoulos M and Touloumi G, eds, Proceedings of the 17th International Workshop on Statistical Modelling, Greece, Chania, pages 391–95.

25.

, Lin

, Muller

(2010) Bayesian in- ference in semiparametric mixed models for longitudinal data. Biometrics , 66, 70–78.

26.

, Tian

, Xiao

, Chen

(2015) Variance change-point detection in panel data models. Economics Letters , 126, 140–43.

27.

Lunn

, Spiegelhalter

, Thomas

, Best

(2009) The BUGS project: Evolution, critique and future directions (with discussion). Statistics in Medicine , 28, 3049–82.

28.

MacEachern

(1994) Estimating normal means with a conjugate style Dirichlet process prior. Communications in Statistics , 23, 727–41.

29.

The MathWorks, Inc. (2012) MATLAB and Statistics Toolbox Release . Natick, MA: Author.

30.

McLachlan

, Peel

(2000) Finite Mixture Models . New York, NY: Wiley.

31.

McNicholas

(2016) Mixture Model-based Classification . Boca Raton, Florida: Chapman and Hall/CRC.

32.

Morgan

, Hodge

, Trepinksi

, Anderson

(2014) The stability of teacher performance and effectiveness: Implications for policies concerning teacher evaluation. Education Policy Analysis Archives , 22, 1–18.

33.

Morse

(1993) Does variability increase with age? An archival study of cognitive measures. Psychology and Aging , 8, 156–64.

34.

Nagin

(1999) Analyzing developmental trajec- tories: A semiparametric, group-based app- roach. Psychological Methods , 4, 139–57.

35.

Nagin

(2005) Group-based Modeling of Development . Cambridge, MA: Harvard University Press.

36.

Nie

, Chen

, Zhang

, Guo

(2010) Credit card customer analysis based on panel data clustering. Procedia Computer Science , 1, 2489–97.

37.

Nigg

(2001) Explaining adolescent exercise behaviour change: A longitudinal application of the transtheoretical model. Annals of Behavioral Medicine , 23, 11–20.

38.

Pennell

, Dunson

(2007) Fitting semiparametric random effects models to large datasets. Biostatistics , 8, 821–34.

39.

Rathod

, Yadav

, Rathore

(2011) Probabilistic modelling of fatigue damage accumulation for reliability prediction. International Journal of Quality, Statistics, and Reliability , 2011, 10. DOI:10.1155/2011/718901.

40.

Rikhtehgaran

, Kazemi

(2016) The determination of uncertainty levels in robust clustering of subjects with longitudinal observations using the Dirichlet process mixture. Advances in Data Analysis and Classification , 10, 541–62. DOI:10.1007/s11634-016-0262-x

41.

Sethuraman

(1994) A constructive definition of Dirichlet priors. Statistica Sinica 4639–50.

42.

Suarez

, Ghosal

(2016) Bayesian clustering of functional data using local features. Bayesian Analysis , 11, 71–98.

43.

Theodore

, Broadbent

, Nagin

, Ambler

, Hogan

, Ramrakha

, Cutfield

, Williams

MJA

, Harrington

, Moffitt

, Caspi

, Milne

, Poulton

(2015) Childhood to Early-midlife systolic blood pressure trajectories: Early-life predictors, effect modifiers, and adult cardiovascular outcomes. Hypertension , 66, 1108–15.

44.

Tobia

, Inauen

(2010) Gathering time-series data for evaluating behavior-change campaigns in developing countries: Reactivity of diaries and interviews. Evaluation Review , 34, 367–90.

45.

Vermunt

(2010) Longitudinal research using mixture models. In vanMontfort K, Oud JHL and Satorra A eds. Longitudinal Research with Latent Variables , Chapter 4, Heidelberg: Springer, pages 119–152.

46.

Vogt

, Linton

(2016) Classification of non-parametric regression functions in longitudinal data models. Journal of the Royal Statistical Society: Series B . DOI: 10.1111/rssb.12155

47.

Wang

, Coit

(2007) Reliability and degradation modelling with random or uncertain failure threshold. In Proceeding of the Annual Reliability and Maintainability Symposium, IEEE Reliability Society, pages 392–97 Orlando, FL.

48.

Wang

, Wang

(2013) Hierarchical Dirichlet process model for gene expression clustering. EURASIP Journal on Bioinformatics and Systems Biology , 2013, 5.

49.

Weisburd

, Bushway

, Lum

, Yang

(2004) Trajectories of crime at places: A longitudinal study of street segments in the city of seattle. Criminology , 42, 283–22.

50.

, Wunsch

(2005) Survey of clustering algorithms. IEEE Transactions on Neural Networks , 16, 645–78.

51.

, Zhong

P-S

, Wang

(2016) Detecting variance change-points for blocked time series and dependent panel data. Journal of Business and Economic Statistics , 34, 213–26.