src/language/stats/roc.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include "roc.h"
  20 #include <data/procedure.h>
  21 #include <language/lexer/variable-parser.h>
  22 #include <language/lexer/value-parser.h>
  23 #include <language/lexer/lexer.h>
  24
  25 #include <data/casegrouper.h>
  26 #include <data/casereader.h>
  27 #include <data/casewriter.h>
  28 #include <data/dictionary.h>
  29 #include <math/sort.h>
  30
  31 #include <libpspp/misc.h>
  32
  33 #include <gsl/gsl_cdf.h>
  34 #include <output/table.h>
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38 #define N_(msgid) msgid
  39
  40 struct cmd_roc
  41 {
  42   size_t n_vars;
  43   const struct variable **vars;
  44
  45   struct variable *state_var ;
  46   union value state_value;
  47
  48   /* Plot the roc curve */
  49   bool curve;
  50   /* Plot the reference line */
  51   bool reference;
  52
  53   double ci;
  54
  55   bool print_coords;
  56   bool print_se;
  57   bool bi_neg_exp; /* True iff the bi-negative exponential critieria
  58                       should be used */
  59   enum mv_class exclude;
  60
  61   bool invert ; /* True iff a smaller test result variable indicates
  62                    a positive result */
  63 };
  64
  65 static int run_roc (struct dataset *ds, struct cmd_roc *roc);
  66
  67 int
  68 cmd_roc (struct lexer *lexer, struct dataset *ds)
  69 {
  70   struct cmd_roc roc ;
  71   const struct dictionary *dict = dataset_dict (ds);
  72
  73   roc.vars = NULL;
  74   roc.n_vars = 0;
  75   roc.print_se = false;
  76   roc.print_coords = false;
  77   roc.exclude = MV_ANY;
  78   roc.curve = true;
  79   roc.reference = false;
  80   roc.ci = 95;
  81   roc.bi_neg_exp = false;
  82   roc.invert = false;
  83
  84   if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars,
  85                               PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
  86     return 2;
  87
  88   if ( ! lex_force_match (lexer, T_BY))
  89     {
  90       return 2;
  91     }
  92
  93   roc.state_var = parse_variable (lexer, dict);
  94
  95   if ( !lex_force_match (lexer, '('))
  96     {
  97       return 2;
  98     }
  99
 100   parse_value (lexer, &roc.state_value, var_get_width (roc.state_var));
 101
 102
 103   if ( !lex_force_match (lexer, ')'))
 104     {
 105       return 2;
 106     }
 107
 108
 109   while (lex_token (lexer) != '.')
 110     {
 111       lex_match (lexer, '/');
 112       if (lex_match_id (lexer, "MISSING"))
 113         {
 114           lex_match (lexer, '=');
 115           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 116             {
 117               if (lex_match_id (lexer, "INCLUDE"))
 118                 {
 119                   roc.exclude = MV_SYSTEM;
 120                 }
 121               else if (lex_match_id (lexer, "EXCLUDE"))
 122                 {
 123                   roc.exclude = MV_ANY;
 124                 }
 125               else
 126                 {
 127                   lex_error (lexer, NULL);
 128                   return 2;
 129                 }
 130             }
 131         }
 132       else if (lex_match_id (lexer, "PLOT"))
 133         {
 134           lex_match (lexer, '=');
 135           if (lex_match_id (lexer, "CURVE"))
 136             {
 137               roc.curve = true;
 138               if (lex_match (lexer, '('))
 139                 {
 140                   roc.reference = true;
 141                   lex_force_match_id (lexer, "REFERENCE");
 142                   lex_force_match (lexer, ')');
 143                 }
 144             }
 145           else if (lex_match_id (lexer, "NONE"))
 146             {
 147               roc.curve = false;
 148             }
 149           else
 150             {
 151               lex_error (lexer, NULL);
 152               return 2;
 153             }
 154         }
 155       else if (lex_match_id (lexer, "PRINT"))
 156         {
 157           lex_match (lexer, '=');
 158           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 159             {
 160               if (lex_match_id (lexer, "SE"))
 161                 {
 162                   roc.print_se = true;
 163                 }
 164               else if (lex_match_id (lexer, "COORDINATES"))
 165                 {
 166                   roc.print_coords = true;
 167                 }
 168               else
 169                 {
 170                   lex_error (lexer, NULL);
 171                   return 2;
 172                 }
 173             }
 174         }
 175       else if (lex_match_id (lexer, "CRITERIA"))
 176         {
 177           lex_match (lexer, '=');
 178           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 179             {
 180               if (lex_match_id (lexer, "CUTOFF"))
 181                 {
 182                   lex_force_match (lexer, '(');
 183                   if (lex_match_id (lexer, "INCLUDE"))
 184                     {
 185                       roc.exclude = MV_SYSTEM;
 186                     }
 187                   else if (lex_match_id (lexer, "EXCLUDE"))
 188                     {
 189                       roc.exclude = MV_USER | MV_SYSTEM;
 190                     }
 191                   else
 192                     {
 193                       lex_error (lexer, NULL);
 194                       return 2;
 195                     }
 196                   lex_force_match (lexer, ')');
 197                 }
 198               else if (lex_match_id (lexer, "TESTPOS"))
 199                 {
 200                   lex_force_match (lexer, '(');
 201                   if (lex_match_id (lexer, "LARGE"))
 202                     {
 203                       roc.invert = false;
 204                     }
 205                   else if (lex_match_id (lexer, "SMALL"))
 206                     {
 207                       roc.invert = true;
 208                     }
 209                   else
 210                     {
 211                       lex_error (lexer, NULL);
 212                       return 2;
 213                     }
 214                   lex_force_match (lexer, ')');
 215                 }
 216               else if (lex_match_id (lexer, "CI"))
 217                 {
 218                   lex_force_match (lexer, '(');
 219                   lex_force_num (lexer);
 220                   roc.ci = lex_number (lexer);
 221                   lex_get (lexer);
 222                   lex_force_match (lexer, ')');
 223                 }
 224               else if (lex_match_id (lexer, "DISTRIBUTION"))
 225                 {
 226                   lex_force_match (lexer, '(');
 227                   if (lex_match_id (lexer, "FREE"))
 228                     {
 229                       roc.bi_neg_exp = false;
 230                     }
 231                   else if (lex_match_id (lexer, "NEGEXPO"))
 232                     {
 233                       roc.bi_neg_exp = true;
 234                     }
 235                   else
 236                     {
 237                       lex_error (lexer, NULL);
 238                       return 2;
 239                     }
 240                   lex_force_match (lexer, ')');
 241                 }
 242               else
 243                 {
 244                   lex_error (lexer, NULL);
 245                   return 2;
 246                 }
 247             }
 248         }
 249       else
 250         {
 251           lex_error (lexer, NULL);
 252           break;
 253         }
 254     }
 255
 256   run_roc (ds, &roc);
 257
 258   return 1;
 259 }
 260
 261
 262
 263
 264 static void
 265 do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict);
 266
 267
 268 static int
 269 run_roc (struct dataset *ds, struct cmd_roc *roc)
 270 {
 271   struct dictionary *dict = dataset_dict (ds);
 272   bool ok;
 273   struct casereader *group;
 274
 275   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
 276   while (casegrouper_get_next_group (grouper, &group))
 277     {
 278       do_roc (roc, group, dataset_dict (ds));
 279     }
 280   ok = casegrouper_destroy (grouper);
 281   ok = proc_commit (ds) && ok;
 282
 283   return ok;
 284 }
 285
 286
 287 static void
 288 dump_casereader (struct casereader *reader)
 289 {
 290   struct ccase *c;
 291   struct casereader *r = casereader_clone (reader);
 292
 293   for ( ; (c = casereader_read (r) ); case_unref (c))
 294     {
 295       int i;
 296       for (i = 0 ; i < case_get_value_cnt (c); ++i)
 297         {
 298           printf ("%g ", case_data_idx (c, i)->f);
 299         }
 300       printf ("\n");
 301     }
 302
 303   casereader_destroy (r);
 304 }
 305
 306 static bool
 307 match_positives (const struct ccase *c, void *aux)
 308 {
 309   struct cmd_roc *roc = aux;
 310
 311   return 0 == value_compare_3way (case_data (c, roc->state_var),
 312                                  &roc->state_value,
 313                                  var_get_width (roc->state_var));
 314 }
 315
 316
 317 #define VALUE  0
 318 #define N_EQ   1
 319 #define N_PRED 2
 320
 321 struct roc_state
 322 {
 323   double auc;
 324
 325   double n1;
 326   double n2;
 327
 328   double q1hat;
 329   double q2hat;
 330 };
 331
 332
 333 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
 334
 335
 336
 337 static struct casereader *
 338 process_group (const struct variable *var, struct casereader *reader,
 339                bool (*pred) (double, double),
 340                const struct dictionary *dict,
 341                double *cc)
 342 {
 343   const struct variable *w = dict_get_weight (dict);
 344   const int weight_idx  = w ? var_get_case_index (w) :
 345     caseproto_get_n_widths (casereader_get_proto (reader)) - 1;
 346
 347   struct casereader *r1 =
 348     casereader_create_distinct (sort_execute_1var (reader, var), var, w);
 349
 350   struct ccase *c1;
 351
 352   struct casereader *rclone = casereader_clone (r1);
 353   struct casewriter *wtr;
 354   struct caseproto *proto = caseproto_create ();
 355
 356   proto = caseproto_add_width (proto, 0);
 357   proto = caseproto_add_width (proto, 0);
 358   proto = caseproto_add_width (proto, 0);
 359
 360   wtr = autopaging_writer_create (proto);
 361
 362   *cc = 0;
 363
 364   for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
 365     {
 366       struct ccase *c2;
 367       struct casereader *r2 = casereader_clone (rclone);
 368
 369       const double weight1 = case_data_idx (c1, weight_idx)->f;
 370       const double d1 = case_data (c1, var)->f;
 371       double n_eq = 0.0;
 372       double n_pred = 0.0;
 373
 374
 375       struct ccase *new_case = case_create (proto);
 376
 377       *cc += weight1;
 378
 379       for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
 380         {
 381           const double d2 = case_data (c2, var)->f;
 382           const double weight2 = case_data_idx (c2, weight_idx)->f;
 383
 384           if ( d1 == d2 )
 385             {
 386               n_eq += weight2;
 387               continue;
 388             }
 389           else  if ( pred (d2, d1))
 390             {
 391               n_pred += weight2;
 392             }
 393         }
 394
 395       case_data_rw_idx (new_case, VALUE)->f = d1;
 396       case_data_rw_idx (new_case, N_EQ)->f = n_eq;
 397       case_data_rw_idx (new_case, N_PRED)->f = n_pred;
 398
 399       casewriter_write (wtr, new_case);
 400
 401       casereader_destroy (r2);
 402     }
 403
 404   casereader_destroy (r1);
 405   casereader_destroy (rclone);
 406
 407   return casewriter_make_reader (wtr);
 408 }
 409
 410 static bool
 411 gt (double d1, double d2)
 412 {
 413   return d1 > d2;
 414 }
 415
 416 static bool
 417 lt (double d1, double d2)
 418 {
 419   return d1 < d2;
 420 }
 421
 422
 423 static void
 424 do_roc (struct cmd_roc *roc, struct casereader *input, struct dictionary *dict)
 425 {
 426   int i;
 427
 428   struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
 429
 430   const struct caseproto *proto = casereader_get_proto (input);
 431
 432   struct casewriter *neg_wtr = autopaging_writer_create (proto);
 433
 434   struct casereader *negatives = NULL;
 435
 436   struct casereader *positives =
 437     casereader_create_filter_func (input,
 438                                    match_positives,
 439                                    NULL,
 440                                    roc,
 441                                    neg_wtr);
 442
 443
 444   for (i = 0 ; i < roc->n_vars; ++i)
 445     {
 446       struct ccase *cpos;
 447       struct casereader *n_neg ;
 448       const struct variable *var = roc->vars[i];
 449
 450       struct casereader *neg ;
 451       struct casereader *pos = casereader_clone (positives);
 452
 453       struct casereader *n_pos = process_group (var, pos, gt, dict, &rs[i].n1);
 454
 455       if ( negatives == NULL)
 456         {
 457           negatives = casewriter_make_reader (neg_wtr);
 458         }
 459
 460       neg = casereader_clone (negatives);
 461
 462       n_neg = process_group (var, neg, lt, dict, &rs[i].n2);
 463
 464       /* Simple join on VALUE */
 465       for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos))
 466         {
 467           struct ccase *cneg = NULL;
 468           double dneg = -DBL_MAX;
 469           const double dpos = case_data_idx (cpos, VALUE)->f;
 470           while (dneg < dpos)
 471             {
 472               if ( cneg )
 473                 case_unref (cneg);
 474
 475               cneg = casereader_read (n_neg);
 476               dneg = case_data_idx (cneg, VALUE)->f;
 477             }
 478
 479           if ( dpos == dneg )
 480             {
 481               double n_pos_eq = case_data_idx (cpos, N_EQ)->f;
 482               double n_neg_eq = case_data_idx (cneg, N_EQ)->f;
 483               double n_pos_gt = case_data_idx (cpos, N_PRED)->f;
 484               double n_neg_lt = case_data_idx (cneg, N_PRED)->f;
 485
 486               rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
 487               rs[i].q1hat +=
 488                 n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
 489               rs[i].q2hat +=
 490                 n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
 491             }
 492
 493           if ( cneg )
 494             case_unref (cneg);
 495         }
 496
 497       rs[i].auc /=  rs[i].n1 * rs[i].n2;
 498       if ( roc->invert )
 499         rs[i].auc = 1 - rs[i].auc;
 500
 501       if ( roc->bi_neg_exp )
 502         {
 503           rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
 504           rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
 505         }
 506       else
 507         {
 508           rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
 509           rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
 510         }
 511     }
 512
 513   casereader_destroy (positives);
 514   casereader_destroy (negatives);
 515
 516   output_roc (rs, roc);
 517
 518   free (rs);
 519 }
 520
 521
 522
 523
 524 static void
 525 show_auc  (struct roc_state *rs, const struct cmd_roc *roc)
 526 {
 527   int i;
 528   const int n_fields = roc->print_se ? 5 : 1;
 529   const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields;
 530   const int n_rows = 2 + roc->n_vars;
 531   struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
 532
 533   if ( roc->n_vars > 1)
 534     tab_title (tbl, _("Area Under the Curve"));
 535   else
 536     tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
 537
 538   tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
 539
 540   tab_dim (tbl, tab_natural_dimensions, NULL);
 541
 542   tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
 543
 544   tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
 545
 546   tab_box (tbl,
 547            TAL_2, TAL_2,
 548            -1, TAL_1,
 549            0, 0,
 550            n_cols - 1,
 551            n_rows - 1);
 552
 553   if ( roc->print_se )
 554     {
 555       tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
 556       tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
 557
 558       tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
 559       tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
 560
 561       tab_joint_text (tbl, n_cols - 2, 0, 4, 0,
 562                       TAT_TITLE | TAB_CENTER | TAT_PRINTF,
 563                       _("Asymp. %g%% Confidence Interval"), roc->ci);
 564       tab_vline (tbl, 0, n_cols - 1, 0, 0);
 565       tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1);
 566     }
 567
 568   if ( roc->n_vars > 1)
 569     tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
 570
 571   if ( roc->n_vars > 1)
 572     tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
 573
 574
 575   for ( i = 0 ; i < roc->n_vars ; ++i )
 576     {
 577       tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
 578
 579       tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
 580
 581       if ( roc->print_se )
 582         {
 583
 584           double se ;
 585           const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
 586                                       (12 * rs[i].n1 * rs[i].n2));
 587           double ci ;
 588           double yy ;
 589
 590           se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (rs[i].q1hat - pow2 (rs[i].auc)) +
 591             (rs[i].n2 - 1) * (rs[i].q2hat - pow2 (rs[i].auc));
 592
 593           se /= rs[i].n1 * rs[i].n2;
 594
 595           se = sqrt (se);
 596
 597           tab_double (tbl, n_cols - 4, 2 + i, 0,
 598                       se,
 599                       NULL);
 600
 601           ci = 1 - roc->ci / 100.0;
 602           yy = gsl_cdf_gaussian_Qinv (ci, se) ;
 603
 604           tab_double (tbl, n_cols - 2, 2 + i, 0,
 605                       rs[i].auc - yy,
 606                       NULL);
 607
 608           tab_double (tbl, n_cols - 1, 2 + i, 0,
 609                       rs[i].auc + yy,
 610                       NULL);
 611
 612           tab_double (tbl, n_cols - 3, 2 + i, 0,
 613                       2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)),
 614                       NULL);
 615         }
 616     }
 617
 618   tab_submit (tbl);
 619 }
 620
 621
 622
 623
 624 static void
 625 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
 626 {
 627 #if 0
 628   show_summary (roc);
 629
 630   if ( roc->curve )
 631     draw_roc (rs, roc);
 632 #endif
 633
 634   show_auc (rs, roc);
 635
 636 #if 0
 637   if ( roc->print_coords )
 638     show_coords (rs, roc);
 639 #endif
 640 }