From: John Darrington Date: Thu, 18 Dec 2014 14:37:22 +0000 (+0100) Subject: FACTOR: Added PROMAX rotation. X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?p=pspp;a=commitdiff_plain;h=3a21788dd4f7a4f10d12a12e83c8890f857b5e98 FACTOR: Added PROMAX rotation. --- diff --git a/doc/statistics.texi b/doc/statistics.texi index 9cc49557ea..f17a9ece9c 100644 --- a/doc/statistics.texi +++ b/doc/statistics.texi @@ -684,7 +684,7 @@ FACTOR VARIABLES=@var{var_list} [ /EXTRACTION=@{PC, PAF@}] - [ /ROTATION=@{VARIMAX, EQUAMAX, QUARTIMAX, NOROTATE@}] + [ /ROTATION=@{VARIMAX, EQUAMAX, QUARTIMAX, PROMAX[(@var{k})], NOROTATE@}] [ /PRINT=[INITIAL] [EXTRACTION] [ROTATION] [UNIVARIATE] [CORRELATION] [COVARIANCE] [DET] [KMO] [SIG] [ALL] [DEFAULT] ] @@ -708,9 +708,13 @@ If @subcmd{PAF} is specified, then Principal Axis Factoring is used. By default Principal Components Analysis will be used. The @subcmd{/ROTATION} subcommand is used to specify the method by which the extracted solution will be rotated. -Three methods are available: @subcmd{VARIMAX} (which is the default), @subcmd{EQUAMAX}, and @subcmd{QUARTIMAX}. -If don't want any rotation to be performed, the word @subcmd{NOROTATE} will prevent the command from performing any -rotation on the data. Oblique rotations are not supported. +Three orthogonal rotation methods are available: +@subcmd{VARIMAX} (which is the default), @subcmd{EQUAMAX}, and @subcmd{QUARTIMAX}. +There is one oblique rotation method, @i{viz}: @subcmd{PROMAX}. +Optionally you may enter the power of the promax rotation @var{k}, which must be enclosed in parentheses. +The default value of @var{k} is 5. +If you don't want any rotation to be performed, the word @subcmd{NOROTATE} will prevent the command from performing any +rotation on the data. The @subcmd{/METHOD} subcommand should be used to determine whether the covariance matrix or the correlation matrix of the data is to be analysed. By default, the correlation matrix is analysed. diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c index 4094a22b43..c3a7e8ee24 100644 --- a/src/language/stats/factor.c +++ b/src/language/stats/factor.c @@ -96,6 +96,7 @@ enum rotation_type ROT_VARIMAX = 0, ROT_EQUAMAX, ROT_QUARTIMAX, + ROT_PROMAX, ROT_NONE }; @@ -131,13 +132,65 @@ quartimax_coefficients (double *x, double *y, *y = c ; } -static const rotation_coefficients rotation_coeff[3] = { +static const rotation_coefficients rotation_coeff[] = { varimax_coefficients, equamax_coefficients, - quartimax_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 { size_t n_vars; @@ -153,6 +206,7 @@ struct cmd_factor enum plot_opts plot; enum rotation_type rotation; int rotation_iterations; + int promax_power; /* Extraction Criteria */ int n_factors; @@ -262,7 +316,7 @@ ssq_od_n (const gsl_matrix *m, int n) -#if 0 +#if 1 static void dump_matrix (const gsl_matrix *m) { @@ -622,7 +676,9 @@ static void rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, const gsl_vector *communalities, gsl_matrix *result, - gsl_vector *rotated_loadings + gsl_vector *rotated_loadings, + gsl_matrix *pattern_matrix, + gsl_matrix *factor_correlation_matrix ) { int j, k; @@ -731,6 +787,125 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, 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); + } + /* reflect negative sums and populate the rotated loadings vector*/ for (i = 0 ; i < result->size2; ++i) @@ -741,7 +916,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, { double s = gsl_matrix_get (result, j, i); ssq += s * s; - sum += gsl_matrix_get (result, j, i); + sum += s; } gsl_vector_set (rotated_loadings, i, ssq); @@ -918,6 +1093,18 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) { 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); + lex_force_match (lexer, T_RPAREN); + } + factor.rotation = ROT_PROMAX; + } else if (lex_match_id (lexer, "NOROTATE")) { factor.rotation = ROT_NONE; @@ -1337,6 +1524,7 @@ static void show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm) { int i; + const int n_factors = idata->n_extractions; const int heading_columns = 1; @@ -1456,7 +1644,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) @@ -1508,18 +1698,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_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 %")); - - tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1); } for (i = 0 ; i < initial_eigenvalues->size; ++i) @@ -1580,8 +1775,11 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, { r_cum += r_percent; tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL, RC_OTHER); - 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); + 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); + } } } } @@ -1591,6 +1789,67 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, } +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, _("%d"), i + 1); + } + + for (i = 0 ; i < fcm->size2; ++i) + { + tab_text_format (t, 0, heading_rows + i, TAB_CENTER | TAT_TITLE, _("%d"), i + 1); + } + + + for (i = 0 ; i < fcm->size1; ++i) + { + for (j = 0 ; j < fcm->size2; ++j) + tab_double (t, heading_columns + i, heading_rows +j, 0, + gsl_matrix_get (fcm, i, j), NULL, RC_OTHER); + } + + tab_submit (t); +} + + static void show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata) { @@ -1938,6 +2197,8 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) { 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; @@ -2004,10 +2265,16 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) { 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); + 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); @@ -2018,16 +2285,27 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) 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->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"), + (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); @@ -2043,4 +2321,3 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) } - diff --git a/tests/automake.mk b/tests/automake.mk index 661bfbb9a5..31a7877d07 100644 --- a/tests/automake.mk +++ b/tests/automake.mk @@ -239,7 +239,8 @@ EXTRA_DIST += \ tests/data/v14.sav \ tests/language/data-io/Book1.gnm.unzipped \ tests/language/data-io/test.ods \ - tests/language/data-io/newone.ods + tests/language/data-io/newone.ods \ + tests/language/stats/llz.zsav CLEANFILES += *.save pspp.* foo* diff --git a/tests/language/stats/factor.at b/tests/language/stats/factor.at index e90e862f9a..59acba9d48 100644 --- a/tests/language/stats/factor.at +++ b/tests/language/stats/factor.at @@ -1921,4 +1921,145 @@ x,.990 y,.808 z,.723 ]) -AT_CLEANUP \ No newline at end of file +AT_CLEANUP + + +AT_SETUP([FACTOR promax]) +AT_DATA([factor-promax.sps], [dnl +set decimal=dot. +set format=F22.3. + +get file='llz.zsav'. + +factor + /variables pz pn ps nz nn ns tz tn ts oz on os sz sn ss zz zn zs + /missing listwise + /print initial extraction rotation + /criteria mineigen(1) iterate(25) + /extraction paf + /method correlation + /rotation promax (5). +]) + +AT_CHECK([ln -s $top_srcdir/tests/language/stats/llz.zsav .], [0], [ignore]) + +AT_CHECK([pspp -O format=csv factor-promax.sps], [0], [dnl +Table: Communalities +,Initial,Extraction +PZ,.191,.375 +PN,.042,.102 +PS,.458,.623 +NZ,.100,.163 +NN,.065,.079 +NS,.129,.148 +TZ,.181,.344 +TN,.102,.142 +TS,.310,.372 +OZ,.097,.158 +ON,.323,.410 +OS,.469,.617 +SZ,.104,.170 +SN,.154,.267 +SS,.081,.180 +ZZ,.123,.192 +ZN,.208,.412 +ZS,.130,.158 + +Table: Total Variance Explained +,Initial Eigenvalues,,,Extraction Sums of Squared Loadings,,,Rotation Sums of Squared Loadings +Factor,Total,% of Variance,Cumulative %,Total,% of Variance,Cumulative %,Total +1,2.968,16.491,16.491,2.411,13.393,13.393,2.355 +2,2.026,11.253,27.744,1.271,7.059,20.452,1.209 +3,1.622,9.011,36.756,.948,5.264,25.716,1.231 +4,1.086,6.032,42.788,.283,1.574,27.290,.770 +5,.996,5.533,48.321,,,, +6,.923,5.130,53.451,,,, +7,.873,4.852,58.303,,,, +8,.856,4.756,63.060,,,, +9,.836,4.644,67.703,,,, +10,.816,4.534,72.237,,,, +11,.785,4.359,76.596,,,, +12,.740,4.110,80.706,,,, +13,.713,3.964,84.670,,,, +14,.653,3.626,88.296,,,, +15,.633,3.519,91.815,,,, +16,.604,3.356,95.171,,,, +17,.484,2.687,97.858,,,, +18,.386,2.142,100.000,,,, + +Table: Factor Matrix +,Factor,,, +,1,2,3,4 +PZ,-.276,.154,.510,.124 +PN,.096,.129,-.091,.261 +PS,.746,-.085,.234,.063 +NZ,-.111,.323,.206,-.058 +NN,.007,.260,-.083,-.069 +NS,.366,.096,.046,.051 +TZ,-.228,.172,.509,.059 +TN,.131,.345,-.074,.029 +TS,.601,-.005,.098,.030 +OZ,-.145,.166,.322,-.081 +ON,.607,.082,.073,-.173 +OS,.757,-.059,.171,-.104 +SZ,-.142,.307,.226,-.066 +SN,.175,.436,-.183,.115 +SS,.199,.206,-.083,.302 +ZZ,-.074,.411,-.080,-.104 +ZN,.015,.580,-.252,-.114 +ZS,.365,.156,-.004,.015 + +Table: Pattern Matrix +,Factor,,, +,1,2,3,4 +PZ,-.063,-.126,.599,.085 +PN,-.035,.000,-.033,.325 +PS,.762,-.175,.058,.081 +NZ,.027,.230,.327,-.044 +NN,.008,.289,.008,-.026 +NS,.344,.044,.015,.091 +TZ,.004,-.074,.589,.020 +TN,.097,.307,.033,.103 +TS,.585,-.043,-.017,.062 +OZ,.046,.067,.382,-.109 +ON,.654,.151,-.029,-.145 +OS,.803,-.037,-.009,-.092 +SZ,.009,.213,.345,-.060 +SN,.065,.376,-.036,.227 +SS,.054,.042,-.013,.388 +ZZ,-.044,.434,.078,-.046 +ZN,-.025,.646,-.041,-.006 +ZS,.337,.133,-.013,.067 + +Table: Structure Matrix +,Factor,,, +,1,2,3,4 +PZ,-.177,-.058,.598,-.022 +PN,.068,.110,-.049,.317 +PS,.771,-.138,-.136,.240 +NZ,-.060,.236,.339,.019 +NN,.000,.281,.027,.076 +NS,.368,.080,-.068,.207 +TZ,-.127,-.028,.582,-.049 +TN,.122,.345,.023,.235 +TS,.607,-.018,-.160,.221 +OZ,-.074,.055,.384,-.101 +ON,.619,.104,-.160,.102 +OS,.778,-.064,-.190,.132 +SZ,-.086,.215,.361,-.009 +SN,.143,.453,-.044,.380 +SS,.171,.176,-.052,.420 +ZZ,-.073,.422,.120,.085 +ZN,-.013,.641,.008,.214 +ZS,.361,.158,-.088,.213 + +Table: Factor Correlation Matrix +Factor,1,2,3,4 +1,1.000,.008,-.232,.294 +2,.008,1.000,.065,.347 +3,-.232,.065,1.000,-.076 +4,.294,.347,-.076,1.000 +]) + + +AT_CLEANUP diff --git a/tests/language/stats/llz.zsav b/tests/language/stats/llz.zsav new file mode 100644 index 0000000000..4477cb88c1 Binary files /dev/null and b/tests/language/stats/llz.zsav differ