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