setwd("/Users/subasishdas1/Dropbox/---- TRB_2017 ----/SHRP2 Inclement")
aa <- read.csv("FINAL_5a.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
## Excellent Visibility:11666 Fatal : 181 Min. :25.00
## Medium Visibility : 9815 Injury :11293 1st Qu.:40.00
## Poor Visibility : 2138 No Injury:12145 Median :45.00
## Mean :45.25
## 3rd Qu.:50.00
## Max. :70.00
##
## RCIAADT SKIDNUMBER AVG_SH_WID AGE_DRPED
## Min. : 0 Min. : 0.00 Min. : 0.000 Min. : 15.00
## 1st Qu.: 26000 1st Qu.:33.00 1st Qu.: 1.000 1st Qu.: 25.00
## Median : 40500 Median :36.00 Median : 3.000 Median : 38.00
## Mean : 60909 Mean :35.29 Mean : 3.844 Mean : 39.57
## 3rd Qu.: 60500 3rd Qu.:40.00 3rd Qu.: 6.000 3rd Qu.: 51.00
## Max. :304000 Max. :62.00 Max. :19.500 Max. :108.00
##
## RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1
## Min. : 0.000 >60 : 2600 >124,999 :3131 > 50 : 250
## 1st Qu.: 3.300 0-30 : 1827 0-9999 : 871 20-30: 1738
## Median : 4.530 30-45:15013 10,000-34,999 :8127 30-40:15920
## Mean : 5.171 45-60: 4179 35,000-54,999 :7226 40-50: 4702
## 3rd Qu.: 6.260 55,000-124,999:4050 NA's : 1009
## Max. :31.770 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 Visibility|Medium Visibility 0.1235 0.1710 0.722
## Medium Visibility|Poor Visibility 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 Visibility|Medium Visibility 0.1277 0.1707 0.748
## Medium Visibility|Poor Visibility 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 Visibility|Medium Visibility
## 0.1235173
## Medium Visibility|Poor Visibility
## 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
## 1 Excellent Visibility No Injury 35 34000 36 4.5
## 2 Excellent Visibility Injury 30 21000 44 1.0
## 3 Excellent Visibility Injury 30 21000 44 1.0
## 4 Excellent Visibility Injury 45 22500 33 7.5
## 5 Excellent Visibility No Injury 45 22500 33 6.0
## 6 Excellent Visibility No Injury 45 22500 33 6.0
## AGE_DRPED RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1
## 1 31 8.35 30-45 10,000-34,999 30-40
## 2 55 8.35 0-30 10,000-34,999 40-50
## 3 47 8.35 0-30 10,000-34,999 40-50
## 4 44 8.33 30-45 10,000-34,999 30-40
## 5 32 8.33 30-45 10,000-34,999 30-40
## 6 40 8.33 30-45 10,000-34,999 30-40
## LIGHTCOND1 DIV_UNDIV1 AGE_DRPED1 Excellent Visibility
## 1 Dark(No Street light) undivided 30-39 0.4976801
## 2 Dark(Street Light) undivided 50-59 0.5535496
## 3 Dark(Street Light) undivided 40-49 0.5424877
## 4 Daylight divided 40-49 0.4862752
## 5 Daylight divided 30-39 0.4711119
## 6 Daylight divided 40-49 0.4822525
## Medium Visibility Poor Visibility
## 1 0.4135601 0.08875980
## 2 0.3742370 0.07221346
## 3 0.3822488 0.07526348
## 4 0.4212042 0.09252061
## 5 0.4311374 0.09775064
## 6 0.4238656 0.09388186
names(newdat)
## [1] "Vis_Score1" "HIGHESTINJ1" "RCIMAXSPD"
## [4] "RCIAADT" "SKIDNUMBER" "AVG_SH_WID"
## [7] "AGE_DRPED" "RCIAVGTFCT" "RCIMAXSPD1"
## [10] "RCIAADT1" "SKIDNUMBER1" "LIGHTCOND1"
## [13] "DIV_UNDIV1" "AGE_DRPED1" "Excellent Visibility"
## [16] "Medium Visibility" "Poor Visibility"
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
## 1 Excellent Visibility No Injury 35 34000 36 4.5
## 2 Excellent Visibility Injury 30 21000 44 1.0
## 3 Excellent Visibility Injury 30 21000 44 1.0
## 4 Excellent Visibility Injury 45 22500 33 7.5
## 5 Excellent Visibility No Injury 45 22500 33 6.0
## 6 Excellent Visibility No Injury 45 22500 33 6.0
## AGE_DRPED RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1
## 1 31 8.35 30-45 10,000-34,999 30-40
## 2 55 8.35 0-30 10,000-34,999 40-50
## 3 47 8.35 0-30 10,000-34,999 40-50
## 4 44 8.33 30-45 10,000-34,999 30-40
## 5 32 8.33 30-45 10,000-34,999 30-40
## 6 40 8.33 30-45 10,000-34,999 30-40
## LIGHTCOND1 DIV_UNDIV1 AGE_DRPED1 Level
## 1 Dark(No Street light) undivided 30-39 Excellent Visibility
## 2 Dark(Street Light) undivided 50-59 Excellent Visibility
## 3 Dark(Street Light) undivided 40-49 Excellent Visibility
## 4 Daylight divided 40-49 Excellent Visibility
## 5 Daylight divided 30-39 Excellent Visibility
## 6 Daylight divided 40-49 Excellent Visibility
## Probability
## 1 0.4976801
## 2 0.5535496
## 3 0.5424877
## 4 0.4862752
## 5 0.4711119
## 6 0.4822525
names(lnewdat)[2] <- "Severity"
head(lnewdat)
## Vis_Score1 Severity RCIMAXSPD RCIAADT SKIDNUMBER AVG_SH_WID
## 1 Excellent Visibility No Injury 35 34000 36 4.5
## 2 Excellent Visibility Injury 30 21000 44 1.0
## 3 Excellent Visibility Injury 30 21000 44 1.0
## 4 Excellent Visibility Injury 45 22500 33 7.5
## 5 Excellent Visibility No Injury 45 22500 33 6.0
## 6 Excellent Visibility No Injury 45 22500 33 6.0
## AGE_DRPED RCIAVGTFCT RCIMAXSPD1 RCIAADT1 SKIDNUMBER1
## 1 31 8.35 30-45 10,000-34,999 30-40
## 2 55 8.35 0-30 10,000-34,999 40-50
## 3 47 8.35 0-30 10,000-34,999 40-50
## 4 44 8.33 30-45 10,000-34,999 30-40
## 5 32 8.33 30-45 10,000-34,999 30-40
## 6 40 8.33 30-45 10,000-34,999 30-40
## LIGHTCOND1 DIV_UNDIV1 AGE_DRPED1 Level
## 1 Dark(No Street light) undivided 30-39 Excellent Visibility
## 2 Dark(Street Light) undivided 50-59 Excellent Visibility
## 3 Dark(Street Light) undivided 40-49 Excellent Visibility
## 4 Daylight divided 40-49 Excellent Visibility
## 5 Daylight divided 30-39 Excellent Visibility
## 6 Daylight divided 40-49 Excellent Visibility
## Probability
## 1 0.4976801
## 2 0.5535496
## 3 0.5424877
## 4 0.4862752
## 5 0.4711119
## 6 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))
setwd("/Users/subasishdas1/Dropbox/---- TRB_2017 ----/SHRP2 Inclement")
aa <- read.csv("FINAL_5a.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"
ex <- subset(aa1, Vis_Score1=="Excellent Visibility")
md <- subset(aa1, Vis_Score1=="Medium Visibility")
pr <- subset(aa1, Vis_Score1=="Poor Visibility")
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$Vis_Score1)
##
## Excellent Visibility Medium Visibility Poor Visibility
## 11666 0 0
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 Visibility")+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 Visibility")+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 Visibility")+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 Visibility")+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 Visibility")+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 Visibility")+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`.
