X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Ffactor.c;h=a1cd333389c7276bc5a6ed7c37b776f2daf583a1;hb=445beb7fedbdee5ddaaa481e9ce77f2e2a5a2d70;hp=29e30c16b4b022b3eb26858507562dedc90f6d0e;hpb=dfd1972f7bcb550a4fc3b05dbe7e71d12334b0a7;p=pspp diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c index 29e30c16b4..a1cd333389 100644 --- a/src/language/stats/factor.c +++ b/src/language/stats/factor.c @@ -1,5 +1,6 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2009 Free Software Foundation, Inc. + Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, + 2016, 2017 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -16,38 +17,38 @@ #include - #include #include #include -#include -#include +#include +#include #include +#include + +#include "data/any-reader.h" +#include "data/casegrouper.h" +#include "data/casereader.h" +#include "data/casewriter.h" +#include "data/dataset.h" +#include "data/dictionary.h" +#include "data/format.h" +#include "data/subcase.h" +#include "language/command.h" +#include "language/lexer/lexer.h" +#include "language/lexer/value-parser.h" +#include "language/lexer/variable-parser.h" +#include "language/data-io/file-handle.h" +#include "language/data-io/matrix-reader.h" +#include "libpspp/cast.h" +#include "libpspp/message.h" +#include "libpspp/misc.h" +#include "math/correlation.h" +#include "math/covariance.h" +#include "math/moments.h" +#include "output/chart-item.h" +#include "output/charts/scree.h" +#include "output/tab.h" -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include -#include #include "gettext.h" #define _(msgid) gettext (msgid) @@ -91,12 +92,111 @@ enum print_opts PRINT_EXTRACTION = 0x0100, PRINT_INITIAL = 0x0200, PRINT_KMO = 0x0400, - PRINT_REPR = 0x0800, + PRINT_REPR = 0x0800, PRINT_FSCORE = 0x1000 }; +enum rotation_type + { + ROT_VARIMAX = 0, + ROT_EQUAMAX, + ROT_QUARTIMAX, + ROT_PROMAX, + ROT_NONE + }; + +typedef void (*rotation_coefficients) (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ); + + +static void +varimax_coefficients (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ) +{ + *x = d - 2 * a * b / loadings->size1; + *y = c - (a * a - b * b) / loadings->size1; +} + +static void +equamax_coefficients (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ) +{ + *x = d - loadings->size2 * a * b / loadings->size1; + *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1); +} + +static void +quartimax_coefficients (double *x, double *y, + double a UNUSED, double b UNUSED, double c, double d, + const gsl_matrix *loadings UNUSED) +{ + *x = d ; + *y = c ; +} + +static const rotation_coefficients rotation_coeff[] = { + varimax_coefficients, + equamax_coefficients, + quartimax_coefficients, + varimax_coefficients /* PROMAX is identical to VARIMAX */ +}; + + +/* return diag (C'C) ^ {-0.5} */ +static gsl_matrix * +diag_rcp_sqrt (const gsl_matrix *C) +{ + int j; + gsl_matrix *d = gsl_matrix_calloc (C->size1, C->size2); + gsl_matrix *r = gsl_matrix_calloc (C->size1, C->size2); + + assert (C->size1 == C->size2); + + gsl_linalg_matmult_mod (C, GSL_LINALG_MOD_TRANSPOSE, + C, GSL_LINALG_MOD_NONE, + d); + + for (j = 0 ; j < d->size2; ++j) + { + double e = gsl_matrix_get (d, j, j); + e = 1.0 / sqrt (e); + gsl_matrix_set (r, j, j, e); + } + + gsl_matrix_free (d); + + return r; +} + + + +/* return diag ((C'C)^-1) ^ {-0.5} */ +static gsl_matrix * +diag_rcp_inv_sqrt (const gsl_matrix *CCinv) +{ + int j; + gsl_matrix *r = gsl_matrix_calloc (CCinv->size1, CCinv->size2); + + assert (CCinv->size1 == CCinv->size2); + + for (j = 0 ; j < CCinv->size2; ++j) + { + double e = gsl_matrix_get (CCinv, j, j); + e = 1.0 / sqrt (e); + gsl_matrix_set (r, j, j, e); + } + + return r; +} + -struct cmd_factor + + + +struct cmd_factor { size_t n_vars; const struct variable **vars; @@ -109,25 +209,30 @@ struct cmd_factor enum print_opts print; enum extraction_method extraction; enum plot_opts plot; + enum rotation_type rotation; + int rotation_iterations; + int promax_power; /* Extraction Criteria */ int n_factors; double min_eigen; double econverge; - int iterations; + int extraction_iterations; + + double rconverge; /* Format */ double blank; bool sort; }; + struct idata { /* Intermediate values used in calculation */ + struct matrix_material mm; - const gsl_matrix *corr ; /* The correlation matrix */ - const gsl_matrix *cov ; /* The covariance matrix */ - const gsl_matrix *n ; /* Matrix of number of samples */ + gsl_matrix *analysis_matrix; /* A pointer to either mm.corr or mm.cov */ gsl_vector *eval ; /* The eigenvalues */ gsl_matrix *evec ; /* The eigenvectors */ @@ -135,6 +240,12 @@ struct idata int n_extractions; gsl_vector *msr ; /* Multiple Squared Regressions */ + + double detR; /* The determinant of the correlation matrix */ + + gsl_matrix *ai_cov; /* The anti-image covariance matrix */ + gsl_matrix *ai_cor; /* The anti-image correlation matrix */ + struct covariance *cvm; }; static struct idata * @@ -157,11 +268,108 @@ idata_free (struct idata *id) gsl_vector_free (id->msr); gsl_vector_free (id->eval); gsl_matrix_free (id->evec); + gsl_matrix_free (id->ai_cov); + gsl_matrix_free (id->ai_cor); free (id); } +/* Return the sum of squares of all the elements in row J excluding column J */ +static double +ssq_row_od_n (const gsl_matrix *m, int j) +{ + int i; + double ss = 0; + assert (m->size1 == m->size2); + + assert (j < m->size1); + + for (i = 0; i < m->size1; ++i) + { + if (i == j ) continue; + ss += pow2 (gsl_matrix_get (m, i, j)); + } + + return ss; +} + +/* Return the sum of all the elements excluding row N */ +static double +ssq_od_n (const gsl_matrix *m, int n) +{ + int i, j; + double ss = 0; + assert (m->size1 == m->size2); + + assert (n < m->size1); + + for (i = 0; i < m->size1; ++i) + { + if (i == n ) continue; + for (j = 0; j < m->size2; ++j) + { + ss += pow2 (gsl_matrix_get (m, i, j)); + } + } + + return ss; +} + + +static gsl_matrix * +anti_image_corr (const gsl_matrix *m, const struct idata *idata) +{ + int i, j; + gsl_matrix *a; + assert (m->size1 == m->size2); + + a = gsl_matrix_alloc (m->size1, m->size2); + + for (i = 0; i < m->size1; ++i) + { + for (j = 0; j < m->size2; ++j) + { + double *p = gsl_matrix_ptr (a, i, j); + *p = gsl_matrix_get (m, i, j); + *p /= sqrt (gsl_matrix_get (m, i, i) * + gsl_matrix_get (m, j, j)); + } + } + + for (i = 0; i < m->size1; ++i) + { + double r = ssq_row_od_n (idata->mm.corr, i); + double u = ssq_row_od_n (a, i); + gsl_matrix_set (a, i, i, r / (r + u)); + } + + return a; +} + +static gsl_matrix * +anti_image_cov (const gsl_matrix *m) +{ + int i, j; + gsl_matrix *a; + assert (m->size1 == m->size2); + + a = gsl_matrix_alloc (m->size1, m->size2); + for (i = 0; i < m->size1; ++i) + { + for (j = 0; j < m->size2; ++j) + { + double *p = gsl_matrix_ptr (a, i, j); + *p = gsl_matrix_get (m, i, j); + *p /= gsl_matrix_get (m, i, i); + *p /= gsl_matrix_get (m, j, j); + } + } + + return a; +} + +#if 0 static void dump_matrix (const gsl_matrix *m) { @@ -175,7 +383,6 @@ dump_matrix (const gsl_matrix *m) } } - static void dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p) { @@ -200,13 +407,14 @@ dump_vector (const gsl_vector *v) } printf ("\n"); } +#endif -static int +static int n_extracted_factors (const struct cmd_factor *factor, struct idata *idata) { int i; - + /* If there is a cached value, then return that. */ if ( idata->n_extractions != 0) return idata->n_extractions; @@ -218,7 +426,7 @@ n_extracted_factors (const struct cmd_factor *factor, struct idata *idata) idata->n_extractions = factor->n_factors; goto finish; } - + /* Use the MIN_EIGEN setting. */ for (i = 0 ; i < idata->eval->size; ++i) { @@ -253,7 +461,7 @@ struct smr_workspace { /* Copy of the subject */ gsl_matrix *m; - + gsl_matrix *inverse; gsl_permutation *perm; @@ -266,7 +474,7 @@ struct smr_workspace static struct smr_workspace *ws_create (const gsl_matrix *input) { struct smr_workspace *ws = xmalloc (sizeof (*ws)); - + ws->m = gsl_matrix_alloc (input->size1, input->size2); ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1); ws->perm = gsl_permutation_alloc (input->size1 - 1); @@ -289,13 +497,13 @@ ws_destroy (struct smr_workspace *ws) } -/* +/* Return the square of the regression coefficient for VAR regressed against all other variables. */ static double squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws) { - /* For an explanation of what this is doing, see + /* For an explanation of what this is doing, see http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm */ @@ -307,7 +515,7 @@ squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_worksp gsl_matrix_swap_rows (ws->m, 0, var); gsl_matrix_swap_columns (ws->m, 0, var); - rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1); + rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1); gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum); @@ -356,7 +564,7 @@ factor_matrix_workspace_alloc (size_t n, size_t nf) ws->eval = gsl_vector_alloc (n); ws->evec = gsl_matrix_alloc (n, n); ws->r = gsl_matrix_alloc (n, n); - + return ws; } @@ -391,10 +599,10 @@ perm_shift_apply (gsl_permutation *target, const gsl_permutation *p, } -/* +/* Indirectly sort the rows of matrix INPUT, storing the sort order in PERM. The sort criteria are as follows: - + Rows are sorted on the first column, until the absolute value of an element in a subsequent column is greater than that of the first column. Thereafter, rows will be sorted on the second column, @@ -427,7 +635,7 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm) } } - while (column_n < m && row_n < n) + while (column_n < m && row_n < n) { gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n); gsl_sort_vector_index (p, &columni.vector); @@ -436,7 +644,7 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm) { gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]); size_t maxindex = gsl_vector_max_index (&row.vector); - + if ( maxindex > column_n ) break; @@ -454,11 +662,342 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm) gsl_permutation_free (p); gsl_matrix_free (mat); - + assert ( 0 == gsl_permutation_valid (perm)); /* We want the biggest value to be first */ - gsl_permutation_reverse (perm); + gsl_permutation_reverse (perm); +} + + +static void +drot_go (double phi, double *l0, double *l1) +{ + double r0 = cos (phi) * *l0 + sin (phi) * *l1; + double r1 = - sin (phi) * *l0 + cos (phi) * *l1; + + *l0 = r0; + *l1 = r1; +} + + +static gsl_matrix * +clone_matrix (const gsl_matrix *m) +{ + int j, k; + gsl_matrix *c = gsl_matrix_calloc (m->size1, m->size2); + + for (j = 0 ; j < c->size1; ++j) + { + for (k = 0 ; k < c->size2; ++k) + { + const double *v = gsl_matrix_const_ptr (m, j, k); + gsl_matrix_set (c, j, k, *v); + } + } + + return c; +} + + +static double +initial_sv (const gsl_matrix *fm) +{ + int j, k; + + double sv = 0.0; + for (j = 0 ; j < fm->size2; ++j) + { + double l4s = 0; + double l2s = 0; + + for (k = j + 1 ; k < fm->size2; ++k) + { + double lambda = gsl_matrix_get (fm, k, j); + double lambda_sq = lambda * lambda; + double lambda_4 = lambda_sq * lambda_sq; + + l4s += lambda_4; + l2s += lambda_sq; + } + sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 ); + } + return sv; +} + +static void +rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, + const gsl_vector *communalities, + gsl_matrix *result, + gsl_vector *rotated_loadings, + gsl_matrix *pattern_matrix, + gsl_matrix *factor_correlation_matrix + ) +{ + int j, k; + int i; + double prev_sv; + + /* First get a normalised version of UNROT */ + gsl_matrix *normalised = gsl_matrix_calloc (unrot->size1, unrot->size2); + gsl_matrix *h_sqrt = gsl_matrix_calloc (communalities->size, communalities->size); + gsl_matrix *h_sqrt_inv ; + + /* H is the diagonal matrix containing the absolute values of the communalities */ + for (i = 0 ; i < communalities->size ; ++i) + { + double *ptr = gsl_matrix_ptr (h_sqrt, i, i); + *ptr = fabs (gsl_vector_get (communalities, i)); + } + + /* Take the square root of the communalities */ + gsl_linalg_cholesky_decomp (h_sqrt); + + + /* Save a copy of h_sqrt and invert it */ + h_sqrt_inv = clone_matrix (h_sqrt); + gsl_linalg_cholesky_decomp (h_sqrt_inv); + gsl_linalg_cholesky_invert (h_sqrt_inv); + + /* normalised vertion is H^{1/2} x UNROT */ + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, h_sqrt_inv, unrot, 0.0, normalised); + + gsl_matrix_free (h_sqrt_inv); + + + /* Now perform the rotation iterations */ + + prev_sv = initial_sv (normalised); + for (i = 0 ; i < cf->rotation_iterations ; ++i) + { + double sv = 0.0; + for (j = 0 ; j < normalised->size2; ++j) + { + /* These variables relate to the convergence criterium */ + double l4s = 0; + double l2s = 0; + + for (k = j + 1 ; k < normalised->size2; ++k) + { + int p; + double a = 0.0; + double b = 0.0; + double c = 0.0; + double d = 0.0; + double x, y; + double phi; + + for (p = 0; p < normalised->size1; ++p) + { + double jv = gsl_matrix_get (normalised, p, j); + double kv = gsl_matrix_get (normalised, p, k); + + double u = jv * jv - kv * kv; + double v = 2 * jv * kv; + a += u; + b += v; + c += u * u - v * v; + d += 2 * u * v; + } + + rotation_coeff [cf->rotation] (&x, &y, a, b, c, d, normalised); + + phi = atan2 (x, y) / 4.0 ; + + /* Don't bother rotating if the angle is small */ + if ( fabs (sin (phi) ) <= pow (10.0, -15.0)) + continue; + + for (p = 0; p < normalised->size1; ++p) + { + double *lambda0 = gsl_matrix_ptr (normalised, p, j); + double *lambda1 = gsl_matrix_ptr (normalised, p, k); + drot_go (phi, lambda0, lambda1); + } + + /* Calculate the convergence criterium */ + { + double lambda = gsl_matrix_get (normalised, k, j); + double lambda_sq = lambda * lambda; + double lambda_4 = lambda_sq * lambda_sq; + + l4s += lambda_4; + l2s += lambda_sq; + } + } + sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 ); + } + + if ( fabs (sv - prev_sv) <= cf->rconverge) + break; + + prev_sv = sv; + } + + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, + h_sqrt, normalised, 0.0, result); + + gsl_matrix_free (h_sqrt); + gsl_matrix_free (normalised); + + if (cf->rotation == ROT_PROMAX) + { + /* general purpose m by m matrix, where m is the number of factors */ + gsl_matrix *mm1 = gsl_matrix_calloc (unrot->size2, unrot->size2); + gsl_matrix *mm2 = gsl_matrix_calloc (unrot->size2, unrot->size2); + + /* general purpose m by p matrix, where p is the number of variables */ + gsl_matrix *mp1 = gsl_matrix_calloc (unrot->size2, unrot->size1); + + gsl_matrix *pm1 = gsl_matrix_calloc (unrot->size1, unrot->size2); + + gsl_permutation *perm = gsl_permutation_alloc (unrot->size2); + + int signum; + + int i, j; + + /* The following variables follow the notation by SPSS Statistical Algorithms + page 342 */ + gsl_matrix *L = gsl_matrix_calloc (unrot->size2, unrot->size2); + gsl_matrix *P = clone_matrix (result); + gsl_matrix *D ; + gsl_matrix *Q ; + + + /* Vector of length p containing (indexed by i) + \Sum^m_j {\lambda^2_{ij}} */ + gsl_vector *rssq = gsl_vector_calloc (unrot->size1); + + for (i = 0; i < P->size1; ++i) + { + double sum = 0; + for (j = 0; j < P->size2; ++j) + { + sum += gsl_matrix_get (result, i, j) + * gsl_matrix_get (result, i, j); + + } + + gsl_vector_set (rssq, i, sqrt (sum)); + } + + for (i = 0; i < P->size1; ++i) + { + for (j = 0; j < P->size2; ++j) + { + double l = gsl_matrix_get (result, i, j); + double r = gsl_vector_get (rssq, i); + gsl_matrix_set (P, i, j, pow (fabs (l / r), cf->promax_power + 1) * r / l); + } + } + + gsl_vector_free (rssq); + + gsl_linalg_matmult_mod (result, + GSL_LINALG_MOD_TRANSPOSE, + result, + GSL_LINALG_MOD_NONE, + mm1); + + gsl_linalg_LU_decomp (mm1, perm, &signum); + gsl_linalg_LU_invert (mm1, perm, mm2); + + gsl_linalg_matmult_mod (mm2, GSL_LINALG_MOD_NONE, + result, GSL_LINALG_MOD_TRANSPOSE, + mp1); + + gsl_linalg_matmult_mod (mp1, GSL_LINALG_MOD_NONE, + P, GSL_LINALG_MOD_NONE, + L); + + D = diag_rcp_sqrt (L); + Q = gsl_matrix_calloc (unrot->size2, unrot->size2); + + gsl_linalg_matmult_mod (L, GSL_LINALG_MOD_NONE, + D, GSL_LINALG_MOD_NONE, + Q); + + gsl_matrix *QQinv = gsl_matrix_calloc (unrot->size2, unrot->size2); + + gsl_linalg_matmult_mod (Q, GSL_LINALG_MOD_TRANSPOSE, + Q, GSL_LINALG_MOD_NONE, + QQinv); + + gsl_linalg_cholesky_decomp (QQinv); + gsl_linalg_cholesky_invert (QQinv); + + + gsl_matrix *C = diag_rcp_inv_sqrt (QQinv); + gsl_matrix *Cinv = clone_matrix (C); + + gsl_linalg_cholesky_decomp (Cinv); + gsl_linalg_cholesky_invert (Cinv); + + + gsl_linalg_matmult_mod (result, GSL_LINALG_MOD_NONE, + Q, GSL_LINALG_MOD_NONE, + pm1); + + gsl_linalg_matmult_mod (pm1, GSL_LINALG_MOD_NONE, + Cinv, GSL_LINALG_MOD_NONE, + pattern_matrix); + + + gsl_linalg_matmult_mod (C, GSL_LINALG_MOD_NONE, + QQinv, GSL_LINALG_MOD_NONE, + mm1); + + gsl_linalg_matmult_mod (mm1, GSL_LINALG_MOD_NONE, + C, GSL_LINALG_MOD_TRANSPOSE, + factor_correlation_matrix); + + gsl_linalg_matmult_mod (pattern_matrix, GSL_LINALG_MOD_NONE, + factor_correlation_matrix, GSL_LINALG_MOD_NONE, + pm1); + + gsl_matrix_memcpy (result, pm1); + + + gsl_matrix_free (QQinv); + gsl_matrix_free (C); + gsl_matrix_free (Cinv); + + gsl_matrix_free (D); + gsl_matrix_free (Q); + gsl_matrix_free (L); + gsl_matrix_free (P); + + gsl_permutation_free (perm); + + gsl_matrix_free (mm1); + gsl_matrix_free (mm2); + gsl_matrix_free (mp1); + gsl_matrix_free (pm1); + } + + + /* reflect negative sums and populate the rotated loadings vector*/ + for (i = 0 ; i < result->size2; ++i) + { + double ssq = 0.0; + double sum = 0.0; + for (j = 0 ; j < result->size1; ++j) + { + double s = gsl_matrix_get (result, j, i); + ssq += s * s; + sum += s; + } + + gsl_vector_set (rotated_loadings, i, ssq); + + if ( sum < 0 ) + for (j = 0 ; j < result->size1; ++j) + { + double *lambda = gsl_matrix_ptr (result, j, i); + *lambda = - *lambda; + } + } } @@ -468,7 +1007,8 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm) WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init. */ static void -iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, struct factor_matrix_workspace *ws) +iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, + struct factor_matrix_workspace *ws) { size_t i; gsl_matrix_view mv ; @@ -502,8 +1042,7 @@ iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matri /* Take the square root of gamma */ gsl_linalg_cholesky_decomp (ws->gamma); - gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, - 1.0, &mv.matrix, ws->gamma, 0.0, factors); + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, &mv.matrix, ws->gamma, 0.0, factors); for (i = 0 ; i < r->size1 ; ++i) { @@ -516,14 +1055,18 @@ iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matri static bool run_factor (struct dataset *ds, const struct cmd_factor *factor); +static void do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata); + + int cmd_factor (struct lexer *lexer, struct dataset *ds) { - bool extraction_seen = false; - const struct dictionary *dict = dataset_dict (ds); - + struct dictionary *dict = NULL; + int n_iterations = 25; struct cmd_factor factor; + factor.n_vars = 0; + factor.vars = NULL; factor.method = METHOD_CORR; factor.missing_type = MISS_LISTWISE; factor.exclude = MV_ANY; @@ -531,38 +1074,118 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) factor.extraction = EXTRACTION_PC; factor.n_factors = 0; factor.min_eigen = SYSMIS; - factor.iterations = 25; + factor.extraction_iterations = 25; + factor.rotation_iterations = 25; factor.econverge = 0.001; + factor.blank = 0; factor.sort = false; factor.plot = 0; + factor.rotation = ROT_VARIMAX; + factor.wv = NULL; - factor.wv = dict_get_weight (dict); - - lex_match (lexer, '/'); - - if (!lex_force_match_id (lexer, "VARIABLES")) - { - goto error; - } + factor.rconverge = 0.0001; - lex_match (lexer, '='); + lex_match (lexer, T_SLASH); - if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars, - PV_NO_DUPLICATE | PV_NUMERIC)) - goto error; + struct matrix_reader *mr = NULL; + struct casereader *matrix_reader = NULL; - while (lex_token (lexer) != '.') + if (lex_match_id (lexer, "VARIABLES")) { - lex_match (lexer, '/'); + lex_match (lexer, T_EQUALS); + dict = dataset_dict (ds); + factor.wv = dict_get_weight (dict); - if (lex_match_id (lexer, "PLOT")) + if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars, + PV_NO_DUPLICATE | PV_NUMERIC)) + goto error; + } + else if (lex_match_id (lexer, "MATRIX")) + { + lex_match (lexer, T_EQUALS); + if (! lex_force_match_id (lexer, "IN")) + goto error; + if (!lex_force_match (lexer, T_LPAREN)) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') - { - if (lex_match_id (lexer, "EIGEN")) - { + goto error; + } + if (lex_match_id (lexer, "CORR")) + { + } + else if (lex_match_id (lexer, "COV")) + { + } + else + { + lex_error (lexer, _("Matrix input for %s must be either COV or CORR"), "FACTOR"); + goto error; + } + if (! lex_force_match (lexer, T_EQUALS)) + goto error; + if (lex_match (lexer, T_ASTERISK)) + { + dict = dataset_dict (ds); + matrix_reader = casereader_clone (dataset_source (ds)); + } + else + { + struct file_handle *fh = fh_parse (lexer, FH_REF_FILE, NULL); + if (fh == NULL) + goto error; + + matrix_reader + = any_reader_open_and_decode (fh, NULL, &dict, NULL); + + if (! (matrix_reader && dict)) + { + goto error; + } + } + + if (! lex_force_match (lexer, T_RPAREN)) + goto error; + + mr = create_matrix_reader_from_case_reader (dict, matrix_reader, + &factor.vars, &factor.n_vars); + } + else + { + goto error; + } + + while (lex_token (lexer) != T_ENDCMD) + { + lex_match (lexer, T_SLASH); + + if (lex_match_id (lexer, "ANALYSIS")) + { + struct const_var_set *vs; + const struct variable **vars; + size_t n_vars; + bool ok; + + lex_match (lexer, T_EQUALS); + + vs = const_var_set_create_from_array (factor.vars, factor.n_vars); + ok = parse_const_var_set_vars (lexer, vs, &vars, &n_vars, + PV_NO_DUPLICATE | PV_NUMERIC); + const_var_set_destroy (vs); + + if (!ok) + goto error; + + free (factor.vars); + factor.vars = vars; + factor.n_vars = n_vars; + } + else if (lex_match_id (lexer, "PLOT")) + { + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) + { + if (lex_match_id (lexer, "EIGEN")) + { factor.plot |= PLOT_SCREE; } #if FACTOR_FULLY_IMPLEMENTED @@ -579,8 +1202,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "METHOD")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "COVARIANCE")) { @@ -597,17 +1220,40 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } } -#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "ROTATION")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { - if (lex_match_id (lexer, "VARIMAX")) + /* VARIMAX and DEFAULT are defaults */ + if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT")) { + factor.rotation = ROT_VARIMAX; } - else if (lex_match_id (lexer, "DEFAULT")) + else if (lex_match_id (lexer, "EQUAMAX")) + { + factor.rotation = ROT_EQUAMAX; + } + else if (lex_match_id (lexer, "QUARTIMAX")) + { + factor.rotation = ROT_QUARTIMAX; + } + else if (lex_match_id (lexer, "PROMAX")) + { + factor.promax_power = 5; + if (lex_match (lexer, T_LPAREN) + && lex_force_int (lexer)) + { + factor.promax_power = lex_integer (lexer); + lex_get (lexer); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; + } + factor.rotation = ROT_PROMAX; + } + else if (lex_match_id (lexer, "NOROTATE")) { + factor.rotation = ROT_NONE; } else { @@ -615,58 +1261,73 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) goto error; } } + factor.rotation_iterations = n_iterations; } -#endif else if (lex_match_id (lexer, "CRITERIA")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "FACTORS")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN) + && lex_force_int (lexer)) { - lex_force_int (lexer); factor.n_factors = lex_integer (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; } } else if (lex_match_id (lexer, "MINEIGEN")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN) + && lex_force_num (lexer)) { - lex_force_num (lexer); factor.min_eigen = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; } } else if (lex_match_id (lexer, "ECONVERGE")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN) + && lex_force_num (lexer)) { - lex_force_num (lexer); factor.econverge = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; } } + else if (lex_match_id (lexer, "RCONVERGE")) + { + if (lex_force_match (lexer, T_LPAREN) + && lex_force_num (lexer)) + { + factor.rconverge = lex_number (lexer); + lex_get (lexer); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; + } + } else if (lex_match_id (lexer, "ITERATE")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN) + && lex_force_int (lexer)) { - lex_force_int (lexer); - factor.iterations = lex_integer (lexer); + n_iterations = lex_integer (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; } } else if (lex_match_id (lexer, "DEFAULT")) { factor.n_factors = 0; factor.min_eigen = 1; - factor.iterations = 25; + n_iterations = 25; } else { @@ -677,9 +1338,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "EXTRACTION")) { - extraction_seen = true; - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "PAF")) { @@ -703,11 +1363,12 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) goto error; } } + factor.extraction_iterations = n_iterations; } else if (lex_match_id (lexer, "FORMAT")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "SORT")) { @@ -715,12 +1376,13 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "BLANK")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN) + && lex_force_num (lexer)) { - lex_force_num (lexer); factor.blank = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + if (! lex_force_match (lexer, T_RPAREN)) + goto error; } } else if (lex_match_id (lexer, "DEFAULT")) @@ -738,8 +1400,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) else if (lex_match_id (lexer, "PRINT")) { factor.print = 0; - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "UNIVARIATE")) { @@ -753,10 +1415,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) else if (lex_match_id (lexer, "INV")) { } +#endif else if (lex_match_id (lexer, "AIC")) { + factor.print |= PRINT_AIC; } -#endif else if (lex_match_id (lexer, "SIG")) { factor.print |= PRINT_SIG; @@ -765,11 +1428,10 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) { factor.print |= PRINT_CORRELATION; } -#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "COVARIANCE")) { + factor.print |= PRINT_COVARIANCE; } -#endif else if (lex_match_id (lexer, "ROTATION")) { factor.print |= PRINT_ROTATION; @@ -782,10 +1444,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) { factor.print |= PRINT_INITIAL; } -#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "KMO")) { + factor.print |= PRINT_KMO; } +#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "REPR")) { } @@ -812,8 +1475,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "MISSING")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "INCLUDE")) { @@ -849,13 +1512,46 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } - if ( ! run_factor (ds, &factor)) - goto error; + if ( factor.rotation == ROT_NONE ) + factor.print &= ~PRINT_ROTATION; + + if (factor.n_vars < 2) + msg (MW, _("Factor analysis on a single variable is not useful.")); + + if (factor.n_vars < 1) + { + msg (ME, _("Factor analysis without variables is not possible.")); + goto error; + } + + if (matrix_reader) + { + struct idata *id = idata_alloc (factor.n_vars); + + while (next_matrix_from_reader (&id->mm, mr, + factor.vars, factor.n_vars)) + { + do_factor_by_matrix (&factor, id); + + gsl_matrix_free (id->mm.corr); + id->mm.corr = NULL; + gsl_matrix_free (id->mm.cov); + id->mm.cov = NULL; + } + + idata_free (id); + } + else + if ( ! run_factor (ds, &factor)) + goto error; + + destroy_matrix_reader (mr); free (factor.vars); return CMD_SUCCESS; error: + destroy_matrix_reader (mr); free (factor.vars); return CMD_FAILURE; } @@ -915,14 +1611,14 @@ the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_fa /* Return the communality of variable N, calculated to N_FACTORS */ static double -communality (struct idata *idata, int n, int n_factors) +communality (const struct idata *idata, int n, int n_factors) { return the_communality (idata->evec, idata->eval, n, n_factors); } static void -show_scree (const struct cmd_factor *f, struct idata *idata) +show_scree (const struct cmd_factor *f, const struct idata *idata) { struct scree *s; const char *label ; @@ -996,10 +1692,10 @@ show_communalities (const struct cmd_factor * factor, tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i])); if (factor->print & PRINT_INITIAL) - tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL); + tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL, RC_OTHER); if (factor->print & PRINT_EXTRACTION) - tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL); + tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL, RC_OTHER); } tab_submit (t); @@ -1007,10 +1703,11 @@ show_communalities (const struct cmd_factor * factor, static void -show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const gsl_matrix *fm) +show_factor_matrix (const struct cmd_factor *factor, const struct idata *idata, const char *title, const gsl_matrix *fm) { int i; - const int n_factors = n_extracted_factors (factor, idata); + + const int n_factors = idata->n_extractions; const int heading_columns = 1; const int heading_rows = 2; @@ -1020,10 +1717,14 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const struct tab_table *t = tab_create (nc, nr); + /* if ( factor->extraction == EXTRACTION_PC ) tab_title (t, _("Component Matrix")); - else + else tab_title (t, _("Factor Matrix")); + */ + + tab_title (t, "%s", title); tab_headers (t, heading_columns, 0, heading_rows, 0); @@ -1084,7 +1785,7 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const if ( fabs (x) < factor->blank) continue; - tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL); + tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL, RC_OTHER); } } @@ -1095,9 +1796,11 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const static void -show_explained_variance (const struct cmd_factor * factor, struct idata *idata, +show_explained_variance (const struct cmd_factor * factor, + const struct idata *idata, const gsl_vector *initial_eigenvalues, - const gsl_vector *extracted_eigenvalues) + const gsl_vector *extracted_eigenvalues, + const gsl_vector *rotated_loadings) { size_t i; int c = 0; @@ -1113,6 +1816,8 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, double e_total = 0.0; double e_cum = 0.0; + double r_cum = 0.0; + int nc = heading_columns; if (factor->print & PRINT_EXTRACTION) @@ -1122,7 +1827,9 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, nc += 3; if (factor->print & PRINT_ROTATION) - nc += 3; + { + nc += factor->rotation == ROT_PROMAX ? 1 : 3; + } /* No point having a table with only headings */ if ( nc <= heading_columns) @@ -1174,17 +1881,23 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, if (factor->print & PRINT_ROTATION) { - tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings")); - c += 3; + const int width = factor->rotation == ROT_PROMAX ? 0 : 2; + tab_joint_text (t, c, 0, c + width, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings")); + c += width + 1; } - for (i = 0; i < (nc - heading_columns) / 3 ; ++i) + for (i = 0; i < (nc - heading_columns + 2) / 3 ; ++i) { tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total")); - tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance")); - tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %")); tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1); + + if (i == 2 && factor->rotation == ROT_PROMAX) + continue; + + /* xgettext:no-c-format */ + tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance")); + tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %")); } for (i = 0 ; i < initial_eigenvalues->size; ++i) @@ -1199,7 +1912,6 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, e_total = i_total; } - for (i = 0 ; i < factor->n_vars; ++i) { const double i_lambda = gsl_vector_get (initial_eigenvalues, i); @@ -1210,7 +1922,7 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, c = 0; - tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%d"), i + 1); + tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1); i_cum += i_percent; e_cum += e_percent; @@ -1218,26 +1930,183 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, /* Initial Eigenvalues */ if (factor->print & PRINT_INITIAL) { - tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL); - tab_double (t, c++, i + heading_rows, 0, i_percent, NULL); - tab_double (t, c++, i + heading_rows, 0, i_cum, NULL); + tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL, RC_OTHER); + tab_double (t, c++, i + heading_rows, 0, i_percent, NULL, RC_OTHER); + tab_double (t, c++, i + heading_rows, 0, i_cum, NULL, RC_OTHER); } + if (factor->print & PRINT_EXTRACTION) { - if ( i < n_extracted_factors (factor, idata)) + if (i < idata->n_extractions) { /* Sums of squared loadings */ - tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL); - tab_double (t, c++, i + heading_rows, 0, e_percent, NULL); - tab_double (t, c++, i + heading_rows, 0, e_cum, NULL); + tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL, RC_OTHER); + tab_double (t, c++, i + heading_rows, 0, e_percent, NULL, RC_OTHER); + tab_double (t, c++, i + heading_rows, 0, e_cum, NULL, RC_OTHER); } } + + if (rotated_loadings != NULL) + { + const double r_lambda = gsl_vector_get (rotated_loadings, i); + double r_percent = 100.0 * r_lambda / e_total ; + + if (factor->print & PRINT_ROTATION) + { + if (i < idata->n_extractions) + { + r_cum += r_percent; + tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL, RC_OTHER); + if (factor->rotation != ROT_PROMAX) + { + tab_double (t, c++, i + heading_rows, 0, r_percent, NULL, RC_OTHER); + tab_double (t, c++, i + heading_rows, 0, r_cum, NULL, RC_OTHER); + } + } + } + } + } + + tab_submit (t); +} + + +static void +show_factor_correlation (const struct cmd_factor * factor, const gsl_matrix *fcm) +{ + size_t i, j; + const int heading_columns = 1; + const int heading_rows = 1; + const int nr = heading_rows + fcm->size2; + const int nc = heading_columns + fcm->size1; + struct tab_table *t = tab_create (nc, nr); + + tab_title (t, _("Factor Correlation Matrix")); + + tab_headers (t, heading_columns, 0, heading_rows, 0); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + /* Vertical lines */ + tab_box (t, + -1, -1, + -1, TAL_1, + heading_columns, 0, + nc - 1, nr - 1); + + tab_hline (t, TAL_1, 0, nc - 1, heading_rows); + tab_hline (t, TAL_1, 1, nc - 1, 1); + + tab_vline (t, TAL_2, heading_columns, 0, nr - 1); + + + if ( factor->extraction == EXTRACTION_PC) + tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("Component")); + else + tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("Factor")); + + for (i = 0 ; i < fcm->size1; ++i) + { + tab_text_format (t, heading_columns + i, 0, TAB_CENTER | TAT_TITLE, _("%zu"), i + 1); + } + + for (i = 0 ; i < fcm->size2; ++i) + { + tab_text_format (t, 0, heading_rows + i, TAB_CENTER | TAT_TITLE, _("%zu"), i + 1); + } + + + for (i = 0 ; i < fcm->size1; ++i) + { + for (j = 0 ; j < fcm->size2; ++j) + tab_double (t, heading_columns + j, heading_rows + i, 0, + gsl_matrix_get (fcm, i, j), NULL, RC_OTHER); } tab_submit (t); } +static void +show_aic (const struct cmd_factor *factor, const struct idata *idata) +{ + struct tab_table *t ; + size_t i; + + const int heading_rows = 1; + const int heading_columns = 2; + + const int nc = heading_columns + factor->n_vars; + const int nr = heading_rows + 2 * factor->n_vars; + + if ((factor->print & PRINT_AIC) == 0) + return; + + t = tab_create (nc, nr); + + tab_title (t, _("Anti-Image Matrices")); + + tab_hline (t, TAL_1, 0, nc - 1, heading_rows); + + tab_headers (t, heading_columns, 0, heading_rows, 0); + + tab_vline (t, TAL_2, 2, 0, nr - 1); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + /* Vertical lines */ + tab_box (t, + -1, -1, + -1, TAL_1, + heading_columns, 0, + nc - 1, nr - 1); + + + for (i = 0; i < factor->n_vars; ++i) + tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i])); + + tab_text (t, 0, heading_rows, TAT_TITLE, _("Anti-image Covariance")); + tab_hline (t, TAL_1, 0, nc - 1, heading_rows + factor->n_vars); + tab_text (t, 0, heading_rows + factor->n_vars, TAT_TITLE, _("Anti-image Correlation")); + + for (i = 0; i < factor->n_vars; ++i) + { + tab_text (t, 1, i + heading_rows, TAT_TITLE, + var_to_string (factor->vars[i])); + + tab_text (t, 1, factor->n_vars + i + heading_rows, TAT_TITLE, + var_to_string (factor->vars[i])); + } + + for (i = 0; i < factor->n_vars; ++i) + { + int j; + for (j = 0; j < factor->n_vars; ++j) + { + tab_double (t, heading_columns + i, heading_rows + j, 0, + gsl_matrix_get (idata->ai_cov, i, j), NULL, RC_OTHER); + } + + + for (j = 0; j < factor->n_vars; ++j) + { + tab_double (t, heading_columns + i, factor->n_vars + heading_rows + j, 0, + gsl_matrix_get (idata->ai_cor, i, j), NULL, RC_OTHER); + } + } + + tab_submit (t); +} static void show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata) @@ -1327,26 +2196,26 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id for (i = 0; i < factor->n_vars; ++i) { for (j = 0; j < factor->n_vars; ++j) - tab_double (t, heading_columns + i, y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL); + tab_double (t, heading_columns + j, y + i, 0, gsl_matrix_get (idata->mm.corr, i, j), NULL, RC_OTHER); } } if (factor->print & PRINT_SIG) { const double y = heading_rows + y_pos_sig * factor->n_vars; - tab_text (t, 0, y, TAT_TITLE, _("Sig. 1-tailed")); + tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)")); for (i = 0; i < factor->n_vars; ++i) { for (j = 0; j < factor->n_vars; ++j) { - double rho = gsl_matrix_get (idata->corr, i, j); - double w = gsl_matrix_get (idata->n, i, j); + double rho = gsl_matrix_get (idata->mm.corr, i, j); + double w = gsl_matrix_get (idata->mm.n, i, j); if (i == j) continue; - tab_double (t, heading_columns + i, y + j, 0, significance_of_correlation (rho, w), NULL); + tab_double (t, heading_columns + j, y + i, 0, significance_of_correlation (rho, w), NULL, RC_PVALUE); } } } @@ -1354,67 +2223,194 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id if (factor->print & PRINT_DETERMINANT) { - int sign = 0; - double det = 0.0; + tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant")); - const int size = idata->corr->size1; - gsl_permutation *p = gsl_permutation_calloc (size); - gsl_matrix *tmp = gsl_matrix_calloc (size, size); - gsl_matrix_memcpy (tmp, idata->corr); + tab_double (t, 1, nr, 0, idata->detR, NULL, RC_OTHER); + } - gsl_linalg_LU_decomp (tmp, p, &sign); - det = gsl_linalg_LU_det (tmp, sign); - gsl_permutation_free (p); - gsl_matrix_free (tmp); + tab_submit (t); +} +static void +show_covariance_matrix (const struct cmd_factor *factor, const struct idata *idata) +{ + struct tab_table *t ; + size_t i, j; + int y_pos_corr = -1; + int suffix_rows = 0; - tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant")); - tab_double (t, 1, nr, 0, det, NULL); + const int heading_rows = 1; + const int heading_columns = 1; + + int nc = heading_columns ; + int nr = heading_rows ; + int n_data_sets = 0; + + if (factor->print & PRINT_COVARIANCE) + { + y_pos_corr = n_data_sets; + n_data_sets++; + nc = heading_columns + factor->n_vars; + } + + nr += n_data_sets * factor->n_vars; + + /* If the table would contain only headings, don't bother rendering it */ + if (nr <= heading_rows && suffix_rows == 0) + return; + + t = tab_create (nc, nr + suffix_rows); + + tab_title (t, _("Covariance Matrix")); + + tab_hline (t, TAL_1, 0, nc - 1, heading_rows); + + if (nr > heading_rows) + { + tab_headers (t, heading_columns, 0, heading_rows, 0); + + tab_vline (t, TAL_2, heading_columns, 0, nr - 1); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + /* Vertical lines */ + tab_box (t, + -1, -1, + -1, TAL_1, + heading_columns, 0, + nc - 1, nr - 1); + + + for (i = 0; i < factor->n_vars; ++i) + tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i])); + + + for (i = 0 ; i < n_data_sets; ++i) + { + int y = heading_rows + i * factor->n_vars; + size_t v; + for (v = 0; v < factor->n_vars; ++v) + tab_text (t, heading_columns -1, y + v, TAT_TITLE, var_to_string (factor->vars[v])); + + tab_hline (t, TAL_1, 0, nc - 1, y); + } + + if (factor->print & PRINT_COVARIANCE) + { + const double y = heading_rows + y_pos_corr; + + for (i = 0; i < factor->n_vars; ++i) + { + for (j = 0; j < factor->n_vars; ++j) + tab_double (t, heading_columns + j, y + i, 0, gsl_matrix_get (idata->mm.cov, i, j), NULL, RC_OTHER); + } + } } tab_submit (t); } - static void do_factor (const struct cmd_factor *factor, struct casereader *r) { struct ccase *c; - const gsl_matrix *var_matrix; - const gsl_matrix *mean_matrix; - - const gsl_matrix *analysis_matrix; struct idata *idata = idata_alloc (factor->n_vars); - struct covariance *cov = covariance_create (factor->n_vars, factor->vars, - factor->wv, factor->exclude); + idata->cvm = covariance_1pass_create (factor->n_vars, factor->vars, + factor->wv, factor->exclude, true); for ( ; (c = casereader_read (r) ); case_unref (c)) { - covariance_accumulate (cov, c); + covariance_accumulate (idata->cvm, c); + } + + idata->mm.cov = covariance_calculate (idata->cvm); + + if (idata->mm.cov == NULL) + { + msg (MW, _("The dataset contains no complete observations. No analysis will be performed.")); + covariance_destroy (idata->cvm); + goto finish; } - idata->cov = covariance_calculate (cov); + idata->mm.var_matrix = covariance_moments (idata->cvm, MOMENT_VARIANCE); + idata->mm.mean_matrix = covariance_moments (idata->cvm, MOMENT_MEAN); + idata->mm.n = covariance_moments (idata->cvm, MOMENT_NONE); + + do_factor_by_matrix (factor, idata); - var_matrix = covariance_moments (cov, MOMENT_VARIANCE); - mean_matrix = covariance_moments (cov, MOMENT_MEAN); - idata->n = covariance_moments (cov, MOMENT_NONE); + finish: + gsl_matrix_free (idata->mm.corr); + gsl_matrix_free (idata->mm.cov); + + idata_free (idata); + casereader_destroy (r); +} - if ( factor->method == METHOD_CORR) +static void +do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata) +{ + if (!idata->mm.cov && !idata->mm.corr) { - idata->corr = correlation_from_covariance (idata->cov, var_matrix); - analysis_matrix = idata->corr; + msg (ME, _("The dataset has no complete covariance or correlation matrix.")); + return; } + + if (idata->mm.cov && !idata->mm.corr) + idata->mm.corr = correlation_from_covariance (idata->mm.cov, idata->mm.var_matrix); + if (idata->mm.corr && !idata->mm.cov) + idata->mm.cov = covariance_from_correlation (idata->mm.corr, idata->mm.var_matrix); + if (factor->method == METHOD_CORR) + idata->analysis_matrix = idata->mm.corr; else - analysis_matrix = idata->cov; + idata->analysis_matrix = idata->mm.cov; + + gsl_matrix *r_inv; + r_inv = clone_matrix (idata->mm.corr); + gsl_linalg_cholesky_decomp (r_inv); + gsl_linalg_cholesky_invert (r_inv); + + idata->ai_cov = anti_image_cov (r_inv); + idata->ai_cor = anti_image_corr (r_inv, idata); + + int i; + double sum_ssq_r = 0; + double sum_ssq_a = 0; + for (i = 0; i < r_inv->size1; ++i) + { + sum_ssq_r += ssq_od_n (r_inv, i); + sum_ssq_a += ssq_od_n (idata->ai_cov, i); + } + + gsl_matrix_free (r_inv); + + if (factor->print & PRINT_DETERMINANT + || factor->print & PRINT_KMO) + { + int sign = 0; + + const int size = idata->mm.corr->size1; + gsl_permutation *p = gsl_permutation_calloc (size); + gsl_matrix *tmp = gsl_matrix_calloc (size, size); + gsl_matrix_memcpy (tmp, idata->mm.corr); + + gsl_linalg_LU_decomp (tmp, p, &sign); + idata->detR = gsl_linalg_LU_det (tmp, sign); + gsl_permutation_free (p); + gsl_matrix_free (tmp); + } if ( factor->print & PRINT_UNIVARIATE) { + const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0; const int nc = 4; int i; - const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0; - const int heading_columns = 1; const int heading_rows = 1; @@ -1422,6 +2418,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) const int nr = heading_rows + factor->n_vars; struct tab_table *t = tab_create (nc, nr); + tab_set_format (t, RC_WEIGHT, wfmt); tab_title (t, _("Descriptive Statistics")); tab_headers (t, heading_columns, 0, heading_rows, 0); @@ -1452,44 +2449,131 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) const struct variable *v = factor->vars[i]; tab_text (t, 0, i + heading_rows, TAB_LEFT | TAT_TITLE, var_to_string (v)); - tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL); - tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL); - tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt); + tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (idata->mm.mean_matrix, i, i), NULL, RC_OTHER); + tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (idata->mm.var_matrix, i, i)), NULL, RC_OTHER); + tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->mm.n, i, i), NULL, RC_WEIGHT); } tab_submit (t); } + if (factor->print & PRINT_KMO) + { + int i; + double df = factor->n_vars * (factor->n_vars - 1) / 2; + + double w = 0; + + + double xsq; + + const int heading_columns = 2; + const int heading_rows = 0; + + const int nr = heading_rows + 4; + const int nc = heading_columns + 1; + + + + struct tab_table *t = tab_create (nc, nr); + tab_title (t, _("KMO and Bartlett's Test")); + + + tab_headers (t, heading_columns, 0, heading_rows, 0); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + tab_vline (t, TAL_2, heading_columns, 0, nr - 1); + + tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy")); + + tab_double (t, 2, 0, 0, sum_ssq_r / (sum_ssq_r + sum_ssq_a), NULL, RC_OTHER); + + tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity")); + + tab_text (t, 1, 1, TAT_TITLE, _("Approx. Chi-Square")); + tab_text (t, 1, 2, TAT_TITLE, _("df")); + tab_text (t, 1, 3, TAT_TITLE, _("Sig.")); + + + /* The literature doesn't say what to do for the value of W when + missing values are involved. The best thing I can think of + is to take the mean average. */ + w = 0; + for (i = 0; i < idata->mm.n->size1; ++i) + w += gsl_matrix_get (idata->mm.n, i, i); + w /= idata->mm.n->size1; + + xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0; + xsq *= -log (idata->detR); + + tab_double (t, 2, 1, 0, xsq, NULL, RC_OTHER); + tab_double (t, 2, 2, 0, df, NULL, RC_INTEGER); + tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL, RC_PVALUE); + + + tab_submit (t); + } + show_correlation_matrix (factor, idata); + show_covariance_matrix (factor, idata); + if (idata->cvm) + covariance_destroy (idata->cvm); -#if 1 { + gsl_matrix *am = matrix_dup (idata->analysis_matrix); gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars); - - gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace); + + gsl_eigen_symmv (am, idata->eval, idata->evec, workspace); gsl_eigen_symmv_free (workspace); + gsl_matrix_free (am); } gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC); -#endif + + idata->n_extractions = n_extracted_factors (factor, idata); + + if (idata->n_extractions == 0) + { + msg (MW, _("The %s criteria result in zero factors extracted. Therefore no analysis will be performed."), "FACTOR"); + return; + } + + if (idata->n_extractions > factor->n_vars) + { + msg (MW, + _("The %s criteria result in more factors than variables, which is not meaningful. No analysis will be performed."), + "FACTOR"); + return; + } { + gsl_matrix *rotated_factors = NULL; + gsl_matrix *pattern_matrix = NULL; + gsl_matrix *fcm = NULL; + gsl_vector *rotated_loadings = NULL; + const gsl_vector *extracted_eigenvalues = NULL; gsl_vector *initial_communalities = gsl_vector_alloc (factor->n_vars); gsl_vector *extracted_communalities = gsl_vector_alloc (factor->n_vars); size_t i; - struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, n_extracted_factors (factor, idata)); + struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions); gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors); if ( factor->extraction == EXTRACTION_PAF) { gsl_vector *diff = gsl_vector_alloc (idata->msr->size); - struct smr_workspace *ws = ws_create (analysis_matrix); + struct smr_workspace *ws = ws_create (idata->analysis_matrix); for (i = 0 ; i < factor->n_vars ; ++i) { - double r2 = squared_multiple_correlation (analysis_matrix, i, ws); + double r2 = squared_multiple_correlation (idata->analysis_matrix, i, ws); gsl_vector_set (idata->msr, i, r2); } @@ -1497,52 +2581,96 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) gsl_vector_memcpy (initial_communalities, idata->msr); - for (i = 0; i < factor->iterations; ++i) + for (i = 0; i < factor->extraction_iterations; ++i) { double min, max; gsl_vector_memcpy (diff, idata->msr); - iterate_factor_matrix (analysis_matrix, idata->msr, factor_matrix, fmw); - + iterate_factor_matrix (idata->analysis_matrix, idata->msr, factor_matrix, fmw); + gsl_vector_sub (diff, idata->msr); gsl_vector_minmax (diff, &min, &max); - + if ( fabs (min) < factor->econverge && fabs (max) < factor->econverge) break; } gsl_vector_free (diff); + + gsl_vector_memcpy (extracted_communalities, idata->msr); extracted_eigenvalues = fmw->eval; } else if (factor->extraction == EXTRACTION_PC) { - for (i = 0 ; i < factor->n_vars; ++i) - { - gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars)); - } + for (i = 0; i < factor->n_vars; ++i) + gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars)); + gsl_vector_memcpy (extracted_communalities, initial_communalities); - iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw); + iterate_factor_matrix (idata->analysis_matrix, extracted_communalities, factor_matrix, fmw); + + extracted_eigenvalues = idata->eval; } + + show_aic (factor, idata); show_communalities (factor, initial_communalities, extracted_communalities); - show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues); + if ( factor->rotation != ROT_NONE) + { + rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2); + rotated_loadings = gsl_vector_calloc (factor_matrix->size2); + if (factor->rotation == ROT_PROMAX) + { + pattern_matrix = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2); + fcm = gsl_matrix_calloc (factor_matrix->size2, factor_matrix->size2); + } + + + rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings, pattern_matrix, fcm); + } + + show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings); factor_matrix_workspace_free (fmw); show_scree (factor, idata); - show_factor_matrix (factor, idata, factor_matrix); + show_factor_matrix (factor, idata, + factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"), + factor_matrix); + + if ( factor->rotation == ROT_PROMAX) + { + show_factor_matrix (factor, idata, _("Pattern Matrix"), pattern_matrix); + gsl_matrix_free (pattern_matrix); + } + + if ( factor->rotation != ROT_NONE) + { + show_factor_matrix (factor, idata, + (factor->rotation == ROT_PROMAX) ? _("Structure Matrix") : + (factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : + _("Rotated Factor Matrix")), + rotated_factors); + + gsl_matrix_free (rotated_factors); + } + + if ( factor->rotation == ROT_PROMAX) + { + show_factor_correlation (factor, fcm); + gsl_matrix_free (fcm); + } + gsl_matrix_free (factor_matrix); + gsl_vector_free (rotated_loadings); gsl_vector_free (initial_communalities); gsl_vector_free (extracted_communalities); } +} - idata_free (idata); - casereader_destroy (r); -}