DiagnosisMed friends,<br><br>Finally we are getting some requests about the functions available in the package. Therefore, It would be nice if we get used to share this informations in the diagnosismed email list. <br><br>
As in earlier emails the ROC function in the packages has a error in the AUC confidence limist. It did come to my attention in late december. This week two emails from users came quastioning about the AUC and its confidence limists. Tura proposed a fix. However the fix gave me some errors when I runned the function. There goes...<br>
<br><span style="font-family: courier new,monospace;"># gold is the vector meneaning the reference standard e test is the test under evaluarion</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> if (is.numeric(gold)==TRUE){</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> # Making test summary, overall, disease only, without disease only</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> test.summary<-round(c(summary(test),sd(test)),digits=5)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> test.summary<-rbind(test.summary,round(c(summary(test[gold==0]),sd(test[gold==0])),digits=5))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> test.summary<-rbind(test.summary,round(c(summary(test[gold==1]),sd(test[gold==1])),digits=5))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> colnames(test.summary)<-c("Min.","1st Qu.","Median","Mean","3rd Qu.","Max.","SD")</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> rownames(test.summary)<-c("Overall summary","Without disease", "With disease")</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> # Estimating the AUC and confidence limits, inspired in auc{PresenceAbsence}</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> # m length(test[gold == 0])</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> # n length(test[gold == 1])</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> AUC <- ((as.double(length(test[gold == 0]))) * (as.double(length(test[gold ==</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> 1]))) + ((as.double(length(test[gold == 0]))) * ((as.double(length(</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> test[gold == 0]))) + 1))/2 - sum(rank(test,ties.method = "average")[</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> gold == 0]))/((as.double(length(test[gold == 0]))) * (as.double(</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> length(test[gold == 1]))))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> AUC[AUC < 0.5] <- 1 - AUC</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> X<-test[gold==1] # pag 209</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> Y<-test[gold==1] # pag 209</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># The same as before just with different codes in the variables </span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> if (is.factor(gold)==TRUE){</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> # The same tests summary but with different reference standard codes</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> test.summary<-round(c(summary(test),sd(test)),digits=5)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> test.summary<-rbind(test.summary,round(c(summary(test[gold=="negative"]),sd(test[gold=="negative"])),digits=5))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> test.summary<-rbind(test.summary,round(c(summary(test[gold=="positive"]),sd(test[gold=="positive"])),digits=5))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> colnames(test.summary)<-c("Min.","1st Qu.","Median","Mean","3rd Qu.","Max.","SD")</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> rownames(test.summary)<-c("Overall summary","Without disease", "With disease")</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> AUC <- ((as.double(length(test[gold =="negative"]))) * (as.double(length(test[gold ==</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> "positive"]))) + ((as.double(length(test[gold =="negative"]))) * ((as.double(length(</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> test[gold =="negative"]))) + 1))/2 - sum(rank(test,ties.method = "average")[</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> gold =="negative"]))/((as.double(length(test[gold =="negative"]))) * (as.double(</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> length(test[gold == "positive"]))))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> AUC[AUC < 0.5] <- 1 - AUC</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> X<-test[gold=="negative"] # pag 209</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> Y<-test[gold=="positive"] # pag 209</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># Here is the fix </span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> m<-length(X) # pag 209</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> n<-length(Y) # pag 209</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> D10X<-function(Xi){(1/n)*sum(Y>=Xi)} # pag 211</span> <u><b> I believe the error is here</b></u> although not sure<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> DO1Y<-function(Yi){(1/m)*sum(Yi>=X)} # pag 211</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> VAR.AUC<-sum((tapply(X,X,"D10X")-AUC)^2)/(m*(m-1))+sum((tapply(Y,Y,"D10X")-AUC)^2)/(n*(n-1)) # pag 211</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> SD.AUC<-sqrt(VAR.AUC)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> alpha<-1-CL</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> AUC.summary<-c(AUC- qnorm(1-alpha/2)*SD.AUC,AUC,AUC+ qnorm(1-alpha/2)*SD.AUC)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> names(AUC.summary)<-c("AUC inf conf limit", "AUC","AUC sup conf limit")</span><br style="font-family: courier new,monospace;" clear="all">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"># here is the output and the error below. the graph comes OK... </span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">> ROC(rocdata$Gold,rocdata$test1)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$pop.prevalence</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Informed disease prevalence - same as sample prevalence if not informed </span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> 0.3445946 </span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$sample.size</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">[1] 148</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$sample.prevalence</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">Observed prevalence by gold standard </span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> 0.3445946 </span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$test.summary</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> Min. 1st Qu. Median Mean 3rd Qu. Max. SD</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">Overall summary 0.056 0.102 0.1315 0.3835 0.6312 1.914 0.39349</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Without disease 0.056 0.098 0.1110 0.1590 0.1330 0.837 0.15833</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">With disease 0.098 0.557 0.7950 0.8106 1.0490 1.914 0.35106</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$AUC.summary</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">AUC inf conf limit AUC AUC sup conf limit </span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> 0.7545496 0.9540125 1.1534755 </span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">$test.best.cutoff</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> best.cutoff Sensitivity <a href="http://Se.inf.cl">Se.inf.cl</a> <a href="http://Se.sup.cl">Se.sup.cl</a> Specificity <a href="http://Sp.inf.cl">Sp.inf.cl</a> <a href="http://Sp.sup.cl">Sp.sup.cl</a> PLR <a href="http://PLR.inf.cl">PLR.inf.cl</a> <a href="http://PLR.sup.cl">PLR.sup.cl</a></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. Accuracy 0.365 0.9216 0.8150 0.9691 0.9278 0.8585 0.9646 12.76 4.74 34.36</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. DOR 0.275 0.9608 0.8678 0.9892 0.8969 0.8205 0.9430 9.32 4.13 21.05</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Min. Error rate 0.365 0.9216 0.8150 0.9691 0.9278 0.8585 0.9646 12.76 4.74 34.36</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. Accuracy area 0.275 0.9608 0.8678 0.9892 0.8969 0.8205 0.9430 9.32 4.13 21.05</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. Sens+Spec 0.275 0.9608 0.8678 0.9892 0.8969 0.8205 0.9430 9.32 4.13 21.05</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. Youden 0.275 0.9608 0.8678 0.9892 0.8969 0.8205 0.9430 9.32 4.13 21.05</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Se=Sp 0.356 0.9216 0.8150 0.9691 0.9175 0.8456 0.9576 11.17 4.44 28.08</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Min. ROC distance 0.365 0.9216 0.8150 0.9691 0.9278 0.8585 0.9646 12.76 4.74 34.36</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Max. Efficiency 0.365 0.9216 0.8150 0.9691 0.9278 0.8585 0.9646 12.76 4.74 34.36</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">Min. MCT 0.365 0.9216 0.8150 0.9691 0.9278 0.8585 0.9646 12.76 4.74 34.36</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">Warning messages:</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">1: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">2: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">3: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">4: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">5: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">6: In Y >= Xi :</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> longer object length is not a multiple of shorter object length</span><br><br># Besides I would expect the interval to come closer to something like 0,907 and 0,982 as some other analsys I did way before and comparing with some other software. It seems weird AUC conf limist wider than the upper limit 1. <br>
<br><br><br>-- <br>Abraço forte e que a força esteja com você,<br>Pedro Emmanuel<br>