1 /* PSPP - linear regression.
2 Copyright (C) 2005 Free Software Foundation, Inc.
3 Written by Jason H Stover <jason@sakla.net>.
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
10 This program is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 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, write to the Free Software
17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
22 #include <gsl/gsl_cdf.h>
23 #include <gsl/gsl_vector.h>
24 #include <gsl/gsl_matrix.h>
29 #include "cat-routines.h"
31 #include "design-matrix.h"
32 #include "dictionary.h"
34 #include "file-handle.h"
37 #include <linreg/pspp_linreg.h>
38 #include "missing-values.h"
39 #include "regression_export.h"
41 #include "value-labels.h"
45 #define REG_LARGE_DATA 1000
50 "REGRESSION" (regression_):
75 static struct cmd_regression cmd;
78 Array holding the subscripts of the independent variables.
83 File where the model will be saved if the EXPORT subcommand
86 struct file_handle *model_file;
89 Return value for the procedure.
91 int pspp_reg_rc = CMD_SUCCESS;
93 static void run_regression (const struct casefile *, void *);
95 STATISTICS subcommand output functions.
97 static void reg_stats_r (pspp_linreg_cache *);
98 static void reg_stats_coeff (pspp_linreg_cache *);
99 static void reg_stats_anova (pspp_linreg_cache *);
100 static void reg_stats_outs (pspp_linreg_cache *);
101 static void reg_stats_zpp (pspp_linreg_cache *);
102 static void reg_stats_label (pspp_linreg_cache *);
103 static void reg_stats_sha (pspp_linreg_cache *);
104 static void reg_stats_ci (pspp_linreg_cache *);
105 static void reg_stats_f (pspp_linreg_cache *);
106 static void reg_stats_bcov (pspp_linreg_cache *);
107 static void reg_stats_ses (pspp_linreg_cache *);
108 static void reg_stats_xtx (pspp_linreg_cache *);
109 static void reg_stats_collin (pspp_linreg_cache *);
110 static void reg_stats_tol (pspp_linreg_cache *);
111 static void reg_stats_selection (pspp_linreg_cache *);
112 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
113 int, pspp_linreg_cache *);
116 reg_stats_r (pspp_linreg_cache * c)
126 rsq = c->ssm / c->sst;
127 adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
128 std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
129 t = tab_create (n_cols, n_rows, 0);
130 tab_dim (t, tab_natural_dimensions);
131 tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
132 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
133 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
134 tab_vline (t, TAL_0, 1, 0, 0);
136 tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
137 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
138 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
139 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
140 tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
141 tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
142 tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
143 tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
144 tab_title (t, 0, _("Model Summary"));
149 Table showing estimated regression coefficients.
152 reg_stats_coeff (pspp_linreg_cache * c)
167 n_rows = c->n_coeffs + 2;
169 t = tab_create (n_cols, n_rows, 0);
170 tab_headers (t, 2, 0, 1, 0);
171 tab_dim (t, tab_natural_dimensions);
172 tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
173 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
174 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
175 tab_vline (t, TAL_0, 1, 0, 0);
177 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
178 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
179 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
180 tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
181 tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
182 tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
183 coeff = c->coeff[0].estimate;
184 tab_float (t, 2, 1, 0, coeff, 10, 2);
185 std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
186 tab_float (t, 3, 1, 0, std_err, 10, 2);
187 beta = coeff / c->depvar_std;
188 tab_float (t, 4, 1, 0, beta, 10, 2);
189 t_stat = coeff / std_err;
190 tab_float (t, 5, 1, 0, t_stat, 10, 2);
191 pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
192 tab_float (t, 6, 1, 0, pval, 10, 2);
193 for (j = 1; j <= c->n_indeps; j++)
196 label = var_to_string (c->coeff[j].v);
197 tab_text (t, 1, j + 1, TAB_CENTER, label);
199 Regression coefficients.
201 coeff = c->coeff[j].estimate;
202 tab_float (t, 2, j + 1, 0, coeff, 10, 2);
204 Standard error of the coefficients.
206 std_err = sqrt (gsl_matrix_get (c->cov, j, j));
207 tab_float (t, 3, j + 1, 0, std_err, 10, 2);
209 'Standardized' coefficient, i.e., regression coefficient
210 if all variables had unit variance.
212 beta = gsl_vector_get (c->indep_std, j);
213 beta *= coeff / c->depvar_std;
214 tab_float (t, 4, j + 1, 0, beta, 10, 2);
217 Test statistic for H0: coefficient is 0.
219 t_stat = coeff / std_err;
220 tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
222 P values for the test statistic above.
224 pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
225 tab_float (t, 6, j + 1, 0, pval, 10, 2);
227 tab_title (t, 0, _("Coefficients"));
232 Display the ANOVA table.
235 reg_stats_anova (pspp_linreg_cache * c)
239 const double msm = c->ssm / c->dfm;
240 const double mse = c->sse / c->dfe;
241 const double F = msm / mse;
242 const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
247 t = tab_create (n_cols, n_rows, 0);
248 tab_headers (t, 2, 0, 1, 0);
249 tab_dim (t, tab_natural_dimensions);
251 tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
253 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
254 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
255 tab_vline (t, TAL_0, 1, 0, 0);
257 tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
258 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
259 tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
260 tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
261 tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
263 tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
264 tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
265 tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
267 /* Sums of Squares */
268 tab_float (t, 2, 1, 0, c->ssm, 10, 2);
269 tab_float (t, 2, 3, 0, c->sst, 10, 2);
270 tab_float (t, 2, 2, 0, c->sse, 10, 2);
273 /* Degrees of freedom */
274 tab_float (t, 3, 1, 0, c->dfm, 4, 0);
275 tab_float (t, 3, 2, 0, c->dfe, 4, 0);
276 tab_float (t, 3, 3, 0, c->dft, 4, 0);
280 tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
281 tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
283 tab_float (t, 5, 1, 0, F, 8, 3);
285 tab_float (t, 6, 1, 0, pval, 8, 3);
287 tab_title (t, 0, _("ANOVA"));
291 reg_stats_outs (pspp_linreg_cache * c)
296 reg_stats_zpp (pspp_linreg_cache * c)
301 reg_stats_label (pspp_linreg_cache * c)
306 reg_stats_sha (pspp_linreg_cache * c)
311 reg_stats_ci (pspp_linreg_cache * c)
316 reg_stats_f (pspp_linreg_cache * c)
321 reg_stats_bcov (pspp_linreg_cache * c)
334 n_cols = c->n_indeps + 1 + 2;
335 n_rows = 2 * (c->n_indeps + 1);
336 t = tab_create (n_cols, n_rows, 0);
337 tab_headers (t, 2, 0, 1, 0);
338 tab_dim (t, tab_natural_dimensions);
339 tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
340 tab_hline (t, TAL_2, 0, n_cols - 1, 1);
341 tab_vline (t, TAL_2, 2, 0, n_rows - 1);
342 tab_vline (t, TAL_0, 1, 0, 0);
343 tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
344 tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
345 for (i = 1; i < c->n_indeps + 1; i++)
347 j = indep_vars[(i - 1)];
348 struct variable *v = cmd.v_variables[j];
349 label = var_to_string (v);
350 tab_text (t, 2, i, TAB_CENTER, label);
351 tab_text (t, i + 2, 0, TAB_CENTER, label);
352 for (k = 1; k < c->n_indeps + 1; k++)
354 col = (i <= k) ? k : i;
355 row = (i <= k) ? i : k;
356 tab_float (t, k + 2, i, TAB_CENTER,
357 gsl_matrix_get (c->cov, row, col), 8, 3);
360 tab_title (t, 0, _("Coefficient Correlations"));
364 reg_stats_ses (pspp_linreg_cache * c)
369 reg_stats_xtx (pspp_linreg_cache * c)
374 reg_stats_collin (pspp_linreg_cache * c)
379 reg_stats_tol (pspp_linreg_cache * c)
384 reg_stats_selection (pspp_linreg_cache * c)
390 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
391 int keyword, pspp_linreg_cache * c)
400 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
403 The order here must match the order in which the STATISTICS
404 keywords appear in the specification section above.
431 Set everything but F.
433 for (i = 0; i < f; i++)
440 for (i = 0; i < all; i++)
448 Default output: ANOVA table, parameter estimates,
449 and statistics for variables not entered into model,
452 if (keywords[defaults] | d)
460 statistics_keyword_output (reg_stats_r, keywords[r], c);
461 statistics_keyword_output (reg_stats_anova, keywords[anova], c);
462 statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
463 statistics_keyword_output (reg_stats_outs, keywords[outs], c);
464 statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
465 statistics_keyword_output (reg_stats_label, keywords[label], c);
466 statistics_keyword_output (reg_stats_sha, keywords[sha], c);
467 statistics_keyword_output (reg_stats_ci, keywords[ci], c);
468 statistics_keyword_output (reg_stats_f, keywords[f], c);
469 statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
470 statistics_keyword_output (reg_stats_ses, keywords[ses], c);
471 statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
472 statistics_keyword_output (reg_stats_collin, keywords[collin], c);
473 statistics_keyword_output (reg_stats_tol, keywords[tol], c);
474 statistics_keyword_output (reg_stats_selection, keywords[selection], c);
477 reg_inserted (struct variable *v, struct variable **varlist, int n_vars)
481 for (i = 0; i < n_vars; i++)
483 if (v->index == varlist[i]->index)
491 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
496 struct variable **varlist;
497 struct pspp_linreg_coeff coeff;
500 fprintf (fp, "%s", reg_export_categorical_encode_1);
502 varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
503 for (i = 1; i < c->n_indeps; i++) /* c->coeff[0] is the intercept. */
506 if (coeff.v->type == ALPHA)
508 if (!reg_inserted (coeff.v, varlist, n_vars))
510 fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
512 varlist[n_vars] = coeff.v;
517 fprintf (fp, "int n_vars = %d;\n\t", n_vars);
518 fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
520 for (i = 0; i < n_vars - 1; i++)
522 fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
524 fprintf (fp, "&%s};\n\t", varlist[i]->name);
526 for (i = 0; i < n_vars; i++)
529 fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
531 fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
532 varlist[i]->obs_vals->n_categories);
534 for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
536 val = cat_subscript_to_value ((const size_t) j, varlist[i]);
537 fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
538 value_to_string (val, varlist[i]));
541 fprintf (fp, "%s", reg_export_categorical_encode_2);
545 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
548 struct pspp_linreg_coeff coeff;
550 fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
551 for (i = 1; i < c->n_indeps; i++)
554 fprintf (fp, "\"%s\",\n\t\t", coeff.v->name);
557 fprintf (fp, "\"%s\"};\n\t", coeff.v->name);
560 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
562 fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
563 fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
564 reg_print_depvars (fp, c);
565 fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
567 "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
568 fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
571 subcommand_export (int export, pspp_linreg_cache * c)
576 int n_quantiles = 100;
579 struct pspp_linreg_coeff coeff;
584 assert (model_file != NULL);
586 fp = fopen (handle_get_filename (model_file), "w");
587 fprintf (fp, "%s", reg_preamble);
588 reg_print_getvar (fp, c);
589 reg_print_categorical_encoding (fp, c);
590 fprintf (fp, "%s", reg_export_t_quantiles_1);
591 increment = 0.5 / (double) increment;
592 for (i = 0; i < n_quantiles - 1; i++)
594 tmp = 0.5 + 0.005 * (double) i;
595 fprintf (fp, "%.15e,\n\t\t",
596 gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
598 fprintf (fp, "%.15e};\n\t",
599 gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
600 fprintf (fp, "%s", reg_export_t_quantiles_2);
601 fprintf (fp, "%s", reg_mean_cmt);
602 fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
603 fprintf (fp, "const char *var_names[])\n{\n\t");
604 fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
605 for (i = 1; i < c->n_indeps; i++)
608 fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
611 fprintf (fp, "%.15e};\n\t", coeff.estimate);
613 fprintf (fp, "double estimate = %.15e;\n\t", coeff.estimate);
614 fprintf (fp, "int i;\n\tint j;\n\n\t");
615 fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
616 fprintf (fp, "%s", reg_getvar);
617 fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs,
619 for (i = 0; i < c->cov->size1 - 1; i++)
622 for (j = 0; j < c->cov->size2 - 1; j++)
624 fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
626 fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
629 for (j = 0; j < c->cov->size2 - 1; j++)
631 fprintf (fp, "%.15e, ",
632 gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
634 fprintf (fp, "%.15e}\n\t",
635 gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
636 fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t",
638 fprintf (fp, "double unshuffled_vals[%d];\n\t", c->n_indeps);
639 fprintf (fp, "%s", reg_variance);
640 fprintf (fp, "%s", reg_export_confidence_interval);
641 tmp = c->mse * c->mse;
642 fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
643 fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
644 fprintf (fp, "%s", reg_export_prediction_interval_3);
646 fp = fopen ("pspp_model_reg.h", "w");
647 fprintf (fp, "%s", reg_header);
652 regression_custom_export (struct cmd_regression *cmd)
654 /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
655 if (!lex_force_match ('('))
662 model_file = fh_parse ();
663 if (model_file == NULL)
667 if (!lex_force_match (')'))
674 cmd_regression (void)
676 if (!parse_regression (&cmd))
680 multipass_procedure_with_splits (run_regression, &cmd);
686 Is variable k one of the dependent variables?
692 for (j = 0; j < cmd.n_dependent; j++)
695 compare_var_names returns 0 if the variable
698 if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
705 Mark missing cases. Return the number of non-missing cases.
708 mark_missing_cases (const struct casefile *cf, struct variable *v,
709 double *is_missing_case, double n_data)
711 struct casereader *r;
716 for (r = casefile_get_reader (cf);
717 casereader_read (r, &c); case_destroy (&c))
719 row = casereader_cnum (r) - 1;
721 val = case_data (&c, v->fv);
722 cat_value_update (v, val);
723 if (mv_is_value_missing (&v->miss, val))
725 if (!is_missing_case[row])
727 /* Now it is missing. */
729 is_missing_case[row] = 1;
733 casereader_destroy (r);
738 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
748 Keep track of the missing cases.
750 int *is_missing_case;
751 const union value *val;
752 struct casereader *r;
755 struct variable *depvar;
756 struct variable **indep_vars;
757 struct design_matrix *X;
759 pspp_linreg_cache *lcache;
760 pspp_linreg_opts lopts;
762 n_data = casefile_get_case_cnt (cf);
764 for (i = 0; i < cmd.n_dependent; i++)
766 if (cmd.v_dependent[i]->type != NUMERIC)
768 msg (SE, gettext ("Dependent variable must be numeric."));
769 pspp_reg_rc = CMD_FAILURE;
774 is_missing_case = xnmalloc (n_data, sizeof (*is_missing_case));
775 for (i = 0; i < n_data; i++)
776 is_missing_case[i] = 0;
778 n_indep = cmd.n_variables - cmd.n_dependent;
779 indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
781 lopts.get_depvar_mean_std = 1;
782 lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
785 Read from the active file. The first pass encodes categorical
786 variables and drops cases with missing values.
789 for (i = 0; i < cmd.n_variables; i++)
793 v = cmd.v_variables[i];
796 if (v->type == ALPHA)
798 /* Make a place to hold the binary vectors
799 corresponding to this variable's values. */
800 cat_stored_values_create (v);
802 n_data = mark_missing_cases (cf, v, is_missing_case, n_data);
807 Drop cases with missing values for any dependent variable.
810 for (i = 0; i < cmd.n_dependent; i++)
812 v = cmd.v_dependent[i];
814 n_data = mark_missing_cases (cf, v, is_missing_case, n_data);
817 for (k = 0; k < cmd.n_dependent; k++)
819 depvar = cmd.v_dependent[k];
820 Y = gsl_vector_alloc (n_data);
823 design_matrix_create (n_indep, (const struct variable **) indep_vars,
825 lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
826 lcache->indep_means = gsl_vector_alloc (X->m->size2);
827 lcache->indep_std = gsl_vector_alloc (X->m->size2);
828 lcache->depvar = (const struct variable *) depvar;
830 For large data sets, use QR decomposition.
832 if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
834 lcache->method = PSPP_LINREG_SVD;
838 The second pass creates the design matrix.
841 for (r = casefile_get_reader (cf); casereader_read (r, &c);
843 /* Iterate over the cases. */
845 case_num = casereader_cnum (r) - 1;
846 if (!is_missing_case[case_num])
848 for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
849 for the current case.
852 v = cmd.v_variables[i];
853 val = case_data (&c, v->fv);
855 Independent/dependent variable separation. The
856 'variables' subcommand specifies a varlist which contains
857 both dependent and independent variables. The dependent
858 variables are specified with the 'dependent'
859 subcommand, and maybe also in the 'variables' subcommand.
860 We need to separate the two.
864 if (v->type == ALPHA)
866 design_matrix_set_categorical (X, row, v, val);
868 else if (v->type == NUMERIC)
870 design_matrix_set_numeric (X, row, v, val);
872 lopts.get_indep_mean_std[i] = 1;
875 val = case_data (&c, depvar->fv);
876 gsl_vector_set (Y, row, val->f);
881 Now that we know the number of coefficients, allocate space
882 and store pointers to the variables that correspond to the
885 lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
886 for (i = 0; i < X->m->size2; i++)
888 j = i + 1; /* The first coeff is the intercept. */
890 (const struct variable *) design_matrix_col_to_var (X, i);
891 assert (lcache->coeff[j].v != NULL);
895 Find the least-squares estimates and other statistics.
897 pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
898 subcommand_statistics (cmd.a_statistics, lcache);
899 subcommand_export (cmd.sbc_export, lcache);
901 design_matrix_destroy (X);
902 pspp_linreg_cache_free (lcache);
903 free (lopts.get_indep_mean_std);
904 casereader_destroy (r);
907 free (is_missing_case);