REGRESSION: Implement /ORIGIN subcommand.
[pspp] / src / language / stats / glm.c
index d2b3a37e862e39ac44587975037d26b38eab067c..74e918b886b21f26d05b1fadf2d245faa60096bb 100644 (file)
@@ -81,7 +81,7 @@ struct glm_workspace
 
   struct categoricals *cats;
 
-  /* 
+  /*
      Sums of squares due to different variables. Element 0 is the SSE
      for the entire model. For i > 0, element i is the SS due to
      variable i.
@@ -155,7 +155,8 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
                              PV_NO_DUPLICATE | PV_NUMERIC))
     goto error;
 
-  lex_force_match (lexer, T_BY);
+  if (! lex_force_match (lexer, T_BY))
+    goto error;
 
   if (!parse_variables_const (lexer, glm.dict,
                              &glm.factor_vars, &glm.n_factor_vars,
@@ -229,7 +230,7 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
                      lex_error (lexer, NULL);
                      goto error;
                    }
-                 
+
                  glm.alpha = lex_number (lexer);
                  lex_get (lexer);
                  if ( ! lex_force_match (lexer, T_RPAREN))
@@ -267,7 +268,7 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
            }
 
          glm.ss_type = lex_integer (lexer);
-         if (1 > glm.ss_type  && 3 < glm.ss_type )
+         if (1 > glm.ss_type  ||  3 < glm.ss_type )
            {
              msg (ME, _("Only types 1, 2 & 3 sums of squares are currently implemented"));
              goto error;
@@ -324,6 +325,9 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
 
   const_var_set_destroy (factors);
   free (glm.factor_vars);
+  for (i = 0 ; i < glm.n_interactions; ++i)
+    interaction_destroy (glm.interactions[i]);
+
   free (glm.interactions);
   free (glm.dep_vars);
 
@@ -356,11 +360,11 @@ fill_submatrix (const gsl_matrix * cov, gsl_matrix * submatrix, bool *dropped_f)
   size_t j;
   size_t n = 0;
   size_t m = 0;
-  
+
   for (i = 0; i < cov->size1; i++)
     {
       if (not_dropped (i, dropped_f))
-       {         
+       {
          m = 0;
          for (j = 0; j < cov->size2; j++)
            {
@@ -369,7 +373,7 @@ fill_submatrix (const gsl_matrix * cov, gsl_matrix * submatrix, bool *dropped_f)
                  gsl_matrix_set (submatrix, n, m,
                                  gsl_matrix_get (cov, i, j));
                  m++;
-               }       
+               }
            }
          n++;
        }
@@ -377,14 +381,14 @@ fill_submatrix (const gsl_matrix * cov, gsl_matrix * submatrix, bool *dropped_f)
 }
 
 
-/* 
+/*
    Type 1 sums of squares.
    Populate SSQ with the Type 1 sums of squares according to COV
  */
 static void
 ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 {
-  gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+  const gsl_matrix *cm = covariance_calculate_unnormalized (cov);
   size_t i;
   size_t k;
   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
@@ -406,7 +410,7 @@ ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
     {
       gsl_matrix *model_cov = NULL;
       gsl_matrix *submodel_cov = NULL;
-      
+
       n_dropped_submodel = n_dropped_model;
       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
        {
@@ -415,7 +419,7 @@ ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 
       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
        {
-         const struct interaction * x = 
+         const struct interaction * x =
            categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
 
          if ( x == cmd->interactions [k])
@@ -444,17 +448,16 @@ ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 
   free (model_dropped);
   free (submodel_dropped);
-  gsl_matrix_free (cm);
 }
 
-/* 
+/*
    Type 2 sums of squares.
    Populate SSQ with the Type 2 sums of squares according to COV
  */
 static void
 ssq_type2 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 {
-  gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+  const gsl_matrix *cm = covariance_calculate_unnormalized (cov);
   size_t i;
   size_t k;
   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
@@ -469,7 +472,7 @@ ssq_type2 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
       size_t n_dropped_submodel = 0;
       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
        {
-         const struct interaction * x = 
+         const struct interaction * x =
            categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
 
          model_dropped[i] = false;
@@ -508,17 +511,16 @@ ssq_type2 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 
   free (model_dropped);
   free (submodel_dropped);
-  gsl_matrix_free (cm);
 }
 
-/* 
+/*
    Type 3 sums of squares.
    Populate SSQ with the Type 2 sums of squares according to COV
  */
 static void
 ssq_type3 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 {
-  gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+  const gsl_matrix *cm = covariance_calculate_unnormalized (cov);
   size_t i;
   size_t k;
   bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
@@ -540,7 +542,7 @@ ssq_type3 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 
       for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
        {
-         const struct interaction * x = 
+         const struct interaction * x =
            categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
 
          model_dropped[i] = false;
@@ -565,8 +567,6 @@ ssq_type3 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
       gsl_matrix_free (model_cov);
     }
   free (model_dropped);
-
-  gsl_matrix_free (cm);
 }
 
 
@@ -587,11 +587,21 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
   struct glm_workspace ws;
   struct covariance *cov;
 
+  input  = casereader_create_filter_missing (input,
+                                            cmd->dep_vars, cmd->n_dep_vars,
+                                            cmd->exclude,
+                                            NULL,  NULL);
+
+  input  = casereader_create_filter_missing (input,
+                                            cmd->factor_vars, cmd->n_factor_vars,
+                                            cmd->exclude,
+                                            NULL,  NULL);
+
   ws.cats = categoricals_create (cmd->interactions, cmd->n_interactions,
                                 cmd->wv, cmd->exclude, MV_ANY);
 
   cov = covariance_2pass_create (cmd->n_dep_vars, cmd->dep_vars,
-                                ws.cats, cmd->wv, cmd->exclude);
+                                ws.cats, cmd->wv, cmd->exclude, true);
 
 
   c = casereader_peek (input, 0);
@@ -654,7 +664,9 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
     }
 
   {
-    gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+    const gsl_matrix *ucm = covariance_calculate_unnormalized (cov);
+    gsl_matrix *cm = gsl_matrix_alloc (ucm->size1, ucm->size2);
+    gsl_matrix_memcpy (cm, ucm);
 
     //    dump_matrix (cm);
 
@@ -683,7 +695,6 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
        break;
       }
     //    dump_matrix (cm);
-
     gsl_matrix_free (cm);
   }
 
@@ -698,7 +709,7 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
   taint_destroy (taint);
 }
 
-static const char *roman[] = 
+static const char *roman[] =
   {
     "", /* The Romans had no concept of zero */
     "I",
@@ -730,8 +741,8 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
   if (cmd->intercept)
     nr += 2;
 
-  msg (MW, "GLM is experimental.  Do not rely on these results.");
   t = tab_create (nc, nr);
+  tab_set_format (t, RC_WEIGHT, wfmt);
   tab_title (t, _("Tests of Between-Subjects Effects"));
 
   tab_headers (t, heading_columns, 0, heading_rows, 0);
@@ -745,7 +756,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 
   /* TRANSLATORS: The parameter is a roman numeral */
   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE,
-                  _("Type %s Sum of Squares"), 
+                  _("Type %s Sum of Squares"),
                   roman[cmd->ss_type]);
   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
@@ -774,12 +785,17 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
       const double df = 1.0;
       const double F = intercept_ssq / df / mse;
       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Intercept"));
-      tab_double (t, 1, r, 0, intercept_ssq, NULL);
-      tab_double (t, 2, r, 0, 1.00, wfmt);
-      tab_double (t, 3, r, 0, intercept_ssq / df, NULL);
-      tab_double (t, 4, r, 0, F, NULL);
-      tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
-                 NULL);
+      /* The intercept for unbalanced models is of limited use and
+        nobody knows how to calculate it properly */
+      if (categoricals_isbalanced (ws->cats))
+       {
+         tab_double (t, 1, r, 0, intercept_ssq, NULL, RC_OTHER);
+         tab_double (t, 2, r, 0, 1.00, NULL, RC_WEIGHT);
+         tab_double (t, 3, r, 0, intercept_ssq / df, NULL, RC_OTHER);
+         tab_double (t, 4, r, 0, F, NULL, RC_OTHER);
+         tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
+                     NULL, RC_PVALUE);
+       }
       r++;
     }
 
@@ -793,7 +809,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 
       ssq_effects += ssq;
 
-      if (! cmd->intercept) 
+      if (! cmd->intercept)
        {
          df++;
          ssq += intercept_ssq;
@@ -804,13 +820,13 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, ds_cstr (&str));
       ds_destroy (&str);
 
-      tab_double (t, 1, r, 0, ssq, NULL);
-      tab_double (t, 2, r, 0, df, wfmt);
-      tab_double (t, 3, r, 0, ssq / df, NULL);
-      tab_double (t, 4, r, 0, F, NULL);
+      tab_double (t, 1, r, 0, ssq, NULL, RC_OTHER);
+      tab_double (t, 2, r, 0, df, NULL, RC_WEIGHT);
+      tab_double (t, 3, r, 0, ssq / df, NULL, RC_OTHER);
+      tab_double (t, 4, r, 0, F, NULL, RC_OTHER);
 
       tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr),
-                 NULL);
+                 NULL, RC_PVALUE);
       r++;
     }
 
@@ -826,13 +842,13 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
       ssq += intercept_ssq;
 
     F = ssq / df / mse;
-    tab_double (t, 1, heading_rows, 0, ssq, NULL);
-    tab_double (t, 2, heading_rows, 0, df, wfmt);
-    tab_double (t, 3, heading_rows, 0, ssq / df, NULL);
-    tab_double (t, 4, heading_rows, 0, F, NULL);
+    tab_double (t, 1, heading_rows, 0, ssq, NULL, RC_OTHER);
+    tab_double (t, 2, heading_rows, 0, df, NULL, RC_WEIGHT);
+    tab_double (t, 3, heading_rows, 0, ssq / df, NULL, RC_OTHER);
+    tab_double (t, 4, heading_rows, 0, F, NULL, RC_OTHER);
 
     tab_double (t, 5, heading_rows, 0,
-               gsl_cdf_fdist_Q (F, df, n_total - df_corr), NULL);
+               gsl_cdf_fdist_Q (F, df, n_total - df_corr), NULL, RC_PVALUE);
   }
 
   {
@@ -840,24 +856,24 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
     const double ssq = gsl_vector_get (ws->ssq, 0);
     const double mse = ssq / df;
     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Error"));
-    tab_double (t, 1, r, 0, ssq, NULL);
-    tab_double (t, 2, r, 0, df, wfmt);
-    tab_double (t, 3, r++, 0, mse, NULL);
+    tab_double (t, 1, r, 0, ssq, NULL, RC_OTHER);
+    tab_double (t, 2, r, 0, df, NULL, RC_WEIGHT);
+    tab_double (t, 3, r++, 0, mse, NULL, RC_OTHER);
   }
 
   {
     tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Total"));
-    tab_double (t, 1, r, 0, ws->total_ssq + intercept_ssq, NULL);
-    tab_double (t, 2, r, 0, n_total, wfmt);
-    
+    tab_double (t, 1, r, 0, ws->total_ssq + intercept_ssq, NULL, RC_OTHER);
+    tab_double (t, 2, r, 0, n_total, NULL, RC_WEIGHT);
+
     r++;
   }
 
   if (cmd->intercept)
     {
       tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Corrected Total"));
-      tab_double (t, 1, r, 0, ws->total_ssq, NULL);
-      tab_double (t, 2, r, 0, n_total - 1.0, wfmt);
+      tab_double (t, 1, r, 0, ws->total_ssq, NULL, RC_OTHER);
+      tab_double (t, 2, r, 0, n_total - 1.0, NULL, RC_WEIGHT);
     }
 
   tab_submit (t);
@@ -883,66 +899,6 @@ dump_matrix (const gsl_matrix * m)
 
 
 \f
-
-/* Match a variable.
-   If the match succeeds, the variable will be placed in VAR.
-   Returns true if successful */
-static bool
-lex_match_variable (struct lexer *lexer, const struct dictionary *dict, const struct variable **var)
-{
-  if (lex_token (lexer) !=  T_ID)
-    return false;
-
-  *var = parse_variable_const  (lexer, dict);
-
-  if ( *var == NULL)
-    return false;
-  return true;
-}
-
-/* An interaction is a variable followed by {*, BY} followed by an interaction */
-static bool
-parse_design_interaction (struct lexer *lexer, const struct dictionary *dict, struct interaction **iact)
-{
-  const struct variable *v = NULL;
-  assert (iact);
-
-  switch  (lex_next_token (lexer, 1))
-    {
-    case T_ENDCMD:
-    case T_SLASH:
-    case T_COMMA:
-    case T_ID:
-    case T_BY:
-    case T_ASTERISK:
-      break;
-    default:
-      return false;
-      break;
-    }
-
-  if (! lex_match_variable (lexer, dict, &v))
-    {
-      interaction_destroy (*iact);
-      *iact = NULL;
-      return false;
-    }
-  
-  assert (v);
-
-  if ( *iact == NULL)
-    *iact = interaction_create (v);
-  else
-    interaction_add_variable (*iact, v);
-
-  if ( lex_match (lexer, T_ASTERISK) || lex_match (lexer, T_BY))
-    {
-      return parse_design_interaction (lexer, dict, iact);
-    }
-
-  return true;
-}
-
 static bool
 parse_nested_variable (struct lexer *lexer, struct glm_spec *glm)
 {
@@ -959,7 +915,7 @@ parse_nested_variable (struct lexer *lexer, struct glm_spec *glm)
        return false;
     }
 
-  lex_error (lexer, "Nested variables are not yet implemented"); return false;  
+  lex_error (lexer, "Nested variables are not yet implemented"); return false;
   return true;
 }