+ struct ccase *cc = case_create (casewriter_get_proto (writer));
+
+ *case_num_rw_idx (cc, ROC_CUTPOINT) = cutpoint;
+ *case_num_rw_idx (cc, ROC_TP) = 0;
+ *case_num_rw_idx (cc, ROC_FN) = 0;
+ *case_num_rw_idx (cc, ROC_TN) = 0;
+ *case_num_rw_idx (cc, ROC_FP) = 0;
+
+ casewriter_write (writer, cc);
+}
+
+
+/*
+ Create and initialise the rs[x].cutpoint_rdr casereaders. That is, the readers will
+ be created with width 5, ready to take the values (cutpoint, ROC_TP, ROC_FN, ROC_TN, ROC_FP), and the
+ reader will be populated with its final number of cases.
+ However on exit from this function, only ROC_CUTPOINT entries will be set to their final
+ value. The other entries will be initialised to zero.
+*/
+static void
+prepare_cutpoints (struct cmd_roc *roc, struct roc_state *rs, struct casereader *input)
+{
+ int i;
+ struct casereader *r = casereader_clone (input);
+ struct ccase *c;
+
+ {
+ struct caseproto *proto = caseproto_create ();
+ struct subcase ordering;
+ subcase_init (&ordering, ROC_CUTPOINT, 0, SC_ASCEND);
+
+ proto = caseproto_add_width (proto, 0); /* cutpoint */
+ proto = caseproto_add_width (proto, 0); /* ROC_TP */
+ proto = caseproto_add_width (proto, 0); /* ROC_FN */
+ proto = caseproto_add_width (proto, 0); /* ROC_TN */
+ proto = caseproto_add_width (proto, 0); /* ROC_FP */
+
+ for (i = 0 ; i < roc->n_vars; ++i)
+ {
+ rs[i].cutpoint_wtr = sort_create_writer (&ordering, proto);
+ rs[i].prev_result = SYSMIS;
+ rs[i].max = -DBL_MAX;
+ rs[i].min = DBL_MAX;
+ }
+
+ caseproto_unref (proto);
+ subcase_uninit (&ordering);
+ }
+
+ for (; (c = casereader_read (r)) != NULL; case_unref (c))
+ {
+ for (i = 0 ; i < roc->n_vars; ++i)
+ {
+ const union value *v = case_data (c, roc->vars[i]);
+ const double result = v->f;
+
+ if (mv_is_value_missing (var_get_missing_values (roc->vars[i]), v)
+ & roc->exclude)
+ continue;
+
+ minimize (&rs[i].min, result);
+ maximize (&rs[i].max, result);
+
+ if (rs[i].prev_result != SYSMIS && rs[i].prev_result != result)
+ {
+ const double mean = (result + rs[i].prev_result) / 2.0;
+ append_cutpoint (rs[i].cutpoint_wtr, mean);
+ }
+
+ rs[i].prev_result = result;
+ }
+ }
+ casereader_destroy (r);
+
+
+ /* Append the min and max cutpoints */
+ for (i = 0 ; i < roc->n_vars; ++i)
+ {
+ append_cutpoint (rs[i].cutpoint_wtr, rs[i].min - 1);
+ append_cutpoint (rs[i].cutpoint_wtr, rs[i].max + 1);
+
+ rs[i].cutpoint_rdr = casewriter_make_reader (rs[i].cutpoint_wtr);
+ }
+}
+
+static void
+do_roc (struct cmd_roc *roc, struct casereader *reader, struct dictionary *dict)
+{
+ int i;
+
+ struct roc_state *rs = XCALLOC (roc->n_vars, struct roc_state);
+
+ 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_num_idx (cpos, VALUE);
+
+ while ((cneg = casereader_read (n_neg_reader)))
+ {
+ struct ccase *nc = case_create (n_proto);
+
+ const double jneg = case_num_idx (cneg, VALUE);
+
+ *case_num_rw_idx (nc, VALUE) = jneg;
+ *case_num_rw_idx (nc, N_POS_EQ) = 0;
+
+ *case_num_rw_idx (nc, N_POS_GT) = 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_num_rw_idx (pos_case, VALUE) = 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_num_rw_idx (pos_case, N_NEG_EQ) = 0;
+ *case_num_rw_idx (pos_case, N_NEG_LT) = 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_num_idx (c, N_POS_GT);
+ struct ccase *nc = case_clone (c);
+
+ if (n_pos_gt == SYSMIS)
+ {
+ n_pos_gt = prev_pos_gt;
+ *case_num_rw_idx (nc, N_POS_GT) = 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_num_idx (c, N_NEG_LT);
+ struct ccase *nc = case_clone (c);
+
+ if (n_neg_lt == SYSMIS)
+ {
+ n_neg_lt = prev_neg_lt;
+ *case_num_rw_idx (nc, N_NEG_LT) = 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_num_idx (c, VALUE);
+ double n_pos_eq = case_num_idx (c, N_POS_EQ);
+ double n_pos_gt = case_num_idx (c, N_POS_GT);
+ double n_neg_eq = case_num_idx (c, N_NEG_EQ);
+ double n_neg_lt = case_num_idx (c, N_NEG_LT);
+
+ if (prev_case && j == case_num_idx (prev_case, VALUE))
+ {
+ if (0 == case_num_idx (c, N_POS_EQ))
+ {
+ n_pos_eq = case_num_idx (prev_case, N_POS_EQ);
+ n_pos_gt = case_num_idx (prev_case, N_POS_GT);
+ }
+
+ if (0 == case_num_idx (c, N_NEG_EQ))
+ {
+ n_neg_eq = case_num_idx (prev_case, N_NEG_EQ);
+ n_neg_lt = case_num_idx (prev_case, N_NEG_LT);
+ }
+ }
+
+ if (NULL == next_case || j != case_num_idx (next_case, VALUE))
+ {
+ 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_uninit (&up_ordering);
+ subcase_uninit (&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)
+{
+ struct pivot_table *table = pivot_table_create (N_("Area Under the Curve"));
+
+ struct pivot_dimension *statistics = pivot_dimension_create (
+ table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+ N_("Area"), PIVOT_RC_OTHER);
+ if (roc->print_se)
+ {
+ pivot_category_create_leaves (
+ statistics->root,
+ N_("Std. Error"), PIVOT_RC_OTHER,
+ N_("Asymptotic Sig."), PIVOT_RC_SIGNIFICANCE);
+ struct pivot_category *interval = pivot_category_create_group__ (
+ statistics->root,
+ pivot_value_new_text_format (N_("Asymp. %g%% Confidence Interval"),
+ roc->ci));
+ pivot_category_create_leaves (interval,
+ N_("Lower Bound"), PIVOT_RC_OTHER,
+ N_("Upper Bound"), PIVOT_RC_OTHER);
+ }
+
+ struct pivot_dimension *variables = pivot_dimension_create (
+ table, PIVOT_AXIS_ROW, N_("Variable under test"));
+ variables->root->show_label = true;
+
+ for (size_t i = 0 ; i < roc->n_vars ; ++i)
+ {
+ int var_idx = pivot_category_create_leaf (
+ variables->root, pivot_value_new_variable (roc->vars[i]));
+
+ pivot_table_put2 (table, 0, var_idx, pivot_value_new_number (rs[i].auc));
+
+ if (roc->print_se)
+ {
+ double 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);
+
+ double ci = 1 - roc->ci / 100.0;
+ double yy = gsl_cdf_gaussian_Qinv (ci, se);
+
+ double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) /
+ (12 * rs[i].n1 * rs[i].n2));
+ double sig = 2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5)
+ / sd_0_5));
+ double entries[] = { se, sig, rs[i].auc - yy, rs[i].auc + yy };
+ for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
+ pivot_table_put2 (table, i + 1, var_idx,
+ pivot_value_new_number (entries[i]));
+ }
+ }
+
+ pivot_table_submit (table);
+}
+
+
+static void
+show_summary (const struct cmd_roc *roc)
+{
+ struct pivot_table *table = pivot_table_create (N_("Case Summary"));
+
+ struct pivot_dimension *statistics = pivot_dimension_create (
+ table, PIVOT_AXIS_COLUMN, N_("Valid N (listwise)"),
+ N_("Unweighted"), PIVOT_RC_INTEGER,
+ N_("Weighted"), PIVOT_RC_OTHER);
+ statistics->root->show_label = true;
+
+ struct pivot_dimension *cases = pivot_dimension_create__ (
+ table, PIVOT_AXIS_ROW, pivot_value_new_variable (roc->state_var));
+ cases->root->show_label = true;
+ pivot_category_create_leaves (cases->root, N_("Positive"), N_("Negative"));
+
+ struct entry
+ {
+ int stat_idx;
+ int case_idx;
+ double x;
+ }
+ entries[] = {
+ { 0, 0, roc->pos },
+ { 0, 1, roc->neg },
+ { 1, 0, roc->pos_weighted },
+ { 1, 1, roc->neg_weighted },
+ };
+ for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
+ {
+ const struct entry *e = &entries[i];
+ pivot_table_put2 (table, e->stat_idx, e->case_idx,
+ pivot_value_new_number (e->x));
+ }
+ pivot_table_submit (table);
+}
+
+static void
+show_coords (struct roc_state *rs, const struct cmd_roc *roc)
+{
+ struct pivot_table *table = pivot_table_create (
+ N_("Coordinates of the Curve"));
+
+ pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+ N_("Positive if greater than or equal to"),
+ N_("Sensitivity"), N_("1 - Specificity"));
+
+ struct pivot_dimension *coordinates = pivot_dimension_create (
+ table, PIVOT_AXIS_ROW, N_("Coordinates"));
+ coordinates->hide_all_labels = true;
+
+ struct pivot_dimension *variables = pivot_dimension_create (
+ table, PIVOT_AXIS_ROW, N_("Test variable"));
+ variables->root->show_label = true;
+
+
+ int n_coords = 0;
+ for (size_t i = 0; i < roc->n_vars; ++i)
+ {
+ struct casereader *r = casereader_clone (rs[i].cutpoint_rdr);
+
+ int var_idx = pivot_category_create_leaf (
+ variables->root, pivot_value_new_variable (roc->vars[i]));
+
+ struct ccase *cc;
+ int coord_idx = 0;
+ for (; (cc = casereader_read (r)) != NULL; case_unref (cc))
+ {
+ const double se = case_num_idx (cc, ROC_TP) /
+ (case_num_idx (cc, ROC_TP) + case_num_idx (cc, ROC_FN));
+
+ const double sp = case_num_idx (cc, ROC_TN) /
+ (case_num_idx (cc, ROC_TN) + case_num_idx (cc, ROC_FP));
+
+ if (coord_idx >= n_coords)
+ {
+ assert (coord_idx == n_coords);
+ pivot_category_create_leaf (
+ coordinates->root, pivot_value_new_integer (++n_coords));
+ }
+
+ pivot_table_put3 (
+ table, 0, coord_idx, var_idx,
+ pivot_value_new_var_value (roc->vars[i],
+ case_data_idx (cc, ROC_CUTPOINT)));
+
+ pivot_table_put3 (table, 1, coord_idx, var_idx,
+ pivot_value_new_number (se));
+ pivot_table_put3 (table, 2, coord_idx, var_idx,
+ pivot_value_new_number (1 - sp));
+ coord_idx++;
+ }
+
+ casereader_destroy (r);
+ }
+
+ pivot_table_submit (table);
+}
+
+
+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);