# Precision diagnostica para detectar lesiones periapicales entre panoramica y periapical con CBCT como gold standard
# acuerdo entre observadores
df <- read.csv("precisiondiagnosticapanoperiapicalcbct.csv", header = TRUE, sep = ",")
df <- subset(df, select = c(1:10)) #hay dos columnas extras, las remuevo
df<- df[complete.cases(df),] # omito todos los NA
#veo los datos
str(df)
## 'data.frame': 42 obs. of 10 variables:
## $ caso : int 1 2 3 4 5 6 7 8 9 10 ...
## $ CBCT : int 1 1 1 1 0 1 1 1 1 0 ...
## $ pano1.obs1: int 1 1 1 1 0 1 0 1 1 0 ...
## $ pano1.obs2: int 1 1 1 1 0 1 1 0 1 0 ...
## $ pa1.obs1 : int 1 1 1 1 0 1 0 1 1 0 ...
## $ pa1.obs2 : int 1 1 1 0 0 1 0 1 1 0 ...
## $ pano2.obs1: int 1 1 1 1 0 1 0 0 1 0 ...
## $ pano2.obs2: int 1 1 1 1 0 1 1 0 1 0 ...
## $ pa2.obs1 : int 1 1 1 1 0 1 1 1 1 0 ...
## $ pa2.obs2 : int 1 1 1 0 0 1 0 1 1 0 ...
head(df)
## caso CBCT pano1.obs1 pano1.obs2 pa1.obs1 pa1.obs2 pano2.obs1 pano2.obs2
## 1 1 1 1 1 1 1 1 1
## 2 2 1 1 1 1 1 1 1
## 3 3 1 1 1 1 1 1 1
## 4 4 1 1 1 1 0 1 1
## 5 5 0 0 0 0 0 0 0
## 6 6 1 1 1 1 1 1 1
## pa2.obs1 pa2.obs2
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 0
## 5 0 0
## 6 1 1
df
## caso CBCT pano1.obs1 pano1.obs2 pa1.obs1 pa1.obs2 pano2.obs1 pano2.obs2
## 1 1 1 1 1 1 1 1 1
## 2 2 1 1 1 1 1 1 1
## 3 3 1 1 1 1 1 1 1
## 4 4 1 1 1 1 0 1 1
## 5 5 0 0 0 0 0 0 0
## 6 6 1 1 1 1 1 1 1
## 7 7 1 0 1 0 0 0 1
## 8 8 1 1 0 1 1 0 0
## 9 9 1 1 1 1 1 1 1
## 10 10 0 0 0 0 0 0 0
## 11 11 1 1 0 0 0 1 0
## 12 12 1 1 1 1 1 1 1
## 13 13 1 1 0 0 0 0 0
## 14 14 1 1 1 1 0 1 1
## 15 15 0 0 0 0 0 0 0
## 16 16 1 0 1 0 0 0 0
## 17 17 1 0 0 0 0 0 0
## 18 18 1 1 1 0 0 1 1
## 19 19 1 1 1 0 0 0 1
## 20 20 0 1 0 0 0 0 0
## 21 21 1 0 0 1 1 0 0
## 22 22 1 1 1 1 1 1 1
## 23 23 1 1 0 1 1 1 1
## 24 24 1 1 0 1 1 1 1
## 25 25 0 0 0 0 0 0 0
## 26 26 1 0 0 0 0 0 0
## 27 27 0 0 1 0 0 0 1
## 28 28 1 0 0 0 0 0 0
## 29 29 1 1 1 1 1 1 1
## 30 30 1 1 0 1 0 1 0
## 31 31 0 0 0 0 0 0 0
## 32 32 1 0 1 1 0 0 0
## 33 33 1 0 0 0 0 0 0
## 34 34 1 1 1 1 1 0 1
## 35 35 1 0 1 0 0 0 0
## 36 36 0 0 0 0 0 0 0
## 37 37 1 0 1 0 0 1 1
## 38 38 0 1 1 0 0 0 1
## 39 39 1 1 1 1 1 1 0
## 40 40 1 1 1 1 1 1 1
## 41 41 1 1 1 0 1 0 1
## 42 42 0 1 0 0 0 0 0
## pa2.obs1 pa2.obs2
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 0
## 5 0 0
## 6 1 1
## 7 1 0
## 8 1 1
## 9 1 1
## 10 0 0
## 11 1 0
## 12 1 1
## 13 1 0
## 14 1 1
## 15 0 1
## 16 0 0
## 17 1 0
## 18 0 1
## 19 0 0
## 20 0 0
## 21 1 0
## 22 1 1
## 23 1 1
## 24 1 1
## 25 0 0
## 26 0 0
## 27 0 0
## 28 0 0
## 29 1 1
## 30 1 0
## 31 1 0
## 32 0 0
## 33 1 0
## 34 1 1
## 35 1 0
## 36 0 0
## 37 0 0
## 38 0 0
## 39 1 1
## 40 1 1
## 41 1 0
## 42 0 0
# para el kappa
install.packages("irr")
## Error in install.packages : Updating loaded packages
library(irr)
# hago las matrices
# para las radiografÃas
pano1<-as.matrix(df[,c(3,4)])
pano2<-as.matrix(df[,c(7,8)])
pa1<-as.matrix(df[,c(5,6)])
pa2<-as.matrix(df[,c(9,10)])
# para los observadores
obs1pano<-as.matrix(df[,c(3,7)])
obs1pa<-as.matrix(df[,c(5,9)])
obs2pano<-as.matrix(df[,c(4,8)])
obs2pa<-as.matrix(df[,c(6,10)])
####### kappa inter
#####################
# Kappa inter para los examenes
# pano
KappaPano1Inter <- kappa2(pano1, "unweighted")
KappaPano2Inter <- kappa2(pano2, "unweighted")
# periapical
KappaPa1Inter <- kappa2(pa1, "unweighted")
KappaPa2Inter <- kappa2(pa2, "unweighted")
####### kappa intra
#####################
# Kappa inter para los observadores
# 1
KappaObs1Pano <- kappa2(obs1pano, "unweighted")
KappaObs1Pa <- kappa2(obs1pa, "unweighted")
# 2
KappaObs2Pano <- kappa2(obs2pano, "unweighted")
KappaObs2Pa <- kappa2(obs2pa, "unweighted")
#reportes de kappa
# kappa inter para pano1
KappaPano1Inter # 0.321
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.321
##
## z = 2.09
## p-value = 0.0366
#kappa inter para pano 2
KappaPano2Inter # 0.571
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.571
##
## z = 3.74
## p-value = 0.000183
# kappa inter para periapical 1
KappaPa1Inter # 0.756
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.756
##
## z = 4.96
## p-value = 7.23e-07
# kappa inter para periapical 2
KappaPa2Inter # 0.408
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.408
##
## z = 2.9
## p-value = 0.00376
# kappa intra para obs 1 para pano
KappaObs1Pano # 0.583
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.583
##
## z = 3.99
## p-value = 6.53e-05
# kappa intra para obs 2 para pano
KappaObs2Pano # 0.714
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.714
##
## z = 4.65
## p-value = 3.32e-06
# kappa intra para obs 1 para periapicales
KappaObs1Pa # 0.581
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.581
##
## z = 3.98
## p-value = 6.82e-05
# kappa intra para obs 1 para periapicales
KappaObs2Pa # 0.751
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 42
## Raters = 2
## Kappa = 0.751
##
## z = 4.87
## p-value = 1.11e-06
### Precisión diagnóstica
######################
install.packages("caret")
## Error in install.packages : Updating loaded packages
library(caret)
str(df)
## 'data.frame': 42 obs. of 10 variables:
## $ caso : int 1 2 3 4 5 6 7 8 9 10 ...
## $ CBCT : int 1 1 1 1 0 1 1 1 1 0 ...
## $ pano1.obs1: int 1 1 1 1 0 1 0 1 1 0 ...
## $ pano1.obs2: int 1 1 1 1 0 1 1 0 1 0 ...
## $ pa1.obs1 : int 1 1 1 1 0 1 0 1 1 0 ...
## $ pa1.obs2 : int 1 1 1 0 0 1 0 1 1 0 ...
## $ pano2.obs1: int 1 1 1 1 0 1 0 0 1 0 ...
## $ pano2.obs2: int 1 1 1 1 0 1 1 0 1 0 ...
## $ pa2.obs1 : int 1 1 1 1 0 1 1 1 1 0 ...
## $ pa2.obs2 : int 1 1 1 0 0 1 0 1 1 0 ...
# primero, como son cuatro pano (dos del obs1 y dos del obs2)
# voy a crear un CBCT x 4
cbct4 <- rep(df$CBCT, 4)
cbct4
## [1] 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1
## [36] 0 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 1
## [71] 1 1 0 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1
## [106] 1 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 0 1 1 1 1
## [141] 0 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1 0 1 0 1 1 1 0
# ahora combino todas las mediciones de cada test
pano <- c(df$pano1.obs1, df$pano1.obs2, df$pano2.obs1, df$pano2.obs2)
pa <- c(df$pa1.obs1, df$pa1.obs2, df$pa2.obs1, df$pa2.obs2)
# hago una tabla para test diagnostico
cbctXpano <- with(df, table(pano, cbct4))
cbctXpa <- with(df, table(pa, cbct4))
cbctXpano
## cbct4
## pano 0 1
## 0 33 48
## 1 7 80
cbctXpa
## cbct4
## pa 0 1
## 0 38 52
## 1 2 76
# Evaluar precisión diagnostica
Panoramica <- confusionMatrix(pano, cbct4, positive="1")
Periapical <- confusionMatrix(pa, cbct4, positive="1")
#Precisión diagnostica
# Panorámicas
Panoramica
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 33 48
## 1 7 80
##
## Accuracy : 0.6726
## 95% CI : (0.5961, 0.7429)
## No Information Rate : 0.7619
## P-Value [Acc > NIR] : 0.9967
##
## Kappa : 0.3328
## Mcnemar's Test P-Value : 6.906e-08
##
## Sensitivity : 0.6250
## Specificity : 0.8250
## Pos Pred Value : 0.9195
## Neg Pred Value : 0.4074
## Prevalence : 0.7619
## Detection Rate : 0.4762
## Detection Prevalence : 0.5179
## Balanced Accuracy : 0.7250
##
## 'Positive' Class : 1
##
# ver en http://araw.mede.uic.edu/cgi-bin/testcalc.pl?DT=80&Dt=48&dT=7&dt=33&2x2=Compute
# Periapical
Periapical
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 38 52
## 1 2 76
##
## Accuracy : 0.6786
## 95% CI : (0.6023, 0.7484)
## No Information Rate : 0.7619
## P-Value [Acc > NIR] : 0.9946
##
## Kappa : 0.3803
## Mcnemar's Test P-Value : 2.592e-11
##
## Sensitivity : 0.5938
## Specificity : 0.9500
## Pos Pred Value : 0.9744
## Neg Pred Value : 0.4222
## Prevalence : 0.7619
## Detection Rate : 0.4524
## Detection Prevalence : 0.4643
## Balanced Accuracy : 0.7719
##
## 'Positive' Class : 1
##
# ver en http://araw.mede.uic.edu/cgi-bin/testcalc.pl?DT=76&Dt=52&dT=1&dt=38&2x2=Compute
# elijo solo pano 2 por su kappa de 0.57 y periapical 1 por kappa 0.756
# hago un cbct2 de 84
cbct2 <- rep(df$CBCT, 2)
mejorpano <- c(df$pano2.obs1, df$pano2.obs2)
mejorperiapical <- c(df$pa1.obs1, df$pa1.obs2)
MejorPanoramica <- confusionMatrix(mejorpano, cbct2, positive="1")
MejorPeriapical <- confusionMatrix(mejorperiapical, cbct2, positive="1")
# precisión diagnostica utilizando las pano y pa con mejor kappa (pano 2 y pa1)
# precisión diagnostica de pano2
MejorPanoramica
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 18 27
## 1 2 37
##
## Accuracy : 0.6548
## 95% CI : (0.5431, 0.7552)
## No Information Rate : 0.7619
## P-Value [Acc > NIR] : 0.9905
##
## Kappa : 0.3344
## Mcnemar's Test P-Value : 8.324e-06
##
## Sensitivity : 0.5781
## Specificity : 0.9000
## Pos Pred Value : 0.9487
## Neg Pred Value : 0.4000
## Prevalence : 0.7619
## Detection Rate : 0.4405
## Detection Prevalence : 0.4643
## Balanced Accuracy : 0.7391
##
## 'Positive' Class : 1
##
# precisión diagnostica de pa1
MejorPeriapical
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 20 29
## 1 0 35
##
## Accuracy : 0.6548
## 95% CI : (0.5431, 0.7552)
## No Information Rate : 0.7619
## P-Value [Acc > NIR] : 0.9905
##
## Kappa : 0.365
## Mcnemar's Test P-Value : 1.999e-07
##
## Sensitivity : 0.5469
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.4082
## Prevalence : 0.7619
## Detection Rate : 0.4167
## Detection Prevalence : 0.4167
## Balanced Accuracy : 0.7734
##
## 'Positive' Class : 1
##
install.packages("pROC")
## Error in install.packages : Updating loaded packages
install.packages("doParallel")
## Error in install.packages : Updating loaded packages
library(pROC)
library(doParallel)
table(pano)
## pano
## 0 1
## 81 87
table(cbct4)
## cbct4
## 0 1
## 40 128
rocPano <- plot.roc(cbct4, pano , percent = TRUE, main="Smoothing")
rocPeriapical <- plot.roc(cbct4, pa , percent = TRUE, main="Smoothing")
# en un solo gráfico
# calculo las áreas bajo la curva
AUCPano <- auc(cbct4, pano)
AUCPano
## Area under the curve: 0.725
AUCPa <- auc(cbct4, pa)
AUCPa
## Area under the curve: 0.7719
# hago el gráfico
par(mfrow=c(1,1)) #n?mero de figuras en el plot;c(filas,columnas)
par(mar=c(5,5,3,2)+0.1) #margenes
rocobj1 <- plot.roc(cbct4, pano , main="Comparación de curvas ROC", percent=TRUE, col="#1c61b6",
xlab="Especificidad (%)", ylab="Sensibilidad (%)",
ylim=c(0,100))
rocobj2 <- lines.roc(cbct4, pa , percent=TRUE, col="#008600")
testobj <- roc.test(rocobj1, rocobj2)
text(50, 50, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0, .5))
text(70, 66, labels=paste("AUC Panorámica =", AUCPano), adj=c(0, .5), col="#1c61b6")
text(90, 79, labels=paste("AUC Periapical =", AUCPa), adj=c(0, .5), col="#008600")
legend("bottomright", legend=c("Panorámica", "Periapical"), col=c("#1c61b6", "#008600"), lwd=2)
