src/examine.q

   1 /* PSPP - EXAMINE data for normality . -*-c-*-
   2
   3 Copyright (C) 2004 Free Software Foundation, Inc.
   4 Author: John Darrington 2004
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License as
   8 published by the Free Software Foundation; either version 2 of the
   9 License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program; if not, write to the Free Software
  18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19 02111-1307, USA. */
  20
  21 #include <config.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include "error.h"
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <math.h>
  27 #include "alloc.h"
  28 #include "str.h"
  29 #include "case.h"
  30 #include "dictionary.h"
  31 #include "command.h"
  32 #include "lexer.h"
  33 #include "error.h"
  34 #include "magic.h"
  35 #include "misc.h"
  36 #include "tab.h"
  37 #include "som.h"
  38 #include "value-labels.h"
  39 #include "var.h"
  40 #include "vfm.h"
  41 #include "hash.h"
  42 #include "casefile.h"
  43 #include "factor_stats.h"
  44 /* (headers) */
  45 #include "chart.h"
  46
  47 /* (specification)
  48    "EXAMINE" (xmn_):
  49    *variables=custom;
  50    +total=custom;
  51    +nototal=custom;
  52    +missing=miss:pairwise/!listwise,
  53    rep:report/!noreport,
  54    incl:include/!exclude;
  55    +compare=cmp:variables/!groups;
  56    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  57    +cinterval=double;
  58    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  59 */
  60
  61 /* (declarations) */
  62
  63 /* (functions) */
  64
  65
  66
  67 static struct cmd_examine cmd;
  68
  69 static struct variable **dependent_vars;
  70
  71 static int n_dependent_vars;
  72
  73
  74 struct factor
  75 {
  76   /* The independent variable */
  77   struct variable *indep_var[2];
  78
  79
  80   /* Hash table of factor stats indexed by 2 values */
  81   struct hsh_table *fstats;
  82
  83   /* The hash table after it has been crunched */
  84   struct factor_statistics **fs;
  85
  86   struct factor *next;
  87
  88 };
  89
  90 /* Linked list of factors */
  91 static struct factor *factors=0;
  92
  93 static struct metrics *totals=0;
  94
  95 void
  96 print_factors(void)
  97 {
  98   struct factor *f = factors;
  99
 100   while (f)
 101     {
 102       struct  factor_statistics **fs = f->fs;
 103
 104       printf("Factor: %s BY %s\n",
 105              var_to_string(f->indep_var[0]),
 106              var_to_string(f->indep_var[1]) );
 107
 108
 109       printf("Contains %d entries\n", hsh_count(f->fstats));
 110
 111
 112       while (*fs)
 113         {
 114           printf("Factor %g; %g\n", (*fs)->id[0].f, (*fs)->id[1].f);
 115
 116           /*
 117              printf("Factor %s; %s\n",
 118              value_to_string(&(*fs)->id[0], f->indep_var[0]),
 119              value_to_string(&(*fs)->id[1], f->indep_var[1]));
 120           */
 121
 122
 123           printf("Sum is %g; ",(*fs)->m[0].sum);
 124           printf("N is %g; ",(*fs)->m[0].n);
 125           printf("Mean is %g\n",(*fs)->m[0].mean);
 126
 127           fs++ ;
 128         }
 129
 130       f = f->next;
 131     }
 132
 133
 134 }
 135
 136
 137 /* Parse the clause specifying the factors */
 138 static int examine_parse_independent_vars(struct cmd_examine *cmd);
 139
 140
 141
 142 /* Output functions */
 143 static void show_summary(struct variable **dependent_var, int n_dep_var,
 144                          const struct factor *f);
 145
 146 static void show_extremes(struct variable **dependent_var,
 147                           int n_dep_var,
 148                           const struct factor *factor,
 149                           int n_extremities);
 150
 151 static void show_descriptives(struct variable **dependent_var,
 152                               int n_dep_var,
 153                               struct factor *factor);
 154
 155
 156 void np_plot(const struct metrics *m, const char *factorname);
 157
 158
 159
 160
 161 /* Per Split function */
 162 static void run_examine(const struct casefile *cf, void *cmd_);
 163
 164 static void output_examine(void);
 165
 166
 167 void factor_calc(struct ccase *c, int case_no,
 168                  double weight, int case_missing);
 169
 170
 171 /* Function to use for testing for missing values */
 172 static is_missing_func value_is_missing;
 173
 174
 175 int
 176 cmd_examine(void)
 177 {
 178
 179   if ( !parse_examine(&cmd) )
 180     return CMD_FAILURE;
 181
 182   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 183   if (cmd.incl == XMN_INCLUDE )
 184     value_is_missing = is_system_missing;
 185   else
 186     value_is_missing = is_missing;
 187
 188   if ( cmd.st_n == SYSMIS )
 189     cmd.st_n = 5;
 190
 191   if ( ! cmd.sbc_cinterval)
 192     cmd.n_cinterval[0] = 95.0;
 193
 194   multipass_procedure_with_splits (run_examine, &cmd);
 195
 196   if ( totals )
 197     free(totals);
 198
 199   return CMD_SUCCESS;
 200 };
 201
 202
 203
 204 /* Show all the appropriate tables */
 205 static void
 206 output_examine(void)
 207 {
 208   struct factor *fctr;
 209
 210   /* Show totals if appropriate */
 211   if ( ! cmd.sbc_nototal || factors == 0 )
 212     {
 213       show_summary(dependent_vars, n_dependent_vars, 0);
 214
 215       if ( cmd.sbc_statistics )
 216         {
 217           if ( cmd.a_statistics[XMN_ST_EXTREME])
 218             show_extremes(dependent_vars, n_dependent_vars, 0, cmd.st_n);
 219
 220           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 221             show_descriptives(dependent_vars, n_dependent_vars, 0);
 222
 223         }
 224
 225       if ( cmd.sbc_plot)
 226         {
 227           if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 228             {
 229               int v;
 230
 231               for ( v = 0 ; v < n_dependent_vars; ++v )
 232                   np_plot(&totals[v], var_to_string(dependent_vars[v]));
 233             }
 234         }
 235
 236
 237     }
 238
 239
 240   /* Show grouped statistics  as appropriate */
 241   fctr = factors;
 242   while ( fctr )
 243     {
 244       show_summary(dependent_vars, n_dependent_vars, fctr);
 245
 246       if ( cmd.sbc_statistics )
 247         {
 248           if ( cmd.a_statistics[XMN_ST_EXTREME])
 249             show_extremes(dependent_vars, n_dependent_vars, fctr, cmd.st_n);
 250
 251           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 252             show_descriptives(dependent_vars, n_dependent_vars, fctr);
 253         }
 254
 255       if ( cmd.sbc_plot)
 256         {
 257           if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 258             {
 259               int v;
 260               for ( v = 0 ; v < n_dependent_vars; ++ v)
 261                 {
 262
 263                   struct factor_statistics **fs = fctr->fs ;
 264                   for ( fs = fctr->fs ; *fs ; ++fs )
 265                     {
 266                       char buf1[100];
 267                       char buf2[100];
 268                       sprintf(buf1, "%s (",
 269                               var_to_string(dependent_vars[v]));
 270
 271                       sprintf(buf2, "%s = %s",
 272                              var_to_string(fctr->indep_var[0]),
 273                              value_to_string(&(*fs)->id[0],fctr->indep_var[0]));
 274
 275                       strcat(buf1, buf2);
 276
 277
 278                       if ( fctr->indep_var[1] )
 279                         {
 280                           sprintf(buf2, "; %s = %s)",
 281                                   var_to_string(fctr->indep_var[1]),
 282                                   value_to_string(&(*fs)->id[1],
 283                                                   fctr->indep_var[1]));
 284                           strcat(buf1, buf2);
 285                         }
 286                       else
 287                         {
 288                           strcat(buf1, ")");
 289                         }
 290
 291                       np_plot(&(*fs)->m[v],buf1);
 292
 293                     }
 294
 295                 }
 296
 297             }
 298         }
 299
 300       fctr = fctr->next;
 301     }
 302
 303 }
 304
 305
 306
 307 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
 308 static int
 309 xmn_custom_total(struct cmd_examine *p)
 310 {
 311   if ( p->sbc_nototal )
 312     {
 313       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 314       return 0;
 315     }
 316
 317   return 1;
 318 }
 319
 320 static int
 321 xmn_custom_nototal(struct cmd_examine *p)
 322 {
 323   if ( p->sbc_total )
 324     {
 325       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 326       return 0;
 327     }
 328
 329   return 1;
 330 }
 331
 332
 333
 334 /* Parser for the variables sub command */
 335 static int
 336 xmn_custom_variables(struct cmd_examine *cmd )
 337 {
 338
 339   lex_match('=');
 340
 341   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 342       && token != T_ALL)
 343     return 2;
 344
 345   if (!parse_variables (default_dict, &dependent_vars, &n_dependent_vars,
 346                         PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
 347     {
 348       free (dependent_vars);
 349       return 0;
 350     }
 351
 352   assert(n_dependent_vars);
 353
 354   totals = xmalloc( sizeof(struct metrics) * n_dependent_vars);
 355
 356   if ( lex_match(T_BY))
 357     {
 358       return examine_parse_independent_vars(cmd);
 359     }
 360
 361   return 1;
 362 }
 363
 364
 365
 366 /* Parse the clause specifying the factors */
 367 static int
 368 examine_parse_independent_vars(struct cmd_examine *cmd)
 369 {
 370
 371   struct factor *sf = xmalloc(sizeof(struct factor));
 372
 373   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 374       && token != T_ALL)
 375     return 2;
 376
 377
 378   sf->indep_var[0] = parse_variable();
 379   sf->indep_var[1] = 0;
 380
 381   if ( token == T_BY )
 382     {
 383
 384       lex_match(T_BY);
 385
 386       if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 387           && token != T_ALL)
 388         return 2;
 389
 390       sf->indep_var[1] = parse_variable();
 391
 392     }
 393
 394
 395   sf->fstats = hsh_create(4,
 396                           (hsh_compare_func *) factor_statistics_compare,
 397                           (hsh_hash_func *) factor_statistics_hash,
 398                           (hsh_free_func *) factor_statistics_free,
 399                           0);
 400
 401   sf->next = factors;
 402   factors = sf;
 403
 404   lex_match(',');
 405
 406   if ( token == '.' || token == '/' )
 407     return 1;
 408
 409   return examine_parse_independent_vars(cmd);
 410 }
 411
 412
 413
 414
 415 void populate_descriptives(struct tab_table *t, int col, int row,
 416                            const struct metrics *fs);
 417
 418 void populate_extremes(struct tab_table *t, int col, int row, int n,
 419                        const struct metrics *m);
 420
 421 void populate_summary(struct tab_table *t, int col, int row,
 422                       const struct metrics *m);
 423
 424
 425
 426
 427 static int bad_weight_warn = 1;
 428
 429
 430 /* Perform calculations for the sub factors */
 431 void
 432 factor_calc(struct ccase *c, int case_no, double weight, int case_missing)
 433 {
 434   int v;
 435   struct factor *fctr = factors;
 436
 437   while ( fctr)
 438     {
 439       union value indep_vals[2] ;
 440
 441       indep_vals[0] = * case_data(c, fctr->indep_var[0]->fv);
 442
 443       if ( fctr->indep_var[1] )
 444         indep_vals[1] = * case_data(c, fctr->indep_var[1]->fv);
 445       else
 446         indep_vals[1].f = SYSMIS;
 447
 448       assert(fctr->fstats);
 449
 450       struct factor_statistics **foo = ( struct factor_statistics ** )
 451         hsh_probe(fctr->fstats, (void *) &indep_vals);
 452
 453       if ( !*foo )
 454         {
 455
 456           *foo = create_factor_statistics(n_dependent_vars,
 457                                           &indep_vals[0],
 458                                           &indep_vals[1]);
 459
 460           for ( v =  0 ; v  < n_dependent_vars ; ++v )
 461             {
 462               metrics_precalc( &(*foo)->m[v] );
 463             }
 464
 465         }
 466
 467       for ( v =  0 ; v  < n_dependent_vars ; ++v )
 468         {
 469           const struct variable *var = dependent_vars[v];
 470           const union value *val = case_data (c, var->fv);
 471
 472           if ( value_is_missing(val,var) || case_missing )
 473             val = 0;
 474
 475           metrics_calc( &(*foo)->m[v], val, weight, case_no );
 476         }
 477
 478       fctr = fctr->next;
 479     }
 480
 481
 482 }
 483
 484
 485
 486
 487 static void
 488 run_examine(const struct casefile *cf, void *cmd_ )
 489 {
 490   struct casereader *r;
 491   struct ccase c;
 492   int v;
 493
 494   const struct cmd_examine *cmd = (struct cmd_examine *) cmd_;
 495
 496   /* Make sure we haven't got rubbish left over from a
 497      previous split */
 498
 499   struct factor *fctr = factors;
 500   while (fctr)
 501     {
 502       struct factor *next = fctr->next;
 503
 504       hsh_clear(fctr->fstats);
 505
 506       fctr->fs = 0;
 507
 508       fctr = next;
 509     }
 510
 511
 512
 513   for ( v = 0 ; v < n_dependent_vars ; ++v )
 514     metrics_precalc(&totals[v]);
 515
 516   for(r = casefile_get_reader (cf);
 517       casereader_read (r, &c) ;
 518       case_destroy (&c) )
 519     {
 520       int case_missing=0;
 521       const int case_no = casereader_cnum(r);
 522
 523       const double weight =
 524         dict_get_case_weight(default_dict, &c, &bad_weight_warn);
 525
 526       if ( cmd->miss == XMN_LISTWISE )
 527         {
 528           for ( v = 0 ; v < n_dependent_vars ; ++v )
 529             {
 530               const struct variable *var = dependent_vars[v];
 531               const union value *val = case_data (&c, var->fv);
 532
 533               if ( value_is_missing(val,var))
 534                 case_missing = 1;
 535
 536             }
 537         }
 538
 539       for ( v = 0 ; v < n_dependent_vars ; ++v )
 540         {
 541           const struct variable *var = dependent_vars[v];
 542           const union value *val = case_data (&c, var->fv);
 543
 544           if ( value_is_missing(val,var) || case_missing )
 545             val = 0;
 546
 547           metrics_calc(&totals[v], val, weight, case_no );
 548
 549         }
 550
 551       factor_calc(&c, case_no, weight, case_missing);
 552
 553     }
 554
 555
 556   for ( v = 0 ; v < n_dependent_vars ; ++v)
 557     {
 558       fctr = factors;
 559       while ( fctr )
 560         {
 561           struct hsh_iterator hi;
 562           struct factor_statistics *fs;
 563
 564           for ( fs = hsh_first(fctr->fstats, &hi);
 565                 fs != 0 ;
 566                 fs = hsh_next(fctr->fstats, &hi))
 567             {
 568               metrics_postcalc(&fs->m[v]);
 569             }
 570
 571           fctr = fctr->next;
 572         }
 573       metrics_postcalc(&totals[v]);
 574     }
 575
 576
 577   /* Make sure that the combination of factors are complete */
 578
 579   fctr = factors;
 580   while ( fctr )
 581     {
 582       struct hsh_iterator hi;
 583       struct hsh_iterator hi0;
 584       struct hsh_iterator hi1;
 585       struct factor_statistics *fs;
 586
 587       struct hsh_table *idh0=0;
 588       struct hsh_table *idh1=0;
 589       union value *val0;
 590       union value *val1;
 591
 592       idh0 = hsh_create(4, (hsh_compare_func *) compare_values,
 593                         (hsh_hash_func *) hash_value,
 594                         0,0);
 595
 596       idh1 = hsh_create(4, (hsh_compare_func *) compare_values,
 597                         (hsh_hash_func *) hash_value,
 598                         0,0);
 599
 600
 601       for ( fs = hsh_first(fctr->fstats, &hi);
 602             fs != 0 ;
 603             fs = hsh_next(fctr->fstats, &hi))
 604         {
 605           hsh_insert(idh0,(void *) &fs->id[0]);
 606           hsh_insert(idh1,(void *) &fs->id[1]);
 607         }
 608
 609       /* Ensure that the factors combination is complete */
 610       for ( val0 = hsh_first(idh0, &hi0);
 611             val0 != 0 ;
 612             val0 = hsh_next(idh0, &hi0))
 613         {
 614           for ( val1 = hsh_first(idh1, &hi1);
 615                 val1 != 0 ;
 616                 val1 = hsh_next(idh1, &hi1))
 617             {
 618               struct factor_statistics **ffs;
 619               union value key[2];
 620               key[0] = *val0;
 621               key[1] = *val1;
 622
 623               ffs = (struct factor_statistics **)
 624                 hsh_probe(fctr->fstats, (void *) &key );
 625
 626               if ( !*ffs ) {
 627                 int i;
 628                 (*ffs) = create_factor_statistics (n_dependent_vars,
 629                                                    &key[0], &key[1]);
 630                 for ( i = 0 ; i < n_dependent_vars ; ++i )
 631                   metrics_precalc( &(*ffs)->m[i]);
 632               }
 633             }
 634         }
 635
 636       hsh_destroy(idh0);
 637       hsh_destroy(idh1);
 638
 639       fctr->fs = (struct factor_statistics **) hsh_sort_copy(fctr->fstats);
 640
 641       fctr = fctr->next;
 642     }
 643
 644   output_examine();
 645
 646   for ( v = 0 ; v < n_dependent_vars ; ++v )
 647     hsh_destroy(totals[v].ordered_data);
 648
 649 }
 650
 651
 652 static void
 653 show_summary(struct variable **dependent_var, int n_dep_var,
 654              const struct factor *fctr)
 655 {
 656   static const char *subtitle[]=
 657     {
 658       N_("Valid"),
 659       N_("Missing"),
 660       N_("Total")
 661     };
 662
 663   int i;
 664   int heading_columns ;
 665   int n_cols;
 666   const int heading_rows = 3;
 667   struct tab_table *tbl;
 668
 669   int n_rows ;
 670   int n_factors = 1;
 671
 672   if ( fctr )
 673     {
 674       heading_columns = 2;
 675       n_factors = hsh_count(fctr->fstats);
 676       n_rows = n_dep_var * n_factors ;
 677
 678       if ( fctr->indep_var[1] )
 679           heading_columns = 3;
 680     }
 681   else
 682     {
 683       heading_columns = 1;
 684       n_rows = n_dep_var;
 685     }
 686
 687   n_rows += heading_rows;
 688
 689   n_cols = heading_columns + 6;
 690
 691   tbl = tab_create (n_cols,n_rows,0);
 692   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 693
 694   tab_dim (tbl, tab_natural_dimensions);
 695
 696   /* Outline the box */
 697   tab_box (tbl,
 698            TAL_2, TAL_2,
 699            -1, -1,
 700            0, 0,
 701            n_cols - 1, n_rows - 1);
 702
 703   /* Vertical lines for the data only */
 704   tab_box (tbl,
 705            -1, -1,
 706            -1, TAL_1,
 707            heading_columns, 0,
 708            n_cols - 1, n_rows - 1);
 709
 710
 711   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
 712   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
 713   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
 714
 715   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
 716
 717
 718   tab_title (tbl, 0, _("Case Processing Summary"));
 719
 720
 721   tab_joint_text(tbl, heading_columns, 0,
 722                  n_cols -1, 0,
 723                  TAB_CENTER | TAT_TITLE,
 724                  _("Cases"));
 725
 726   /* Remove lines ... */
 727   tab_box (tbl,
 728            -1, -1,
 729            TAL_0, TAL_0,
 730            heading_columns, 0,
 731            n_cols - 1, 0);
 732
 733   for ( i = 0 ; i < 3 ; ++i )
 734     {
 735       tab_text (tbl, heading_columns + i*2 , 2, TAB_CENTER | TAT_TITLE,
 736                 _("N"));
 737
 738       tab_text (tbl, heading_columns + i*2 + 1, 2, TAB_CENTER | TAT_TITLE,
 739                 _("Percent"));
 740
 741       tab_joint_text(tbl, heading_columns + i*2 , 1,
 742                      heading_columns + i*2 + 1, 1,
 743                      TAB_CENTER | TAT_TITLE,
 744                      subtitle[i]);
 745
 746       tab_box (tbl, -1, -1,
 747                TAL_0, TAL_0,
 748                heading_columns + i*2, 1,
 749                heading_columns + i*2 + 1, 1);
 750
 751     }
 752
 753
 754   /* Titles for the independent variables */
 755   if ( fctr )
 756     {
 757       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 758                 var_to_string(fctr->indep_var[0]));
 759
 760       if ( fctr->indep_var[1] )
 761         {
 762           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 763                     var_to_string(fctr->indep_var[1]));
 764         }
 765
 766     }
 767
 768
 769   for ( i = 0 ; i < n_dep_var ; ++i )
 770     {
 771       int n_factors = 1;
 772       if ( fctr )
 773         n_factors = hsh_count(fctr->fstats);
 774
 775
 776       if ( i > 0 )
 777         tab_hline(tbl, TAL_1, 0, n_cols -1 , i * n_factors + heading_rows);
 778
 779       tab_text (tbl,
 780                 0, i * n_factors + heading_rows,
 781                 TAB_LEFT | TAT_TITLE,
 782                 var_to_string(dependent_var[i])
 783                 );
 784
 785
 786       if ( !fctr )
 787         populate_summary(tbl, heading_columns,
 788                          (i * n_factors) + heading_rows,
 789                          &totals[i]);
 790
 791
 792       else
 793         {
 794           struct factor_statistics **fs = fctr->fs;
 795           int count = 0 ;
 796
 797           while (*fs)
 798             {
 799               static union value prev;
 800
 801               if ( 0 != compare_values(&prev, &(*fs)->id[0],
 802                                        fctr->indep_var[0]->width))
 803                 {
 804                    tab_text (tbl,
 805                              1,
 806                              (i * n_factors ) + count +
 807                              heading_rows,
 808                              TAB_LEFT | TAT_TITLE,
 809                              value_to_string(&(*fs)->id[0], fctr->indep_var[0])
 810                              );
 811
 812                    if (fctr->indep_var[1] && count > 0 )
 813                      tab_hline(tbl, TAL_1, 1, n_cols - 1,
 814                                (i * n_factors ) + count + heading_rows);
 815
 816                 }
 817
 818               prev = (*fs)->id[0];
 819
 820
 821               if ( fctr->indep_var[1])
 822                 tab_text (tbl,
 823                           2,
 824                           (i * n_factors ) + count +
 825                           heading_rows,
 826                           TAB_LEFT | TAT_TITLE,
 827                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
 828                           );
 829
 830               populate_summary(tbl, heading_columns,
 831                                (i * n_factors) + count
 832                                + heading_rows,
 833                                &(*fs)->m[i]);
 834
 835               count++ ;
 836               fs++;
 837             }
 838         }
 839     }
 840
 841   tab_submit (tbl);
 842 }
 843
 844
 845 void
 846 populate_summary(struct tab_table *t, int col, int row,
 847                  const struct metrics *m)
 848
 849 {
 850   const double total = m->n + m->n_missing ;
 851
 852   tab_float(t, col + 0, row + 0, TAB_RIGHT, m->n, 8, 0);
 853   tab_float(t, col + 2, row + 0, TAB_RIGHT, m->n_missing, 8, 0);
 854   tab_float(t, col + 4, row + 0, TAB_RIGHT, total, 8, 0);
 855
 856
 857   if ( total > 0 ) {
 858     tab_text (t, col + 1, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
 859               100.0 * m->n / total );
 860
 861     tab_text (t, col + 3, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
 862               100.0 * m->n_missing / total );
 863
 864     /* This seems a bit pointless !!! */
 865     tab_text (t, col + 5, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
 866               100.0 * total / total );
 867
 868
 869   }
 870
 871
 872 }
 873
 874
 875
 876 static void
 877 show_extremes(struct variable **dependent_var, int n_dep_var,
 878               const struct factor *fctr, int n_extremities)
 879 {
 880   int i;
 881   int heading_columns ;
 882   int n_cols;
 883   const int heading_rows = 1;
 884   struct tab_table *tbl;
 885
 886   int n_factors = 1;
 887   int n_rows ;
 888
 889   if ( fctr )
 890     {
 891       heading_columns = 2;
 892       n_factors = hsh_count(fctr->fstats);
 893
 894       n_rows = n_dep_var * 2 * n_extremities * n_factors;
 895
 896       if ( fctr->indep_var[1] )
 897           heading_columns = 3;
 898     }
 899   else
 900     {
 901       heading_columns = 1;
 902       n_rows = n_dep_var * 2 * n_extremities;
 903     }
 904
 905   n_rows += heading_rows;
 906
 907   heading_columns += 2;
 908   n_cols = heading_columns + 2;
 909
 910   tbl = tab_create (n_cols,n_rows,0);
 911   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 912
 913   tab_dim (tbl, tab_natural_dimensions);
 914
 915   /* Outline the box, No internal lines*/
 916   tab_box (tbl,
 917            TAL_2, TAL_2,
 918            -1, -1,
 919            0, 0,
 920            n_cols - 1, n_rows - 1);
 921
 922   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
 923
 924   tab_title (tbl, 0, _("Extreme Values"));
 925
 926
 927   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows -1);
 928   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows -1);
 929
 930   if ( fctr )
 931     {
 932       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 933                 var_to_string(fctr->indep_var[0]));
 934
 935       if ( fctr->indep_var[1] )
 936         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 937                   var_to_string(fctr->indep_var[1]));
 938     }
 939
 940   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Value"));
 941   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Case Number"));
 942
 943
 944
 945
 946   for ( i = 0 ; i < n_dep_var ; ++i )
 947     {
 948
 949       if ( i > 0 )
 950         tab_hline(tbl, TAL_1, 0, n_cols -1 ,
 951                   i * 2 * n_extremities * n_factors + heading_rows);
 952
 953       tab_text (tbl, 0,
 954                 i * 2 * n_extremities * n_factors  + heading_rows,
 955                 TAB_LEFT | TAT_TITLE,
 956                 var_to_string(dependent_var[i])
 957                 );
 958
 959
 960       if ( !fctr )
 961         populate_extremes(tbl, heading_columns - 2,
 962                           i * 2 * n_extremities * n_factors  + heading_rows,
 963                           n_extremities, &totals[i]);
 964
 965       else
 966         {
 967           struct factor_statistics **fs = fctr->fs;
 968           int count = 0 ;
 969
 970           while (*fs)
 971             {
 972               static union value prev ;
 973
 974               const int row = heading_rows + ( 2 * n_extremities )  *
 975                 ( ( i  * n_factors  ) +  count );
 976
 977
 978               if ( 0 != compare_values(&prev, &(*fs)->id[0],
 979                                        fctr->indep_var[0]->width))
 980                 {
 981
 982                   if ( count > 0 )
 983                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
 984
 985                   tab_text (tbl,
 986                             1, row,
 987                             TAB_LEFT | TAT_TITLE,
 988                             value_to_string(&(*fs)->id[0], fctr->indep_var[0])
 989                             );
 990                 }
 991
 992               prev = (*fs)->id[0];
 993
 994               if (fctr->indep_var[1] && count > 0 )
 995                 tab_hline(tbl, TAL_1, 2, n_cols - 1, row);
 996
 997               if ( fctr->indep_var[1])
 998                 tab_text (tbl, 2, row,
 999                           TAB_LEFT | TAT_TITLE,
1000                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
1001                           );
1002
1003               populate_extremes(tbl, heading_columns - 2,
1004                                 row, n_extremities,
1005                                 &(*fs)->m[i]);
1006
1007               count++ ;
1008               fs++;
1009             }
1010         }
1011     }
1012
1013   tab_submit(tbl);
1014 }
1015
1016
1017
1018 /* Fill in the extremities table */
1019 void
1020 populate_extremes(struct tab_table *t,
1021                   int col, int row, int n, const struct metrics *m)
1022 {
1023   int extremity;
1024   int idx=0;
1025
1026
1027   tab_text(t, col, row,
1028            TAB_RIGHT | TAT_TITLE ,
1029            _("Highest")
1030            );
1031
1032   tab_text(t, col, row + n ,
1033            TAB_RIGHT | TAT_TITLE ,
1034            _("Lowest")
1035            );
1036
1037
1038   tab_hline(t, TAL_1, col, col + 3, row + n );
1039
1040   for (extremity = 0; extremity < n ; ++extremity )
1041     {
1042       /* Highest */
1043       tab_float(t, col + 1, row + extremity,
1044                 TAB_RIGHT,
1045                 extremity + 1, 8, 0);
1046
1047
1048       /* Lowest */
1049       tab_float(t, col + 1, row + extremity + n,
1050                 TAB_RIGHT,
1051                 extremity + 1, 8, 0);
1052
1053     }
1054
1055
1056   /* Lowest */
1057   for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx )
1058     {
1059       int j;
1060       const struct weighted_value *wv = m->wvp[idx];
1061       struct case_node *cn = wv->case_nos;
1062
1063
1064       for (j = 0 ; j < wv->w ; ++j  )
1065         {
1066           if ( extremity + j >= n )
1067             break ;
1068
1069           tab_float(t, col + 3, row + extremity + j  + n,
1070                     TAB_RIGHT,
1071                     wv->v.f, 8, 2);
1072
1073           tab_float(t, col + 2, row + extremity + j  + n,
1074                     TAB_RIGHT,
1075                     cn->num, 8, 0);
1076
1077           if ( cn->next )
1078               cn = cn->next;
1079
1080         }
1081
1082       extremity +=  wv->w ;
1083     }
1084
1085
1086   /* Highest */
1087   for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx )
1088     {
1089       int j;
1090       const struct weighted_value *wv = m->wvp[idx];
1091       struct case_node *cn = wv->case_nos;
1092
1093       for (j = 0 ; j < wv->w ; ++j  )
1094         {
1095           if ( extremity + j >= n )
1096             break ;
1097
1098           tab_float(t, col + 3, row + extremity + j,
1099                     TAB_RIGHT,
1100                     wv->v.f, 8, 2);
1101
1102           tab_float(t, col + 2, row + extremity + j,
1103                     TAB_RIGHT,
1104                     cn->num, 8, 0);
1105
1106           if ( cn->next )
1107               cn = cn->next;
1108
1109         }
1110
1111       extremity +=  wv->w ;
1112     }
1113 }
1114
1115
1116 /* Show the descriptives table */
1117 void
1118 show_descriptives(struct variable **dependent_var,
1119                   int n_dep_var,
1120                   struct factor *fctr)
1121 {
1122   int i;
1123   int heading_columns ;
1124   int n_cols;
1125   const int n_stat_rows = 13;
1126
1127   const int heading_rows = 1;
1128
1129   struct tab_table *tbl;
1130
1131   int n_factors = 1;
1132   int n_rows ;
1133
1134   if ( fctr )
1135     {
1136       heading_columns = 4;
1137       n_factors = hsh_count(fctr->fstats);
1138
1139       n_rows = n_dep_var * n_stat_rows * n_factors;
1140
1141       if ( fctr->indep_var[1] )
1142           heading_columns = 5;
1143     }
1144   else
1145     {
1146       heading_columns = 3;
1147       n_rows = n_dep_var * n_stat_rows;
1148     }
1149
1150   n_rows += heading_rows;
1151
1152   n_cols = heading_columns + 2;
1153
1154
1155   tbl = tab_create (n_cols, n_rows, 0);
1156
1157   tab_headers (tbl, heading_columns + 1, 0, heading_rows, 0);
1158
1159   tab_dim (tbl, tab_natural_dimensions);
1160
1161   /* Outline the box and have no internal lines*/
1162   tab_box (tbl,
1163            TAL_2, TAL_2,
1164            -1, -1,
1165            0, 0,
1166            n_cols - 1, n_rows - 1);
1167
1168   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1169
1170   tab_vline (tbl, TAL_1, heading_columns, 0, n_rows - 1);
1171   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows - 1);
1172   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1173
1174   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Statistic"));
1175   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
1176
1177   tab_title (tbl, 0, _("Descriptives"));
1178
1179
1180   for ( i = 0 ; i < n_dep_var ; ++i )
1181     {
1182       const int row = heading_rows + i * n_stat_rows * n_factors ;
1183
1184       if ( i > 0 )
1185         tab_hline(tbl, TAL_1, 0, n_cols - 1, row );
1186
1187       tab_text (tbl, 0,
1188                 i * n_stat_rows * n_factors  + heading_rows,
1189                 TAB_LEFT | TAT_TITLE,
1190                 var_to_string(dependent_var[i])
1191                 );
1192
1193
1194       if ( fctr  )
1195         {
1196           struct factor_statistics **fs = fctr->fs;
1197           int count = 0;
1198
1199           tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1200                     var_to_string(fctr->indep_var[0]));
1201
1202
1203           if ( fctr->indep_var[1])
1204             tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1205                       var_to_string(fctr->indep_var[1]));
1206
1207           while( *fs )
1208             {
1209
1210               static union value prev ;
1211
1212               const int row = heading_rows + n_stat_rows  *
1213                 ( ( i  * n_factors  ) +  count );
1214
1215
1216               if ( 0 != compare_values(&prev, &(*fs)->id[0],
1217                                        fctr->indep_var[0]->width))
1218                 {
1219
1220                   if ( count > 0 )
1221                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1222
1223                   tab_text (tbl,
1224                             1, row,
1225                             TAB_LEFT | TAT_TITLE,
1226                             value_to_string(&(*fs)->id[0], fctr->indep_var[0])
1227                             );
1228                 }
1229
1230               prev = (*fs)->id[0];
1231
1232               if (fctr->indep_var[1] && count > 0 )
1233                 tab_hline(tbl, TAL_1, 2, n_cols - 1, row);
1234
1235               if ( fctr->indep_var[1])
1236                 tab_text (tbl, 2, row,
1237                           TAB_LEFT | TAT_TITLE,
1238                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
1239                           );
1240
1241               populate_descriptives(tbl, heading_columns - 2,
1242                                 row, &(*fs)->m[i]);
1243
1244               count++ ;
1245               fs++;
1246             }
1247
1248         }
1249
1250       else
1251         {
1252
1253           populate_descriptives(tbl, heading_columns - 2,
1254                                 i * n_stat_rows * n_factors  + heading_rows,
1255                                 &totals[i]);
1256         }
1257     }
1258
1259   tab_submit(tbl);
1260
1261 }
1262
1263
1264
1265
1266
1267
1268 /* Fill in the descriptives data */
1269 void
1270 populate_descriptives(struct tab_table *tbl, int col, int row,
1271                       const struct metrics *m)
1272 {
1273
1274   const double t = gsl_cdf_tdist_Qinv(1 - cmd.n_cinterval[0]/100.0/2.0, \
1275                                       m->n -1);
1276
1277
1278   tab_text (tbl, col,
1279             row,
1280             TAB_LEFT | TAT_TITLE,
1281             _("Mean"));
1282
1283   tab_float (tbl, col + 2,
1284              row,
1285              TAB_CENTER,
1286              m->mean,
1287              8,2);
1288
1289   tab_float (tbl, col + 3,
1290              row,
1291              TAB_CENTER,
1292              m->stderr,
1293              8,3);
1294
1295
1296   tab_text (tbl, col,
1297             row + 1,
1298             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1299             _("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1300
1301
1302   tab_text (tbl, col + 1,
1303             row  + 1,
1304             TAB_LEFT | TAT_TITLE,
1305             _("Lower Bound"));
1306
1307   tab_float (tbl, col + 2,
1308              row + 1,
1309              TAB_CENTER,
1310              m->mean - t * m->stderr,
1311              8,3);
1312
1313   tab_text (tbl, col + 1,
1314             row + 2,
1315             TAB_LEFT | TAT_TITLE,
1316             _("Upper Bound"));
1317
1318
1319   tab_float (tbl, col + 2,
1320              row + 2,
1321              TAB_CENTER,
1322              m->mean + t * m->stderr,
1323              8,3);
1324
1325   tab_text (tbl, col,
1326             row + 3,
1327             TAB_LEFT | TAT_TITLE,
1328             _("5% Trimmed Mean"));
1329
1330   tab_float (tbl, col + 2,
1331              row + 3,
1332              TAB_CENTER,
1333              m->trimmed_mean,
1334              8,2);
1335
1336   tab_text (tbl, col,
1337             row + 4,
1338             TAB_LEFT | TAT_TITLE,
1339             _("Median"));
1340
1341   tab_text (tbl, col,
1342             row + 5,
1343             TAB_LEFT | TAT_TITLE,
1344             _("Variance"));
1345
1346   tab_float (tbl, col + 2,
1347              row + 5,
1348              TAB_CENTER,
1349              m->var,
1350              8,3);
1351
1352
1353   tab_text (tbl, col,
1354             row + 6,
1355             TAB_LEFT | TAT_TITLE,
1356             _("Std. Deviation"));
1357
1358
1359   tab_float (tbl, col + 2,
1360              row + 6,
1361              TAB_CENTER,
1362              m->stddev,
1363              8,3);
1364
1365
1366   tab_text (tbl, col,
1367             row + 7,
1368             TAB_LEFT | TAT_TITLE,
1369             _("Minimum"));
1370
1371   tab_float (tbl, col + 2,
1372              row + 7,
1373              TAB_CENTER,
1374              m->min,
1375              8,3);
1376
1377   tab_text (tbl, col,
1378             row + 8,
1379             TAB_LEFT | TAT_TITLE,
1380             _("Maximum"));
1381
1382   tab_float (tbl, col + 2,
1383              row + 8,
1384              TAB_CENTER,
1385              m->max,
1386              8,3);
1387
1388
1389   tab_text (tbl, col,
1390             row + 9,
1391             TAB_LEFT | TAT_TITLE,
1392             _("Range"));
1393
1394
1395   tab_float (tbl, col + 2,
1396              row + 9,
1397              TAB_CENTER,
1398              m->max - m->min,
1399              8,3);
1400
1401   tab_text (tbl, col,
1402             row + 10,
1403             TAB_LEFT | TAT_TITLE,
1404             _("Interquartile Range"));
1405
1406   tab_text (tbl, col,
1407             row + 11,
1408             TAB_LEFT | TAT_TITLE,
1409             _("Skewness"));
1410
1411   tab_text (tbl, col,
1412             row + 12,
1413             TAB_LEFT | TAT_TITLE,
1414             _("Kurtosis"));
1415 }
1416
1417
1418
1419
1420
1421 /* Plot the normal and detrended normal plots for m
1422    Label the plots with factorname */
1423 void
1424 np_plot(const struct metrics *m, const char *factorname)
1425 {
1426   int i;
1427   double yfirst=0, ylast=0;
1428
1429   /* Normal Plot */
1430   struct chart np_chart;
1431
1432   /* Detrended Normal Plot */
1433   struct chart dnp_chart;
1434
1435   const struct weighted_value *wv = *(m->wvp);
1436
1437
1438   /* The slope and intercept of the ideal normal probability line */
1439   const double slope = 1.0 / m->stddev;
1440   const double intercept = - m->mean / m->stddev;
1441
1442   /* Cowardly refuse to plot an empty data set */
1443   if ( m->n_data == 0 )
1444     return ;
1445
1446   chart_initialise(&np_chart);
1447   chart_write_title(&np_chart, _("Normal Q-Q Plot of %s"), factorname);
1448   chart_write_xlabel(&np_chart, _("Observed Value"));
1449   chart_write_ylabel(&np_chart, _("Expected Normal"));
1450
1451   chart_initialise(&dnp_chart);
1452   chart_write_title(&dnp_chart, _("Detrended Normal Q-Q Plot of %s"),
1453                     factorname);
1454   chart_write_xlabel(&dnp_chart, _("Observed Value"));
1455   chart_write_ylabel(&dnp_chart, _("Dev from Normal"));
1456
1457   yfirst = gsl_cdf_ugaussian_Pinv (wv[0].rank / ( m->n + 1));
1458   ylast =  gsl_cdf_ugaussian_Pinv (wv[m->n_data-1].rank / ( m->n + 1));
1459
1460   {
1461     /* Need to make sure that both the scatter plot and the ideal fit into the
1462        plot */
1463     double x_lower = min(m->min, (yfirst - intercept) / slope) ;
1464     double x_upper = max(m->max, (ylast  - intercept) / slope) ;
1465     double slack = (x_upper - x_lower)  * 0.05 ;
1466
1467     chart_write_xscale(&np_chart, x_lower  - slack, x_upper + slack,
1468                        chart_rounded_tick((m->max - m->min) / 5.0));
1469
1470
1471     chart_write_xscale(&dnp_chart, m->min, m->max,
1472                        chart_rounded_tick((m->max - m->min) / 5.0));
1473
1474   }
1475
1476   chart_write_yscale(&np_chart, yfirst, ylast,
1477                      chart_rounded_tick((ylast - yfirst)/5.0) );
1478
1479   {
1480   /* We have to cache the detrended data, beacause we need to
1481      find its limits before we can plot it */
1482   double *d_data;
1483   d_data = xmalloc (m->n_data * sizeof(double));
1484   double d_max = -DBL_MAX;
1485   double d_min = DBL_MAX;
1486   for ( i = 0 ; i < m->n_data; ++i )
1487     {
1488       const double ns = gsl_cdf_ugaussian_Pinv (wv[i].rank / ( m->n + 1));
1489
1490       chart_datum(&np_chart, 0, wv[i].v.f, ns);
1491
1492       d_data[i] = (wv[i].v.f - m->mean) / m->stddev  - ns;
1493
1494       if ( d_data[i] < d_min ) d_min = d_data[i];
1495       if ( d_data[i] > d_max ) d_max = d_data[i];
1496     }
1497
1498   chart_write_yscale(&dnp_chart, d_min, d_max,
1499                      chart_rounded_tick((d_max - d_min) / 5.0));
1500
1501   for ( i = 0 ; i < m->n_data; ++i )
1502       chart_datum(&dnp_chart, 0, wv[i].v.f, d_data[i]);
1503
1504   free(d_data);
1505   }
1506
1507   chart_line(&np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1508   chart_line(&dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1509
1510   chart_finalise(&np_chart);
1511   chart_finalise(&dnp_chart);
1512
1513 }