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