ROT_VARIMAX = 0,
ROT_EQUAMAX,
ROT_QUARTIMAX,
+ ROT_PROMAX,
ROT_NONE
};
*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;
enum plot_opts plot;
enum rotation_type rotation;
int rotation_iterations;
+ int promax_power;
/* Extraction Criteria */
int n_factors;
-#if 0
+#if 1
static void
dump_matrix (const gsl_matrix *m)
{
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;
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)
{
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);
{
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;
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;
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)
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)
{
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);
+ }
}
}
}
}
+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)
{
{
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;
{
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);
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);
}
-
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
projects / pspp / commitdiff