1
/* pROC: Tools Receiver operating characteristic (ROC curves) with
2
   (partial) area under the curve, confidence intervals and comparison. 
3
   Copyright (C) 2016 Xavier Robin, Stefan Siegert <stefan_siegert@gmx.de>
4

5
   This program is free software: you can redistribute it and/or modify
6
   it under the terms of the GNU General Public License as published by
7
   the Free Software Foundation, either version 3 of the License, or
8
   (at your option) any later version.
9
  
10
   This program is distributed in the hope that it will be useful,
11
   but WITHOUT ANY WARRANTY; without even the implied warranty of
12
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
   GNU General Public License for more details.
14
  
15
   You should have received a copy of the GNU General Public License
16
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
17
*/
18

19
#include <Rcpp.h> 
20
using namespace Rcpp;
21

22 29
bool _cmp(std::pair<int, double> l, std::pair<int, double> r) {
23 29
  return l.second < r.second;
24
}
25

26
// [[Rcpp::export]]
27 29
List delongPlacementsCpp(List roc) {
28

29
  int i, j, k, m, n, mdupl, ndupl, L;
30

31 29
  std::vector<double> cases = roc["cases"];
32 29
  std::vector<double> controls = roc["controls"];
33 29
  std::string direction = roc["direction"];
34 29
  m = cases.size();
35 29
  n = controls.size();
36 29
  L = m + n;
37
  
38
  // For direction ">" we must reverse the data
39 29
  if (direction == ">") {
40 29
    for (i = 0; i < m; i++) {
41 29
      cases[i] = -cases[i];
42
    }
43 29
    for (i = 0; i < n; i++) {
44 29
      controls[i] = -controls[i];
45
    }
46
  }
47

48
  // concatenate cases and controls into a vector of L pairs of the form
49
  // (index, value), also save class labels (1 for cases, 0 for controls)
50 29
  std::vector< std::pair<int, double> > Z;
51 29
  std::vector< bool > labels;
52 29
  for (i = 0; i < m; i++) {
53 29
    Z.push_back(std::pair<int, double>(i, cases.at(i)));
54 29
    labels.push_back(true);
55
  }
56 29
  Rcpp::checkUserInterrupt();
57 29
  for (j = 0; j < n; j++) {
58 29
    Z.push_back(std::pair<int, double>(m+j, controls.at(j)));
59 29
    labels.push_back(false);
60
  }
61 29
  Rcpp::checkUserInterrupt();
62

63
  // sort Z from smallest to largest value, so Z holds the order indices and
64
  // order statistics of all classifiers
65 29
  std::sort(Z.begin(), Z.end(), _cmp);
66 29
  Rcpp::checkUserInterrupt();
67

68
  // the following calculates the "Delong-placements" X and Y in a single pass
69
  // over the vector Z, instead of having to double loop over all pairs of
70
  // (X_i, Y_j)
71 29
  std::vector< double > XY(L, 0.0);  // vector to hold the unnormalised X and Y values
72 29
  std::vector< int > X_inds, Y_inds; // temporary vectors to save indices of duplicates
73 29
  m = n = i = 0;                     // initialisation
74 29
  while (i < L) {
75 29
    X_inds.clear();
76 29
    Y_inds.clear();
77 29
    mdupl = ndupl = 0;
78 29
    if (i % 10000 == 0) Rcpp::checkUserInterrupt();
79 29
    while(1) {
80 29
      j = Z.at(i).first;
81 29
      if (labels.at(j)) {
82 29
        mdupl++;
83 29
        X_inds.push_back(j);
84
      } else {
85 29
        ndupl++;
86 29
        Y_inds.push_back(j);
87
      }
88 29
      if (i == L-1) {
89 29
        break;
90
      }
91 29
      if (Z.at(i).second != Z.at(i+1).second) {
92 29
        break;
93
      }
94 29
      i++;
95
    }
96 29
    for (k = 0; k < mdupl; k++) {
97 29
      XY.at(X_inds.at(k)) = n + ndupl/2.0;
98
    }
99 29
    for (k = 0; k < ndupl; k++) {
100 29
      XY.at(Y_inds.at(k)) = m + mdupl/2.0;
101
    }
102 29
    n += ndupl;
103 29
    m += mdupl;
104 29
    i++;
105
  }
106

107 29
  double sum = 0.0;
108 29
  std::vector<double> X, Y;
109 29
  Rcpp::checkUserInterrupt();
110

111 29
  for (i = 0; i < L; i++) {
112 29
    if (labels.at(i)) {
113 29
      sum += XY.at(i);
114 29
      X.push_back(XY.at(i) / n);
115
    } else {
116 29
      Y.push_back(1.0 - XY.at(i) / m);
117
    }
118
  }
119

120 29
  List ret;
121 29
  ret["theta"] = sum / m / n;
122 29
  ret["X"] = X;
123 29
  ret["Y"] = Y;
124 29
  return(ret);
125
}

Read our documentation on viewing source code .

Loading