setwd("C:/Users/koho0/Desktop")
require(Epi)
## Loading required package: Epi
require(pROC)
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
require(ztable)
## Loading required package: ztable
## Welcome to package ztable ver 0.2.0
require(moonBook)
## Loading required package: moonBook
source("ROC_sub.R")
crps <-read.csv("crps.csv")
crps$result <-as.factor(crps$result)
crps$body_t <-as.factor(crps$body_t)
crps$sudomotor <-as.factor(crps$sudomotor)
crps$wbbs <-as.factor(crps$wbbs)
str(crps)
## 'data.frame': 32 obs. of 5 variables:
## $ result : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp : num 1.29 0.98 0.77 1.58 0.8 2.17 0.63 0.67 3.5 0.3 ...
## $ body_t : Factor w/ 2 levels "1","2": 1 2 2 1 2 1 2 2 1 2 ...
## $ sudomotor: Factor w/ 2 levels "1","2": 2 2 1 1 2 2 2 2 2 2 ...
## $ wbbs : Factor w/ 2 levels "1","2": 2 2 2 2 1 2 2 2 2 2 ...
head(crps)
## result temp body_t sudomotor wbbs
## 1 0 1.29 1 2 2
## 2 0 0.98 2 2 2
## 3 0 0.77 2 1 2
## 4 0 1.58 1 1 2
## 5 0 0.80 2 2 1
## 6 0 2.17 1 2 2
tsw <- glm(result ~ temp + sudomotor + wbbs, data=crps, family="binomial")
fitted(tsw)
## 1 2 3 4 5 6 7
## 0.13865573 0.10431215 0.44079797 0.64745775 0.68904508 0.28749152 0.07476648
## 8 9 10 11 12 13 14
## 0.07770768 0.61803030 0.05415315 0.96921101 0.99245943 0.96559234 0.68459202
## 15 16 17 18 19 20 21
## 0.99214033 0.15430921 0.22232186 0.96524373 0.59845742 0.80404033 0.96761425
## 22 23 24 25 26 27 28
## 0.96661807 0.97865245 0.99704021 0.97971639 0.81213621 0.89694676 0.99146140
## 29 30 31 32
## 0.99930515 0.96761425 0.96303029 0.99907905
crps$pred_tsw=fitted(tsw)
a1=ROC(form=result~pred_tsw,data=crps,plot="ROC")
a1
## $res
## sens spec pvp pvn lr.eta
## 1.00000000 0.0 NaN 0.31250000 -Inf
## 0.0821257466746746 1.00000000 0.1 0.0000000 0.29032258 0.08212575
## 0.0914494613541141 1.00000000 0.2 0.0000000 0.26666667 0.09144946
## 0.09285495190196 1.00000000 0.3 0.0000000 0.24137931 0.09285495
## 0.106471852972547 1.00000000 0.4 0.0000000 0.21428571 0.10647185
## 0.126626195193206 1.00000000 0.5 0.0000000 0.18518519 0.12662620
## 0.136849117383786 0.95454545 0.5 0.1666667 0.19230769 0.13684912
## 0.189506873389417 0.90909091 0.5 0.2857143 0.20000000 0.18950687
## 0.253323637040478 0.90909091 0.6 0.2500000 0.16666667 0.25332364
## 0.448862673376277 0.90909091 0.7 0.2222222 0.13043478 0.44886267
## 0.66714370038935 0.86363636 0.7 0.3000000 0.13636364 0.66714370
## 0.691490391795274 0.86363636 0.8 0.2727273 0.09523810 0.69149039
## 0.726152347258012 0.86363636 0.9 0.2500000 0.05000000 0.72615235
## 0.766258077475893 0.81818182 0.9 0.3076923 0.05263158 0.76625808
## 0.77078297256912 0.81818182 1.0 0.2857143 0.00000000 0.77078297
## 0.866411799100658 0.77272727 1.0 0.3333333 0.00000000 0.86641180
## 0.871674110357614 0.72727273 1.0 0.3750000 0.00000000 0.87167411
## 0.916849130687502 0.68181818 1.0 0.4117647 0.00000000 0.91684913
## 0.941462115977832 0.63636364 1.0 0.4444444 0.00000000 0.94146212
## 0.942155025369667 0.59090909 1.0 0.4736842 0.00000000 0.94215503
## 0.942263454752985 0.54545455 1.0 0.5000000 0.00000000 0.94226345
## 0.942581379577532 0.50000000 1.0 0.5238095 0.00000000 0.94258138
## 0.942888570363306 0.40909091 1.0 0.5652174 0.00000000 0.94288857
## 0.943377744759295 0.36363636 1.0 0.5833333 0.00000000 0.94337774
## 0.946190546741504 0.31818182 1.0 0.6000000 0.00000000 0.94619055
## 0.946499128323634 0.27272727 1.0 0.6153846 0.00000000 0.94649913
## 0.94979637708728 0.22727273 1.0 0.6296296 0.00000000 0.94979638
## 0.949980973831789 0.18181818 1.0 0.6428571 0.00000000 0.94998097
## 0.950067512964902 0.13636364 1.0 0.6551724 0.00000000 0.95006751
## 0.951294279635132 0.09090909 1.0 0.6666667 0.00000000 0.95129428
## 0.951831051244294 0.04545455 1.0 0.6774194 0.00000000 0.95183105
## 0.951890230079601 0.00000000 1.0 0.6875000 NaN 0.95189023
##
## $AUC
## [1] 0.9363636
##
## $lr
##
## Call: glm(formula = form, family = binomial, data = data)
##
## Coefficients:
## (Intercept) pred_tsw
## -2.723 5.712
##
## Degrees of Freedom: 31 Total (i.e. Null); 30 Residual
## Null Deviance: 39.75
## Residual Deviance: 21.75 AIC: 25.75
ts <- glm(result ~ temp + sudomotor, data=crps, family="binomial")
fitted(ts)
## 1 2 3 4 5 6 7 8
## 0.3285762 0.2673359 0.7710740 0.8788305 0.2352987 0.5296329 0.2075745 0.2138738
## 9 10 11 12 13 14 15 16
## 0.7986846 0.1608258 0.8476465 0.9531830 0.8337024 0.8940611 0.9514634 0.3541150
## 17 18 19 20 21 22 23 24
## 0.4516967 0.8323854 0.8571742 0.3497955 0.8414309 0.8376037 0.8866688 0.9794836
## 25 26 27 28 29 30 31 32
## 0.8913401 0.3606391 0.9674641 0.9478429 0.9944148 0.8414309 0.7419513 0.9927994
crps$pred_ts=fitted(ts)
a2=ROC(form=result~pred_ts,data=crps,plot="ROC")
a2
## $res
## sens spec pvp pvn lr.eta
## 1.00000000 0.0 NaN 0.31250000 -Inf
## 0.173275069573077 1.00000000 0.1 0.0000000 0.29032258 0.1732751
## 0.209239759642206 1.00000000 0.2 0.0000000 0.26666667 0.2092398
## 0.214483791429557 1.00000000 0.3 0.0000000 0.24137931 0.2144838
## 0.23302870154027 1.00000000 0.4 0.0000000 0.21428571 0.2330287
## 0.262782171456349 1.00000000 0.5 0.0000000 0.18518519 0.2627822
## 0.326022860155418 1.00000000 0.6 0.0000000 0.15384615 0.3260229
## 0.349682635343628 0.95454545 0.6 0.1428571 0.16000000 0.3496826
## 0.354595786817145 0.90909091 0.6 0.2500000 0.16666667 0.3545958
## 0.362074586149775 0.86363636 0.6 0.3333333 0.17391304 0.3620746
## 0.471934535420382 0.81818182 0.6 0.4000000 0.18181818 0.4719345
## 0.56861215173662 0.81818182 0.7 0.3636364 0.14285714 0.5686122
## 0.791626817458375 0.77272727 0.7 0.4166667 0.15000000 0.7916268
## 0.814566205549124 0.77272727 0.8 0.3846154 0.10526316 0.8145662
## 0.834466311808699 0.77272727 0.9 0.3571429 0.05555556 0.8344663
## 0.856392123946376 0.72727273 0.9 0.4000000 0.05882353 0.8563921
## 0.857197787077499 0.68181818 0.9 0.4375000 0.06250000 0.8571978
## 0.859562280273369 0.63636364 0.9 0.4705882 0.06666667 0.8595623
## 0.861849938574947 0.54545455 0.9 0.5263158 0.07692308 0.8618499
## 0.865498478461376 0.50000000 0.9 0.5500000 0.08333333 0.8654985
## 0.87093290859367 0.45454545 0.9 0.5714286 0.09090909 0.8709329
## 0.882591796596249 0.45454545 1.0 0.5454545 0.00000000 0.8825918
## 0.886581200882974 0.40909091 1.0 0.5652174 0.00000000 0.8865812
## 0.888902145578829 0.36363636 1.0 0.5833333 0.00000000 0.8889021
## 0.890234814457609 0.31818182 1.0 0.6000000 0.00000000 0.8902348
## 0.913827038521329 0.27272727 1.0 0.6153846 0.00000000 0.9138270
## 0.915237923659119 0.22727273 1.0 0.6296296 0.00000000 0.9152379
## 0.915900669771847 0.18181818 1.0 0.6428571 0.00000000 0.9159007
## 0.921225203611029 0.13636364 1.0 0.6551724 0.00000000 0.9212252
## 0.925465709196979 0.09090909 1.0 0.6666667 0.00000000 0.9254657
## 0.929917790556103 0.04545455 1.0 0.6774194 0.00000000 0.9299178
## 0.930440863252602 0.00000000 1.0 0.6875000 NaN 0.9304409
##
## $AUC
## [1] 0.8818182
##
## $lr
##
## Call: glm(formula = form, family = binomial, data = data)
##
## Coefficients:
## (Intercept) pred_ts
## -2.364 4.986
##
## Degrees of Freedom: 31 Total (i.e. Null); 30 Residual
## Null Deviance: 39.75
## Residual Deviance: 28.3 AIC: 32.3
tw <- glm(result ~ temp + wbbs, data=crps, family="binomial")
fitted(tw)
## 1 2 3 4 5 6 7 8
## 0.2880471 0.2257393 0.1893135 0.3547131 0.8692218 0.5062965 0.1676404 0.1736226
## 9 10 11 12 13 14 15 16
## 0.8070366 0.1244202 0.9185306 0.9795779 0.9093965 0.3942980 0.9787148 0.3147383
## 17 18 19 20 21 22 23 24
## 0.4198009 0.9085219 0.3079403 0.9253712 0.9144876 0.9119753 0.9428395 0.9920114
## 25 26 27 28 29 30 31 32
## 0.9456217 0.9289396 0.7263933 0.9768802 0.9981502 0.9144876 0.9877316 0.9975408
crps$pred_tw=fitted(tw)
a3=ROC(form=result~pred_tw,data=crps,plot="ROC")
a3
## $res
## sens spec pvp pvn lr.eta
## 1.00000000 0.0 NaN 0.31250000 -Inf
## 0.142859326208696 1.00000000 0.1 0.0000000 0.29032258 0.1428593
## 0.173546638254637 1.00000000 0.2 0.0000000 0.26666667 0.1735466
## 0.178181417309418 1.00000000 0.3 0.0000000 0.24137931 0.1781814
## 0.190797297660939 1.00000000 0.4 0.0000000 0.21428571 0.1907973
## 0.222680425778328 1.00000000 0.5 0.0000000 0.18518519 0.2226804
## 0.285563139317172 1.00000000 0.6 0.0000000 0.15384615 0.2855631
## 0.30774403425921 0.95454545 0.6 0.1428571 0.16000000 0.3077440
## 0.315539557571566 0.90909091 0.6 0.2500000 0.16666667 0.3155396
## 0.363396801739424 0.90909091 0.7 0.2222222 0.13043478 0.3633968
## 0.413615254259296 0.86363636 0.7 0.3000000 0.13636364 0.4136153
## 0.447022681811059 0.81818182 0.7 0.3636364 0.14285714 0.4470227
## 0.562096980976581 0.81818182 0.8 0.3333333 0.10000000 0.5620970
## 0.806318187506068 0.77272727 0.8 0.3846154 0.10526316 0.8063182
## 0.864989085380419 0.77272727 0.9 0.3571429 0.05555556 0.8649891
## 0.899328765046921 0.77272727 1.0 0.3333333 0.00000000 0.8993288
## 0.916817934760421 0.72727273 1.0 0.3750000 0.00000000 0.9168179
## 0.917173802540877 0.68181818 1.0 0.4117647 0.00000000 0.9171738
## 0.918215038603007 0.63636364 1.0 0.4444444 0.00000000 0.9182150
## 0.919217944921909 0.54545455 1.0 0.5000000 0.00000000 0.9192179
## 0.920808366515454 0.50000000 1.0 0.5238095 0.00000000 0.9208084
## 0.923434246921541 0.45454545 1.0 0.5454545 0.00000000 0.9234342
## 0.924772127906817 0.40909091 1.0 0.5652174 0.00000000 0.9247721
## 0.929781045650915 0.36363636 1.0 0.5833333 0.00000000 0.9297810
## 0.930745903818899 0.31818182 1.0 0.6000000 0.00000000 0.9307459
## 0.940771898691554 0.27272727 1.0 0.6153846 0.00000000 0.9407719
## 0.941315998566666 0.22727273 1.0 0.6296296 0.00000000 0.9413160
## 0.941570355991776 0.18181818 1.0 0.6428571 0.00000000 0.9415704
## 0.943922720663524 0.13636364 1.0 0.6551724 0.00000000 0.9439227
## 0.945121512982558 0.09090909 1.0 0.6666667 0.00000000 0.9451215
## 0.946634616441523 0.04545455 1.0 0.6774194 0.00000000 0.9466346
## 0.946798947019223 0.00000000 1.0 0.6875000 NaN 0.9467989
##
## $AUC
## [1] 0.9272727
##
## $lr
##
## Call: glm(formula = form, family = binomial, data = data)
##
## Coefficients:
## (Intercept) pred_tw
## -2.457 5.346
##
## Degrees of Freedom: 31 Total (i.e. Null); 30 Residual
## Null Deviance: 39.75
## Residual Deviance: 24.74 AIC: 28.74
optimal_lr.eta=function(x){
no=which.max(x$res$sens+x$res$spec)[1]
result=x$res$lr.eta[no]
result
}
optimal_lr.eta(a1)
## [1] 0.770783
optimal_lr.eta(a2)
## [1] 0.8344663
optimal_lr.eta(a3)
## [1] 0.8993288
optimal_cutpoint=function(x){
y=optimal_lr.eta(x)
b0=unname(x$lr$coeff[1])
b1=unname(x$lr$coeff[2])
result=(-log(1/y-1)-b0)/b1
result
}
optimal_cutpoint(a1)
## [1] 0.6890451
optimal_cutpoint(a2)
## [1] 0.7986846
optimal_cutpoint(a3)
## [1] 0.8692218
plot(a1$res$sens~I(1-a1$res$spec),type="l")
par(new=TRUE)
plot(a2$res$sens~I(1-a2$res$spec),type="l",axes=F,xlab="",ylab="",col="red")
par(new=TRUE)
plot(a3$res$sens~I(1-a2$res$spec),type="l",axes=F,xlab="",ylab="",col="blue")
b1=roc(result~pred_tsw,crps,ci=T,percent=T)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
b1
##
## Call:
## roc.formula(formula = result ~ pred_tsw, data = crps, ci = T, percent = T)
##
## Data: pred_tsw in 10 controls (result 0) < 22 cases (result 1).
## Area under the curve: 93.64%
## 95% CI: 85.75%-100% (DeLong)
b2=roc(result~pred_ts,crps,ci=T,percent=T)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
b2
##
## Call:
## roc.formula(formula = result ~ pred_ts, data = crps, ci = T, percent = T)
##
## Data: pred_ts in 10 controls (result 0) < 22 cases (result 1).
## Area under the curve: 88.18%
## 95% CI: 75.28%-100% (DeLong)
b3=roc(result~pred_tw,crps,ci=T,percent=T)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
b3
##
## Call:
## roc.formula(formula = result ~ pred_tw, data = crps, ci = T, percent = T)
##
## Data: pred_tw in 10 controls (result 0) < 22 cases (result 1).
## Area under the curve: 92.73%
## 95% CI: 84.1%-100% (DeLong)
plot(b1)
plot(b2,add=TRUE,col="red")
plot(b3,add=TRUE,col="blue")
roc.test(b1,b2,plot=T)
##
## DeLong's test for two correlated ROC curves
##
## data: b1 and b2
## Z = 0.94756, p-value = 0.3434
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 93.63636 88.18182
roc.test(b1,b3,plot=T)
##
## DeLong's test for two correlated ROC curves
##
## data: b1 and b3
## Z = 0.36769, p-value = 0.7131
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 93.63636 92.72727