1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2005, 2009, 2010, 2011, 2012, 2013, 2014,
3 2016, 2017, 2019 Free Software Foundation, Inc.
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>. */
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_cdf.h>
25 #include <gsl/gsl_matrix.h>
27 #include <data/dataset.h>
28 #include <data/casewriter.h>
30 #include "language/command.h"
31 #include "language/lexer/lexer.h"
32 #include "language/lexer/value-parser.h"
33 #include "language/lexer/variable-parser.h"
36 #include "data/casegrouper.h"
37 #include "data/casereader.h"
38 #include "data/dictionary.h"
40 #include "math/covariance.h"
41 #include "math/linreg.h"
42 #include "math/moments.h"
44 #include "libpspp/message.h"
45 #include "libpspp/taint.h"
47 #include "output/pivot-table.h"
49 #include "gl/intprops.h"
50 #include "gl/minmax.h"
53 #define _(msgid) gettext (msgid)
54 #define N_(msgid) msgid
65 #define STATS_DEFAULT (STATS_R | STATS_COEFF | STATS_ANOVA | STATS_OUTS)
73 const struct variable **vars;
76 const struct variable **dep_vars;
88 struct regression_workspace
90 /* The new variables which will be introduced by /SAVE */
91 const struct variable **predvars;
92 const struct variable **residvars;
94 /* A reader/writer pair to temporarily hold the
95 values of the new variables */
96 struct casewriter *writer;
97 struct casereader *reader;
99 /* Indeces of the new values in the reader/writer (-1 if not applicable) */
103 /* 0, 1 or 2 depending on what new variables are to be created */
107 static void run_regression (const struct regression *cmd,
108 struct regression_workspace *ws,
109 struct casereader *input);
112 /* Return a string based on PREFIX which may be used as the name
113 of a new variable in DICT */
115 reg_get_name (const struct dictionary *dict, const char *prefix)
120 /* XXX handle too-long prefixes */
121 name = xmalloc (strlen (prefix) + INT_BUFSIZE_BOUND (i) + 1);
124 sprintf (name, "%s%d", prefix, i);
125 if (dict_lookup_var (dict, name) == NULL)
131 static const struct variable *
132 create_aux_var (struct dataset *ds, const char *prefix)
134 struct variable *var;
135 struct dictionary *dict = dataset_dict (ds);
136 char *name = reg_get_name (dict, prefix);
137 var = dict_create_var_assert (dict, name, 0);
142 /* Auxiliary data for transformation when /SAVE is entered */
143 struct save_trans_data
146 struct regression_workspace *ws;
150 save_trans_free (void *aux)
152 struct save_trans_data *save_trans_data = aux;
153 free (save_trans_data->ws->predvars);
154 free (save_trans_data->ws->residvars);
156 casereader_destroy (save_trans_data->ws->reader);
157 free (save_trans_data->ws);
158 free (save_trans_data);
162 static enum trns_result
163 save_trans_func (void *aux, struct ccase **c, casenumber x UNUSED)
165 struct save_trans_data *save_trans_data = aux;
166 struct regression_workspace *ws = save_trans_data->ws;
167 struct ccase *in = casereader_read (ws->reader);
172 *c = case_unshare (*c);
174 for (k = 0; k < save_trans_data->n_dep_vars; ++k)
176 if (ws->pred_idx != -1)
178 double pred = case_num_idx (in, ws->extras * k + ws->pred_idx);
179 *case_num_rw (*c, ws->predvars[k]) = pred;
182 if (ws->res_idx != -1)
184 double resid = case_num_idx (in, ws->extras * k + ws->res_idx);
185 *case_num_rw (*c, ws->residvars[k]) = resid;
191 return TRNS_CONTINUE;
195 cmd_regression (struct lexer *lexer, struct dataset *ds)
197 struct regression_workspace workspace;
198 struct regression regression;
199 const struct dictionary *dict = dataset_dict (ds);
202 memset (®ression, 0, sizeof (struct regression));
204 regression.ci = 0.95;
205 regression.stats = STATS_DEFAULT;
206 regression.pred = false;
207 regression.resid = false;
210 regression.origin = false;
212 bool variables_seen = false;
213 bool method_seen = false;
214 bool dependent_seen = false;
217 while (lex_token (lexer) != T_ENDCMD)
219 lex_match (lexer, T_SLASH);
221 if (lex_match_id (lexer, "VARIABLES"))
225 lex_next_error (lexer, -1, -1,
226 _("VARIABLES may not appear after %s"), "METHOD");
231 lex_next_error (lexer, -1, -1,
232 _("VARIABLES may not appear after %s"), "DEPENDENT");
235 variables_seen = true;
236 lex_match (lexer, T_EQUALS);
238 if (!parse_variables_const (lexer, dict,
239 ®ression.vars, ®ression.n_vars,
240 PV_NO_DUPLICATE | PV_NUMERIC))
243 else if (lex_match_id (lexer, "DEPENDENT"))
245 dependent_seen = true;
246 lex_match (lexer, T_EQUALS);
248 free (regression.dep_vars);
249 regression.n_dep_vars = 0;
251 if (!parse_variables_const (lexer, dict,
252 ®ression.dep_vars,
253 ®ression.n_dep_vars,
254 PV_NO_DUPLICATE | PV_NUMERIC))
257 else if (lex_match_id (lexer, "ORIGIN"))
259 regression.origin = true;
261 else if (lex_match_id (lexer, "NOORIGIN"))
263 regression.origin = false;
265 else if (lex_match_id (lexer, "METHOD"))
268 lex_match (lexer, T_EQUALS);
270 if (!lex_force_match_id (lexer, "ENTER"))
275 if (! variables_seen)
277 if (!parse_variables_const (lexer, dict,
278 ®ression.vars, ®ression.n_vars,
279 PV_NO_DUPLICATE | PV_NUMERIC))
283 else if (lex_match_id (lexer, "STATISTICS"))
285 unsigned long statistics = 0;
286 lex_match (lexer, T_EQUALS);
288 while (lex_token (lexer) != T_ENDCMD
289 && lex_token (lexer) != T_SLASH)
291 if (lex_match (lexer, T_ALL))
295 else if (lex_match_id (lexer, "DEFAULTS"))
297 statistics |= STATS_DEFAULT;
299 else if (lex_match_id (lexer, "R"))
301 statistics |= STATS_R;
303 else if (lex_match_id (lexer, "COEFF"))
305 statistics |= STATS_COEFF;
307 else if (lex_match_id (lexer, "ANOVA"))
309 statistics |= STATS_ANOVA;
311 else if (lex_match_id (lexer, "BCOV"))
313 statistics |= STATS_BCOV;
315 else if (lex_match_id (lexer, "TOL"))
317 statistics |= STATS_TOL;
319 else if (lex_match_id (lexer, "CI"))
321 statistics |= STATS_CI;
323 if (lex_match (lexer, T_LPAREN) &&
324 lex_force_num (lexer))
326 regression.ci = lex_number (lexer) / 100.0;
328 if (! lex_force_match (lexer, T_RPAREN))
334 lex_error (lexer, NULL);
340 regression.stats = statistics;
343 else if (lex_match_id (lexer, "SAVE"))
345 save_start = lex_ofs (lexer) - 1;
346 lex_match (lexer, T_EQUALS);
348 while (lex_token (lexer) != T_ENDCMD
349 && lex_token (lexer) != T_SLASH)
351 if (lex_match_id (lexer, "PRED"))
353 regression.pred = true;
355 else if (lex_match_id (lexer, "RESID"))
357 regression.resid = true;
361 lex_error (lexer, NULL);
365 save_end = lex_ofs (lexer) - 1;
369 lex_error (lexer, NULL);
374 if (!regression.vars)
376 dict_get_vars (dict, ®ression.vars, ®ression.n_vars, 0);
379 save = regression.pred || regression.resid;
380 workspace.extras = 0;
381 workspace.res_idx = -1;
382 workspace.pred_idx = -1;
383 workspace.writer = NULL;
384 workspace.reader = NULL;
385 workspace.residvars = NULL;
386 workspace.predvars = NULL;
390 struct caseproto *proto = caseproto_create ();
392 if (regression.resid)
394 workspace.res_idx = workspace.extras ++;
395 workspace.residvars = xcalloc (regression.n_dep_vars, sizeof (*workspace.residvars));
397 for (i = 0; i < regression.n_dep_vars; ++i)
399 workspace.residvars[i] = create_aux_var (ds, "RES");
400 proto = caseproto_add_width (proto, 0);
406 workspace.pred_idx = workspace.extras ++;
407 workspace.predvars = xcalloc (regression.n_dep_vars, sizeof (*workspace.predvars));
409 for (i = 0; i < regression.n_dep_vars; ++i)
411 workspace.predvars[i] = create_aux_var (ds, "PRED");
412 proto = caseproto_add_width (proto, 0);
416 if (proc_make_temporary_transformations_permanent (ds))
417 lex_ofs_msg (lexer, SW, save_start, save_end,
418 _("REGRESSION with SAVE ignores TEMPORARY. "
419 "Temporary transformations will be made permanent."));
421 if (dict_get_filter (dict))
422 lex_ofs_msg (lexer, SW, save_start, save_end,
423 _("REGRESSION with SAVE ignores FILTER. "
424 "All cases will be processed."));
426 workspace.writer = autopaging_writer_create (proto);
427 caseproto_unref (proto);
432 struct casegrouper *grouper;
433 struct casereader *group;
436 grouper = casegrouper_create_splits (proc_open_filtering (ds, !save), dict);
439 while (casegrouper_get_next_group (grouper, &group))
441 run_regression (®ression,
446 ok = casegrouper_destroy (grouper);
447 ok = proc_commit (ds) && ok;
450 if (workspace.writer)
452 struct save_trans_data *save_trans_data = xmalloc (sizeof *save_trans_data);
453 struct casereader *r = casewriter_make_reader (workspace.writer);
454 workspace.writer = NULL;
455 workspace.reader = r;
456 save_trans_data->ws = xmalloc (sizeof (workspace));
457 memcpy (save_trans_data->ws, &workspace, sizeof (workspace));
458 save_trans_data->n_dep_vars = regression.n_dep_vars;
460 static const struct trns_class trns_class = {
461 .name = "REGRESSION",
462 .execute = save_trans_func,
463 .destroy = save_trans_free,
465 add_transformation (ds, &trns_class, save_trans_data);
469 free (regression.vars);
470 free (regression.dep_vars);
475 free (regression.vars);
476 free (regression.dep_vars);
480 /* Return the size of the union of dependent and independent variables */
482 get_n_all_vars (const struct regression *cmd)
484 size_t result = cmd->n_vars;
488 result += cmd->n_dep_vars;
489 for (i = 0; i < cmd->n_dep_vars; i++)
491 for (j = 0; j < cmd->n_vars; j++)
493 if (cmd->vars[j] == cmd->dep_vars[i])
502 /* Fill VARS with the union of dependent and independent variables */
504 fill_all_vars (const struct variable **vars, const struct regression *cmd)
508 for (i = 0; i < cmd->n_vars; i++)
510 vars[i] = cmd->vars[i];
513 for (i = 0; i < cmd->n_dep_vars; i++)
517 for (j = 0; j < cmd->n_vars; j++)
519 if (cmd->dep_vars[i] == cmd->vars[j])
527 vars[cmd->n_vars + x++] = cmd->dep_vars[i];
533 /* Fill the array VARS, with all the predictor variables from CMD, except
536 fill_predictor_x (const struct variable **vars, const struct variable *x, const struct regression *cmd)
541 for (i = 0; i < cmd->n_vars; i++)
543 if (cmd->vars[i] == x)
546 vars[n++] = cmd->vars[i];
551 Is variable k the dependent variable?
554 is_depvar (const struct regression *cmd, size_t k, const struct variable *v)
556 return v == cmd->vars[k];
560 /* Identify the explanatory variables in v_variables. Returns
561 the number of independent variables. */
563 identify_indep_vars (const struct regression *cmd,
564 const struct variable **indep_vars,
565 const struct variable *depvar)
567 int n_indep_vars = 0;
570 for (i = 0; i < cmd->n_vars; i++)
571 if (!is_depvar (cmd, i, depvar))
572 indep_vars[n_indep_vars++] = cmd->vars[i];
573 if ((n_indep_vars < 1) && is_depvar (cmd, 0, depvar))
576 There is only one independent variable, and it is the same
577 as the dependent variable. Print a warning and continue.
581 ("The dependent variable is equal to the independent variable. "
582 "The least squares line is therefore Y=X. "
583 "Standard errors and related statistics may be meaningless."));
585 indep_vars[0] = cmd->vars[0];
591 fill_covariance (gsl_matrix * cov, struct covariance *all_cov,
592 const struct variable **vars,
593 size_t n_vars, const struct variable *dep_var,
594 const struct variable **all_vars, size_t n_all_vars,
599 size_t dep_subscript = SIZE_MAX;
601 const gsl_matrix *ssizes;
602 const gsl_matrix *mean_matrix;
603 const gsl_matrix *ssize_matrix;
606 const gsl_matrix *cm = covariance_calculate_unnormalized (all_cov);
611 rows = xnmalloc (cov->size1 - 1, sizeof (*rows));
613 for (i = 0; i < n_all_vars; i++)
615 for (j = 0; j < n_vars; j++)
617 if (vars[j] == all_vars[i])
622 if (all_vars[i] == dep_var)
627 assert (dep_subscript != SIZE_MAX);
629 mean_matrix = covariance_moments (all_cov, MOMENT_MEAN);
630 ssize_matrix = covariance_moments (all_cov, MOMENT_NONE);
631 for (i = 0; i < cov->size1 - 1; i++)
633 means[i] = gsl_matrix_get (mean_matrix, rows[i], 0)
634 / gsl_matrix_get (ssize_matrix, rows[i], 0);
635 for (j = 0; j < cov->size2 - 1; j++)
637 gsl_matrix_set (cov, i, j, gsl_matrix_get (cm, rows[i], rows[j]));
638 gsl_matrix_set (cov, j, i, gsl_matrix_get (cm, rows[j], rows[i]));
641 means[cov->size1 - 1] = gsl_matrix_get (mean_matrix, dep_subscript, 0)
642 / gsl_matrix_get (ssize_matrix, dep_subscript, 0);
643 ssizes = covariance_moments (all_cov, MOMENT_NONE);
644 result = gsl_matrix_get (ssizes, dep_subscript, rows[0]);
645 for (i = 0; i < cov->size1 - 1; i++)
647 gsl_matrix_set (cov, i, cov->size1 - 1,
648 gsl_matrix_get (cm, rows[i], dep_subscript));
649 gsl_matrix_set (cov, cov->size1 - 1, i,
650 gsl_matrix_get (cm, rows[i], dep_subscript));
651 if (result > gsl_matrix_get (ssizes, rows[i], dep_subscript))
653 result = gsl_matrix_get (ssizes, rows[i], dep_subscript);
656 gsl_matrix_set (cov, cov->size1 - 1, cov->size1 - 1,
657 gsl_matrix_get (cm, dep_subscript, dep_subscript));
664 struct model_container
666 struct linreg **models;
670 STATISTICS subcommand output functions.
672 static void reg_stats_r (const struct linreg *, const struct variable *);
673 static void reg_stats_coeff (const struct regression *, const struct linreg *,
674 const struct model_container *, const gsl_matrix *,
675 const struct variable *);
676 static void reg_stats_anova (const struct linreg *, const struct variable *);
677 static void reg_stats_bcov (const struct linreg *, const struct variable *);
680 static struct linreg **
681 run_regression_get_models (const struct regression *cmd,
682 struct casereader *input,
686 struct model_container *model_container = XCALLOC (cmd->n_vars, struct model_container);
689 struct covariance *cov;
690 struct casereader *reader;
692 if (cmd->stats & STATS_TOL)
694 for (i = 0; i < cmd->n_vars; i++)
696 struct regression subreg;
697 subreg.origin = cmd->origin;
699 subreg.n_vars = cmd->n_vars - 1;
700 subreg.n_dep_vars = 1;
701 subreg.vars = xmalloc (sizeof (*subreg.vars) * cmd->n_vars - 1);
702 subreg.dep_vars = xmalloc (sizeof (*subreg.dep_vars));
703 fill_predictor_x (subreg.vars, cmd->vars[i], cmd);
704 subreg.dep_vars[0] = cmd->vars[i];
705 subreg.stats = STATS_R;
707 subreg.resid = false;
710 model_container[i].models =
711 run_regression_get_models (&subreg, input, false);
713 free (subreg.dep_vars);
717 size_t n_all_vars = get_n_all_vars (cmd);
718 const struct variable **all_vars = xnmalloc (n_all_vars, sizeof (*all_vars));
720 /* In the (rather pointless) case where the dependent variable is
721 the independent variable, n_all_vars == 1.
722 However this would result in a buffer overflow so we must
723 over-allocate the space required in this malloc call.
725 double *means = xnmalloc (n_all_vars <= 1 ? 2 : n_all_vars,
727 fill_all_vars (all_vars, cmd);
728 cov = covariance_1pass_create (n_all_vars, all_vars,
729 dict_get_weight (dataset_dict (cmd->ds)),
730 MV_ANY, cmd->origin == false);
732 reader = casereader_clone (input);
733 reader = casereader_create_filter_missing (reader, all_vars, n_all_vars,
736 struct casereader *r = casereader_clone (reader);
738 for (; (c = casereader_read (r)) != NULL; case_unref (c))
740 covariance_accumulate (cov, c);
742 casereader_destroy (r);
745 struct linreg **models = XCALLOC (cmd->n_dep_vars, struct linreg*);
747 for (int k = 0; k < cmd->n_dep_vars; k++)
749 const struct variable **vars = xnmalloc (cmd->n_vars, sizeof (*vars));
750 const struct variable *dep_var = cmd->dep_vars[k];
751 int n_indep = identify_indep_vars (cmd, vars, dep_var);
752 gsl_matrix *cov_matrix = gsl_matrix_alloc (n_indep + 1, n_indep + 1);
753 double n_data = fill_covariance (cov_matrix, cov, vars, n_indep,
754 dep_var, all_vars, n_all_vars, means);
755 models[k] = linreg_alloc (dep_var, vars, n_data, n_indep, cmd->origin);
756 for (i = 0; i < n_indep; i++)
758 linreg_set_indep_variable_mean (models[k], i, means[i]);
760 linreg_set_depvar_mean (models[k], means[i]);
763 linreg_fit (cov_matrix, models[k]);
765 if (output && !taint_has_tainted_successor (casereader_get_taint (input)))
768 Find the least-squares estimates and other statistics.
770 if (cmd->stats & STATS_R)
771 reg_stats_r (models[k], dep_var);
773 if (cmd->stats & STATS_ANOVA)
774 reg_stats_anova (models[k], dep_var);
776 if (cmd->stats & STATS_COEFF)
777 reg_stats_coeff (cmd, models[k],
779 cov_matrix, dep_var);
781 if (cmd->stats & STATS_BCOV)
782 reg_stats_bcov (models[k], dep_var);
787 msg (SE, _("No valid data found. This command was skipped."));
790 gsl_matrix_free (cov_matrix);
793 casereader_destroy (reader);
795 for (int i = 0; i < cmd->n_vars; i++)
797 if (model_container[i].models)
799 linreg_unref (model_container[i].models[0]);
801 free (model_container[i].models);
803 free (model_container);
807 covariance_destroy (cov);
812 run_regression (const struct regression *cmd,
813 struct regression_workspace *ws,
814 struct casereader *input)
816 struct linreg **models = run_regression_get_models (cmd, input, true);
821 struct casereader *r = casereader_clone (input);
823 for (; (c = casereader_read (r)) != NULL; case_unref (c))
825 struct ccase *outc = case_create (casewriter_get_proto (ws->writer));
826 for (int k = 0; k < cmd->n_dep_vars; k++)
828 const struct variable **vars = xnmalloc (cmd->n_vars, sizeof (*vars));
829 const struct variable *dep_var = cmd->dep_vars[k];
830 int n_indep = identify_indep_vars (cmd, vars, dep_var);
831 double *vals = xnmalloc (n_indep, sizeof (*vals));
832 for (int i = 0; i < n_indep; i++)
834 const union value *tmp = case_data (c, vars[i]);
840 double pred = linreg_predict (models[k], vals, n_indep);
841 *case_num_rw_idx (outc, k * ws->extras + ws->pred_idx) = pred;
846 double obs = case_num (c, linreg_dep_var (models[k]));
847 double res = linreg_residual (models[k], obs, vals, n_indep);
848 *case_num_rw_idx (outc, k * ws->extras + ws->res_idx) = res;
853 casewriter_write (ws->writer, outc);
855 casereader_destroy (r);
858 for (int k = 0; k < cmd->n_dep_vars; k++)
860 linreg_unref (models[k]);
864 casereader_destroy (input);
871 reg_stats_r (const struct linreg * c, const struct variable *var)
873 struct pivot_table *table = pivot_table_create__ (
874 pivot_value_new_text_format (N_("Model Summary (%s)"),
875 var_to_string (var)),
878 pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
879 N_("R"), N_("R Square"), N_("Adjusted R Square"),
880 N_("Std. Error of the Estimate"));
882 double rsq = linreg_ssreg (c) / linreg_sst (c);
883 double adjrsq = (rsq -
884 (1.0 - rsq) * linreg_n_coeffs (c)
885 / (linreg_n_obs (c) - linreg_n_coeffs (c) - 1));
886 double std_error = sqrt (linreg_mse (c));
889 sqrt (rsq), rsq, adjrsq, std_error
891 for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
892 pivot_table_put1 (table, i, pivot_value_new_number (entries[i]));
894 pivot_table_submit (table);
898 Table showing estimated regression coefficients.
901 reg_stats_coeff (const struct regression *cmd, const struct linreg *c,
902 const struct model_container *mc, const gsl_matrix *cov,
903 const struct variable *var)
905 struct pivot_table *table = pivot_table_create__ (
906 pivot_value_new_text_format (N_("Coefficients (%s)"), var_to_string (var)),
909 struct pivot_dimension *statistics = pivot_dimension_create (
910 table, PIVOT_AXIS_COLUMN, N_("Statistics"));
911 pivot_category_create_group (statistics->root,
912 N_("Unstandardized Coefficients"),
913 N_("B"), N_("Std. Error"));
914 pivot_category_create_group (statistics->root,
915 N_("Standardized Coefficients"), N_("Beta"));
916 pivot_category_create_leaves (statistics->root, N_("t"),
917 N_("Sig."), PIVOT_RC_SIGNIFICANCE);
918 if (cmd->stats & STATS_CI)
920 struct pivot_category *interval = pivot_category_create_group__ (
921 statistics->root, pivot_value_new_text_format (
922 N_("%g%% Confidence Interval for B"),
924 pivot_category_create_leaves (interval, N_("Lower Bound"),
928 if (cmd->stats & STATS_TOL)
929 pivot_category_create_group (statistics->root,
930 N_("Collinearity Statistics"),
931 N_("Tolerance"), N_("VIF"));
934 struct pivot_dimension *variables = pivot_dimension_create (
935 table, PIVOT_AXIS_ROW, N_("Variables"));
937 double df = linreg_n_obs (c) - linreg_n_coeffs (c) - 1;
938 double q = (1 - cmd->ci) / 2.0; /* 2-tailed test */
939 double tval = gsl_cdf_tdist_Qinv (q, df);
943 int var_idx = pivot_category_create_leaf (
944 variables->root, pivot_value_new_text (N_("(Constant)")));
946 double std_err = sqrt (gsl_matrix_get (linreg_cov (c), 0, 0));
947 double t_stat = linreg_intercept (c) / std_err;
948 double base_entries[] = {
949 linreg_intercept (c),
953 2.0 * gsl_cdf_tdist_Q (fabs (t_stat),
954 linreg_n_obs (c) - linreg_n_coeffs (c)),
958 for (size_t i = 0; i < sizeof base_entries / sizeof *base_entries; i++)
959 pivot_table_put2 (table, col++, var_idx,
960 pivot_value_new_number (base_entries[i]));
962 if (cmd->stats & STATS_CI)
964 double interval_entries[] = {
965 linreg_intercept (c) - tval * std_err,
966 linreg_intercept (c) + tval * std_err,
969 for (size_t i = 0; i < sizeof interval_entries / sizeof *interval_entries; i++)
970 pivot_table_put2 (table, col++, var_idx,
971 pivot_value_new_number (interval_entries[i]));
975 for (size_t j = 0; j < linreg_n_coeffs (c); j++)
977 const struct variable *v = linreg_indep_var (c, j);
978 int var_idx = pivot_category_create_leaf (
979 variables->root, pivot_value_new_variable (v));
981 double std_err = sqrt (gsl_matrix_get (linreg_cov (c), j + 1, j + 1));
982 double t_stat = linreg_coeff (c, j) / std_err;
983 double base_entries[] = {
985 sqrt (gsl_matrix_get (linreg_cov (c), j + 1, j + 1)),
986 (sqrt (gsl_matrix_get (cov, j, j)) * linreg_coeff (c, j) /
987 sqrt (gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1))),
989 2 * gsl_cdf_tdist_Q (fabs (t_stat), df)
993 for (size_t i = 0; i < sizeof base_entries / sizeof *base_entries; i++)
994 pivot_table_put2 (table, col++, var_idx,
995 pivot_value_new_number (base_entries[i]));
997 if (cmd->stats & STATS_CI)
999 double interval_entries[] = {
1000 linreg_coeff (c, j) - tval * std_err,
1001 linreg_coeff (c, j) + tval * std_err,
1005 for (size_t i = 0; i < sizeof interval_entries / sizeof *interval_entries; i++)
1006 pivot_table_put2 (table, col++, var_idx,
1007 pivot_value_new_number (interval_entries[i]));
1010 if (cmd->stats & STATS_TOL)
1013 struct linreg *m = mc[j].models[0];
1014 double rsq = linreg_ssreg (m) / linreg_sst (m);
1015 pivot_table_put2 (table, col++, var_idx, pivot_value_new_number (1.0 - rsq));
1016 pivot_table_put2 (table, col++, var_idx, pivot_value_new_number (1.0 / (1.0 - rsq)));
1021 pivot_table_submit (table);
1025 Display the ANOVA table.
1028 reg_stats_anova (const struct linreg * c, const struct variable *var)
1030 struct pivot_table *table = pivot_table_create__ (
1031 pivot_value_new_text_format (N_("ANOVA (%s)"), var_to_string (var)),
1034 pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
1035 N_("Sum of Squares"), PIVOT_RC_OTHER,
1036 N_("df"), PIVOT_RC_INTEGER,
1037 N_("Mean Square"), PIVOT_RC_OTHER,
1038 N_("F"), PIVOT_RC_OTHER,
1039 N_("Sig."), PIVOT_RC_SIGNIFICANCE);
1041 pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Source"),
1042 N_("Regression"), N_("Residual"), N_("Total"));
1044 double msm = linreg_ssreg (c) / linreg_dfmodel (c);
1045 double mse = linreg_mse (c);
1046 double F = msm / mse;
1055 /* Sums of Squares. */
1056 { 0, 0, linreg_ssreg (c) },
1057 { 0, 1, linreg_sse (c) },
1058 { 0, 2, linreg_sst (c) },
1059 /* Degrees of freedom. */
1060 { 1, 0, linreg_dfmodel (c) },
1061 { 1, 1, linreg_dferror (c) },
1062 { 1, 2, linreg_dftotal (c) },
1069 { 4, 0, gsl_cdf_fdist_Q (F, linreg_dfmodel (c), linreg_dferror (c)) },
1071 for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
1073 const struct entry *e = &entries[i];
1074 pivot_table_put2 (table, e->stat_idx, e->source_idx,
1075 pivot_value_new_number (e->x));
1078 pivot_table_submit (table);
1083 reg_stats_bcov (const struct linreg * c, const struct variable *var)
1085 struct pivot_table *table = pivot_table_create__ (
1086 pivot_value_new_text_format (N_("Coefficient Correlations (%s)"),
1087 var_to_string (var)),
1088 "Coefficient Correlations");
1090 for (size_t i = 0; i < 2; i++)
1092 struct pivot_dimension *models = pivot_dimension_create (
1093 table, i ? PIVOT_AXIS_ROW : PIVOT_AXIS_COLUMN, N_("Models"));
1094 for (size_t j = 0; j < linreg_n_coeffs (c); j++)
1095 pivot_category_create_leaf (
1096 models->root, pivot_value_new_variable (
1097 linreg_indep_var (c, j)));
1100 pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
1103 for (size_t i = 0; i < linreg_n_coeffs (c); i++)
1104 for (size_t k = 0; k < linreg_n_coeffs (c); k++)
1106 double cov = gsl_matrix_get (linreg_cov (c), MIN (i, k), MAX (i, k));
1107 pivot_table_put3 (table, k, i, 0, pivot_value_new_number (cov));
1110 pivot_table_submit (table);