Washington Rural Multi-lane Divided

Subasish Das and Choalun Ma

2018-11-12

# Step 1: Load R Packages 
### options(repos='http://cran.rstudio.com/')
#install.packages("astsa")
#install.packages('ggplot2')
#install.packages('forecast')
#install.packages('tseries')
#install.packages("data.table")

library(astsa)
library(forecast)
library(tseries)
library(zoo)
library(tseries)

library(data.table)
library(dplyr)
library(tidyr)
library(naniar)
library(stringr)
library(ggplot2)
library(DT)
library(lubridate)
library(ggpubr)


setwd("/scratch/user/cma16/Task4_Deliverable2/Process4/AllCrash/FacilityBased/")
load("./multi-lane_divided_WA_reduce_withCrash.rData")
mytype = 'RMD'
setwd(paste0("/scratch/user/cma16/Task4_Deliverable2/Process4/AllCrash/FacilityBased/",mytype))

df_RMD <- W_mun_med
dim(df_RMD)
## [1] 3784320      64
df_RMD$spd_av = 3600*df_RMD$DISTANCE/df_RMD$Travel_TIME_ALL_VEHICLES
df_RMD$spd_pv = 3600*df_RMD$DISTANCE/df_RMD$Travel_TIME_PASSENGER_VEHICLES
df_RMD$spd_ft = 3600*df_RMD$DISTANCE/df_RMD$Travel_TIME_FREIGHT_TRUCKS

### Month, Day
df_RMD$date <- as.character(df_RMD$DATE)
df_RMD$date <- str_pad(df_RMD$DATE, 8, pad = "0")
df_RMD$Month <- substr(df_RMD$date, start = 1, stop = 2)
df_RMD$Day   <- substr(df_RMD$date, start = 3, stop = 4)
df_RMD$Year  <- substr(df_RMD$date, start = 5, stop = 8)
df_RMD$MonthDay <- paste0(df_RMD$Month,"_", df_RMD$Day)
head(df_RMD)
##            TimeStamp       TMC     V1    DATE EPOCH15
## 1:  114N05428_0101_0 114N05428 425569 1012015       0
## 2:  114N05428_0101_1 114N05428 425570 1012015       1
## 3: 114N05428_0101_10 114N05428 425579 1012015      10
## 4: 114N05428_0101_11 114N05428 425580 1012015      11
## 5: 114N05428_0101_12 114N05428 425581 1012015      12
## 6: 114N05428_0101_13 114N05428 425582 1012015      13
##    Travel_TIME_ALL_VEHICLES Travel_TIME_PASSENGER_VEHICLES
## 1:                       NA                             NA
## 2:                       NA                             NA
## 3:                       NA                             NA
## 4:                       NA                             NA
## 5:                       NA                             NA
## 6:                       NA                             NA
##    Travel_TIME_FREIGHT_TRUCKS NP ADMIN_LEVE ADMIN_LE_1 ADMIN_LE_2 DISTANCE
## 1:                         NA  N        USA Washington       King   8.3187
## 2:                         NA  N        USA Washington       King   8.3187
## 3:                         NA  N        USA Washington       King   8.3187
## 4:                         NA  N        USA Washington       King   8.3187
## 5:                         NA  N        USA Washington       King   8.3187
## 6:                         NA  N        USA Washington       King   8.3187
##    ROAD_NUMBE ROAD_NAME LATITUDE LONGITUDE ROAD_DIREC  ORN_FID    FID_1
## 1:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
## 2:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
## 3:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
## 4:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
## 5:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
## 6:       US-2           47.71281 -121.2902  Westbound 2638.708 1096.573
##    ACCESS LSHL_TY2 LSHL_TYP MED_TYPE NHS_IND PRK_ZNE RSHL_TY2 RSHL_TYP
## 1:      P                 A                Y                         A
## 2:      P                 A                Y                         A
## 3:      P                 A                Y                         A
## 4:      P                 A                Y                         A
## 5:      P                 A                Y                         A
## 6:      P                 A                Y                         A
##    SURF_TYP SURF_TY2 TERRAIN COMP_DIR COUNTY FUNC_CLS MEDBARTY ST_FUNC
## 1:        A                M        E     17       43               R1
## 2:        A                M        E     17       43               R1
## 3:        A                M        E     17       43               R1
## 4:        A                M        E     17       43               R1
## 5:        A                M        E     17       43               R1
## 6:        A                M        E     17       43               R1
##    RTE_NBR          HPMS ROAD_INV SPD_LIMT BEGMP ENDMP LSHLDWID    MEDWID
## 1:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
## 2:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
## 3:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
## 4:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
## 5:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
## 6:       2 5002172008272      002 55.07084 52.16 60.54 6.073922 0.8364139
##    NO_LANE1 NO_LANE2 NO_LANES RSHLDWID  RSHL_WD2   SEG_LNG lanewid
## 1:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
## 2:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
## 3:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
## 4:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
## 5:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
## 6:  1.80421        1        3 4.867425 0.1539141 0.2689276  12.047
##    rdwy_wd1 rdwy_wd2 rdwy_wid     AADT      mvmt rodwycls ORN_FID_1 Total
## 1: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
## 2: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
## 3: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
## 4: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
## 5: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
## 6: 33.41593 0.489148 33.90508 5086.424 0.4998998        7  2638.708    16
##    Fatal Injury PDO DAYMTH Crash spd_av spd_pv spd_ft     date Month Day
## 1:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
## 2:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
## 3:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
## 4:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
## 5:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
## 6:     1      5  10   0101     0     NA     NA     NA 01012015    01  01
##    Year MonthDay
## 1: 2015    01_01
## 2: 2015    01_01
## 3: 2015    01_01
## 4: 2015    01_01
## 5: 2015    01_01
## 6: 2015    01_01
day1<- df_RMD[,-c(1)] %>% group_by(MonthDay) %>% summarize(Speed_All_Mean=mean(spd_av, na.rm=TRUE))
day1
## # A tibble: 365 x 2
##    MonthDay Speed_All_Mean
##    <chr>             <dbl>
##  1 01_01              58.7
##  2 01_02              56.4
##  3 01_03              57.6
##  4 01_04              54.7
##  5 01_05              55.9
##  6 01_06              56.7
##  7 01_07              56.5
##  8 01_08              56.8
##  9 01_09              56.9
## 10 01_10              57.6
## # ... with 355 more rows
# Step 2: Examine Data
speed_clean <- tsclean(ts(day1$Speed_All_Mean))
plot.ts(speed_clean)

# ggplot() + geom_line(data = Q1, aes(x = TimeStamp, y = speed_clean)) + ylab('Cleaned Speed Records')

day1$cnt_ma = ma(speed_clean, order=7) # using the clean count with no outliers
day1$cnt_ma30 = ma(speed_clean, order=30)
# Step 3: Decompose Your Data
count_ma = ts(na.omit(speed_clean), frequency=30)
decomp = stl(count_ma, s.window="periodic")
deseasonal_cnt <- seasadj(decomp)
plot(decomp)

# Step 4: Stationarity
# statinary test
adf.test(count_ma, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  count_ma
## Dickey-Fuller = -2.8345, Lag order = 7, p-value = 0.2247
## alternative hypothesis: stationary
adf.test(deseasonal_cnt, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  deseasonal_cnt
## Dickey-Fuller = -2.6891, Lag order = 7, p-value = 0.2861
## alternative hypothesis: stationary
d1 = diff(deseasonal_cnt)
adf.test(d1, alternative = "stationary")
## Warning in adf.test(d1, alternative = "stationary"): p-value smaller than
## printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  d1
## Dickey-Fuller = -10.33, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
# Step 5: Autocorrelations and Choosing Model Order
# check ACF and PACF
acf2(count_ma)

##          ACF  PACF
##   [1,]  0.49  0.49
##   [2,]  0.13 -0.15
##   [3,]  0.07  0.10
##   [4,]  0.06  0.00
##   [5,]  0.13  0.13
##   [6,]  0.39  0.37
##   [7,]  0.66  0.48
##   [8,]  0.36 -0.16
##   [9,]  0.06 -0.07
##  [10,]  0.03  0.00
##  [11,]  0.00 -0.10
##  [12,]  0.04 -0.08
##  [13,]  0.33  0.18
##  [14,]  0.62  0.30
##  [15,]  0.35 -0.06
##  [16,]  0.05 -0.01
##  [17,] -0.03 -0.11
##  [18,] -0.01  0.06
##  [19,]  0.05 -0.01
##  [20,]  0.31  0.03
##  [21,]  0.58  0.15
##  [22,]  0.26 -0.18
##  [23,]  0.00  0.01
##  [24,] -0.03  0.01
##  [25,] -0.06 -0.05
##  [26,] -0.01 -0.04
##  [27,]  0.26  0.08
##  [28,]  0.54  0.13
##  [29,]  0.27 -0.03
##  [30,] -0.04 -0.09
##  [31,] -0.11 -0.09
##  [32,] -0.10  0.00
##  [33,] -0.03  0.02
##  [34,]  0.25  0.02
##  [35,]  0.49  0.06
##  [36,]  0.22 -0.03
##  [37,] -0.06 -0.01
##  [38,] -0.12 -0.04
##  [39,] -0.11  0.00
##  [40,] -0.06 -0.05
##  [41,]  0.20  0.04
##  [42,]  0.44  0.01
##  [43,]  0.19 -0.04
##  [44,] -0.08 -0.02
##  [45,] -0.13  0.04
##  [46,] -0.12  0.05
##  [47,] -0.07 -0.03
##  [48,]  0.19  0.00
##  [49,]  0.43  0.05
##  [50,]  0.18 -0.02
##  [51,] -0.09  0.02
##  [52,] -0.14 -0.05
##  [53,] -0.12  0.02
##  [54,] -0.08 -0.02
##  [55,]  0.15 -0.05
##  [56,]  0.41  0.08
##  [57,]  0.15 -0.07
##  [58,] -0.12  0.00
##  [59,] -0.17  0.00
##  [60,] -0.16 -0.02
##  [61,] -0.10 -0.03
##  [62,]  0.13  0.00
##  [63,]  0.38  0.03
##  [64,]  0.13 -0.05
##  [65,] -0.15 -0.02
##  [66,] -0.18  0.00
##  [67,] -0.18 -0.04
##  [68,] -0.13  0.00
##  [69,]  0.12  0.01
##  [70,]  0.35  0.01
##  [71,]  0.10 -0.02
##  [72,] -0.14  0.05
##  [73,] -0.20 -0.04
##  [74,] -0.20 -0.04
##  [75,] -0.13  0.02
##  [76,]  0.11  0.02
##  [77,]  0.33 -0.02
##  [78,]  0.09 -0.02
##  [79,] -0.15  0.03
##  [80,] -0.19  0.00
##  [81,] -0.20 -0.02
##  [82,] -0.14 -0.04
##  [83,]  0.08 -0.01
##  [84,]  0.29 -0.01
##  [85,]  0.07  0.00
##  [86,] -0.17  0.00
##  [87,] -0.20  0.00
##  [88,] -0.21  0.01
##  [89,] -0.17 -0.05
##  [90,]  0.06 -0.01
##  [91,]  0.28  0.00
##  [92,]  0.05 -0.01
##  [93,] -0.17  0.03
##  [94,] -0.22 -0.02
##  [95,] -0.20  0.04
##  [96,] -0.16 -0.03
##  [97,]  0.07  0.03
##  [98,]  0.26 -0.03
##  [99,]  0.05 -0.01
## [100,] -0.17 -0.01
## [101,] -0.22 -0.02
## [102,] -0.21  0.03
## [103,] -0.16 -0.01
## [104,]  0.05 -0.02
## [105,]  0.25  0.00
## [106,]  0.04  0.00
## [107,] -0.18 -0.01
## [108,] -0.21  0.03
## [109,] -0.19  0.01
## [110,] -0.14  0.03
## [111,]  0.04 -0.06
## [112,]  0.23 -0.01
## [113,]  0.04 -0.01
## [114,] -0.17 -0.01
## [115,] -0.22 -0.02
## [116,] -0.20 -0.01
## [117,] -0.15  0.02
## [118,]  0.03 -0.02
## [119,]  0.22 -0.01
## [120,]  0.03  0.00
acf2(deseasonal_cnt)

##          ACF  PACF
##   [1,]  0.50  0.50
##   [2,]  0.13 -0.16
##   [3,]  0.08  0.11
##   [4,]  0.07 -0.01
##   [5,]  0.13  0.13
##   [6,]  0.40  0.38
##   [7,]  0.68  0.49
##   [8,]  0.37 -0.17
##   [9,]  0.06 -0.09
##  [10,]  0.03  0.01
##  [11,]  0.01 -0.10
##  [12,]  0.04 -0.07
##  [13,]  0.35  0.19
##  [14,]  0.63  0.29
##  [15,]  0.35 -0.08
##  [16,]  0.05 -0.01
##  [17,] -0.03 -0.12
##  [18,] -0.01  0.05
##  [19,]  0.05  0.00
##  [20,]  0.31  0.02
##  [21,]  0.59  0.16
##  [22,]  0.27 -0.18
##  [23,] -0.01  0.00
##  [24,] -0.03  0.01
##  [25,] -0.05 -0.04
##  [26,]  0.00 -0.04
##  [27,]  0.27  0.08
##  [28,]  0.55  0.13
##  [29,]  0.28 -0.04
##  [30,] -0.05 -0.10
##  [31,] -0.10 -0.09
##  [32,] -0.10 -0.01
##  [33,] -0.03  0.02
##  [34,]  0.25  0.01
##  [35,]  0.49  0.06
##  [36,]  0.22 -0.01
##  [37,] -0.06  0.01
##  [38,] -0.12 -0.04
##  [39,] -0.11  0.00
##  [40,] -0.06 -0.05
##  [41,]  0.20  0.04
##  [42,]  0.45  0.02
##  [43,]  0.20 -0.02
##  [44,] -0.08 -0.01
##  [45,] -0.13  0.04
##  [46,] -0.12  0.03
##  [47,] -0.07 -0.04
##  [48,]  0.19  0.00
##  [49,]  0.44  0.05
##  [50,]  0.18 -0.03
##  [51,] -0.09  0.04
##  [52,] -0.14 -0.05
##  [53,] -0.13  0.01
##  [54,] -0.08 -0.03
##  [55,]  0.15 -0.04
##  [56,]  0.41  0.09
##  [57,]  0.16 -0.07
##  [58,] -0.12  0.01
##  [59,] -0.16 -0.01
##  [60,] -0.16 -0.03
##  [61,] -0.11 -0.04
##  [62,]  0.13  0.00
##  [63,]  0.38  0.04
##  [64,]  0.13 -0.06
##  [65,] -0.16 -0.02
##  [66,] -0.18  0.01
##  [67,] -0.18 -0.02
##  [68,] -0.13  0.00
##  [69,]  0.11  0.01
##  [70,]  0.35  0.01
##  [71,]  0.11 -0.02
##  [72,] -0.14  0.06
##  [73,] -0.19 -0.03
##  [74,] -0.19 -0.03
##  [75,] -0.13  0.02
##  [76,]  0.10  0.01
##  [77,]  0.33 -0.04
##  [78,]  0.09 -0.01
##  [79,] -0.15  0.04
##  [80,] -0.19 -0.01
##  [81,] -0.20  0.01
##  [82,] -0.15 -0.04
##  [83,]  0.07 -0.01
##  [84,]  0.29 -0.02
##  [85,]  0.08  0.00
##  [86,] -0.17  0.00
##  [87,] -0.20  0.00
##  [88,] -0.21  0.01
##  [89,] -0.18 -0.07
##  [90,]  0.05 -0.01
##  [91,]  0.27  0.00
##  [92,]  0.05 -0.02
##  [93,] -0.17  0.04
##  [94,] -0.22 -0.03
##  [95,] -0.20  0.05
##  [96,] -0.16 -0.03
##  [97,]  0.06  0.05
##  [98,]  0.26 -0.05
##  [99,]  0.05 -0.01
## [100,] -0.17 -0.01
## [101,] -0.22 -0.02
## [102,] -0.20  0.05
## [103,] -0.15  0.01
## [104,]  0.05 -0.02
## [105,]  0.25 -0.01
## [106,]  0.04 -0.01
## [107,] -0.18 -0.02
## [108,] -0.21  0.04
## [109,] -0.19  0.02
## [110,] -0.15  0.01
## [111,]  0.04 -0.03
## [112,]  0.24 -0.02
## [113,]  0.04 -0.01
## [114,] -0.18 -0.02
## [115,] -0.21 -0.01
## [116,] -0.20 -0.01
## [117,] -0.15  0.03
## [118,]  0.03 -0.01
## [119,]  0.22 -0.02
## [120,]  0.02 -0.01

Seasonility Not in Consideration

# Step 6: Fitting an ARIMA model
auto.arima(deseasonal_cnt, seasonal=FALSE)
## Series: deseasonal_cnt 
## ARIMA(1,1,5) 
## 
## Coefficients:
##           ar1    ma1      ma2      ma3     ma4     ma5
##       -0.5514  0.015  -0.9499  -0.4531  0.1934  0.3637
## s.e.   0.1155  0.114   0.0689   0.1084  0.0468  0.0679
## 
## sigma^2 estimated as 0.4451:  log likelihood=-367.57
## AIC=749.14   AICc=749.45   BIC=776.42
# Step 7: Evaluate and Iterate
# (try different model)
fit<-auto.arima(deseasonal_cnt, seasonal=FALSE)
tsdisplay(residuals(fit), lag.max=45, main='Model Residuals [Seasonality not considered]')

# step 8 forcast
fcast <- forecast(fit, h=30)
plot(fcast)

Seasonility in Consideration

# Step 6: Fitting an ARIMA model
auto.arima(deseasonal_cnt, seasonal=TRUE)
## Series: deseasonal_cnt 
## ARIMA(1,1,3)(0,0,2)[30] 
## 
## Coefficients:
##           ar1     ma1      ma2      ma3     sma1     sma2
##       -0.7027  0.1690  -0.7765  -0.2002  -0.3791  -0.2691
## s.e.   0.1942  0.2113   0.0833   0.1163   0.0641   0.0753
## 
## sigma^2 estimated as 0.4219:  log likelihood=-364.12
## AIC=742.24   AICc=742.56   BIC=769.52
# Step 7: Evaluate and Iterate
# (try different model)
fit<-auto.arima(deseasonal_cnt, seasonal=TRUE)
tsdisplay(residuals(fit), lag.max=45, main='Model Residuals [Seasonality considered]')

# step 8 forcast
fcast <- forecast(fit, h=30)
plot(fcast)