1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_cdf.h>
20 #include <gsl/gsl_matrix.h>
23 #include "data/case.h"
24 #include "data/casegrouper.h"
25 #include "data/casereader.h"
26 #include "data/dataset.h"
27 #include "data/dictionary.h"
28 #include "data/format.h"
29 #include "data/value.h"
30 #include "language/command.h"
31 #include "language/dictionary/split-file.h"
32 #include "language/lexer/lexer.h"
33 #include "language/lexer/value-parser.h"
34 #include "language/lexer/variable-parser.h"
35 #include "libpspp/ll.h"
36 #include "libpspp/message.h"
37 #include "libpspp/misc.h"
38 #include "libpspp/taint.h"
39 #include "linreg/sweep.h"
40 #include "math/categoricals.h"
41 #include "math/covariance.h"
42 #include "math/levene.h"
43 #include "math/moments.h"
44 #include "output/tab.h"
47 #define _(msgid) gettext (msgid)
58 STATS_DESCRIPTIVES = 0x0001,
59 STATS_HOMOGENEITY = 0x0002
72 struct ll_list coefficient_list;
74 bool bad_count; /* True if the number of coefficients does not equal the number of groups */
80 const struct variable **vars;
82 const struct variable *indep_var;
84 enum statistics stats;
86 enum missing_type missing_type;
87 enum mv_class exclude;
89 /* List of contrasts */
90 struct ll_list contrast_list;
92 /* The weight variable */
93 const struct variable *wv;
96 /* Per category data */
97 struct descriptive_data
99 const struct variable *var;
100 struct moments1 *mom;
106 /* Workspace variable for each dependent variable */
109 struct categoricals *cat;
110 struct covariance *cov;
122 struct oneway_workspace
124 /* The number of distinct values of the independent variable, when all
125 missing values are disregarded */
126 int actual_number_of_groups;
128 struct per_var_ws *vws;
130 /* An array of descriptive data. One for each dependent variable */
131 struct descriptive_data **dd_total;
134 /* Routines to show the output tables */
135 static void show_anova_table (const struct oneway_spec *, const struct oneway_workspace *);
136 static void show_descriptives (const struct oneway_spec *, const struct oneway_workspace *);
137 static void show_homogeneity (const struct oneway_spec *, const struct oneway_workspace *);
139 static void output_oneway (const struct oneway_spec *, struct oneway_workspace *ws);
140 static void run_oneway (const struct oneway_spec *cmd, struct casereader *input, const struct dataset *ds);
143 cmd_oneway (struct lexer *lexer, struct dataset *ds)
145 const struct dictionary *dict = dataset_dict (ds);
146 struct oneway_spec oneway ;
149 oneway.indep_var = NULL;
151 oneway.missing_type = MISS_ANALYSIS;
152 oneway.exclude = MV_ANY;
153 oneway.wv = dict_get_weight (dict);
155 ll_init (&oneway.contrast_list);
158 if ( lex_match (lexer, T_SLASH))
160 if (!lex_force_match_id (lexer, "VARIABLES"))
164 lex_match (lexer, T_EQUALS);
167 if (!parse_variables_const (lexer, dict,
168 &oneway.vars, &oneway.n_vars,
169 PV_NO_DUPLICATE | PV_NUMERIC))
172 lex_force_match (lexer, T_BY);
174 oneway.indep_var = parse_variable_const (lexer, dict);
176 while (lex_token (lexer) != T_ENDCMD)
178 lex_match (lexer, T_SLASH);
180 if (lex_match_id (lexer, "STATISTICS"))
182 lex_match (lexer, T_EQUALS);
183 while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
185 if (lex_match_id (lexer, "DESCRIPTIVES"))
187 oneway.stats |= STATS_DESCRIPTIVES;
189 else if (lex_match_id (lexer, "HOMOGENEITY"))
191 oneway.stats |= STATS_HOMOGENEITY;
195 lex_error (lexer, NULL);
200 else if (lex_match_id (lexer, "CONTRAST"))
202 struct contrasts_node *cl = xzalloc (sizeof *cl);
204 struct ll_list *coefficient_list = &cl->coefficient_list;
205 lex_match (lexer, T_EQUALS);
207 ll_init (coefficient_list);
209 while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
211 if ( lex_is_number (lexer))
213 struct coeff_node *cc = xmalloc (sizeof *cc);
214 cc->coeff = lex_number (lexer);
216 ll_push_tail (coefficient_list, &cc->ll);
221 lex_error (lexer, NULL);
226 ll_push_tail (&oneway.contrast_list, &cl->ll);
228 else if (lex_match_id (lexer, "MISSING"))
230 lex_match (lexer, T_EQUALS);
231 while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
233 if (lex_match_id (lexer, "INCLUDE"))
235 oneway.exclude = MV_SYSTEM;
237 else if (lex_match_id (lexer, "EXCLUDE"))
239 oneway.exclude = MV_ANY;
241 else if (lex_match_id (lexer, "LISTWISE"))
243 oneway.missing_type = MISS_LISTWISE;
245 else if (lex_match_id (lexer, "ANALYSIS"))
247 oneway.missing_type = MISS_ANALYSIS;
251 lex_error (lexer, NULL);
258 lex_error (lexer, NULL);
265 struct casegrouper *grouper;
266 struct casereader *group;
269 grouper = casegrouper_create_splits (proc_open (ds), dict);
270 while (casegrouper_get_next_group (grouper, &group))
271 run_oneway (&oneway, group, ds);
272 ok = casegrouper_destroy (grouper);
273 ok = proc_commit (ds) && ok;
288 static struct descriptive_data *
289 dd_create (const struct variable *var)
291 struct descriptive_data *dd = xmalloc (sizeof *dd);
293 dd->mom = moments1_create (MOMENT_VARIANCE);
294 dd->minimum = DBL_MAX;
295 dd->maximum = -DBL_MAX;
302 dd_destroy (struct descriptive_data *dd)
304 moments1_destroy (dd->mom);
309 makeit (void *aux1, void *aux2 UNUSED)
311 const struct variable *var = aux1;
313 struct descriptive_data *dd = dd_create (var);
319 updateit (void *user_data,
320 enum mv_class exclude,
321 const struct variable *wv,
322 const struct variable *catvar UNUSED,
323 const struct ccase *c,
324 void *aux1, void *aux2)
326 struct descriptive_data *dd = user_data;
328 const struct variable *varp = aux1;
330 const union value *valx = case_data (c, varp);
332 struct descriptive_data *dd_total = aux2;
336 if ( var_is_value_missing (varp, valx, exclude))
339 weight = wv != NULL ? case_data (c, wv)->f : 1.0;
341 moments1_add (dd->mom, valx->f, weight);
342 if (valx->f < dd->minimum)
343 dd->minimum = valx->f;
345 if (valx->f > dd->maximum)
346 dd->maximum = valx->f;
349 const struct variable *var = dd_total->var;
350 const union value *val = case_data (c, var);
352 moments1_add (dd_total->mom,
356 if (val->f < dd_total->minimum)
357 dd_total->minimum = val->f;
359 if (val->f > dd_total->maximum)
360 dd_total->maximum = val->f;
365 run_oneway (const struct oneway_spec *cmd,
366 struct casereader *input,
367 const struct dataset *ds)
371 struct dictionary *dict = dataset_dict (ds);
372 struct casereader *reader;
375 struct oneway_workspace ws;
377 ws.actual_number_of_groups = 0;
378 ws.vws = xzalloc (cmd->n_vars * sizeof (*ws.vws));
379 ws.dd_total = xmalloc (sizeof (struct descriptive_data) * cmd->n_vars);
381 for (v = 0 ; v < cmd->n_vars; ++v)
382 ws.dd_total[v] = dd_create (cmd->vars[v]);
384 for (v = 0; v < cmd->n_vars; ++v)
386 ws.vws[v].cat = categoricals_create (&cmd->indep_var, 1, cmd->wv,
387 cmd->exclude, makeit, updateit,
388 CONST_CAST (struct variable *,
392 ws.vws[v].cov = covariance_2pass_create (1, &cmd->vars[v],
394 cmd->wv, cmd->exclude);
395 ws.vws[v].nl = levene_create (var_get_width (cmd->indep_var), NULL);
398 c = casereader_peek (input, 0);
401 casereader_destroy (input);
404 output_split_file_values (ds, c);
407 taint = taint_clone (casereader_get_taint (input));
409 input = casereader_create_filter_missing (input, &cmd->indep_var, 1,
410 cmd->exclude, NULL, NULL);
411 if (cmd->missing_type == MISS_LISTWISE)
412 input = casereader_create_filter_missing (input, cmd->vars, cmd->n_vars,
413 cmd->exclude, NULL, NULL);
414 input = casereader_create_filter_weight (input, dict, NULL, NULL);
416 reader = casereader_clone (input);
417 for (; (c = casereader_read (reader)) != NULL; case_unref (c))
420 double w = dict_get_case_weight (dict, c, NULL);
422 for (i = 0; i < cmd->n_vars; ++i)
424 struct per_var_ws *pvw = &ws.vws[i];
425 const struct variable *v = cmd->vars[i];
426 const union value *val = case_data (c, v);
428 if ( MISS_ANALYSIS == cmd->missing_type)
430 if ( var_is_value_missing (v, val, cmd->exclude))
434 covariance_accumulate_pass1 (pvw->cov, c);
435 levene_pass_one (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
438 casereader_destroy (reader);
440 reader = casereader_clone (input);
441 for ( ; (c = casereader_read (reader) ); case_unref (c))
444 double w = dict_get_case_weight (dict, c, NULL);
445 for (i = 0; i < cmd->n_vars; ++i)
447 struct per_var_ws *pvw = &ws.vws[i];
448 const struct variable *v = cmd->vars[i];
449 const union value *val = case_data (c, v);
451 if ( MISS_ANALYSIS == cmd->missing_type)
453 if ( var_is_value_missing (v, val, cmd->exclude))
457 covariance_accumulate_pass2 (pvw->cov, c);
458 levene_pass_two (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
461 casereader_destroy (reader);
463 reader = casereader_clone (input);
464 for ( ; (c = casereader_read (reader) ); case_unref (c))
467 double w = dict_get_case_weight (dict, c, NULL);
469 for (i = 0; i < cmd->n_vars; ++i)
471 struct per_var_ws *pvw = &ws.vws[i];
472 const struct variable *v = cmd->vars[i];
473 const union value *val = case_data (c, v);
475 if ( MISS_ANALYSIS == cmd->missing_type)
477 if ( var_is_value_missing (v, val, cmd->exclude))
481 levene_pass_three (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
484 casereader_destroy (reader);
487 for (v = 0; v < cmd->n_vars; ++v)
489 struct per_var_ws *pvw = &ws.vws[v];
490 gsl_matrix *cm = covariance_calculate_unnormalized (pvw->cov);
491 const struct categoricals *cats = covariance_get_categoricals (pvw->cov);
494 moments1_calculate (ws.dd_total[v]->mom, &n, NULL, NULL, NULL, NULL);
496 pvw->sst = gsl_matrix_get (cm, 0, 0);
500 pvw->sse = gsl_matrix_get (cm, 0, 0);
502 pvw->ssa = pvw->sst - pvw->sse;
504 pvw->n_groups = categoricals_total (cats);
506 pvw->mse = (pvw->sst - pvw->ssa) / (n - pvw->n_groups);
508 gsl_matrix_free (cm);
511 for (v = 0; v < cmd->n_vars; ++v)
513 const struct categoricals *cats = covariance_get_categoricals (ws.vws[v].cov);
515 categoricals_done (cats);
517 if (categoricals_total (cats) > ws.actual_number_of_groups)
518 ws.actual_number_of_groups = categoricals_total (cats);
521 casereader_destroy (input);
523 if (!taint_has_tainted_successor (taint))
524 output_oneway (cmd, &ws);
526 taint_destroy (taint);
529 for (v = 0; v < cmd->n_vars; ++v)
531 covariance_destroy (ws.vws[v].cov);
532 levene_destroy (ws.vws[v].nl);
533 dd_destroy (ws.dd_total[v]);
539 static void show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspace *ws);
540 static void show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspace *ws);
543 output_oneway (const struct oneway_spec *cmd, struct oneway_workspace *ws)
547 /* Check the sanity of the given contrast values */
548 struct contrasts_node *coeff_list = NULL;
549 ll_for_each (coeff_list, struct contrasts_node, ll, &cmd->contrast_list)
551 struct coeff_node *cn = NULL;
553 struct ll_list *cl = &coeff_list->coefficient_list;
556 if (ll_count (cl) != ws->actual_number_of_groups)
559 _("Number of contrast coefficients must equal the number of groups"));
560 coeff_list->bad_count = true;
564 ll_for_each (cn, struct coeff_node, ll, cl)
568 msg (SW, _("Coefficients for contrast %zu do not total zero"), i);
571 if (cmd->stats & STATS_DESCRIPTIVES)
572 show_descriptives (cmd, ws);
574 if (cmd->stats & STATS_HOMOGENEITY)
575 show_homogeneity (cmd, ws);
577 show_anova_table (cmd, ws);
579 if (ll_count (&cmd->contrast_list) > 0)
581 show_contrast_coeffs (cmd, ws);
582 show_contrast_tests (cmd, ws);
587 /* Show the ANOVA table */
589 show_anova_table (const struct oneway_spec *cmd, const struct oneway_workspace *ws)
593 size_t n_rows = cmd->n_vars * 3 + 1;
595 struct tab_table *t = tab_create (n_cols, n_rows);
597 tab_headers (t, 2, 0, 1, 0);
603 n_cols - 1, n_rows - 1);
605 tab_hline (t, TAL_2, 0, n_cols - 1, 1 );
606 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
607 tab_vline (t, TAL_0, 1, 0, 0);
609 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
610 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
611 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
612 tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
613 tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
616 for (i = 0; i < cmd->n_vars; ++i)
621 const char *s = var_to_string (cmd->vars[i]);
622 const struct per_var_ws *pvw = &ws->vws[i];
624 moments1_calculate (ws->dd_total[i]->mom, &n, NULL, NULL, NULL, NULL);
626 df1 = pvw->n_groups - 1;
627 df2 = n - pvw->n_groups;
628 msa = pvw->ssa / df1;
630 tab_text (t, 0, i * 3 + 1, TAB_LEFT | TAT_TITLE, s);
631 tab_text (t, 1, i * 3 + 1, TAB_LEFT | TAT_TITLE, _("Between Groups"));
632 tab_text (t, 1, i * 3 + 2, TAB_LEFT | TAT_TITLE, _("Within Groups"));
633 tab_text (t, 1, i * 3 + 3, TAB_LEFT | TAT_TITLE, _("Total"));
636 tab_hline (t, TAL_1, 0, n_cols - 1, i * 3 + 1);
639 /* Sums of Squares */
640 tab_double (t, 2, i * 3 + 1, 0, pvw->ssa, NULL);
641 tab_double (t, 2, i * 3 + 3, 0, pvw->sst, NULL);
642 tab_double (t, 2, i * 3 + 2, 0, pvw->sse, NULL);
645 /* Degrees of freedom */
646 tab_fixed (t, 3, i * 3 + 1, 0, df1, 4, 0);
647 tab_fixed (t, 3, i * 3 + 2, 0, df2, 4, 0);
648 tab_fixed (t, 3, i * 3 + 3, 0, n - 1, 4, 0);
651 tab_double (t, 4, i * 3 + 1, TAB_RIGHT, msa, NULL);
652 tab_double (t, 4, i * 3 + 2, TAB_RIGHT, pvw->mse, NULL);
655 const double F = msa / pvw->mse ;
658 tab_double (t, 5, i * 3 + 1, 0, F, NULL);
660 /* The significance */
661 tab_double (t, 6, i * 3 + 1, 0, gsl_cdf_fdist_Q (F, df1, df2), NULL);
665 tab_title (t, _("ANOVA"));
670 /* Show the descriptives table */
672 show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace *ws)
679 const double confidence = 0.95;
680 const double q = (1.0 - confidence) / 2.0;
682 const struct fmt_spec *wfmt = cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
686 for (v = 0; v < cmd->n_vars; ++v)
687 n_rows += ws->actual_number_of_groups + 1;
689 t = tab_create (n_cols, n_rows);
690 tab_headers (t, 2, 0, 2, 0);
692 /* Put a frame around the entire box, and vertical lines inside */
697 n_cols - 1, n_rows - 1);
699 /* Underline headers */
700 tab_hline (t, TAL_2, 0, n_cols - 1, 2);
701 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
703 tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("N"));
704 tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Mean"));
705 tab_text (t, 4, 1, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
706 tab_text (t, 5, 1, TAB_CENTER | TAT_TITLE, _("Std. Error"));
709 tab_vline (t, TAL_0, 7, 0, 0);
710 tab_hline (t, TAL_1, 6, 7, 1);
711 tab_joint_text_format (t, 6, 0, 7, 0, TAB_CENTER | TAT_TITLE,
712 _("%g%% Confidence Interval for Mean"),
715 tab_text (t, 6, 1, TAB_CENTER | TAT_TITLE, _("Lower Bound"));
716 tab_text (t, 7, 1, TAB_CENTER | TAT_TITLE, _("Upper Bound"));
718 tab_text (t, 8, 1, TAB_CENTER | TAT_TITLE, _("Minimum"));
719 tab_text (t, 9, 1, TAB_CENTER | TAT_TITLE, _("Maximum"));
721 tab_title (t, _("Descriptives"));
724 for (v = 0; v < cmd->n_vars; ++v)
726 const char *s = var_to_string (cmd->vars[v]);
727 const struct fmt_spec *fmt = var_get_print_format (cmd->vars[v]);
731 struct per_var_ws *pvw = &ws->vws[v];
732 const struct categoricals *cats = covariance_get_categoricals (pvw->cov);
734 tab_text (t, 0, row, TAB_LEFT | TAT_TITLE, s);
736 tab_hline (t, TAL_1, 0, n_cols - 1, row);
738 for (count = 0; count < categoricals_total (cats); ++count)
741 double n, mean, variance;
742 double std_dev, std_error ;
746 const union value *gval = categoricals_get_value_by_category (cats, count);
747 const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, count);
749 moments1_calculate (dd->mom, &n, &mean, &variance, NULL, NULL);
751 std_dev = sqrt (variance);
752 std_error = std_dev / sqrt (n) ;
754 ds_init_empty (&vstr);
756 var_append_value_name (cmd->indep_var, gval, &vstr);
758 tab_text (t, 1, row + count,
759 TAB_LEFT | TAT_TITLE,
764 /* Now fill in the numbers ... */
766 tab_double (t, 2, row + count, 0, n, wfmt);
768 tab_double (t, 3, row + count, 0, mean, NULL);
770 tab_double (t, 4, row + count, 0, std_dev, NULL);
773 tab_double (t, 5, row + count, 0, std_error, NULL);
775 /* Now the confidence interval */
777 T = gsl_cdf_tdist_Qinv (q, n - 1);
779 tab_double (t, 6, row + count, 0,
780 mean - T * std_error, NULL);
782 tab_double (t, 7, row + count, 0,
783 mean + T * std_error, NULL);
787 tab_double (t, 8, row + count, 0, dd->minimum, fmt);
788 tab_double (t, 9, row + count, 0, dd->maximum, fmt);
793 double n, mean, variance;
797 moments1_calculate (ws->dd_total[v]->mom, &n, &mean, &variance, NULL, NULL);
799 std_dev = sqrt (variance);
800 std_error = std_dev / sqrt (n) ;
802 tab_text (t, 1, row + count,
803 TAB_LEFT | TAT_TITLE, _("Total"));
805 tab_double (t, 2, row + count, 0, n, wfmt);
807 tab_double (t, 3, row + count, 0, mean, NULL);
809 tab_double (t, 4, row + count, 0, std_dev, NULL);
811 tab_double (t, 5, row + count, 0, std_error, NULL);
813 /* Now the confidence interval */
814 T = gsl_cdf_tdist_Qinv (q, n - 1);
816 tab_double (t, 6, row + count, 0,
817 mean - T * std_error, NULL);
819 tab_double (t, 7, row + count, 0,
820 mean + T * std_error, NULL);
823 tab_double (t, 8, row + count, 0, ws->dd_total[v]->minimum, fmt);
824 tab_double (t, 9, row + count, 0, ws->dd_total[v]->maximum, fmt);
827 row += categoricals_total (cats) + 1;
833 /* Show the homogeneity table */
835 show_homogeneity (const struct oneway_spec *cmd, const struct oneway_workspace *ws)
839 size_t n_rows = cmd->n_vars + 1;
841 struct tab_table *t = tab_create (n_cols, n_rows);
842 tab_headers (t, 1, 0, 1, 0);
844 /* Put a frame around the entire box, and vertical lines inside */
849 n_cols - 1, n_rows - 1);
852 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
853 tab_vline (t, TAL_2, 1, 0, n_rows - 1);
855 tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Levene Statistic"));
856 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df1"));
857 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df2"));
858 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
860 tab_title (t, _("Test of Homogeneity of Variances"));
862 for (v = 0; v < cmd->n_vars; ++v)
865 const struct per_var_ws *pvw = &ws->vws[v];
866 double F = levene_calculate (pvw->nl);
868 const struct variable *var = cmd->vars[v];
869 const char *s = var_to_string (var);
872 moments1_calculate (ws->dd_total[v]->mom, &n, NULL, NULL, NULL, NULL);
874 df1 = pvw->n_groups - 1;
875 df2 = n - pvw->n_groups;
877 tab_text (t, 0, v + 1, TAB_LEFT | TAT_TITLE, s);
879 tab_double (t, 1, v + 1, TAB_RIGHT, F, NULL);
880 tab_fixed (t, 2, v + 1, TAB_RIGHT, df1, 8, 0);
881 tab_fixed (t, 3, v + 1, TAB_RIGHT, df2, 8, 0);
883 /* Now the significance */
884 tab_double (t, 4, v + 1, TAB_RIGHT, gsl_cdf_fdist_Q (F, df1, df2), NULL);
891 /* Show the contrast coefficients table */
893 show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspace *ws)
898 int n_contrasts = ll_count (&cmd->contrast_list);
899 int n_cols = 2 + ws->actual_number_of_groups;
900 int n_rows = 2 + n_contrasts;
904 const struct covariance *cov = ws->vws[0].cov ;
906 t = tab_create (n_cols, n_rows);
907 tab_headers (t, 2, 0, 2, 0);
909 /* Put a frame around the entire box, and vertical lines inside */
914 n_cols - 1, n_rows - 1);
928 tab_hline (t, TAL_1, 2, n_cols - 1, 1);
929 tab_hline (t, TAL_2, 0, n_cols - 1, 2);
931 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
933 tab_title (t, _("Contrast Coefficients"));
935 tab_text (t, 0, 2, TAB_LEFT | TAT_TITLE, _("Contrast"));
938 tab_joint_text (t, 2, 0, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
939 var_to_string (cmd->indep_var));
941 for ( cli = ll_head (&cmd->contrast_list);
942 cli != ll_null (&cmd->contrast_list);
946 struct contrasts_node *cn = ll_data (cli, struct contrasts_node, ll);
949 tab_text_format (t, 1, c_num + 2, TAB_CENTER, "%d", c_num + 1);
951 for (coeffi = ll_head (&cn->coefficient_list);
952 coeffi != ll_null (&cn->coefficient_list);
953 ++count, coeffi = ll_next (coeffi))
955 const struct categoricals *cats = covariance_get_categoricals (cov);
956 const union value *val = categoricals_get_value_by_category (cats, count);
959 ds_init_empty (&vstr);
961 var_append_value_name (cmd->indep_var, val, &vstr);
963 tab_text (t, count + 2, 1, TAB_CENTER | TAT_TITLE, ds_cstr (&vstr));
968 tab_text (t, count + 2, c_num + 2, TAB_RIGHT, "?" );
971 struct coeff_node *coeffn = ll_data (coeffi, struct coeff_node, ll);
973 tab_text_format (t, count + 2, c_num + 2, TAB_RIGHT, "%g", coeffn->coeff);
983 /* Show the results of the contrast tests */
985 show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspace *ws)
987 int n_contrasts = ll_count (&cmd->contrast_list);
990 size_t n_rows = 1 + cmd->n_vars * 2 * n_contrasts;
994 t = tab_create (n_cols, n_rows);
995 tab_headers (t, 3, 0, 1, 0);
997 /* Put a frame around the entire box, and vertical lines inside */
1002 n_cols - 1, n_rows - 1);
1010 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
1011 tab_vline (t, TAL_2, 3, 0, n_rows - 1);
1013 tab_title (t, _("Contrast Tests"));
1015 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Contrast"));
1016 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Value of Contrast"));
1017 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
1018 tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
1019 tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("df"));
1020 tab_text (t, 7, 0, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1022 for (v = 0; v < cmd->n_vars; ++v)
1024 const struct per_var_ws *pvw = &ws->vws[v];
1025 const struct categoricals *cats = covariance_get_categoricals (pvw->cov);
1028 int lines_per_variable = 2 * n_contrasts;
1030 tab_text (t, 0, (v * lines_per_variable) + 1, TAB_LEFT | TAT_TITLE,
1031 var_to_string (cmd->vars[v]));
1033 for ( cli = ll_head (&cmd->contrast_list);
1034 cli != ll_null (&cmd->contrast_list);
1035 ++i, cli = ll_next (cli))
1037 struct contrasts_node *cn = ll_data (cli, struct contrasts_node, ll);
1040 double contrast_value = 0.0;
1041 double coef_msq = 0.0;
1044 double std_error_contrast;
1046 double sec_vneq = 0.0;
1048 /* Note: The calculation of the degrees of freedom in the
1049 "variances not equal" case is painfull!!
1050 The following formula may help to understand it:
1051 \frac{\left (\sum_{i=1}^k{c_i^2\frac{s_i^2}{n_i}}\right)^2}
1054 \frac{\left (c_i^2\frac{s_i^2}{n_i}\right)^2} {n_i-1}
1059 double df_denominator = 0.0;
1060 double df_numerator = 0.0;
1063 moments1_calculate (ws->dd_total[v]->mom, &grand_n, NULL, NULL, NULL, NULL);
1064 df = grand_n - pvw->n_groups;
1068 tab_text (t, 1, (v * lines_per_variable) + i + 1,
1069 TAB_LEFT | TAT_TITLE,
1070 _("Assume equal variances"));
1072 tab_text (t, 1, (v * lines_per_variable) + i + 1 + n_contrasts,
1073 TAB_LEFT | TAT_TITLE,
1074 _("Does not assume equal"));
1077 tab_text_format (t, 2, (v * lines_per_variable) + i + 1,
1078 TAB_CENTER | TAT_TITLE, "%d", i + 1);
1081 tab_text_format (t, 2,
1082 (v * lines_per_variable) + i + 1 + n_contrasts,
1083 TAB_CENTER | TAT_TITLE, "%d", i + 1);
1088 for (coeffi = ll_head (&cn->coefficient_list);
1089 coeffi != ll_null (&cn->coefficient_list);
1090 ++ci, coeffi = ll_next (coeffi))
1092 double n, mean, variance;
1093 const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, ci);
1094 struct coeff_node *cn = ll_data (coeffi, struct coeff_node, ll);
1095 const double coef = cn->coeff;
1098 moments1_calculate (dd->mom, &n, &mean, &variance, NULL, NULL);
1100 winv = variance / n;
1102 contrast_value += coef * mean;
1104 coef_msq += (pow2 (coef)) / n;
1106 sec_vneq += (pow2 (coef)) * variance / n;
1108 df_numerator += (pow2 (coef)) * winv;
1109 df_denominator += pow2((pow2 (coef)) * winv) / (n - 1);
1112 sec_vneq = sqrt (sec_vneq);
1114 df_numerator = pow2 (df_numerator);
1116 tab_double (t, 3, (v * lines_per_variable) + i + 1,
1117 TAB_RIGHT, contrast_value, NULL);
1119 tab_double (t, 3, (v * lines_per_variable) + i + 1 +
1121 TAB_RIGHT, contrast_value, NULL);
1123 std_error_contrast = sqrt (pvw->mse * coef_msq);
1126 tab_double (t, 4, (v * lines_per_variable) + i + 1,
1127 TAB_RIGHT, std_error_contrast,
1130 T = fabs (contrast_value / std_error_contrast);
1134 tab_double (t, 5, (v * lines_per_variable) + i + 1,
1139 /* Degrees of Freedom */
1140 tab_fixed (t, 6, (v * lines_per_variable) + i + 1,
1145 /* Significance TWO TAILED !!*/
1146 tab_double (t, 7, (v * lines_per_variable) + i + 1,
1147 TAB_RIGHT, 2 * gsl_cdf_tdist_Q (T, df),
1150 /* Now for the Variances NOT Equal case */
1154 (v * lines_per_variable) + i + 1 + n_contrasts,
1155 TAB_RIGHT, sec_vneq,
1158 T = contrast_value / sec_vneq;
1160 (v * lines_per_variable) + i + 1 + n_contrasts,
1164 df = df_numerator / df_denominator;
1167 (v * lines_per_variable) + i + 1 + n_contrasts,
1171 /* The Significance */
1172 tab_double (t, 7, (v * lines_per_variable) + i + 1 + n_contrasts,
1173 TAB_RIGHT, 2 * gsl_cdf_tdist_Q (T,df),
1178 tab_hline (t, TAL_1, 0, n_cols - 1, (v * lines_per_variable) + 1);