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