CROSSTABS: Calculate significance for lambda and Somers' d.

[pspp] / src / language / stats / crosstabs.q

diff --git a/src/language/stats/crosstabs.q b/src/language/stats/crosstabs.q
index e3f5aee1d573ec8dc105d085538cc677090fcd26..aa7e2457ee24ae9d41d40bd989da985b2c822fe3 100644 (file)
--- a/src/language/stats/crosstabs.q
+++ b/src/language/stats/crosstabs.q
@@ -16,11 +16,10 @@
 
 /* FIXME:
 
-   - How to calculate significance of symmetric and directional measures?
-   - How to calculate ASE for asymmetric lambda?
-   - How to calculate ASE for symmetric Somers' d?
+   - 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?
-   - Approx. T of uncertainty coefficient is wrong.
+   - How to calculate approx. T of symmetric uncertainty coefficient?
 
 */
 
@@ -1975,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
@@ -2042,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);
@@ -2087,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);
     }
 
@@ -2721,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. */
@@ -2800,8 +2800,17 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
       v[1] = (sum_fmj - rm) / (pt->total - rm);
       v[2] = (sum_fim - cm) / (pt->total - cm);
 
-      /* XXX We don't have a working formula for ASE1. */
-      ase[2] = SYSMIS;
+      /* ASE1 for Y given PT. */
+      {
+        double accum;
+
+        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. */
       {
@@ -2814,8 +2823,17 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
        t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / pt->total) / (pt->total - cm));
       }
 
-      /* XXX We don't have a working formula for ASE1. */
-      ase[1] = SYSMIS;
+      /* ASE1 for PT given Y. */
+      {
+        double accum;
+
+        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. */
       {
@@ -2847,11 +2865,15 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
                       / (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;
@@ -2920,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);
@@ -2951,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]));
             }
         }
     }