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 <data/format.h>
  30 #include <math/sort.h>
  31 #include <data/subcase.h>
  32
  33
  34 #include <libpspp/misc.h>
  35
  36 #include <gsl/gsl_cdf.h>
  37 #include <output/table.h>
  38
  39 #include <output/charts/plot-chart.h>
  40 #include <output/charts/cartesian.h>
  41
  42 #include "gettext.h"
  43 #define _(msgid) gettext (msgid)
  44 #define N_(msgid) msgid
  45
  46 struct cmd_roc
  47 {
  48   size_t n_vars;
  49   const struct variable **vars;
  50
  51   struct variable *state_var ;
  52   union value state_value;
  53
  54   /* Plot the roc curve */
  55   bool curve;
  56   /* Plot the reference line */
  57   bool reference;
  58
  59   double ci;
  60
  61   bool print_coords;
  62   bool print_se;
  63   bool bi_neg_exp; /* True iff the bi-negative exponential critieria
  64                       should be used */
  65   enum mv_class exclude;
  66
  67   bool invert ; /* True iff a smaller test result variable indicates
  68                    a positive result */
  69 };
  70
  71 static int run_roc (struct dataset *ds, struct cmd_roc *roc);
  72
  73 int
  74 cmd_roc (struct lexer *lexer, struct dataset *ds)
  75 {
  76   struct cmd_roc roc ;
  77   const struct dictionary *dict = dataset_dict (ds);
  78
  79   roc.vars = NULL;
  80   roc.n_vars = 0;
  81   roc.print_se = false;
  82   roc.print_coords = false;
  83   roc.exclude = MV_ANY;
  84   roc.curve = true;
  85   roc.reference = false;
  86   roc.ci = 95;
  87   roc.bi_neg_exp = false;
  88   roc.invert = false;
  89
  90   if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars,
  91                               PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
  92     return 2;
  93
  94   if ( ! lex_force_match (lexer, T_BY))
  95     {
  96       return 2;
  97     }
  98
  99   roc.state_var = parse_variable (lexer, dict);
 100
 101   if ( !lex_force_match (lexer, '('))
 102     {
 103       return 2;
 104     }
 105
 106   parse_value (lexer, &roc.state_value, var_get_width (roc.state_var));
 107
 108
 109   if ( !lex_force_match (lexer, ')'))
 110     {
 111       return 2;
 112     }
 113
 114
 115   while (lex_token (lexer) != '.')
 116     {
 117       lex_match (lexer, '/');
 118       if (lex_match_id (lexer, "MISSING"))
 119         {
 120           lex_match (lexer, '=');
 121           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 122             {
 123               if (lex_match_id (lexer, "INCLUDE"))
 124                 {
 125                   roc.exclude = MV_SYSTEM;
 126                 }
 127               else if (lex_match_id (lexer, "EXCLUDE"))
 128                 {
 129                   roc.exclude = MV_ANY;
 130                 }
 131               else
 132                 {
 133                   lex_error (lexer, NULL);
 134                   return 2;
 135                 }
 136             }
 137         }
 138       else if (lex_match_id (lexer, "PLOT"))
 139         {
 140           lex_match (lexer, '=');
 141           if (lex_match_id (lexer, "CURVE"))
 142             {
 143               roc.curve = true;
 144               if (lex_match (lexer, '('))
 145                 {
 146                   roc.reference = true;
 147                   lex_force_match_id (lexer, "REFERENCE");
 148                   lex_force_match (lexer, ')');
 149                 }
 150             }
 151           else if (lex_match_id (lexer, "NONE"))
 152             {
 153               roc.curve = false;
 154             }
 155           else
 156             {
 157               lex_error (lexer, NULL);
 158               return 2;
 159             }
 160         }
 161       else if (lex_match_id (lexer, "PRINT"))
 162         {
 163           lex_match (lexer, '=');
 164           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 165             {
 166               if (lex_match_id (lexer, "SE"))
 167                 {
 168                   roc.print_se = true;
 169                 }
 170               else if (lex_match_id (lexer, "COORDINATES"))
 171                 {
 172                   roc.print_coords = true;
 173                 }
 174               else
 175                 {
 176                   lex_error (lexer, NULL);
 177                   return 2;
 178                 }
 179             }
 180         }
 181       else if (lex_match_id (lexer, "CRITERIA"))
 182         {
 183           lex_match (lexer, '=');
 184           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 185             {
 186               if (lex_match_id (lexer, "CUTOFF"))
 187                 {
 188                   lex_force_match (lexer, '(');
 189                   if (lex_match_id (lexer, "INCLUDE"))
 190                     {
 191                       roc.exclude = MV_SYSTEM;
 192                     }
 193                   else if (lex_match_id (lexer, "EXCLUDE"))
 194                     {
 195                       roc.exclude = MV_USER | MV_SYSTEM;
 196                     }
 197                   else
 198                     {
 199                       lex_error (lexer, NULL);
 200                       return 2;
 201                     }
 202                   lex_force_match (lexer, ')');
 203                 }
 204               else if (lex_match_id (lexer, "TESTPOS"))
 205                 {
 206                   lex_force_match (lexer, '(');
 207                   if (lex_match_id (lexer, "LARGE"))
 208                     {
 209                       roc.invert = false;
 210                     }
 211                   else if (lex_match_id (lexer, "SMALL"))
 212                     {
 213                       roc.invert = true;
 214                     }
 215                   else
 216                     {
 217                       lex_error (lexer, NULL);
 218                       return 2;
 219                     }
 220                   lex_force_match (lexer, ')');
 221                 }
 222               else if (lex_match_id (lexer, "CI"))
 223                 {
 224                   lex_force_match (lexer, '(');
 225                   lex_force_num (lexer);
 226                   roc.ci = lex_number (lexer);
 227                   lex_get (lexer);
 228                   lex_force_match (lexer, ')');
 229                 }
 230               else if (lex_match_id (lexer, "DISTRIBUTION"))
 231                 {
 232                   lex_force_match (lexer, '(');
 233                   if (lex_match_id (lexer, "FREE"))
 234                     {
 235                       roc.bi_neg_exp = false;
 236                     }
 237                   else if (lex_match_id (lexer, "NEGEXPO"))
 238                     {
 239                       roc.bi_neg_exp = true;
 240                     }
 241                   else
 242                     {
 243                       lex_error (lexer, NULL);
 244                       return 2;
 245                     }
 246                   lex_force_match (lexer, ')');
 247                 }
 248               else
 249                 {
 250                   lex_error (lexer, NULL);
 251                   return 2;
 252                 }
 253             }
 254         }
 255       else
 256         {
 257           lex_error (lexer, NULL);
 258           break;
 259         }
 260     }
 261
 262   run_roc (ds, &roc);
 263
 264   return 1;
 265 }
 266
 267
 268
 269
 270 static void
 271 do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict);
 272
 273
 274 static int
 275 run_roc (struct dataset *ds, struct cmd_roc *roc)
 276 {
 277   struct dictionary *dict = dataset_dict (ds);
 278   bool ok;
 279   struct casereader *group;
 280
 281   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
 282   while (casegrouper_get_next_group (grouper, &group))
 283     {
 284       do_roc (roc, group, dataset_dict (ds));
 285     }
 286   ok = casegrouper_destroy (grouper);
 287   ok = proc_commit (ds) && ok;
 288
 289   return ok;
 290 }
 291
 292 #if 0
 293 static void
 294 dump_casereader (struct casereader *reader)
 295 {
 296   struct ccase *c;
 297   struct casereader *r = casereader_clone (reader);
 298
 299   for ( ; (c = casereader_read (r) ); case_unref (c))
 300     {
 301       int i;
 302       for (i = 0 ; i < case_get_value_cnt (c); ++i)
 303         {
 304           printf ("%g ", case_data_idx (c, i)->f);
 305         }
 306       printf ("\n");
 307     }
 308
 309   casereader_destroy (r);
 310 }
 311 #endif
 312
 313 static bool
 314 match_positives (const struct ccase *c, void *aux)
 315 {
 316   struct cmd_roc *roc = aux;
 317
 318   return 0 == value_compare_3way (case_data (c, roc->state_var),
 319                                   &roc->state_value,
 320                                   var_get_width (roc->state_var));
 321 }
 322
 323
 324 #define VALUE  0
 325 #define N_EQ   1
 326 #define N_PRED 2
 327
 328 struct roc_state
 329 {
 330   double auc;
 331
 332   double n1;
 333   double n2;
 334
 335   double q1hat;
 336   double q2hat;
 337
 338   struct casewriter *cutpoint_wtr;
 339   struct casereader *cutpoint_rdr;
 340   double prev_result;
 341   double min;
 342   double max;
 343 };
 344
 345
 346
 347 #define CUTPOINT 0
 348 #define TP 1
 349 #define FN 2
 350 #define TN 3
 351 #define FP 4
 352
 353
 354 static struct casereader *
 355 accumulate_counts (struct casereader *cutpoint_rdr,
 356                    double result, double weight,
 357                    bool (*pos_cond) (double, double),
 358                    int true_index, int false_index)
 359 {
 360   const struct caseproto *proto = casereader_get_proto (cutpoint_rdr);
 361   struct casewriter *w =
 362     autopaging_writer_create (proto);
 363   struct casereader *r = casereader_clone (cutpoint_rdr);
 364   struct ccase *cpc;
 365   double prev_cp = SYSMIS;
 366
 367
 368   for ( ; (cpc = casereader_read (r) ); case_unref (cpc))
 369     {
 370       struct ccase *new_case;
 371       const double cp = case_data_idx (cpc, CUTPOINT)->f;
 372
 373       /* We don't want duplicates here */
 374       if ( cp == prev_cp )
 375         continue;
 376
 377       new_case = case_clone (cpc);
 378
 379       if ( pos_cond (result, cp))
 380         {
 381           case_data_rw_idx (new_case, true_index)->f += weight;
 382         }
 383       else
 384         {
 385           case_data_rw_idx (new_case, false_index)->f += weight;
 386         }
 387
 388       prev_cp = cp;
 389
 390       casewriter_write (w, new_case);
 391     }
 392   casereader_destroy (r);
 393
 394   return casewriter_make_reader (w);
 395 }
 396
 397
 398
 399 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
 400
 401
 402 static struct casereader *
 403 process_group (const struct variable *var, struct casereader *reader,
 404                bool (*pred) (double, double),
 405                const struct dictionary *dict,
 406                double *cc,
 407                struct casereader **cutpoint_rdr,
 408                bool (*pos_cond) (double, double),
 409                int true_index,
 410                int false_index
 411                )
 412 {
 413   const struct variable *w = dict_get_weight (dict);
 414   struct casereader *r1 =
 415     casereader_create_distinct (sort_execute_1var (reader, var), var, w);
 416
 417   const int weight_idx  = w ? var_get_case_index (w) :
 418     caseproto_get_n_widths (casereader_get_proto (r1)) - 1;
 419
 420   struct ccase *c1;
 421
 422   struct casereader *rclone = casereader_clone (r1);
 423   struct casewriter *wtr;
 424   struct caseproto *proto = caseproto_create ();
 425
 426   proto = caseproto_add_width (proto, 0);
 427   proto = caseproto_add_width (proto, 0);
 428   proto = caseproto_add_width (proto, 0);
 429
 430   wtr = autopaging_writer_create (proto);
 431
 432   *cc = 0;
 433
 434   for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
 435     {
 436       struct ccase *new_case = case_create (proto);
 437       struct ccase *c2;
 438       struct casereader *r2 = casereader_clone (rclone);
 439
 440       const double weight1 = case_data_idx (c1, weight_idx)->f;
 441       const double d1 = case_data (c1, var)->f;
 442       double n_eq = 0.0;
 443       double n_pred = 0.0;
 444
 445       *cutpoint_rdr = accumulate_counts (*cutpoint_rdr, d1, weight1,
 446                                          pos_cond,
 447                                          true_index, false_index);
 448
 449       *cc += weight1;
 450
 451       for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
 452         {
 453           const double d2 = case_data (c2, var)->f;
 454           const double weight2 = case_data_idx (c2, weight_idx)->f;
 455
 456           if ( d1 == d2 )
 457             {
 458               n_eq += weight2;
 459               continue;
 460             }
 461           else  if ( pred (d2, d1))
 462             {
 463               n_pred += weight2;
 464             }
 465         }
 466
 467       case_data_rw_idx (new_case, VALUE)->f = d1;
 468       case_data_rw_idx (new_case, N_EQ)->f = n_eq;
 469       case_data_rw_idx (new_case, N_PRED)->f = n_pred;
 470
 471       casewriter_write (wtr, new_case);
 472
 473       casereader_destroy (r2);
 474     }
 475
 476   casereader_destroy (r1);
 477   casereader_destroy (rclone);
 478
 479   return casewriter_make_reader (wtr);
 480 }
 481
 482 static bool
 483 gt (double d1, double d2)
 484 {
 485   return d1 > d2;
 486 }
 487
 488
 489 static bool
 490 ge (double d1, double d2)
 491 {
 492   return d1 > d2;
 493 }
 494
 495 static bool
 496 lt (double d1, double d2)
 497 {
 498   return d1 < d2;
 499 }
 500
 501 static struct casereader *
 502 process_positive_group (const struct variable *var, struct casereader *reader,
 503                         const struct dictionary *dict,
 504                         struct roc_state *rs)
 505 {
 506   return process_group (var, reader, gt, dict, &rs->n1,
 507                         &rs->cutpoint_rdr,
 508                         ge,
 509                         TP, FN);
 510 }
 511
 512
 513 static struct casereader *
 514 process_negative_group (const struct variable *var, struct casereader *reader,
 515                         const struct dictionary *dict,
 516                         struct roc_state *rs)
 517 {
 518   return process_group (var, reader, lt, dict, &rs->n2,
 519                         &rs->cutpoint_rdr,
 520                         lt,
 521                         TN, FP);
 522 }
 523
 524
 525
 526
 527 static void
 528 append_cutpoint (struct casewriter *writer, double cutpoint)
 529 {
 530   struct ccase *cc = case_create (casewriter_get_proto (writer));
 531
 532   case_data_rw_idx (cc, CUTPOINT)->f = cutpoint;
 533   case_data_rw_idx (cc, TP)->f = 0;
 534   case_data_rw_idx (cc, FN)->f = 0;
 535   case_data_rw_idx (cc, TN)->f = 0;
 536   case_data_rw_idx (cc, FP)->f = 0;
 537
 538
 539   casewriter_write (writer, cc);
 540 }
 541
 542
 543 /* Prepare the cutpoints */
 544 static void
 545 prepare_cutpoints (struct cmd_roc *roc, struct roc_state *rs, struct casereader *input)
 546 {
 547   int i;
 548   struct casereader *r = casereader_clone (input);
 549   struct ccase *c;
 550   struct caseproto *proto = caseproto_create ();
 551
 552   struct subcase ordering;
 553   subcase_init (&ordering, CUTPOINT, 0, SC_ASCEND);
 554
 555   proto = caseproto_add_width (proto, 0); /* cutpoint */
 556   proto = caseproto_add_width (proto, 0); /* TP */
 557   proto = caseproto_add_width (proto, 0); /* FN */
 558   proto = caseproto_add_width (proto, 0); /* TN */
 559   proto = caseproto_add_width (proto, 0); /* FP */
 560
 561   for (i = 0 ; i < roc->n_vars; ++i)
 562     {
 563       rs[i].cutpoint_wtr = sort_create_writer (&ordering, proto);
 564       rs[i].prev_result = SYSMIS;
 565       rs[i].max = -DBL_MAX;
 566       rs[i].min = DBL_MAX;
 567     }
 568
 569   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 570     {
 571       for (i = 0 ; i < roc->n_vars; ++i)
 572         {
 573           const double result = case_data (c, roc->vars[i])->f;
 574
 575           minimize (&rs[i].min, result);
 576           maximize (&rs[i].max, result);
 577
 578           if ( rs[i].prev_result != SYSMIS && rs[i].prev_result != result )
 579             {
 580               const double mean = (result + rs[i].prev_result ) / 2.0;
 581               append_cutpoint (rs[i].cutpoint_wtr, mean);
 582             }
 583
 584           rs[i].prev_result = result;
 585         }
 586     }
 587   casereader_destroy (r);
 588
 589
 590   /* Append the min and max cutpoints */
 591   for (i = 0 ; i < roc->n_vars; ++i)
 592     {
 593       append_cutpoint (rs[i].cutpoint_wtr, rs[i].min - 1);
 594       append_cutpoint (rs[i].cutpoint_wtr, rs[i].max + 1);
 595
 596       rs[i].cutpoint_rdr = casewriter_make_reader (rs[i].cutpoint_wtr);
 597     }
 598 }
 599
 600 static void
 601 do_roc (struct cmd_roc *roc, struct casereader *input, struct dictionary *dict)
 602 {
 603   int i;
 604
 605   struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
 606
 607   struct casewriter *neg_wtr = autopaging_writer_create (casereader_get_proto (input));
 608
 609   struct casereader *negatives = NULL;
 610   struct casereader *positives = NULL;
 611
 612   struct caseproto *n_proto = caseproto_create ();
 613
 614   struct subcase up_ordering;
 615   struct subcase down_ordering;
 616
 617   prepare_cutpoints (roc, rs, input);
 618
 619   positives =
 620     casereader_create_filter_func (input,
 621                                    match_positives,
 622                                    NULL,
 623                                    roc,
 624                                    neg_wtr);
 625
 626   n_proto = caseproto_create ();
 627
 628   n_proto = caseproto_add_width (n_proto, 0);
 629   n_proto = caseproto_add_width (n_proto, 0);
 630   n_proto = caseproto_add_width (n_proto, 0);
 631   n_proto = caseproto_add_width (n_proto, 0);
 632   n_proto = caseproto_add_width (n_proto, 0);
 633
 634   subcase_init (&up_ordering, VALUE, 0, SC_ASCEND);
 635   subcase_init (&down_ordering, VALUE, 0, SC_DESCEND);
 636
 637   for (i = 0 ; i < roc->n_vars; ++i)
 638     {
 639       struct casewriter *w = NULL;
 640       struct casereader *r = NULL;
 641
 642       struct ccase *c;
 643
 644       struct ccase *cpos;
 645       struct casereader *n_neg ;
 646       const struct variable *var = roc->vars[i];
 647
 648       struct casereader *neg ;
 649       struct casereader *pos = casereader_clone (positives);
 650
 651       struct casereader *n_pos =
 652         process_positive_group (var, pos, dict, &rs[i]);
 653
 654       if ( negatives == NULL)
 655         {
 656           negatives = casewriter_make_reader (neg_wtr);
 657         }
 658
 659       neg = casereader_clone (negatives);
 660
 661       n_neg = process_negative_group (var, neg, dict, &rs[i]);
 662
 663       w = sort_create_writer (&up_ordering, n_proto);
 664       for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos))
 665         {
 666           struct ccase *pos_case = case_create (n_proto);
 667           struct ccase *cneg;
 668           const double jpos = case_data_idx (cpos, VALUE)->f;
 669
 670           while ((cneg = casereader_read (n_neg)))
 671             {
 672               struct ccase *nc = case_create (n_proto);
 673
 674               const double jneg = case_data_idx (cneg, VALUE)->f;
 675
 676               case_data_rw_idx (nc, VALUE)->f = jneg;
 677               case_data_rw_idx (nc, N_EQ)->f = 0;
 678
 679               case_data_rw_idx (nc, N_PRED)->f = SYSMIS;
 680
 681               *case_data_rw_idx (nc, 3) = *case_data_idx (cneg, N_EQ);
 682               *case_data_rw_idx (nc, 4) = *case_data_idx (cneg, N_PRED);
 683
 684               casewriter_write (w, nc);
 685
 686               case_unref (cneg);
 687               if ( jneg > jpos)
 688                 break;
 689             }
 690
 691           case_data_rw_idx (pos_case, VALUE)->f = jpos;
 692           *case_data_rw_idx (pos_case, N_EQ) = *case_data_idx (cpos, N_EQ);
 693           *case_data_rw_idx (pos_case, N_PRED) = *case_data_idx (cpos, N_PRED);
 694           case_data_rw_idx (pos_case, 3)->f = 0;
 695           case_data_rw_idx (pos_case, 4)->f = SYSMIS;
 696
 697           casewriter_write (w, pos_case);
 698         }
 699
 700       r = casewriter_make_reader (w);
 701
 702       {
 703         double prev_pos_gt = rs[i].n1;
 704         w = sort_create_writer (&down_ordering, n_proto);
 705
 706         for ( ; (c = casereader_read (r) ); case_unref (c))
 707           {
 708             double n_pos_gt = case_data_idx (c, N_PRED)->f;
 709             struct ccase *nc = case_clone (c);
 710
 711             if ( n_pos_gt == SYSMIS)
 712               {
 713                 n_pos_gt = prev_pos_gt;
 714                 case_data_rw_idx (nc, N_PRED)->f = n_pos_gt;
 715               }
 716
 717             casewriter_write (w, nc);
 718             prev_pos_gt = n_pos_gt;
 719           }
 720
 721         r = casewriter_make_reader (w);
 722       }
 723
 724       {
 725         double prev_neg_lt = rs[i].n2;
 726         w = sort_create_writer (&up_ordering, n_proto);
 727
 728         for ( ; (c = casereader_read (r) ); case_unref (c))
 729           {
 730             double n_neg_lt = case_data_idx (c, 4)->f;
 731             struct ccase *nc = case_clone (c);
 732
 733             if ( n_neg_lt == SYSMIS)
 734               {
 735                 n_neg_lt = prev_neg_lt;
 736                 case_data_rw_idx (nc, 4)->f = n_neg_lt;
 737               }
 738
 739             casewriter_write (w, nc);
 740             prev_neg_lt = n_neg_lt;
 741           }
 742
 743         r = casewriter_make_reader (w);
 744       }
 745
 746       {
 747         struct ccase *prev_case = NULL;
 748         for ( ; (c = casereader_read (r) ); case_unref (c))
 749           {
 750             const struct ccase *next_case = casereader_peek (r, 0);
 751
 752             const double j = case_data_idx (c, VALUE)->f;
 753             double n_pos_eq = case_data_idx (c, N_EQ)->f;
 754             double n_pos_gt = case_data_idx (c, N_PRED)->f;
 755             double n_neg_eq = case_data_idx (c, 3)->f;
 756             double n_neg_lt = case_data_idx (c, 4)->f;
 757
 758             if ( prev_case && j == case_data_idx (prev_case, VALUE)->f)
 759               {
 760                 if ( 0 ==  case_data_idx (c, N_EQ)->f)
 761                   {
 762                     n_pos_eq = case_data_idx (prev_case, N_EQ)->f;
 763                     n_pos_gt = case_data_idx (prev_case, N_PRED)->f;
 764                   }
 765
 766                 if ( 0 ==  case_data_idx (c, 3)->f)
 767                   {
 768                     n_neg_eq = case_data_idx (prev_case, 3)->f;
 769                     n_neg_lt = case_data_idx (prev_case, 4)->f;
 770                   }
 771               }
 772
 773             if ( NULL == next_case || j != case_data_idx (next_case, VALUE)->f)
 774               {
 775                 rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
 776
 777                 rs[i].q1hat +=
 778                   n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
 779                 rs[i].q2hat +=
 780                   n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
 781
 782               }
 783
 784             case_unref (prev_case);
 785             prev_case = case_clone (c);
 786           }
 787
 788         rs[i].auc /=  rs[i].n1 * rs[i].n2;
 789         if ( roc->invert )
 790           rs[i].auc = 1 - rs[i].auc;
 791
 792         if ( roc->bi_neg_exp )
 793           {
 794             rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
 795             rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
 796           }
 797         else
 798           {
 799             rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
 800             rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
 801           }
 802       }
 803     }
 804
 805   casereader_destroy (positives);
 806   casereader_destroy (negatives);
 807
 808   output_roc (rs, roc);
 809
 810   free (rs);
 811 }
 812
 813 static void
 814 show_auc  (struct roc_state *rs, const struct cmd_roc *roc)
 815 {
 816   int i;
 817   const int n_fields = roc->print_se ? 5 : 1;
 818   const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields;
 819   const int n_rows = 2 + roc->n_vars;
 820   struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
 821
 822   if ( roc->n_vars > 1)
 823     tab_title (tbl, _("Area Under the Curve"));
 824   else
 825     tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
 826
 827   tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
 828
 829   tab_dim (tbl, tab_natural_dimensions, NULL);
 830
 831   tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
 832
 833   tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
 834
 835   tab_box (tbl,
 836            TAL_2, TAL_2,
 837            -1, TAL_1,
 838            0, 0,
 839            n_cols - 1,
 840            n_rows - 1);
 841
 842   if ( roc->print_se )
 843     {
 844       tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
 845       tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
 846
 847       tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
 848       tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
 849
 850       tab_joint_text (tbl, n_cols - 2, 0, 4, 0,
 851                       TAT_TITLE | TAB_CENTER | TAT_PRINTF,
 852                       _("Asymp. %g%% Confidence Interval"), roc->ci);
 853       tab_vline (tbl, 0, n_cols - 1, 0, 0);
 854       tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1);
 855     }
 856
 857   if ( roc->n_vars > 1)
 858     tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
 859
 860   if ( roc->n_vars > 1)
 861     tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
 862
 863
 864   for ( i = 0 ; i < roc->n_vars ; ++i )
 865     {
 866       tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
 867
 868       tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
 869
 870       if ( roc->print_se )
 871         {
 872           double se ;
 873           const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
 874                                       (12 * rs[i].n1 * rs[i].n2));
 875           double ci ;
 876           double yy ;
 877
 878           se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (rs[i].q1hat - pow2 (rs[i].auc)) +
 879             (rs[i].n2 - 1) * (rs[i].q2hat - pow2 (rs[i].auc));
 880
 881           se /= rs[i].n1 * rs[i].n2;
 882
 883           se = sqrt (se);
 884
 885           tab_double (tbl, n_cols - 4, 2 + i, 0,
 886                       se,
 887                       NULL);
 888
 889           ci = 1 - roc->ci / 100.0;
 890           yy = gsl_cdf_gaussian_Qinv (ci, se) ;
 891
 892           tab_double (tbl, n_cols - 2, 2 + i, 0,
 893                       rs[i].auc - yy,
 894                       NULL);
 895
 896           tab_double (tbl, n_cols - 1, 2 + i, 0,
 897                       rs[i].auc + yy,
 898                       NULL);
 899
 900           tab_double (tbl, n_cols - 3, 2 + i, 0,
 901                       2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)),
 902                       NULL);
 903         }
 904     }
 905
 906   tab_submit (tbl);
 907 }
 908
 909
 910 static void
 911 show_summary (const struct cmd_roc *roc)
 912 {
 913   const int n_cols = 3;
 914   const int n_rows = 4;
 915   struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
 916
 917   tab_title (tbl, _("Case Summary"));
 918
 919   tab_headers (tbl, 1, 0, 2, 0);
 920
 921   tab_dim (tbl, tab_natural_dimensions, NULL);
 922
 923   tab_box (tbl,
 924            TAL_2, TAL_2,
 925            -1, -1,
 926            0, 0,
 927            n_cols - 1,
 928            n_rows - 1);
 929
 930   tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
 931   tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
 932
 933
 934   tab_hline (tbl, TAL_2, 1, n_cols - 1, 1);
 935   tab_vline (tbl, TAL_1, 2, 1, n_rows - 1);
 936
 937
 938   tab_text (tbl, 0, 1, TAT_TITLE | TAB_LEFT, var_to_string (roc->state_var));
 939   tab_text (tbl, 1, 1, TAT_TITLE, _("Unweighted"));
 940   tab_text (tbl, 2, 1, TAT_TITLE, _("Weighted"));
 941
 942   tab_joint_text (tbl, 1, 0, 2, 0,
 943                   TAT_TITLE | TAB_CENTER,
 944                   _("Valid N (listwise)"));
 945
 946
 947   tab_text (tbl, 0, 2, TAB_LEFT, _("Positive"));
 948   tab_text (tbl, 0, 3, TAB_LEFT, _("Negative"));
 949
 950
 951 #if 0
 952   tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0);
 953   tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0);
 954
 955   tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0);
 956   tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0);
 957 #endif
 958
 959   tab_submit (tbl);
 960 }
 961
 962
 963 static void
 964 show_coords (struct roc_state *rs, const struct cmd_roc *roc)
 965 {
 966   int x = 1;
 967   int i;
 968   const int n_cols = roc->n_vars > 1 ? 4 : 3;
 969   int n_rows = 1;
 970   struct tab_table *tbl ;
 971
 972   for (i = 0; i < roc->n_vars; ++i)
 973     n_rows += casereader_count_cases (rs[i].cutpoint_rdr);
 974
 975   tbl = tab_create (n_cols, n_rows, 0);
 976
 977   if ( roc->n_vars > 1)
 978     tab_title (tbl, _("Coordinates of the Curve"));
 979   else
 980     tab_title (tbl, _("Coordinates of the Curve (%s)"), var_to_string (roc->vars[0]));
 981
 982
 983   tab_headers (tbl, 1, 0, 1, 0);
 984
 985   tab_dim (tbl, tab_natural_dimensions, NULL);
 986
 987   tab_hline (tbl, TAL_2, 0, n_cols - 1, 1);
 988
 989   if ( roc->n_vars > 1)
 990     tab_text (tbl, 0, 0, TAT_TITLE, _("Test variable"));
 991
 992   tab_text (tbl, n_cols - 3, 0, TAT_TITLE, _("Positive if greater than or equal to"));
 993   tab_text (tbl, n_cols - 2, 0, TAT_TITLE, _("Sensitivity"));
 994   tab_text (tbl, n_cols - 1, 0, TAT_TITLE, _("1 - Specificity"));
 995
 996   tab_box (tbl,
 997            TAL_2, TAL_2,
 998            -1, TAL_1,
 999            0, 0,
1000            n_cols - 1,
1001            n_rows - 1);
1002
1003   if ( roc->n_vars > 1)
1004     tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
1005
1006   for (i = 0; i < roc->n_vars; ++i)
1007     {
1008       struct ccase *cc;
1009       struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
1010
1011       if ( roc->n_vars > 1)
1012         tab_text (tbl, 0, x, TAT_TITLE, var_to_string (roc->vars[i]));
1013
1014       if ( i > 0)
1015         tab_hline (tbl, TAL_1, 0, n_cols - 1, x);
1016
1017
1018       for (; (cc = casereader_read (r)) != NULL;
1019            case_unref (cc), x++)
1020         {
1021           const double se = case_data_idx (cc, TP)->f /
1022             (
1023              case_data_idx (cc, TP)->f
1024              +
1025              case_data_idx (cc, FN)->f
1026              );
1027
1028           const double sp = case_data_idx (cc, TN)->f /
1029             (
1030              case_data_idx (cc, TN)->f
1031              +
1032              case_data_idx (cc, FP)->f
1033              );
1034
1035           tab_double (tbl, n_cols - 3, x, 0, case_data_idx (cc, CUTPOINT)->f,
1036                       var_get_print_format (roc->vars[i]));
1037
1038           tab_double (tbl, n_cols - 2, x, 0, se, NULL);
1039           tab_double (tbl, n_cols - 1, x, 0, 1 - sp, NULL);
1040         }
1041
1042       casereader_destroy (r);
1043     }
1044
1045   tab_submit (tbl);
1046 }
1047
1048
1049 static void
1050 draw_roc (struct roc_state *rs, const struct cmd_roc *roc)
1051 {
1052   int i;
1053
1054   struct chart *roc_chart = chart_create ();
1055
1056   chart_write_title (roc_chart, _("ROC Curve"));
1057   chart_write_xlabel (roc_chart, _("1 - Specificity"));
1058   chart_write_ylabel (roc_chart, _("Sensitivity"));
1059
1060   chart_write_xscale (roc_chart, 0, 1, 5);
1061   chart_write_yscale (roc_chart, 0, 1, 5);
1062
1063   if ( roc->reference )
1064     {
1065       chart_line (roc_chart, 1.0, 0,
1066                   0.0, 1.0,
1067                   CHART_DIM_X);
1068     }
1069
1070   for (i = 0; i < roc->n_vars; ++i)
1071     {
1072       struct ccase *cc;
1073       struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
1074
1075       chart_vector_start (roc_chart, var_get_name (roc->vars[i]));
1076       for (; (cc = casereader_read (r)) != NULL;
1077            case_unref (cc))
1078         {
1079           double se = case_data_idx (cc, TP)->f;
1080           double sp = case_data_idx (cc, TN)->f;
1081
1082           se /= case_data_idx (cc, FN)->f +
1083             case_data_idx (cc, TP)->f ;
1084
1085           sp /= case_data_idx (cc, TN)->f +
1086             case_data_idx (cc, FP)->f ;
1087
1088           chart_vector (roc_chart, 1 - sp, se);
1089         }
1090       chart_vector_end (roc_chart);
1091       casereader_destroy (r);
1092     }
1093
1094   chart_write_legend (roc_chart);
1095
1096   chart_submit (roc_chart);
1097 }
1098
1099
1100 static void
1101 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
1102 {
1103   show_summary (roc);
1104
1105   if ( roc->curve )
1106     draw_roc (rs, roc);
1107
1108   show_auc (rs, roc);
1109
1110
1111   if ( roc->print_coords )
1112     show_coords (rs, roc);
1113 }
1114