[Diagnosismed-commits] r16 - in pkg/DiagnosisMed: . R chm html man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 18 05:43:48 CEST 2009
Author: pedrobrasil
Date: 2009-06-18 05:43:46 +0200 (Thu, 18 Jun 2009)
New Revision: 16
Added:
pkg/DiagnosisMed/R/plot.TGROC.r
Removed:
pkg/DiagnosisMed/chm/00Index.html
pkg/DiagnosisMed/chm/DiagnosisMed.hhp
pkg/DiagnosisMed/chm/DiagnosisMed.toc
pkg/DiagnosisMed/chm/ROC.html
pkg/DiagnosisMed/chm/Rchm.css
pkg/DiagnosisMed/chm/diagnosis.html
pkg/DiagnosisMed/chm/diagtest.chm
pkg/DiagnosisMed/chm/diagtest.hhp
pkg/DiagnosisMed/chm/diagtest.toc
pkg/DiagnosisMed/chm/interact.ROC.html
pkg/DiagnosisMed/chm/logo.jpg
pkg/DiagnosisMed/chm/plot.diag.html
pkg/DiagnosisMed/chm/rocdata.html
pkg/DiagnosisMed/chm/tutorial.html
pkg/DiagnosisMed/html/ROC.html
pkg/DiagnosisMed/html/diagtest.HTML
pkg/DiagnosisMed/html/plot.diag.HTML
pkg/DiagnosisMed/html/print.diag.HTML
Modified:
pkg/DiagnosisMed/DESCRIPTION
pkg/DiagnosisMed/R/ROC.r
pkg/DiagnosisMed/R/TGROC.r
pkg/DiagnosisMed/R/diagnosis.r
pkg/DiagnosisMed/R/diagnosisI.r
pkg/DiagnosisMed/R/plot.ROC.r
pkg/DiagnosisMed/R/print.ROC.r
pkg/DiagnosisMed/man/LRgrgaph.Rd
pkg/DiagnosisMed/man/ROC.Rd
pkg/DiagnosisMed/man/TGROC.Rd
pkg/DiagnosisMed/man/diagnosis.Rd
pkg/DiagnosisMed/man/interact.ROC.Rd
pkg/DiagnosisMed/man/plot.diag.Rd
Log:
Modified: pkg/DiagnosisMed/DESCRIPTION
===================================================================
--- pkg/DiagnosisMed/DESCRIPTION 2009-05-26 18:03:34 UTC (rev 15)
+++ pkg/DiagnosisMed/DESCRIPTION 2009-06-18 03:43:46 UTC (rev 16)
@@ -1,9 +1,9 @@
Package: DiagnosisMed
-Version: 0.2.2
-Date: 2009-05-26
+Version: 0.2.2.1
+Date: 2009-06-17
Author: Pedro Brasil <pedro.brasil at ipec.fiocruz.br>
Maintainer: Pedro Brasil <pedro.brasil at ipec.fiocruz.br>
-Depends: R (>= 2.7.2),epitools,epicalc, TeachingDemos, tcltk, AMORE
+Depends: R (>= 2.7.2),epitools, epicalc, TeachingDemos, tcltk, AMORE
Title: Diagnostic test accuracy evaluation for medical professionals.
Description: DiagnosisMed is a package to analyze data from diagnostic test accuracy evaluating health conditions. It is being built to be used by health professionals. This package is - or soon will be - able to estimate sample size for common situations in diagnostic test accuracy, estimate sensitivity and specificity from categorical and continuous test results including some evaluations of indeterminate results, or compare different analysis strategies into measures commonly used by health professionals.
License: GPL (>= 2)
Modified: pkg/DiagnosisMed/R/ROC.r
===================================================================
--- pkg/DiagnosisMed/R/ROC.r 2009-05-26 18:03:34 UTC (rev 15)
+++ pkg/DiagnosisMed/R/ROC.r 2009-06-18 03:43:46 UTC (rev 16)
@@ -55,15 +55,15 @@
rownames(test.summary)<-c("Overall summary","Without disease", "With disease")
AUC <- ((as.double(length(test[gold =="negative"]))) * (as.double(length(test[gold =="positive"]))) + ((as.double(length(test[gold =="negative"]))) * ((as.double(length(test[gold =="negative"]))) + 1))/2 - sum(rank(test,ties.method = "average")[gold =="negative"]))/((as.double(length(test[gold =="negative"]))) * (as.double(length(test[gold == "positive"]))))
AUC[AUC < 0.5] <- 1 - AUC
- X<-test[gold=="negative"] # pag 209
- Y<-test[gold=="positive"] # pag 209
+ X<-test[gold=="negative"]
+ Y<-test[gold=="positive"]
}
- m<-length(X) # pag 209
- n<-length(Y) # pag 209
- D10X<-function(Xi){(1/n)*sum(Y>=Xi)} # pag 211
- D01Y<-function(Yi){(1/m)*sum(Yi>=X)} # pag 211
- VAR.AUC<-sum((tapply(X,X,"D10X")-AUC)^2)/(m*(m-1))+sum((tapply(Y,Y,"D01Y")-AUC)^2)/(n*(n-1)) # pag 211
+ m<-length(X)
+ n<-length(Y)
+ D10X<-function(Xi){(1/n)*sum(Y>=Xi)}
+ D01Y<-function(Yi){(1/m)*sum(Yi>=X)}
+ VAR.AUC<-sum((tapply(X,X,"D10X")-AUC)^2)/(m*(m-1))+sum((tapply(Y,Y,"D01Y")-AUC)^2)/(n*(n-1))
SD.AUC<-sqrt(VAR.AUC)
alpha<-1-CL
AUC.summary<-c(AUC- qnorm(1-alpha/2)*SD.AUC,AUC,AUC+ qnorm(1-alpha/2)*SD.AUC)
@@ -73,88 +73,68 @@
#FP sum(test.table[i:nrow(test.table),1])
#TN sum(test.table[1:i-1,1])
#FN sum(test.table[1:i-1,2])
+ D<-sum(test.table[,2])
+ ND<-sum(test.table[,1])
# Taking the rownames of the test.table to be results first column
test.values<-(as.numeric(rownames(unclass(test.table))))
test.diag.table<-as.data.frame(test.values)
# Making a table with Se Sp PLR NLR PPV NPV and its confidence limits for each cut-off
for (i in 1:nrow(test.diag.table)) {
- test.diag.table$Sensitivity[i] <- round(sum(test.table[(i:nrow(test.table)),2])/sum(test.table[,2]),digits=4)
- test.diag.table$Se.inf.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),sum(test.table[,2]),conf.level=CL)[4]),digits=4)
- test.diag.table$Se.sup.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),sum(test.table[,2]),conf.level=CL)[5]),digits=4)
- test.diag.table$Specificity[i] <- round((sum(test.table[(1:i-1),1]))/(sum(test.table[,1])),digits=4)
- test.diag.table$Sp.inf.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[(1:i-1),1]),sum(test.table[,1]),conf.level=CL)[4]),digits=4)
- test.diag.table$Sp.sup.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[(1:i-1),1]),sum(test.table[,1]),conf.level=CL)[5]),digits=4)
- test.diag.table$PPV[i]<-round((sum(test.table[i:nrow(test.table),2]))/((sum(test.table[i:nrow(test.table),2]))+(sum(test.table[i:nrow(test.table),1]))),digits=4)
- test.diag.table$PPV.inf.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),((sum(test.table[i:nrow(test.table),2]))+(sum(
- test.table[i:nrow(test.table),1]))),conf.level=CL)[4]),digits=4)
- test.diag.table$PPV.sup.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),((sum(test.table[i:nrow(test.table),2]))+(sum(
- test.table[i:nrow(test.table),1]))),conf.level=CL)[5]),digits=4)
- test.diag.table$NPV[i]<-round((sum(test.table[1:i-1,1]))/((sum(test.table[1:i-1,1]))+sum(test.table[1:i-1,2])),digits=4)
- test.diag.table$NPV.inf.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[1:i-1,1]),((sum(test.table[1:i-1,1]))+(sum(test.table[1:i-1,])))
- ,conf.level=CL)[4]),digits=4)
- test.diag.table$NPV.sup.cl[i]<-round(as.numeric(binom.wilson(sum(test.table[1:i-1,1]),((sum(test.table[1:i-1,1]))+(sum(test.table[1:i-1,])))
- ,conf.level=CL)[5]),digits=4)
- }
+ test.diag.table$TP[i] <- sum(test.table[i:nrow(test.table),2])
+ test.diag.table$FN[i] <- sum(test.table[1:i-1,2])
+ test.diag.table$FP[i] <- sum(test.table[i:nrow(test.table),1])
+ test.diag.table$TN[i] <- sum(test.table[1:i-1,1])
+ }
+ test.diag.table$Sensitivity <- round(test.diag.table$TP/D,digits=4)
+ test.diag.table$Se.inf.cl <- round(binom.wilson(test.diag.table$TP,D,conf.level=CL)[4]$lower,digits=4)
+ test.diag.table$Se.sup.cl <- round(binom.wilson(test.diag.table$TP,D,conf.level=CL)[5]$upper,digits=4)
+ test.diag.table$Specificity <- round(test.diag.table$TN/ND,digits=4)
+ test.diag.table$Sp.inf.cl <- round(binom.wilson(test.diag.table$TN,ND,conf.level=CL)[4]$lower,digits=4)
+ test.diag.table$Sp.sup.cl <- round(binom.wilson(test.diag.table$TN,ND,conf.level=CL)[5]$upper,digits=4)
+ test.diag.table$PPV <- round(test.diag.table$TP/(test.diag.table$TP + test.diag.table$FP),digits=4)
+ test.diag.table$PPV.inf.cl <- round(binom.wilson(test.diag.table$TP,(test.diag.table$TP + test.diag.table$TP),conf.level=CL)[4]$lower,digits=4)
+ test.diag.table$PPV.sup.cl <- round(binom.wilson(test.diag.table$TP,(test.diag.table$TP + test.diag.table$FN),conf.level=CL)[5]$upper,digits=4)
+ test.diag.table$NPV <- round(test.diag.table$TN/(test.diag.table$TN + test.diag.table$FN),digits=4)
+ test.diag.table$NPV.inf.cl <- round(binom.wilson(test.diag.table$TN,(test.diag.table$TN + test.diag.table$FN),conf.level=CL)[4]$lower,digits=4)
+ test.diag.table$NPV.sup.cl <- round(binom.wilson(test.diag.table$TN,(test.diag.table$TN + test.diag.table$FN),conf.level=CL)[5]$upper,digits=4)
test.diag.table$PLR<-round(test.diag.table$Sensitivity/(1-test.diag.table$Specificity),digits=2)
- test.diag.table$PLR.inf.cl<-round(exp(log(test.diag.table$PLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-test.diag.table$Sensitivity)/
- ((sum(test.diag.table[,2]))*test.diag.table$Specificity)+(test.diag.table$Specificity)/((sum(test.diag.table[,1]))*(1-test.diag.table$Specificity)))),digits=2)
- test.diag.table$PLR.sup.cl<-round(exp(log(test.diag.table$PLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-test.diag.table$Sensitivity)/((sum(
- test.diag.table[,2]))*test.diag.table$Specificity)+(test.diag.table$Specificity)/((sum(test.diag.table[,1]))*(1-test.diag.table$Specificity)))),digits=2)
+ test.diag.table$PLR.inf.cl<-round(exp(log(test.diag.table$PLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-test.diag.table$Sensitivity)/(
+ (D)*test.diag.table$Specificity)+(test.diag.table$Specificity)/((ND)*(1-test.diag.table$Specificity)))),digits=2)
+ test.diag.table$PLR.sup.cl<-round(exp(log(test.diag.table$PLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-test.diag.table$Sensitivity)/(
+ (D)*test.diag.table$Specificity)+(test.diag.table$Specificity)/((ND)*(1-test.diag.table$Specificity)))),digits=2)
test.diag.table$NLR<-round((1-test.diag.table$Sensitivity)/test.diag.table$Specificity,digits=2)
- test.diag.table$NLR.inf.cl<-round(exp(log(test.diag.table$NLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((test.diag.table$Sensitivity)/((sum(test.diag.table[,2]))*(1-test.diag.table$Sensitivity))+(1-test.diag.table$Specificity)/((sum(test.diag.table[,1]))*(test.diag.table$Specificity)))),digits=2)
- test.diag.table$NLR.sup.cl<-round(exp(log(test.diag.table$NLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((test.diag.table$Sensitivity)/((sum(test.diag.table[,2]))*(1-test.diag.table$Sensitivity))+(1-test.diag.table$Specificity)/((sum(test.diag.table[,1]))*(test.diag.table$Specificity)))),digits=2)
+ test.diag.table$NLR.inf.cl<-round(exp(log(test.diag.table$NLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((test.diag.table$Sensitivity)/((D)*(1-test.diag.table$Sensitivity))+(1-test.diag.table$Specificity)/((ND)*(test.diag.table$Specificity)))),digits=2)
+ test.diag.table$NLR.sup.cl<-round(exp(log(test.diag.table$NLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((test.diag.table$Sensitivity)/((D)*(1-test.diag.table$Sensitivity))+(1-test.diag.table$Specificity)/((ND)*(test.diag.table$Specificity)))),digits=2)
+ test.diag.table$Accuracy <- (test.diag.table$TN + test.diag.table$TP)/sample.size
+ test.diag.table$DOR <- ((test.diag.table$TN)*(test.diag.table$TP))/((test.diag.table$FP)*(test.diag.table$FN))
+ test.diag.table$DOR<-ifelse(test.diag.table$DOR==Inf,NA,test.diag.table$DOR)
+ test.diag.table$Error.rate <- ((test.diag.table$FP)+(test.diag.table$FN))/sample.size
+ test.diag.table$Accuracy.area <- ((test.diag.table$TP)*(test.diag.table$TN))/(D*ND)
+ test.diag.table$Max.Se.Sp <- test.diag.table$Sensitivity + test.diag.table$Specificity
+ test.diag.table$Youden <- test.diag.table$Sensitivity + test.diag.table$Specificity - 1
+ test.diag.table$Se.equals.Sp <- abs(test.diag.table$Specificity-test.diag.table$Sensitivity)
+ test.diag.table$MinRocDist <- (test.diag.table$Specificity-1)^2+(1-test.diag.table$Sensitivity)^2
+ test.diag.table$Efficiency<-(test.diag.table$Sensitivity*(pop.prevalence))+((1-(pop.prevalence))*test.diag.table$Specificity)
+ test.diag.table$MCT<-(1-(pop.prevalence))*(1-test.diag.table$Specificity)+(cost*(pop.prevalence))*(1-test.diag.table$Sensitivity)
- # The following is not result but will help to find the best cut-offs
- test.cutoff.table<-as.data.frame(test.values)
- for(i in 1:nrow(test.cutoff.table)) {
- #Accuracy is (TP+TN)/sample.size
- test.cutoff.table$Accuracy[i]<-(sum(test.table[i:nrow(test.table),2])+sum(test.table[1:i-1,1]))/sample.size
- test.cutoff.table$DOR[i]<-(((sum(test.table[i:nrow(test.table),2])*(sum(test.table[1:i-1,1])))/
- ((sum(test.table[i:nrow(test.table),1]))*(sum(test.table[1:i-1,2])))))
- test.cutoff.table$Error.rate[i]<-((sum(test.table[1:i-1,2]))+(sum(test.table[i:nrow(test.table),1])))/sample.size
- #Error.trade.off - error out of limits, only zeros can be used with negative subscripts
- #Error.trade.off[i]<-(sum(test.table[1:i-1,2])-sum(test.table[1:i-2,2]))
- # (sum(test.table[i:nrow(test.table),1])-sum(test.table[i+
- # 1:nrow(test.table),1]))
- test.cutoff.table$Accuracy.area[i]<-round((sum(test.table[i:nrow(test.table),2])*sum(test.table[1:i-1,1]))/((sample.size-pop.prevalence)*pop.prevalence),digits=4)
- }
- test.cutoff.table$DOR<-ifelse(test.cutoff.table$DOR==Inf,NA,test.cutoff.table$DOR)
- test.cutoff.table$Max.Se.Sp<-(test.diag.table$Sensitivity + test.diag.table$Specificity)
- test.cutoff.table$Youden<-test.diag.table$Sensitivity+test.diag.table$Specificity - 1
- test.cutoff.table$Se.equals.Sp<-abs(test.diag.table$Specificity-test.diag.table$Sensitivity)
- test.cutoff.table$MinRocDist<-(test.diag.table$Specificity-1)^2+(1-test.diag.table$Sensitivity)^2
- # Efficiency= Se*prevalence+(1-prevalence)*Sp
- test.cutoff.table$Efficiency<-(test.diag.table$Sensitivity*(pop.prevalence))+((1-(pop.prevalence))*test.diag.table$Specificity)
- # MissclassificatioCostTerm(MCT)=(1-prevalence)(1-Sp)+r*prevalence(1-Se) - r=cost(FN)/cost(FP)
- test.cutoff.table$MCT<-(1-(pop.prevalence))*(1-test.diag.table$Specificity)+(Cost*(pop.prevalence))*(1-test.diag.table$Sensitivity)
# Making a table with the test result for each best cut-off and attaching validity measures
-
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$Accuracy)]
- test.best.cutoff<-test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)]
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$DOR)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.min(test.cutoff.table$Error.rate)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$Accuracy.area)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$Max.Se.Sp)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$Youden)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.min(test.cutoff.table$Se.equals.Sp)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.min(test.cutoff.table$MinRocDist)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.max(test.cutoff.table$Efficiency)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
- best.cutoff<-test.cutoff.table$test.values[which.min(test.cutoff.table$MCT)]
- test.best.cutoff<-rbind(test.best.cutoff,test.diag.table[test.diag.table$test.values==best.cutoff,c(1:7,14:16)])
+ test.best.cutoff <- as.data.frame(rbind(
+ test.diag.table[which.max(test.diag.table$Accuracy),c(1,6:11,18:20)],
+ test.diag.table[which.max(test.diag.table$DOR),c(1,6:11,18:20)],
+ test.diag.table[which.min(test.diag.table$Error.rate),c(1,6:11,18:20)],
+ test.diag.table[which.max(test.diag.table$Accuracy.area),c(1,6:11,18:20)],
+ test.diag.table[which.max(test.diag.table$Max.Se.Sp),c(1,6:11,18:20)],
+ test.diag.table[which.max(test.diag.table$Youden),c(1,6:11,18:20)],
+ test.diag.table[which.min(test.diag.table$Se.equals.Sp),c(1,6:11,18:20)],
+ test.diag.table[which.min(test.diag.table$MinRocDist),c(1,6:11,18:20)],
+ test.diag.table[which.max(test.diag.table$Efficiency),c(1,6:11,18:20)],
+ test.diag.table[which.min(test.diag.table$MCT),c(1,6:11,18:20)]
+ ))
rownames(test.best.cutoff)<- c("Max. Accuracy", "Max. DOR","Min. Error rate",
"Max. Accuracy area","Max. Sens+Spec","Max. Youden","Se=Sp","Min. ROC distance",
"Max. Efficiency", "Min. MCT")
- rm(best.cutoff)
-
+
#names(pop.prevalence)<-c("Informed disease prevalence - same as sample prevalence if not informed")
#names(sample.prevalence)<-c("Observed prevalence by gold standard")
reteval<-list(pop.prevalence=pop.prevalence,
@@ -166,8 +146,7 @@
test.best.cutoff=test.best.cutoff,
test.diag.table=test.diag.table,
CL=CL,
- cost=cost,
- test.cutoff.table=test.cutoff.table)
+ cost=cost)
class(reteval)<-"ROC"
if(Print==TRUE){
Modified: pkg/DiagnosisMed/R/TGROC.r
===================================================================
--- pkg/DiagnosisMed/R/TGROC.r 2009-05-26 18:03:34 UTC (rev 15)
+++ pkg/DiagnosisMed/R/TGROC.r 2009-06-18 03:43:46 UTC (rev 16)
@@ -39,12 +39,6 @@
if(is.null(precision)||!is.numeric(precision)){
stop("Precision must be set to a numeric value!")
}
- if(Plot!="Both" & Plot!="Non-parametric" & Plot!="Parametric" & Plot!="None"){
- stop("Plot must be set either to 'None','Both','Non-parametric' or 'Parametric'!")
- }
- if(Plot.cutoff!="Min.MCT" & Plot.cutoff!="Se=Sp" & Plot.cutoff!="Max.Efficiency" & Plot.cutoff!="None"){
- stop("Plot.cutoff must be set either to 'None','Max.Efficiency','Min.MCT' or 'Se=Sp'!")
- }
sample.prevalence<-(sum(test.table[,2]))/(sum(test.table))
if (Prevalence==0){
pop.prevalence<-sample.prevalence
@@ -64,28 +58,34 @@
conf.limit<-CL
inc<-Inconclusive
names(inc)<-"Inconclusive tolerance level"
+
+ D<-sum(test.table[,2])
+ ND<-sum(test.table[,1])
+
+ # Taking the rownames of the test.table to be results first column
test.values<-(as.numeric(rownames(unclass(test.table))))
non.parametric<-as.data.frame(test.values)
-
+ # Making a table with Se Sp PLR NLR PPV NPV and its confidence limits for each cut-off
for (i in 1:nrow(non.parametric)) {
- non.parametric$Sensitivity[i] <- round(sum(test.table[(i:nrow(test.table)),2])/sum(test.table[,2]),digits=4)
- non.parametric$Se.inf.Cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),sum(test.table[,2]),
- conf.level=CL)[4]),digits=4)
- non.parametric$Se.sup.Cl[i]<-round(as.numeric(binom.wilson(sum(test.table[i:nrow(test.table),2]),sum(test.table[,2]),
- conf.level=CL)[5]),digits=4)
- non.parametric$Specificity[i] <- round((sum(test.table[(1:i-1),1]))/(sum(test.table[,1])),digits=4)
- non.parametric$Sp.inf.Cl[i]<-round(as.numeric(binom.wilson(sum(test.table[(1:i-1),1]),sum(test.table[,1]),
- conf.level=CL)[4]),digits=4)
- non.parametric$Sp.sup.Cl[i]<-round(as.numeric(binom.wilson(sum(test.table[(1:i-1),1]),sum(test.table[,1]),
- conf.level=CL)[5]),digits=4)
- }
+ non.parametric$TP[i] <- sum(test.table[i:nrow(test.table),2])
+ non.parametric$FN[i] <- sum(test.table[1:i-1,2])
+ non.parametric$FP[i] <- sum(test.table[i:nrow(test.table),1])
+ non.parametric$TN[i] <- sum(test.table[1:i-1,1])
+ }
+
+ non.parametric$Sensitivity <- round(non.parametric$TP/D,digits=4)
+ non.parametric$Se.inf.cl <- round(binom.wilson(non.parametric$TP,D,conf.level=CL)[4]$lower,digits=4)
+ non.parametric$Se.sup.cl <- round(binom.wilson(non.parametric$TP,D,conf.level=CL)[5]$upper,digits=4)
+ non.parametric$Specificity <- round(non.parametric$TN/ND,digits=4)
+ non.parametric$Sp.inf.cl <- round(binom.wilson(non.parametric$TN,ND,conf.level=CL)[4]$lower,digits=4)
+ non.parametric$Sp.sup.cl <- round(binom.wilson(non.parametric$TN,ND,conf.level=CL)[5]$upper,digits=4)
+
non.parametric$PLR<-round(non.parametric$Sensitivity/(1-non.parametric$Specificity),digits=2)
- non.parametric$PLR.inf.cl<-round(exp(log(non.parametric$PLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-non.parametric$Sensitivity)/
- ((sum(non.parametric[,2]))*non.parametric$Specificity)+(non.parametric$Specificity)/((sum(non.parametric[,1]))*
- (1-non.parametric$Specificity)))),digits=2)
- non.parametric$PLR.sup.cl<-round(exp(log(non.parametric$PLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt(
- (1-non.parametric$Sensitivity)/((sum(non.parametric[,2]))*non.parametric$Specificity)+(
- non.parametric$Specificity)/((sum(non.parametric[,1]))*(1-non.parametric$Specificity)))),digits=2)
+ non.parametric$PLR.inf.cl<-round(exp(log(non.parametric$PLR)-(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-non.parametric$Sensitivity)/(
+ (D)*non.parametric$Specificity)+(non.parametric$Specificity)/((ND)*(1-non.parametric$Specificity)))),digits=2)
+ non.parametric$PLR.sup.cl<-round(exp(log(non.parametric$PLR)+(qnorm(1-((1-CL)/2),mean=0,sd=1))*sqrt((1-non.parametric$Sensitivity)/(
+ (D)*non.parametric$Specificity)+(non.parametric$Specificity)/((ND)*(1-non.parametric$Specificity)))),digits=2)
+
# Se=Sp cut-off
non.parametric$Se.equals.Sp<-abs(non.parametric$Specificity-non.parametric$Sensitivity)
# Efficiency= Se*prevalence+(1-prevalence)*Se
@@ -101,9 +101,9 @@
# )) Does Not work... somethig worng with the lines or columns selection
np.test.best.cutoff<-as.data.frame(rbind(
- non.parametric[which.min(non.parametric$Se.equals.Sp),1:10],
- non.parametric[which.max(non.parametric$Efficiency),1:10],
- non.parametric[which.min(non.parametric$MCT),1:10]))
+ non.parametric[which.min(non.parametric$Se.equals.Sp),c(1,6:14)],
+ non.parametric[which.max(non.parametric$Efficiency),c(1,6:14)],
+ non.parametric[which.min(non.parametric$MCT),c(1,6:14)]))
#best.cutoff,non.parametric[which.min(non.parametric$Se.equals.Sp),1:10])
#best.cutoff<-non.parametric$test.values[which.max(non.parametric$Efficiency)]
@@ -115,8 +115,8 @@
rownames(np.test.best.cutoff)<- c("Se=Sp","Max. Efficiency","Min. MCT")
non.parametric.inconclusive<-as.data.frame(rbind(
- non.parametric[which.min(abs(Inconclusive-non.parametric$Sensitivity)),1:10],
- non.parametric[which.min(abs(Inconclusive-non.parametric$Specificity)),1:10]))
+ non.parametric[which.min(abs(Inconclusive-non.parametric$Sensitivity)),c(1,6:14)],
+ non.parametric[which.min(abs(Inconclusive-non.parametric$Specificity)),c(1,6:14)]))
rownames(non.parametric.inconclusive)<-c("Lower inconclusive","Upper inconclusive")
if(is.null(t.max)){
@@ -135,27 +135,27 @@
test.values<-seq(t.min,t.max,precision)
parametric<-as.data.frame(test.values)
parametric$Sensitivity <- as.numeric(sim(net.Se$net, test.values))
- parametric$Se.inf.Cl<-parametric$Sensitivity - qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Sensitivity *
+ parametric$Se.inf.cl<-parametric$Sensitivity - qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Sensitivity *
(1 - parametric$Sensitivity))/(sample.size * sample.prevalence)))
- parametric$Se.sup.Cl<-parametric$Sensitivity + qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Sensitivity *
+ parametric$Se.sup.cl<-parametric$Sensitivity + qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Sensitivity *
(1 - parametric$Sensitivity))/(sample.size * sample.prevalence)))
parametric$Specificity <- as.numeric(sim(net.Sp$net, test.values))
- parametric$Sp.inf.Cl <- parametric$Specificity - qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Specificity *
+ parametric$Sp.inf.cl <- parametric$Specificity - qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Specificity *
(1 - parametric$Specificity))/(sample.size * (1-sample.prevalence))))
- parametric$Sp.sup.Cl <- parametric$Specificity + qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Specificity *
+ parametric$Sp.sup.cl <- parametric$Specificity + qnorm(1 - (1-conf.limit)/2) * sqrt(((parametric$Specificity *
(1 - parametric$Specificity))/(sample.size * (1-sample.prevalence))))
parametric$PLR <- parametric$Sensitivity/(1-parametric$Specificity)
- parametric$PLR.inf.Cl<-exp(log(parametric$PLR)-(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((1-parametric$Sensitivity)/
+ parametric$PLR.inf.cl<-exp(log(parametric$PLR)-(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((1-parametric$Sensitivity)/
((sample.size * sample.prevalence)*parametric$Specificity)+(parametric$Specificity)/((sample.size *
(1-sample.prevalence))*(1-parametric$Specificity))))
- parametric$PLR.sup.Cl<-exp(log(parametric$PLR)+(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((1-parametric$Sensitivity)/
+ parametric$PLR.sup.cl<-exp(log(parametric$PLR)+(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((1-parametric$Sensitivity)/
((sample.size * sample.prevalence)*parametric$Specificity)+(parametric$Specificity)/((sample.size *
(1-sample.prevalence))*(1-parametric$Specificity))))
#parametric$NLR <- (1-parametric$Specificity)/parametric$Sensitivity
- #parametric$NLR.inf.Cl <- exp(log(parametric$NLR)-(qnorm(1-((1-(1-conf.limit))/2),mean=0,sd=1))*
+ #parametric$NLR.inf.cl <- exp(log(parametric$NLR)-(qnorm(1-((1-(1-conf.limit))/2),mean=0,sd=1))*
# sqrt((parametric$Sensitivity)/((sample.size * sample.prevalence)*(1-parametric$Sensitivity))+(1-parametric$Specificity)/
# ((sample.size * (1-sample.prevalence))*(parametric$Specificity))))
- #parametric$NLR.sup.Cl <- exp(log(parametric$NLR)+(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((parametric$Sensitivity)/
+ #parametric$NLR.sup.cl <- exp(log(parametric$NLR)+(qnorm(1-((1-conf.limit)/2),mean=0,sd=1))*sqrt((parametric$Sensitivity)/
# ((sample.size * sample.prevalence)*(1-parametric$Sensitivity))+(1-parametric$Specificity)/((sample.size *
# (1-sample.prevalence))*(parametric$Specificity))))
parametric$Se.equals.Sp<-abs(parametric$Sensitivity-parametric$Specificity)
@@ -173,125 +173,9 @@
parametric[which.max(parametric$Efficiency),1:10],
parametric[which.min(parametric$MCT),1:10]
))
- rownames(par.test.best.cutoff)<- c("Se=Sp","Max. Efficiency","Min. MCT")
-
-# Plot Commands
- if(Plot!="None"){
- if(Plot=="Parametric"|Plot=="Both"){
- plot(parametric$test.values,parametric$Sensitivity,ylim=c(0,1),type="l",col=2, xlab="Test scale",
- ylab="Sensitivity & Specificity",cex=cex,lty=1)
- lines(parametric$test.values,parametric$Specificity,col=4,type="l",lty=2,cex=cex)
- if(Plot=="Both"){
- lines(non.parametric$test.values,non.parametric$Sensitivity,col=2,type="o",lty=1,cex=cex)
- lines(non.parametric$test.values,non.parametric$Specificity,col=4,type="o",lty=2,cex=cex)
- }
- leg.txt<-c("Se", "Sp")
- fill.col<-c(2,4)
- line.type<-c(1,2)
- subtitle<-""
- }
- if(Plot=="Non-parametric"){
- plot(non.parametric$test.values,non.parametric$Sensitivity,
- ylim=c(0,1),type="o",col=2, xlab="test scale",ylab="Sensitivity & Specificity",
- lty=1,cex=cex)
- lines(non.parametric$test.values,non.parametric$Specificity,col=4,type="o",lty=2,cex=cex)
- leg.txt<-c("Se", "Sp")
- fill.col<-c(2,4)
- line.type<-c(1,2)
- subtitle<-""
- }
- if(Plot.inc.range==TRUE){
- abline(h=Inconclusive,col="lightgray",lty=4)
- if(Plot=="Parametric"|Plot=="Both"){
- abline(v=(parametric.inconclusive[1,1]),col="lightgray",lty=4)
- abline(v=(parametric.inconclusive[2,1]),col="lightgray",lty=4)
- subtitle<-paste("Parametric inconclusive limits at",formatC(inc),"level:",formatC(parametric.inconclusive[1,1]),
- "-",formatC(parametric.inconclusive[2,1]),".")
- }
- if(Plot=="Non-parametric"){
- abline(v=(non.parametric.inconclusive[1,1]),col="lightgray",lty=4)
- abline(v=(non.parametric.inconclusive[2,1]),col="lightgray",lty=4)
- subtitle<-paste("Non-parametric inconclusive limits at",formatC(inc),"level:",formatC(non.parametric.inconclusive[1,1]),
- "-",formatC(non.parametric.inconclusive[2,1]),".")
- }
- leg.txt<-c(leg.txt,c("Inc limits"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,4)
- }
- if(Plot.Cl==TRUE){
- if(Plot=="Both"|Plot=="Parametric"){
- lines(parametric$test.values,parametric$Se.inf.Cl,lty=5,col=2)
- lines(parametric$test.values,parametric$Se.sup.Cl,lty=5,col=2)
- lines(parametric$test.values,parametric$Sp.inf.Cl,lty=3,col=4)
- lines(parametric$test.values,parametric$Sp.sup.Cl,lty=3,col=4)
- }
- if(Plot=="Non-parametric"){
- lines(non.parametric$test.values,non.parametric$Se.inf.Cl,lty=5,col=2)
- lines(non.parametric$test.values,non.parametric$Se.sup.Cl,lty=5,col=2)
- lines(non.parametric$test.values,non.parametric$Sp.inf.Cl,lty=3,col=4)
- lines(non.parametric$test.values,non.parametric$Sp.sup.Cl,lty=3,col=4)
- }
- leg.txt<-c(leg.txt,c("Se conf. band","Sp conf. band"))
- fill.col<-c(fill.col,c(2,4))
- line.type<-c(line.type,5,3)
- }
- if(Plot.cutoff=="Se=Sp"){
- if(Plot=="Both"|Plot=="Parametric"){
- abline(v=(par.test.best.cutoff[1,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by parametric Se=Sp:",formatC(par.test.best.cutoff[1,1])))
- }
- if(Plot=="Non-parametric"){
- abline(v=(np.test.best.cutoff[1,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by Se=Sp:",formatC(np.test.best.cutoff[1,1])))
- }
+ rownames(par.test.best.cutoff)<- c("Se=Sp","Max. Efficiency","Min. MCT")
- }
- if(Plot.cutoff=="Max.Efficiency"){
- if(Plot=="Both"|Plot=="Parametric"){
- abline(v=(par.test.best.cutoff[2,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by parametric Max. Efficiency:",formatC(par.test.best.cutoff[2,1]),"."))
- #"Pop. prevalence:",formatC(pop.prevalence)))) Does not fit in the graph
- }
- if(Plot=="Non-parametric"){
- abline(v=(np.test.best.cutoff[2,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by Max. Efficiency:",formatC(np.test.best.cutoff[2,1]),"."))
- #"Pop. prevalence:",formatC(pop.prevalence)))) Does not fit in the graph
- }
- }
- if(Plot.cutoff=="Min.MCT"){
- if(Plot=="Both"|Plot=="Parametric"){
- abline(v=(par.test.best.cutoff[3,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by minimizing paramaetric MCT:",formatC(np.test.best.cutoff[3,1],".")))
- #,"Pop. prevalence:",formatC(pop.prevalence),"Cost FN/FP:",formatC(cost))) Does not fit in the graph
- }
- if(Plot=="Non-parametric"){
- abline(v=(np.test.best.cutoff[3,1]),col="lightgray",lty=6)
- leg.txt<-c(leg.txt,c("Best cut-off"))
- fill.col<-c(fill.col,c("lightgray"))
- line.type<-c(line.type,6)
- subtitle<-paste(subtitle,paste("Cut-off estimated by Minimizing MCT:",formatC(np.test.best.cutoff[3,1],".")))
- #,"Pop. prevalence:",formatC(pop.prevalence),"Cost FN/FP:",formatC(cost))) Does not fit in the graph
- }
- }
- legend("right",legend=leg.txt,col=fill.col,lty=line.type, bty="n")
- title(sub=subtitle,cex.sub=cex.sub)
- }
- rm(test.values)
+# rm(test.values)
if(non.parametric.inconclusive[1,1]>non.parametric.inconclusive[2,1]){
warning("Non-parametric lower inconclusive limit is higher than upper inconclusive limit.")
}
@@ -301,22 +185,22 @@
if(np.test.best.cutoff[1,1]>non.parametric.inconclusive[2,1]|
np.test.best.cutoff[2,1]>non.parametric.inconclusive[2,1]|
np.test.best.cutoff[3,1]>non.parametric.inconclusive[2,1]){
- warning("At least one fo the non-parametric best cut-off is higher then upper inconclusive limit.")
+ warning("At least one of the non-parametric best cut-off is higher then upper inconclusive limit.")
}
if(np.test.best.cutoff[1,1]<non.parametric.inconclusive[1,1]|
np.test.best.cutoff[2,1]<non.parametric.inconclusive[1,1]|
np.test.best.cutoff[3,1]<non.parametric.inconclusive[1,1]){
- warning("At least one fo the non-parametric best cut-off is lower then lower inconclusive limit.")
+ warning("At least one of the non-parametric best cut-off is lower then lower inconclusive limit.")
}
if(par.test.best.cutoff[1,1]>parametric.inconclusive[2,1]|
par.test.best.cutoff[2,1]>parametric.inconclusive[2,1]|
par.test.best.cutoff[3,1]>parametric.inconclusive[2,1]){
- warning("At least one fo the parametric best cut-off is higher then upper inconclusive limit.")
+ warning("At least one of the parametric best cut-off is higher then upper inconclusive limit.")
}
if(par.test.best.cutoff[1,1]<parametric.inconclusive[1,1]|
par.test.best.cutoff[2,1]<parametric.inconclusive[1,1]|
par.test.best.cutoff[3,1]<parametric.inconclusive[1,1]){
- warning("At least one fo the parametric best cut-off is lower then lower inconclusive limit.")
+ warning("At least one of the parametric best cut-off is lower then lower inconclusive limit.")
}
reteval<-list(sample.size = sample.size,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/diagnosismed -r 16
More information about the Diagnosismed-commits
mailing list