+ int i;
+
+ struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs);
+
+ struct casereader *negatives = NULL;
+ struct casereader *positives = NULL;
+
+ struct caseproto *n_proto = NULL;
+
+ struct subcase up_ordering;
+ struct subcase down_ordering;
+
+ struct casewriter *neg_wtr = NULL;
+
+ struct casereader *input = casereader_create_filter_missing (reader,
+ roc->vars, roc->n_vars,
+ roc->exclude,
+ NULL,
+ NULL);
+
+ input = casereader_create_filter_missing (input,
+ &roc->state_var, 1,
+ roc->exclude,
+ NULL,
+ NULL);
+
+ neg_wtr = autopaging_writer_create (casereader_get_proto (input));
+
+ prepare_cutpoints (roc, rs, input);
+
+
+ /* Separate the positive actual state cases from the negative ones */
+ positives =
+ casereader_create_filter_func (input,
+ match_positives,
+ NULL,
+ roc,
+ neg_wtr);
+
+ n_proto = caseproto_create ();
+
+ n_proto = caseproto_add_width (n_proto, 0);
+ n_proto = caseproto_add_width (n_proto, 0);
+ n_proto = caseproto_add_width (n_proto, 0);
+ n_proto = caseproto_add_width (n_proto, 0);
+ n_proto = caseproto_add_width (n_proto, 0);
+
+ subcase_init (&up_ordering, VALUE, 0, SC_ASCEND);
+ subcase_init (&down_ordering, VALUE, 0, SC_DESCEND);
+
+ for (i = 0 ; i < roc->n_vars; ++i)
+ {
+ struct casewriter *w = NULL;
+ struct casereader *r = NULL;
+
+ struct ccase *c;
+
+ struct ccase *cpos;
+ struct casereader *n_neg_reader ;
+ const struct variable *var = roc->vars[i];
+
+ struct casereader *neg ;
+ struct casereader *pos = casereader_clone (positives);
+
+ struct casereader *n_pos_reader =
+ process_positive_group (var, pos, dict, &rs[i]);
+
+ if ( negatives == NULL)
+ {
+ negatives = casewriter_make_reader (neg_wtr);
+ }
+
+ neg = casereader_clone (negatives);
+
+ n_neg_reader = process_negative_group (var, neg, dict, &rs[i]);
+
+ /* Merge the n_pos and n_neg casereaders */
+ w = sort_create_writer (&up_ordering, n_proto);
+ for ( ; (cpos = casereader_read (n_pos_reader) ); case_unref (cpos))
+ {
+ struct ccase *pos_case = case_create (n_proto);
+ struct ccase *cneg;
+ const double jpos = case_data_idx (cpos, VALUE)->f;
+
+ while ((cneg = casereader_read (n_neg_reader)))
+ {
+ struct ccase *nc = case_create (n_proto);
+
+ const double jneg = case_data_idx (cneg, VALUE)->f;
+
+ case_data_rw_idx (nc, VALUE)->f = jneg;
+ case_data_rw_idx (nc, N_POS_EQ)->f = 0;
+
+ case_data_rw_idx (nc, N_POS_GT)->f = SYSMIS;
+
+ *case_data_rw_idx (nc, N_NEG_EQ) = *case_data_idx (cneg, N_EQ);
+ *case_data_rw_idx (nc, N_NEG_LT) = *case_data_idx (cneg, N_PRED);
+
+ casewriter_write (w, nc);
+
+ case_unref (cneg);
+ if ( jneg > jpos)
+ break;
+ }
+
+ case_data_rw_idx (pos_case, VALUE)->f = jpos;
+ *case_data_rw_idx (pos_case, N_POS_EQ) = *case_data_idx (cpos, N_EQ);
+ *case_data_rw_idx (pos_case, N_POS_GT) = *case_data_idx (cpos, N_PRED);
+ case_data_rw_idx (pos_case, N_NEG_EQ)->f = 0;
+ case_data_rw_idx (pos_case, N_NEG_LT)->f = SYSMIS;
+
+ casewriter_write (w, pos_case);
+ }
+
+ casereader_destroy (n_pos_reader);
+ casereader_destroy (n_neg_reader);
+
+/* These aren't used anymore */
+#undef N_EQ
+#undef N_PRED
+
+ r = casewriter_make_reader (w);
+
+ /* Propagate the N_POS_GT values from the positive cases
+ to the negative ones */
+ {
+ double prev_pos_gt = rs[i].n1;
+ w = sort_create_writer (&down_ordering, n_proto);
+
+ for ( ; (c = casereader_read (r) ); case_unref (c))
+ {
+ double n_pos_gt = case_data_idx (c, N_POS_GT)->f;
+ struct ccase *nc = case_clone (c);
+
+ if ( n_pos_gt == SYSMIS)
+ {
+ n_pos_gt = prev_pos_gt;
+ case_data_rw_idx (nc, N_POS_GT)->f = n_pos_gt;
+ }
+
+ casewriter_write (w, nc);
+ prev_pos_gt = n_pos_gt;
+ }
+
+ casereader_destroy (r);
+ r = casewriter_make_reader (w);
+ }
+
+ /* Propagate the N_NEG_LT values from the negative cases
+ to the positive ones */
+ {
+ double prev_neg_lt = rs[i].n2;
+ w = sort_create_writer (&up_ordering, n_proto);
+
+ for ( ; (c = casereader_read (r) ); case_unref (c))
+ {
+ double n_neg_lt = case_data_idx (c, N_NEG_LT)->f;
+ struct ccase *nc = case_clone (c);
+
+ if ( n_neg_lt == SYSMIS)
+ {
+ n_neg_lt = prev_neg_lt;
+ case_data_rw_idx (nc, N_NEG_LT)->f = n_neg_lt;
+ }
+
+ casewriter_write (w, nc);
+ prev_neg_lt = n_neg_lt;
+ }
+
+ casereader_destroy (r);
+ r = casewriter_make_reader (w);
+ }
+
+ {
+ struct ccase *prev_case = NULL;
+ for ( ; (c = casereader_read (r) ); case_unref (c))
+ {
+ struct ccase *next_case = casereader_peek (r, 0);
+
+ const double j = case_data_idx (c, VALUE)->f;
+ double n_pos_eq = case_data_idx (c, N_POS_EQ)->f;
+ double n_pos_gt = case_data_idx (c, N_POS_GT)->f;
+ double n_neg_eq = case_data_idx (c, N_NEG_EQ)->f;
+ double n_neg_lt = case_data_idx (c, N_NEG_LT)->f;
+
+ if ( prev_case && j == case_data_idx (prev_case, VALUE)->f)
+ {
+ if ( 0 == case_data_idx (c, N_POS_EQ)->f)
+ {
+ n_pos_eq = case_data_idx (prev_case, N_POS_EQ)->f;
+ n_pos_gt = case_data_idx (prev_case, N_POS_GT)->f;
+ }
+
+ if ( 0 == case_data_idx (c, N_NEG_EQ)->f)
+ {
+ n_neg_eq = case_data_idx (prev_case, N_NEG_EQ)->f;
+ n_neg_lt = case_data_idx (prev_case, N_NEG_LT)->f;
+ }
+ }
+
+ if ( NULL == next_case || j != case_data_idx (next_case, VALUE)->f)
+ {
+ rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0;
+
+ rs[i].q1hat +=
+ n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0);
+ rs[i].q2hat +=
+ n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0);
+
+ }
+
+ case_unref (next_case);
+ case_unref (prev_case);
+ prev_case = case_clone (c);
+ }
+ casereader_destroy (r);
+ case_unref (prev_case);
+
+ rs[i].auc /= rs[i].n1 * rs[i].n2;
+ if ( roc->invert )
+ rs[i].auc = 1 - rs[i].auc;
+
+ if ( roc->bi_neg_exp )
+ {
+ rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc);
+ rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc);
+ }
+ else
+ {
+ rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1);
+ rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2);
+ }
+ }
+ }
+
+ casereader_destroy (positives);
+ casereader_destroy (negatives);
+
+ caseproto_unref (n_proto);
+ subcase_destroy (&up_ordering);
+ subcase_destroy (&down_ordering);
+
+ output_roc (rs, roc);
+
+ for (i = 0 ; i < roc->n_vars; ++i)
+ casereader_destroy (rs[i].cutpoint_rdr);
+
+ free (rs);
+}
+
+static void
+show_auc (struct roc_state *rs, const struct cmd_roc *roc)
+{
+ int i;
+ const int n_fields = roc->print_se ? 5 : 1;
+ const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields;
+ const int n_rows = 2 + roc->n_vars;
+ struct tab_table *tbl = tab_create (n_cols, n_rows);
+
+ if ( roc->n_vars > 1)
+ tab_title (tbl, _("Area Under the Curve"));
+ else
+ tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0]));
+
+ tab_headers (tbl, n_cols - n_fields, 0, 1, 0);
+
+
+ tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area"));
+
+ tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
+
+ tab_box (tbl,
+ TAL_2, TAL_2,
+ -1, TAL_1,
+ 0, 0,
+ n_cols - 1,
+ n_rows - 1);
+
+ if ( roc->print_se )
+ {
+ tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error"));
+ tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig."));
+
+ tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound"));
+ tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound"));
+
+ tab_joint_text_format (tbl, n_cols - 2, 0, 4, 0,
+ TAT_TITLE | TAB_CENTER,
+ _("Asymp. %g%% Confidence Interval"), roc->ci);
+ tab_vline (tbl, 0, n_cols - 1, 0, 0);
+ tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1);
+ }
+
+ if ( roc->n_vars > 1)
+ tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test"));
+
+ if ( roc->n_vars > 1)
+ tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
+
+
+ for ( i = 0 ; i < roc->n_vars ; ++i )
+ {
+ tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i]));
+
+ tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL);
+
+ if ( roc->print_se )
+ {
+ double se ;
+ const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
+ (12 * rs[i].n1 * rs[i].n2));
+ double ci ;
+ double yy ;
+
+ se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (rs[i].q1hat - pow2 (rs[i].auc)) +
+ (rs[i].n2 - 1) * (rs[i].q2hat - pow2 (rs[i].auc));
+
+ se /= rs[i].n1 * rs[i].n2;
+
+ se = sqrt (se);
+
+ tab_double (tbl, n_cols - 4, 2 + i, 0,
+ se,
+ NULL);
+
+ ci = 1 - roc->ci / 100.0;
+ yy = gsl_cdf_gaussian_Qinv (ci, se) ;
+
+ tab_double (tbl, n_cols - 2, 2 + i, 0,
+ rs[i].auc - yy,
+ NULL);
+
+ tab_double (tbl, n_cols - 1, 2 + i, 0,
+ rs[i].auc + yy,
+ NULL);
+
+ tab_double (tbl, n_cols - 3, 2 + i, 0,
+ 2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)),
+ NULL);
+ }
+ }
+
+ tab_submit (tbl);
+}
+
+
+static void
+show_summary (const struct cmd_roc *roc)
+{
+ const int n_cols = 3;
+ const int n_rows = 4;
+ struct tab_table *tbl = tab_create (n_cols, n_rows);
+
+ tab_title (tbl, _("Case Summary"));
+
+ tab_headers (tbl, 1, 0, 2, 0);
+
+ tab_box (tbl,
+ TAL_2, TAL_2,
+ -1, -1,
+ 0, 0,
+ n_cols - 1,
+ n_rows - 1);
+
+ tab_hline (tbl, TAL_2, 0, n_cols - 1, 2);
+ tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
+
+
+ tab_hline (tbl, TAL_2, 1, n_cols - 1, 1);
+ tab_vline (tbl, TAL_1, 2, 1, n_rows - 1);
+
+
+ tab_text (tbl, 0, 1, TAT_TITLE | TAB_LEFT, var_to_string (roc->state_var));
+ tab_text (tbl, 1, 1, TAT_TITLE, _("Unweighted"));
+ tab_text (tbl, 2, 1, TAT_TITLE, _("Weighted"));
+
+ tab_joint_text (tbl, 1, 0, 2, 0,
+ TAT_TITLE | TAB_CENTER,
+ _("Valid N (listwise)"));
+
+
+ tab_text (tbl, 0, 2, TAB_LEFT, _("Positive"));
+ tab_text (tbl, 0, 3, TAB_LEFT, _("Negative"));
+
+
+ tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0);
+ tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0);
+
+ tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0);
+ tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0);
+
+ tab_submit (tbl);
+}
+
+
+static void
+show_coords (struct roc_state *rs, const struct cmd_roc *roc)
+{
+ int x = 1;
+ int i;
+ const int n_cols = roc->n_vars > 1 ? 4 : 3;
+ int n_rows = 1;
+ struct tab_table *tbl ;
+
+ for (i = 0; i < roc->n_vars; ++i)
+ n_rows += casereader_count_cases (rs[i].cutpoint_rdr);
+
+ tbl = tab_create (n_cols, n_rows);
+
+ if ( roc->n_vars > 1)
+ tab_title (tbl, _("Coordinates of the Curve"));
+ else
+ tab_title (tbl, _("Coordinates of the Curve (%s)"), var_to_string (roc->vars[0]));
+
+
+ tab_headers (tbl, 1, 0, 1, 0);
+
+ tab_hline (tbl, TAL_2, 0, n_cols - 1, 1);
+
+ if ( roc->n_vars > 1)
+ tab_text (tbl, 0, 0, TAT_TITLE, _("Test variable"));
+
+ tab_text (tbl, n_cols - 3, 0, TAT_TITLE, _("Positive if greater than or equal to"));
+ tab_text (tbl, n_cols - 2, 0, TAT_TITLE, _("Sensitivity"));
+ tab_text (tbl, n_cols - 1, 0, TAT_TITLE, _("1 - Specificity"));
+
+ tab_box (tbl,
+ TAL_2, TAL_2,
+ -1, TAL_1,
+ 0, 0,
+ n_cols - 1,
+ n_rows - 1);
+
+ if ( roc->n_vars > 1)
+ tab_vline (tbl, TAL_2, 1, 0, n_rows - 1);
+
+ for (i = 0; i < roc->n_vars; ++i)
+ {
+ struct ccase *cc;
+ struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
+
+ if ( roc->n_vars > 1)
+ tab_text (tbl, 0, x, TAT_TITLE, var_to_string (roc->vars[i]));
+
+ if ( i > 0)
+ tab_hline (tbl, TAL_1, 0, n_cols - 1, x);
+
+
+ for (; (cc = casereader_read (r)) != NULL;
+ case_unref (cc), x++)
+ {
+ const double se = case_data_idx (cc, ROC_TP)->f /
+ (
+ case_data_idx (cc, ROC_TP)->f
+ +
+ case_data_idx (cc, ROC_FN)->f
+ );
+
+ const double sp = case_data_idx (cc, ROC_TN)->f /
+ (
+ case_data_idx (cc, ROC_TN)->f
+ +
+ case_data_idx (cc, ROC_FP)->f
+ );
+
+ tab_double (tbl, n_cols - 3, x, 0, case_data_idx (cc, ROC_CUTPOINT)->f,
+ var_get_print_format (roc->vars[i]));
+
+ tab_double (tbl, n_cols - 2, x, 0, se, NULL);
+ tab_double (tbl, n_cols - 1, x, 0, 1 - sp, NULL);
+ }
+
+ casereader_destroy (r);
+ }
+
+ tab_submit (tbl);
+}
+
+
+static void
+output_roc (struct roc_state *rs, const struct cmd_roc *roc)
+{
+ show_summary (roc);
+
+ if ( roc->curve )
+ {
+ struct roc_chart *rc;
+ size_t i;
+
+ rc = roc_chart_create (roc->reference);
+ for (i = 0; i < roc->n_vars; i++)
+ roc_chart_add_var (rc, var_get_name (roc->vars[i]),
+ rs[i].cutpoint_rdr);
+ roc_chart_submit (rc);
+ }
+
+ show_auc (rs, roc);
+
+ if ( roc->print_coords )
+ show_coords (rs, roc);