setwd("/Users/subasishdas1/Dropbox/---- TRB_2017 ----/SHRP2 Inclement")
aa <- read.csv("FINAL_4a.csv")
names(aa)
## [1] "CR.A.AMSHNUM_AG" "ID_A" "CRASHNUM1"
## [4] "DataType" "HIGHESTINJ" "Category"
## [7] "Category1" "Vis_Score" "SKIDNUMBER"
## [10] "RCISLDWTH1" "RCISLDWTH2" "AVG_SH_WID"
## [13] "RCIAADT" "RCIAVGTFCT" "RCIMAXSPD"
## [16] "LIGHTCOND" "WEATHCOND" "DIV_UNDIV"
## [19] "AGE_DRPED" "RDSURFCOND" "YEAR"
## [22] "CNTOFVEH" "ID" "Airport_Na"
## [25] "Airport_Nu" "LATITUDE" "LONGITUDE"
## [28] "CRSHCAUSE1" "CNTOFINJ" "CNTOFFATL"
## [31] "CNTOFSVINJ" "RCIFUNCLAS" "RCIACC"
## [34] "CR_LATITUDE" "CR_LONGITUDE" "RCISLDWTH3"
## [37] "RCIMEDWDTH" "StartDate" "Start_Time"
## [40] "End_Date" "End_Time" "Duration"
## [43] "HIGHESTINJ1" "Vis_Score1" "SKIDNUMBER1"
## [46] "RCIAADT1" "RCIAVGTFCT1" "RCIMAXSPD1"
## [49] "LIGHTCOND1" "WEATHCOND1" "DIV_UNDIV1"
## [52] "AGE_DRPED1" "RDSURFCOND1" "CNTOFVEH1"
## [55] "RCIFUNCLAS1" "RCIMEDWDTH1"
aa1 <- aa[c(44, 43, 15, 13, 9, 12, 19, 14, 48, 46, 45, 49, 51, 52)]
names(aa1)
## [1] "Vis_Score1" "HIGHESTINJ1" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "AVG_SH_WID" "AGE_DRPED" "RCIAVGTFCT" "RCIMAXSPD1" "RCIAADT1"
## [11] "SKIDNUMBER1" "LIGHTCOND1" "DIV_UNDIV1" "AGE_DRPED1"
summary(aa1)
## Vis_Score1 HIGHESTINJ1 RCIMAXSPD RCIAADT
## Excellent:11666 Fatal : 181 Min. :25.00 Min. : 0
## Medium : 9815 Injury :11293 1st Qu.:40.00 1st Qu.: 26000
## Poor : 2138 No Injury:12145 Median :45.00 Median : 40500
## Mean :45.25 Mean : 60909
## 3rd Qu.:50.00 3rd Qu.: 60500
## Max. :70.00 Max. :304000
##
## SKIDNUMBER AVG_SH_WID AGE_DRPED RCIAVGTFCT
## Min. : 0.00 Min. : 0.000 Min. : 15.00 Min. : 0.000
## 1st Qu.:33.00 1st Qu.: 1.000 1st Qu.: 25.00 1st Qu.: 3.300
## Median :36.00 Median : 3.000 Median : 38.00 Median : 4.530
## Mean :35.29 Mean : 3.844 Mean : 39.57 Mean : 5.171
## 3rd Qu.:40.00 3rd Qu.: 6.000 3rd Qu.: 51.00 3rd Qu.: 6.260
## Max. :62.00 Max. :19.500 Max. :108.00 Max. :31.770
##
## RCIMAXSPD1 RCIAADT1 SKIDNUMBER1
## >60 : 2600 >124,999 :3131 > 50 : 250
## 0-30 : 1827 0-9999 : 871 20-30: 1738
## 30-45:15013 10,000-34,999 :8127 30-40:15920
## 45-60: 4179 35,000-54,999 :7226 40-50: 4702
## 55,000-124,999:4050 NA's : 1009
## NA's : 214
##
## LIGHTCOND1 DIV_UNDIV1 AGE_DRPED1
## Dark(No Street light): 818 divided :16111 20-29 :6190
## Dark(Street Light) : 6417 undivided: 7192 30-39 :4555
## Dawn/Dusk : 1195 unknown : 316 40-49 :4358
## Daylight :15189 50-59 :3455
## 60-69 :1914
## (Other):2686
## NA's : 461
aa2 <- aa1[c(1:8)]
names(aa2)
## [1] "Vis_Score1" "HIGHESTINJ1" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "AVG_SH_WID" "AGE_DRPED" "RCIAVGTFCT"
library(ordinal)
m1 <- clm(Vis_Score1 ~HIGHESTINJ1+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCIAVGTFCT+
AVG_SH_WID+ AGE_DRPED, data = aa2)
## Warning: (2) Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## In addition: Absolute and relative convergence criteria were met
summary(m1)
## formula:
## Vis_Score1 ~ HIGHESTINJ1 + RCIMAXSPD + RCIAADT + SKIDNUMBER + RCIAVGTFCT + AVG_SH_WID + AGE_DRPED
## data: aa2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 23619 -21934.08 43888.17 7(0) 6.38e-07 4.6e+12
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## HIGHESTINJ1Injury 2.651e-01 1.515e-01 1.750 0.0801 .
## HIGHESTINJ1No Injury 2.630e-01 1.514e-01 1.736 0.0825 .
## RCIMAXSPD 8.029e-03 1.984e-03 4.048 5.17e-05 ***
## RCIAADT -7.073e-07 2.789e-07 -2.536 0.0112 *
## SKIDNUMBER -6.498e-03 1.544e-03 -4.209 2.56e-05 ***
## RCIAVGTFCT 9.086e-04 4.623e-03 0.197 0.8442
## AVG_SH_WID 2.722e-03 5.076e-03 0.536 0.5918
## AGE_DRPED -5.583e-03 7.675e-04 -7.273 3.50e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## Excellent|Medium 0.1235 0.1710 0.722
## Medium|Poor 2.4617 0.1721 14.305
anova(m1, update(m1, scale = ~.-(RCIAVGTFCT+AVG_SH_WID)))
## Warning: (-2) Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## In addition: maximum number of consecutive Newton modifications reached
## Likelihood ratio tests of cumulative link models:
##
## formula:
## m1 Vis_Score1 ~ HIGHESTINJ1 + RCIMAXSPD + RCIAADT + SKIDNUMBER + RCIAVGTFCT + AVG_SH_WID + AGE_DRPED
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) Vis_Score1 ~ HIGHESTINJ1 + RCIMAXSPD + RCIAADT + SKIDNUMBER + RCIAVGTFCT + AVG_SH_WID + AGE_DRPED
## scale:
## m1 ~1
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) ~. - (RCIAVGTFCT + AVG_SH_WID)
## link: threshold:
## m1 logit flexible
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) logit flexible
##
## no.par AIC logLik
## m1 10 43888 -21934
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) 18 34803 -17383
## LR.stat df Pr(>Chisq)
## m1
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) 9101.6 8 < 2.2e-16
##
## m1
## update(m1, scale = ~. - (RCIAVGTFCT + AVG_SH_WID)) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- clm(Vis_Score1 ~HIGHESTINJ1+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ AGE_DRPED, data = aa2)
## Warning: (2) Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## In addition: Absolute and relative convergence criteria were met
summary(m2)
## formula:
## Vis_Score1 ~ HIGHESTINJ1 + RCIMAXSPD + RCIAADT + SKIDNUMBER + AGE_DRPED
## data: aa2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 23619 -21934.26 43884.53 6(0) 6.71e-07 4.6e+12
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## HIGHESTINJ1Injury 2.650e-01 1.515e-01 1.750 0.08017 .
## HIGHESTINJ1No Injury 2.626e-01 1.514e-01 1.734 0.08289 .
## RCIMAXSPD 8.688e-03 1.638e-03 5.305 1.13e-07 ***
## RCIAADT -7.113e-07 2.727e-07 -2.609 0.00909 **
## SKIDNUMBER -6.786e-03 1.467e-03 -4.627 3.72e-06 ***
## AGE_DRPED -5.578e-03 7.675e-04 -7.268 3.65e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## Excellent|Medium 0.1277 0.1707 0.748
## Medium|Poor 2.4658 0.1718 14.352
library(MASS)
m <- polr(Vis_Score1 ~HIGHESTINJ1+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCIAVGTFCT+ AVG_SH_WID+ AGE_DRPED, data = aa2, na.action = na.omit, Hess=TRUE)
m
## Call:
## polr(formula = Vis_Score1 ~ HIGHESTINJ1 + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCIAVGTFCT + AVG_SH_WID + AGE_DRPED, data = aa2,
## na.action = na.omit, Hess = TRUE)
##
## Coefficients:
## HIGHESTINJ1Injury HIGHESTINJ1No Injury RCIMAXSPD
## 2.651082e-01 2.629668e-01 8.029317e-03
## RCIAADT SKIDNUMBER RCIAVGTFCT
## -7.072602e-07 -6.497780e-03 9.085730e-04
## AVG_SH_WID AGE_DRPED
## 2.721568e-03 -5.582666e-03
##
## Intercepts:
## Excellent|Medium Medium|Poor
## 0.1235173 2.4616696
##
## Residual Deviance: 43868.17
## AIC: 43888.17
newdat <- cbind(aa1, predict(m, aa2[c(2:8)], type = "probs"))
head(newdat)
## Vis_Score1 HIGHESTINJ1 RCIMAXSPD RCIAADT SKIDNUMBER AVG_SH_WID AGE_DRPED
## 1 Excellent No Injury 35 34000 36 4.5 31
## 2 Excellent Injury 30 21000 44 1.0 55
## 3 Excellent Injury 30 21000 44 1.0 47
## 4 Excellent Injury 45 22500 33 7.5 44
## 5 Excellent No Injury 45 22500 33 6.0 32
## 6 Excellent No Injury 45 22500 33 6.0 40
## RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1 LIGHTCOND1
## 1 8.35 30-45 10,000-34,999 30-40 Dark(No Street light)
## 2 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 3 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 4 8.33 30-45 10,000-34,999 30-40 Daylight
## 5 8.33 30-45 10,000-34,999 30-40 Daylight
## 6 8.33 30-45 10,000-34,999 30-40 Daylight
## DIV_UNDIV1 AGE_DRPED1 Excellent Medium Poor
## 1 undivided 30-39 0.4976801 0.4135601 0.08875980
## 2 undivided 50-59 0.5535496 0.3742370 0.07221346
## 3 undivided 40-49 0.5424877 0.3822488 0.07526348
## 4 divided 40-49 0.4862752 0.4212042 0.09252061
## 5 divided 30-39 0.4711119 0.4311374 0.09775064
## 6 divided 40-49 0.4822525 0.4238656 0.09388186
names(newdat)
## [1] "Vis_Score1" "HIGHESTINJ1" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "AVG_SH_WID" "AGE_DRPED" "RCIAVGTFCT" "RCIMAXSPD1" "RCIAADT1"
## [11] "SKIDNUMBER1" "LIGHTCOND1" "DIV_UNDIV1" "AGE_DRPED1" "Excellent"
## [16] "Medium" "Poor"
library(reshape2)
library(ggplot2)
lnewdat <- melt(newdat, id.vars = c("Vis_Score1", "HIGHESTINJ1", "RCIMAXSPD", "RCIAADT", "SKIDNUMBER", "AVG_SH_WID", "AGE_DRPED", "RCIAVGTFCT", "RCIMAXSPD1", "RCIAADT1", "SKIDNUMBER1", "LIGHTCOND1", "DIV_UNDIV1", "AGE_DRPED1"), variable.name = "Level", value.name="Probability")
head(lnewdat)
## Vis_Score1 HIGHESTINJ1 RCIMAXSPD RCIAADT SKIDNUMBER AVG_SH_WID AGE_DRPED
## 1 Excellent No Injury 35 34000 36 4.5 31
## 2 Excellent Injury 30 21000 44 1.0 55
## 3 Excellent Injury 30 21000 44 1.0 47
## 4 Excellent Injury 45 22500 33 7.5 44
## 5 Excellent No Injury 45 22500 33 6.0 32
## 6 Excellent No Injury 45 22500 33 6.0 40
## RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1 LIGHTCOND1
## 1 8.35 30-45 10,000-34,999 30-40 Dark(No Street light)
## 2 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 3 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 4 8.33 30-45 10,000-34,999 30-40 Daylight
## 5 8.33 30-45 10,000-34,999 30-40 Daylight
## 6 8.33 30-45 10,000-34,999 30-40 Daylight
## DIV_UNDIV1 AGE_DRPED1 Level Probability
## 1 undivided 30-39 Excellent 0.4976801
## 2 undivided 50-59 Excellent 0.5535496
## 3 undivided 40-49 Excellent 0.5424877
## 4 divided 40-49 Excellent 0.4862752
## 5 divided 30-39 Excellent 0.4711119
## 6 divided 40-49 Excellent 0.4822525
names(lnewdat)[2] <- "Severity"
head(lnewdat)
## Vis_Score1 Severity RCIMAXSPD RCIAADT SKIDNUMBER AVG_SH_WID AGE_DRPED
## 1 Excellent No Injury 35 34000 36 4.5 31
## 2 Excellent Injury 30 21000 44 1.0 55
## 3 Excellent Injury 30 21000 44 1.0 47
## 4 Excellent Injury 45 22500 33 7.5 44
## 5 Excellent No Injury 45 22500 33 6.0 32
## 6 Excellent No Injury 45 22500 33 6.0 40
## RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1 LIGHTCOND1
## 1 8.35 30-45 10,000-34,999 30-40 Dark(No Street light)
## 2 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 3 8.35 0-30 10,000-34,999 40-50 Dark(Street Light)
## 4 8.33 30-45 10,000-34,999 30-40 Daylight
## 5 8.33 30-45 10,000-34,999 30-40 Daylight
## 6 8.33 30-45 10,000-34,999 30-40 Daylight
## DIV_UNDIV1 AGE_DRPED1 Level Probability
## 1 undivided 30-39 Excellent 0.4976801
## 2 undivided 50-59 Excellent 0.5535496
## 3 undivided 40-49 Excellent 0.5424877
## 4 divided 40-49 Excellent 0.4862752
## 5 divided 30-39 Excellent 0.4711119
## 6 divided 40-49 Excellent 0.4822525
## Max Speed
library(ggplot2)
m1 <- ggplot(lnewdat, aes(x = RCIMAXSPD, y = Probability, colour = Level)) +
geom_point(alpha = 1/50, size=1) + facet_grid(~ Severity, labeller="label_both")+theme_bw()+
aes(colour = Level) + stat_smooth()+labs(x="Max. Speed (mph) \n (a)")+
theme(legend.position="none")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
## AADT
library(ggplot2)
m2 <- ggplot(lnewdat, aes(x = RCIAADT, y = Probability, colour = Level)) +
geom_point(alpha = 1/50, size=1) + facet_grid(~ Severity, labeller="label_both")+theme_bw()+
aes(colour = Level) + stat_smooth()+labs(x="AADT (vpd) \n (b)")+
theme(legend.position="none")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
## SKIDNUMBER
library(ggplot2)
m3 <- ggplot(lnewdat, aes(x = SKIDNUMBER, y = Probability, colour = Level)) +
geom_point(alpha = 1/50, size=1) + facet_grid(~ Severity, labeller="label_both")+theme_bw()+
aes(colour = Level) + stat_smooth()+ labs(x="Skid Number \n (c)")+
theme(legend.position="none")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
## Age
library(ggplot2)
m4 <- ggplot(lnewdat, aes(x = AGE_DRPED, y = Probability, colour = Level)) +
geom_point(alpha = 1/50, size=1) + facet_grid(~ Severity, labeller="label_both")+theme_bw()+
aes(colour = Level) + stat_smooth()+ labs(x="Driver Age \n (d)")+
theme(legend.position="bottom", legend.text=element_text(size=20))+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(4, 1)))
vplayout <- function(x, y)
viewport(layout.pos.row = x, layout.pos.col = y)
print(m1, vp = vplayout(1, 1))
## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.
print(m2, vp = vplayout(2, 1))
print(m3, vp = vplayout(3, 1))
print(m4, vp = vplayout(4, 1))
ex <- subset(aa1, Vis_Score1=="Excellent")
md <- subset(aa1, Vis_Score1=="Medium")
pr <- subset(aa1, Vis_Score1=="Poor")
dim(ex)
## [1] 11666 14
dim(md)
## [1] 9815 14
dim(pr)
## [1] 2138 14
names(ex)
## [1] "Vis_Score1" "HIGHESTINJ1" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "AVG_SH_WID" "AGE_DRPED" "RCIAVGTFCT" "RCIMAXSPD1" "RCIAADT1"
## [11] "SKIDNUMBER1" "LIGHTCOND1" "DIV_UNDIV1" "AGE_DRPED1"
table(ex$RCIAADT)
##
## 0 550 1100 1300 1365 1500 1600 1800 2000 2035
## 107 1 2 5 1 1 1 2 1 4
## 2100 2400 2700 2800 2900 3100 3200 3300 3400 3500
## 5 5 1 2 2 7 2 3 2 1
## 3600 3681 3700 3734 3800 3900 3950 4000 4100 4200
## 5 4 5 1 1 8 2 4 3 2
## 4300 4600 4700 4800 4900 5000 5100 5200 5300 5400
## 14 5 10 3 5 9 2 4 2 4
## 5450 5500 5600 5800 5900 6000 6080 6100 6139 6200
## 3 8 1 1 3 6 4 8 1 1
## 6229 6300 6366 6400 6500 6600 6735 6744 6800 6900
## 2 4 4 1 1 7 2 3 1 5
## 7000 7100 7200 7300 7400 7433 7500 7600 7700 7800
## 1 7 9 3 3 1 2 2 7 4
## 7900 8000 8100 8160 8200 8300 8400 8420 8500 8600
## 3 8 5 1 11 16 2 2 19 4
## 8700 8800 8900 9000 9100 9200 9300 9400 9500 9600
## 3 17 9 4 11 8 7 13 6 20
## 9632 9700 9746 9800 9900 10000 10100 10200 10300 10400
## 3 4 2 2 13 18 12 9 4 11
## 10500 10600 10700 10800 10900 11000 11100 11200 11300 11400
## 29 3 6 7 6 16 1 2 5 9
## 11500 11600 11624 11700 11800 11848 11900 12000 12100 12300
## 21 2 2 8 7 1 8 24 4 8
## 12400 12415 12500 12600 12700 12800 12900 12928 13000 13015
## 1 3 48 13 7 8 8 1 31 1
## 13100 13138 13188 13200 13300 13400 13500 13600 13700 13800
## 5 4 2 2 26 5 27 19 7 24
## 13900 14000 14100 14200 14300 14400 14500 14600 14700 14800
## 3 34 6 5 2 5 23 7 13 3
## 14900 15000 15100 15200 15300 15400 15500 15600 15700 15900
## 11 31 14 11 11 7 24 13 14 12
## 16000 16100 16200 16300 16400 16500 16600 16700 16800 16900
## 18 9 5 14 20 62 14 14 18 9
## 17000 17100 17200 17277 17300 17400 17500 17600 17700 17800
## 38 5 11 2 14 10 45 30 8 40
## 17900 18000 18100 18200 18300 18400 18500 18600 18700 18800
## 8 13 24 22 16 7 19 15 7 15
## 18900 19000 19100 19200 19300 19400 19499 19500 19600 19700
## 16 25 13 7 14 14 2 19 21 14
## 19800 19900 20000 20100 20200 20300 20400 20500 20600 20800
## 17 22 23 13 7 13 5 39 2 2
## 20900 21000 21200 21400 21500 21938 22000 22182 22200 22301
## 2 85 10 3 107 2 76 2 2 2
## 22500 22935 23000 23180 23300 23500 24000 24112 24500 24604
## 81 2 80 2 1 84 94 1 76 1
## 24867 25000 25300 25400 25500 25533 26000 26400 26500 26513
## 4 72 2 2 92 2 69 1 73 2
## 26626 26893 26997 27000 27141 27167 27500 27598 27729 27781
## 1 4 1 106 4 1 59 5 1 1
## 28000 28500 29000 29065 29398 29500 30000 30500 31000 31500
## 57 91 76 2 1 94 100 72 72 80
## 32000 32500 32728 32742 32975 33000 33012 33500 33628 33860
## 107 101 2 6 3 102 3 104 2 1
## 34000 34325 34500 35000 35500 36000 36394 36500 36564 36579
## 160 4 104 118 141 71 6 112 2 8
## 37000 37345 37500 38000 38017 38035 38102 38500 39000 39500
## 94 9 125 45 2 1 4 101 125 93
## 40000 40500 41000 41500 41792 42000 42500 43000 43300 43500
## 111 121 99 86 2 62 123 119 2 90
## 43700 44000 44500 44907 44967 45000 45500 46000 46278 46449
## 1 143 82 2 2 135 85 66 2 14
## 46470 46500 46773 47000 47500 47571 48000 48475 48500 49000
## 3 49 2 95 74 5 67 2 93 55
## 49500 50000 50065 50500 50877 51000 51300 51500 51700 51900
## 69 88 3 66 6 58 8 41 1 1
## 52000 52500 53000 53040 53500 54000 54400 54500 55000 55500
## 82 68 88 2 54 46 2 75 81 46
## 56000 56500 56966 57000 57073 57500 57596 58000 58426 58500
## 57 23 9 81 2 61 5 55 2 50
## 58590 59000 59133 59500 60000 60440 60500 61000 61500 62000
## 7 25 11 39 19 1 38 46 1 25
## 62500 63000 63500 64000 64205 64500 64973 65000 65054 65125
## 34 74 32 53 3 23 3 19 3 2
## 65500 66000 66139 66500 66952 67000 67500 68000 68500 69000
## 16 14 1 38 4 6 8 37 22 10
## 69100 69500 70000 70500 71500 72000 72500 73000 73300 73500
## 1 14 21 5 21 12 17 10 2 1
## 74000 74500 74915 75000 75500 76000 77000 77500 78000 78500
## 19 17 1 13 3 19 22 7 5 28
## 79000 79500 80000 80500 81000 81500 82000 82500 83000 83500
## 14 7 9 4 1 16 25 1 23 8
## 84000 84500 86000 86500 87000 88000 88900 89000 89500 90000
## 6 16 19 12 13 3 3 4 6 7
## 90500 91000 91500 92000 93500 94000 94500 95500 97000 98000
## 5 8 8 11 12 6 7 11 6 6
## 98500 99500 99800 100000 100500 100800 101000 102000 102500 103000
## 1 19 2 5 5 6 2 9 2 3
## 104000 104500 105000 105500 106000 106200 107000 108000 108500 109000
## 3 7 4 1 2 2 27 20 3 1
## 109500 110000 110500 110600 111000 111700 112000 112500 113000 113500
## 2 12 3 2 24 1 27 16 8 10
## 114500 114800 115779 116500 117000 117716 118470 118776 119000 119500
## 5 1 4 11 4 10 2 3 2 3
## 120000 121000 122457 122500 123000 123500 123696 123731 125000 125427
## 2 7 3 10 7 6 2 9 1 2
## 125500 126500 127500 129000 130325 130575 131000 131500 132000 132500
## 1 1 23 2 2 6 3 3 8 3
## 133500 134000 134500 134764 135000 136407 136500 137000 138000 138500
## 2 6 3 2 7 1 13 6 27 4
## 139000 139500 140000 141000 141500 142000 142500 143500 144000 147000
## 16 8 11 8 28 11 2 10 1 2
## 147500 148000 148500 149000 149500 150000 150500 151000 151500 152000
## 11 4 6 10 8 1 13 2 2 2
## 152500 153000 154500 155500 156500 157000 158000 158500 160000 161000
## 3 2 3 9 25 6 7 9 22 4
## 163000 163500 164000 164500 165000 165028 166000 166198 167000 168500
## 9 6 7 2 6 2 2 5 2 10
## 169000 170500 171000 171500 172000 172500 173000 173500 174500 175000
## 7 15 2 2 10 4 7 8 1 9
## 175500 176500 177000 177500 178500 179500 180500 181500 183000 184000
## 17 11 7 7 5 10 13 5 33 17
## 184500 185000 186000 187500 189500 190000 190500 191500 192000 192443
## 16 4 4 1 10 7 52 16 1 1
## 192500 193000 193478 193500 194000 194500 194657 195000 195500 196000
## 3 9 14 16 2 15 7 4 5 4
## 197000 197500 198000 198500 199000 200000 203000 204000 205000 205500
## 32 16 19 6 2 14 5 3 14 5
## 206000 207000 209000 210000 213000 214000 217000 219000 222000 224000
## 2 18 15 6 3 3 20 12 6 14
## 225000 227000 229000 230000 231000 235000 236000 238626 239000 239351
## 6 3 5 20 17 21 6 12 5 2
## 240000 241000 244000 245000 247000 248000 252000 253000 254000 255000
## 4 2 8 15 2 10 18 16 4 7
## 256000 258000 261000 262000 263000 264000 265000 266000 268000 270000
## 10 13 21 17 24 25 12 2 24 7
## 273000 274000 275000 277000 278000 279000 281000 287000 294000 296000
## 7 10 31 6 11 16 26 4 1 6
## 299000 301000 304000
## 15 4 5
hist(ex$RCIAADT)
library(scales)
library(ggplot2)
m1 <-ggplot(ex, aes(x=SKIDNUMBER)) + geom_histogram(colour = "#4E2A1E", fill = "#046C9A")+theme_bw()+scale_y_continuous(limits=c(0,2200))+labs(title="Excellent")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m2 <- ggplot(md, aes(x=SKIDNUMBER)) + geom_histogram(colour = "#4E2A1E", fill = "#5BBCD6")+theme_bw()+scale_y_continuous(limits=c(0,2200))+labs(title="Medium")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m3 <- ggplot(pr, aes(x=SKIDNUMBER)) + geom_histogram(colour = "#4E2A1E", fill = "#3F5151")+theme_bw()+scale_y_continuous(limits=c(0,2200))+labs(title="Poor")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
ex0 <- subset(ex, RCIAADT > 5000)
dim(ex0)
## [1] 11430 14
dim(ex)
## [1] 11666 14
md0 <- subset(md, RCIAADT > 500)
pr0 <- subset(pr, RCIAADT > 500)
ex1 <- ggplot(ex, aes(x=RCIAADT))
md1 <- ggplot(md, aes(x=RCIAADT))
pr1 <- ggplot(pr, aes(x=RCIAADT))
m01 <- ex1+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A", breaks=c(0, 50000, 100000, 150000, 200000, 250000, 300000))+ scale_y_continuous(labels=percent, limits=c(0,0.7))+
labs(x=" \n (d)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m02 <- md1+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6", breaks=c(0, 50000, 100000, 150000, 200000, 250000, 300000))+ scale_y_continuous(labels=percent, limits=c(0,0.7))+
labs(x="AADT (vpd) \n (e)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m03 <- pr1+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312", breaks=c(0, 50000, 100000, 150000, 200000, 250000, 300000))+ scale_y_continuous(labels=percent, limits=c(0,0.7))+
labs(x=" \n (f)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
ex1 <- ggplot(ex, aes(x=RCIMAXSPD))
md1 <- ggplot(md, aes(x=RCIMAXSPD))
pr1 <- ggplot(pr, aes(x=RCIMAXSPD))
m1 <- ex1+ geom_bar(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A")+ scale_y_continuous(labels=percent, limits=c(0,0.15))+
labs(x=" \n (a)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(title="Excellent")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m2 <- md1+ geom_bar(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6")+ scale_y_continuous(labels=percent, limits=c(0,0.15))+
labs(x="Max. Speed (mph) \n (b)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+labs(title="Medium")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m3 <- pr1+ geom_bar(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312")+ scale_y_continuous(labels=percent, limits=c(0,0.15))+
labs(x=" \n (c)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+labs(title="Poor")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
exx <- subset(ex, SKIDNUMBER >20)
mdd <- subset(md, SKIDNUMBER >20)
prr <- subset(pr, SKIDNUMBER >20)
ex3 <- ggplot(exx, aes(x=SKIDNUMBER))
md3 <- ggplot(mdd, aes(x=SKIDNUMBER))
pr3 <- ggplot(prr, aes(x=SKIDNUMBER))
m4 <- ex3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A")+ scale_y_continuous(labels=percent, limits=c(0,0.20))+
labs(x=" \n (g)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m5 <- md3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6")+ scale_y_continuous(labels=percent, limits=c(0,0.20))+
labs(x="Skid Number \n (h)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m6 <- pr3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312")+ scale_y_continuous(labels=percent, limits=c(0,0.20))+
labs(x=" \n (i)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
ex3 <- ggplot(ex, aes(x=AGE_DRPED))
md3 <- ggplot(md, aes(x=AGE_DRPED))
pr3 <- ggplot(pr, aes(x=AGE_DRPED))
m4a <- ex3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (j)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m5a <- md3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x="Driver Age \n (k)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m6a <- pr3+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (l)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
ex4 <- ggplot(ex, aes(x=RCIAVGTFCT))
md4 <- ggplot(md, aes(x=RCIAVGTFCT))
pr4 <- ggplot(pr, aes(x=RCIAVGTFCT))
m4b <- ex4+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (m)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m5b <- md4+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x="Percentage of Trucks \n (n)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m6b <- pr4+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (o)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
ex4t <- ggplot(ex, aes(x=AVG_SH_WID))
md4t <- ggplot(md, aes(x=AVG_SH_WID))
pr4t <- ggplot(pr, aes(x=AVG_SH_WID))
m4c <- ex4t+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#046C9A")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (p)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m5c <- md4t+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#5BBCD6")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x="Avg. Shoulder Width (ft.) \n (q)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
m6c <- pr4t+ geom_histogram(aes(y = (..count..)/sum(..count..)), colour = "#4E2A1E", fill = "#C93312")+ scale_y_continuous(labels=percent, limits=c(0,0.125))+
labs(x=" \n (r)", y="Percentage") +
theme(text = element_text(size=30)) + theme_bw()+labs(y=" ")+theme(text = element_text(size=16), axis.text = element_text(size = 16), strip.text = element_text(size = 17, colour = "#550307"))
library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(6, 3)))
vplayout <- function(x, y)
viewport(layout.pos.row = x, layout.pos.col = y)
print(m1, vp = vplayout(1, 1))
print(m2, vp = vplayout(1, 2))
print(m3, vp = vplayout(1, 3))
print(m01, vp = vplayout(2, 1))
print(m02, vp = vplayout(2, 2))
print(m03, vp = vplayout(2, 3))
print(m4, vp = vplayout(3, 1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m5, vp = vplayout(3, 2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m6, vp = vplayout(3, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m4a, vp = vplayout(4, 1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m5a, vp = vplayout(4, 2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m6a, vp = vplayout(4, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m4b, vp = vplayout(5, 1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m5b, vp = vplayout(5, 2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m6b, vp = vplayout(5, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m4c, vp = vplayout(6, 1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m5c, vp = vplayout(6, 2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
print(m6c, vp = vplayout(6, 3))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
