# Simultaneous Generation of Optimum Pavement Clusters and Associated Performance Models

Authors: Mukesh Khadka; Alexander Paz; Cristian Arteaga; David K. Hale
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

With regard to developing pavement performance models (PPMs), the existing state-of-the-art proposes Clusterwise Linear Regression (CLR) to determine the pavement clusters and associated PPMs simultaneously. However, the approach does not determine optimal clustering to minimize error; that is, the number of clusters and explanatory variables are prespecified to determine the corresponding coefficients of the PPMs. In addition, existing formulations do no address issues associated with overfitting as there is no limit to include parameters in the model. In order to address this limitation, this paper proposes a mathematical program within the CLR approach to determine simultaneously(1) an optimal number of clusters, (2) assignment of segments into clusters, and (3) regression coefficients for all prespecified explanatory variables required to minimize the estimation error. The Bayesian Information Criteria is proposed to limit the number of optimal clusters. A simulated annealing coupled with ordinary least squares was used to solve the mathematical program.


---

## Body

## 1. Introduction

Typically, pavement performance models (PPMs) are developed using a two-step approach. First, pavement segments with similar characteristics are grouped into clusters using a few critical factors, such as pavement type, age, and traffic volume. Then, PPMs for each of the clusters are developed using statistical techniques. The objective of clustering is to group the pavement segments that perform similarly over time. However, in practice, the performances of pavement segments within a cluster differ significantly because clusters are formed using only a few critical factors [1, 2]. A major challenge is to select characteristics that define clusters and the corresponding segments associated with them [3]. If inappropriate characteristics are used, clusters may include homogeneous segments with different performance behavior or heterogeneous segments with similar performance behavior [4]. The prediction accuracy of PPMs can be improved by subdividing the pavement segments into more uniform clusters. However, this subdivision is not always possible due to limited information [1].Figures1(a) and 1(b) provide an example of heterogeneous performance behavior for two segments, each grouped within the same cluster (the Prioritization Category), using the two-step approach and actual data from the Pavement Management System (PMS) of the Nevada Department of Transportation (NDOT). Segments “SR445 (SB), MP: 40-39” and “SR445 (NB), MP: 36-37” were assigned to one cluster, Prioritization Category 4. Segments “SR156 (WB), MP: 3-2” and “SR892 (SB), MP: 25-24” were assigned into Prioritization Category 5.Figure 1
(a) and (b) Heterogeneous performance behavior of pavement segments from the same Prioritization Category; (c) and (d) potential clusters with pavement segments having homogeneous performance behavior. NB: north bound; SB: south bound; WB: west bound; MP: milepost.
(a) (b) (c) (d)In contrast, Figure1(c) illustrates that segments “SR445 (SB), MP: 40-39” and “SR156 (WB), MP: 3-2” had homogeneous performance behavior. Similarly, segments “SR892 (SB), MP: 25-24” and “SR445 (NB), MP: 36-37” showed a consistent behavior (see Figure 1(d)). This suggests that factors other than the Prioritization Category are critical in causing the differences in performance behavior. Influencing factors could include subgrade type, traffic loading characteristics, or any hidden effects.To address the limitations of the two-step approach, Spath [5] proposed using Clusterwise Linear Regression (CLR) to determine clusters and associated regression models simultaneously. CLR assigns pavement segments that have similar regression effects on the pavement performance such that the overall sum of squared errors is minimal. Hence, CLR minimizes the overall prediction error by simultaneously determining the explanatory variables’ coefficients and assigning each pavement segment into an appropriate cluster.In the context of pavement management, CLR first was used by Luo and Chou [1] to model the deterioration of pavement conditions. First, the pavement network was clustered by using a few critical pavement characteristics. The subdivision was performed at the data-point level; that is, data points collected over various years for a pavement segment could be assigned to multiple clusters. Hence, there was a chance of pavement segments being associated with multiple performance models. An additional step was proposed to predict performance using the results from multiple models. Later, Luo and Yin [2] expanded their research, using CLR to formulate the development of pavement distresses. Both studies [1, 2] used pavement age as the only explanatory variable.To address some of the limitations present in the studies completed by Luo and Chou [1] and Luo and Yin [2], Zhang and Durango-Cohen [6] developed CLR models using multiple explanatory variables. Data used in the study were collected during the AASHO Road Test [7], conducted in the late 1950s in Ottawa, Illinois. The data were collected from a single site in a relatively controlled environment; clearly, the site characteristics were not representative of all other locations. For example, the data did not represent the varieties of soil and climatic conditions. In addition, construction techniques and materials have changed substantially since this test was conducted. In this study, the number of clusters was determined manually using an “elbow” criterion. Experiments were run only for a number of clusters equal to 1, 5, 10, 15, 20, and 25. One possible reason for not examining all the possible number of clusters might have been large amount of required computational time given that an exchange algorithm was utilized with 100 instances with various initial assignments.The previous CLR methods used to estimate PPMs do not test the explanatory power of variables in both cluster and regression analyses. That is, all the explanatory variables used in regression models are assumed to be significant. The effects of any insignificant explanatory variables on the dependent variable are also accounted for during assignment of pavement segments into clusters and estimation of the regression coefficients. The presence of any insignificant explanatory variables distorts the underlying regression effects of other significant explanatory variables. This may lead to an incorrect assignment of pavement segments into clusters and estimation of PPMs. Khadka and Paz [8] developed PPMs using CLR while testing the significance of the explanatory variables and obtained superior results relative to the case when significance was assumed.Another key limitation of the existing CLR is the need to prespecify the number of clusters. In order to avoid prespecifying the number of required clusters, this study proposes a mathematical program to simultaneously determine an optimal number of clusters, the assignment of segments into clusters, and regression coefficients for all explanatory variables. This mathematical program is flexible enough to handle multiple explanatory variables, multiple observations per pavement segments, and user-defined constraints on cluster characteristics.In previous studies using CLR for pavement management [1, 2, 6], the objective function was to minimize the sum of squared errors of prediction (SSE). It is intuitive that SSE decreases monotonically as the number of clusters increases. That is, for a given dataset, the optimum number of clusters always is the total number of data points [9]. An optimum number of clusters are equal to the total number of pavement segments, and each pavement segment is the sole member of its own cluster. Such clustering structure is unlikely to provide statistical reliable models. In addition, SSE always decreases when a new explanatory variable is added to the model [10]. Usually, this leads to an overfitting problem [11]. Therefore, SSE is not the best objective function to use for searching an optimal number of clusters.Even though SSE decreases as the number of clusters increases, the rate of improvement diminishes significantly after a relative optimal [12] number of clusters known as the elbow point. Increasing the number of clusters beyond the elbow point provides a very small reduction in SSE. An SSE versus the number of cluster curve might not exhibit an elbow point distinctly in all cases. Hence, it could be very challenging to choose the right number of clusters.To address these limitations, this study extended the existing CLR framework by simultaneously (i) incorporating the Bayesian Information Criteria (BIC) [13] as the objective function, (ii) finding the optimal number of clusters, and (iii) finding the maximum number of clusters as required for model estimation. The BIC penalizes more for the inclusion of additional parameters relative to the Akaike Information Criteria (AIC) [14]. The BIC selects simpler models than the AIC when the sample size is greater than eight [15]. Hence, the optimal number of clusters found using BIC provides a balance between model complexity and goodness of fit [16]. On the other hand, several studies showed that the number of parameters in a model selected using AIC was overestimated [14, 17–19]. The search process for the best model specification using BIC has the property of consistency, which asymptotically selects this model with a probability of “1” [20–23]. Galimberti et al. [24] used BIC to propose a unified framework for model-based clustering, linear regression, and multiple cluster structure detection.In this study, the data limitations in the existing literature were addressed using actual field data collected across a variety of environmental, traffic, design, construction, and maintenance conditions. Pavement data collected for the past 12 years over the entire State of Nevada were used. This data included significant variations across a large range of characteristics, e.g., pavement segments exposed to either extreme desert heat or to very low winter temperatures in the mountains.


## 2. Methodology

### 2.1. Problem Formulation

This section includes the definition of a pavement sample, notation and terms, proposed mathematical program, and a procedure to find upper bound of the range of the feasible number of clusters.


#### 2.1.1. Definition of a Pavement Sample

The condition of a pavement segment improves when intervention occurs by applying an M&R treatment. Such intervention alters the physical characteristics of the pavement. Hence, the performance of a pavement before and after the intervention differs, even though all other contributing factors remain constant. In this circumstance, the same pavement segment before and after intervention should be treated as two different samples. Given that the physical location of a pavement segment is the same, a different identifier is required to distinguish the set of consecutive observations before and after the intervention. In this study, the term pavement sample is used as an identifier to uniquely represent the set of consecutive observations that accounts for historical interventions made on a pavement segment. Figure2 provides a simplified depiction of how consecutive observations of a pavement segment are divided into two pavement samples.Figure 2
A typical pavement performance curve and a simplified depiction of how the observations of a pavement segments are divided into two samples.In this study, a pavement sample, instead of pavement segment, was considered as a single entity during cluster analysis. Hence, if a pavement segment consists of two or more pavement samples, these samples could be assigned to different clusters.


#### 2.1.2. Notation and Terms

The following notation and terms are used in describing the proposed model:i = Subscript for a pavement sample in the networkI = Number of pavement samples in the network, indexed 1≤i≤It = Subscript for an observation period for a pavement sampleTi = Number of observation periods for a pavement sample i, indexed 1≤t≤TiO = Total number of observations = ∑iITi∀i∈Ij = Subscript for an explanatory variableJ = Number of explanatory variables, indexed 1≤j≤Jk = Subscript for a clusterKmax = Maximum number of clustersK = Number of optimum clustersxijtk = Measurement of an explanatory variable j(1≤j≤J) for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)yitk = Measurement of dependent variable for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)n = Minimum number of observations required in a clusterCk = Set of pavement samples that are assigned to cluster k(1≤k≤Kmax)pik = Cluster membership of a pavement sample i to a cluster k,∀1≤i≤I, and 1≤k≤Kmax.δk = Intercept for regression model in cluster k∀1≤k≤Kmaxβjk = Slope coefficient for explanatory variable j(1≤j≤J) in cluster k(1≤k≤Kmax)


#### 2.1.3. Mathematical Program

A major aspect of model development is to identify variables that explain an actual physical process. Pavement performance models (PPMs) should only include explanatory variables that affect pavement deterioration. In this study, explanatory variables from the existing literature were selected. For example, pavement age, traffic and environmental conditions, and structural and material properties of pavement were considered as factors affecting pavement performance [25–28].PSI was chosen as the dependent variable,y. PSI serves as a unified standard and is widely accepted for evaluating pavement performance, especially in terms of ride quality [29–31]. In addition, PSI reflects the human rider’s response and is understood by highway users and legislators [32]. The adopted functional form for the regression model is expressed as(1)yit=β0k+∑j=1Jβjk∗xijtThis study proposes a mixed-integer, nonlinear mathematical program to determine an optimal number of clusters, assignment of segments into clusters, and regression coefficients for all prespecified explanatory variables. The problem was defined by the optimum number of clusters,K; the number of predefined explanatory variables, J; the number of pavement samples to be clustered, I; and the number of observation periods, Ti, associated with each pavement sample. The formulation partitions pavement samples into an optimum number of clusters, with a PPM model fit to each cluster.The objective function involved minimization of the BIC acrossK clusters as expressed by (2)Min.BIC=O+O∗ln⁡2π+O∗ln⁡SSEO+JK+2K-1∗ln⁡OGiven that BIC is an increasing function of SSE and free parameters (the number of clusters and regression coefficients) to be estimated, a clustering with the lowest BIC provides the optimal solution. Minimizing BIC reduces unexplained variation in the dependent variable [13]. The total SSE was calculated using(3)SSE=∑k=1K∑i=1I∑t=1Tiδk+∑j=1Jβjk∗xijtk-yitk2∗pik∀i,j,t,kEach cluster was associated with a linear regression model with predefined explanatory variables. Deviations of predicted statistics from actual data were calculated separately for each cluster and summed to obtain the total SSE.Decision variables to be determined were the coefficients for all prespecified explanatory variables,δk and βjk; the optimum number of clusters, K; and the cluster membership, pik.The following constraints were imposed to describe the proposed problem:(4)∑kpik=1∀i∈i,k(5)pik=1,ifsampleiisassignedtoclusterk;0,otherwise∀i,k(6)∑i∈CkTi≥n∀Ck(7)1≤k≤Kmax(8)Kmax=FI,Ti,nEquations (4) and (5) ensure that each pavement sample was assigned to exactly one cluster. The indicator pik equals 1 if and only if a pavement sample i belongs to cluster k. Otherwise, it takes a value of zero.Equation (6) describes a minimum-size constraint, which is imposed to ensure sufficient observations in each cluster for a statistically reliable estimation of coefficients. The total number of observations in a cluster is required to be no less than the minimum number of observations, n.Equations (7) and (8) are imposed to determine the maximum number of potential clusters. If a pavement sample has more than n observations, regression over these observations could generate statistically reliable estimates of the coefficients. Hence, a cluster could be formed with only one pavement sample that has more than n observations. If all pavement samples have more than n observations, each pavement sample could form a cluster. In this case, the maximum number of clusters would be the total number of pavement samples, I. However, in reality, it is possible that none of the pavement samples would have more than n observations. In this case, samples would be grouped to form a cluster at the sample level, but not at the observation level; that is, observations of a sample must not be assigned to more than one cluster. Equation (8) determines the maximum number of potential clusters in both cases. Function F in this constraint represents the following algorithm to determine Kmax. A flowchart for this algorithm is provided in Figure 3.Figure 3
Algorithm utilized to calculate the maximum number of potential clusters.


#### 2.1.4. Proposed Algorithm to DetermineKmax

Step 1.
If the total number of observations,O, is less than the minimum number of observations required to form a cluster, n, then set Kmax= zero and go to Step 6. Otherwise, create a matrix, M, of size (τmax x 2) with the following elements, where τmax is maximum number of observations of a pavement segment(s) in the dataset:(a)
The first column ofM includes all integer values from 1 to τmax.(b)
The second column includes the number of segments associated with the number of observations.(c)
If no segments have a particular number of observations in the dataset, then set the second column of the matrix to zero.Step 2.
If any segment has a number of observations greater than or equal ton(mτ,1≥n), then (a)
calculateKmax=∑τ≥nmτ,2(b)
updatemτ,2 with 0 for τ≥n
Otherwise, go to Step 3 to find the maximum number of clusters that could be formed.Step 3.
If the matrixM has all zeros in its second column ∑τmτ,2=0, then return Kmax=Kmax and go to Step 6. Otherwise,(a)
updateM by removing all rows that have the number of segments equal to zero (mτ,2=0)(b)
initialize two indices asω = ϑ = number of rows inM(c)
make a copy ofM and let it be represented by M′(d)
if the remaining total number of observations(∑τ=1ωmτ,1∗mτ,2) is less than n, then, Kmax=Kmax, and go to Step 6. Otherwise, initialize S with the value of mω,1 and mω,2=mω,2-1Step 4.
Repeat the following steps untilS=n.Step  4.1. If (mϑ,2=0), then ϑ=ϑ-1. Otherwise, go to Step 4.3.Step  4.2. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setn=n+1(iv)
setS=0(v)
go to Step 4.6Otherwise, go to Step 4.3.Step  4.3. If (S>n), then(i)
S=S-mϑ,1(ii)
mϑ,2=mϑ,2+1(iii)
Setϑ=ϑ-1Otherwise, go to Step 4.6.Step  4.4. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setω=ω-1(iv)
setS=0Otherwise, go to Step 4.6.Step  4.5. If (ω=0), then(i)
update both indicesω and ϑ with the number of rows of M(ii)
setn=n+1(iii)
go to Step 4.6Otherwise,(i)
setS=mω,1(ii)
setmω,2=mω,2-1(iii)
go to Step 4.6Step  4.6. Update S with S+mϑ,1 and mϑ,2 with mϑ,2-1.Step 5.
SetKmax=Kmax+1, and go to Step 3.Step 6.
Return the current value ofKmax and stop.


### 2.2. Solution to the Mathematical Program

Simulated annealing (SA) coupled with an ordinary least square (OLS) algorithm was implemented to solve the above mathematical program. SA was used to cluster the dataset estimate,pik. For each accepted neighborhood cluster, OLS was utilized to estimate the regression coefficients, δk and βjk. The fitting linear models (lm) function, available in the statistical software, R, was used to estimate these coefficients [33]. DeSarbo et al. [34] successfully implemented such an algorithm to solve the CLR problem.The algorithm utilized to solve the clusterwise multiple linear regression is described as follows, and illustrated in Figure4.Figure 4
Algorithm utilized to solve the clusterwise multiple linear regression.Step 1.
Initialization:
Step  1.1. Set K=2, and BICmin=infinity.
Step  1.2. Set values of initial temperature (θ0), final minimum temperature (θmin), cooling rate (λ), and the maximum number of neighbors to be generated (Nmax) at each temperature level. Set the iterator N=1.Step 2.
Calculate maximum number of potential clusters,Kmax, utilizing function F as described above, as part of Constraint 8.Step 3.
Initial estimation of regression coefficients:
Step  3.1. For a given number of clusters, K, randomly assign cluster memberships to all pavement samples.
Step  3.2. Count the number of observations of all pavement samples assigned to each cluster. If all clusters have at least n observations, then go to Step 4; otherwise, reassign the cluster memberships until all clusters have at least n observations. Let CKN be the valid initial clusters.
Step  3.3. Estimate δk and βjk for all K clusters using OLS.Step 4.
Evaluate objective function,BICKN using (2).Step 5.
Generate a set of neighborhood clusters near to the previous one, using the following steps.
Step  5.1. Randomly select a prespecified number of pavement samples (Nps) to change their memberships.
Step  5.2. For each of the samples selected, assign a new membership by generating a random number u~U1,K. If the new membership is same as the previous one, regenerate a random number u′~U1,K until it is different. Repeat this process until the memberships of all the selected pavement samples are different from those that were previously assigned.
Step  5.3. Count the total number of observations of all pavement samples assigned to each cluster.
Step  5.4. If all clusters have at least n observations, go to Step 6; otherwise, repeat Steps 5.1., 5.2., and 5.3. until all clusters have at least n observations. Let CKN+1 be a new set of valid neighborhood clusters.Step 6.
Search for a solution:
Step  6.1. For CK+1N, estimate new δk and βjk for all K clusters using OLS.
Step  6.2. Evaluate BICKN+1 using (1).
Step  6.3. Calculate ΔBIC=BICKN+1-BICKN.
Step  6.4. Check the following two conditions:(a)
IfΔBIC<0, accept the current set of clusters, CN+1K, and corresponding δk and βjk. Go to Step 7; otherwise, go to Step (b).(b)
Generate a random numberu′′~U0,1. Calculate the acceptance probability, paccept=exp-ΔBIC/B∗T, where B is a Boltzmann’s constant. If paccept>u′′, accept the current set of clusters, CK+1N, and corresponding δk and βjk. Go to Step 7; otherwise, return to Step 5.Step 7.
Counter and temperature update:
Step  7.1. Repeat Steps 5 and 6 for Nmax times.
Step  7.2. If θ<θmin, stop the algorithm. Otherwise, reduce the temperature by multiplying the current temperature by the prespecified cooling rate, λ, set N=1, and go to Step 5.Step 8.
Stopping criteria:
Step  8.1. Update BICmin with the smallest between the one obtained in Step 7 and the current BICmin. Set Koptimal equal to K.
Step  8.2. Repeat Steps 3 to 7 for Kmax-1 times.This algorithm seeks solutions using a probabilistic approach. The algorithm starts with a high temperature,θ, and a high probability of accepting a worse solution, paccept. This enables occasional “uphill” moves, which help escape from the local minima. The algorithm builds up a rough view of the search space by moving with large step lengths. As θ drops, paccept decreases to behave more closely as a greedy algorithm, with small step lengths slowly focusing on the most promising solution space. Theoretical studies have shown that, with infinitely slow cooling, the algorithm converges to a global minimum [35].


## 2.1. Problem Formulation

This section includes the definition of a pavement sample, notation and terms, proposed mathematical program, and a procedure to find upper bound of the range of the feasible number of clusters.


### 2.1.1. Definition of a Pavement Sample

The condition of a pavement segment improves when intervention occurs by applying an M&R treatment. Such intervention alters the physical characteristics of the pavement. Hence, the performance of a pavement before and after the intervention differs, even though all other contributing factors remain constant. In this circumstance, the same pavement segment before and after intervention should be treated as two different samples. Given that the physical location of a pavement segment is the same, a different identifier is required to distinguish the set of consecutive observations before and after the intervention. In this study, the term pavement sample is used as an identifier to uniquely represent the set of consecutive observations that accounts for historical interventions made on a pavement segment. Figure2 provides a simplified depiction of how consecutive observations of a pavement segment are divided into two pavement samples.Figure 2
A typical pavement performance curve and a simplified depiction of how the observations of a pavement segments are divided into two samples.In this study, a pavement sample, instead of pavement segment, was considered as a single entity during cluster analysis. Hence, if a pavement segment consists of two or more pavement samples, these samples could be assigned to different clusters.


### 2.1.2. Notation and Terms

The following notation and terms are used in describing the proposed model:i = Subscript for a pavement sample in the networkI = Number of pavement samples in the network, indexed 1≤i≤It = Subscript for an observation period for a pavement sampleTi = Number of observation periods for a pavement sample i, indexed 1≤t≤TiO = Total number of observations = ∑iITi∀i∈Ij = Subscript for an explanatory variableJ = Number of explanatory variables, indexed 1≤j≤Jk = Subscript for a clusterKmax = Maximum number of clustersK = Number of optimum clustersxijtk = Measurement of an explanatory variable j(1≤j≤J) for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)yitk = Measurement of dependent variable for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)n = Minimum number of observations required in a clusterCk = Set of pavement samples that are assigned to cluster k(1≤k≤Kmax)pik = Cluster membership of a pavement sample i to a cluster k,∀1≤i≤I, and 1≤k≤Kmax.δk = Intercept for regression model in cluster k∀1≤k≤Kmaxβjk = Slope coefficient for explanatory variable j(1≤j≤J) in cluster k(1≤k≤Kmax)


### 2.1.3. Mathematical Program

A major aspect of model development is to identify variables that explain an actual physical process. Pavement performance models (PPMs) should only include explanatory variables that affect pavement deterioration. In this study, explanatory variables from the existing literature were selected. For example, pavement age, traffic and environmental conditions, and structural and material properties of pavement were considered as factors affecting pavement performance [25–28].PSI was chosen as the dependent variable,y. PSI serves as a unified standard and is widely accepted for evaluating pavement performance, especially in terms of ride quality [29–31]. In addition, PSI reflects the human rider’s response and is understood by highway users and legislators [32]. The adopted functional form for the regression model is expressed as(1)yit=β0k+∑j=1Jβjk∗xijtThis study proposes a mixed-integer, nonlinear mathematical program to determine an optimal number of clusters, assignment of segments into clusters, and regression coefficients for all prespecified explanatory variables. The problem was defined by the optimum number of clusters,K; the number of predefined explanatory variables, J; the number of pavement samples to be clustered, I; and the number of observation periods, Ti, associated with each pavement sample. The formulation partitions pavement samples into an optimum number of clusters, with a PPM model fit to each cluster.The objective function involved minimization of the BIC acrossK clusters as expressed by (2)Min.BIC=O+O∗ln⁡2π+O∗ln⁡SSEO+JK+2K-1∗ln⁡OGiven that BIC is an increasing function of SSE and free parameters (the number of clusters and regression coefficients) to be estimated, a clustering with the lowest BIC provides the optimal solution. Minimizing BIC reduces unexplained variation in the dependent variable [13]. The total SSE was calculated using(3)SSE=∑k=1K∑i=1I∑t=1Tiδk+∑j=1Jβjk∗xijtk-yitk2∗pik∀i,j,t,kEach cluster was associated with a linear regression model with predefined explanatory variables. Deviations of predicted statistics from actual data were calculated separately for each cluster and summed to obtain the total SSE.Decision variables to be determined were the coefficients for all prespecified explanatory variables,δk and βjk; the optimum number of clusters, K; and the cluster membership, pik.The following constraints were imposed to describe the proposed problem:(4)∑kpik=1∀i∈i,k(5)pik=1,ifsampleiisassignedtoclusterk;0,otherwise∀i,k(6)∑i∈CkTi≥n∀Ck(7)1≤k≤Kmax(8)Kmax=FI,Ti,nEquations (4) and (5) ensure that each pavement sample was assigned to exactly one cluster. The indicator pik equals 1 if and only if a pavement sample i belongs to cluster k. Otherwise, it takes a value of zero.Equation (6) describes a minimum-size constraint, which is imposed to ensure sufficient observations in each cluster for a statistically reliable estimation of coefficients. The total number of observations in a cluster is required to be no less than the minimum number of observations, n.Equations (7) and (8) are imposed to determine the maximum number of potential clusters. If a pavement sample has more than n observations, regression over these observations could generate statistically reliable estimates of the coefficients. Hence, a cluster could be formed with only one pavement sample that has more than n observations. If all pavement samples have more than n observations, each pavement sample could form a cluster. In this case, the maximum number of clusters would be the total number of pavement samples, I. However, in reality, it is possible that none of the pavement samples would have more than n observations. In this case, samples would be grouped to form a cluster at the sample level, but not at the observation level; that is, observations of a sample must not be assigned to more than one cluster. Equation (8) determines the maximum number of potential clusters in both cases. Function F in this constraint represents the following algorithm to determine Kmax. A flowchart for this algorithm is provided in Figure 3.Figure 3
Algorithm utilized to calculate the maximum number of potential clusters.


### 2.1.4. Proposed Algorithm to DetermineKmax

Step 1.
If the total number of observations,O, is less than the minimum number of observations required to form a cluster, n, then set Kmax= zero and go to Step 6. Otherwise, create a matrix, M, of size (τmax x 2) with the following elements, where τmax is maximum number of observations of a pavement segment(s) in the dataset:(a)
The first column ofM includes all integer values from 1 to τmax.(b)
The second column includes the number of segments associated with the number of observations.(c)
If no segments have a particular number of observations in the dataset, then set the second column of the matrix to zero.Step 2.
If any segment has a number of observations greater than or equal ton(mτ,1≥n), then (a)
calculateKmax=∑τ≥nmτ,2(b)
updatemτ,2 with 0 for τ≥n
Otherwise, go to Step 3 to find the maximum number of clusters that could be formed.Step 3.
If the matrixM has all zeros in its second column ∑τmτ,2=0, then return Kmax=Kmax and go to Step 6. Otherwise,(a)
updateM by removing all rows that have the number of segments equal to zero (mτ,2=0)(b)
initialize two indices asω = ϑ = number of rows inM(c)
make a copy ofM and let it be represented by M′(d)
if the remaining total number of observations(∑τ=1ωmτ,1∗mτ,2) is less than n, then, Kmax=Kmax, and go to Step 6. Otherwise, initialize S with the value of mω,1 and mω,2=mω,2-1Step 4.
Repeat the following steps untilS=n.Step  4.1. If (mϑ,2=0), then ϑ=ϑ-1. Otherwise, go to Step 4.3.Step  4.2. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setn=n+1(iv)
setS=0(v)
go to Step 4.6Otherwise, go to Step 4.3.Step  4.3. If (S>n), then(i)
S=S-mϑ,1(ii)
mϑ,2=mϑ,2+1(iii)
Setϑ=ϑ-1Otherwise, go to Step 4.6.Step  4.4. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setω=ω-1(iv)
setS=0Otherwise, go to Step 4.6.Step  4.5. If (ω=0), then(i)
update both indicesω and ϑ with the number of rows of M(ii)
setn=n+1(iii)
go to Step 4.6Otherwise,(i)
setS=mω,1(ii)
setmω,2=mω,2-1(iii)
go to Step 4.6Step  4.6. Update S with S+mϑ,1 and mϑ,2 with mϑ,2-1.Step 5.
SetKmax=Kmax+1, and go to Step 3.Step 6.
Return the current value ofKmax and stop.


## 2.1.1. Definition of a Pavement Sample

The condition of a pavement segment improves when intervention occurs by applying an M&R treatment. Such intervention alters the physical characteristics of the pavement. Hence, the performance of a pavement before and after the intervention differs, even though all other contributing factors remain constant. In this circumstance, the same pavement segment before and after intervention should be treated as two different samples. Given that the physical location of a pavement segment is the same, a different identifier is required to distinguish the set of consecutive observations before and after the intervention. In this study, the term pavement sample is used as an identifier to uniquely represent the set of consecutive observations that accounts for historical interventions made on a pavement segment. Figure2 provides a simplified depiction of how consecutive observations of a pavement segment are divided into two pavement samples.Figure 2
A typical pavement performance curve and a simplified depiction of how the observations of a pavement segments are divided into two samples.In this study, a pavement sample, instead of pavement segment, was considered as a single entity during cluster analysis. Hence, if a pavement segment consists of two or more pavement samples, these samples could be assigned to different clusters.


## 2.1.2. Notation and Terms

The following notation and terms are used in describing the proposed model:i = Subscript for a pavement sample in the networkI = Number of pavement samples in the network, indexed 1≤i≤It = Subscript for an observation period for a pavement sampleTi = Number of observation periods for a pavement sample i, indexed 1≤t≤TiO = Total number of observations = ∑iITi∀i∈Ij = Subscript for an explanatory variableJ = Number of explanatory variables, indexed 1≤j≤Jk = Subscript for a clusterKmax = Maximum number of clustersK = Number of optimum clustersxijtk = Measurement of an explanatory variable j(1≤j≤J) for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)yitk = Measurement of dependent variable for a sample i(1≤i≤I) at observation period t(1≤t≤Ti) that is assigned to a cluster k(1≤k≤Kmax)n = Minimum number of observations required in a clusterCk = Set of pavement samples that are assigned to cluster k(1≤k≤Kmax)pik = Cluster membership of a pavement sample i to a cluster k,∀1≤i≤I, and 1≤k≤Kmax.δk = Intercept for regression model in cluster k∀1≤k≤Kmaxβjk = Slope coefficient for explanatory variable j(1≤j≤J) in cluster k(1≤k≤Kmax)


## 2.1.3. Mathematical Program

A major aspect of model development is to identify variables that explain an actual physical process. Pavement performance models (PPMs) should only include explanatory variables that affect pavement deterioration. In this study, explanatory variables from the existing literature were selected. For example, pavement age, traffic and environmental conditions, and structural and material properties of pavement were considered as factors affecting pavement performance [25–28].PSI was chosen as the dependent variable,y. PSI serves as a unified standard and is widely accepted for evaluating pavement performance, especially in terms of ride quality [29–31]. In addition, PSI reflects the human rider’s response and is understood by highway users and legislators [32]. The adopted functional form for the regression model is expressed as(1)yit=β0k+∑j=1Jβjk∗xijtThis study proposes a mixed-integer, nonlinear mathematical program to determine an optimal number of clusters, assignment of segments into clusters, and regression coefficients for all prespecified explanatory variables. The problem was defined by the optimum number of clusters,K; the number of predefined explanatory variables, J; the number of pavement samples to be clustered, I; and the number of observation periods, Ti, associated with each pavement sample. The formulation partitions pavement samples into an optimum number of clusters, with a PPM model fit to each cluster.The objective function involved minimization of the BIC acrossK clusters as expressed by (2)Min.BIC=O+O∗ln⁡2π+O∗ln⁡SSEO+JK+2K-1∗ln⁡OGiven that BIC is an increasing function of SSE and free parameters (the number of clusters and regression coefficients) to be estimated, a clustering with the lowest BIC provides the optimal solution. Minimizing BIC reduces unexplained variation in the dependent variable [13]. The total SSE was calculated using(3)SSE=∑k=1K∑i=1I∑t=1Tiδk+∑j=1Jβjk∗xijtk-yitk2∗pik∀i,j,t,kEach cluster was associated with a linear regression model with predefined explanatory variables. Deviations of predicted statistics from actual data were calculated separately for each cluster and summed to obtain the total SSE.Decision variables to be determined were the coefficients for all prespecified explanatory variables,δk and βjk; the optimum number of clusters, K; and the cluster membership, pik.The following constraints were imposed to describe the proposed problem:(4)∑kpik=1∀i∈i,k(5)pik=1,ifsampleiisassignedtoclusterk;0,otherwise∀i,k(6)∑i∈CkTi≥n∀Ck(7)1≤k≤Kmax(8)Kmax=FI,Ti,nEquations (4) and (5) ensure that each pavement sample was assigned to exactly one cluster. The indicator pik equals 1 if and only if a pavement sample i belongs to cluster k. Otherwise, it takes a value of zero.Equation (6) describes a minimum-size constraint, which is imposed to ensure sufficient observations in each cluster for a statistically reliable estimation of coefficients. The total number of observations in a cluster is required to be no less than the minimum number of observations, n.Equations (7) and (8) are imposed to determine the maximum number of potential clusters. If a pavement sample has more than n observations, regression over these observations could generate statistically reliable estimates of the coefficients. Hence, a cluster could be formed with only one pavement sample that has more than n observations. If all pavement samples have more than n observations, each pavement sample could form a cluster. In this case, the maximum number of clusters would be the total number of pavement samples, I. However, in reality, it is possible that none of the pavement samples would have more than n observations. In this case, samples would be grouped to form a cluster at the sample level, but not at the observation level; that is, observations of a sample must not be assigned to more than one cluster. Equation (8) determines the maximum number of potential clusters in both cases. Function F in this constraint represents the following algorithm to determine Kmax. A flowchart for this algorithm is provided in Figure 3.Figure 3
Algorithm utilized to calculate the maximum number of potential clusters.


## 2.1.4. Proposed Algorithm to DetermineKmax

Step 1.
If the total number of observations,O, is less than the minimum number of observations required to form a cluster, n, then set Kmax= zero and go to Step 6. Otherwise, create a matrix, M, of size (τmax x 2) with the following elements, where τmax is maximum number of observations of a pavement segment(s) in the dataset:(a)
The first column ofM includes all integer values from 1 to τmax.(b)
The second column includes the number of segments associated with the number of observations.(c)
If no segments have a particular number of observations in the dataset, then set the second column of the matrix to zero.Step 2.
If any segment has a number of observations greater than or equal ton(mτ,1≥n), then (a)
calculateKmax=∑τ≥nmτ,2(b)
updatemτ,2 with 0 for τ≥n
Otherwise, go to Step 3 to find the maximum number of clusters that could be formed.Step 3.
If the matrixM has all zeros in its second column ∑τmτ,2=0, then return Kmax=Kmax and go to Step 6. Otherwise,(a)
updateM by removing all rows that have the number of segments equal to zero (mτ,2=0)(b)
initialize two indices asω = ϑ = number of rows inM(c)
make a copy ofM and let it be represented by M′(d)
if the remaining total number of observations(∑τ=1ωmτ,1∗mτ,2) is less than n, then, Kmax=Kmax, and go to Step 6. Otherwise, initialize S with the value of mω,1 and mω,2=mω,2-1Step 4.
Repeat the following steps untilS=n.Step  4.1. If (mϑ,2=0), then ϑ=ϑ-1. Otherwise, go to Step 4.3.Step  4.2. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setn=n+1(iv)
setS=0(v)
go to Step 4.6Otherwise, go to Step 4.3.Step  4.3. If (S>n), then(i)
S=S-mϑ,1(ii)
mϑ,2=mϑ,2+1(iii)
Setϑ=ϑ-1Otherwise, go to Step 4.6.Step  4.4. If (ϑ = 0), then(i)
M=M′(ii)
setϑ = number of rows ofM(iii)
setω=ω-1(iv)
setS=0Otherwise, go to Step 4.6.Step  4.5. If (ω=0), then(i)
update both indicesω and ϑ with the number of rows of M(ii)
setn=n+1(iii)
go to Step 4.6Otherwise,(i)
setS=mω,1(ii)
setmω,2=mω,2-1(iii)
go to Step 4.6Step  4.6. Update S with S+mϑ,1 and mϑ,2 with mϑ,2-1.Step 5.
SetKmax=Kmax+1, and go to Step 3.Step 6.
Return the current value ofKmax and stop.


## 2.2. Solution to the Mathematical Program

Simulated annealing (SA) coupled with an ordinary least square (OLS) algorithm was implemented to solve the above mathematical program. SA was used to cluster the dataset estimate,pik. For each accepted neighborhood cluster, OLS was utilized to estimate the regression coefficients, δk and βjk. The fitting linear models (lm) function, available in the statistical software, R, was used to estimate these coefficients [33]. DeSarbo et al. [34] successfully implemented such an algorithm to solve the CLR problem.The algorithm utilized to solve the clusterwise multiple linear regression is described as follows, and illustrated in Figure4.Figure 4
Algorithm utilized to solve the clusterwise multiple linear regression.Step 1.
Initialization:
Step  1.1. Set K=2, and BICmin=infinity.
Step  1.2. Set values of initial temperature (θ0), final minimum temperature (θmin), cooling rate (λ), and the maximum number of neighbors to be generated (Nmax) at each temperature level. Set the iterator N=1.Step 2.
Calculate maximum number of potential clusters,Kmax, utilizing function F as described above, as part of Constraint 8.Step 3.
Initial estimation of regression coefficients:
Step  3.1. For a given number of clusters, K, randomly assign cluster memberships to all pavement samples.
Step  3.2. Count the number of observations of all pavement samples assigned to each cluster. If all clusters have at least n observations, then go to Step 4; otherwise, reassign the cluster memberships until all clusters have at least n observations. Let CKN be the valid initial clusters.
Step  3.3. Estimate δk and βjk for all K clusters using OLS.Step 4.
Evaluate objective function,BICKN using (2).Step 5.
Generate a set of neighborhood clusters near to the previous one, using the following steps.
Step  5.1. Randomly select a prespecified number of pavement samples (Nps) to change their memberships.
Step  5.2. For each of the samples selected, assign a new membership by generating a random number u~U1,K. If the new membership is same as the previous one, regenerate a random number u′~U1,K until it is different. Repeat this process until the memberships of all the selected pavement samples are different from those that were previously assigned.
Step  5.3. Count the total number of observations of all pavement samples assigned to each cluster.
Step  5.4. If all clusters have at least n observations, go to Step 6; otherwise, repeat Steps 5.1., 5.2., and 5.3. until all clusters have at least n observations. Let CKN+1 be a new set of valid neighborhood clusters.Step 6.
Search for a solution:
Step  6.1. For CK+1N, estimate new δk and βjk for all K clusters using OLS.
Step  6.2. Evaluate BICKN+1 using (1).
Step  6.3. Calculate ΔBIC=BICKN+1-BICKN.
Step  6.4. Check the following two conditions:(a)
IfΔBIC<0, accept the current set of clusters, CN+1K, and corresponding δk and βjk. Go to Step 7; otherwise, go to Step (b).(b)
Generate a random numberu′′~U0,1. Calculate the acceptance probability, paccept=exp-ΔBIC/B∗T, where B is a Boltzmann’s constant. If paccept>u′′, accept the current set of clusters, CK+1N, and corresponding δk and βjk. Go to Step 7; otherwise, return to Step 5.Step 7.
Counter and temperature update:
Step  7.1. Repeat Steps 5 and 6 for Nmax times.
Step  7.2. If θ<θmin, stop the algorithm. Otherwise, reduce the temperature by multiplying the current temperature by the prespecified cooling rate, λ, set N=1, and go to Step 5.Step 8.
Stopping criteria:
Step  8.1. Update BICmin with the smallest between the one obtained in Step 7 and the current BICmin. Set Koptimal equal to K.
Step  8.2. Repeat Steps 3 to 7 for Kmax-1 times.This algorithm seeks solutions using a probabilistic approach. The algorithm starts with a high temperature,θ, and a high probability of accepting a worse solution, paccept. This enables occasional “uphill” moves, which help escape from the local minima. The algorithm builds up a rough view of the search space by moving with large step lengths. As θ drops, paccept decreases to behave more closely as a greedy algorithm, with small step lengths slowly focusing on the most promising solution space. Theoretical studies have shown that, with infinitely slow cooling, the algorithm converges to a global minimum [35].


## 3. Experiment and Results

### 3.1. Data Resources

Data used in this study were extracted from the PMS database of NDOT [36]. The data consisted of various classes, location data, segment data, contract data, environmental data, traffic data, and pavement condition data, collected for flexible pavement in the entire State of Nevada.A detailed data analysis was performed to check for inconsistent and missing information in the dataset, and some were found. Some of the missing data were synthesized based on associated information available in the dataset. In preparing the PMS data, the following filters were applied:(i)
Only one-mile segments were selected for consistency.(ii)
Only pavement segments with the most recent maintenance contracts awarded in 2001 or later were used in the study.(iii)
PSI of a pavement should deteriorate over time if no M&R treatment occurred. If PSI of a segment in any year increased by 0.1 or more points from the previous year without any M&R treatment, all observations for that year were excluded from the analysis. However, if an increase in PSI in any year was less than 0.1 from the previous year, it was assumed to be a random error during the process of pavement evaluation or data processing. Therefore, those observations were included in the analysis.(iv)
If the PSI of any year decreased by one or more points from the previous year, all observations for that year were excluded from the analysis.(v)
In practice, the PSI range is between 4.5 and 1.5. Therefore, if a pavement segment had a PSI beyond these limits in any year, it was considered an outlier, and all observations for that year were excluded.(vi)
Only PSI values used were within the interval of the mean, minus three standard deviations to the mean plus three standard deviations.(vii)
Pavement samples that did not consist of data regarding conditions for at least two consecutive years were excluded.(viii)
Data analysis showed that an improvement in PSI was seen one or two years after the contract award date. Hence, the age of the pavement sample was set to 0 when the actual improvement occurred rather than when the contract was awarded.After data preparation was completed, 4,138 flexible pavement samples with 17,642 observations were available. For CLR modelling, 14,637 observations, collected from 2001 to 2010, were used; the remaining 3,005 observations, collected in 2011 and 2012 (about 17% of the total number), were used as test dataset to check the accuracy of the CLR models. Tables1 and 2 provide the descriptive statistics of continuous and categorical variables used in this paper, respectively.Table 1
Descriptive statistics for the continuous variables.
Variable Minimum Maximum Mean Std. deviation psi 1.60 4.57 4.01 0.41 age 0.00 8.00 2.24 2.01 adt 20.00 132000.00 4844.45 9812.57 trucks 1.00 7731.00 862.29 1082.20 elevation 228.60 2667.00 1368.25 415.19 precip 3.94 89.28 19.33 10.10 min_temp -6.67 13.33 3.20 4.00 max_temp 7.78 31.67 20.31 4.22 wet_days 11.00 81.00 42.14 15.67 freeze_thaw 0.00 230.00 136.75 51.51 rut_depth 0.00 1.60 0.14 0.14Table 2
Descriptive statistics for the continuous variables.
Variable Category Dummy variable Number of observations Percent System ID IR - 4,622 31.6 NHS nhs 5,063 34.6 STP stp 4,952 33.8 Number of Lanes 1 - 8,450 57.7 2 lane=2 5,612 38.3 ≥ 3 lane≥3 575 3.9 Prioritization Category 1 - 5,004 34.2 2 category =2 3,181 21.7 3 category =3 3,118 21.3 4 category =4 1,558 10.6 5 category =5 1,776 12.1 Functional Class 1 - 4,695 32.1 2 f_class = 2 107 0.7 3 f_class = 3 5,106 34.9 4 f_class = 4 2,604 17.8 5 f_class = 5 1,830 12.5 6 f_class = 6 251 1.7 7 f_class = 7 44 0.3 Note: IR, interstate route; NHS, national highway system; STP, Surface Transportation Program.


### 3.2. Parameters of the Algorithm

Performance of the SA algorithm generally depends on the values of the optimization parameters utilized for a given problem. To ensure proper initialization and search for optimal solutions, selection of the most appropriate parameter values is critical [37, 38]. A body of literature exists regarding various methodologies for finding the most appropriate values for annealing parameters in SA [37, 39–43].If an SA algorithm is allowed to run for a sufficiently long time by setting a high initial temperature with a slow cooling rate, the algorithm performs well, as shown by Anily and Federgruen [44]. In such a cooling scheme, the selection of the most appropriate parameter values may not be critical. However, computation time cannot always be ignored. Hence, the algorithm has to find a good solution in a reasonable amount of time [39].Effective values to be assigned to the optimization parameters depend on the type and complexity of the problem. These values may not be obvious to determine, but rather might be determined by trial and error for a given problem [40]. In this study, values assigned to the optimization parameters were determined using experience gained from previous research [45–51] that involved SA and other comparable algorithms. Table 3 lists the parameter values used in this study.Table 3
Setup parameters for implementation of the proposed algorithm.
Parameter Value Remarks θ 0 10 Initial temperature θ m i n 10e-17 Final minimum temperature B 3000 Boltzmann constant λ 0.97 Cooling rate N m a x 10 Number of neighborhood solutions generated at each temperature level n 800 Minimum number of observations required in a cluster N p s 25 Number of pavement samples, which memberships were changed to generate a neighborhood clusterThe minimum number of observations required in a cluster,n, is one of the parameters to be defined by the analyst for each dataset and application as it is required for any other statistical analysis regarding the minimum sample size. That is, the proposed methodology and contributions are not restricted or affected by the sample size. The analyst must have available sufficient data for each cluster to be able to obtain reliable estimates. A too small n will result in statistically unreliable models and a potentially time-consuming search process as the maximum number of feasible clusters; Kmax will be very large. In contrast, a too large n will result in insufficient number of clusters required to provide the optimum goodness of fit. Hence, sensitivity analysis is required to achieve balance between reliable estimates and Kmax.


### 3.3. Results and Discussion

Given the constraints for feasible partitions defined in the problem formulation and the minimum number of observations required in a cluster,n=800, the proposed algorithm determined 16 as the maximum number of potential clusters. The algorithm searched for the optimum number of clusters from 2 to 16. Seven-cluster CLR models provided the optimum solution with the lowest BIC. The estimated regression coefficients for the CLR models are presented in Table 4.Table 4
Estimated model parameters using the proposed CLR approach.
Parameters C 1  (2,279) C 2  (1,959) C 3  (2,169) C 4  (2,094) C 5  (1,883) C 6  (1,936) C 7  (2,317) Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. intercept 4.145 0.312 6.242 0.369 2.7910 0.311 6.7280 0.358 7.7810 0.394 12.1400 0.401 3.8730 0.331 age -0.035 0.002 -0.040 0.003 -0.0350 0.003 -0.0327 0.003 -0.0464 0.003 -0.0392 0.004 -0.0498 0.003 adt† -0.006 <0.001 -0.004 <0.001 -0.0262 <0.001 -0.0028 <0.001 -0.0078 <0.001 -0.0053 <0.001 -0.0334 <0.001 t r u c k s † 0.0002‡ <0.001 0.0205‡ <0.001 0.0190‡ <0.001 -0.0151‡ <0.001 0.0306 <0.001 0.0557 <0.001 0.0752 <0.001 e l e v a t i o n † 0.0066‡ <0.001 0.0182‡ <0.001 -0.0352‡ <0.001 -0.1079 <0.001 -0.1418 <0.001 -0.0131‡ <0.001 0.1060 <0.001 precip -0.0037‡ 0.008 -0.0118‡ 0.009 -0.0037‡ 0.008 -0.0248 0.009 0.0094‡ 0.011 -0.0518 0.012 -0.0252 0.008 min_temp -0.030 0.009 0.0129‡ 0.010 -0.0092‡ 0.009 0.0497 0.010 -0.0321 0.011 -0.0447 0.013 0.0532 0.009 max_temp 0.025 0.007 -0.030 0.008 0.0249 0.007 -0.0568 0.008 -0.0124‡ 0.009 -0.0554 0.010 -0.0311 0.007 wet_days 0.005 0.002 -0.010 0.002 0.0115 0.002 0.0031‡ 0.002 -0.0028‡ 0.002 0.0004‡ 0.003 -0.0061 0.002 freeze_ t h a w † -1.697 <0.001 1.513 <0.001 -0.2029‡ <0.001 1.8550 <0.001 -1.4100 0.001 -13.7900 0.001 4.2370 <0.001 rut_depth -0.632 0.115 -1.002 0.115 -1.1060 0.095 -0.5614 0.117 -0.8408 0.119 -0.3999 0.143 -0.9900 0.099 lane=2 -0.371 0.029 -0.166 0.027 0.0121‡ 0.027 -0.1216 0.028 -0.5574 0.032 -0.2745 0.032 0.0258‡ 0.026 lane≥3 -0.325 0.044 -0.195 0.046 -0.1713 0.052 -0.3215 0.045 -0.3974 0.052 -0.2511 0.054 0.0391‡ 0.046 nhs -0.442 0.079 -0.407 0.172 0.7454 0.138 -0.2811 0.072 -0.2639 0.150 -0.4300 0.085 1.4320 0.138 stp -0.802 0.094 -0.487 0.180 0.6121 0.153 -0.3452 0.091 -0.1346‡ 0.131 -0.3379 0.104 1.1350 0.151 f_class=2 0.521 0.115 0.0940‡ 0.165 -0.9383 0.136 0.4118 0.091 0.7050 0.140 0.9474 0.115 -1.5090 0.149 f_class=3 0.382 0.100 0.3377‡ 0.174 -0.8551 0.139 0.3270 0.080 0.3107 0.148 0.4025 0.101 -1.3770 0.136 f_class=4 0.619 0.112 0.2964‡ 0.181 -0.6925 0.151 0.2945 0.094 0.0189‡ 0.137 0.5468 0.117 -1.0330 0.146 f_class=5 0.618 0.113 -0.0703‡ 0.182 -0.6471 0.152 -0.7763 0.098 -0.6250 0.143 0.3015 0.119 -1.5590 0.148 f_class=6 0.472 0.118 -0.419 0.188 -1.2170 0.159 -0.6566 0.108 -0.7955 0.145 0.0582‡ 0.124 -1.7880 0.154 f_class=7 0.554 0.204 0.416 0.199 -1.3840 0.173 -0.1375‡ 0.171 NA 0.151 0.7614 0.144 -1.3640 0.169 category=2 -0.306 0.070 -0.0444‡ 0.051 0.0762‡ 0.056 -0.1300 0.045 -0.6077 0.051 -0.1175‡ 0.065 0.0246‡ 0.035 category=3 -0.319 0.073 -0.0156‡ 0.056 -0.0421‡ 0.060 -0.2457 0.051 -0.6000 0.056 -0.2317 0.069 -0.0080‡ 0.040 category=4 -0.468 0.077 -0.298 0.064 -0.3665 0.065 -0.1535 0.062 -0.6331 0.066 -0.5050 0.077 -0.3982 0.052 category=5 -0.463 0.077 -0.356 0.064 -0.7165 0.065 -0.1899 0.063 -0.8446 0.067 -0.6001 0.076 -0.4145 0.052  BIC 136 338 238 216 496 857 271 Note: the quantity included in parentheses represents the total number of observations in a cluster. † variable value in thousands. ‡ coefficient with p value > 0.05. NA = not applicable.Figure5(a) shows the smallest BIC for each of the clusters (K=2 to 16) considered in this experiment. Figure 5(b) shows the trajectory of the objective function, BIC, when the CLR models were used. The initial value of BIC was 8,502. After 1,360 iterations, the BIC decreased to 3,008. This change was equivalent to an improvement of 65%.Figure 5
BIC trend over the number of clusters (a); trajectory of the BIC during optimization for seven-cluster models (b); SSE trend over the number of clusters (c); and trajectory of the SSE during optimization for eleven-cluster models (d).
(a) (b) (c) (d)It was observed that not all coefficients had associated p values less than 0.05. In this study, the significance level was considered to be 5%. As expected, coefficients differed in magnitude and sign across the clusters, which indicated that the deterioration patterns of pavement samples varied among the clusters. However, seven explanatory variables had the same sign across all clusters.Different clusters had different numbers of significant explanatory variables. For example, Cluster 2 had 10 insignificant explanatory variables. In addition, among all seven clusters, five variables,age,adt,rut_depth,category=4, andcategory=5, were significant. However, four variables,trucks,elevation,precip, andcategory=2, were not significant in four different clusters.Trucks play a key role in pavement deterioration because they transfer heavy loads to the pavement [52]. Hence, it is expected that variable “trucks” be significant. There could be various reasons precluding the statistical significance of this variable. Two of them include (i) multicollinearity effects, which are not addressed in the existing CLR literature and motivates the expansion of the framework as recommend in this paper under future research and (ii) different lanes mistakenly used to collect pavement performance and truck traffic data. The data used in this study does not include the lane used to collect the information.The performance of the proposed CLR approach was compared with that of the existing CLR for pavement management. That is, in this section models estimated using the proposed approach were compared with those estimated using the existing CLR [6] described in the Introduction. Experiments using the existing CLR approach were run for all feasible clusters (K=2 to 16). Figure 5(c) shows the smallest SEE for each of these clusters. As expected, SSE decreased with an increasing number of clusters, but at a very small rate after K=11. In this case, Figure 5(c) does not exhibit a clear elbow point. Hence, an optimum number of clusters needed to be decided by visual inspection while considering the trade-off between goodness of fit and model complexity (i.e., of the number of models and the explanatory variables). This inherent subjectivity when choosing an optimum number of clusters is a major drawback for the existing state-of-the-art CLR approach.After careful assessment, 11-cluster CLR models were selected as the optimum solution. Figure5(d) shows the trajectory of SSE, when 11-cluster CLR models were used, and Table 5 provides the corresponding regression coefficients. Similar to the results obtained from the proposed CLR approach, the coefficients differed in terms of magnitude and sign. In addition, some coefficients had p values larger than 0.05.Table 5
Coefficients obtained using the existing state-of-the-art CLR approach.
Parameters C 1  (1,229) C 2  (1,365) C 3  (1,413) C 4  (1,091) C 5  (1,423) C 6  (1,412) C 7  (1,430) C 8  (1,321) C 9  (1,272) C 10  (1,360) C 11  (1,321) Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. intercept 3.92 0.40 4.12 0.40 4.95 0.44 14.43 0.60 7.41 0.45 3.70 0.36 4.61 0.40 3.80 0.47 7.95 0.48 6.50 0.42 7.54 0.43 age -0.04 <0.01 -0.04 <0.01 -0.05 <0.01 -0.04 <0.01 -0.05 <0.01 -0.04 <0.01 -0.04 <0.01 -0.05 <0.01 -0.05 <0.01 -0.04 <0.01 -0.03 <0.01 adt † -0.04 <0.01 -0.02 <0.01 -0.02 <0.01 0.01 <0.01 -0.04 <0.01 0.00 <0.01 -0.01‡ <0.01 -0.04 <0.01 -0.01 <0.01 -0.01‡ <0.01 0.00 <0.01 trucks † 0.11 0.02 0.03 0.02 0.04 0.02 0.13 0.02 0.05 0.02 0.01‡ <0.01 -0.03‡ <0.01 -0.01‡ 0.02 0.05 <0.01 0.01‡ 0.02 0.03 <0.01 elevation † -0.02‡ <0.01 0.18 <0.01 -0.05‡ <0.01 -0.05‡ <0.01 -0.22 <0.01 0.01‡ <0.01 -0.11 <0.01 -0.17 <0.01 0.18 <0.01 -0.05‡ <0.01 0.01‡ <0.01 precip -0.01‡ <0.01 -0.02‡ <0.01 -0.01‡ <0.01 -0.13 <0.01 -0.01‡ <0.01 -0.01‡ <0.01 0.02‡ <0.01 0.04 <0.01 -0.05 <0.01 -0.03 <0.01 -0.05 <0.01 min_temp 0.01‡ <0.01 -0.02‡ <0.01 0.09 <0.01 -0.09 0.02 0.01‡ <0.01 -0.02 <0.01 -0.03 <0.01 -0.02‡ <0.01 0.12 <0.01 -0.01‡ <0.01 0.05 <0.01 max_temp -0.01‡ <0.01 0.01‡ <0.01 -0.06 <0.01 -0.05 <0.01 -0.03 <0.01 0.02 <0.01 0.02‡ <0.01 0.03 <0.01 -0.12 <0.01 -0.01‡ <0.01 -0.06 <0.01 wet_days 0.01‡ <0.01 -0.01 <0.01 -0.01‡ <0.01 0.02 <0.01 0.01‡ <0.01 0.01 <0.01 -0.01‡ <0.01 0.01‡ <0.01 -0.02 <0.01 0.01‡ <0.01 0.01 <0.01 freeze_thaw † 1.00‡ 0.56 -0.86‡ 0.50 6.59 0.55 -20.89 0.85 0.24‡ 0.59 -1.38 0.44 0.78‡ 0.52 1.45 0.59 5.22 0.71 -2.19 0.59 -1.93 0.59 rut_depth -0.64 0.15 -1.00 0.14 -1.22 0.13 -1.03 0.17 -0.27 0.13 -0.46 0.13 -0.76 0.13 -1.02 0.13 -1.80 0.15 -0.88 0.14 -1.20 0.15 lane=2 0.02‡ 0.04 -0.16 0.03 -0.35 0.04 -0.10 0.04 0.03‡ 0.04 -0.06‡ 0.03 -0.26 0.03 -0.22 0.03 -0.53 0.04 -0.28 0.04 -0.26 0.04 lane≥3 0.22 0.06 -0.03‡ 0.06 -0.43 0.06 -0.31 0.06 0.31 0.08 -0.43 0.06 -0.48 0.05 0.01‡ 0.07 -0.33 0.06 -0.27 0.06 -0.39 0.06 nhs 1.39 0.15 0.28‡ 0.17 0.51 0.16 -0.12‡ 0.16 1.20 0.14 -0.25 0.08 -0.41 0.08 0.80 0.18 -0.70 0.15 -0.50 0.17 0.44 0.15 stp 1.33 0.17 0.18‡ 0.19 0.42 0.17 -0.27‡ 0.18 0.67 0.17 -0.38 0.11 -0.50 0.10 0.82 0.19 -0.64 0.15 -0.65 0.19 0.65 0.14 f_class=2 -1.42 0.23 -0.39‡ 0.25 -0.24‡ 0.17 0.44 0.18 -0.83 0.24 0.32 0.09 0.52 0.10 -0.93 0.17 1.43 0.20 -0.25‡ 0.17 -0.71 0.17 f_class=3 -1.36 0.18 -0.55 0.18 -0.38 0.15 0.12‡ 0.18 -0.88 0.15 0.15‡ 0.09 0.33 0.12 -0.99 0.17 0.93 0.14 0.43 0.16 -0.55 0.13 f_class=4 -1.25 0.18 -0.49 0.20 -0.39 0.16 0.22‡ 0.20 -0.56 0.18 0.27 0.11 0.41 0.13 -1.12 0.18 0.89 0.14 0.69 0.18 -0.59 0.12 f_class=5 -1.85 0.19 -0.51 0.20 -0.91 0.16 0.03‡ 0.20 -1.43 0.18 0.21‡ 0.11 0.18‡ 0.13 -2.18 0.18 0.60 0.14 0.43 0.18 -0.47 0.12 f_class=6 -2.07 0.20 -0.71 0.20 -0.77 0.17 -0.72 0.21 -1.71 0.19 0.14‡ 0.12 -0.66 0.14 -2.24 0.19 0.24‡ 0.15 0.56 0.19 -0.71 0.13 f_class=7 NA -0.63 0.22 -1.53 0.25 0.42‡ 0.29 -0.68 0.19 0.09‡ 0.16 0.18‡ 0.16 -2.28 0.25 NA 0.72 0.25 NA category=2 0.02‡ 0.09 0.07‡ 0.06 -0.36 0.05 0.22 0.10 -0.27 0.07 0.06‡ 0.05 -0.20 0.09 -0.07‡ 0.06 -0.66 0.06 -0.13 0.05 -0.14‡ 0.07 category=3 0.02‡ 0.10 0.01‡ 0.07 -0.43 0.06 0.14‡ 0.11 -0.17 0.08 -0.05‡ 0.06 -0.29 0.10 -0.08‡ 0.06 -0.76 0.07 -0.20 0.06 -0.46 0.08 category=4 -0.24 0.11 -0.41 0.08 -0.69 0.07 0.02‡ 0.12 -0.22 0.09 -0.13 0.06 -0.42 0.10 -0.15 0.07 -1.01 0.08 -0.74 0.07 -0.68 0.09 category=5 -0.38 0.11 -0.33 0.07 -0.80 0.07 -0.07‡ 0.12 -0.18 0.09 -0.08‡ 0.06 -0.38 0.10 -0.18 0.07 -1.15 0.08 -1.29 0.07 -1.14 0.09  BIC 125 97 420 553 314 226 88 103 337 233 276 Note: the quantity included in the parenthesis represents the total number of observations in a cluster. † variable value in thousands. ‡ coefficient with p value > 0.05. NA = not applicable.TheBIC for these models are provided in Tables 4 and 5. To compare the goodness of fit, overallBIC values were calculated. The overallBIC for the seven-cluster models obtained from the proposed approach was 3,008, whereas theBIC for 11-cluster models, obtained from the existing approach, was 3,171. This difference was the result of similar or better explanatory power provided by the proposed approach with seven clusters versus 11. That is, the more clusters, the more coefficients for explanatory variables needed to be estimated for a similar goodness of fit; thus, theBIC increased.It is observed that the individual BIC for the 11-cluster are slightly smaller than those for the seven-cluster models. Similar to SSE, the more the clusters are used, the smaller the individual model BIC is. However, this does not necessarily translate into a smaller overall BIC.


### 3.4. Model Assessment and Validation

#### 3.4.1. Potential Overfitting

Brusco et al. [53] noted that clusterwise regression models have great potential for overfitting. Often, a variation in the response variable is governed by clustering. Hence, they recommend investigating the potential presence of overfitting in CLR models. This study adopted a procedure proposed by Brusco et al. [53] and the test dataset to diagnose overfitting. For the optimum seven-cluster models, the total sum of squares (TSS) was 2,419, the between-clusters sum of squares (BCSS) was 30, the within-clusters sum of squares (WCSS) was 2,389, the sum of squares due to regression (SSR) was 1,456, and the SSE was 933. The BCSS was around 1% of TSS and SSR was 62% of WCSS. These results indicated that there was no overfitting, as most of the variation in PSI was explained by within-cluster regressions. SSE accounted for 38% of TSS, which suggests that the models still have a relatively high rate of errors. A nonlinear functional form should be investigated to reduce the existing errors.


#### 3.4.2. Model Accuracy

The accuracy of the models obtained from both approaches was assessed by calculating the overall root-mean-square error (RMSE), as follows:(9)RMSE=∑1ηyitk-y^itk2ηwhere, yitk = the observed PSIy^itk = the predicted PSIη = the number of predictionsBoth models were applied to the test dataset. Memberships of the pavement samples were assigned by mapping the sample IDs and memberships determined by the CLR models. Associated regression models and observed data were used to estimate the PSIs. Predicted PSIs then were compared with the observed PSIs, as shown in Figure6. Results indicate that both CLR models overestimated the PSI. A possible reason might be the existence of multicollinearity among explanatory variables.Figure 6
PSI predictions using the test dataset with (a) the proposed CLR approach and (b) the state-of-the-art approach.
(a) (b)The RMSE for 2011 and 2012 were calculated for models obtained using both approaches. The RMSE were, respectively, 0.429 and 0.439 for the 11-cluster models estimated using the existing state-of-the-art and the seven-cluster models estimated using the proposed CLR approach. The prediction accuracy of models estimated by the existing state-of-the-art was slightly higher than that of the model estimated by the proposed approach. However, this difference in accuracy is very small. The additional four clusters required 100 parameters corresponding to twenty-five explanatory variables including the intercept. Hence, seven-cluster models estimated using the proposed approach were more parsimonious and preferred over the 11-cluster models.


## 3.1. Data Resources

Data used in this study were extracted from the PMS database of NDOT [36]. The data consisted of various classes, location data, segment data, contract data, environmental data, traffic data, and pavement condition data, collected for flexible pavement in the entire State of Nevada.A detailed data analysis was performed to check for inconsistent and missing information in the dataset, and some were found. Some of the missing data were synthesized based on associated information available in the dataset. In preparing the PMS data, the following filters were applied:(i)
Only one-mile segments were selected for consistency.(ii)
Only pavement segments with the most recent maintenance contracts awarded in 2001 or later were used in the study.(iii)
PSI of a pavement should deteriorate over time if no M&R treatment occurred. If PSI of a segment in any year increased by 0.1 or more points from the previous year without any M&R treatment, all observations for that year were excluded from the analysis. However, if an increase in PSI in any year was less than 0.1 from the previous year, it was assumed to be a random error during the process of pavement evaluation or data processing. Therefore, those observations were included in the analysis.(iv)
If the PSI of any year decreased by one or more points from the previous year, all observations for that year were excluded from the analysis.(v)
In practice, the PSI range is between 4.5 and 1.5. Therefore, if a pavement segment had a PSI beyond these limits in any year, it was considered an outlier, and all observations for that year were excluded.(vi)
Only PSI values used were within the interval of the mean, minus three standard deviations to the mean plus three standard deviations.(vii)
Pavement samples that did not consist of data regarding conditions for at least two consecutive years were excluded.(viii)
Data analysis showed that an improvement in PSI was seen one or two years after the contract award date. Hence, the age of the pavement sample was set to 0 when the actual improvement occurred rather than when the contract was awarded.After data preparation was completed, 4,138 flexible pavement samples with 17,642 observations were available. For CLR modelling, 14,637 observations, collected from 2001 to 2010, were used; the remaining 3,005 observations, collected in 2011 and 2012 (about 17% of the total number), were used as test dataset to check the accuracy of the CLR models. Tables1 and 2 provide the descriptive statistics of continuous and categorical variables used in this paper, respectively.Table 1
Descriptive statistics for the continuous variables.
Variable Minimum Maximum Mean Std. deviation psi 1.60 4.57 4.01 0.41 age 0.00 8.00 2.24 2.01 adt 20.00 132000.00 4844.45 9812.57 trucks 1.00 7731.00 862.29 1082.20 elevation 228.60 2667.00 1368.25 415.19 precip 3.94 89.28 19.33 10.10 min_temp -6.67 13.33 3.20 4.00 max_temp 7.78 31.67 20.31 4.22 wet_days 11.00 81.00 42.14 15.67 freeze_thaw 0.00 230.00 136.75 51.51 rut_depth 0.00 1.60 0.14 0.14Table 2
Descriptive statistics for the continuous variables.
Variable Category Dummy variable Number of observations Percent System ID IR - 4,622 31.6 NHS nhs 5,063 34.6 STP stp 4,952 33.8 Number of Lanes 1 - 8,450 57.7 2 lane=2 5,612 38.3 ≥ 3 lane≥3 575 3.9 Prioritization Category 1 - 5,004 34.2 2 category =2 3,181 21.7 3 category =3 3,118 21.3 4 category =4 1,558 10.6 5 category =5 1,776 12.1 Functional Class 1 - 4,695 32.1 2 f_class = 2 107 0.7 3 f_class = 3 5,106 34.9 4 f_class = 4 2,604 17.8 5 f_class = 5 1,830 12.5 6 f_class = 6 251 1.7 7 f_class = 7 44 0.3 Note: IR, interstate route; NHS, national highway system; STP, Surface Transportation Program.


## 3.2. Parameters of the Algorithm

Performance of the SA algorithm generally depends on the values of the optimization parameters utilized for a given problem. To ensure proper initialization and search for optimal solutions, selection of the most appropriate parameter values is critical [37, 38]. A body of literature exists regarding various methodologies for finding the most appropriate values for annealing parameters in SA [37, 39–43].If an SA algorithm is allowed to run for a sufficiently long time by setting a high initial temperature with a slow cooling rate, the algorithm performs well, as shown by Anily and Federgruen [44]. In such a cooling scheme, the selection of the most appropriate parameter values may not be critical. However, computation time cannot always be ignored. Hence, the algorithm has to find a good solution in a reasonable amount of time [39].Effective values to be assigned to the optimization parameters depend on the type and complexity of the problem. These values may not be obvious to determine, but rather might be determined by trial and error for a given problem [40]. In this study, values assigned to the optimization parameters were determined using experience gained from previous research [45–51] that involved SA and other comparable algorithms. Table 3 lists the parameter values used in this study.Table 3
Setup parameters for implementation of the proposed algorithm.
Parameter Value Remarks θ 0 10 Initial temperature θ m i n 10e-17 Final minimum temperature B 3000 Boltzmann constant λ 0.97 Cooling rate N m a x 10 Number of neighborhood solutions generated at each temperature level n 800 Minimum number of observations required in a cluster N p s 25 Number of pavement samples, which memberships were changed to generate a neighborhood clusterThe minimum number of observations required in a cluster,n, is one of the parameters to be defined by the analyst for each dataset and application as it is required for any other statistical analysis regarding the minimum sample size. That is, the proposed methodology and contributions are not restricted or affected by the sample size. The analyst must have available sufficient data for each cluster to be able to obtain reliable estimates. A too small n will result in statistically unreliable models and a potentially time-consuming search process as the maximum number of feasible clusters; Kmax will be very large. In contrast, a too large n will result in insufficient number of clusters required to provide the optimum goodness of fit. Hence, sensitivity analysis is required to achieve balance between reliable estimates and Kmax.


## 3.3. Results and Discussion

Given the constraints for feasible partitions defined in the problem formulation and the minimum number of observations required in a cluster,n=800, the proposed algorithm determined 16 as the maximum number of potential clusters. The algorithm searched for the optimum number of clusters from 2 to 16. Seven-cluster CLR models provided the optimum solution with the lowest BIC. The estimated regression coefficients for the CLR models are presented in Table 4.Table 4
Estimated model parameters using the proposed CLR approach.
Parameters C 1  (2,279) C 2  (1,959) C 3  (2,169) C 4  (2,094) C 5  (1,883) C 6  (1,936) C 7  (2,317) Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. Coeff. Std. err. intercept 4.145 0.312 6.242 0.369 2.7910 0.311 6.7280 0.358 7.7810 0.394 12.1400 0.401 3.8730 0.331 age -0.035 0.002 -0.040 0.003 -0.0350 0.003 -0.0327 0.003 -0.0464 0.003 -0.0392 0.004 -0.0498 0.003 adt† -0.006 <0.001 -0.004 <0.001 -0.0262 <0.001 -0.0028 <0.001 -0.0078 <0.001 -0.0053 <0.001 -0.0334 <0.001 t r u c k s † 0.0002‡ <0.001 0.0205‡ <0.001 0.0190‡ <0.001 -0.0151‡ <0.001 0.0306 <0.001 0.0557 <0.001 0.0752 <0.001 e l e v a t i o n † 0.0066‡ <0.001 0.0182‡ <0.001 -0.0352‡ <0.001 -0.1079 <0.001 -0.1418 <0.001 -0.0131‡ <0.001 0.1060 <0.001 precip -0.0037‡ 0.008 -0.0118‡ 0.009 -0.0037‡ 0.008 -0.0248 0.009 0.0094‡ 0.011 -0.0518 0.012 -0.0252 0.008 min_temp -0.030 0.009 0.0129‡ 0.010 -0.0092‡ 0.009 0.0497 0.010 -0.0321 0.011 -0.0447 0.013 0.0532 0.009 max_temp 0.025 0.007 -0.030 0.008 0.0249 0.007 -0.0568 0.008 -0.0124‡ 0.009 -0.0554 0.010 -0.0311 0.007 wet_days 0.005 0.002 -0.010 0.002 0.0115 0.002 0.0031‡ 0.002 -0.0028‡ 0.002 0.0004‡ 0.003 -0.0061 0.002 freeze_ t h a w † -1.697 <0.001 1.513 <0.001 -0.2029‡ <0.001 1.8550 <0.001 -1.4100 0.001 -13.7900 0.001 4.2370 <0.001 rut_depth -0.632 0.115 -1.002 0.115 -1.1060 0.095 -0.5614 0.117 -0.8408 0.119 -0.3999 0.143 -0.9900 0.099 lane=2 -0.371 0.029 -0.166 0.027 0.0121‡ 0.027 -0.1216 0.028 -0.5574 0.032 -0.2745 0.032 0.0258‡ 0.026 lane≥3 -0.325 0.044 -0.195 0.046 -0.1713 0.052 -0.3215 0.045 -0.3974 0.052 -0.2511 0.054 0.0391‡ 0.046 nhs -0.442 0.079 -0.407 0.172 0.7454 0.138 -0.2811 0.072 -0.2639 0.150 -0.4300 0.085 1.4320 0.138 stp -0.802 0.094 -0.487 0.180 0.6121 0.153 -0.3452 0.091 -0.1346‡ 0.131 -0.3379 0.104 1.1350 0.151 f_class=2 0.521 0.115 0.0940‡ 0.165 -0.9383 0.136 0.4118 0.091 0.7050 0.140 0.9474 0.115 -1.5090 0.149 f_class=3 0.382 0.100 0.3377‡ 0.174 -0.8551 0.139 0.3270 0.080 0.3107 0.148 0.4025 0.101 -1.3770 0.136 f_class=4 0.619 0.112 0.2964‡ 0.181 -0.6925 0.151 0.2945 0.094 0.0189‡ 0.137 0.5468 0.117 -1.0330 0.146 f_class=5 0.618 0.113 -0.0703‡ 0.182 -0.6471 0.152 -0.7763 0.098 -0.6250 0.143 0.3015 0.119 -1.5590 0.148 f_class=6 0.472 0.118 -0.419 0.188 -1.2170 0.159 -0.6566 0.108 -0.7955 0.145 0.0582‡ 0.124 -1.7880 0.154 f_class=7 0.554 0.204 0.416 0.199 -1.3840 0.173 -0.1375‡ 0.171 NA 0.151 0.7614 0.144 -1.3640 0.169 category=2 -0.306 0.070 -0.0444‡ 0.051 0.0762‡ 0.056 -0.1300 0.045 -0.6077 0.051 -0.1175‡ 0.065 0.0246‡ 0.035 category=3 -0.319 0.073 -0.0156‡ 0.056 -0.0421‡ 0.060 -0.2457 0.051 -0.6000 0.056 -0.2317 0.069 -0.0080‡ 0.040 category=4 -0.468 0.077 -0.298 0.064 -0.3665 0.065 -0.1535 0.062 -0.6331 0.066 -0.5050 0.077 -0.3982 0.052 category=5 -0.463 0.077 -0.356 0.064 -0.7165 0.065 -0.1899 0.063 -0.8446 0.067 -0.6001 0.076 -0.4145 0.052  BIC 136 338 238 216 496 857 271 Note: the quantity included in parentheses represents the total number of observations in a cluster. † variable value in thousands. ‡ coefficient with p value > 0.05. NA = not applicable.Figure5(a) shows the smallest BIC for each of the clusters (K=2 to 16) considered in this experiment. Figure 5(b) shows the trajectory of the objective function, BIC, when the CLR models were used. The initial value of BIC was 8,502. After 1,360 iterations, the BIC decreased to 3,008. This change was equivalent to an improvement of 65%.Figure 5
BIC trend over the number of clusters (a); trajectory of the BIC during optimization for seven-cluster models (b); SSE trend over the number of clusters (c); and trajectory of the SSE during optimization for eleven-cluster models (d).
(a) (b) (c) (d)It was observed that not all coefficients had associated p values less than 0.05. In this study, the significance level was considered to be 5%. As expected, coefficients differed in magnitude and sign across the clusters, which indicated that the deterioration patterns of pavement samples varied among the clusters. However, seven explanatory variables had the same sign across all clusters.Different clusters had different numbers of significant explanatory variables. For example, Cluster 2 had 10 insignificant explanatory variables. In addition, among all seven clusters, five variables,age,adt,rut_depth,category=4, andcategory=5, were significant. However, four variables,trucks,elevation,precip, andcategory=2, were not significant in four different clusters.Trucks play a key role in pavement deterioration because they transfer heavy loads to the pavement [52]. Hence, it is expected that variable “trucks” be significant. There could be various reasons precluding the statistical significance of this variable. Two of them include (i) multicollinearity effects, which are not addressed in the existing CLR literature and motivates the expansion of the framework as recommend in this paper under future research and (ii) different lanes mistakenly used to collect pavement performance and truck traffic data. The data used in this study does not include the lane used to collect the information.The performance of the proposed CLR approach was compared with that of the existing CLR for pavement management. That is, in this section models estimated using the proposed approach were compared with those estimated using the existing CLR [6] described in the Introduction. Experiments using the existing CLR approach were run for all feasible clusters (K=2 to 16). Figure 5(c) shows the smallest SEE for each of these clusters. As expected, SSE decreased with an increasing number of clusters, but at a very small rate after K=11. In this case, Figure 5(c) does not exhibit a clear elbow point. Hence, an optimum number of clusters needed to be decided by visual inspection while considering the trade-off between goodness of fit and model complexity (i.e., of the number of models and the explanatory variables). This inherent subjectivity when choosing an optimum number of clusters is a major drawback for the existing state-of-the-art CLR approach.After careful assessment, 11-cluster CLR models were selected as the optimum solution. Figure5(d) shows the trajectory of SSE, when 11-cluster CLR models were used, and Table 5 provides the corresponding regression coefficients. Similar to the results obtained from the proposed CLR approach, the coefficients differed in terms of magnitude and sign. In addition, some coefficients had p values larger than 0.05.Table 5
Coefficients obtained using the existing state-of-the-art CLR approach.
Parameters C 1  (1,229) C 2  (1,365) C 3  (1,413) C 4  (1,091) C 5  (1,423) C 6  (1,412) C 7  (1,430) C 8  (1,321) C 9  (1,272) C 10  (1,360) C 11  (1,321) Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. Coeff. Std. Err. intercept 3.92 0.40 4.12 0.40 4.95 0.44 14.43 0.60 7.41 0.45 3.70 0.36 4.61 0.40 3.80 0.47 7.95 0.48 6.50 0.42 7.54 0.43 age -0.04 <0.01 -0.04 <0.01 -0.05 <0.01 -0.04 <0.01 -0.05 <0.01 -0.04 <0.01 -0.04 <0.01 -0.05 <0.01 -0.05 <0.01 -0.04 <0.01 -0.03 <0.01 adt † -0.04 <0.01 -0.02 <0.01 -0.02 <0.01 0.01 <0.01 -0.04 <0.01 0.00 <0.01 -0.01‡ <0.01 -0.04 <0.01 -0.01 <0.01 -0.01‡ <0.01 0.00 <0.01 trucks † 0.11 0.02 0.03 0.02 0.04 0.02 0.13 0.02 0.05 0.02 0.01‡ <0.01 -0.03‡ <0.01 -0.01‡ 0.02 0.05 <0.01 0.01‡ 0.02 0.03 <0.01 elevation † -0.02‡ <0.01 0.18 <0.01 -0.05‡ <0.01 -0.05‡ <0.01 -0.22 <0.01 0.01‡ <0.01 -0.11 <0.01 -0.17 <0.01 0.18 <0.01 -0.05‡ <0.01 0.01‡ <0.01 precip -0.01‡ <0.01 -0.02‡ <0.01 -0.01‡ <0.01 -0.13 <0.01 -0.01‡ <0.01 -0.01‡ <0.01 0.02‡ <0.01 0.04 <0.01 -0.05 <0.01 -0.03 <0.01 -0.05 <0.01 min_temp 0.01‡ <0.01 -0.02‡ <0.01 0.09 <0.01 -0.09 0.02 0.01‡ <0.01 -0.02 <0.01 -0.03 <0.01 -0.02‡ <0.01 0.12 <0.01 -0.01‡ <0.01 0.05 <0.01 max_temp -0.01‡ <0.01 0.01‡ <0.01 -0.06 <0.01 -0.05 <0.01 -0.03 <0.01 0.02 <0.01 0.02‡ <0.01 0.03 <0.01 -0.12 <0.01 -0.01‡ <0.01 -0.06 <0.01 wet_days 0.01‡ <0.01 -0.01 <0.01 -0.01‡ <0.01 0.02 <0.01 0.01‡ <0.01 0.01 <0.01 -0.01‡ <0.01 0.01‡ <0.01 -0.02 <0.01 0.01‡ <0.01 0.01 <0.01 freeze_thaw † 1.00‡ 0.56 -0.86‡ 0.50 6.59 0.55 -20.89 0.85 0.24‡ 0.59 -1.38 0.44 0.78‡ 0.52 1.45 0.59 5.22 0.71 -2.19 0.59 -1.93 0.59 rut_depth -0.64 0.15 -1.00 0.14 -1.22 0.13 -1.03 0.17 -0.27 0.13 -0.46 0.13 -0.76 0.13 -1.02 0.13 -1.80 0.15 -0.88 0.14 -1.20 0.15 lane=2 0.02‡ 0.04 -0.16 0.03 -0.35 0.04 -0.10 0.04 0.03‡ 0.04 -0.06‡ 0.03 -0.26 0.03 -0.22 0.03 -0.53 0.04 -0.28 0.04 -0.26 0.04 lane≥3 0.22 0.06 -0.03‡ 0.06 -0.43 0.06 -0.31 0.06 0.31 0.08 -0.43 0.06 -0.48 0.05 0.01‡ 0.07 -0.33 0.06 -0.27 0.06 -0.39 0.06 nhs 1.39 0.15 0.28‡ 0.17 0.51 0.16 -0.12‡ 0.16 1.20 0.14 -0.25 0.08 -0.41 0.08 0.80 0.18 -0.70 0.15 -0.50 0.17 0.44 0.15 stp 1.33 0.17 0.18‡ 0.19 0.42 0.17 -0.27‡ 0.18 0.67 0.17 -0.38 0.11 -0.50 0.10 0.82 0.19 -0.64 0.15 -0.65 0.19 0.65 0.14 f_class=2 -1.42 0.23 -0.39‡ 0.25 -0.24‡ 0.17 0.44 0.18 -0.83 0.24 0.32 0.09 0.52 0.10 -0.93 0.17 1.43 0.20 -0.25‡ 0.17 -0.71 0.17 f_class=3 -1.36 0.18 -0.55 0.18 -0.38 0.15 0.12‡ 0.18 -0.88 0.15 0.15‡ 0.09 0.33 0.12 -0.99 0.17 0.93 0.14 0.43 0.16 -0.55 0.13 f_class=4 -1.25 0.18 -0.49 0.20 -0.39 0.16 0.22‡ 0.20 -0.56 0.18 0.27 0.11 0.41 0.13 -1.12 0.18 0.89 0.14 0.69 0.18 -0.59 0.12 f_class=5 -1.85 0.19 -0.51 0.20 -0.91 0.16 0.03‡ 0.20 -1.43 0.18 0.21‡ 0.11 0.18‡ 0.13 -2.18 0.18 0.60 0.14 0.43 0.18 -0.47 0.12 f_class=6 -2.07 0.20 -0.71 0.20 -0.77 0.17 -0.72 0.21 -1.71 0.19 0.14‡ 0.12 -0.66 0.14 -2.24 0.19 0.24‡ 0.15 0.56 0.19 -0.71 0.13 f_class=7 NA -0.63 0.22 -1.53 0.25 0.42‡ 0.29 -0.68 0.19 0.09‡ 0.16 0.18‡ 0.16 -2.28 0.25 NA 0.72 0.25 NA category=2 0.02‡ 0.09 0.07‡ 0.06 -0.36 0.05 0.22 0.10 -0.27 0.07 0.06‡ 0.05 -0.20 0.09 -0.07‡ 0.06 -0.66 0.06 -0.13 0.05 -0.14‡ 0.07 category=3 0.02‡ 0.10 0.01‡ 0.07 -0.43 0.06 0.14‡ 0.11 -0.17 0.08 -0.05‡ 0.06 -0.29 0.10 -0.08‡ 0.06 -0.76 0.07 -0.20 0.06 -0.46 0.08 category=4 -0.24 0.11 -0.41 0.08 -0.69 0.07 0.02‡ 0.12 -0.22 0.09 -0.13 0.06 -0.42 0.10 -0.15 0.07 -1.01 0.08 -0.74 0.07 -0.68 0.09 category=5 -0.38 0.11 -0.33 0.07 -0.80 0.07 -0.07‡ 0.12 -0.18 0.09 -0.08‡ 0.06 -0.38 0.10 -0.18 0.07 -1.15 0.08 -1.29 0.07 -1.14 0.09  BIC 125 97 420 553 314 226 88 103 337 233 276 Note: the quantity included in the parenthesis represents the total number of observations in a cluster. † variable value in thousands. ‡ coefficient with p value > 0.05. NA = not applicable.TheBIC for these models are provided in Tables 4 and 5. To compare the goodness of fit, overallBIC values were calculated. The overallBIC for the seven-cluster models obtained from the proposed approach was 3,008, whereas theBIC for 11-cluster models, obtained from the existing approach, was 3,171. This difference was the result of similar or better explanatory power provided by the proposed approach with seven clusters versus 11. That is, the more clusters, the more coefficients for explanatory variables needed to be estimated for a similar goodness of fit; thus, theBIC increased.It is observed that the individual BIC for the 11-cluster are slightly smaller than those for the seven-cluster models. Similar to SSE, the more the clusters are used, the smaller the individual model BIC is. However, this does not necessarily translate into a smaller overall BIC.


## 3.4. Model Assessment and Validation

### 3.4.1. Potential Overfitting

Brusco et al. [53] noted that clusterwise regression models have great potential for overfitting. Often, a variation in the response variable is governed by clustering. Hence, they recommend investigating the potential presence of overfitting in CLR models. This study adopted a procedure proposed by Brusco et al. [53] and the test dataset to diagnose overfitting. For the optimum seven-cluster models, the total sum of squares (TSS) was 2,419, the between-clusters sum of squares (BCSS) was 30, the within-clusters sum of squares (WCSS) was 2,389, the sum of squares due to regression (SSR) was 1,456, and the SSE was 933. The BCSS was around 1% of TSS and SSR was 62% of WCSS. These results indicated that there was no overfitting, as most of the variation in PSI was explained by within-cluster regressions. SSE accounted for 38% of TSS, which suggests that the models still have a relatively high rate of errors. A nonlinear functional form should be investigated to reduce the existing errors.


### 3.4.2. Model Accuracy

The accuracy of the models obtained from both approaches was assessed by calculating the overall root-mean-square error (RMSE), as follows:(9)RMSE=∑1ηyitk-y^itk2ηwhere, yitk = the observed PSIy^itk = the predicted PSIη = the number of predictionsBoth models were applied to the test dataset. Memberships of the pavement samples were assigned by mapping the sample IDs and memberships determined by the CLR models. Associated regression models and observed data were used to estimate the PSIs. Predicted PSIs then were compared with the observed PSIs, as shown in Figure6. Results indicate that both CLR models overestimated the PSI. A possible reason might be the existence of multicollinearity among explanatory variables.Figure 6
PSI predictions using the test dataset with (a) the proposed CLR approach and (b) the state-of-the-art approach.
(a) (b)The RMSE for 2011 and 2012 were calculated for models obtained using both approaches. The RMSE were, respectively, 0.429 and 0.439 for the 11-cluster models estimated using the existing state-of-the-art and the seven-cluster models estimated using the proposed CLR approach. The prediction accuracy of models estimated by the existing state-of-the-art was slightly higher than that of the model estimated by the proposed approach. However, this difference in accuracy is very small. The additional four clusters required 100 parameters corresponding to twenty-five explanatory variables including the intercept. Hence, seven-cluster models estimated using the proposed approach were more parsimonious and preferred over the 11-cluster models.


## 3.4.1. Potential Overfitting

Brusco et al. [53] noted that clusterwise regression models have great potential for overfitting. Often, a variation in the response variable is governed by clustering. Hence, they recommend investigating the potential presence of overfitting in CLR models. This study adopted a procedure proposed by Brusco et al. [53] and the test dataset to diagnose overfitting. For the optimum seven-cluster models, the total sum of squares (TSS) was 2,419, the between-clusters sum of squares (BCSS) was 30, the within-clusters sum of squares (WCSS) was 2,389, the sum of squares due to regression (SSR) was 1,456, and the SSE was 933. The BCSS was around 1% of TSS and SSR was 62% of WCSS. These results indicated that there was no overfitting, as most of the variation in PSI was explained by within-cluster regressions. SSE accounted for 38% of TSS, which suggests that the models still have a relatively high rate of errors. A nonlinear functional form should be investigated to reduce the existing errors.


## 3.4.2. Model Accuracy

The accuracy of the models obtained from both approaches was assessed by calculating the overall root-mean-square error (RMSE), as follows:(9)RMSE=∑1ηyitk-y^itk2ηwhere, yitk = the observed PSIy^itk = the predicted PSIη = the number of predictionsBoth models were applied to the test dataset. Memberships of the pavement samples were assigned by mapping the sample IDs and memberships determined by the CLR models. Associated regression models and observed data were used to estimate the PSIs. Predicted PSIs then were compared with the observed PSIs, as shown in Figure6. Results indicate that both CLR models overestimated the PSI. A possible reason might be the existence of multicollinearity among explanatory variables.Figure 6
PSI predictions using the test dataset with (a) the proposed CLR approach and (b) the state-of-the-art approach.
(a) (b)The RMSE for 2011 and 2012 were calculated for models obtained using both approaches. The RMSE were, respectively, 0.429 and 0.439 for the 11-cluster models estimated using the existing state-of-the-art and the seven-cluster models estimated using the proposed CLR approach. The prediction accuracy of models estimated by the existing state-of-the-art was slightly higher than that of the model estimated by the proposed approach. However, this difference in accuracy is very small. The additional four clusters required 100 parameters corresponding to twenty-five explanatory variables including the intercept. Hence, seven-cluster models estimated using the proposed approach were more parsimonious and preferred over the 11-cluster models.


## 4. Conclusions

This study proposed and implemented a clusterwise multiple linear regression to develop pavement performance models. A mixed-integer nonlinear mathematical program was formulated to explain the problem. The CLR approach simultaneously divided pavement samples into an optimum number of clusters, and estimated a PPM for each cluster.In the experiments, various environmental factors were considered as potential explanatory variables, including elevation, annual precipitation, average minimum and maximum temperatures, the number of wet days, and freeze and thaw cycles. The proposed approach enabled consideration of other types of variables, such as economic and social factors. Formulation of the mathematical program developed in this study supports a number of explanatory variables, multiple observations per pavement segment, and user-defined constraints on cluster characteristics.A simulated annealing coupled with OLS was used to solve the mathematical program. For the data used in the experiments, the algorithm found that 7-cluster models provided the optimum solution. Results obtained from the proposed CLR models were compared with results obtained from the state-of-the-art approach. This comparison showed that the proposed CLR approach performed better than the state-of-the-art approach in predicting the PSI of pavement samples.The analysis showed that overfitting was not an issue for the resulting clusters and regression models. As expected, the use of the BIC as an objective function to determine the best model specification provided a more parsimonious structure compared with that obtained using SSE. This was a consequence of the consistency property of the BIC [10, 20–23].


## 5. Future Work

This study did not address all limitations of the state-of-the-art CLR approach that were discussed in the Introduction of this paper. The following are potential extensions.One limitation in this study is related to the error that likely is caused by including insignificant explanatory variables during clustering and regression analyses. The proposed formulation needs to be extended to include only significant variables. Hence, the CLR models would not be restricted only to prespecified explanatory variables. Instead, the models could include cluster-specific significant explanatory variables. In addition, the multicollinearity among explanatory variables should be investigated to exclude highly correlated variables during the model estimation process.The proposed mathematical formulation was limited to a linear functional form for PPMs. Luo and Chau [1] implemented a CLR approach using an exponential form. However, their study used only pavement age as an explanatory variable. An interesting aspect worthy of investigation would be to explore utilizing the proposed CLR approach by allowing nonlinear relationships between pavement performance measures and multiple explanatory variables.The proposed SA with the OLS algorithm was designed to search for a global minimum. However, a large amount of computational time was required. Hence, another avenue for future research would be to develop faster and more efficient combinatorial algorithms that could guarantee global optimality.



---

*Source: 2159865-2018-12-12.xml*