Bayesian Detection of Piecewise Linear Trends in Replicated Time-Series with Application to Growth Data Modelling

2.1 Prior assumptions

For the mean parameters we assume that

independent for n = 1, …, N; t = 1, …, T. The quantities ν0>0 and μ0t∈R; t = 1, …, T, correspond to fixed hyper-parameters. The following default values are considered: μ0t=1NR∑n=1N∑r=1Rxntr, for t = 1, …, T. It is suggested that the parameter ν0 should be sufficiently small so that the prior distribution in (7) has large variability around the global mean. This is particularly important in cases where the multiple time series of an experiment exhibit strong heterogeneity. Values between 5×10−2⩽ν0⩽5×10−1 performed reasonably well in our setup. The variance σnt2 may be fixed (see Section 3.2) or not (see Section 3.3).

In order to specify the prior distribution of locations for a given number of change-points (ℓn), we are taking into account the prior assumption that later time-points are more likely to contain changes than earlier time-points. The growth levels of ℓn consecutive time-points during the end of the observation period are more likely to break collinearity than earlier stages. For example, all time-series consist of an initial period with very small and almost constant growth level in which no change is expected to occur. On the other hand, the last part of the observation period may exhibit larger heterogeneity. A simple and efficient way to incorporate such a prior information while supporting (with positive probability) all possible configurations with respect to constraint (2) is the following.

Let g(i;ii,i2)=1i2−i1+1I(i1⩽i⩽i2) denote the probability mass function of the discrete uniform distribution defined over the finite set of integers i such that i1⩽i⩽i2, where I(⋅) denotes the indicator function. For j = 1 we assume that τn1∼g(⋅;2,T−ℓn). For j⩾2 and conditionally on the event (τn1=t1,…,τn;j−1=tj−1), we assume that τnj follows a (discrete) uniform prior distribution defined over tj−1+1⩽τnj⩽T−ℓn+j−1. Thus, the prior distribution f(τn|ℓn>0) for a specific realization (t1,…,tℓn) of τn is defined as

Note that according to eq. (8), P(τn1=T−ℓn,τn2=T−ℓn+1,…,τnℓn=T−1|ℓn)≈1/T, while P(τn1=2,τn2=3,…,τnℓn=ℓn+1|ℓn)≈1/Tℓn, which satisfies our prior expectation discussed in the previous paragraph. We assume a-priori independence of τn for n = 1, …, N. Obviously, τn makes sense only in the case that the total number of change-points is strictly positive, thus eq. (8) is defined conditionally on the event ℓn>0.

Finally, the prior distribution of the number of change-points should be defined. Recall that the distribution in eq. (1) assigns a distinct mean parameter per time-point. However, given ℓn, only a small subset of θ’s will influence the likelihood in eq. (4) and the rest of them will affect the posterior solely due to their contribution to the prior distribution. Our motivation is based on the fact that we are trying to find very minimal models with few change points because it makes the data easier to interpret.

For a given number of change-points (ℓ), the number of possible models is (T∗ℓ), a number which grows rapidly with ℓ, where T∗=T−2. In order to penalize model complexity we consider prior distributions which favour sparse configurations. A natural way to incorporate such information on our prior distribution is to assume that the prior probability of ℓ change-points is inversely proportional to (T∗ℓ), that is, the number of models of size ℓ. This approach is also justified by the work of [10], where they consider that the number of change-points follows an exponentially decreasing prior distribution.

The prior P(⋅) has exponential decrease if, for some constants C > 0 and D < 1,

where ℓ∗ denotes the true value of ℓn. In the context of multivariate normal mean models with an underlying sparse true mean vector, it has been shown [10] that asymptotically, priors satisfying (9) lead to posterior distributions that concentrate on the sparse underlying true generative model. Members of the family of the so-called “complexity priors”, defined as,

have exponential decrease (9) for b>1+e [10]. In our applications we consider the choices b = 3.72 and α = 2 as default values, but we also report various prior sensitivity checks in the Appendix. As noted by Castillo and van der Vaart [10] it holds that eℓlog⁡(T∗/ℓ)⩽(T∗ℓ)⩽eℓlog⁡(eT∗/ℓ) implying that (10) is inversely proportional to the number of models of size ℓ. Thus, this choice is suited to the purpose of penalizing model complexity. Note that the right hand side of eq. (10) for ℓ=0 should be perceived as the limit limℓ↓0e−αℓlog⁡(bT∗/ℓ)=1.

Finally, we mention that typical prior assumptions on the number of change-points (for example a truncated Poisson or a uniform distribution over a pre-specified set of non-negative integer values) tend to overfit the number of change-points. This behaviour is demonstrated in Section 4.1 using simulated data with a known number of change-points, as well as in Section C of the Appendix using our real dataset.

3.1 Metropolis–Hastings MCMC Sampler

Assume first that the variances σ={σnt2;n=1,…,N,t=1,…,T} are given. This assumption will be relaxed at the end of this section. Observe that conditionally on the vector σ2, {(θn,ℓn,τn);n=1,…,N} are a-posteriori independent. Therefore, the inferential procedure breaks down to N independent tasks. The posterior distribution is written as

Although analytical evaluation of the marginal posterior distribution f(τ,ℓ|x,σ2)=∫Θf(θ,τ,ℓ|x,σ2)dθ is possible since we use conjugate prior assumptions, the integration cannot be carried out independently within each segment due to the continuity constraint. This fact makes the use of standard dynamic programming approaches not applicable in our set-up and the discrete nature of the sampling problem will make the computation of the involved expressions a time consuming task since the possible combinations of change-points increases rapidly with T. Thus, in order to make inference we sample from the target posterior distribution using a Metropolis-Hastings MCMC sampler that updates (θ,ℓ,τ).

At each step, the state of the chain is updated using four move-types: move 1 updates the number of change-points, move 2 updates the mean parameters by using a random walk proposal centered at the current values, move 3 updates the position of change-points and move 4 updates the subset of mean parameters that are not allocated to a change-point. In each case the proposed move is accepted according to the usual Metropolis-Hastings acceptance ratio, that is, Rn=min{1,αn} where

In eq. (13), (θn,ℓn,τn) denotes the current state of the n-th chain and (θn′,ℓn′,τn′) denotes the candidate state. Moreover, we use the notation Pprop(x→y) to denote the probability of proposing state y when the current state of the chain is x.

Move 1 This move updates the number of change-points, while keeping the mean parameters constant. We introduce two move types which propose to update the total number of change-points by 1. These move-types are complementary to each other: addition/deletion of a change-point. In the following, {τn∪t} denotes the resulting ordered set when a new change-point t is added to the current configuration τn. In a similar fashion, {τn∖t} denotes the remaining set when a specific member t of τn is removed from the c urrent configuration.

At a given state consisting of ℓn change-points, we propose addition/deletion with probabilities pa(ℓn) and pd(ℓn)=1−pa(ℓn), respectively. The addition probabilities are defined as

where L denotes the maximum number of change-points (L⩽T−2). In case of addition, we propose to add a randomly drawn change-point between two successive ones. The probability of proposing the addition of change-point t∗ such that τn;j−1<t∗<τnj; for some j=1,…,ℓn+1, is equal to pa(ℓn)1τnj−τn;j−1−11ℓn+1. In case that τnj−τn;j−1=1 the proposed move is immediately rejected. In the reverse move, a the previously added change-point is selected with probability pd(ℓn+1)1ℓn+1 and is deleted from τn∪t∗. Thus, the acceptace probability for an addition move is equal to

In the case of proposing deletion of a change-point τnj, the corresponding acceptance ratio term is equal to

At this point we underline that using an overfitted set of model parameters (one mean θnt per time-point t = 1, …, T) allows us to use the standard Metropolis-Hastings ratio for proposing additions/deletions of change-points. This would not be true if the number of mean parameters was defined conditionally on ℓn: in such a case the Reversible Jump algorithm or integration of mean parameters is required.

Move 2 In this move the update of mean parameters θ is proposed, while all other parameters remain unchanged. For this purpose a random walk centered at the current values of the chain is used. For subject n = 1, …, N, let θnt denote the current value of the mean parameters at time-point t = 1, …, T. Then a new state is proposed according to

independent for all t and n, for some constant c > 0. Recall that σnt2 denotes the variance of the random sample (Xnt1,…,XntR) which is assumed known. Note that the ratio of the proposal distribution for the transitions θn→θn′ and θn′→θn is 1. Thus, the Metropolis-Hastings acceptance ratio (13) simplifies to

Move 3.a The candidate state is generated by using a proposal distribution which will jointly update the change-points τn, while the total number of change-points ℓn and mean parameters θn are kept constant. Let ε=(ε1,…,εℓn) and

where d1>0 denotes a pre-specified positive integer and g(⋅;−d1,d1) denotes the discrete uniform distribution over {−d1,−d1+1,…,d1−1,d1}. Then, the proposed state is generated as

independently for i=1,…,ℓn, while ℓn′=ℓn as well as θn′=θn. In this case, the proposal ratio in eq. (13) is equal to 1 so the acceptance ratio is written as the posterior probability ratio.

where θτn as in (5). A small value for d1 will be capable of achieving optimal acceptance rates. Thus, this move is oriented towards the local exploration of the posterior surface, given the current state.

Move 3.b This is a similar proposal to Move 3.b, but instead of proposing the simultaneous update of all cut-points, just one entry is modified and the rest remain the same. As in Move 3.a, both the number of change-points as well as the values of mean parameters remain the same. Thus, let i∗ denote a randomly drawn index from the set {1,…,ℓn} and

where d2>0 denotes a pre-specified positive integer. The proposed state is generated as

The Metropolis-Hastings acceptance probability simplifies to eq. (18). In this case, a sufficiently large value for d2 will propose moves that are more likely to be accepted compared to Move 3.a, since only one entry is changed. Thus, move 3.b will be used as complementary to move 3.a, in order to facilitate the ability of escaping from local modes of the posterior distribution.

Move 4 Let θn[−τn] denote the mean of those time-points that do not correspond to change-points for time series n. In this case it can be easily seen that the full conditional distribution of (θn[−τn]|τn,xn,ℓn) is the prior distribution in eq. (7). Hence, a draw from the prior distribution will perfom a Gibbs update to θn[−τn].

Note that both moves 3.a and 3.b are able to propose states that have zero prior probability in eq. (8). Although this is not a frequent event, in this case the proposed state is immediately rejected, since the prior probability ratio is equal to zero.

Recall that the previous MCMC steps are defined conditionally on σ2. The next subsection deals with the case where the variance is estimated at a pre-processing stage and plugged into the previously described MCMC sampler. For this purpose, we consider that the variance can be either shared between different time series or not and two estimators are derived. We assume that we have enough data in order to obtain a robust estimate of variance, which is the case in our application. The more general case where the variance is treated as a random variable is discussed in Subsection 3.3. In this case, the variance is updated by the full conditional distribution using a Gibbs sampling step.

3.2 Variance estimation at a pre-processing stage

In this case the variance is considered known and in practice it should be estimated at a pre-processing stage. We use the posterior mean arising from a multivariate normal–inverse gamma model as a point estimate. For this purpose we ignore the piecewise linear parameterization of the mean function and use the same likelihood as in eq. (1) and the same prior assumptions for θnt as in eq. (7).

We will assume two parameterizations: the full model where the variances are a-priori distributed as:

independent for n = 1, …, N; t = 1, …, T, where IG(α,β) denotes the inverse Gamma distribution. The second model parameterization imposes the restriction of common variance across different time series and replicates, that is,

Under (21), a-priori it is assumed that

independent for t = 1, …, T. The quantities α0>0 and β0>0 correspond to fixed hyper-parameters.

Let us define the function

for t = 1, …, T; n = 1, …, N. It easily follows that under (20):

independent for n = 1, …, N; t = 1, …, N. In the case of the restricted parameterization in (21), the corresponding posterior distribution is

independent for t = 1, …, T. Then, we use the posterior means as the plug-in point estimates of the variance per time-point, that is,

provided that α0+R/2>1 so that both posterior expectations exist. The default values we use for the constants in eq. (22) and (23) are α0=1 and β0=1, but we also report some prior sensitivity checks in the simulation section. We will use the labels “s1” and “s2” to refer to the MCMC sampler using the plug-in estimates (22) and (23), respectively.

3.3 Updating the variance

Here we discuss the case where the variance in eq. (1) is unknown. In the case where all variances are unrestricted (eq. (20)), the full conditional distributions are:

independent for n = 1, …, N, t = 1, …, T. Note that in this case, the full conditional distribution depends on the piecewise linear mean function μ(t;θn,τn). Under restriction (21), it follows that:

independent for t = 1, …, T.

We will use the labels “s3” and “s4” to refer to the MCMC sampler using the Gibbs steps (20) and (21), respectively. Since the full conditional distribution of σnt in eq. (20) depends only on the quantities (xn,θn,τn), the joint posterior distribution factorizes over n = 1, …, N. Thus, we can split the MCMC sampling into N independent samplers, as also done to the case where the variance is fixed. However, this does not hold for eq. (21) where the state of the parameters μ(t;θn,τn) of all time-series (n = 1, …, N) is required. This imposes a large computational burden since N can be quite large. Therefore, in the examples of the next section we have only applied the first three MCMC samplers (s1, s2 and s3).

4.1 Simulation study

We considered simulated datasets of length T = 1000 time-points consisting of N = 1000 independent multivariate observations, while the number of replicates (R) is equal to R = 3. For n = 1, …, N, the number of change-points ℓn is drawn uniformly at random from the set {0, 1, …, 9}. A detailed description of the simulation mechanism is given in Section B of the Appendix.

All three samplers were applied using the following set of hyper-parameter values: α0=0.1, β0=0.1, ν0=0.005. The prior distribution on the number of change-points corresponds to the complexity prior distribution in (10), with parameters α = 2 and b = 3.72. The fixed-variance samplers (s1 and s2) are initialized from a state with 1 change-point with parameters randomly generated from the prior distribution. On the other hand, we observed the sampler s3 remains trapped in areas of low posterior probability when using random starting values. In order to deal with this issue, we initialized sampler s3 using a run of sampler s1 (with 30000 iterations). The inference is based on 50000 iterations following a burn-in period of 20000 iterations.

Our results are benchmarked against the not package available in the Comprehensive R Archive Network, which implements the narrowest-over-threshold method of Baranowski et al. [22]. The number of change-points is inferred using a strengthened version of the Bayesian Information Criterion [26, 38]. For the problem of detecting changes in the slope, this approach only considers univariate time-series with constant variance, thus, it is not suitable for multivariate time-series where the variance varies with time (which is the case in our datasets). In order to apply this method we averaged across the replicates of the time-series.

Figure 2 displays the difference of the estimated number of change-points (ℓ^) from the true number (ℓ), stratified according to the true number of change-points used to generate the data. We conclude that in all cases the estimates are centered in the true value of the number of change-points. At the lower panel of Figure 2 we calculated the Mean Absolute Error (MAE) of the resulting estimates, with error bars corresponding to the standard error. For smaller number of change-points the MAE is consistently higher for sampler s1. Although the simulation scenario does not assume the same variance per time-series, note that the sampler s2 (which uses the pooled variance estimate in eq. (21)) gives slightly better results.

Figure 2:

Benchmarking the estimation of the number of change-points arising from the narrowest-over-threshold (“not”) method of [22] and our MCMC sampler with fixed different variance (s1), fixed shared variance (s2) and unknown different variance (s3) per time-series, under the complexity prior distribution. The number after the name of each sampler indicates the true number of change-points and for each value 100 synthetic datasets were simulated. Each time-series consists of three replicates and 1000 time-points.

We conclude that in many cases the method of [22] tends to overfit the number of change-points, resulting to worse performance than our method. Clearly, benchmarking against this simpler approach should not be perceived as a fair comparison between the two methods, but as a means of highlighting the usefulness of our modelling.

Figure 3 displays the output of the MCMC sampler s1 (lower panel) and s2 (upper panel) for a time-series where the true number of change-points equals to 8. The output of sampler s3 is almost identical to the lower panel of Figure 3. The posterior distribution of the number of change-points is shown at the left panels, considering both the complexity prior distribution (gray trace) as well as a Poisson(1) distribution truncated on the set {0, 1, …, 30}. Observe that the first choice quickly converges to a state with 8 change-points (that is, the true number). This is not the case under the Poisson prior distribution, which supports larger values than the true number of change-points. This behaviour is typical in other simulated datasets we tried, therefore, all results in the remaining sections are based on the complexity prior. The right panels display the posterior distribution of change-point locations (under the complexity prior) and it evident that the method is able to accurately infer all change-point locations.

Figure 3:

Example of a simulated time-series with 8 change-points. Replicates are shown in red, blue and green color. Left: MCMC trace of the sampled number of change-points considering both the complexity and a truncated Poisson(1) prior distributions (every 20th iteration is displayed). Right: output of the MCMC sampler conditionally on the estimated MAP of changepoints (which is equal to 8) according to the complexity prior. The upper and lower panels correspond to the MCMC sampler using the variance estimates arising from eq. (23) (sampler s2) and (22) (sampler s1), respectively. The gray lines correspond to the posterior mean estimates of the piecewise linear mean function and the boxplots display the posterior marginal distribution of each change-point. Vertical dotted lines correspond to the central line of each boxplot (median). The coloured outer regions correspond to two estimated standard deviations from the mean.

Additional results based on synthetic data are provided in Section B of the Appendix, including prior sensitivity checks as well as alternative simulation scenarios. Further comparisons against not are provided in Section D of the Appendix.

4.2 Phase detection in parallel time-series analysis of fungal growth

The filamentous fungal pathogen Aspergillus fumigatus is a major pathogen of the human lung causing more deaths per annum than tuberculosis or malaria [39]. A time series study of fungal growth was performed in liquid culture by analysing, in parallel, the growth characteristics of 411 independent transcription factor gene deletion mutants. The mutant strains were cultivated in a microtiter plate containing 200 µL of a fungal culture medium and incubated at 37 °C. Optical density (at 600 nm) was measured at 10 minute intervals for a total period of 48 hours. The growth analysis was performed on three separate occasions.

The observed data consists of N × R × T growth levels for R = 3 replicates of N = 411 objects (mutants) measured every 10 minutes for T = 289 time-points. Figure 1 displays the observed time series for four mutants. Visual inspection reveals that describing growth with a piecewise linear mean function with an unknown number of segments is a reasonable assumption for the observed data. Regarding the fixed hyper-parameter values, we considered that α = 2 (eq. (10)) and ν0=10−1 (eq. (7) and (23)). After estimating the variance per time-point using the estimator in (23), the MCMC sampler s2 ran for m = 50000 iterations, following a burn-in period of 20000.

The boxplots in Figure 1 correspond to the estimate of the marginal posterior distribution of each change-point for specific subset of four mutants, conditionally on the mode of the posterior distribution of the number of change-points. Figure 4 displays the averaged profile per mutant (mean of three replicates) coloured according to the most probable number of change-points for each of N = 411 subjects. We conclude that the majority of the mutants (343) consist of three growth phases. It is clear that mutants with a smaller number of change-points are also the more slowly growing mutants, which is reasonable since these mutants most likely have not been able to reach the later growth phases in the time of the experiment. In particular the method inferred 35 mutants with only 2 growth phases and slow growth behaviour while 12 mutants have a single phase and very slow growth behaviour. Finally, 21 mutants consist of 4 growth phases and some of them exhibit a faster growth rate at later observation stages (t > 220).

Figure 4:

Visualization of the growth dataset with respect to the estimated MAP number of phases. For plotting convenience, each curve corresponds to the average growth time series across the three replicates.

Amongst 12 fungal mutants identified as having a single change-point during growth curve analysis, and therefore exhibiting severely retarded growth kinetics, seven mutants had previously been characterised [40, 41, 42, 43, 44, 45]. Without exception previously characterised mutants had been reported as having various morphological defects, and these morphogenesis mutants were correctly identified in our analysis. The remaining five mutants have, until now, remained uncharacterised and therefore provide promising candidates to investigate further in order to establish their roles in fungal morphogenesis. Detailed phenotypical characterisations of the identified mutants will be described elsewhere.

a

We have also considered a Poisson(1) prior distribution on the number of change-points. As already demonstrated in Section 4.1, in this case the sampler selects a larger number of change-points which are less interpretable. The reader is referred to Section C of the Appendix.