3
$\begingroup$

Problem

We have an array $\underline{X}$ of order $$items \times items \times people$$

(that is, $J \times J \times K$) where the cells are Pearson's correlation coefficients, each across some observations of some item pair. (See this related answer on why this correlation matrix is the only comparable data we have; we can't go back to rawer data.) These correlation coefficients can be interpreted as similarity, where high values are very similar, low values are dissimilar and 0is unrelated.

We now want to reduce the dimensionality, in all three ranks, i.e., the rank of $people$ and $items$, so that we get lower people dimensions, each of which has some lower item dimensions.

Being the Principal Components Analysis (PCA) fanboy that I am (mostly because I understand it somewhat, and, more appropriately, because we really just want dimensionality reduction), I turned to its n-mode generalisations, Candecomp/Parafac (CP), and Tucker (aka. multilinear SVD).

Which of the two (if any) is appropriate for our data?

We tried both of them (see below), but found them both a bit lacking / uncertain, so we're between that rock and hard place they talk about.


Attempt 1: Candecomp/Parafac

We're aware that our data isn't $I \times J \times K$, as would be the case for most most multi-modal arrays, but only $J \times J \times K$, so we actually have three-way, two-mode data. I had naively assumed, by analogy to PCA, that this wouldn't be a problem, since you always just pass a correlation matrix to a PCA instead of a data matrix, as in prcomp(x) == prcomp(cor(x)), except, of course that the latter can't give you scores.

Turns out, we're lucky here, and Giordani and Kiers 2014: 17ff actually describe a similar use case in their vignette for ThreeWay:

[...] The data array is such that the A- and B-mode entities (the kinship terms) coincide. It corresponds to a three-way proximity array (three-way, two-mode data), namely a collection of proximity data matrices, one for each group of students. This particular array can still be analyzed by means of CP. In this case, the CP model is usually referred to as the Indscal model to highlight that the model represents a three-way extension of the classical (two-way) multidimensional scaling. In their original paper, Carroll and Chang (1970) introduce such a model as a tool for summarizing a three-way array of scalar products through underlying components. Thus, the Indscal model can be understood as the CP model for arrays with symmetric slabs. When the CP model is fitted to proximity arrays, the symmetry of the slabs implies that two component matrices (in our case the matrices A and B) are the same. [...]

So, here we go.

Apologies for the not-quite-reprex:

  • $\underline{X}$ aka cora is too large for dput-action, but can be found here.
  • ThreeWay::CP can only be run interactively AFAIK, so I can only paste the console output.

Call:

npca_res <- CP(data = cora,
               laba = dimnames(cora)[[1]],
               labb = dimnames(cora)[[2]],
               labc = dimnames(cora)[[3]])

Messy console interaction:

WELCOME to the interactive CANDECOMP/PARAFAC analysis program
Warning: If you insert an object of mode CHARACTER when not requested, an error occurs and the program stops! 


To see ANOVA results, specify 1:
1: 
Read 0 items

How do you want to center your array?
0 = none (default)
1 = across A-mode
2 = across B-mode
3 = across C-mode
12 = across A-mode and across B-mode
13 = across A-mode and across C-mode
23 = across B-mode and across C-mode
1: 0
Read 1 item
Data have not been centered

How do you want to normalize your array?
0 = none (default)
1 = within A-mode
2 = within B-mode
3 = within C-mode
1: 0
Read 1 item
Data have not been normalized

Note: The preprocessed data are now available in Xprep.
Xprep can be used for analyses outside the Candecomp/Parafac program

To see PCA of MEAN, specify '1':
1: 0
Read 1 item

You are now about to do a Candecomp/Parafac analysis.

You have to specify the dimensionality to use.
To search this, you are advised to run a series of (unconstrained)
Candecomp/Parafac analyses using different dimensionalities
The results can next be used to find useful dimensionalities by the scree test.
To reduce the computation time, note that the maximum allowed number of iterations is 5000,
the convergence criterion is 1e-6 and the rational start via eigendecompositions is considered.
If you want to do this for choosing your dimensionality, specify '1':
1: 0
Read 1 item

You can now do a CANDECOMP/PARAFAC ANALYSIS

How many components do you want to use?
1: 3
Read 1 item

Do you want to use constraints? If so, enter '1':
1: 
Read 0 items
No constraints imposed

Specify convergence criterion (default=1e-6)
1: 
Read 0 items

By default, only a rationally started analysis run will be carried out.
To decrease the chance of missing the optimal solution, you may use additional, randomly started runs.
If you want additional runs, specify how many (e.g., 4):
1: 4
Read 1 item

Specify the maximum number of iterations you allow (default=10000).
1: 50000
Read 1 item
Run no. 1
Candecomp/Parafac function value at Start is  14265.9792386173
Candecomp/Parafac function value is 11248.0847401712 after 15 iterations
Fit percentage is 20.8684707651605 %
Procedure used 0.14 seconds
Run no. 2
Candecomp/Parafac function value at Start is  14215.4293443188
Candecomp/Parafac function value is 11284.7457771764 after 26 iterations
Fit percentage is 20.6105562856232 %
Procedure used 0.23 seconds
Run no. 3
Candecomp/Parafac function value at Start is  14219.5279522587
f= 11229.5642472746 after 50 iters; diff.= 18.7124782152951
Candecomp/Parafac function value is 11208.8164256207 after 60 iterations
Fit percentage is 21.1447277326922 %
Procedure used 0.56 seconds
Run no. 4
Candecomp/Parafac function value at Start is  14219.946679173
Candecomp/Parafac function value is 11312.4618534379 after 16 iterations
Fit percentage is 20.4155706014269 %
Procedure used 0.11 seconds
Run no. 5
Candecomp/Parafac function value at Start is  14216.8408202143
Candecomp/Parafac function value is 11298.2532050982 after 21 iterations
Fit percentage is 20.5155300254045 %
Procedure used 0.14 seconds

Fit (%) values from all runs:
Start n.1 Start n.2 Start n.3 Start n.4 Start n.5 
    20.87     20.61     21.14     20.42     20.52 

Candecomp/Parafac analysis with  3  components, gave a fit of  21.14 %
Simple check on degeneracy: inspect matrix of triple congruences
       Comp.1 Comp.2 Comp.3
Comp.1 1.0000 0.0125 0.0028
Comp.2 0.0125 1.0000 0.0004
Comp.3 0.0028 0.0004 1.0000

It is sometimes useful to SCALE solution, e.g., scale two matrices so that
they have unit sum of squares compensating this scale in remaining matrix.
If you want to scale components, specify '1':
1: 1
Read 1 item
What modes do you want to scale?
1 = B- and C-modes (scaling compensated in A) 
2 = A- and C-modes (scaling compensated in B) 
3 = A- and B-modes (scaling compensated in C) 
1: 3
Read 1 item

Solution for A, B and C in detail.
A
                     Comp.1 Comp.2 Comp.3
but-how               -0.17  -0.12  -0.15
encyclopedia           0.17  -0.05   0.22
alien                  0.17  -0.25  -0.21
language-of-bees      -0.17   0.12   0.19
bad-hen                0.18  -0.26  -0.20
correspondence         0.16   0.13  -0.17
bamboozable           -0.17  -0.17   0.14
inventions             0.16   0.27   0.19
gray-hair              0.17  -0.06   0.20
i-we                  -0.17   0.15  -0.12
the-same              -0.17   0.20  -0.14
the-better            -0.17   0.17   0.23
cockatoo              -0.17   0.05   0.23
caravan                0.16   0.05  -0.08
comma                 -0.17   0.24  -0.13
headturner             0.17  -0.11   0.21
crocodile             -0.16  -0.29   0.04
countries              0.17  -0.06   0.27
level                  0.16  -0.26  -0.11
morning               -0.17   0.29  -0.14
idiom                  0.17  -0.29  -0.18
lullaby                0.17  -0.04   0.18
riddle                 0.17  -0.07   0.23
eating-grandpa        -0.17   0.16  -0.14
easter-bunny          -0.17  -0.14   0.15
alphabet-of-swearing   0.17   0.17  -0.09
guilt                  0.17  -0.10  -0.07
trouble               -0.18   0.08  -0.11
eating-animals         0.17   0.05  -0.19
random-poetry         -0.16   0.12  -0.08
comparatively         -0.16  -0.15  -0.09
one-letter             0.17  -0.07   0.21
resistance            -0.18   0.21  -0.17
conjugation            0.17   0.00   0.22
censored              -0.17  -0.22  -0.12
B
                     Comp.1 Comp.2 Comp.3
but-how                0.17  -0.12  -0.15
encyclopedia          -0.17  -0.05   0.22
alien                 -0.17  -0.25  -0.21
language-of-bees       0.17   0.12   0.19
bad-hen               -0.18  -0.26  -0.20
correspondence        -0.16   0.13  -0.17
bamboozable            0.17  -0.17   0.14
inventions            -0.16   0.27   0.19
gray-hair             -0.17  -0.06   0.20
i-we                   0.17   0.15  -0.12
the-same               0.17   0.20  -0.14
the-better             0.17   0.17   0.23
cockatoo               0.17   0.05   0.23
caravan               -0.16   0.05  -0.08
comma                  0.17   0.24  -0.13
headturner            -0.17  -0.11   0.21
crocodile              0.16  -0.29   0.04
countries             -0.17  -0.06   0.27
level                 -0.16  -0.26  -0.11
morning                0.17   0.29  -0.14
idiom                 -0.17  -0.29  -0.18
lullaby               -0.17  -0.04   0.18
riddle                -0.17  -0.07   0.23
eating-grandpa         0.17   0.16  -0.14
easter-bunny           0.17  -0.14   0.15
alphabet-of-swearing  -0.17   0.17  -0.09
guilt                 -0.17  -0.10  -0.07
trouble                0.18   0.08  -0.11
eating-animals        -0.17   0.05  -0.19
random-poetry          0.16   0.12  -0.08
comparatively          0.16  -0.15  -0.09
one-letter            -0.17  -0.07   0.21
resistance             0.18   0.21  -0.17
conjugation           -0.17   0.00   0.22
censored               0.17  -0.22  -0.12
C
          Comp.1 Comp.2 Comp.3
Willy      -0.26   2.39   1.16
Karen      -0.31   3.17   1.09
Kristina   -0.41   9.25   0.53
Stefan     -0.92   1.29   1.85
Niklas      0.02   1.70   3.01
Fiona      -1.03   2.34   1.98
Finley     -1.14   3.67   1.82
Zoe        -0.54   3.12   2.04
Annika     -1.27   0.53   2.27
Ian        -1.18   2.06   1.12
Frank      -0.54   2.29   0.99
William    -1.23   2.12   0.76
Jana       -0.92   3.42   2.24
Babsi      -0.07   4.01   2.14
Collin     -0.80   1.27   2.54
Patrizia   -3.57   3.90   1.18
Valentine  -0.33   1.19   0.96
Jennifer  -34.99  -0.19  -0.01
Julius     -2.46   2.65   1.75
Thomas     -0.66   2.33   1.22
Kim        -1.20   3.80   1.59
Tamara     -1.16   3.89   2.20
Caroline   -1.18   4.82   5.20
Dennis     -1.43   3.71   1.78
Shae       -0.43   1.50   1.04
Knut       -0.27   2.31   6.99
Angelina   -0.71   2.78   2.93
Victor     -0.96   4.09   0.92
Fred       -1.80   3.80   1.33
Nora       -0.58   5.53   2.25
Carmen     -1.06   2.89   3.04
Tashina    -0.97   1.13   1.82
Tessa      -0.69   9.30   0.65
Adrian     -0.62   3.56   1.68
Constanze   0.66   8.07   2.88
Justin     -0.32   5.63   1.89
Wanda      -0.56   3.69   1.21
Irene      -0.85   8.40   2.07
Anastasia  -0.06   2.33  22.67
Sabrina    -0.57   3.49   2.03
Jasmine    -1.15   6.08   6.53
Dana       -2.61   3.84   2.33
Amanda      0.04   0.84   0.18
Jalena     -2.01   5.25   1.90
Marcel     -0.28   4.14   1.30
Leandro    -0.72   4.02   3.28
Justus     -0.58   2.54   3.09
Hagen      -1.57   6.63   1.47
Ralf       -1.30   4.06   3.31
Shaana     -1.07   0.48   2.18
Nhome      -1.17   2.09   3.81
Azra       -2.51   3.66   1.25

You can now manually PERMUTE and REFLECT columns of solution
If you want to reflect/permute columns, specify '1':
1: 
Read 0 items

If you want to carry out a STABILITY CHECK on current or different solution, specify '1':
1: 
Read 0 items

If you want to carry out a BOOTSTRAP procedure for
computing confidence intervals for the current solution, specify '1':
1: 
Read 0 items

If you want to carry out a FITPARTITIONING on current solution, specify '1':
1: 
Read 0 items

Press 'return' to conclude the analysis
1: 
Read 0 items
The component matrices are normalized such that A and B have unit sums of squares
To see component matrices and fit and other results, digits: $Xprep, $A, $B, $C, $fit, $fitA, $fitB, $fitC, ...

The fit isn't great, but that might just be our crappy / preliminary /small dataset:

npca_res$fit
#> Fit (%) 
#21.14473 

Happily, A-mode matrices ($Items$) and B-mode matrices ($Items$), are same-ish (computational artefact, I'm guessing, plus random reflection?)

all(abs(round(x = npca_res$A, digits = 1)) == abs(round(x = npca_res$B, digits = 1)))
#>[1] TRUE

Worries and problems with Candecomp/Parafac

  1. Giordani and Kiers (ibid.) recommend this to an array of proximity matrices, not correlation matrices. Can we even do this? Is this ok?
  2. CP, by definition, as I understand, always extracts the same number of components on each mode, but we have no reasonable expectation that that would/should be so, empirically. In fact, it seems to us there should be more $item$ components than $people$ components.
  3. Any thoughts on what would be appropriate constraint, scaling, centering and transformation options? (we're thinking no on the latter three, b/c correlation coefficients already are centered).

Attempt 2: Tucker

So, to allow us to extract different (numbers of) components for each mode, we turn to Tucker, in flagrant disregard of actually having only an $J \times J \times K$ array, not $I \times J \times K$. Weirdly, it works.

Same apologies for non-reprex apply.

Call:

tucker_res <- T3(data = cora,
                 laba = dimnames(cora)[[1]],
                 labb = dimnames(cora)[[2]],
                 labc = dimnames(cora)[[3]])

Messy Terminal Interaction is here because too long for SO.

As advised by Giordani and Kiers (2014: 8) we then do some post-processing to make the core matrix an identity matrix (I think?):

tucker_res$A <- tucker_res$A %*% tucker_res$core
tucker_res$core <- solve(tucker_res$core) %*% tucker_res$core

This ominously fails with:

#> Error in solve.default(tucker_res$core) : 'a' (3 x 6) must be square

$Items$ and $Items$ mode matrices are the same:

all(round(tucker_res$A, digits = 1) == round(tucker_res$B, digits = 1))
#> [1] TRUE

but the core ain't an identity matrix, making this weird/hard to interpret:

tucker_res$core
#>         B1xC1       B2xC1       B3xC1      B1xC2      B2xC2      B3xC2
#> A1 -0.2351972 -1.59675721  6.47329544 36.2976853 -0.6195908  2.0033419
#> A2 -1.5995439 25.46542681 -0.05905521 -0.6205366 -0.4781704 -1.1461151
#> A3  6.4725537 -0.03896677 24.04316413  2.0032077 -1.1465854 -0.3047064

Worries and problems with Tucker

  1. Plugging in an $items \times items \times people$ matrix, can we just do this?!? Are the results interpretable or is this just garbage-in-garbage-out?
  2. How do we interpret a non-identity core matrix?
  3. Is this, overall, appropriate?
$\endgroup$
2
  • $\begingroup$ incidentally, if someone could sponsor a nPCA, multiway, Tucker or Candecomp/Parafac tag, that would be great – couldn't find anything remotely applicable from the available tags. $\endgroup$
    – maxheld
    Aug 24 '16 at 17:30
  • 1
    $\begingroup$ One note (relevant or not- you to decide). JxJ dim. red. task is MDS task. The JxK dim. red. task is MDU task. Also, (JxJ)xK could be conceptualized as INDSCAL or similar version of MDS, which meas that one can first do the MDS obtaining the K maps and then may find some approach to summarise the K to a lower "dimension". Just speculating. $\endgroup$
    – ttnphns
    Aug 29 '16 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.