1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2009 Free Software Foundation, Inc.
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.
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.
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/>. */
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>
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>
34 #include <libpspp/misc.h>
36 #include <gsl/gsl_cdf.h>
37 #include <output/table.h>
39 #include <output/charts/plot-chart.h>
40 #include <output/charts/cartesian.h>
43 #define _(msgid) gettext (msgid)
44 #define N_(msgid) msgid
49 const struct variable **vars;
50 const struct dictionary *dict;
52 const struct variable *state_var ;
53 union value state_value;
55 /* Plot the roc curve */
57 /* Plot the reference line */
64 bool bi_neg_exp; /* True iff the bi-negative exponential critieria
66 enum mv_class exclude;
68 bool invert ; /* True iff a smaller test result variable indicates
77 static int run_roc (struct dataset *ds, struct cmd_roc *roc);
80 cmd_roc (struct lexer *lexer, struct dataset *ds)
83 const struct dictionary *dict = dataset_dict (ds);
88 roc.print_coords = false;
91 roc.reference = false;
93 roc.bi_neg_exp = false;
95 roc.pos = roc.pos_weighted = 0;
96 roc.neg = roc.neg_weighted = 0;
97 roc.dict = dataset_dict (ds);
99 if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars,
100 PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
103 if ( ! lex_force_match (lexer, T_BY))
108 roc.state_var = parse_variable (lexer, dict);
110 if ( !lex_force_match (lexer, '('))
115 parse_value (lexer, &roc.state_value, var_get_width (roc.state_var));
118 if ( !lex_force_match (lexer, ')'))
124 while (lex_token (lexer) != '.')
126 lex_match (lexer, '/');
127 if (lex_match_id (lexer, "MISSING"))
129 lex_match (lexer, '=');
130 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
132 if (lex_match_id (lexer, "INCLUDE"))
134 roc.exclude = MV_SYSTEM;
136 else if (lex_match_id (lexer, "EXCLUDE"))
138 roc.exclude = MV_ANY;
142 lex_error (lexer, NULL);
147 else if (lex_match_id (lexer, "PLOT"))
149 lex_match (lexer, '=');
150 if (lex_match_id (lexer, "CURVE"))
153 if (lex_match (lexer, '('))
155 roc.reference = true;
156 lex_force_match_id (lexer, "REFERENCE");
157 lex_force_match (lexer, ')');
160 else if (lex_match_id (lexer, "NONE"))
166 lex_error (lexer, NULL);
170 else if (lex_match_id (lexer, "PRINT"))
172 lex_match (lexer, '=');
173 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
175 if (lex_match_id (lexer, "SE"))
179 else if (lex_match_id (lexer, "COORDINATES"))
181 roc.print_coords = true;
185 lex_error (lexer, NULL);
190 else if (lex_match_id (lexer, "CRITERIA"))
192 lex_match (lexer, '=');
193 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
195 if (lex_match_id (lexer, "CUTOFF"))
197 lex_force_match (lexer, '(');
198 if (lex_match_id (lexer, "INCLUDE"))
200 roc.exclude = MV_SYSTEM;
202 else if (lex_match_id (lexer, "EXCLUDE"))
204 roc.exclude = MV_USER | MV_SYSTEM;
208 lex_error (lexer, NULL);
211 lex_force_match (lexer, ')');
213 else if (lex_match_id (lexer, "TESTPOS"))
215 lex_force_match (lexer, '(');
216 if (lex_match_id (lexer, "LARGE"))
220 else if (lex_match_id (lexer, "SMALL"))
226 lex_error (lexer, NULL);
229 lex_force_match (lexer, ')');
231 else if (lex_match_id (lexer, "CI"))
233 lex_force_match (lexer, '(');
234 lex_force_num (lexer);
235 roc.ci = lex_number (lexer);
237 lex_force_match (lexer, ')');
239 else if (lex_match_id (lexer, "DISTRIBUTION"))
241 lex_force_match (lexer, '(');
242 if (lex_match_id (lexer, "FREE"))
244 roc.bi_neg_exp = false;
246 else if (lex_match_id (lexer, "NEGEXPO"))
248 roc.bi_neg_exp = true;
252 lex_error (lexer, NULL);
255 lex_force_match (lexer, ')');
259 lex_error (lexer, NULL);
266 lex_error (lexer, NULL);
271 if ( ! run_roc (ds, &roc))
285 do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict);
289 run_roc (struct dataset *ds, struct cmd_roc *roc)
291 struct dictionary *dict = dataset_dict (ds);
293 struct casereader *group;
295 struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
296 while (casegrouper_get_next_group (grouper, &group))
298 do_roc (roc, group, dataset_dict (ds));
300 ok = casegrouper_destroy (grouper);
301 ok = proc_commit (ds) && ok;
308 dump_casereader (struct casereader *reader)
311 struct casereader *r = casereader_clone (reader);
313 for ( ; (c = casereader_read (r) ); case_unref (c))
316 for (i = 0 ; i < case_get_value_cnt (c); ++i)
318 printf ("%g ", case_data_idx (c, i)->f);
323 casereader_destroy (r);
328 match_positives (const struct ccase *c, void *aux)
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;
334 bool positive = ( 0 == value_compare_3way (case_data (c, roc->state_var),
336 var_get_width (roc->state_var)));
341 roc->pos_weighted += weight;
346 roc->neg_weighted += weight;
367 struct casewriter *cutpoint_wtr;
368 struct casereader *cutpoint_rdr;
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)
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);
394 double prev_cp = SYSMIS;
396 for ( ; (cpc = casereader_read (r) ); case_unref (cpc))
398 struct ccase *new_case;
399 const double cp = case_data_idx (cpc, CUTPOINT)->f;
401 assert (cp != SYSMIS);
403 /* We don't want duplicates here */
407 new_case = case_clone (cpc);
409 if ( pos_cond (result, cp))
411 case_data_rw_idx (new_case, true_index)->f += weight;
415 case_data_rw_idx (new_case, false_index)->f += weight;
420 casewriter_write (w, new_case);
422 casereader_destroy (r);
424 return casewriter_make_reader (w);
429 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
432 static struct casereader *
433 process_group (const struct variable *var, struct casereader *reader,
434 bool (*pred) (double, double),
435 const struct dictionary *dict,
437 struct casereader **cutpoint_rdr,
438 bool (*pos_cond) (double, double),
442 const struct variable *w = dict_get_weight (dict);
444 struct casereader *r1 =
445 casereader_create_distinct (sort_execute_1var (reader, var), var, w);
447 const int weight_idx = w ? var_get_case_index (w) :
448 caseproto_get_n_widths (casereader_get_proto (r1)) - 1;
452 struct casereader *rclone = casereader_clone (r1);
453 struct casewriter *wtr;
454 struct caseproto *proto = caseproto_create ();
456 proto = caseproto_add_width (proto, 0);
457 proto = caseproto_add_width (proto, 0);
458 proto = caseproto_add_width (proto, 0);
460 wtr = autopaging_writer_create (proto);
464 for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
466 struct ccase *new_case = case_create (proto);
468 struct casereader *r2 = casereader_clone (rclone);
470 const double weight1 = case_data_idx (c1, weight_idx)->f;
471 const double d1 = case_data (c1, var)->f;
475 *cutpoint_rdr = accumulate_counts (*cutpoint_rdr, d1, weight1,
477 true_index, false_index);
481 for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
483 const double d2 = case_data (c2, var)->f;
484 const double weight2 = case_data_idx (c2, weight_idx)->f;
491 else if ( pred (d2, d1))
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;
501 casewriter_write (wtr, new_case);
503 casereader_destroy (r2);
506 casereader_destroy (r1);
507 casereader_destroy (rclone);
509 return casewriter_make_reader (wtr);
513 gt (double d1, double d2)
520 ge (double d1, double d2)
526 lt (double d1, double d2)
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)
536 return process_group (var, reader, gt, dict, &rs->n1,
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)
548 return process_group (var, reader, lt, dict, &rs->n2,
558 append_cutpoint (struct casewriter *writer, double cutpoint)
560 struct ccase *cc = case_create (casewriter_get_proto (writer));
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;
569 casewriter_write (writer, cc);
573 /* Prepare the cutpoints */
575 prepare_cutpoints (struct cmd_roc *roc, struct roc_state *rs, struct casereader *input)
578 struct casereader *r = casereader_clone (input);
580 struct caseproto *proto = caseproto_create ();
582 struct subcase ordering;
583 subcase_init (&ordering, CUTPOINT, 0, SC_ASCEND);
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 */
591 for (i = 0 ; i < roc->n_vars; ++i)
593 rs[i].cutpoint_wtr = sort_create_writer (&ordering, proto);
594 rs[i].prev_result = SYSMIS;
595 rs[i].max = -DBL_MAX;
599 for (; (c = casereader_read (r)) != NULL; case_unref (c))
601 for (i = 0 ; i < roc->n_vars; ++i)
603 const union value *v = case_data (c, roc->vars[i]);
604 const double result = v->f;
606 if ( mv_is_value_missing (var_get_missing_values (roc->vars[i]), v, roc->exclude))
609 minimize (&rs[i].min, result);
610 maximize (&rs[i].max, result);
612 if ( rs[i].prev_result != SYSMIS && rs[i].prev_result != result )
614 const double mean = (result + rs[i].prev_result ) / 2.0;
615 append_cutpoint (rs[i].cutpoint_wtr, mean);
618 rs[i].prev_result = result;
621 casereader_destroy (r);
624 /* Append the min and max cutpoints */
625 for (i = 0 ; i < roc->n_vars; ++i)
627 append_cutpoint (rs[i].cutpoint_wtr, rs[i].min - 1);
628 append_cutpoint (rs[i].cutpoint_wtr, rs[i].max + 1);
630 rs[i].cutpoint_rdr = casewriter_make_reader (rs[i].cutpoint_wtr);
635 do_roc (struct cmd_roc *roc, struct casereader *reader, struct dictionary *dict)
639 struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
641 struct casereader *negatives = NULL;
642 struct casereader *positives = NULL;
644 struct caseproto *n_proto = caseproto_create ();
646 struct subcase up_ordering;
647 struct subcase down_ordering;
649 struct casewriter *neg_wtr = NULL;
651 struct casereader *input = casereader_create_filter_missing (reader,
652 roc->vars, roc->n_vars,
657 input = casereader_create_filter_missing (input,
663 neg_wtr = autopaging_writer_create (casereader_get_proto (input));
665 prepare_cutpoints (roc, rs, input);
668 casereader_create_filter_func (input,
674 n_proto = caseproto_create ();
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);
682 subcase_init (&up_ordering, VALUE, 0, SC_ASCEND);
683 subcase_init (&down_ordering, VALUE, 0, SC_DESCEND);
685 for (i = 0 ; i < roc->n_vars; ++i)
687 struct casewriter *w = NULL;
688 struct casereader *r = NULL;
693 struct casereader *n_neg ;
694 const struct variable *var = roc->vars[i];
696 struct casereader *neg ;
697 struct casereader *pos = casereader_clone (positives);
699 struct casereader *n_pos =
700 process_positive_group (var, pos, dict, &rs[i]);
702 if ( negatives == NULL)
704 negatives = casewriter_make_reader (neg_wtr);
707 neg = casereader_clone (negatives);
709 n_neg = process_negative_group (var, neg, dict, &rs[i]);
711 w = sort_create_writer (&up_ordering, n_proto);
712 for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos))
714 struct ccase *pos_case = case_create (n_proto);
716 const double jpos = case_data_idx (cpos, VALUE)->f;
718 while ((cneg = casereader_read (n_neg)))
720 struct ccase *nc = case_create (n_proto);
722 const double jneg = case_data_idx (cneg, VALUE)->f;
724 case_data_rw_idx (nc, VALUE)->f = jneg;
725 case_data_rw_idx (nc, N_EQ)->f = 0;
727 case_data_rw_idx (nc, N_PRED)->f = SYSMIS;
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);
732 casewriter_write (w, nc);
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;
745 casewriter_write (w, pos_case);
748 r = casewriter_make_reader (w);
751 double prev_pos_gt = rs[i].n1;
752 w = sort_create_writer (&down_ordering, n_proto);
754 for ( ; (c = casereader_read (r) ); case_unref (c))
756 double n_pos_gt = case_data_idx (c, N_PRED)->f;
757 struct ccase *nc = case_clone (c);
759 if ( n_pos_gt == SYSMIS)
761 n_pos_gt = prev_pos_gt;
762 case_data_rw_idx (nc, N_PRED)->f = n_pos_gt;
765 casewriter_write (w, nc);
766 prev_pos_gt = n_pos_gt;
769 r = casewriter_make_reader (w);
773 double prev_neg_lt = rs[i].n2;
774 w = sort_create_writer (&up_ordering, n_proto);
776 for ( ; (c = casereader_read (r) ); case_unref (c))
778 double n_neg_lt = case_data_idx (c, 4)->f;
779 struct ccase *nc = case_clone (c);
781 if ( n_neg_lt == SYSMIS)
783 n_neg_lt = prev_neg_lt;
784 case_data_rw_idx (nc, 4)->f = n_neg_lt;
787 casewriter_write (w, nc);
788 prev_neg_lt = n_neg_lt;
791 r = casewriter_make_reader (w);
795 struct ccase *prev_case = NULL;
796 for ( ; (c = casereader_read (r) ); case_unref (c))
798 const struct ccase *next_case = casereader_peek (r, 0);
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;
806 if ( prev_case && j == case_data_idx (prev_case, VALUE)->f)
808 if ( 0 == case_data_idx (c, N_EQ)->f)
810 n_pos_eq = case_data_idx (prev_case, N_EQ)->f;
811 n_pos_gt = case_data_idx (prev_case, N_PRED)->f;
814 if ( 0 == case_data_idx (c, 3)->f)
816 n_neg_eq = case_data_idx (prev_case, 3)->f;
817 n_neg_lt = case_data_idx (prev_case, 4)->f;
821 if ( NULL == next_case || j != case_data_idx (next_case, VALUE)->f)
823 rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
826 n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
828 n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
832 case_unref (prev_case);
833 prev_case = case_clone (c);
836 rs[i].auc /= rs[i].n1 * rs[i].n2;
838 rs[i].auc = 1 - rs[i].auc;
840 if ( roc->bi_neg_exp )
842 rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
843 rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
847 rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
848 rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
853 casereader_destroy (positives);
854 casereader_destroy (negatives);
856 output_roc (rs, roc);
862 show_auc (struct roc_state *rs, const struct cmd_roc *roc)
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);
870 if ( roc->n_vars > 1)
871 tab_title (tbl, _("Area Under the Curve"));
873 tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
875 tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
877 tab_dim (tbl, tab_natural_dimensions, NULL);
879 tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
881 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
892 tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
893 tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
895 tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
896 tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
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);
905 if ( roc->n_vars > 1)
906 tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
908 if ( roc->n_vars > 1)
909 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
912 for ( i = 0 ; i < roc->n_vars ; ++i )
914 tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
916 tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
921 const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
922 (12 * rs[i].n1 * rs[i].n2));
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));
929 se /= rs[i].n1 * rs[i].n2;
933 tab_double (tbl, n_cols - 4, 2 + i, 0,
937 ci = 1 - roc->ci / 100.0;
938 yy = gsl_cdf_gaussian_Qinv (ci, se) ;
940 tab_double (tbl, n_cols - 2, 2 + i, 0,
944 tab_double (tbl, n_cols - 1, 2 + i, 0,
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)),
959 show_summary (const struct cmd_roc *roc)
961 const int n_cols = 3;
962 const int n_rows = 4;
963 struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
965 tab_title (tbl, _("Case Summary"));
967 tab_headers (tbl, 1, 0, 2, 0);
969 tab_dim (tbl, tab_natural_dimensions, NULL);
978 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
979 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
982 tab_hline (tbl, TAL_2, 1, n_cols - 1, 1);
983 tab_vline (tbl, TAL_1, 2, 1, n_rows - 1);
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"));
990 tab_joint_text (tbl, 1, 0, 2, 0,
991 TAT_TITLE | TAB_CENTER,
992 _("Valid N (listwise)"));
995 tab_text (tbl, 0, 2, TAB_LEFT, _("Positive"));
996 tab_text (tbl, 0, 3, TAB_LEFT, _("Negative"));
999 tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0);
1000 tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0);
1002 tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0);
1003 tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0);
1010 show_coords (struct roc_state *rs, const struct cmd_roc *roc)
1014 const int n_cols = roc->n_vars > 1 ? 4 : 3;
1016 struct tab_table *tbl ;
1018 for (i = 0; i < roc->n_vars; ++i)
1019 n_rows += casereader_count_cases (rs[i].cutpoint_rdr);
1021 tbl = tab_create (n_cols, n_rows, 0);
1023 if ( roc->n_vars > 1)
1024 tab_title (tbl, _("Coordinates of the Curve"));
1026 tab_title (tbl, _("Coordinates of the Curve (%s)"), var_to_string (roc->vars[0]));
1029 tab_headers (tbl, 1, 0, 1, 0);
1031 tab_dim (tbl, tab_natural_dimensions, NULL);
1033 tab_hline (tbl, TAL_2, 0, n_cols - 1, 1);
1035 if ( roc->n_vars > 1)
1036 tab_text (tbl, 0, 0, TAT_TITLE, _("Test variable"));
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"));
1049 if ( roc->n_vars > 1)
1050 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
1052 for (i = 0; i < roc->n_vars; ++i)
1055 struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
1057 if ( roc->n_vars > 1)
1058 tab_text (tbl, 0, x, TAT_TITLE, var_to_string (roc->vars[i]));
1061 tab_hline (tbl, TAL_1, 0, n_cols - 1, x);
1064 for (; (cc = casereader_read (r)) != NULL;
1065 case_unref (cc), x++)
1067 const double se = case_data_idx (cc, TP)->f /
1069 case_data_idx (cc, TP)->f
1071 case_data_idx (cc, FN)->f
1074 const double sp = case_data_idx (cc, TN)->f /
1076 case_data_idx (cc, TN)->f
1078 case_data_idx (cc, FP)->f
1081 tab_double (tbl, n_cols - 3, x, 0, case_data_idx (cc, CUTPOINT)->f,
1082 var_get_print_format (roc->vars[i]));
1084 tab_double (tbl, n_cols - 2, x, 0, se, NULL);
1085 tab_double (tbl, n_cols - 1, x, 0, 1 - sp, NULL);
1088 casereader_destroy (r);
1096 draw_roc (struct roc_state *rs, const struct cmd_roc *roc)
1100 struct chart *roc_chart = chart_create ();
1102 chart_write_title (roc_chart, _("ROC Curve"));
1103 chart_write_xlabel (roc_chart, _("1 - Specificity"));
1104 chart_write_ylabel (roc_chart, _("Sensitivity"));
1106 chart_write_xscale (roc_chart, 0, 1, 5);
1107 chart_write_yscale (roc_chart, 0, 1, 5);
1109 if ( roc->reference )
1111 chart_line (roc_chart, 1.0, 0,
1116 for (i = 0; i < roc->n_vars; ++i)
1119 struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
1121 chart_vector_start (roc_chart, var_get_name (roc->vars[i]));
1122 for (; (cc = casereader_read (r)) != NULL;
1125 double se = case_data_idx (cc, TP)->f;
1126 double sp = case_data_idx (cc, TN)->f;
1128 se /= case_data_idx (cc, FN)->f +
1129 case_data_idx (cc, TP)->f ;
1131 sp /= case_data_idx (cc, TN)->f +
1132 case_data_idx (cc, FP)->f ;
1134 chart_vector (roc_chart, 1 - sp, se);
1136 chart_vector_end (roc_chart);
1137 casereader_destroy (r);
1140 chart_write_legend (roc_chart);
1142 chart_submit (roc_chart);
1147 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
1157 if ( roc->print_coords )
1158 show_coords (rs, roc);