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>
32 #include <libpspp/misc.h>
34 #include <gsl/gsl_cdf.h>
35 #include <output/table.h>
38 #define _(msgid) gettext (msgid)
39 #define N_(msgid) msgid
44 const struct variable **vars;
46 struct variable *state_var ;
47 union value state_value;
49 /* Plot the roc curve */
51 /* Plot the reference line */
58 bool bi_neg_exp; /* True iff the bi-negative exponential critieria
60 enum mv_class exclude;
62 bool invert ; /* True iff a smaller test result variable indicates
66 static int run_roc (struct dataset *ds, struct cmd_roc *roc);
69 cmd_roc (struct lexer *lexer, struct dataset *ds)
72 const struct dictionary *dict = dataset_dict (ds);
77 roc.print_coords = false;
80 roc.reference = false;
82 roc.bi_neg_exp = false;
85 if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars,
86 PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
89 if ( ! lex_force_match (lexer, T_BY))
94 roc.state_var = parse_variable (lexer, dict);
96 if ( !lex_force_match (lexer, '('))
101 parse_value (lexer, &roc.state_value, var_get_width (roc.state_var));
104 if ( !lex_force_match (lexer, ')'))
110 while (lex_token (lexer) != '.')
112 lex_match (lexer, '/');
113 if (lex_match_id (lexer, "MISSING"))
115 lex_match (lexer, '=');
116 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
118 if (lex_match_id (lexer, "INCLUDE"))
120 roc.exclude = MV_SYSTEM;
122 else if (lex_match_id (lexer, "EXCLUDE"))
124 roc.exclude = MV_ANY;
128 lex_error (lexer, NULL);
133 else if (lex_match_id (lexer, "PLOT"))
135 lex_match (lexer, '=');
136 if (lex_match_id (lexer, "CURVE"))
139 if (lex_match (lexer, '('))
141 roc.reference = true;
142 lex_force_match_id (lexer, "REFERENCE");
143 lex_force_match (lexer, ')');
146 else if (lex_match_id (lexer, "NONE"))
152 lex_error (lexer, NULL);
156 else if (lex_match_id (lexer, "PRINT"))
158 lex_match (lexer, '=');
159 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
161 if (lex_match_id (lexer, "SE"))
165 else if (lex_match_id (lexer, "COORDINATES"))
167 roc.print_coords = true;
171 lex_error (lexer, NULL);
176 else if (lex_match_id (lexer, "CRITERIA"))
178 lex_match (lexer, '=');
179 while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
181 if (lex_match_id (lexer, "CUTOFF"))
183 lex_force_match (lexer, '(');
184 if (lex_match_id (lexer, "INCLUDE"))
186 roc.exclude = MV_SYSTEM;
188 else if (lex_match_id (lexer, "EXCLUDE"))
190 roc.exclude = MV_USER | MV_SYSTEM;
194 lex_error (lexer, NULL);
197 lex_force_match (lexer, ')');
199 else if (lex_match_id (lexer, "TESTPOS"))
201 lex_force_match (lexer, '(');
202 if (lex_match_id (lexer, "LARGE"))
206 else if (lex_match_id (lexer, "SMALL"))
212 lex_error (lexer, NULL);
215 lex_force_match (lexer, ')');
217 else if (lex_match_id (lexer, "CI"))
219 lex_force_match (lexer, '(');
220 lex_force_num (lexer);
221 roc.ci = lex_number (lexer);
223 lex_force_match (lexer, ')');
225 else if (lex_match_id (lexer, "DISTRIBUTION"))
227 lex_force_match (lexer, '(');
228 if (lex_match_id (lexer, "FREE"))
230 roc.bi_neg_exp = false;
232 else if (lex_match_id (lexer, "NEGEXPO"))
234 roc.bi_neg_exp = true;
238 lex_error (lexer, NULL);
241 lex_force_match (lexer, ')');
245 lex_error (lexer, NULL);
252 lex_error (lexer, NULL);
266 do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict);
270 run_roc (struct dataset *ds, struct cmd_roc *roc)
272 struct dictionary *dict = dataset_dict (ds);
274 struct casereader *group;
276 struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
277 while (casegrouper_get_next_group (grouper, &group))
279 do_roc (roc, group, dataset_dict (ds));
281 ok = casegrouper_destroy (grouper);
282 ok = proc_commit (ds) && ok;
289 dump_casereader (struct casereader *reader)
292 struct casereader *r = casereader_clone (reader);
294 for ( ; (c = casereader_read (r) ); case_unref (c))
297 for (i = 0 ; i < case_get_value_cnt (c); ++i)
299 printf ("%g ", case_data_idx (c, i)->f);
304 casereader_destroy (r);
308 match_positives (const struct ccase *c, void *aux)
310 struct cmd_roc *roc = aux;
312 return 0 == value_compare_3way (case_data (c, roc->state_var),
314 var_get_width (roc->state_var));
334 static void output_roc (struct roc_state *rs, const struct cmd_roc *roc);
338 static struct casereader *
339 process_group (const struct variable *var, struct casereader *reader,
340 bool (*pred) (double, double),
341 const struct dictionary *dict,
344 const struct variable *w = dict_get_weight (dict);
345 const int weight_idx = w ? var_get_case_index (w) :
346 caseproto_get_n_widths (casereader_get_proto (reader)) - 1;
348 struct casereader *r1 =
349 casereader_create_distinct (sort_execute_1var (reader, var), var, w);
353 struct casereader *rclone = casereader_clone (r1);
354 struct casewriter *wtr;
355 struct caseproto *proto = caseproto_create ();
357 proto = caseproto_add_width (proto, 0);
358 proto = caseproto_add_width (proto, 0);
359 proto = caseproto_add_width (proto, 0);
361 wtr = autopaging_writer_create (proto);
365 for ( ; (c1 = casereader_read (r1) ); case_unref (c1))
368 struct casereader *r2 = casereader_clone (rclone);
370 const double weight1 = case_data_idx (c1, weight_idx)->f;
371 const double d1 = case_data (c1, var)->f;
376 struct ccase *new_case = case_create (proto);
380 for ( ; (c2 = casereader_read (r2) ); case_unref (c2))
382 const double d2 = case_data (c2, var)->f;
383 const double weight2 = case_data_idx (c2, weight_idx)->f;
390 else if ( pred (d2, d1))
396 case_data_rw_idx (new_case, VALUE)->f = d1;
397 case_data_rw_idx (new_case, N_EQ)->f = n_eq;
398 case_data_rw_idx (new_case, N_PRED)->f = n_pred;
400 casewriter_write (wtr, new_case);
402 casereader_destroy (r2);
405 casereader_destroy (r1);
406 casereader_destroy (rclone);
408 return casewriter_make_reader (wtr);
412 gt (double d1, double d2)
418 lt (double d1, double d2)
425 do_roc (struct cmd_roc *roc, struct casereader *input, struct dictionary *dict)
429 struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
431 const struct caseproto *proto = casereader_get_proto (input);
433 struct casewriter *neg_wtr = autopaging_writer_create (proto);
435 struct casereader *negatives = NULL;
437 struct casereader *positives =
438 casereader_create_filter_func (input,
445 for (i = 0 ; i < roc->n_vars; ++i)
448 struct casereader *n_neg ;
449 const struct variable *var = roc->vars[i];
451 struct casereader *neg ;
452 struct casereader *pos = casereader_clone (positives);
454 struct casereader *n_pos = process_group (var, pos, gt, dict, &rs[i].n1);
456 if ( negatives == NULL)
458 negatives = casewriter_make_reader (neg_wtr);
461 neg = casereader_clone (negatives);
463 n_neg = process_group (var, neg, lt, dict, &rs[i].n2);
465 /* Simple join on VALUE */
466 for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos))
468 struct ccase *cneg = NULL;
469 double dneg = -DBL_MAX;
470 const double dpos = case_data_idx (cpos, VALUE)->f;
476 cneg = casereader_read (n_neg);
479 dneg = case_data_idx (cneg, VALUE)->f;
484 double n_pos_eq = case_data_idx (cpos, N_EQ)->f;
485 double n_neg_eq = case_data_idx (cneg, N_EQ)->f;
486 double n_pos_gt = case_data_idx (cpos, N_PRED)->f;
487 double n_neg_lt = case_data_idx (cneg, N_PRED)->f;
489 rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
491 n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
493 n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
500 rs[i].auc /= rs[i].n1 * rs[i].n2;
502 rs[i].auc = 1 - rs[i].auc;
504 if ( roc->bi_neg_exp )
506 rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
507 rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
511 rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
512 rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
516 casereader_destroy (positives);
517 casereader_destroy (negatives);
519 output_roc (rs, roc);
528 show_auc (struct roc_state *rs, const struct cmd_roc *roc)
531 const int n_fields = roc->print_se ? 5 : 1;
532 const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields;
533 const int n_rows = 2 + roc->n_vars;
534 struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
536 if ( roc->n_vars > 1)
537 tab_title (tbl, _("Area Under the Curve"));
539 tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
541 tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
543 tab_dim (tbl, tab_natural_dimensions, NULL);
545 tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
547 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
558 tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
559 tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
561 tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
562 tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
564 tab_joint_text (tbl, n_cols - 2, 0, 4, 0,
565 TAT_TITLE | TAB_CENTER | TAT_PRINTF,
566 _("Asymp. %g%% Confidence Interval"), roc->ci);
567 tab_vline (tbl, 0, n_cols - 1, 0, 0);
568 tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1);
571 if ( roc->n_vars > 1)
572 tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
574 if ( roc->n_vars > 1)
575 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
578 for ( i = 0 ; i < roc->n_vars ; ++i )
580 tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
582 tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
588 const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
589 (12 * rs[i].n1 * rs[i].n2));
593 se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (rs[i].q1hat - pow2 (rs[i].auc)) +
594 (rs[i].n2 - 1) * (rs[i].q2hat - pow2 (rs[i].auc));
596 se /= rs[i].n1 * rs[i].n2;
600 tab_double (tbl, n_cols - 4, 2 + i, 0,
604 ci = 1 - roc->ci / 100.0;
605 yy = gsl_cdf_gaussian_Qinv (ci, se) ;
607 tab_double (tbl, n_cols - 2, 2 + i, 0,
611 tab_double (tbl, n_cols - 1, 2 + i, 0,
615 tab_double (tbl, n_cols - 3, 2 + i, 0,
616 2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)),
626 show_summary (const struct cmd_roc *roc)
628 const int n_cols = 3;
629 const int n_rows = 4;
630 struct tab_table *tbl = tab_create (n_cols, n_rows, 0);
632 tab_title (tbl, _("Case Summary"));
634 tab_headers (tbl, 1, 0, 2, 0);
636 tab_dim (tbl, tab_natural_dimensions, NULL);
645 tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
646 tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
649 tab_hline (tbl, TAL_2, 1, n_cols - 1, 1);
650 tab_vline (tbl, TAL_1, 2, 1, n_rows - 1);
653 tab_text (tbl, 0, 1, TAT_TITLE | TAB_LEFT, var_to_string (roc->state_var));
654 tab_text (tbl, 1, 1, TAT_TITLE, _("Unweighted"));
655 tab_text (tbl, 2, 1, TAT_TITLE, _("Weighted"));
657 tab_joint_text (tbl, 1, 0, 2, 0,
658 TAT_TITLE | TAB_CENTER,
659 _("Valid N (listwise)"));
662 tab_text (tbl, 0, 2, TAB_LEFT, _("Positive"));
663 tab_text (tbl, 0, 3, TAB_LEFT, _("Negative"));
667 tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0);
668 tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0);
670 tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0);
671 tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0);
679 output_roc (struct roc_state *rs, const struct cmd_roc *roc)
692 if ( roc->print_coords )
693 show_coords (rs, roc);