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/>. */
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>
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>
40 #define _(msgid) gettext (msgid)
41 #define N_(msgid) msgid
46 const struct variable **vars;
48 struct variable *state_var ;
49 union value state_value;
51 /* Plot the roc curve */
53 /* Plot the reference line */
60 bool bi_neg_exp; /* True iff the bi-negative exponential critieria
62 enum mv_class exclude;
64 bool invert ; /* True iff a smaller test result variable indicates
68 static int run_roc (struct dataset *ds, struct cmd_roc *roc);
71 cmd_roc (struct lexer *lexer, struct dataset *ds)
74 const struct dictionary *dict = dataset_dict (ds);
79 roc.print_coords = false;
82 roc.reference = false;
84 roc.bi_neg_exp = false;
87 if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars,
88 PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
91 if ( ! lex_force_match (lexer, T_BY))
96 roc.state_var = parse_variable (lexer, dict);
98 if ( !lex_force_match (lexer, '('))
103 parse_value (lexer, &roc.state_value, var_get_width (roc.state_var));
106 if ( !lex_force_match (lexer, ')'))
112 while (lex_token (lexer) != '.')
114 lex_match (lexer, '/');
115 if (lex_match_id (lexer, "MISSING"))
117 lex_match (lexer, '=');
118 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
120 if (lex_match_id (lexer, "INCLUDE"))
122 roc.exclude = MV_SYSTEM;
124 else if (lex_match_id (lexer, "EXCLUDE"))
126 roc.exclude = MV_ANY;
130 lex_error (lexer, NULL);
135 else if (lex_match_id (lexer, "PLOT"))
137 lex_match (lexer, '=');
138 if (lex_match_id (lexer, "CURVE"))
141 if (lex_match (lexer, '('))
143 roc.reference = true;
144 lex_force_match_id (lexer, "REFERENCE");
145 lex_force_match (lexer, ')');
148 else if (lex_match_id (lexer, "NONE"))
154 lex_error (lexer, NULL);
158 else if (lex_match_id (lexer, "PRINT"))
160 lex_match (lexer, '=');
161 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
163 if (lex_match_id (lexer, "SE"))
167 else if (lex_match_id (lexer, "COORDINATES"))
169 roc.print_coords = true;
173 lex_error (lexer, NULL);
178 else if (lex_match_id (lexer, "CRITERIA"))
180 lex_match (lexer, '=');
181 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
183 if (lex_match_id (lexer, "CUTOFF"))
185 lex_force_match (lexer, '(');
186 if (lex_match_id (lexer, "INCLUDE"))
188 roc.exclude = MV_SYSTEM;
190 else if (lex_match_id (lexer, "EXCLUDE"))
192 roc.exclude = MV_USER | MV_SYSTEM;
196 lex_error (lexer, NULL);
199 lex_force_match (lexer, ')');
201 else if (lex_match_id (lexer, "TESTPOS"))
203 lex_force_match (lexer, '(');
204 if (lex_match_id (lexer, "LARGE"))
208 else if (lex_match_id (lexer, "SMALL"))
214 lex_error (lexer, NULL);
217 lex_force_match (lexer, ')');
219 else if (lex_match_id (lexer, "CI"))
221 lex_force_match (lexer, '(');
222 lex_force_num (lexer);
223 roc.ci = lex_number (lexer);
225 lex_force_match (lexer, ')');
227 else if (lex_match_id (lexer, "DISTRIBUTION"))
229 lex_force_match (lexer, '(');
230 if (lex_match_id (lexer, "FREE"))
232 roc.bi_neg_exp = false;
234 else if (lex_match_id (lexer, "NEGEXPO"))
236 roc.bi_neg_exp = true;
240 lex_error (lexer, NULL);
243 lex_force_match (lexer, ')');
247 lex_error (lexer, NULL);
254 lex_error (lexer, NULL);
268 do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict);
272 run_roc (struct dataset *ds, struct cmd_roc *roc)
274 struct dictionary *dict = dataset_dict (ds);
276 struct casereader *group;
278 struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
279 while (casegrouper_get_next_group (grouper, &group))
281 do_roc (roc, group, dataset_dict (ds));
283 ok = casegrouper_destroy (grouper);
284 ok = proc_commit (ds) && ok;
291 dump_casereader (struct casereader *reader)
294 struct casereader *r = casereader_clone (reader);
296 for ( ; (c = casereader_read (r) ); case_unref (c))
299 for (i = 0 ; i < case_get_value_cnt (c); ++i)
301 printf ("%g ", case_data_idx (c, i)->f);
306 casereader_destroy (r);
310 match_positives (const struct ccase *c, void *aux)
312 struct cmd_roc *roc = aux;
314 return 0 == value_compare_3way (case_data (c, roc->state_var),
316 var_get_width (roc->state_var));
334 struct casewriter *cutpoint_wtr;
335 struct casereader *cutpoint_rdr;
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)
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);
361 double prev_cp = SYSMIS;
364 for ( ; (cpc = casereader_read (r) ); case_unref (cpc))
366 struct ccase *new_case;
367 const double cp = case_data_idx (cpc, CUTPOINT)->f;
369 /* We don't want duplicates here */
373 new_case = case_clone (cpc);
375 if ( pos_cond (result, cp))
377 case_data_rw_idx (new_case, true_index)->f += weight;
381 case_data_rw_idx (new_case, false_index)->f += weight;
386 casewriter_write (w, new_case);
388 casereader_destroy (r);
390 return casewriter_make_reader (w);
395 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
398 static struct casereader *
399 process_group (const struct variable *var, struct casereader *reader,
400 bool (*pred) (double, double),
401 const struct dictionary *dict,
403 struct casereader **cutpoint_rdr,
404 bool (*pos_cond) (double, double),
409 const struct variable *w = dict_get_weight (dict);
410 struct casereader *r1 =
411 casereader_create_distinct (sort_execute_1var (reader, var), var, w);
413 const int weight_idx = w ? var_get_case_index (w) :
414 caseproto_get_n_widths (casereader_get_proto (r1)) - 1;
418 struct casereader *rclone = casereader_clone (r1);
419 struct casewriter *wtr;
420 struct caseproto *proto = caseproto_create ();
422 proto = caseproto_add_width (proto, 0);
423 proto = caseproto_add_width (proto, 0);
424 proto = caseproto_add_width (proto, 0);
426 wtr = autopaging_writer_create (proto);
430 for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
433 struct casereader *r2 = casereader_clone (rclone);
435 const double weight1 = case_data_idx (c1, weight_idx)->f;
436 const double d1 = case_data (c1, var)->f;
440 *cutpoint_rdr = accumulate_counts (*cutpoint_rdr, d1, weight1,
442 true_index, false_index);
444 struct ccase *new_case = case_create (proto);
448 for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
450 const double d2 = case_data (c2, var)->f;
451 const double weight2 = case_data_idx (c2, weight_idx)->f;
458 else if ( pred (d2, d1))
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;
468 casewriter_write (wtr, new_case);
470 casereader_destroy (r2);
473 casereader_destroy (r1);
474 casereader_destroy (rclone);
476 return casewriter_make_reader (wtr);
480 gt (double d1, double d2)
487 ge (double d1, double d2)
493 lt (double d1, double d2)
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)
503 return process_group (var, reader, gt, dict, &rs->n1,
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)
515 return process_group (var, reader, lt, dict, &rs->n2,
525 append_cutpoint (struct casewriter *writer, double cutpoint)
527 struct ccase *cc = case_create (casewriter_get_proto (writer));
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;
536 casewriter_write (writer, cc);
541 do_roc (struct cmd_roc *roc, struct casereader *input, struct dictionary *dict)
545 struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
547 struct casewriter *neg_wtr = autopaging_writer_create (casereader_get_proto (input));
549 struct casereader *negatives = NULL;
550 struct casereader *positives = NULL;
553 /* Prepare the cutpoints */
555 struct casereader *r = casereader_clone (input);
557 struct caseproto *proto = caseproto_create ();
559 struct subcase ordering;
560 struct variable *iv = var_create_internal (CUTPOINT);
561 subcase_init_var (&ordering, iv, SC_ASCEND);
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 */
571 for (i = 0 ; i < roc->n_vars; ++i)
573 rs[i].cutpoint_wtr = sort_create_writer (&ordering, proto);
574 rs[i].prev_result = SYSMIS;
575 rs[i].max = -DBL_MAX;
579 for (; (c = casereader_read (r)) != NULL; case_unref (c))
581 const double weight = dict_get_case_weight (dict, c, NULL);
582 for (i = 0 ; i < roc->n_vars; ++i)
584 const double result = case_data (c, roc->vars[i])->f;
586 minimize (&rs[i].min, result);
587 maximize (&rs[i].max, result);
589 if ( rs[i].prev_result != SYSMIS && rs[i].prev_result != result )
591 const double mean = (result + rs[i].prev_result ) / 2.0;
592 append_cutpoint (rs[i].cutpoint_wtr, mean);
595 rs[i].prev_result = result;
598 casereader_destroy (r);
601 /* Append the min and max cutpoints */
602 for (i = 0 ; i < roc->n_vars; ++i)
604 append_cutpoint (rs[i].cutpoint_wtr, rs[i].min - 1);
605 append_cutpoint (rs[i].cutpoint_wtr, rs[i].max + 1);
607 rs[i].cutpoint_rdr = casewriter_make_reader (rs[i].cutpoint_wtr);
612 casereader_create_filter_func (input,
619 for (i = 0 ; i < roc->n_vars; ++i)
622 struct casereader *n_neg ;
623 const struct variable *var = roc->vars[i];
625 struct casereader *neg ;
626 struct casereader *pos = casereader_clone (positives);
628 struct casereader *n_pos =
629 process_positive_group (var, pos, dict, &rs[i]);
631 if ( negatives == NULL)
633 negatives = casewriter_make_reader (neg_wtr);
636 neg = casereader_clone (negatives);
638 n_neg = process_negative_group (var, neg, dict, &rs[i]);
641 /* Simple join on VALUE */
642 for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos))
644 struct ccase *cneg = NULL;
645 double dneg = -DBL_MAX;
646 const double dpos = case_data_idx (cpos, VALUE)->f;
652 cneg = casereader_read (n_neg);
655 dneg = case_data_idx (cneg, VALUE)->f;
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;
665 rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
667 n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
669 n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
676 rs[i].auc /= rs[i].n1 * rs[i].n2;
678 rs[i].auc = 1 - rs[i].auc;
680 if ( roc->bi_neg_exp )
682 rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
683 rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
687 rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
688 rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
692 casereader_destroy (positives);
693 casereader_destroy (negatives);
695 output_roc (rs, roc);
704 show_auc (struct roc_state *rs, const struct cmd_roc *roc)
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);
712 if ( roc->n_vars > 1)
713 tab_title (tbl, _("Area Under the Curve"));
715 tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
717 tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
719 tab_dim (tbl, tab_natural_dimensions, NULL);
721 tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
723 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
734 tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
735 tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
737 tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
738 tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
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);
747 if ( roc->n_vars > 1)
748 tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
750 if ( roc->n_vars > 1)
751 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
754 for ( i = 0 ; i < roc->n_vars ; ++i )
756 tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
758 tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
763 const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
764 (12 * rs[i].n1 * rs[i].n2));
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));
771 se /= rs[i].n1 * rs[i].n2;
775 tab_double (tbl, n_cols - 4, 2 + i, 0,
779 ci = 1 - roc->ci / 100.0;
780 yy = gsl_cdf_gaussian_Qinv (ci, se) ;
782 tab_double (tbl, n_cols - 2, 2 + i, 0,
786 tab_double (tbl, n_cols - 1, 2 + i, 0,
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)),
801 show_summary (const struct cmd_roc *roc)
803 const int n_cols = 3;
804 const int n_rows = 4;
805 struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
807 tab_title (tbl, _("Case Summary"));
809 tab_headers (tbl, 1, 0, 2, 0);
811 tab_dim (tbl, tab_natural_dimensions, NULL);
820 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
821 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
824 tab_hline (tbl, TAL_2, 1, n_cols - 1, 1);
825 tab_vline (tbl, TAL_1, 2, 1, n_rows - 1);
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"));
832 tab_joint_text (tbl, 1, 0, 2, 0,
833 TAT_TITLE | TAB_CENTER,
834 _("Valid N (listwise)"));
837 tab_text (tbl, 0, 2, TAB_LEFT, _("Positive"));
838 tab_text (tbl, 0, 3, TAB_LEFT, _("Negative"));
842 tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0);
843 tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0);
845 tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0);
846 tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0);
854 show_coords (struct roc_state *rs, const struct cmd_roc *roc)
858 const int n_cols = roc->n_vars > 1 ? 4 : 3;
860 struct tab_table *tbl ;
862 for (i = 0; i < roc->n_vars; ++i)
863 n_rows += casereader_count_cases (rs[i].cutpoint_rdr);
865 tbl = tab_create (n_cols, n_rows, 0);
867 if ( roc->n_vars > 1)
868 tab_title (tbl, _("Coordinates of the Curve"));
870 tab_title (tbl, _("Coordinates of the Curve (%s)"), var_to_string (roc->vars[0]));
873 tab_headers (tbl, 1, 0, 1, 0);
875 tab_dim (tbl, tab_natural_dimensions, NULL);
877 tab_hline (tbl, TAL_2, 0, n_cols - 1, 1);
879 if ( roc->n_vars > 1)
880 tab_text (tbl, 0, 0, TAT_TITLE, _("Test variable"));
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"));
893 if ( roc->n_vars > 1)
894 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
896 for (i = 0; i < roc->n_vars; ++i)
899 struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
901 if ( roc->n_vars > 1)
902 tab_text (tbl, 0, x, TAT_TITLE, var_to_string (roc->vars[i]));
905 tab_hline (tbl, TAL_1, 0, n_cols - 1, x);
908 for (; (cc = casereader_read (r)) != NULL;
909 case_unref (cc), x++)
911 const double se = case_data_idx (cc, TP)->f /
913 case_data_idx (cc, TP)->f
915 case_data_idx (cc, FN)->f
918 const double sp = case_data_idx (cc, TN)->f /
920 case_data_idx (cc, TN)->f
922 case_data_idx (cc, FP)->f
925 tab_double (tbl, n_cols - 3, x, 0, case_data_idx (cc, CUTPOINT)->f,
926 var_get_print_format (roc->vars[i]));
928 tab_double (tbl, n_cols - 2, x, 0, se, NULL);
929 tab_double (tbl, n_cols - 1, x, 0, 1 - sp, NULL);
932 casereader_destroy (r);
941 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
954 if ( roc->print_coords )
955 show_coords (rs, roc);