1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
# Accuracy measures, modified from the WGCNA version.
# Helper function: contingency table of 2 variables that will also include rows/columns for levels that do
# not appear in x or y.
.table2.allLevels = function(x, y, levels.x = sort(unique(x)), levels.y = sort(unique(y)), setNames = FALSE)
{
nx = length(levels.x);
ny = length(levels.y);
t = table(x, y);
out = matrix(0, nx, ny);
if (setNames)
{
rownames(out) = levels.x;
colnames(out) = levels.y;
}
out[ match(rownames(t), levels.x), match(colnames(t), levels.y) ] = t;
out;
}
# accuracy measures
accuracyMeasures = function(predicted, observed = NULL, type = c("auto", "binary", "quantitative"),
levels = if (isTRUE(all.equal(dim(predicted), c(2,2)))) colnames(predicted)
else if (is.factor(predicted))
sort(unique(c(as.character(predicted), as.character(observed))))
else sort(unique(c(observed, predicted))),
negativeLevel = levels[2], positiveLevel = levels[1] )
{
type = match.arg(type);
if (type=="auto")
{
if (!is.null(dim(predicted)))
{
if (isTRUE(all.equal(dim(predicted), c(2,2))))
{
type = "binary"
} else
stop("If supplying a matrix in 'predicted', it must be a 2x2 contingency table.");
} else {
if (is.null(observed))
stop("When 'predicted' is a vector, 'observed' must be given and have the same length as 'predicted'.");
if (length(levels)==2)
{
type = "binary"
} else
type = "quantitative"
}
}
if (type=="binary")
{
if (is.null(dim(predicted)))
{
if (is.null(observed))
stop("When 'predicted' is a vector, 'observed' must be given and have the same length as 'predicted'.");
if ( length(predicted)!=length(observed) )
stop("When both 'predicted' and 'observed' are given, they must be vectors of the same length.");
if (length(levels)!=2)
stop("'levels' must contain 2 entries (the possible values of the binary variables\n",
" 'predicted' and 'observed').");
tab = .table2.allLevels(predicted, observed, levels.x = levels, levels.y = levels, setNames = TRUE);
} else {
tab = predicted;
if (is.null(colnames(tab)) | is.null(rownames(tab)))
stop("When 'predicted' is a contingency table, it must have valid colnames and rownames.");
}
if ( ncol(tab) !=2 | nrow(tab) !=2 )
stop("The input table must be a 2x2 table. ")
if (negativeLevel==positiveLevel)
stop("'negativeLevel' and 'positiveLevel' cannot be the same.");
neg = match(negativeLevel, colnames(tab));
if (is.na(neg))
stop(.spaste("Cannot find the negative level ", negativeLevel,
" among the colnames of the contingency table.\n Please check the input and try again."))
pos = match(positiveLevel, colnames(tab));
if (is.na(pos))
stop(.spaste("Cannot find the positive level ", positiveLevel,
" among the colnames of the contingency table.\n Please check the input and try again."))
if ( sum(is.na(tab) ) )
warning("Missing data should not be present in input.\n",
" Suggestion: check whether NA should be coded as 0.")
is.wholenumber =function(x, tol = .Machine$double.eps^0.5) { abs(x - round(x)) < tol }
if ( sum( !is.wholenumber(tab), na.rm=T ) >0)
warning("STRONG WARNING: The input table contains non-integers, which does not make sense.")
if ( sum( tab<0, na.rm=T ) >0)
stop("The input table cannot contain negative numbers.");
num1=sum(diag(tab),na.rm=T)
denom1=sum(tab,na.rm=T)
if (denom1==0)
warning("The input table has zero observations (sum of all cells is zero).")
TP=tab[pos, pos]
FP=tab[pos, neg]
FN=tab[neg, pos]
TN=tab[neg, neg]
error.rate= ifelse(denom1==0,NA, 1-num1/denom1)
Accuracy= ifelse(denom1==0,NA, num1/denom1 )
Specificity= ifelse(FP + TN==0, NA, TN / (FP + TN) )
Sensitivity= ifelse(TP + FN==0, NA, TP / (TP + FN) )
NegativePredictiveValue= ifelse(FN + TN==0,NA, TN / (FN + TN) )
PositivePredictiveValue=ifelse(TP + FP==0,NA, TP / (TP + FP) )
FalsePositiveRate = 1 - Specificity
FalseNegativeRate = 1 - Sensitivity
Power = Sensitivity
LikelihoodRatioPositive = ifelse(1 - Specificity==0,NA, Sensitivity / (1 - Specificity) )
LikelihoodRatioNegative = ifelse(Specificity==0, NA, (1 - Sensitivity) / Specificity )
NaiveErrorRate = ifelse(denom1==0,NA,
min(c(tab[pos, pos]+ tab[neg, pos] , tab[pos, neg]+ tab[neg, neg] ))/denom1 )
out=data.frame(
Measure= c("Error.Rate","Accuracy", "Specificity","Sensitivity","NegativePredictiveValue",
"PositivePredictiveValue","FalsePositiveRate","FalseNegativeRate","Power",
"LikelihoodRatioPositive","LikelihoodRatioNegative", "NaiveErrorRate", "NegativeLevel",
"PositiveLevel"),
Value=c(error.rate,Accuracy, Specificity,Sensitivity,NegativePredictiveValue,
PositivePredictiveValue,FalsePositiveRate,FalseNegativeRate,Power,
LikelihoodRatioPositive,LikelihoodRatioNegative,NaiveErrorRate, negativeLevel,
positiveLevel));
} else if (type=="quantitative")
{
if (!is.null(dim(predicted)))
stop("When 'type' is \"quantitative\", 'predicted' cannot be a 2-dimensional matrix.");
if (length(predicted)!=length(observed))
stop("'predicted' and 'observed' must be vectors of the same length.");
cr = cor(predicted, observed, use = 'p');
out = data.frame(
Measure = c("Cor", "R.squared", "MeanSquareError", "MedianAbsoluteError", "Cindex"),
Value = c(cr,
cr^2,
mean( (predicted-observed)^2,na.rm=TRUE),
median((predicted-observed)^2,na.rm=TRUE),
rcorr.cens(predicted,observed,outx=TRUE)[[1]]));
}
out;
}
|
