CROSSTABS: Calculate significance for lambda and Somers' d.
[pspp] / src / language / stats / crosstabs.q
index 5095b7471551bd57fad15da0db9eb16ff6a8e182..aa7e2457ee24ae9d41d40bd989da985b2c822fe3 100644 (file)
 
 /* FIXME:
 
-   - T values for Spearman's R and Pearson's R are wrong.
-   - How to calculate significance of symmetric and directional measures?
-   - Asymmetric ASEs and T values for lambda are wrong.
-   - ASE of Goodman and Kruskal's tau is not calculated.
-   - ASE of symmetric somers' d is wrong.
-   - Approx. T of uncertainty coefficient is wrong.
+   - How to calculate significance of some symmetric and directional measures?
+   - How to calculate ASE for symmetric Somers ' d?
+   - How to calculate ASE for Goodman and Kruskal's tau?
+   - How to calculate approx. T of symmetric uncertainty coefficient?
 
 */
 
@@ -1976,7 +1974,7 @@ display_risk (struct pivot_table *pt, struct tab_table *risk)
 
 static int calc_directional (struct crosstabs_proc *, struct pivot_table *,
                              double[N_DIRECTIONAL], double[N_DIRECTIONAL],
-                            double[N_DIRECTIONAL]);
+                            double[N_DIRECTIONAL], double[N_DIRECTIONAL]);
 
 /* Display directional measures. */
 static void
@@ -2043,10 +2041,11 @@ display_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
   double direct_v[N_DIRECTIONAL];
   double direct_ase[N_DIRECTIONAL];
   double direct_t[N_DIRECTIONAL];
+  double sig[N_DIRECTIONAL];
 
   int i;
 
-  if (!calc_directional (proc, pt, direct_v, direct_ase, direct_t))
+  if (!calc_directional (proc, pt, direct_v, direct_ase, direct_t, sig))
     return;
 
   tab_offset (direct, pt->n_consts + pt->n_vars - 2, -1);
@@ -2088,7 +2087,7 @@ display_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
        tab_double (direct, 4, 0, TAB_RIGHT, direct_ase[i], NULL, RC_OTHER);
       if (direct_t[i] != SYSMIS)
        tab_double (direct, 5, 0, TAB_RIGHT, direct_t[i], NULL, RC_OTHER);
-      /*tab_double (direct, 6, 0, TAB_RIGHT, normal_sig (direct_v[i]), NULL, RC_PVALUE);*/
+      tab_double (direct, 6, 0, TAB_RIGHT, sig[i], NULL, RC_PVALUE);
       tab_next_row (direct);
     }
 
@@ -2265,12 +2264,12 @@ calc_chisq (struct pivot_table *pt,
     }
 }
 
-/* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
-   ASE_1, and ase_0 into ASE_0.  The row and column values must be
+/* Calculate the value of Pearson's r.  r is stored into R, its T value into
+   T, and standard error into ERROR.  The row and column values must be
    passed in PT and Y. */
 static void
 calc_r (struct pivot_table *pt,
-        double *PT, double *Y, double *r, double *ase_0, double *ase_1)
+        double *PT, double *Y, double *r, double *t, double *error)
 {
   double SX, SY, S, T;
   double Xbar, Ybar;
@@ -2308,7 +2307,7 @@ calc_r (struct pivot_table *pt,
   SY = sum_Y2c - pow2 (sum_Yc) / pt->total;
   T = sqrt (SX * SY);
   *r = S / T;
-  *ase_0 = sqrt ((sum_X2Y2f - pow2 (sum_XYf) / pt->total) / (sum_X2r * sum_Y2c));
+  *t = *r / sqrt (1 - pow2 (*r)) * sqrt (pt->total - 2);
 
   {
     double s, c, y, t;
@@ -2329,7 +2328,7 @@ calc_r (struct pivot_table *pt,
          c = (t - s) - y;
          s = t;
        }
-    *ase_1 = sqrt (s) / (T * T);
+    *error = sqrt (s) / (T * T);
   }
 }
 
@@ -2536,7 +2535,7 @@ calc_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
       if (proc->statistics & (1u << CRS_ST_D))
        {
          somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
-         somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
+         somers_d_ase[0] = SYSMIS;
          somers_d_t[0] = (somers_d_v[0]
                           / (4 / (Dc + Dr)
                              * sqrt (ctau_cum - pow2 (P - Q) / pt->total)));
@@ -2596,13 +2595,11 @@ calc_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
       }
 
       calc_r (pt, R, C, &v[6], &t[6], &ase[6]);
-      t[6] = v[6] / t[6];
 
       free (R);
       free (C);
 
       calc_r (pt, (double *) pt->rows, (double *) pt->cols, &v[7], &t[7], &ase[7]);
-      t[7] = v[7] / t[7];
     }
 
   /* Cohen's kappa. */
@@ -2724,13 +2721,13 @@ calc_risk (struct pivot_table *pt,
 static int
 calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
                   double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
-                 double t[N_DIRECTIONAL])
+                 double t[N_DIRECTIONAL], double sig[N_DIRECTIONAL])
 {
   {
     int i;
 
     for (i = 0; i < N_DIRECTIONAL; i++)
-      v[i] = ase[i] = t[i] = SYSMIS;
+      v[i] = ase[i] = t[i] = sig[i] = SYSMIS;
   }
 
   /* Lambda. */
@@ -2805,19 +2802,14 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       /* ASE1 for Y given PT. */
       {
-       double accum;
-
-       for (accum = 0., i = 0; i < pt->n_rows; i++)
-         for (j = 0; j < pt->n_cols; j++)
-           {
-             const int deltaj = j == cm_index;
-             accum += (pt->mat[j + i * pt->n_cols]
-                       * pow2 ((j == fim_index[i])
-                              - deltaj
-                              + v[0] * deltaj));
-           }
+        double accum;
 
-       ase[2] = sqrt (accum - pt->total * v[0]) / (pt->total - cm);
+        accum = 0.;
+       for (i = 0; i < pt->n_rows; i++)
+          if (cm_index == fim_index[i])
+            accum += fim[i];
+        ase[2] = sqrt ((pt->total - sum_fim) * (sum_fim + cm - 2. * accum)
+                       / pow3 (pt->total - cm));
       }
 
       /* ASE0 for Y given PT. */
@@ -2833,19 +2825,14 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       /* ASE1 for PT given Y. */
       {
-       double accum;
+        double accum;
 
-       for (accum = 0., i = 0; i < pt->n_rows; i++)
-         for (j = 0; j < pt->n_cols; j++)
-           {
-             const int deltaj = i == rm_index;
-             accum += (pt->mat[j + i * pt->n_cols]
-                       * pow2 ((i == fmj_index[j])
-                              - deltaj
-                              + v[0] * deltaj));
-           }
-
-       ase[1] = sqrt (accum - pt->total * v[0]) / (pt->total - rm);
+        accum = 0.;
+       for (j = 0; j < pt->n_cols; j++)
+          if (rm_index == fmj_index[j])
+            accum += fmj[j];
+        ase[1] = sqrt ((pt->total - sum_fmj) * (sum_fmj + rm - 2. * accum)
+                       / pow3 (pt->total - rm));
       }
 
       /* ASE0 for PT given Y. */
@@ -2874,15 +2861,19 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
                         * pow2 (temp0 + (v[0] - 1.) * temp1));
            }
        ase[0] = sqrt (accum1 - 4. * pt->total * v[0] * v[0]) / (2. * pt->total - rm - cm);
-       t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / pt->total))
+       t[0] = v[0] / (sqrt (accum0 - pow2 (sum_fim + sum_fmj - cm - rm) / pt->total)
                       / (2. * pt->total - rm - cm));
       }
 
+      for (i = 0; i < 3; i++)
+        sig[i] = 2 * gsl_cdf_ugaussian_Q (t[i]);
+
       free (fim);
       free (fim_index);
       free (fmj);
       free (fmj_index);
 
+      /* Tau. */
       {
        double sum_fij2_ri, sum_fij2_ci;
        double sum_ri2, sum_cj2;
@@ -2951,8 +2942,7 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
       ase[5] = (2. / (pt->total * pow2 (UX + UY))) * sqrt (ase1_sym);
-      t[5] = v[5] / ((2. / (pt->total * (UX + UY)))
-                    * sqrt (P - pow2 (UX + UY - UXY) / pt->total));
+      t[5] = SYSMIS;
 
       v[6] = (UX + UY - UXY) / UX;
       ase[6] = sqrt (ase1_xy) / (pt->total * UX * UX);
@@ -2982,6 +2972,7 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
               v[8 + i] = somers_d_v[i];
               ase[8 + i] = somers_d_ase[i];
               t[8 + i] = somers_d_t[i];
+              sig[8 + i] = 2 * gsl_cdf_ugaussian_Q (fabs (somers_d_t[i]));
             }
         }
     }