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