Reading the Data
setwd("C:/Users/s-das/Syncplicity Folders/SHRP2 Visibility Study (TOG-CTS) (Brad Brimley)/FL_Buffer_5Miles/Test Set removing Holidays+Extreme Weather")
##fl_all <- read.csv("FL_Main_Control1.csv")
##dim(fl_all)
##names(fl_all)
##fl_all_01 <- fl_all[c(6, 67, 136, 123,120, 112,116,117,119,254, 280,155,129, 66,126 , 256, 257, 64, 65, 125, 106,
##56, 46, 47, 75, 281)]
##names(fl_all_01)
##a2 <- fl_all_01[c(11, 14, 15, 2, 4, 5, 6, 13, 21, 26, 12, 8, 9, 1, 3, 10, 15, 18:25, 16, 17)]
##names(a2)
## write.csv(a2, "FL_Selec.csv")
a03 <- read.csv("FL_Selec1.csv")
dim(a03)
## [1] 43274 25
a04 <- a03[!duplicated(a03[,1]),]
dim(a04)
## [1] 13156 25
names(a04)
## [1] "CRASHNUM1" "DataType" "Visibility" "Severity"
## [5] "Vis_Score" "RCIMAXSPD" "RCIAADT" "RCISURFWTH"
## [9] "LIGHTCOND" "SKIDNUMBER" "MatchID" "Category"
## [13] "Category1" "HIGHESTINJ" "ID" "AP_LATITUDE"
## [17] "AP_LONGITUDE" "CR_LATITUDE" "CR_LONGITUDE" "WORKZONE"
## [21] "RCISLDWTH2" "RCIMEDWDTH" "FULLNAME" "RDSURFTYPE"
## [25] "AGE_DRPED"
a05 <- a04[c(2:10)]
summary(a05)
## DataType Visibility Severity Vis_Score
## Control:6654 Excellent:6654 Fatal : 110 Min. : 0.500
## Main :6502 Medium :5328 Injury :6032 1st Qu.: 2.250
## Poor :1174 No Injury:7014 Median :10.000
## Mean : 6.148
## 3rd Qu.:10.000
## Max. :10.000
## RCIMAXSPD RCIAADT RCISURFWTH LIGHTCOND
## Min. :25.00 Min. : 500 Min. : 9.00 Min. : 1.000
## 1st Qu.:40.00 1st Qu.: 27500 1st Qu.:24.00 1st Qu.: 1.000
## Median :45.00 Median : 41500 Median :33.00 Median : 1.000
## Mean :46.08 Mean : 61907 Mean :32.24 Mean : 2.467
## 3rd Qu.:50.00 3rd Qu.: 62000 3rd Qu.:36.00 3rd Qu.: 4.000
## Max. :70.00 Max. :304000 Max. :84.00 Max. :88.000
## SKIDNUMBER
## Min. : 0.00
## 1st Qu.:33.00
## Median :36.00
## Mean :36.55
## 3rd Qu.:40.00
## Max. :62.00
ftable(DataType~Visibility+Severity , data=a05)
## DataType Control Main
## Visibility Severity
## Excellent Fatal 65 0
## Injury 3027 0
## No Injury 3562 0
## Medium Fatal 0 33
## Injury 0 2428
## No Injury 0 2867
## Poor Fatal 0 12
## Injury 0 577
## No Injury 0 585
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
aa <- a05 %>%
group_by(Visibility, Severity) %>%
summarise (n = n()) %>%
mutate(freq = n / sum(n))
aa
## Source: local data frame [9 x 4]
## Groups: Visibility [3]
##
## Visibility Severity n freq
## (fctr) (fctr) (int) (dbl)
## 1 Excellent Fatal 65 0.009768560
## 2 Excellent Injury 3027 0.454914337
## 3 Excellent No Injury 3562 0.535317102
## 4 Medium Fatal 33 0.006193694
## 5 Medium Injury 2428 0.455705706
## 6 Medium No Injury 2867 0.538100601
## 7 Poor Fatal 12 0.010221465
## 8 Poor Injury 577 0.491482112
## 9 Poor No Injury 585 0.498296422
names(a05)
## [1] "DataType" "Visibility" "Severity" "Vis_Score" "RCIMAXSPD"
## [6] "RCIAADT" "RCISURFWTH" "LIGHTCOND" "SKIDNUMBER"
a06 <- a05[c(2, 3, 5:6, 9, 7)]
names(a06)
## [1] "Visibility" "Severity" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "RCISURFWTH"
Response Variable= Visibility
### MULTINOM
library(nnet)
test <- multinom(Visibility ~Severity+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05)
## # weights: 24 (14 variable)
## initial value 14453.343270
## iter 10 value 12538.175524
## iter 20 value 12160.243304
## final value 12160.242548
## converged
summary(test)
## Call:
## multinom(formula = Visibility ~ Severity + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, data = a05)
##
## Coefficients:
## (Intercept) SeverityInjury SeverityNo Injury RCIMAXSPD
## Medium -0.9871889 0.48917664 0.4983746 0.012534852
## Poor -2.3374533 0.03366072 -0.1039777 0.006956064
## RCIAADT SKIDNUMBER RCISURFWTH
## Medium -1.672057e-07 -0.005188727 -0.003221014
## Poor -1.771266e-06 0.005486080 0.007046785
##
## Std. Errors:
## (Intercept) SeverityInjury SeverityNo Injury RCIMAXSPD
## Medium 7.132128e-06 3.304684e-06 3.775045e-06 0.0002943583
## Poor 1.212562e-05 6.091315e-06 5.887973e-06 0.0004992284
## RCIAADT SKIDNUMBER RCISURFWTH
## Medium 3.406609e-07 0.0002742943 0.0001726231
## Poor 6.037767e-07 0.0004717143 0.0002990539
##
## Residual Deviance: 24320.49
## AIC: 24348.49
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
m <- polr(Visibility ~Severity+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH,
data = a05, Hess=TRUE)
newdat <- cbind(a06, predict(m, a06[c(2:6)], type = "probs"))
head(newdat)
## Visibility Severity RCIMAXSPD RCIAADT SKIDNUMBER RCISURFWTH Excellent
## 1 Excellent No Injury 30 11900 44 22 0.5443798
## 3 Excellent No Injury 35 34000 36 19 0.5325164
## 5 Excellent Injury 30 21000 44 44 0.5359791
## 9 Excellent Injury 45 22500 33 24 0.4958473
## 11 Excellent No Injury 45 22500 33 24 0.5047310
## 15 Excellent Injury 45 19100 45 36 0.4997306
## Medium Poor
## 1 0.3784140 0.07720612
## 3 0.3868070 0.08067665
## 5 0.3843704 0.07965052
## 9 0.4118911 0.09226161
## 11 0.4059407 0.08932830
## 15 0.4093005 0.09096888
library(reshape2)
lnewdat <- melt(newdat, id.vars = c("Visibility", "Severity", "RCIMAXSPD", "RCIAADT", "SKIDNUMBER", "RCISURFWTH"),
variable.name = "Level", value.name="Probability")
head(lnewdat)
## Visibility Severity RCIMAXSPD RCIAADT SKIDNUMBER RCISURFWTH Level
## 1 Excellent No Injury 30 11900 44 22 Excellent
## 2 Excellent No Injury 35 34000 36 19 Excellent
## 3 Excellent Injury 30 21000 44 44 Excellent
## 4 Excellent Injury 45 22500 33 24 Excellent
## 5 Excellent No Injury 45 22500 33 24 Excellent
## 6 Excellent Injury 45 19100 45 36 Excellent
## Probability
## 1 0.5443798
## 2 0.5325164
## 3 0.5359791
## 4 0.4958473
## 5 0.5047310
## 6 0.4997306
library(ggplot2)
ggplot(lnewdat, aes(x = RCIMAXSPD, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = RCIAADT, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = RCISURFWTH, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = SKIDNUMBER, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
Response Variable= Severity
### MULTINOM
library(nnet)
test <- multinom(Severity ~Visibility+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05)
## # weights: 24 (14 variable)
## initial value 14453.343270
## iter 10 value 10397.269273
## iter 20 value 9590.032134
## final value 9588.350048
## converged
summary(test)
## Call:
## multinom(formula = Severity ~ Visibility + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, data = a05)
##
## Coefficients:
## (Intercept) VisibilityMedium VisibilityPoor RCIMAXSPD
## Injury 2.896103 0.4770137 0.03917657 -0.02672621
## No Injury 4.268856 0.4861324 -0.09860219 -0.04115792
## RCIAADT SKIDNUMBER RCISURFWTH
## Injury 7.692452e-06 0.05171408 -0.003089616
## No Injury 1.048614e-05 0.03958956 -0.011693828
##
## Std. Errors:
## (Intercept) VisibilityMedium VisibilityPoor RCIMAXSPD
## Injury 3.39709e-06 1.345086e-06 3.169906e-07 0.0001400554
## No Injury 3.40056e-06 1.346232e-06 3.169006e-07 0.0001401717
## RCIAADT SKIDNUMBER RCISURFWTH
## Injury 1.406955e-06 0.0001312772 8.270150e-05
## No Injury 1.405715e-06 0.0001313758 8.277304e-05
##
## Residual Deviance: 19176.7
## AIC: 19204.7
library(MASS)
m <- polr(Severity ~Visibility+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05, Hess=TRUE)
m
## Call:
## polr(formula = Severity ~ Visibility + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, data = a05, Hess = TRUE)
##
## Coefficients:
## VisibilityMedium VisibilityPoor RCIMAXSPD RCIAADT
## 2.411532e-02 -1.361726e-01 -1.543723e-02 3.027336e-06
## SKIDNUMBER RCISURFWTH
## -9.664274e-03 -8.566029e-03
##
## Intercepts:
## Fatal|Injury Injury|No Injury
## -5.943349 -1.289152
##
## Residual Deviance: 19207.10
## AIC: 19223.10
a06 <- a05[c(3, 2, 5:6, 9, 7)]
names(a06)
## [1] "Severity" "Visibility" "RCIMAXSPD" "RCIAADT" "SKIDNUMBER"
## [6] "RCISURFWTH"
newdat <- cbind(a06, predict(m, a06[c(2:6)], type = "probs"))
head(newdat)
## Severity Visibility RCIMAXSPD RCIAADT SKIDNUMBER RCISURFWTH
## 1 No Injury Excellent 30 11900 44 22
## 3 No Injury Excellent 35 34000 36 19
## 5 Injury Excellent 30 21000 44 44
## 9 Injury Excellent 45 22500 33 24
## 11 No Injury Excellent 45 22500 33 24
## 15 Injury Excellent 45 19100 45 36
## Fatal Injury No Injury
## 1 0.007372659 0.4308496 0.5617777
## 3 0.006724122 0.4088157 0.5844602
## 5 0.008648588 0.4694963 0.5218551
## 9 0.008225284 0.4573079 0.5344669
## 11 0.008225284 0.4573079 0.5344669
## 15 0.010320728 0.5124061 0.4772731
library(reshape2)
lnewdat <- melt(newdat, id.vars = c("Severity", "Visibility", "RCIMAXSPD", "RCIAADT", "SKIDNUMBER", "RCISURFWTH"),
variable.name = "Level", value.name="Probability")
head(lnewdat)
## Severity Visibility RCIMAXSPD RCIAADT SKIDNUMBER RCISURFWTH Level
## 1 No Injury Excellent 30 11900 44 22 Fatal
## 2 No Injury Excellent 35 34000 36 19 Fatal
## 3 Injury Excellent 30 21000 44 44 Fatal
## 4 Injury Excellent 45 22500 33 24 Fatal
## 5 No Injury Excellent 45 22500 33 24 Fatal
## 6 Injury Excellent 45 19100 45 36 Fatal
## Probability
## 1 0.007372659
## 2 0.006724122
## 3 0.008648588
## 4 0.008225284
## 5 0.008225284
## 6 0.010320728
library(ggplot2)
ggplot(lnewdat, aes(x = RCIMAXSPD, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Severity~ Visibility, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = RCIAADT, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Severity~ Visibility, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = RCISURFWTH, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
ggplot(lnewdat, aes(x = SKIDNUMBER, y = Probability, colour = Level)) +
geom_point(alpha = 1/10, size=1) + facet_grid(Visibility ~ Severity, labeller="label_both")+
aes(colour = Level) + stat_summary(fun.y = mean, geom="line")+theme_bw()
######## VGAM
library(VGAM)
## Warning: package 'VGAM' was built under R version 3.2.5
## Loading required package: stats4
## Loading required package: splines
mod2 <- vglm(Severity ~Visibility+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05,
family = multinomial(refLevel = "Fatal"))
summary(mod2)
##
## Call:
## vglm(formula = Severity ~ Visibility + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, family = multinomial(refLevel = "Fatal"),
## data = a05)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,2]/mu[,1]) -15.06 -0.6006 -0.545 0.8228 1.303
## log(mu[,3]/mu[,1]) -15.19 -0.6929 0.635 0.7133 0.985
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 2.896e+00 7.212e-01 4.016 5.93e-05 ***
## (Intercept):2 4.269e+00 7.185e-01 5.941 2.83e-09 ***
## VisibilityMedium:1 4.773e-01 2.161e-01 2.209 0.027193 *
## VisibilityMedium:2 4.864e-01 2.159e-01 2.253 0.024229 *
## VisibilityPoor:1 4.019e-02 3.183e-01 0.126 0.899517
## VisibilityPoor:2 -9.759e-02 3.182e-01 -0.307 0.759107
## RCIMAXSPD:1 -2.674e-02 1.101e-02 -2.429 0.015161 *
## RCIMAXSPD:2 -4.118e-02 1.101e-02 -3.740 0.000184 ***
## RCIAADT:1 7.692e-06 3.239e-06 2.375 0.017561 *
## RCIAADT:2 1.049e-05 3.237e-06 3.240 0.001196 **
## SKIDNUMBER:1 5.172e-02 1.108e-02 4.667 3.06e-06 ***
## SKIDNUMBER:2 3.960e-02 1.100e-02 3.601 0.000317 ***
## RCISURFWTH:1 -3.072e-03 1.368e-02 -0.225 0.822330
## RCISURFWTH:2 -1.168e-02 1.367e-02 -0.854 0.392953
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
##
## Dispersion Parameter for multinomial family: 1
##
## Residual deviance: 19176.7 on 26298 degrees of freedom
##
## Log-likelihood: -9588.35 on 26298 degrees of freedom
##
## Number of iterations: 8
##
## Reference group is level 1 of the response
exp(VGAM::coef(mod2))
## (Intercept):1 (Intercept):2 VisibilityMedium:1
## 18.1023313 71.4334729 1.6117614
## VisibilityMedium:2 VisibilityPoor:1 VisibilityPoor:2
## 1.6265306 1.0410134 0.9070251
## RCIMAXSPD:1 RCIMAXSPD:2 RCIAADT:1
## 0.9736101 0.9596603 1.0000077
## RCIAADT:2 SKIDNUMBER:1 SKIDNUMBER:2
## 1.0000105 1.0530820 1.0403915
## RCISURFWTH:1 RCISURFWTH:2
## 0.9969332 0.9883923
mod3 <- vglm(Visibility ~Severity+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05,
family = multinomial(refLevel = "Excellent"))
summary(mod3)
##
## Call:
## vglm(formula = Visibility ~ Severity + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, family = multinomial(refLevel = "Excellent"),
## data = a05)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,2]/mu[,1]) -1.121 -0.8666 -0.7705 1.1958 1.769
## log(mu[,3]/mu[,1]) -0.548 -0.4607 -0.4166 -0.1121 3.780
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -9.872e-01 2.693e-01 -3.666 0.000246 ***
## (Intercept):2 -2.338e+00 4.272e-01 -5.472 4.44e-08 ***
## SeverityInjury:1 4.891e-01 2.161e-01 2.263 0.023611 *
## SeverityInjury:2 3.387e-02 3.180e-01 0.107 0.915155
## SeverityNo Injury:1 4.983e-01 2.159e-01 2.309 0.020967 *
## SeverityNo Injury:2 -1.038e-01 3.179e-01 -0.327 0.744014
## RCIMAXSPD:1 1.254e-02 2.407e-03 5.208 1.91e-07 ***
## RCIMAXSPD:2 6.955e-03 4.176e-03 1.665 0.095865 .
## RCIAADT:1 -1.673e-07 5.422e-07 -0.308 0.757702
## RCIAADT:2 -1.770e-06 9.482e-07 -1.867 0.061883 .
## SKIDNUMBER:1 -5.188e-03 3.146e-03 -1.649 0.099110 .
## SKIDNUMBER:2 5.485e-03 5.556e-03 0.987 0.323469
## RCISURFWTH:1 -3.220e-03 2.592e-03 -1.242 0.214159
## RCISURFWTH:2 7.045e-03 4.406e-03 1.599 0.109821
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
##
## Dispersion Parameter for multinomial family: 1
##
## Residual deviance: 24320.49 on 26298 degrees of freedom
##
## Log-likelihood: -12160.24 on 26298 degrees of freedom
##
## Number of iterations: 5
##
## Reference group is level 1 of the response
exp(VGAM::coef(mod3))
## (Intercept):1 (Intercept):2 SeverityInjury:1
## 0.37261658 0.09656466 1.63092452
## SeverityInjury:2 SeverityNo Injury:1 SeverityNo Injury:2
## 1.03445483 1.64598480 0.90140100
## RCIMAXSPD:1 RCIMAXSPD:2 RCIAADT:1
## 1.01261417 1.00697878 0.99999983
## RCIAADT:2 SKIDNUMBER:1 SKIDNUMBER:2
## 0.99999823 0.99482508 1.00550040
## RCISURFWTH:1 RCISURFWTH:2
## 0.99678481 1.00706952
mod4 <- vglm(Visibility ~Severity+ RCIMAXSPD+ RCIAADT+ SKIDNUMBER+ RCISURFWTH, data = a05,
family = multinomial(refLevel = "Poor"))
summary(mod4)
##
## Call:
## vglm(formula = Visibility ~ Severity + RCIMAXSPD + RCIAADT +
## SKIDNUMBER + RCISURFWTH, family = multinomial(refLevel = "Poor"),
## data = a05)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,1]/mu[,3]) -2.711 -0.5879 0.7928 0.8662 1.074
## log(mu[,2]/mu[,3]) -2.719 -0.4711 -0.4404 1.0518 1.603
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 2.338e+00 4.272e-01 5.472 4.44e-08 ***
## (Intercept):2 1.350e+00 4.472e-01 3.020 0.00253 **
## SeverityInjury:1 -3.387e-02 3.180e-01 -0.107 0.91515
## SeverityInjury:2 4.553e-01 3.410e-01 1.335 0.18180
## SeverityNo Injury:1 1.038e-01 3.179e-01 0.327 0.74401
## SeverityNo Injury:2 6.021e-01 3.409e-01 1.766 0.07733 .
## RCIMAXSPD:1 -6.955e-03 4.176e-03 -1.665 0.09587 .
## RCIMAXSPD:2 5.581e-03 4.220e-03 1.322 0.18601
## RCIAADT:1 1.770e-06 9.482e-07 1.867 0.06188 .
## RCIAADT:2 1.603e-06 9.625e-07 1.666 0.09579 .
## SKIDNUMBER:1 -5.485e-03 5.556e-03 -0.987 0.32347
## SKIDNUMBER:2 -1.067e-02 5.646e-03 -1.890 0.05870 .
## RCISURFWTH:1 -7.045e-03 4.406e-03 -1.599 0.10982
## RCISURFWTH:2 -1.027e-02 4.495e-03 -2.283 0.02241 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
##
## Dispersion Parameter for multinomial family: 1
##
## Residual deviance: 24320.49 on 26298 degrees of freedom
##
## Log-likelihood: -12160.24 on 26298 degrees of freedom
##
## Number of iterations: 5
##
## Reference group is level 3 of the response
exp(VGAM::coef(mod4))
## (Intercept):1 (Intercept):2 SeverityInjury:1
## 10.3557549 3.8587260 0.9666928
## SeverityInjury:2 SeverityNo Injury:1 SeverityNo Injury:2
## 1.5766029 1.1093842 1.8260295
## RCIMAXSPD:1 RCIMAXSPD:2 RCIAADT:1
## 0.9930696 1.0055963 1.0000018
## RCIAADT:2 SKIDNUMBER:1 SKIDNUMBER:2
## 1.0000016 0.9945297 0.9893831
## RCISURFWTH:1 RCISURFWTH:2
## 0.9929801 0.9897875