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
 329
 330 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
 331
 332
 333
 334 static struct casereader *
 335 process_group (const struct variable *var, struct casereader *reader,
 336                bool (*pred) (double, double),
 337                const struct dictionary *dict,
 338                double *cc)
 339 {
 340   const struct variable *w = dict_get_weight (dict);
 341   const int weight_idx  = w ? var_get_case_index (w) :
 342     caseproto_get_n_widths (casereader_get_proto (reader)) - 1;
 343
 344   struct casereader *r1 =
 345     casereader_create_distinct (sort_execute_1var (reader, var), var, w);
 346
 347   struct ccase *c1;
 348
 349   struct casereader *rclone = casereader_clone (r1);
 350   struct casewriter *wtr;
 351   struct caseproto *proto = caseproto_create ();
 352
 353   proto = caseproto_add_width (proto, 0);
 354   proto = caseproto_add_width (proto, 0);
 355   proto = caseproto_add_width (proto, 0);
 356
 357   wtr = autopaging_writer_create (proto);
 358
 359   *cc = 0;
 360
 361   for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
 362     {
 363       struct ccase *c2;
 364       struct casereader *r2 = casereader_clone (rclone);
 365
 366       const double weight1 = case_data_idx (c1, weight_idx)->f;
 367       const double d1 = case_data (c1, var)->f;
 368       double n_eq = 0.0;
 369       double n_pred = 0.0;
 370
 371
 372       struct ccase *new_case = case_create (proto);
 373
 374       *cc += weight1;
 375
 376       for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
 377         {
 378           const double d2 = case_data (c2, var)->f;
 379           const double weight2 = case_data_idx (c2, weight_idx)->f;
 380
 381           if ( d1 == d2 )
 382             {
 383               n_eq += weight2;
 384               continue;
 385             }
 386           else  if ( pred (d2, d1))
 387             {
 388               n_pred += weight2;
 389             }
 390         }
 391
 392       case_data_rw_idx (new_case, VALUE)->f = d1;
 393       case_data_rw_idx (new_case, N_EQ)->f = n_eq;
 394       case_data_rw_idx (new_case, N_PRED)->f = n_pred;
 395
 396       casewriter_write (wtr, new_case);
 397
 398       casereader_destroy (r2);
 399     }
 400
 401   casereader_destroy (r1);
 402   casereader_destroy (rclone);
 403
 404   return casewriter_make_reader (wtr);
 405 }
 406
 407 static bool
 408 gt (double d1, double d2)
 409 {
 410   return d1 > d2;
 411 }
 412
 413 static bool
 414 lt (double d1, double d2)
 415 {
 416   return d1 < d2;
 417 }
 418
 419
 420 static void
 421 do_roc (struct cmd_roc *roc, struct casereader *input, struct dictionary *dict)
 422 {
 423   int i;
 424
 425   struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
 426
 427   const struct caseproto *proto = casereader_get_proto (input);
 428
 429   struct casewriter *neg_wtr = autopaging_writer_create (proto);
 430
 431   struct casereader *negatives = NULL;
 432
 433   struct casereader *positives =
 434     casereader_create_filter_func (input,
 435                                    match_positives,
 436                                    NULL,
 437                                    roc,
 438                                    neg_wtr);
 439
 440
 441   for (i = 0 ; i < roc->n_vars; ++i)
 442     {
 443       double q1hat = 0;
 444       double q2hat = 0;
 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               q1hat += n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
 488               q2hat += n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
 489             }
 490
 491           if ( cneg )
 492             case_unref (cneg);
 493         }
 494
 495       rs[i].auc /=  rs[i].n1 * rs[i].n2;
 496       if ( roc->invert )
 497         rs[i].auc = 1 - rs[i].auc;
 498
 499
 500     }
 501
 502   casereader_destroy (positives);
 503   casereader_destroy (negatives);
 504
 505   output_roc (rs, roc);
 506
 507   free (rs);
 508 }
 509
 510
 511
 512
 513 static void
 514 show_auc  (struct roc_state *rs, const struct cmd_roc *roc)
 515 {
 516   int i;
 517   const int n_fields = roc->print_se ? 5 : 1;
 518   const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields;
 519   const int n_rows = 2 + roc->n_vars;
 520   struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
 521
 522   if ( roc->n_vars > 1)
 523     tab_title (tbl, _("Area Under the Curve"));
 524   else
 525     tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
 526
 527   tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
 528
 529   tab_dim (tbl, tab_natural_dimensions, NULL);
 530
 531   tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
 532
 533   tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
 534
 535   tab_box (tbl,
 536            TAL_2, TAL_2,
 537            -1, TAL_1,
 538            0, 0,
 539            n_cols - 1,
 540            n_rows - 1);
 541
 542   if ( roc->print_se )
 543     {
 544       tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
 545       tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
 546
 547       tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
 548       tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
 549
 550       tab_joint_text (tbl, n_cols - 2, 0, 4, 0,
 551                       TAT_TITLE | TAB_CENTER | TAT_PRINTF,
 552                       _("Asymp. %g%% Confidence Interval"), roc->ci);
 553       tab_vline (tbl, 0, n_cols - 1, 0, 0);
 554       tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1);
 555     }
 556
 557   if ( roc->n_vars > 1)
 558     tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
 559
 560   if ( roc->n_vars > 1)
 561     tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
 562
 563
 564   for ( i = 0 ; i < roc->n_vars ; ++i )
 565     {
 566       tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
 567
 568       tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
 569
 570       if ( roc->print_se )
 571         {
 572
 573           double se ;
 574           const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
 575                                       (12 * rs[i].n1 * rs[i].n2));
 576           double ci ;
 577           double yy ;
 578
 579           double q1 = rs[i].auc / ( 2 - rs[i].auc);
 580           double q2 = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
 581
 582           se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (q1 - pow2 (rs[i].auc)) +
 583             (rs[i].n2 - 1) * (q2 - pow2 (rs[i].auc));
 584
 585           se /= rs[i].n1 * rs[i].n2;
 586
 587           se = sqrt (se);
 588
 589           tab_double (tbl, n_cols - 4, 2 + i, 0,
 590                       se,
 591                       NULL);
 592
 593           ci = 1 - roc->ci / 100.0;
 594           yy = gsl_cdf_gaussian_Qinv (ci, se) ;
 595
 596           tab_double (tbl, n_cols - 2, 2 + i, 0,
 597                       rs[i].auc - yy,
 598                       NULL);
 599
 600           tab_double (tbl, n_cols - 1, 2 + i, 0,
 601                       rs[i].auc + yy,
 602                       NULL);
 603
 604           tab_double (tbl, n_cols - 3, 2 + i, 0,
 605                       2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)),
 606                       NULL);
 607         }
 608     }
 609
 610   tab_submit (tbl);
 611 }
 612
 613
 614
 615
 616 static void
 617 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
 618 {
 619 #if 0
 620   show_summary (roc);
 621
 622   if ( roc->curve )
 623     draw_roc (rs, roc);
 624 #endif
 625
 626   show_auc (rs, roc);
 627
 628 #if 0
 629   if ( roc->print_coords )
 630     show_coords (rs, roc);
 631 #endif
 632 }