A novel weighted multi-output neural network (NN) model is proposed for 3 predicting the deterioration of rigid pavements based on Iowa pavement management 4 system data. This first-of-a-kind model simultaneously predicts four pavement 5 condition metrics concerning rigid pavements, including IRI, faulting, longitudinal 6 crack and transverse crack. It provides an opportunity to efficiently evaluate pavement 7 conditions and to make treatment decisions based on multi-condition metrics, such as 8 the pavement condition index (PCI) for budget allocation models. Compared to 9 traditional single-output NN models, this multi-output model is capable of 10 incorporating correlations among different condition metrics. During model training, 11 each condition metric is assigned a weight to reflect its relative importance. When the 12 weights equal to those in the formula for the multi-condition metric, the prediction 13 performance for PCI is optimal (13% lower MSE than optimal, single-output models). 14 The multi-output model improves on the prediction performance for three of the four 15 individual condition metrics compared to optimal single-output models. Results show 16 that the consideration of correlations could improve the prediction performance for 17 single and multi-condition metrics. Finally, variable weighting is critical for achieving 18 the optimal balance of prediction performance among the various metrics as dictated by 19 the needs of the decisionmaker. 20 Keywords: multi-output; neural network; deterioration; rigid pavement; correlation 21 1. Introduction 22 Performance-based planning (PBP) is an important tool to mitigate the pervasive 23 problem of inadequate budget faced by transportation agencies (Baker et al. 2006). A key 24 element for implementing PBP is efficient prediction of future pavement conditions. This 25 depends on a robust deterioration prediction model. Over the last decades, many deterioration 26 models have been described in the literature. Given the complex physical and chemical 27 processes involved in pavement deterioration, most existing deterioration models are data-28 driven using many statistical methods including Markov chains (Kobayashi et al. 2012, 29 Lethanh and Adey 2013), linear/nonlinear regression (Prozzi and Madanat 2004, Mahmood et 30 al. 2018), neural network (Choi et al. 2004, Hossain et al. 2019), and decision trees (Inkoom 31 et al. 2019). 32 In practice, PBP-based decisions are based on multiple metrics of pavement 33 condition. These can include metrics of pavement roughness (e.g., the international roughness 34 index (IRI)), rutting, faulting, and various forms of cracking (American Association of State 35 Highway and Transportation Officials 2008, Bektas, Fatih; Smadi, Omar G.; and Al-Zoubi 36 2014, FHWA 2019). To date, deterioration models described in the literature predict single 37 metrics (there are papers that describe multiple models for multiple metrics, but each metric 38 is projected using a distinct single model). As will be shown in this paper, by examining 39 panel data from the Iowa Pavement Management System (PMS) (which will be discussed in 40 Section 3), the authors have observed that several of these metrics evolve in a coordinated 41 manner. Failure to account for this correlation may lead to misestimation of future pavement 42 condition, particularly in a stochastic context. 43 In making PBP-based decisions, transportation agencies may also use an overall 1 condition metric – that is a weighted composite metric – to evaluate pavement performance. 2 For example, the state of Iowa uses the pavement condition index (PCI), which is a weighted 3 function of IRI, faulting, longitudinal and transverse cracking for rigid pavements (Bektas, 4 Fatih; Smadi, Omar G.; and Al-Zoubi 2014). As will be shown in this paper, the nature of this 5 weighting should be considered in the development of a correlation-aware prediction model. 6 In this paper, a new weighted multi-output neural network model for rigid pavements 7 is proposed and evaluated using empirical data from the Iowa Pavement Management 8 System. Considering the importance of data cleaning for the training of a data-driven model, 9 a detailed description for the data cleaning process is presented. This model simultaneously 10 predicts four individual metrics of pavement condition and a composite PCI. Different from 11 multiple single-output models, this multi-output model can and does incorporate correlations 12 among different condition metrics. Results show that after incorporating correlations among 13 different condition metrics based on the weighed multi-output model, the prediction 14 performance of the overall condition metric has been improved. In addition, the prediction 15 performances of individual condition metrics are equal to or better than the performance of 16 single-output models. Finally, by adding weights for each condition metric during the model 17 training process, their importance levels can be tailored to the specific decision-making 18 context. 19 This paper is organized as follows. Firstly, a comprehensive literature review is 20 provided. Secondly, in the Data Preparation section, the processes of data collection and 21 cleaning are described in detail. Next, the weighted multi-output neural network is 22 introduced, including the selection of model input parameters, the normalization of training 23 data, the model framework for rigid pavements and the model training process. Then the 24 results about the comparison among single-output models and multi-output models with 25 different weights are discussed. Finally, the main conclusions are summarized, and 26 recommended future research activities are presented. 27 2. Literature Review 28 For data-driven models, the first methodological consideration is selection of the 29 model form and the method to estimate model parameters in the deterioration model. As 30 suggested in (Swei et al. 2018), deterioration modeling methods can be divided into two 31 categories, Markov chains and regression models. Models based on Markov chains focus on 32 the probability that pavement condition evolves from one state to another state (Hong and 33 Wang 2003, Kobayashi et al. 2012, Lethanh and Adey 2013, Abaza 2016). This type of 34 models usually assumes the variation of pavement deterioration to be aleatory. Since the 35 Markov property assumes that future deterioration rates only depend on the current pavement 36 condition (i.e. pavement state), historical dependence is ignored during the prediction 37 process. 38 There are several types of regression models, among which the most commonly 39 applied in the deterioration literature is linear regression (including both explicitly linear 40 forms (Mahmood et al. 2018) and linear forms by transformation1 (Prozzi and Madanat 2003, 1 2004, Swei et al. 2018)). These models can be readily solved by ordinary least squares 2 (OLS). In order to consider the heterogeneity of individual pavements in the same pavement 3 group, Zhang et al. proposed a clusterwise linear regression model (Zhang and Durango-4 Cohen 2013), while Yu et al. proposed linear mixed effects models (Yu et al. 2007). When 5 the prediction is concerned with classification or the probability of performance failure 6 occurrence, logistic regression can be applied. Heidari et al. applied a logistic regression to 7 predict the subclasses of pothole, rutting and protrusion for forest roads (Heidari et al. 2018). 8 Chen et al. evaluated the failure probability of asphalt preventive treatments by both logistic 9 and Bayesian logistic models (Chen et al. 2019). 10 The second type of regression models are non-linear models which are described by 11 the sum of several nonlinear variables or by the exponential form. AASHTO models are 12 examples of this type that were once widely applied and cited. These models were estimated 13 from the AASHO road test data and updated for several versions (Highway Research Board 14 1962, American Association of State Highway and Transportation Officials 1993). Based on 15 the first version of the AASHTO model, Ben-Akiva et al. developed a latent model for 16 pavement deterioration, which assumed that the pavement performance condition is a latent 17 variable (Ben-Akiva and Ramaswamy 1993). Abaza et al. proposed a deterministic prediction 18 model based on the AASHTO serviceability concept and an incremental solution of the 19 AASHTO basic design equation (Abaza 2004). Following the AASHTO model form, Hong 20 et al. developed a nonlinear model that incorporated heterogeneity using Bayesian analysis 21 (Hong and Prozzi 2006, 2010). The Mechanistic-Empirical Pavement Design Guide 22 (MEPDG) is a collection of models that were developed as an update to AASHTO through 23 the National Cooperative Highway Research Program (NCHRP) Project 1-37A (American 24 Association of State Highway and Transportation Officials 2008). Considering the fact that 25 the input parameters in the MEPDG model could be estimated by separate models and thus 26 cause bias, Aguiar-Moya et al. introduced instrumental variables to update the MEPDG 27 model (Aguiar-Moya et al. 2011). Other nonlinear models include the exponential form 28 (Lamptey et al. 2008, Rosa et al. 2017) and multivariate adaptive regression splines (MARS) 29 proposed by (Attoh-Okine et al. 2003). 30 Another class of regression models that have been applied to predict pavement 31 performance are neural networks. A neural network model is an example of machine learning 32 where a prediction algorithm is developed through iterations among a training set of values. 33 The ultimate prediction algorithm combines both linear and non-linear elements weighted to 34 produce a best prediction. Recent years have witnessed rapid expansion in the application of 35 neural network models due to their efficiency to describe nonlinear relationships in many 36 phenomena including pavements’ deterioration (Nii O. Attoh-Okine 1994, Owusu-Ababio 37 1998, Roberts and Attoh-Okine 1998, Lou et al. 2001, Lee and Lee 2004, Thube 2012, 38 Mazari and Rodriguez 2016, Gong et al. 2019, Inkoom et al. 2019). Neural network models 39 usually consist of input layer, hidden layer and output layer. With the increase in model 40 complexity, such as the increase of hidden layers and the increase of the neurons in each 41 hidden layer, the neural network could fit the training data perfectly but may lack the 42 1 Models like can be transformed into explicitly linear models by log-transformation: fx()=∏iaixibilogfx()=∑ilogai+bi⋅logxi() capacity to predict the future deterioration at the same time. This phenomenon is called 1 ‘overfitting’. There are several ways to detect and prevent it, such as out-of-sample 2 prediction, cross-validation, regularization, and early-stopping, etc. However, some exiting 3 models may lack the analysis to check if the neural network model is overfitting or not (Nii 4 O. Attoh-Okine 1994, Owusu-Ababio 1998, Roberts and Attoh-Okine 1998, Mazari and 5 Rodriguez 2016). These models focus more on the fitting process. Attoh-Okine (1994) and 6 Roberts and Attoh-Okine (1998) mainly use the training error as the evaluation criteria for 7 model performance, Owusu-Ababio (1998) and Mazari and Rodriguez (2016) just use 1-fold 8 validation instead of multi folds (i.e. cross validation). 9 There are also some other types of regression models. Inkoom et al. proposed a 10 classification and regression tree (CART) model to predict crack condition (Inkoom et al. 11 2019). Pan et al. used fuzzy regression to predict PSI by considering visual inspection data as 12 fuzzy data (Pan et al. 2011). Bianchini et al. proposed a neuro-fuzzy model with the 13 consideration of IF-THEN fuzzy rules (Bianchini and Bandini 2010). 14 Existing models mainly focus on the prediction of a single condition metric – a single 15 overall metric like PCI, or a specific metric like IRI. This type of single-output models lacks 16 the consideration of correlations among different condition metrics. They could also 17 potentially increase the computational cost when generating multi specific metrics to obtain 18 an overall metric, and lead to some inconvenience for a pavement management system whose 19 treatment decisions are based on multi metrics. 20 The second methodological consideration is the selection of predictor (input) 21 variables for deterioration models, which is strongly related with the available parameters in 22 the PMS dataset available to the model developer. Common input variables include traffic 23 volumes/loads, pavement structure (pavement type, thickness), pavement age, pavement 24 conditions, and maintenance history. Environmental factors can also have a significant 25 influence on the pavement deterioration, such as ambient temperature, precipitation, freeze-26 thaw cycles, and freeze index. Incorporating all parameters in regression models, on one hand 27 could increase the model complexity and lead to overfitting, on the other hand it prevents us 28 from understanding the role of each parameter. Several studies have conducted sensitivity 29 analysis to explore the significance of input variables (Choi et al. 2004, Kargah-Ostadi et al. 30 2010, Karlaftis and Badr 2015, Hossain et al. 2019, Yao et al. 2019). 31 Models are inherently limited to the variables described in available data. 32 Unfortunately, several current PMS datasets do not fully describe the maintenance history of 33 a pavement segment, limiting the ability to model this effect. For example, many models use 34 total pavement thickness as an input variable. However, two pavement segments with the 35 same total thickness are very probable to have different deterioration rates: one asphalt 36 segment has an original 12-inch (304.8 mm) thickness and does not have an overlay, the 37 other one has an original 8-inch (203.2 mm) thickness and also has a 4-inch (101.6 mm) 38 asphalt overlay. 39 Another important methodological consideration is the data “cleaning” process, which 40 usually takes data scientists about 60% of the total time for a project (Press 2016). One 41 common problem is the identification and processing of outliers in the dataset. For example, 42 when there are no treatment actions, the pavement condition should get worse due to the 43 traffic and environmental influence. However, in much real-world PMS data, like LTPP, 44 overall pavement condition is sometimes recorded as getting better from one measurement to 45 the next. There are several reasons for the existence of outliers, such as measurement error or 46 missing treatment records. Given the same raw dataset, with different cleaning criteria, 1 different training datasets are obtained. When datasets differ, even the same training method 2 can generate different models with different evaluation performance. This makes it 3 challenging to compare different methods when cleaning methods are not reproducible. Most 4 existing papers in the deterioration modeling literature have an obscure description about the 5 data cleaning process. 6 To summarize, there are several gaps for existing pavement deterioration models 7 described in the literature. The first gap is about the correlations among condition metrics. 8 Current models are single-output ones that focus on either an overall or a specific condition 9 metric. They lack the consideration of correlations among different condition metrics. The 10 second gap is about input variable selection. Variables in most existing models may not 11 reflect the maintenance history of a pavement segment. The last gap is the unclear description 12 about the data cleaning process, leading to the difficulty of comparing different models. 13 To bridge current gaps in pavement deterioration models, an innovative weighted 14 multi-output neural network model is proposed. Most existing models only focus on flexible 15 pavements. According to the FHWA road statistics (FHWA 2017), rigid pavements account 16 for over 25% in terms of interstate system. This paper uses rigid pavements as an example to 17 illustrate this new multi-output neural network model. A detailed description of data cleaning 18 process is presented. This new model incorporates correlations among different condition 19 metrics during the model training process. Considering the fact that different metrics may 20 have different importance levels for transportation agencies, each output condition metric is 21 assigned a weight number instead of being uniformly weighted during the training process. 22 Model evaluation results show that after incorporating correlations among different condition 23 metrics, the prediction of the overall condition metric could be improved. In addition, the 24 prediction performance of a single condition metric is better or equal to the performance of 25 single-output models. 26 3. Data Preparation 27 The proposed deterioration model is trained on data from the Iowa PMS database and 28 climate data from LTPP. In the following section, a detailed description is presented for the 29 data preparation, which is concerned with data extraction, transformation and cleaning in 30 terms of input and output variables. At the end of this section, the selection of input variables 31 is discussed based on a correlation analysis. 32 3.1 Iowa PMS data 33 The Iowa PMS database contains detailed records for around 4,000 pavements 34 segments on 26 years of pavement condition for the years 1992 to 2017. It covers three 35 systems: 1,566 miles (2,520 km) for interstate system, 4,626 miles (7,445 km) for U.S route 36 system, and 4,891 miles (7,871 km) for Iowa route system. Approximately, 15.6% of total 37 miles are asphalt pavements, 32.9% are concrete pavements, and 51.5% are composite 38 pavements. 39 This data can be divided into three groups based on the units used and the presence of 40 crack information as shown in Table 1. Because the early period of data (Group 1) contains 41 different information about the state of cracking, it is not considered in this analysis. In the 42 dataset, there are two variables called ‘IRIDAT’ and ‘CAPDAT’ which record the dates of 43 measurement for the IRI and crack information, respectively. For each segment in the 44 network, the performance data is collected biennially. For the dataset of year 2013, 1 ‘CAPDAT’ is missing, and for the dataset of year 2016 and 2017, ‘IRIDAT’ is missing. 2 Therefore, the PMS data that is used for training the deterioration model is based on year 3 2000-2012, year 2014 and 2015. 4 5 Table 1. Iowa PMS data groups 6 Group Period Units Crack info 1 1992-1999 1992-1995: English, 1996-1999: Metric DCRACK, HCRACK, LCRACK, TCRACK 2 2000-2012 Metric ACRACK, LCRACK, LWCRACK, TCRACK 3 2013-2017 English ACRACK, LCRACK, LWCRACK, TCRACK * ACRACK: alligator crack, DCRACK: durability crack, HCRACK: half-crack, LCRACK: longitudinal 7 crack, LWCRACK: longitudinal wheel-path crack, TCRACK: transverse crack 8 9 The detailed information for parameters in the selected PMS data can be found in 10 (Iowa DOT Open Data 2019). Based on suggestions from engineering experts, four types of 11 parameters have been selected in this analysis as listed in Table 2. Parameters in the Basic 12 Info group include the segment name (ORIGKEY), the year for the database (PMISYR), 13 system type (SYSTEM, including interstate, US route, and Iowa route), the construction year 14 (CONSYR) and the year for last resurfacing (RESYR). Group Traffic includes two types of 15 traffic statistics, AADT and AADTT. Different kinds of condition metrics and corresponding 16 measurement dates are in Group Distress. The last group Maintenance History records 17 treatment information for each layer. 18 19 20 Table 2. Description of selected parameters. 21 Data type Name Description Basic Info ORIGKEY Original smart key PMISYR Pavement management year SYSTEM Roadway system PAVTYP Pavement type CONYR Year of construction or reconstruction RESYR Year of last resurfacing Traffic AADT Average annual daily traffic TRUCKS Average annual daily truck traffic Distress IRI International roughness index RUT Rut depth FAULT Average faulting only on faulted joints in a segment ACARCKH, ACRACKM, ACRACKL High/moderate/low severity alligator cracks LCARCKH, LCRACKM, LCRACKL High/moderate/low severity longitudinal cracks LWCARCKH, LWCRACKM, LWCRACKL High/moderate/low severity longitudinal wheel-path cracks TCARCKH, TCRACKM, TCRACKL High/moderate/low severity transverse cracks IRIDAT IRI test year CAPDAT Crack & patch collected by vender test year Maintenance History TREATMENT Surface treatment type LAYR (1-8) Layer year # (1-8) SURTYP (1-8) Surface type # (1-8) SURTHK (1-8) Surface thickness # (1-8) BASTYP (1-8) Base type # (1-8) BASTHK (1-8) Base thickness # (1-8) SUBTYP (1-8) Subbase type # (1-8) SUBTHK (1-8) Subbase thickness # (1-8) RMVTYP (1-8) Removal type # (1-8) RMVTHK (1-8) Removal thickness # (1-8) 1 As for Basic Info parameters, based on ‘PMISYR’, ‘CONYR’ and ‘RESYR’, two 2 new variables ‘AGECON’ and ‘AGERES’ were generated by equation (1) and (2). These two 3 variables were used to describe the period since construction/reconstruction and last 4 resurfacing, respectively. 5 (1) 6 (2) 7 Distress parameters in year 2014 and 2015 were transformed into metric units. Each 8 type of crack distress is described by three sub-categories, which can be converted to a single 9 crack distress metric by equation (3) (Bektas, Fatih; Smadi, Omar G.; and Al-Zoubi 2014). 10 Using this, four new variables can be obtained, namely: ‘ACRACK’, ‘LCRACK’, 11 ‘LWCRACK’ and ‘TCRACK’. 12 (3) 13 According to the Iowa Department of Transportation, overall pavement condition is 14 not evaluated directly from condition metric (e.g., IRI). Instead, condition metrics are 15 transformed into indexes including the cracking index, riding index (IRI), and faulting index 16 (Bektas, Fatih; Smadi, Omar G.; and Al-Zoubi 2014). Table 3 lists the threshold values for 17 the calculation of indexes for rigid pavements. Each metric has two threshold values. When 18 the metric is less than or equal to the small threshold value , the index equals 100; 19 When the metric is larger than or equal to the large threshold value , then the index equals 20 0. Other index values can be obtained by linear interpolation. Following this, for any metric 21 , the corresponding index, , is defined formally as: 22 (4) 23 An integrated cracking index was obtained as a weighted sum of indexes of LCRACK and 24 TCRACK with the corresponding weights listed in Table 4. 25 26 Table 3. Thresholds for condition metrics 27 IRI (m/km) FAULT (mm) LCRACK (m/km) TCRACK (count/km) Threshold 0.5&4 0&12 0&250 0&150 AGECON=PMIYR−CONYRAGERES=PMIYR−RESYRCRACK=CRACKL+1.5CRACKM+2CRACKHzζ*ζ*zIzIζ=100ζ≤ζ*ζ−ζ*ζ*−ζ*⋅100ζ*<ζ<ζ*0ζ*≤ζ⎧⎨⎪⎪⎩⎪⎪ Table 4. Cracking sub-index weights 1 sub-index LCRACK TCRACK weights 0.4 0.6 2 Based on these indexes, an overall condition metric called PCI was obtained by 3 equation (5) for rigid pavements. PCI ranges from 0 to 100, and 100 represents that pavement 4 is in perfect condition. In this expression, each of the condition indexes have different 5 coefficients which may represent the significance of that index in the eyes of the Iowa 6 transportation agency. 7 (5) 8 The last group of parameters describe the maintenance history for each segment, from 9 which historical structure information could be generated, including construction year, 10 resurface year, thickness and pavement type. Table 5 presents an example for the 11 maintenance history for a rigid pavement segment. Its ORIGKEY is “08012183 67192 12 8279”. This segment was constructed in 1964 and maintained with an overlay of Portland 13 cement concrete (PCC) in 1984. Based on Table 5, structure information could be obtained as 14 shown in Table 6. From year 1964 to year 1983, the segment structure stayed the same, 15 AGECON (equation (1)) changes each year and AGERES remained 0. After year 1984, due 16 to the overlay action, the segment structure was changed. The new surface thickness 17 SURTHK is 102mm, and previous surface PCC thickness become PCCTHK. TOTTHK is the 18 total thickness, which is equal to the sum of SURTHK and PCCTHK (RMVTHK should be 19 considered when it exists). MRRNUM is the number of treatment actions. In this case, 20 MRRNUM became 1 since 1984. Instead of a single total thickness value as shown in the 21 original PMS dataset, the introduction of SURTHK, PCCTHCK and MRRNUM could reflect 22 the maintenance history of a pavement segment. (If the pavement type is asphalt or composite 23 pavements, a variable called ‘HMATHK’ is introduced similar to PCCTHK.) 24 25 Table 5. Maintenance history for a rigid pavement segment 26 LAYR SURTYP SURTHK BASTYP BASTHK 1964 PCC 254 GSB 102 1984 PCC 102 - - *GSB: granular subbase 27 28 Table 6. Structure information for a rigid pavement segment 29 PMISYR CONYR RESYR PAVTYP TOTTHK SURTHK PCCTHK MRRNUM 1964~1983 1964 - PCC 254 254 0 0 >=1984 1964 1984 PCC 356 102 254 1 30 After generating parameters for model training, the total dataset is filtered based on 31 ‘PAVTYP’ to obtain the data for rigid pavements. Next, input and output variables are 32 determined with the consideration of time series analysis and treatment actions. In terms of 33 the time series analysis, time lag = 1 is suggested in (Chu and Durango-Cohen 2008, Swei et 34 al. 2018), which means that the condition prediction at year t+1 is only related with the 35 information at year t. Since the condition metrics are collected every two years in Iowa, time 36 lag is selected as 2 years here. For example, the information at year 2000 can be considered 37 as inputs and the condition metrics at year 2002 are outputs. 38 0.4 0.4 0.2Crack Ride FaultPCI I I I=+ + During the process of paring input and output variables, attention should be paid to 1 pavement treatments. If there is a treatment action at year 2000, then it is not appropriate to 2 use year 2000 as input and 2002 as output. This kind of pairing data should be excluded from 3 the training dataset. In the dataset, the improvement of pavement condition may be postponed 4 after a treatment. For example, there is a treatment at year 2002, but the condition 5 improvement is reflected at year 2004. Therefore, during the pairing process, when there 6 exists a treatment action, its adjacent years are checked to decide which year should be 7 considered as the treatment year. 8 After this step, input variables include: 1). Basic info: AGECON, AGERES; 2). 9 Traffic: ADT and TRUCKS; 3). Distress: , , , ; 4). 10 Maintenance history: TOTTHK, SURTHK, PCCTHK, MRRNUM. Output variables include 11 , , , . The subscription t and t+2 represent 12 measurement years of condition metrics. 13 14 3.2 Iowa PMS data cleaning 15 After generating the initial input and output variables, the next step is to deal with 16 data outliers. For this case, outliers are defined as values or pairs of values that indicate 17 physically unreasonable pavement condition values. In terms of the overall pavement 18 condition (i.e. PCI in this analysis), it should decrease when there is no treatment action due 19 to running traffic and environmental effects (American Association of State Highway and 20 Transportation Officials 2008). But for other single performance metrics (IRI, faulting and 21 cracks), their performances may improve due to some negative correlations. Differences for 22 all performance metrics between input and output are calculated, 23 (6) 24 (7) 25 (8) 26 (9) 27 (10) 28 When is larger than 0, the data point is considered as an outlier and is deleted from the 29 training dataset. In terms of other single performance metrics, interquartile range (IRQ) is 30 applied, which can be calculated by equation (11), 31 (11) 32 Where and represent values of 25th and 75th percentiles, respectively. Usually, if a 33 data point is below or above , it is considered as an outlier (Rice 34 2006). In this case, all lower quartile boundaries are negative values, indicating unexplained 35 performance improvement. These lower bounds have been used to identify outliers. There is 36 less theoretical support to suggest that pavements could not deteriorate rapidly. Therefore, an 37 upper bound of one half of the threshold values has been adopted for use in Table 3, namely 2 38 for IRI. Table 7 shows the outlier threshold for condition metrics based on IQR and the half 39 IRItFAULTtLCRACKtTCRACKtIRIt+2FAULTt+2LCRACKt+2TCRACKt+2 ΔPCI=PCIt+2−PCItIRI_diff=IRIt+2−IRItFAULT_diff=FAULTt+2−FAULTtLCRACK_diff=LCRACKt+2−LCRACKtTCRACK_diff=TCRACKt+2−TCRACKtΔPCIIQR=Q3−Q1Q1Q3Q1−1.5⋅IQRQ3+1.5⋅IQR of index approach. Values outside the lower and upper bound were eliminated from the 1 training dataset. 2 3 Table 7. Outlier threshold for condition metrics 4 -0.3 -4.9 -23.0 -13.5 0.4 6.0 34.4 18.5 Half of index threshold 2.0 6.0 125.0 75.0 5 3.3 Climate data 6 Climate data was primarily extracted from LTPP, including annual average 7 temperature (TEMP), precipitation (PRECIP), and freeze index (FREEZE). For a given year, 8 there are several records for each climate variable. Annual averages of these were used. 9 Climate data was matched to the other training data (which are generated in section 3.1 and 10 3.2) by condition metric measurement year, IRIDAT and CAPDAT. 11 4. Weighted multi-output neural network model 12 4.1 Selection of input parameters 13 Based on the previous analysis, a candidate set of input variables were generated. 14 However, some variables may not contribute to the pavement deterioration process to an 15 extent that is observable in the current data. To select the most relevant input parameters, a 16 correlation analysis was conducted. This approach is similar to the linear regression analysis 17 described by (Kargah-Ostadi et al. 2010). 18 Fig 1 shows the heatmap for the correlations among input and output variables. The 19 output variables are the variations of condition metrics, including IRI_diff, FAULT_diff, 20 LCRACK_diff, and TCRACK_diff (equations (7)-(10)). All other listed variables are 21 candidate inputs. In this analysis, an input variable would be excluded if either the absolute 22 value of the correlation coefficient is less than 0.05 for all four outputs or if it is less than 23 0.05 for three outputs and less than 0.1 for the fourth. Using these criteria for this dataset, 24 AGERES, PCCTHK, MRRNUM and TEMP were excluded. One main reason that output 25 response is not correlated with AGERES, PCCTHK, MRRNUM is that for most rigid 26 pavements in Iowa the historical overlay number is zero across the dataset. When overlay 27 number is zero, AGERES, PCCTHK, MRRNUM are all equal to zero, and TOTTHK is equal 28 to TOPTHK for most rigid pavements. To simply the problem, AGERES, PCCTHK, and 29 MRRNUM are deleted from the set of input variables. But to be noted that, for asphalt and 30 composite pavements, MRRNUM is usually larger than 0, these variables should be 31 considered. ADT is insignificant for both variations of IRI (IRI_diff) and longitudinal crack 32 (LCRACK_diff). In addition, it is also strongly correlated with TRUCKS. Hence, ADT is 33 also deleted. 34 IRI_diffFault_diffLCRACK_diffTCRACK_diffQ1−1.5⋅IQRQ3+1.5⋅IQR 1 Figure 1. Heatmap for variables’ correlation, green color represents positive correlation, blue 2 represents negative correlation, and grey color represents weak correlation 3 4 Among the four outputs, several interesting correlation trends can be observed. First, 5 the variation of IRI (IRI_diff) is positively correlated with both the level of IRI (collection 6 coefficient, ρ = 0.13) and the level of transverse cracking (TCRACK, ρ = 0.08). For this 7 modeling, a positive correlation suggests that IRI increases more rapidly when values for IRI 8 or TCRACK are higher. Additionally, it has been noticed that the variation of FAULT 9 (FAULT_diff) is strongly negatively corelated with the level of faulting (FAULT, ρ = -0.36) 10 and longitudinal cracking (LCRACK, ρ = -0.12). This suggests that faulting grows more 11 slowly when faulting or longitudinal cracking is more severe. The variation of LCRACK 12 (LCRACK_diff) and TCRACK (TCRACK_diff) are both strongly positively correlated with 13 IRI, LCRACK and TCRACK. From this correlation analysis, it can be concluded that the 14 deterioration of each single condition metric is statistically dependent on at least some of the 15 other condition metrics. This observation implies that it’s better to incorporate the state of all 16 condition metrics during the prediction of single condition metrics. 17 Table 8 summarizes the distribution of input variables for the neural network model. 18 As for four condition metrics, when their standard deviation and mean values are compared, 19 IRI has the smallest ratio between standard deviation and mean value, representing the IRI 20 distribution is more centralized. TCRACK has the largest ratio and its distribution is the 21 sparsest, which might be caused by the measurement quality. 22 Table 8. Statistical summary for input variables 1 AGECON TOTTHK TRUCKS IRI FAULT LCRACK TCRACK PRECIP FREEZE mean 22.2 248.8 1660 2.2 4.6 19.6 13.4 820.8 286.8 std 14.9 32.5 2559.7 0.8 2.6 37.8 26.8 112.2 156.6 min 0 150 1 0.7 0 0 0 650.5 317.8 max 71 559 14293 5 12 350 240.5 1077.2 840.8 *std: standard deviation 2 3 4.2 Normalization of training data 4 As listed in Table 8, each variable has a different magnitude. To avoid fitting 5 problems due to this, all input and output variables are normalized. Considering that both 6 input and output variables have condition metrics but for different years, all condition metrics 7 are normalized together. Suppose X is the input matrix with the dimension is , is 8 the number of input variables, and is the number of data points. is the output matrix. Its 9 dimension is , and is the number of output variables. The normalization process is 10 shown as follows: 11 (12) 12 (13) 13 (14) 14 Where and are normalized input and output values, where i represents the index of 15 data point, j is the input variable index and k is the output variable index. and are 16 mean values for jth input and kth output. and are standard deviations for jth input and 17 kth output. 18 19 4.3 Model framework 20 In recent years, feed-forward neural network (NN) models have been widely applied 21 for pavement deterioration prediction due to their superior performance when dealing with 22 strongly nonlinear relationships (Karlaftis and Badr 2015, Mazari and Rodriguez 2016, 23 Heidari et al. 2018, Hossain et al. 2019, Inkoom et al. 2019, Yao et al. 2019). A NN model 24 usually consists of an input layer, one or several hidden layers, and an output layer. Each 25 layer has several neurons. Hidden layers are usually used to describe the nonlinear 26 relationships via activation functions (also known as transfer functions). The connection 27 between two layers are usually through a linear mapping. A NN model can be mathematically 28 described as equation (15) 29 dx×ndxnYdy×ndyxj,i=Xj,i−Xj,σXj,j∉IRI,FAULT,LCRACK,TCRACK{}()xj,i=Xj,i−Yj,σYj,j∈IRI,FAULT,LCRACK,TCRACK{}()yk,i=Yk,i−Yk,σYk,k∈IRI,FAULT,LCRACK,TCRACK{}()xj,iyk,iXj,Yk,σXj,σYk, (15) 1 whererepresents the prediction values. is the linear transformation 2 matrix, and is the activation function. Common nonlinear activation 3 functions include ReLU, sigmoid, and Tanh. Their ouptut sets are , and . 4 ReLU is often used for regression problems while the other two are more commonly applied 5 for classification problems (Bishop 2006). In terms of the current regression problem, ReLU 6 (equation (16)) is applied for the hidden layers and a linear function is applied for the output 7 layer. 8 (16) 9 The goal of training a NN model is to determine an optimal structure that can 10 minimize the loss funcion. A loss function is used to describe the relationship between the 11 predicted value and the true value . In this regression analysis, mean squared error is 12 chosen as the loss function (equation (17)), where is the weight value for each output 13 variable. 14 (17) 15 After normalization process introduced in section 4.2, the training problem for a NN 16 model (equation (15) and (17)) can be reformulated as, 17 (18) 18 (19) 19 In this case, the prediction value is , which is in the normalization space. Equation (20) is 20 applied to obtain the prediction value with a real magnitude. 21 , and (20) 22 When the dimension of the output variable is equal to 1, the NN is a single-output 23 model, which is very common in existing literatures. When is larger than 1, and 24 for all output variables, then this is a uniformly weighted multi-output model. When is 25 larger than 1, and are different for all output variables, it is a weighted multi-output NN 26 model, which is the focus of this analysis. 27 In order to incorporate weights for each output variable, insert equation (14) to 28 equation (19), then, 29 ˆY=σH+1WH+1σHWHσH−1WH−1...σ2W2σ1W1X()()()()()ˆYWhh=1,…H+1()σhh=1,…H+1()0,∞⎡⎣⎤⎦0,1()−1,1()f(x)=xx>00otherwise⎧⎨⎪⎩⎪ˆYYλk2min:LˆY,Y()=12k=1dy∑λk2ˆYk−Yk22=12k=1dy∑i=1n∑λk2ˆYk,i−Yk,i()2ˆy=σH+1WH+1σHWHσH−1WH−1...σ2W2σ1W1x()()()()()min:Lˆy,y()=12k=1dy∑λk2ˆyk−yk22=12k=1dy∑i=1n∑λk2ˆyk,i−yk,i()2ˆyˆyk,i=ˆYk,i−Yk,σYk,ˆYk,i=Yk,+ˆyk,i⋅σYk,dydyλk=1dyλk (21) 1 By changing the stadard deviaiton during the normalization process, weights in the loss 2 function can be incorporated. 3 4 4.4 Weighted multi-output NN model for rigid pavements 5 In terms of the NN model structure for rigid pavement deterioration, the number of 6 hidden layers can be determined first. As mentioned by Heaton (2008), one hidden layer can 7 approximate any function that contains a continuous mapping from one finite space to 8 another. With the increase of hidden layers, the NN model could have an excellent fitting 9 performance, but it may decrease the prediction capacity due to overfitting. Since the dataset 10 applied in this work contains only around 2,500 training data points, one hidden layer should 11 work well enough and minimize the risk of overfitting (Bishop 2006). In this case, the 12 training problem (equation (18)) can be rewritten as, 13 (22) 14 Where, the activation function for the output layer and hidden layer is chosen as a 15 linear function and ReLU, respectively. 16 The determination of weight matrixes is usually based on a back-propagation 17 algorithm. In terms of equation (22), the partial derivative for the loss function L to can 18 be expressed as, 19 (23) 20 , where (24) 21 where and is the number of neurons in the hidden layer. 22 When , namely, the model has a single output, each value in the linear transition 23 matrix is only influenced by the single output. However, when , are 24 influenced by all outputs. Therefore, the correlations among different output variables could 25 be incorporated during the training process. 26 27 4.5 Model training 28 In addition to the number of hidden layers, there are three other main parameters 29 concerning the training of NN models: 1). the number of neurons M in the hidden layer; 2). 30 The epoch number: one epoch represents one forward pass and one backward pass of all the 31 min:Lˆy,y()=12k=1dy∑i=1n∑λk2ˆyk,i−yk,i()2=12k=1dy∑i=1n∑λk2ˆYk,i−Yk,σYk,−Yk,i−Yk,σYk,⎛⎝⎜⎜⎞⎠⎟⎟2=12k=1dy∑i=1n∑ˆYk,i−Yk,σYk,λk−Yk,i−Yk,σYk,λk⎛⎝⎜⎜⎞⎠⎟⎟2 ˆy=σ2W2σ1W1x()()σ2σ1W1∂L∂W1=∂L∂ˆy∂ˆy∂f2∂f2∂W2f1∂f1∂W1x∂W1x∂W1=W2T∂L∂ˆy∂f1∂W1xxT∂L∂Wm,j1=Wm,2∂L∂ˆy∂f1∂W1xm,xj,∂L∂ˆyk,i=λk2ˆyk,i−ˆyk,i()m∈MMk=1Wm,j1k>1Wm,j1 training examples; 3). Batch size: the number of training examples in one forward/backward 1 pass. 2 As suggested by (Heaton 2008), the number of hidden neurons should be 2/3 the size 3 of the input layer, plus the size of the output layer, which equals to 10 in this analysis. But in 4 reality, considering that each dataset has different characteristics, 6 neuron numbers are 5 selected as candidates including 10, 16, 32, 48, 64 and 80. As for the epoch number, when it 6 is small, the NN model may not be well trained. When it is large, the NN model has the 7 potential to be over-trained. Here, three epoch numbers are considered, 50, 100, and 150. In 8 addition, early-stopping is incorporated in order to avoid overfitting. Batch size determines 9 the number of data points that can be trained for each time. When it is small, the NN model 10 has the potential to be over-trained. Here, three batch sizes (8, 32 and 64) are considered. 11 The optimal combination of these three parameters is determined by grid search based 12 on a 5-fold cross-validation. For each fold, all training data are randomly divided into three 13 groups, including 70% training data, 20% validation data, and 10% test data. The training 14 data and validation data are used to train the model and the test data is used to check the out-15 of-sample performance of the NN model. Theoretically, the dataset should be divided so that 16 the three partitions – training, validation, and test – are large enough that we have a 17 reasonable expectation that they are all representative of the same distribution. Practically, 18 the ratio selection is mainly based on the data amount. For example, when the number of total 19 data points is very large, 50% of the total data can be selected as the test data. When the 20 number of data points is very small, there is only train and validation data, and no test data at 21 all. Because, in this case, it is better to use as many training data as possible to fully train the 22 model. As for Iowa dataset, the number of total data points for rigid pavements is around 23 2,700, so 10% of the test data is selected to make sure the model can be fully trained. By 24 using the 5-fold cross-validation, the model can be tested by 5 different sets of test data. 25 The mean squared error (MSE) of the test data is chosen as the model performance 26 metric, which is the same as the loss function (equation (19)) for the training process. After a 27 5-fold cross-validation, the average MSE is considered to compare NN models with different 28 structure parameters. The model with the lowest average MSE was selected as the 29 representative NN model. 30 5. Results 31 In order to explore the influence of incorporating multi-output variables, and various 32 weights for each output variable on the model prediction performance, four single-output 33 models were trained for each variable (i.e. IRI, FAULT, LCRACK and TCRACK) along with 34 ten multi-output models with different sets of weights. The first multi-output model (listed as 35 NN0 in Table 9) uses a set of weights that are uniform, [1, 1, 1, 1]. The second model (NN1) 36 uses a set of weights equal to those used in the Iowa DOT PCI equation (equation (5)), i.e., 37 [10, 5, 4, 6]. The other eight weight sets (for the other eight models, NN2 – NN10) are shown 38 in Table 9 and explore how disparity of weighing among the four variables impacts model 39 results. 40 The training speed for single-output and multi-output models is related with the 41 training parameters, including the number of neurons, and the epoch number and the batch 42 number. When both types of models use the same parameters: Neuron=64, Epoch=50, 43 Batch=32, the training time for single-output and multi-output models are 26.3s and 29.8s, 44 respectively (the training time is based on 5-fold cross-validation). The difference between 1 these two models in terms of training speed is not large. 2 Tables 9 lists the optimal set of model parameters and the smallest MSEs found 3 during the grid search for these 4 single-output and 10 multi-output models, respectively. For 4 the multi-output models, the optimal NN structure is selected based on average MSE of PCI 5 prediction across a 5-fold cross-validation. Three things are immediately notable from this 6 table. First, we see that, the multi-output models have better prediction performances 7 (smallest MSEs) for IRI, FAULT, and TCRACK (see values with an asterisk in Table 9) 8 compared to the corresponding single variable models. Both types of models have similar 9 prediction performance for LCRACK. The improvement in prediction performance is 10 particularly notable for the prediction of FAULT (single variable model MSE = 1.069, model 11 NN1 MSE for FAULT=0.778, a 27% reduction) and for the highly influential prediction of 12 IRI (single variable model MSE = 0.044, NN1 MSE for IRI = 0.033, a 25% reduction). 13 14 Table 9 Optimal single-output model structure and prediction performance (first four rows) 15 and Optimal multi-output model structure and prediction performance (last 10 rows) 16 NN structure Mean Squared Error (MSE) Name Neuron Epoch Batch PCI IRI FAULT LCRACK TCRACK IRI 32 150 32 11.78 0.044 FAULT 80 150 8 1.069 LCRACK 80 50 64 468.9 TCRACK 32 50 8 137.9 Name Weights Neuron Epoch Batch PCI IRI FAULT LCRACK TCRACK NN0 [1, 1, 1, 1] 48 150 36 10.25* 0.034 0.764 482.3 133.1 NN1 [10, 5, 4, 6] 64 50 32 10.24* (0.035) 0.033* (2e-5) 0.778 (4e-4) 491.3 (0.3) 132.4* (0.04) NN2 [10, 1, 1, 1] 80 100 8 10.37 0.033* 0.800 489.2 141.4 NN3 [50, 1, 1, 1] 48 150 8 10.53 0.033* 0.859 511.6 153.6 NN4 [1, 10, 1, 1] 80 50 8 11.09 0.038 0.703 501.0 140.1 NN5 [1, 50, 1, 1] 64 150 8 12.72 0.051 0.690* 519.1 151.7 NN6 [1, 1, 10, 1] 48 100 8 11.03 0.038 0.827 468.1* 136.4 NN7 [1, 1, 50, 1] 80 150 8 12.62 0.047 0.926 470.2* 158.1 NN8 [1, 1, 1, 10] 48 50 8 10.68 0.039 0.835 500.8 132.1* NN9 [1, 1, 1, 50] 80 100 8 12.04 0.052 0.947 540.5 132.3* *Smallest observed MSE for that variable 17 18 Secondly, and more importantly, from Table 9 it is clear that the PCI prediction 19 performance associated with either multi-output model NN0 or NN1 (MSE = 10.24) is better 20 for than for the PCI prediction derived from the four optimal single output models (MSE = 21 11.78). As mentioned in section 4.4, the training process for multi-output models 22 incorporates the correlations among variables. This correlation could improve the prediction 23 performance for a single variable and should improve the prediction of a multi-variable 24 metric like PCI. 25 Finally, although the multi-output models can outperform single metric models in 26 predicting PCI, because PCI is a composite performance metric, when it is optimal, it does 27 not mean all other single performance metrics are optimal. 28 Due to the existence of uncertainty in a multi-fold validation, the values listed in 29 Tables 9 could change if the same model was regenerated. In order to compare different 30 models statistically, a bootstrap analysis was conducted to estimate the standard deviation for 31 the average MSE obtained by 5-fold cross-validation. Because this uncertainty in MSE 32 derives solely from the random partitioning of training and validation data, all 10 multi-1 output models have the same source of uncertainty. Therefore, a bootstrap analysis was 2 conducted only for the NN1 model. The bootstrap sample size was 100. The standard 3 deviation in MSE for each variable is listed in Table 9 (i.e. values in parentheses). Using 4 these standard deviation values, statistically significant differences among the lowest 5 observed MSE values was tested using a two-sample t-test assuming equal variance and a 6 confidence level of 95%. The MSE values highlighted in Table 9 (asterisk) for PCI, IRI, 7 FAULT, LCRACK and TCRACK are the lowest values. When multiple values in a given 8 column have an asterisk, they are not statistically different. 9 In terms of PCI prediction, normally, the NN1 model would be expected to have the 10 best prediction performance since its weights resemble the PCI formula. However, both NN0 11 and NN1 models have the smallest MSE. This is because in the current PCI formula, the 12 coefficients (weightings) for all four variables are similar in magnitude. Figure 2 shows the 13 comparison between true PCI and predicted PCI based on out-of-sample data for NN1 model, 14 which presents a good prediction performance. 15 16 17 Figure 2. Comparison between true PCI and predicted PCI based on out-of-sample 18 data when weights are [10, 5, 4, 6] (NN1 model). 19 20 The models NN2 to NN9 have two larger weights (10x or 50x) for IRI, FAULT, 21 LCRACK, TCRACK, respectively. For each of these models, the prediction performance of 22 the variable with the largest weight (highlighted in bold) and the PCI prediction performance 23 are emphasized. Generally, when an output is more heavily weighted, the training process 24 emphasizes this output. As shown in equation (21), when the output is heavily weighted, 25 the minimization process will concentrate more on improving the prediction performance on 26 to minimize the total weighed loss value. When this occurs, prediction performance for 27 that focal output would be expected to be higher than for an NN model with a smaller weight 28 on that same output. It is important to note that models NN2 to NN9 are not expected to or 29 intended to improve on the prediction of PCI compared to NN1. Instead the role of these 30 models is to see how the multi-output model algorithm behaves as weighting shifts. In this 31 case, where correlation is present, these also serve to demonstrate both the ability of a multi-32 output model to outperform single-output models and the fact that as weights become more 33 YkYk disparate there may be an important tradeoff between performance in predicting the 1 composite metric (PCI) and the individual metrics. 2 As for IRI, when its weight is 2 times larger than other variables (NN1), its MSE 3 value is already the smallest. With an increase of IRI weight (models NN2 and NN3), IRI 4 prediction performance stays almost the same while MSE for PCI prediction increases very 5 slightly. The behavior of TCRACK prediction is similar to IRI. The NN1 model has a low 6 MSE for TCRACK that is not statistically resolvable from models NN9 and NN10. However, 7 TCRACK is influential in the prediction of PCI (only IRI is more influential). Therefore, as 8 the TCRACK weight grows relative to the other variables (NN9 and NN10), PCI prediction 9 performance deteriorates. For both FAULT and LCRACK, when they have larger weights 10 (i.e. 10 and 50), their prediction performances are better than NN0 and NN1. Hence, with the 11 increase of the weight for a single variable, its own prediction performance could be 12 improved. However, other variables’ prediction performances would decrease, as well as the 13 PCI prediction performance. To this last point, it can be seen that the difference in model 14 performance between a uniformly trained option and a weighted option increases as the 15 heterogeneity of the weights grows. Figure 3 presents the comparison between true and 16 predicted values for IRI, FAULT, LCRACK, and TCRACK. Generally, all four plots 17 demonstrate good model performance without significant systematic deviations in 18 correspondence between the predicted and actual values. Both LCRACK’s and TCRACK’s 19 plots show a small number of low range values (~50 to 100 m/km LCRACK and ~50 to 100 20 count/km TCRACK) that are underestimated by the model. The authors believe this occurs 21 for two reasons: (1). Data quality for these two variables is poor with standard deviations 22 roughly two times the mean for the samples of both; (2). Some crack-related maintenance 23 records may be missing. The record number of surface treatment (e.g. crack sealing) is 24 extremely small. Hence, some input and output parings may not be appropriate when there 25 exists a surface treatment but without a record. Underestimates at these low values should not 26 cause significant problems at these are well below current treatment threshold triggers. 27 28 (a). IRI (b). FAULT (c). LCRACK (d). TCRACK 1 Figure 3. Comparison between true and predicted values based on out-of-sample data 2 for (a). IRI (NN3 model); (b). FAULT (NN5 model); (c). LCRACK (NN7 model); (d). 3 TCRACK (NN9 model) 4 6. Conclusions 5 A new weighted multi-output neutral network model for the deterioration of rigid 6 pavements has been proposed based on the Iowa pavement management system database. 7 This new model simultaneously predicts IRI, faulting (FAULT), longitudinal cracking 8 (LCRACK) and transverse cracking (TCRACK) for pavements, providing convenience for 9 pavement management systems whose treatment decisions are based on composite, multi-10 condition metrics such as the pavement condition index (PCI). Considering the fact that 11 different condition metrics may have different importance levels, each metric to be predicted 12 is given a weight during the model training process. 13 Compared to single-output models, an appropriately weighted, multi-output model 14 (NN1) performs better at estimating PCI (13% lower MSE) than an estimate derived from 15 four optimal, single-output models. Furthermore, the multi-output model generates better 16 estimates than corresponding single-output models for three of the four individual metrics 17 considered – IRI (25% lower MSE), FAULT (27% lower) and TCRACK (4% lower). For the 18 fourth metric, LCRACK, the multi-output model has a 5% higher MSE than the single-output 19 model. The observed improvements derive directly from incorporating correlation 20 relationships among different variables within the multi-output model. The deterioration in 21 LCRACK prediction performance demonstrates the tradeoff that is made when 22 simultaneously estimating correlated metrics. 23 The results presented in this paper make it clear, that multi-output models can 24 improve prediction performance in cases where correlation exists. Furthermore, variable 25 weighting is important to achieve the optimal balance of prediction performance among the 26 various metrics. Analysis of extremely weighted models (NN2 to NN9) suggests that the 27 relevance of understanding the nature of this tradeoff will grow as weights become more 28 heterogeneous. 29 As for the future work to improve the prediction performance for pavement 30 deterioration, one aspect is to improve data quality, which is a common problem for data-31 driven deterioration models. The Iowa PMS provides one of the best pavement databases that 32 are available. However, there still exist some problems, such as large measurement errors, 1 missing maintenance history, and data entry errors. On one hand, these problems increase the 2 difficulty to generate the data for training the model, on the other hand, a large ratio of data 3 points has to be removed. To improve the model performance, it is necessary to obtain data 4 with a good quality. Another aspect is about the time series analysis. Currently, it is assumed 5 that the time step is just two years due to limited data. With the increase of available records, 6 it is necessary to test different time steps. As for the model structure, one potential drawback 7 for neural network model is its ‘black box’ feature, which means that the relationship 8 between input and output variables are not clear. In some research fields such as choice 9 models in transportation, some attempts have been made to integrate theoretical model and 10 data-driven models. In the future, this idea could also be applied in pavement deterioration 11 analysis. 12 13 Acknowledgements 14 This research was conducted as part of the Concrete Sustainability Hub at the 15 Massachusetts Institute of Technology, which is supported by the Portland Cement 16 Association and the Ready Mixed Concrete Research and Education Foundation. 17 18 Declarations of interest: none 19 20 21