Multivariable model to predict PONV after cesarean
Gabriel M N Guimaraes e Helga B G Silva
30 December 2018
library("bnlearn")##
## Attaching package: 'bnlearn'## The following object is masked from 'package:stats':
##
## sigmalibrary(openxlsx)
library(BayesFactor)## Loading required package: coda## Loading required package: Matrix## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************##
## Attaching package: 'BayesFactor'## The following object is masked from 'package:bnlearn':
##
## comparelibrary(e1071)##
## Attaching package: 'e1071'## The following object is masked from 'package:bnlearn':
##
## imputelibrary(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varlibrary(ResourceSelection)## ResourceSelection 0.3-2 2017-02-28read.xlsx2 = function(bdata){
bdata=read.xlsx(bdata)
nv=length(bdata)
ischar=rep(FALSE,nv)
for(i in c(1:nv)){
if(is.character(bdata[,i][0])){
bdata[,i]=factor(bdata[,i])
}
}
bdata
}
sm = function (var,group){
print("Summary statistics for ")
print("Median (min - max)")
print(paste(round(median(var),2),"(",round(min(var),2),"-",round(max(var),2),")"))
print("Mean(sd)")
print(paste(round(mean(var),2),"(",round(sd(var),2),")"))
}
TwoByTwoBFc = function (x,y){
a=sum(x & y,na.rm=TRUE)
b=sum(x & !y,na.rm=TRUE)
c=sum(!x & y,na.rm=TRUE)
d=sum(!x & !y,na.rm=TRUE)
table=as.matrix(data.frame(yes=c(a,c),no=c(b,d)))
bf=extractBF(contingencyTableBF(table, sampleType = "indepMulti", fixedMargin = "cols"))[1]$bf
rownames(table)=c("PONV","NO PONV")
r1=a/(a+c)
r2=b/(b+d)
print(paste("Risco no grupo sim",r1))
print(paste("Risco no grupo não",r2))
print(paste("Risco relativo",r1/r2))
print("Fator bayesiano:")
print(bf)
print("p-valor, teste exato de Fisher:")
print(fisher.test(table)[1]$`p.value`)
print(substitute(y))
table
}
TwoByTwoBFd = function (x,y){
a=sum(x & y,na.rm=TRUE)
b=sum(x & !y,na.rm=TRUE)
c=sum(!x & y,na.rm=TRUE)
d=sum(!x & !y,na.rm=TRUE)
table=as.matrix(data.frame(yes=c(a,c),no=c(b,d)))
bf=extractBF(contingencyTableBF(table, sampleType = "indepMulti", fixedMargin = "cols"))[1]$bf
bf
}Load train and test data sets
traini=read.xlsx2("D:\\Doutorado Gabriel\\train.xlsx")
testi=read.xlsx2("D:\\Doutorado Gabriel\\test.xlsx")
train=droplevels(subset(traini,traini$excluded==FALSE))
test=droplevels(subset(testi,testi$excluded==FALSE))
train$younger=train$age<25
train$morphine.get.80=train$morphine>79
train$neversmoked=train$smoking=="5.never.smoked"
train$preterm=train$gestational.age<38
train$previous.surgery.and.ponv=train$previous.ponv=="previous.surgery.and.ponv"
test$younger=test$age<25
test$morphine.get.80=test$morphine>79
test$neversmoked=test$smoking=="5.never.smoked"
test$preterm=test$gestational.age<38
test$previous.surgery.and.ponv=test$previous.ponv=="previous.surgery.and.ponv"
train$apfel=2+(train$previous.surgery.and.ponv | train$cinetosis)+train$neversmoked
test$apfel=2+(test$previous.surgery.and.ponv | test$cinetosis)+test$neversmoked
test$apfelp=test$apfel*0.2
train$apfelp=train$apfel*0.2
train$lowbupi=train$bupivacaine<13
test$lowbupi=test$bupivacaine<13Descriptive statistics
Variables’ detailed description
| varname | description |
|---|---|
| id | unique case identifier - HIPAA compatible |
| excluded | TRUE if case was excluded |
| exclusion.reason | Reason for exclusion |
| age | Patient Age in Years |
| previous.ponv | If PONV after previous surgery |
| gestational.age | Gestational duration in weeks |
| cinetosis | TRUE if patient refers cinetosis |
| smoking | Cathegorical smoking status from never smoked to smokes |
| nausea1trim | If patient refers significant nausea during the first gestational trimester |
| nausea3trim | If patient refers significant nausea during the third gestational trimester |
| map.basal | Mean arterial pressure before spinal anesthesia in mmHg |
| anesthesia | Anesthesia technique - spinal or other |
| bupivacaine | Heavy Bupivacaine spinal dose in mg |
| sufentanil | Spinal sufentanil in mcg |
| fentanyl | Spinal fentanyl in mcg |
| morphine | Spinal Morphine in mcg |
| intraoperative.nausea | If the patient experienced nausea during the cesarean |
| intraoperative.vomiting | If the patient experienced vomiting during the cesarean |
| lower.intraoperative.map | lower measured intraoperative MAP in mmHg |
| received.ephedrine | If the patient received intravenous ephedrine and if it was rescue or not |
| pacu.nausea | If the patient experienced nausea in the PACU |
| pacu.vomiting | If the patient vomited in the PACU |
| pacu.nausea.intensity | Nausea intensity in PACU |
| pacu.vomiting.number | Number of vomiting episodes in PACU |
| nausea.2.24h | If the patient experienced nausea 2 to 24h after cesarean |
| maximum.nausea.intensity | Maximum nausea intensity 2-24h |
| vomiting.2.24h | If the patient vomited 2-24h after cesarean |
| vomiting.number.2.24h | Number of vomiting Episodes 2-24h after cesarean |
| map.drop | Absolute MAP drop after spinal anesthesia |
| map.drop.proportion | (Basal MAP-lower MAP)/Basal MAP |
| ponv | Either nausea or vomiting from PACU to 24h after cesarean - not IONV |
| ionv | Nausea or vomiting during the cesarean |
| younger | Patients <25 years |
| morphine.get.80 | Spinal Morphine >79mcg |
| neversmoked | Patients who never smoked |
| preterm | Gestational weeks <38 |
| previous.surgery.and.ponv | Patients who had a previous surgery and vomited |
| apfel | Number of Apfel’s PONV risk factors |
| apfelp | Apfel’s predicted probability |
| lowbupi | Spinal Bupivacaine <13mg |
Development and test samples
tsis=length(traini$id)
tsfs=length(train$id)
ttsis=length(testi$id)
ttsfs=length(test$id)
print("Summary")## [1] "Summary"print(paste("Train sample initial size",tsis,"Train sample final size=",tsfs,"Excluded=",tsis-tsfs,"Proportion excluded=",(tsis-tsfs)/tsis))## [1] "Train sample initial size 260 Train sample final size= 250 Excluded= 10 Proportion excluded= 0.0384615384615385"print(paste("Test sample initial size",ttsis,"Test sample final size=",ttsfs,"Excluded=",ttsis-ttsfs,"Proportion excluded=",(ttsis-ttsfs)/ttsis))## [1] "Test sample initial size 100 Test sample final size= 96 Excluded= 4 Proportion excluded= 0.04"Variables description
Numerical variables
meanstrain=c(mean(train$age),mean(train$gestational.age),mean(train$map.basal),mean(train$bupivacaine),mean(train$sufentanil),mean(train$fentanyl),mean(train$morphine),mean(train$lower.intraoperative.map),mean(train$pacu.nausea.intensity),mean(train$pacu.vomiting.number),mean(train$maximum.nausea.intensity,na.rm=TRUE),mean(train$vomiting.number.2.24h),mean(train$map.drop),mean(train$map.drop.proportion))
meanstest=c(mean(test$age),mean(test$gestational.age),mean(test$map.basal),mean(test$bupivacaine),mean(test$sufentanil),mean(test$fentanyl),mean(test$morphine),mean(test$lower.intraoperative.map),mean(test$pacu.nausea.intensity),mean(test$pacu.vomiting.number),mean(test$maximum.nausea.intensity,na.rm=TRUE),mean(test$vomiting.number.2.24h),mean(test$map.drop),mean(test$map.drop.proportion))
sdstrain=c(sd(train$age),sd(train$gestational.age),sd(train$map.basal),sd(train$bupivacaine),sd(train$sufentanil),sd(train$fentanyl),sd(train$morphine),sd(train$lower.intraoperative.map),sd(train$pacu.nausea.intensity),sd(train$pacu.vomiting.number),sd(train$maximum.nausea.intensity,na.rm=TRUE),sd(train$vomiting.number.2.24h),sd(train$map.drop),sd(train$map.drop.proportion))
sdstest=c(sd(test$age),sd(test$gestational.age),sd(test$map.basal),sd(test$bupivacaine),sd(test$sufentanil),sd(test$fentanyl),sd(test$morphine),sd(test$lower.intraoperative.map),sd(test$pacu.nausea.intensity),sd(test$pacu.vomiting.number),sd(test$maximum.nausea.intensity),sd(test$vomiting.number.2.24h),sd(test$map.drop),sd(test$map.drop.proportion))
shapirotrain=c(unname(shapiro.test(train$age)[2]$`p.value`),shapiro.test(train$gestational.age)[2]$`p.value`,shapiro.test(train$map.basal)[2]$`p.value`,shapiro.test(train$bupivacaine)[2]$`p.value`,shapiro.test(train$sufentanil)[2]$`p.value`,shapiro.test(train$fentanyl)[2]$`p.value`,shapiro.test(train$morphine)[2]$`p.value`,shapiro.test(train$lower.intraoperative.map)[2]$`p.value`,shapiro.test(train$pacu.nausea.intensity)[2]$`p.value`,shapiro.test(train$pacu.vomiting.number)[2]$`p.value`,shapiro.test(train$maximum.nausea.intensity)[2]$`p.value`,shapiro.test(train$vomiting.number.2.24h)[2]$`p.value`,shapiro.test(train$map.drop)[2]$`p.value`,shapiro.test(train$map.drop.proportion)[2]$`p.value`)
shapirotest=c(shapiro.test(test$age)[2]$`p.value`,shapiro.test(test$gestational.age)[2]$`p.value`,shapiro.test(test$map.basal)[2]$`p.value`,shapiro.test(test$bupivacaine)[2]$`p.value`,shapiro.test(test$sufentanil)[2]$`p.value`,shapiro.test(test$fentanyl)[2]$`p.value`,shapiro.test(test$morphine)[2]$`p.value`,shapiro.test(test$lower.intraoperative.map)[2]$`p.value`,shapiro.test(test$pacu.nausea.intensity)[2]$`p.value`,shapiro.test(test$pacu.vomiting.number)[2]$`p.value`,shapiro.test(test$maximum.nausea.intensity)[2]$`p.value`,shapiro.test(test$vomiting.number.2.24h)[2]$`p.value`,shapiro.test(test$map.drop)[2]$`p.value`,shapiro.test(test$map.drop.proportion)[2]$`p.value`)
wilcoxtest=c(wilcox.test(train$age,test$age)[3]$`p.value`,wilcox.test(train$gestational.age,test$gestational.age)[3]$`p.value`,wilcox.test(train$map.basal,test$map.basal)[3]$`p.value`,wilcox.test(train$bupivacaine,test$bupivacaine)[3]$`p.value`,wilcox.test(train$sufentanil,test$sufentanil)[3]$`p.value`,wilcox.test(train$fentanyl,test$fentanyl)[3]$`p.value`,wilcox.test(train$morphine,test$morphine)[3]$`p.value`,wilcox.test(train$lower.intraoperative.map,test$lower.intraoperative.map)[3]$`p.value`,wilcox.test(train$pacu.nausea.intensity,test$pacu.nausea.intensity)[3]$`p.value`,wilcox.test(train$pacu.vomiting.number,test$pacu.vomiting.number)[3]$`p.value`,wilcox.test(train$maximum.nausea.intensity,test$maximum.nausea.intensity)[3]$`p.value`,wilcox.test(train$vomiting.number.2.24h,test$vomiting.number.2.24h)[3]$`p.value`,wilcox.test(train$map.drop,test$map.drop)[3]$`p.value`,wilcox.test(train$map.drop.proportion,test$map.drop.proportion)[3]$`p.value`)
parametric=data.frame(Train.mean=round(meanstrain,digits=2),Train.sd=round(sdstrain,digits=2),Train.Shapiro=shapirotrain,Test.mean=round(meanstest,digits=2),Test.sd=round(sdstest,digits=2))
mediantrain=c(median(train$age),median(train$gestational.age),median(train$map.basal),median(train$bupivacaine),median(train$sufentanil),median(train$fentanyl),median(train$morphine),median(train$lower.intraoperative.map),median(train$pacu.nausea.intensity),median(train$pacu.vomiting.number),median(train$maximum.nausea.intensity,na.rm=TRUE),median(train$vomiting.number.2.24h),median(train$map.drop),median(train$map.drop.proportion))
mediantest=c(median(test$age),median(test$gestational.age),median(test$map.basal),median(test$bupivacaine),median(test$sufentanil),median(test$fentanyl),median(test$morphine),median(test$lower.intraoperative.map),median(test$pacu.nausea.intensity),median(test$pacu.vomiting.number),median(test$maximum.nausea.intensity,na.rm=TRUE),median(test$vomiting.number.2.24h),median(test$map.drop),median(test$map.drop.proportion))
mintrain=c(min(train$age),min(train$gestational.age),min(train$map.basal),min(train$bupivacaine),min(train$sufentanil),min(train$fentanyl),min(train$morphine),min(train$lower.intraoperative.map),min(train$pacu.nausea.intensity),min(train$pacu.vomiting.number),min(train$maximum.nausea.intensity,na.rm=TRUE),min(train$vomiting.number.2.24h),min(train$map.drop),min(train$map.drop.proportion))
mintest=c(min(test$age),min(test$gestational.age),min(test$map.basal),min(test$bupivacaine),min(test$sufentanil),min(test$fentanyl),min(test$morphine),min(test$lower.intraoperative.map),min(test$pacu.nausea.intensity),min(test$pacu.vomiting.number),min(test$maximum.nausea.intensity,na.rm=TRUE),min(test$vomiting.number.2.24h),min(test$map.drop),min(test$map.drop.proportion))
maxtrain=c(max(train$age),max(train$gestational.age),max(train$map.basal),max(train$bupivacaine),max(train$sufentanil),max(train$fentanyl),max(train$morphine),max(train$lower.intraoperative.map),max(train$pacu.nausea.intensity),max(train$pacu.vomiting.number),max(train$maximum.nausea.intensity,na.rm=TRUE),max(train$vomiting.number.2.24h),max(train$map.drop),max(train$map.drop.proportion))
maxtest=c(max(test$age),max(test$gestational.age),max(test$map.basal),max(test$bupivacaine),max(test$sufentanil),max(test$fentanyl),max(test$morphine),max(test$lower.intraoperative.map),max(test$pacu.nausea.intensity),max(test$pacu.vomiting.number),max(test$maximum.nausea.intensity,na.rm=TRUE),max(test$vomiting.number.2.24h),max(test$map.drop),max(test$map.drop.proportion))
IQRtrain=c(IQR(train$age),IQR(train$gestational.age),IQR(train$map.basal),IQR(train$bupivacaine),IQR(train$sufentanil),IQR(train$fentanyl),IQR(train$morphine),IQR(train$lower.intraoperative.map),IQR(train$pacu.nausea.intensity),IQR(train$pacu.vomiting.number),IQR(train$maximum.nausea.intensity,na.rm=TRUE),IQR(train$vomiting.number.2.24h),IQR(train$map.drop),IQR(train$map.drop.proportion))
IQRtest=c(IQR(test$age),IQR(test$gestational.age),IQR(test$map.basal),IQR(test$bupivacaine),IQR(test$sufentanil),IQR(test$fentanyl),IQR(test$morphine),IQR(test$lower.intraoperative.map),IQR(test$pacu.nausea.intensity),IQR(test$pacu.vomiting.number),IQR(test$maximum.nausea.intensity,na.rm=TRUE),IQR(test$vomiting.number.2.24h),IQR(test$map.drop),IQR(test$map.drop.proportion))
non.parametric=data.frame(Train.median=round(mediantrain,digits=1),Train.min=round(mintrain,digits=1),Train.max=round(maxtrain,digits=1),Train.iqr=round(IQRtrain,digits=1),Test.median=round(mediantest,digits=1),Test.min=round(mintest,digits=1),Test.max=round(maxtest,digits=1),Test.iqr=round(IQRtest,digits=1))
row.names(parametric)=c("age","gestational.age","map.basal","bupivacaine","sufentanil","fentanyl","morphine","lower.intraoperative.map","pacu.nausea.intensity","pacu.vomiting.number","maximum.nausea.intensity","vomiting.number.2.24h","map.drop","map.drop.proportion")
row.names(non.parametric)=c("age","gestational.age","map.basal","bupivacaine","sufentanil","fentanyl","morphine","lower.intraoperative.map","pacu.nausea.intensity","pacu.vomiting.number","maximum.nausea.intensity","vomiting.number.2.24h","map.drop","map.drop.proportion")
#x_htmlnp <- knitr::kable(non.parametric, "html")
#x_htmlp <- knitr::kable(parametric, "html")
x_htmlnp <- knitr::kable(non.parametric, "pandoc")
x_htmlp <- knitr::kable(parametric, "pandoc")
#kable_styling(kable_input = x_htmlp, "striped")
#kable_styling(kable_input = x_htmlnp, "striped")| Train.mean | Train.sd | Train.Shapiro | Test.mean | Test.sd | |
|---|---|---|---|---|---|
| age | 29.07 | 7.15 | 0.0029885 | 29.58 | 6.50 |
| gestational.age | 38.18 | 2.27 | 0.0000000 | 38.23 | 2.21 |
| map.basal | 92.28 | 14.81 | 0.0002926 | 91.91 | 16.27 |
| bupivacaine | 11.54 | 1.61 | 0.0000000 | 11.50 | 1.60 |
| sufentanil | 1.40 | 1.20 | 0.0000000 | 1.47 | 1.18 |
| fentanyl | 8.46 | 13.15 | 0.0000000 | 8.12 | 13.42 |
| morphine | 89.04 | 15.31 | 0.0000000 | 89.69 | 14.83 |
| lower.intraoperative.map | 66.61 | 16.67 | 0.0003169 | 65.75 | 17.06 |
| pacu.nausea.intensity | 1.08 | 3.04 | 0.0000000 | 1.25 | 3.24 |
| pacu.vomiting.number | 0.08 | 0.39 | 0.0000000 | 0.10 | 0.45 |
| maximum.nausea.intensity | 0.93 | 2.68 | 0.0000000 | 1.05 | 2.80 |
| vomiting.number.2.24h | 0.25 | 1.04 | 0.0000000 | 0.31 | 1.18 |
| map.drop | 25.67 | 16.64 | 0.0000898 | 26.16 | 20.57 |
| map.drop.proportion | 0.27 | 0.16 | 0.0038197 | 0.27 | 0.20 |
| Train.median | Train.min | Train.max | Train.iqr | Test.median | Test.min | Test.max | Test.iqr | |
|---|---|---|---|---|---|---|---|---|
| age | 29.0 | 14.0 | 47.0 | 12.0 | 30.0 | 16.0 | 43.0 | 10.0 |
| gestational.age | 38.0 | 30.0 | 42.0 | 3.0 | 39.0 | 31.0 | 42.0 | 3.0 |
| map.basal | 90.5 | 59.0 | 142.0 | 17.0 | 90.5 | 56.0 | 129.0 | 25.2 |
| bupivacaine | 12.0 | 7.0 | 15.0 | 2.5 | 12.0 | 7.0 | 15.0 | 2.5 |
| sufentanil | 2.0 | 0.0 | 3.0 | 2.5 | 2.0 | 0.0 | 3.0 | 2.5 |
| fentanyl | 0.0 | 0.0 | 50.0 | 25.0 | 0.0 | 0.0 | 50.0 | 25.0 |
| morphine | 100.0 | 0.0 | 100.0 | 20.0 | 100.0 | 0.0 | 100.0 | 20.0 |
| lower.intraoperative.map | 65.0 | 29.0 | 125.0 | 23.0 | 66.0 | 27.0 | 118.0 | 21.2 |
| pacu.nausea.intensity | 0.0 | 0.0 | 10.0 | 0.0 | 0.0 | 0.0 | 10.0 | 0.0 |
| pacu.vomiting.number | 0.0 | 0.0 | 3.0 | 0.0 | 0.0 | 0.0 | 3.0 | 0.0 |
| maximum.nausea.intensity | 0.0 | 0.0 | 10.0 | 0.0 | 0.0 | 0.0 | 10.0 | 0.0 |
| vomiting.number.2.24h | 0.0 | 0.0 | 7.0 | 0.0 | 0.0 | 0.0 | 7.0 | 0.0 |
| map.drop | 23.0 | -9.0 | 70.0 | 24.0 | 25.0 | -14.0 | 79.0 | 29.2 |
| map.drop.proportion | 0.3 | -0.1 | 0.7 | 0.3 | 0.3 | -0.2 | 0.7 | 0.3 |
Categorical variables
presenttrain=c(sum(train$cinetosis),sum(train$nausea1trim),sum(train$nausea3trim),sum(train$ionv),sum(train$younger),sum(train$preterm),sum(train$neversmoked),sum(train$lowbupi))
totalTrain=c(length(train$cinetosis),length(train$nausea1trim),length(train$nausea3trim),length(train$ionv),length(train$younger),length(train$preterm),length(train$neversmoked),length(train$lowbupi))
absentTrain=totalTrain-presenttrain
proportionPresentTrain=presenttrain/totalTrain
presenttest=c(sum(test$cinetosis),sum(test$nausea1trim),sum(test$nausea3trim),sum(test$ionv),sum(test$younger),sum(test$preterm),sum(test$neversmoked),sum(test$lowbupi))
totaltest=c(length(test$cinetosis),length(test$nausea1trim),length(test$nausea3trim),length(test$ionv),length(test$younger),length(test$preterm),length(test$neversmoked),length(test$lowbupi))
absenttest=totaltest-presenttest
proportionPresenttest=presenttest/totaltest
train$sample=rep("train",length(train$id))
test$sample=rep("test",length(test$id))
db=rbind(train,test)
pvaluelist2=c(fisher.test(db$cinetosis,db$sample)[1]$`p.value`,fisher.test(db$nausea1trim,db$sample)[1]$`p.value`,fisher.test(db$nausea3trim,db$sample)[1]$`p.value`,fisher.test(db$ionv,db$sample)[1]$`p.value`,fisher.test(db$younger,db$sample)[1]$`p.value`,fisher.test(db$preterm,db$sample)[1]$`p.value`,fisher.test(db$neversmoked,db$sample)[1]$`p.value`,fisher.test(db$lowbupi,db$sample)[1]$`p.value`)
dichotomic=data.frame(yes.train=presenttrain,no.train=absentTrain,prop=proportionPresentTrain,yes.test=presenttest,no.test=absenttest,prop.test=proportionPresenttest,p.value=pvaluelist2)
row.names(dichotomic)=c("Cinetosis","Nausea 1 Trimester","Nausea 3 Trimester","Intraop.Nausea","Age<25","Gest<38w","Never.smoked","Low.Bupi")
#x_html <- knitr::kable(dichotomic, "html")
x_html <- knitr::kable(dichotomic, "pandoc")
kable_styling(kable_input = x_html, "striped", full_width = FALSE)| yes.train | no.train | prop | yes.test | no.test | prop.test | p.value | |
|---|---|---|---|---|---|---|---|
| Cinetosis | 66 | 184 | 0.264 | 26 | 70 | 0.2708333 | 0.8926667 |
| Nausea 1 Trimester | 179 | 71 | 0.716 | 69 | 27 | 0.7187500 | 1.0000000 |
| Nausea 3 Trimester | 99 | 151 | 0.396 | 41 | 55 | 0.4270833 | 0.6257416 |
| Intraop.Nausea | 110 | 140 | 0.440 | 49 | 47 | 0.5104167 | 0.2783978 |
| Age<25 | 77 | 173 | 0.308 | 23 | 73 | 0.2395833 | 0.2345593 |
| Gest<38w | 75 | 175 | 0.300 | 31 | 65 | 0.3229167 | 0.6971843 |
| Never.smoked | 203 | 47 | 0.812 | 78 | 18 | 0.8125000 | 1.0000000 |
| Low.Bupi | 216 | 34 | 0.864 | 85 | 11 | 0.8854167 | 0.7217297 |
table(db$sample,db$previous.ponv)##
## no.previous.surgery previous.surgery.and.ponv
## test 66 19
## train 172 45
##
## previous.surgery.no.ponv
## test 11
## train 33table(db$sample,db$smoking)##
## 1.smokes 3.stopped.1.3.months 4.stopped.6.months 5.never.smoked
## test 5 1 12 78
## train 13 2 32 203table(db$sample,db$received.ephedrine)##
## Not prophylactic yes.rescue
## test 9 30 57
## train 25 79 146Risk factors for PONV (Training set)
Dichotomic risk factos analysis
presentandponv=c(sum(train$cinetosis & train$ponv),sum(train$nausea1trim & train$ponv),sum(train$nausea3trim & train$ponv),sum(train$ionv & train$ponv),sum(train$younger & train$ponv),sum(train$preterm & train$ponv),sum(train$neversmoked & train$ponv),sum(train$lowbupi & train$ponv))
presentnoponv=c(sum(train$cinetosis & !train$ponv),sum(train$nausea1trim & !train$ponv),sum(train$nausea3trim & !train$ponv),sum(train$ionv & !train$ponv),sum(train$younger & !train$ponv),sum(train$preterm & !train$ponv),sum(train$neversmoked & !train$ponv),sum(train$lowbupi & !train$ponv))
absentandponv=c(sum(!train$cinetosis & train$ponv),sum(!train$nausea1trim & train$ponv),sum(!train$nausea3trim & train$ponv),sum(!train$ionv & train$ponv),sum(!train$younger & train$ponv),sum(!train$preterm & train$ponv),sum(!train$neversmoked & train$ponv), sum(!train$lowbupi & train$ponv))
absentandnoponv=c(sum(!train$cinetosis & !train$ponv),sum(!train$nausea1trim & !train$ponv),sum(!train$nausea3trim & !train$ponv),sum(!train$ionv & !train$ponv),sum(!train$younger & !train$ponv),sum(!train$preterm & !train$ponv),sum(!train$neversmoked & !train$ponv), sum(!train$lowbupi & !train$ponv))
totalTrain=c(length(train$cinetosis),length(train$nausea1trim),length(train$nausea3trim),length(train$ionv),length(train$younger),length(train$preterm),length(train$neversmoked),length(train$lowbupi))
risk=presentandponv/(presentandponv+presentnoponv)
norisk=absentandponv/(absentandponv+absentandnoponv)
rr=risk/norisk
presenttest=presenttest=c(sum(test$cinetosis),sum(test$nausea1trim),sum(test$nausea3trim),sum(test$ionv),sum(test$younger),sum(test$preterm),sum(test$neversmoked),sum(test$lowbupi))
totaltest=c(length(test$cinetosis),length(test$nausea1trim),length(test$nausea3trim),length(test$ionv),length(test$younger),length(test$preterm),length(test$neversmoked),length(test$lowbupi))
absenttest=totaltest-presenttest
proportionPresenttest=presenttest/totaltest
pvaluelist3=c(fisher.test(train$cinetosis,train$ponv)[1]$`p.value`,fisher.test(train$nausea1trim,train$ponv)[1]$`p.value`,fisher.test(train$nausea3trim,train$ponv)[1]$`p.value`,fisher.test(train$ionv,train$ponv)[1]$`p.value`,fisher.test(train$younger,train$ponv)[1]$`p.value`,fisher.test(train$preterm,train$ponv)[1]$`p.value`,fisher.test(train$neversmoked,train$ponv)[1]$`p.value`,fisher.test(train$lowbupi,train$ponv)[1]$`p.value`)
bflist1=c(TwoByTwoBFd(train$cinetosis,train$ponv),TwoByTwoBFd(train$nausea1trim,train$ponv),TwoByTwoBFd(train$nausea3trim,train$ponv),TwoByTwoBFd(train$ionv,train$ponv),TwoByTwoBFd(train$younger,train$ponv),TwoByTwoBFd(train$preterm,train$ponv),TwoByTwoBFd(train$neversmoked,train$ponv),TwoByTwoBFd(train$lowbupi,train$ponv))
odds=c(unname(fisher.test(train$cinetosis,train$ponv)[3]$`estimate`),unname(fisher.test(train$nausea1trim,train$ponv)[3]$`estimate`),unname(fisher.test(train$nausea3trim,train$ponv)[3]$`estimate`),unname(fisher.test(train$ionv,train$ponv)[3]$`estimate`),unname(fisher.test(train$younger,train$ponv)[3]$`estimate`),unname(fisher.test(train$preterm,train$ponv)[3]$`estimate`),unname(fisher.test(train$neversmoked,train$ponv)[3]$`estimate`),unname(fisher.test(train$lowbupi,train$ponv)[3]$`estimate`))
odds_ll=c(fisher.test(train$cinetosis,train$ponv)[2]$`conf.int`[1],fisher.test(train$nausea1trim,train$ponv)[2]$`conf.int`[1],fisher.test(train$nausea3trim,train$ponv)[2]$`conf.int`[1],fisher.test(train$ionv,train$ponv)[2]$`conf.int`[1],fisher.test(train$younger,train$ponv)[2]$`conf.int`[1],fisher.test(train$preterm,train$ponv)[2]$`conf.int`[1],fisher.test(train$neversmoked,train$ponv)[2]$`conf.int`[1],fisher.test(train$lowbupi,train$ponv)[2]$`conf.int`[1])
odds_ul=c(fisher.test(train$cinetosis,train$ponv)[2]$`conf.int`[2],fisher.test(train$nausea1trim,train$ponv)[2]$`conf.int`[2],fisher.test(train$nausea3trim,train$ponv)[2]$`conf.int`[2],fisher.test(train$ionv,train$ponv)[2]$`conf.int`[2],fisher.test(train$younger,train$ponv)[2]$`conf.int`[2],fisher.test(train$preterm,train$ponv)[2]$`conf.int`[2],fisher.test(train$neversmoked,train$ponv)[2]$`conf.int`[2],fisher.test(train$lowbupi,train$ponv)[2]$`conf.int`[2])
odds=round(odds,digits=2)
odds_ll=round(odds_ll,digits=2)
odds_ul=round(odds_ul,digits=2)
rr=round(rr,digits=2)
risk=round(risk,digits=2)
odds2=paste(odds,"(",odds_ll,"-",odds_ul,")")
dichotomic=data.frame(ponv=presentandponv,noponv=presentnoponv,p.ponv=risk,RR=rr, odds.95percent.CI=odds2,p.value=pvaluelist3,BF=bflist1)
row.names(dichotomic)=c("Cinetosis","Nausea 1 Trimester","Nausea 3 Trimester","Intraop.Nausea","Age<25","Gest<38w","Never.smoked","Low.bupi")
#x_html <- knitr::kable(dichotomic, "html")
x_html <- knitr::kable(dichotomic, "pandoc")
kable_styling(kable_input = x_html, "striped", full_width = FALSE)| ponv | noponv | p.ponv | RR | odds.95percent.CI | p.value | BF | |
|---|---|---|---|---|---|---|---|
| Cinetosis | 21 | 45 | 0.32 | 2.09 | 2.59 ( 1.27 - 5.25 ) | 0.0061468 | 9.562329e+00 |
| Nausea 1 Trimester | 25 | 154 | 0.14 | 0.41 | 0.32 ( 0.16 - 0.64 ) | 0.0006816 | 6.750437e+01 |
| Nausea 3 Trimester | 24 | 75 | 0.24 | 1.46 | 1.61 ( 0.82 - 3.17 ) | 0.1452093 | 5.861174e-01 |
| Intraop.Nausea | 40 | 70 | 0.36 | 5.66 | 8.24 ( 3.67 - 20.47 ) | 0.0000000 | 1.170844e+07 |
| Age<25 | 25 | 52 | 0.32 | 2.34 | 2.97 ( 1.49 - 5.96 ) | 0.0009826 | 4.594967e+01 |
| Gest<38w | 21 | 54 | 0.28 | 1.75 | 2.04 ( 1.01 - 4.08 ) | 0.0366478 | 1.855045e+00 |
| Never.smoked | 45 | 158 | 0.22 | 2.60 | 3.05 ( 1.03 - 12.32 ) | 0.0402170 | 1.593157e+00 |
| Low.bupi | 49 | 167 | 0.23 | Inf | Inf ( 2.43 - Inf ) | 0.0006875 | 7.605094e+01 |
Numerical risk factors analysis
meanstrainponv=c(mean(train$age[train$ponv==TRUE]),mean(train$gestational.age[train$ponv==TRUE]),mean(train$map.basal[train$ponv==TRUE]),mean(train$bupivacaine[train$ponv==TRUE]),mean(train$sufentanil[train$ponv==TRUE]),mean(train$fentanyl[train$ponv==TRUE]),mean(train$morphine[train$ponv==TRUE]),mean(train$lower.intraoperative.map[train$ponv==TRUE]),mean(train$map.drop[train$ponv==TRUE]),mean(train$map.drop.proportion[train$ponv==TRUE]))
meanstrainnoponv=c(mean(train$age[train$ponv==FALSE]),mean(train$gestational.age[train$ponv==FALSE]),mean(train$map.basal[train$ponv==FALSE]),mean(train$bupivacaine[train$ponv==FALSE]),mean(train$sufentanil[train$ponv==FALSE]),mean(train$fentanyl[train$ponv==FALSE]),mean(train$morphine[train$ponv==FALSE]),mean(train$lower.intraoperative.map[train$ponv==FALSE]),mean(train$map.drop[train$ponv==FALSE]),mean(train$map.drop.proportion[train$ponv==FALSE]))
sdtrainponv=meanstrain=c(sd(train$age[train$ponv==TRUE]),sd(train$gestational.age[train$ponv==TRUE]),sd(train$map.basal[train$ponv==TRUE]),sd(train$bupivacaine[train$ponv==TRUE]),sd(train$sufentanil[train$ponv==TRUE]),sd(train$fentanyl[train$ponv==TRUE]),sd(train$morphine[train$ponv==TRUE]),sd(train$lower.intraoperative.map[train$ponv==TRUE]),sd(train$map.drop[train$ponv==TRUE]),sd(train$map.drop.proportion[train$ponv==TRUE]))
sdtrainnoponv=c(sd(train$age[train$ponv==FALSE]),sd(train$gestational.age[train$ponv==FALSE]),sd(train$map.basal[train$ponv==FALSE]),sd(train$bupivacaine[train$ponv==FALSE]),sd(train$sufentanil[train$ponv==FALSE]),sd(train$fentanyl[train$ponv==FALSE]),sd(train$morphine[train$ponv==FALSE]),sd(train$lower.intraoperative.map[train$ponv==FALSE]),sd(train$map.drop[train$ponv==FALSE]),sd(train$map.drop.proportion[train$ponv==FALSE]))
plist5=c(wilcox.test(train$age~train$ponv)[3]$`p.value`,wilcox.test(train$gestational.age~train$ponv)[3]$`p.value`,wilcox.test(train$map.basal~train$ponv)[3]$`p.value`,wilcox.test(train$bupivacaine~train$ponv)[3]$`p.value`,wilcox.test(train$sufentanil~train$ponv)[3]$`p.value`,wilcox.test(train$fentanyl~train$ponv)[3]$`p.value`,wilcox.test(train$morphine~train$ponv)[3]$`p.value`,wilcox.test(train$lower.intraoperative.map~train$ponv)[3]$`p.value`,wilcox.test(train$map.drop~train$ponv)[3]$`p.value`,wilcox.test(train$map.drop.proportion~train$ponv)[3]$`p.value`)
ponv=subset(train,train$ponv==TRUE)
nponv=subset(train,train$ponv==FALSE)
bfnum=c(extractBF(ttestBF(ponv$age,nponv$age))$bf,extractBF(ttestBF(ponv$gestational.age,nponv$gestational.age))$bf,extractBF(ttestBF(ponv$map.basal,nponv$map.basal))$bf,extractBF(ttestBF(ponv$bupivacaine,nponv$bupivacaine))$bf,extractBF(ttestBF(ponv$sufentanil,nponv$sufentanil))$bf,extractBF(ttestBF(ponv$fentanyl,nponv$fentanyl))$bf,extractBF(ttestBF(ponv$morphine,nponv$morphine))$bf,extractBF(ttestBF(ponv$lower.intraoperative.map,nponv$lower.intraoperative.map))$bf,extractBF(ttestBF(ponv$map.drop,nponv$map.drop))$bf,extractBF(ttestBF(ponv$map.drop.proportion,nponv$map.drop.proportion))$bf)
numerical=data.frame(mean.ponv=meanstrainponv,sd.ponv=sdtrainponv,mean.no.ponv=meanstrainnoponv,sd.no.ponv=sdtrainnoponv,p.value=plist5,BF=bfnum)
row.names(numerical)=c("age","gestational.age","map.basal","bupivacaine","sufentanil","fentanyl","morphine","lower.intraoperative.map","map.drop","map.drop.proportion")
#x_html <- knitr::kable(numerical, "html")
x_html <- knitr::kable(numerical, "pandoc")
kable_styling(kable_input = x_html, "striped", full_width = FALSE)| mean.ponv | sd.ponv | mean.no.ponv | sd.no.ponv | p.value | BF | |
|---|---|---|---|---|---|---|
| age | 26.1224490 | 6.811854 | 29.7910448 | 7.0573451 | 0.0008506 | 23.5214438 |
| gestational.age | 37.8979592 | 2.608040 | 38.2537313 | 2.1795180 | 0.4448352 | 0.2682916 |
| map.basal | 96.4081633 | 17.643363 | 91.2736318 | 13.9043069 | 0.1557567 | 1.5567891 |
| bupivacaine | 10.7346939 | 2.049141 | 11.7363184 | 1.4213108 | 0.0413993 | 263.7028568 |
| sufentanil | 1.3673469 | 1.206646 | 1.4104478 | 1.2008082 | 0.8228957 | 0.1758287 |
| fentanyl | 10.2040816 | 16.071429 | 8.0348259 | 12.3488372 | 0.6278850 | 0.2813088 |
| morphine | 93.0612245 | 9.618576 | 88.0597015 | 16.2701081 | 0.0486619 | 1.2131889 |
| lower.intraoperative.map | 65.1224490 | 23.624345 | 66.9701493 | 14.5313146 | 0.3053587 | 0.2145840 |
| map.drop | 31.2857143 | 16.904142 | 24.3034826 | 16.3215329 | 0.0086703 | 4.4495144 |
| map.drop.proportion | 0.3306173 | 0.182935 | 0.2575136 | 0.1566235 | 0.0186188 | 6.7189929 |
Feature selection (cathegorical)
BayesFactor::contingencyTableBF(table(train$previous.ponv,train$ponv),sampleType = "indepMulti", fixedMargin = "cols")## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 0.07977785 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, independent multinomialBayesFactor::contingencyTableBF(table(train$received.ephedrine,train$ponv),sampleType = "indepMulti", fixedMargin = "cols")## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 0.4712684 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, independent multinomialBayesFactor::contingencyTableBF(table(train$smoking,train$ponv),sampleType = "indepMulti", fixedMargin = "cols")## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 0.0275246 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, independent multinomialbdbn=data.frame(neversmoked=train$neversmoked,cinetosis=train$cinetosis,nausea1trim=train$nausea1trim, ponv=train$ponv,ionv=train$ionv)
bdbn$younger=as.factor(train$younger)
bdbn$neversmoked=as.factor(bdbn$neversmoked)
bdbn$ionv=as.factor(bdbn$ionv)
bdbn$cinetosis=as.factor(bdbn$cinetosis)
bdbn$nausea1trim=as.factor(bdbn$nausea1trim)
bdbn$ponv=as.factor(bdbn$ponv)
bdbn$preterm=as.factor(train$preterm)
bdbn$ephedrine=as.factor(train$received.ephedrine)
bdbn$hypotension=as.factor(train$map.drop.proportion>0.1)
bdbn$previousponv.and.surgery=as.factor(train$previous.surgery.and.ponv)
bdbn$lowbupi=as.factor(train$lowbupi)
#bdbn
# Blacklist
i=1;
from="";to=""
n=(length(names(bdbn)))
## PONV cannot cause any past variable
for(i in c(1:n)){
#from[i]="ponv";to[i]=names(bdbn)[i]}
from[i]=names(bdbn)[i];to[i]="ponv"}
## Age cannot be caused by any other variable
for(i in c(1:n)){
from[i+n]="younger";to[i+n]=names(bdbn)[i]}
#to[i+n]="younger";from[i+n]=names(bdbn)[i]
to[2*n+1]="ephedrine";from[2*n+1]="hypotension"
to[2*n+2]="ephedrine";from[2*n+2]="neversmoked"
blacklist=data.frame(from=from,to=to)
plot(iamb(bdbn,blacklist=blacklist))
## Step 2
After removing variables excluded by mutual conditional independence
bdbn=data.frame(nausea1trim=train$nausea1trim, ponv=train$ponv,ionv=train$ionv,preterm=train$preterm)
bdbn$ionv=as.factor(bdbn$ionv)
bdbn$preterm=as.factor(bdbn$preterm)
bdbn$nausea1trim=as.factor(bdbn$nausea1trim)
bdbn$ponv=as.factor(bdbn$ponv)
bdbn$ephedrine=as.factor(train$received.ephedrine)
bdbn$hypotension=as.factor(train$map.drop.proportion>0.1)
# Blacklist - temporal and technical knowledge Helps reducing universe of possible models
i=1;
from="";to=""
n=(length(names(bdbn)))
## PONV cannot cause any past variable
for(i in c(1:n)){
to[i]="ponv";from[i]=names(bdbn)[i]}
from[1+n]="ionv";to[1+n]="ephedrine"
to[n+2]="ephedrine";from[n+2]="hypotension"
blacklist=data.frame(from=from,to=to)
BayesNetwork=iamb(bdbn,blacklist=blacklist)
plot(BayesNetwork)
This Bayesian Network is very similar to a NaïveBayes. For classification purposes, it is identical.
NaïveBayes classifier
bdbnt=data.frame(ponv=test$ponv)
bdbnt$ionv=as.factor(test$ionv)
bdbnt$nausea1trim=as.factor(test$nausea1trim)
bdbnt$ponv=as.factor(test$ponv)
bdbnt$preterm=as.factor(test$preterm)
bdbn=data.frame(ponv=test$ponv)
bdbn$ionv=as.factor(test$ionv)
bdbn$nausea1trim=as.factor(test$nausea1trim)
bdbn$ponv=as.factor(test$ponv)
bdbn$preterm=as.factor(test$preterm)
fitnb <-naiveBayes(ponv ~ ., data=bdbn)
nb1=naive.bayes(bdbn,"ponv")
fit=bn.fit(nb1,bdbn)
pred=predict(fit,bdbnt)
table(pred,bdbnt$ponv)##
## pred FALSE TRUE
## FALSE 69 13
## TRUE 4 10tp=sum(pred ==TRUE & bdbnt$ponv==TRUE)
tn=sum(pred ==FALSE & bdbnt$ponv==FALSE)
fp=sum(pred ==TRUE & bdbnt$ponv==FALSE)
fn=sum(pred ==FALSE & bdbnt$ponv==TRUE)
m2=data.frame(naive.p=(predict(fitnb,bdbnt,type="raw")[,2]),observed=bdbnt$ponv)
table(m2)## observed
## naive.p FALSE TRUE
## 0.0315388272614466 23 0
## 0.0584072969116142 10 1
## 0.120777837376385 6 1
## 0.207390399097118 5 1
## 0.2477823063 19 5
## 0.385535426883853 6 5
## 0.581498758763674 4 7
## 0.725773982146764 0 3#m2$naive.p=as.factor(m2$naive.p)
sensitivity.nb=tp/(tp+fn)
specificity.nb=tn/(tn+fp)
fitnb##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## FALSE TRUE
## 0.7604167 0.2395833
##
## Conditional probabilities:
## ionv
## Y FALSE TRUE
## FALSE 0.6027397 0.3972603
## TRUE 0.1304348 0.8695652
##
## nausea1trim
## Y FALSE TRUE
## FALSE 0.2054795 0.7945205
## TRUE 0.5217391 0.4782609
##
## preterm
## Y FALSE TRUE
## FALSE 0.7123288 0.2876712
## TRUE 0.5652174 0.4347826plot.roc(m2$observed,m2$naive.p,percent=TRUE,print.auc=TRUE)
text(x=30,y=20,labels="ROC for NaïveBayes Classifier")
rocnb=roc(m2$observed,m2$naive.p)
ci(roc(m2$observed,m2$naive.p))## 95% CI: 0.7605-0.9297 (DeLong)ci(roc(m2$observed,m2$naive.p),of="sp")## 95% CI (2000 stratified bootstrap replicates):
## se sp.low sp.median sp.high
## 0.0 1.0000 1.0000 1.0000
## 0.1 0.9822 1.0000 1.0000
## 0.2 0.9466 0.9901 1.0000
## 0.3 0.9060 0.9732 1.0000
## 0.4 0.8654 0.9528 0.9937
## 0.5 0.8065 0.9231 0.9874
## 0.6 0.7178 0.8838 0.9690
## 0.7 0.6019 0.8101 0.9380
## 0.8 0.4915 0.6959 0.8890
## 0.9 0.3801 0.5699 0.7628
## 1.0 0.2192 0.3562 0.6164ci(roc(m2$observed,m2$naive.p),of="se")## 95% CI (2000 stratified bootstrap replicates):
## sp se.low se.median se.high
## 0.0 1.0000 1.0000 1.0000
## 0.1 1.0000 1.0000 1.0000
## 0.2 1.0000 1.0000 1.0000
## 0.3 0.9580 1.0000 1.0000
## 0.4 0.8696 0.9791 1.0000
## 0.5 0.7950 0.9372 1.0000
## 0.6 0.7188 0.8717 1.0000
## 0.7 0.6283 0.7946 0.9325
## 0.8 0.5086 0.7075 0.8678
## 0.9 0.3157 0.5600 0.7883
## 1.0 0.0000 0.1304 0.3043score(nb1,bdbnt,type="aic")## [1] -230.3588Multiple logistic regression classifier
mlr1=(glm(ponv~age+previous.surgery.and.ponv+cinetosis+nausea1trim+received.ephedrine+ionv+map.drop.proportion,data=train,family="binomial"))
summary(mlr1)##
## Call:
## glm(formula = ponv ~ age + previous.surgery.and.ponv + cinetosis +
## nausea1trim + received.ephedrine + ionv + map.drop.proportion,
## family = "binomial", data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.59765 -0.50786 -0.17663 -0.04729 2.58399
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.94003 1.23773 2.375 0.01753 *
## age -0.16067 0.03883 -4.138 3.50e-05 ***
## previous.surgery.and.ponvTRUE 1.49376 0.62710 2.382 0.01722 *
## cinetosisTRUE 1.08006 0.45635 2.367 0.01795 *
## nausea1trimTRUE -1.45728 0.45400 -3.210 0.00133 **
## received.ephedrineprophylactic -2.06982 0.75611 -2.737 0.00619 **
## received.ephedrineyes.rescue -3.66654 0.80529 -4.553 5.29e-06 ***
## ionvTRUE 3.63277 0.65432 5.552 2.82e-08 ***
## map.drop.proportion 2.78001 1.39144 1.998 0.04572 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 247.40 on 249 degrees of freedom
## Residual deviance: 154.46 on 241 degrees of freedom
## AIC: 172.46
##
## Number of Fisher Scoring iterations: 6plot.roc(test$ponv,unname(predict(mlr1,test)),percent=TRUE,print.auc=TRUE)
text(x=30,y=20,labels="ROC for Multiple Logistic Regression")
roclogistic=roc(test$ponv,unname(predict(mlr1,test)))
ci(roc(test$ponv,unname(predict(mlr1,test))))## 95% CI: 0.8141-0.9668 (DeLong)ci(roc(test$ponv,unname(predict(mlr1,test))),of="se")## 95% CI (2000 stratified bootstrap replicates):
## sp se.low se.median se.high
## 0.0 1.0000 1.0000 1.0000
## 0.1 1.0000 1.0000 1.0000
## 0.2 1.0000 1.0000 1.0000
## 0.3 1.0000 1.0000 1.0000
## 0.4 0.9565 1.0000 1.0000
## 0.5 0.7826 1.0000 1.0000
## 0.6 0.6946 0.8696 1.0000
## 0.7 0.6522 0.8261 0.9565
## 0.8 0.5652 0.7826 0.9565
## 0.9 0.4783 0.6957 0.8696
## 1.0 0.2609 0.4783 0.7391ci(roc(test$ponv,unname(predict(mlr1,test))),of="sp")## 95% CI (2000 stratified bootstrap replicates):
## se sp.low sp.median sp.high
## 0.0 1.0000 1.0000 1.0000
## 0.1 1.0000 1.0000 1.0000
## 0.2 1.0000 1.0000 1.0000
## 0.3 0.9863 1.0000 1.0000
## 0.4 0.9452 1.0000 1.0000
## 0.5 0.8767 0.9863 1.0000
## 0.6 0.7945 0.9726 1.0000
## 0.7 0.5476 0.8767 1.0000
## 0.8 0.4795 0.7945 0.9863
## 0.9 0.4384 0.5890 0.9041
## 1.0 0.4106 0.5205 0.6575table(test$ponv,unname(predict(mlr1,test)>0.5))##
## FALSE TRUE
## FALSE 73 0
## TRUE 15 8mlr1$aic## [1] 172.4634Apfel’s heuristic
# Female sex = all = 1
# Plan of postoperative opoid (neuraxial?) = same procedure = 1
# Previous PONV or cinetosis
# Smoking ... nerver smoked only?
table(train$apfel,train$ponv)##
## FALSE TRUE
## 2 25 0
## 3 124 28
## 4 52 21roc(test$ponv,test$apfelp)##
## Call:
## roc.default(response = test$ponv, predictor = test$apfelp)
##
## Data: test$apfelp in 73 controls (test$ponv FALSE) < 23 cases (test$ponv TRUE).
## Area under the curve: 0.5989ci(test$ponv,test$apfelp)## 95% CI: 0.4904-0.7073 (DeLong)Apfel’s model’s discrimination: either sensitivity of 42.8% and specificity of 74.1% or sensitivity of 100% and specificity of 12.4%.
Models comparison and model selection
roc.test(rocnb,roc(test$ponv,test$apfelp),method="delong")##
## DeLong's test for two ROC curves
##
## data: rocnb and roc(test$ponv, test$apfelp)
## D = 3.5097, df = 179.39, p-value = 0.0005671
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8451459 0.5988684roc.test(roclogistic,roc(test$ponv,test$apfelp),method="delong")##
## DeLong's test for two correlated ROC curves
##
## data: roclogistic and roc(test$ponv, test$apfelp)
## Z = 4.0918, p-value = 4.281e-05
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8904110 0.5988684roc.test(rocnb,roclogistic,method="delong")##
## DeLong's test for two ROC curves
##
## data: rocnb and roclogistic
## D = -0.77847, df = 188.04, p-value = 0.4373
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8451459 0.8904110roc.test(rocnb,roclogistic,method="bootstrap")##
## Bootstrap test for two ROC curves
##
## data: rocnb and roclogistic
## D = -0.79395, boot.n = 2000, boot.stratified = 1, p-value = 0.4272
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8451459 0.8904110hoslem.test(test$ponv,m2$naive.p)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: test$ponv, m2$naive.p
## X-squared = 1.8854, df = 8, p-value = 0.9843hoslem.test(test$ponv,unname(predict(mlr1,test)))##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: test$ponv, unname(predict(mlr1, test))
## X-squared = -18.767, df = 8, p-value = 1par(mfrow=c(1,3))
plot.roc(m2$observed,m2$naive.p,percent=TRUE,print.auc=TRUE)
text(x=50,y=-20,labels="ROC for NaïveBayes Classifier")
plot.roc(test$ponv,unname(predict(mlr1,test)),percent=TRUE,print.auc=TRUE)
text(x=50,y=-20,labels="Multiple Logistic Regression")
plot.roc(test$ponv,test$apfelp,percent=TRUE,print.auc=TRUE)
text(x=50,y=-20,labels="Simplified Apfel's Model")
We cannot rule out the null hypothesis for the discrimination power difference between models using either Delong’s or bootstrapping methods.
Hosmer-Lemeshow tests high p-values are additional evidence of good model fitness. If p-values < alpha level were obtained, this would be evidence of poor fitness.
More discrimination comparison
mlrp=unname(predict(mlr1,test))
modelcomp=data.frame(observed=test$ponv,logistic=mlrp>0.5,naive=m2$naive.p>0.5)
table(modelcomp)## , , naive = FALSE
##
## logistic
## observed FALSE TRUE
## FALSE 69 0
## TRUE 11 2
##
## , , naive = TRUE
##
## logistic
## observed FALSE TRUE
## FALSE 4 0
## TRUE 4 6modelcomp=data.frame(observed=test$ponv,logistic=mlrp,naive=m2$naive.p)
#Comparisons for 0.5 threshold
logistic=c(8,15)
naive=c(10,13)
naive=c(69,4)
logisti=c(73,0)
## Sensitivity
### Training Sample -
kp = 20 #Naive
np = 49 # total
klp = 18 # Logistic
kpriori = 10 #Naive
npriori = 49 # total
### Test Sample
k = 10+kp #Naive
n = 23+np # total
kl = 8+klp # Logistic
## x-axis for plotting
numSteps = 200
x = seq(0, 1, 1 / numSteps)
## Likelihood (training) function
P = x^kp * (1 - x)^(np - kp)
Pl = x^klp * (1 - x)^(np - klp)
## Likelihood function (test sample)
L = x^k * (1 - x)^(n - k)
Ll = x^kl * (1 - x)^(n - kl)
## Normalization
L = L / sum(L) * numSteps
Ll = Ll / sum(Ll) * numSteps
P = P / sum(P) * numSteps
Pl = Pl / sum(Pl) * numSteps
## Plot prior sensitivity = prevalence = 0.2
plot(x, dbeta(x,4,20), type="l",lty = 3, lwd = 3, col = "blue",ylim = c(0,8),
main = "Sensitivity calibration for each model",
xlab = expression(theta), ylab = "pdf")
## Plot Training
lines(x, P, type="l",lty = 1, lwd = 3, col = "green")
lines(x, Pl, type="l",lty = 1, lwd = 3, col = "pink")
## Plot Likelihood
lines(x, L, type = 'l', lwd = 3,col="yellow")
lines(x, Ll, type = 'l', lwd = 3,col="grey")
## Plot posterior
lines(x, dbeta(x, k + 1, n - k + 1), lty = 3, lwd = 3, col = "red")
lines(x, dbeta(x, kl + 1, n - kl + 1), lty = 3, lwd = 3, col = "black")
legend("topright", c("Prior", "Training-Naïve","Training-Logistic", "Test-Naïve","Test-Logistic","Posterior-Naïve","Posterior-Logistic"), lty = c(3, 1, 1,3,3,3,3), lwd = 3, col = c("blue", "green", "pink","yellow","grey","red","black"))
## Specificity
### Training Sample -
kp = 187 #Naive
np = 201 # total
klp = 194 # Logistic
kpriori = 10 #Naive
npriori = 49 # total
### Test Sample
k = 69+kp #Naive
n = 73+np # total
kl = 72+klp # Logistic
## x-axis for plotting
numSteps = 200
x = seq(0, 1, 1 / numSteps)
## Likelihood (training) function
P = x^kp * (1 - x)^(np - kp)
Pl = x^klp * (1 - x)^(np - klp)
## Likelihood function (test sample)
L = x^k * (1 - x)^(n - k)
Ll = x^kl * (1 - x)^(n - kl)
## Normalization
L = L / sum(L) * numSteps
Ll = Ll / sum(Ll) * numSteps
P = P / sum(P) * numSteps
Pl = Pl / sum(Pl) * numSteps
## Plot prior sensitivity = prevalence = 0.2
plot(x, dbeta(x,40,2), type="l",lty = 3, lwd = 3, col = "blue",ylim = c(0,40), xlim=c(0.8,1),
main = "Specificity calibration for each model",
xlab = expression(theta), ylab = "probability density function")
## Plot Training
lines(x, P, type="l",lty = 1, lwd = 3, col = "green")
lines(x, Pl, type="l",lty = 1, lwd = 3, col = "pink")
## Plot Likelihood
lines(x, L, type = 'l', lwd = 3,col="yellow")
lines(x, Ll, type = 'l', lwd = 3,col="grey")
## Plot posterior
lines(x, dbeta(x, k + 1, n - k + 1), lty = 3, lwd = 3, col = "red")
lines(x, dbeta(x, kl + 1, n - kl + 1), lty = 3, lwd = 3, col = "black")
legend("topleft", c("Prior", "Training-Naïve","Training-Logistic", "Test-Naïve","Test-Logistic","Posterior-Naïve","Posterior-Logistic"), lty = c(3, 1, 1,3,3,3,3), lwd = 3, col = c("blue", "green", "pink","yellow","grey","red","black"))
Cost function
Using a classifier makes no sense without considering the cost function of each prediction error. The default cost function is false positive equals false negative. For this reason, we implemented all the algorithms in an android app where fp/fn cost is an essencial parameter.
Multiple comparison p-values adjustments
testnames=c("age(DxV)","gestational.age(DxV)","map.basal(DxV)","bupivacaine(DxV)","sufentanil(DxV)","fentanyl(DxV)","morphine(DxV)","lower.intraoperative.map(DxV)","pacu.nausea.intensity(DxV)","pacu.vomiting.number(DxV)","maximum.nausea.intensity(DxV)","vomiting.number.2.24h(DxV)","map.drop(DxV)","map.drop.proportion(DxV)","Cinetosis(DxV)","Nausea 1 Trimester(DxV)","Nausea 3 Trimester(DxV)","Intraop.Nausea(DxV)","Age<25(DxV)","Gest<38w(DxV)","Never.smoked(DxV)","Low.Bupi(DxV)","Cinetosis(PONV)","Nausea 1 Trimester(PONV)","Nausea 3 Trimester(PONV)","Intraop.Nausea(PONV)","Age<25(PONV)","Gest<38w(PONV)","Never.smoked(PONV)","Low.Bupi(PONV)","age(PONV)","gestational.age(PONV)","map.basal(PONV)","bupivacaine(PONV)","sufentanil(PONV)","fentanyl(PONV)","morphine(PONV)","lower.intraoperative.map(PONV)","map.drop(PONV)","map.drop.proportion(PONV)","Naïve Hoslem","Multiple Logistic Regression Hoslem","Apfel Hoslem")
testvalues=c(wilcoxtest,pvaluelist2,pvaluelist3,plist5,hoslem.test(test$ponv,m2$naive.p)[3]$`p.value`,hoslem.test(test$ponv,unname(predict(mlr1,test)))[3]$`p.value`,hoslem.test(test$ponv,test$apfelp)[3]$`p.value`)
tests=data.frame(tests=testnames,p=testvalues)
#x_html <- knitr::kable(data.frame(tests,p=p.adjust(testvalues,method="BH")), "html")
x_html <- knitr::kable(data.frame(tests,FDR_adjusted_p=p.adjust(testvalues,method="BH")), "pandoc")
kable_styling(kable_input = x_html, "striped", full_width = FALSE)## Warning in kable_styling(kable_input = x_html, "striped", full_width =
## FALSE): Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for
## details.| tests | p | FDR_adjusted_p |
|---|---|---|
| age(DxV) | 0.4929829 | 1.0000000 |
| gestational.age(DxV) | 0.8869211 | 1.0000000 |
| map.basal(DxV) | 0.9148928 | 1.0000000 |
| bupivacaine(DxV) | 0.9250634 | 1.0000000 |
| sufentanil(DxV) | 0.7042008 | 1.0000000 |
| fentanyl(DxV) | 0.7525312 | 1.0000000 |
| morphine(DxV) | 0.7094460 | 1.0000000 |
| lower.intraoperative.map(DxV) | 0.7897966 | 1.0000000 |
| pacu.nausea.intensity(DxV) | 0.6265914 | 1.0000000 |
| pacu.vomiting.number(DxV) | 0.5883437 | 1.0000000 |
| maximum.nausea.intensity(DxV) | 0.6383733 | 1.0000000 |
| vomiting.number.2.24h(DxV) | 0.7575587 | 1.0000000 |
| map.drop(DxV) | 0.9918572 | 1.0000000 |
| map.drop.proportion(DxV) | 0.9894648 | 1.0000000 |
| Cinetosis(DxV) | 0.8926667 | 1.0000000 |
| Nausea 1 Trimester(DxV) | 1.0000000 | 1.0000000 |
| Nausea 3 Trimester(DxV) | 0.6257416 | 1.0000000 |
| Intraop.Nausea(DxV) | 0.2783978 | 0.7041826 |
| Age<25(DxV) | 0.2345593 | 0.6303780 |
| Gest<38w(DxV) | 0.6971843 | 1.0000000 |
| Never.smoked(DxV) | 1.0000000 | 1.0000000 |
| Low.Bupi(DxV) | 0.7217297 | 1.0000000 |
| Cinetosis(PONV) | 0.0061468 | 0.0377589 |
| Nausea 1 Trimester(PONV) | 0.0006816 | 0.0070421 |
| Nausea 3 Trimester(PONV) | 0.1452093 | 0.4460001 |
| Intraop.Nausea(PONV) | 0.0000000 | 0.0000001 |
| Age<25(PONV) | 0.0009826 | 0.0070421 |
| Gest<38w(PONV) | 0.0366478 | 0.1483475 |
| Never.smoked(PONV) | 0.0402170 | 0.1483475 |
| Low.Bupi(PONV) | 0.0006875 | 0.0070421 |
| age(PONV) | 0.0008506 | 0.0070421 |
| gestational.age(PONV) | 0.4448352 | 1.0000000 |
| map.basal(PONV) | 0.1557567 | 0.4465025 |
| bupivacaine(PONV) | 0.0413993 | 0.1483475 |
| sufentanil(PONV) | 0.8228957 | 1.0000000 |
| fentanyl(PONV) | 0.6278850 | 1.0000000 |
| morphine(PONV) | 0.0486619 | 0.1609587 |
| lower.intraoperative.map(PONV) | 0.3053587 | 0.7294680 |
| map.drop(PONV) | 0.0086703 | 0.0466027 |
| map.drop.proportion(PONV) | 0.0186188 | 0.0889566 |
| Naïve Hoslem | 0.9843251 | 1.0000000 |
| Multiple Logistic Regression Hoslem | 1.0000000 | 1.0000000 |
| Apfel Hoslem | 0.0000000 | 0.0000000 |
bdbn=data.frame(neversmoked=train\(neversmoked,cinetosis=train\)cinetosis,nausea1trim=train\(nausea1trim,ionv=train\)ionv)
bdbn\(younger=as.factor(train\)younger) bdbn\(neversmoked=as.factor(bdbn\)neversmoked) bdbn\(ionv=as.factor(bdbn\)ionv) bdbn\(cinetosis=as.factor(bdbn\)cinetosis) bdbn\(nausea1trim=as.factor(bdbn\)nausea1trim) bdbn\(preterm=as.factor(train\)preterm) bdbn\(ephedrine=as.factor(train\)received.ephedrine) bdbn\(hypotension=as.factor(train\)map.drop.proportion>0.1) bdbn\(previousponv.and.surgery=as.factor(train\)previous.surgery.and.ponv) bdbn\(lowbupi=as.factor(train\)lowbupi) #bdbn # Blacklist i=1; from=“”;to=“” n=(length(names(bdbn))) ## Age cannot be caused by any other variable for(i in c(1:n)){ from[i]=“younger”;to[i]=names(bdbn)[i]} #to[i+n]=“younger”;from[i+n]=names(bdbn)[i]
from[n+1]=“ephedrine”;to[n+1]=“hypotension” from[n+2]=“ephedrine”;to[n+2]=“neversmoked”
blacklist=data.frame(from=from,to=to)
plot(iamb(bdbn,blacklist=blacklist))