Shrinkage plot of PCA component loadings

[dmbates]

@palday and I were discussing in email the idea of creating a shrinkage plot on the scale of the PCA component loadings. That would show that the last one or two components often have very little or no variability after shrinkage.

One issue that this brings up is whether to use a correlation matrix or a covariance matrix to generate the loadings. I have a vague sense that the scale should be determined by the diagonal elements in each block of the leading blocks of A.

If the experimental factors are a two-level factorial with -1/+1 encoding then the diagonal elements in each block for each subject or item are constant if the data are balanced.

julia> using MixedModels

julia> mrk17 = MixedModels.dataset(:mrk17_exp1);

julia> contr = merge(Dict(nm => EffectsCoding() for nm in (:F, :P, :Q, :lQ, :lT)), Dict(nm => Grouping() for nm in (:subj, :item)));

julia> m6form = @formula 1000 / rt ~ 1 + F*P*Q*lQ*lT + (1+F+P+Q+lQ+lT|subj) + (1+P+Q+lQ+lT|item);

julia> m6 = fit(MixedModel, m6form, mrk17; contrasts=contr);

julia> view(m6.A[1].data, :, :, 1)
5×5 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
  54.0  -10.0    2.0    8.0    4.0
 -10.0   54.0   -2.0   -4.0    4.0
   2.0   -2.0   54.0  -16.0  -12.0
   8.0   -4.0  -16.0   54.0   -6.0
   4.0    4.0  -12.0   -6.0   54.0

This means that the columns of Z are the same length within blocks determined by the levels of the grouping factor. I would take this as an indication that the coefficients are "a priori" on the same scale so that one could use the SVD of λ to determine the coefficient loadings.

julia> svd(m6.λ[1])
SVD{Float64, Float64, Matrix{Float64}}
U factor:
5×5 Matrix{Float64}:
 -0.995229   -0.0491655  -0.0629712  -0.0544327   0.0131547
  0.0138492  -0.673998   -0.0936185   0.302618   -0.667228
  0.0840198  -0.220143   -0.912225   -0.221178    0.2518
  0.0393347  -0.0581186   0.136058   -0.914559   -0.37436
  0.0268372  -0.701049    0.369599   -0.141851    0.592525
singular values:
5-element Vector{Float64}:
 0.19417949889806668
 0.05647856144891754
 0.04143641584593457
 0.01915808430179067
 0.0
Vt factor:
5×5 Matrix{Float64}:
 -0.999791     0.0146879   0.013155    0.00537781  -0.0
 -0.0143657   -0.983511    0.0735867  -0.164576     0.0
 -0.0134992   -0.100885   -0.98095     0.165459     0.0
 -0.00539504   0.149375   -0.1793     -0.972373    -0.0
  0.0          0.0         0.0         0.0          1.0

julia> svd(m6.λ[2])
SVD{Float64, Float64, Matrix{Float64}}
U factor:
6×6 Matrix{Float64}:
 -0.996135    -0.0768322   -0.0377731   0.0174735   0.00647847   0.00616525
  0.0137053    0.00193099  -0.146545    0.402011   -0.515128     0.742538
  0.0228672   -0.0556326   -0.345166    0.639659   -0.293233    -0.618138
  0.0675055   -0.489241    -0.771432   -0.151566    0.334395     0.16182
  0.00410365  -0.218233    -0.102201   -0.604129   -0.732544    -0.200796
 -0.0493025    0.839056    -0.5024     -0.202419   -0.0132126   -1.3417e-17
singular values:
6-element Vector{Float64}:
 0.5999594688347276
 0.11604256361820424
 0.08470275955425068
 0.02923672643093188
 0.024456777172460642
 4.392629796786411e-19
Vt factor:
6×6 Matrix{Float64}:
 -0.999874      0.00337994   0.00944645    0.0100762     0.000837818  -0.00700455
 -0.0149165     0.00481949  -0.394519     -0.677365     -0.0737185     0.616321
 -0.00535286   -0.585802    -0.621029     -0.108937      0.0604721    -0.505575
  0.000854701   0.553073     0.0673235    -0.579936      0.0705873    -0.590142
  0.00026508   -0.592377     0.673837     -0.439183      0.00550799   -0.0460494
  0.0           7.3454e-16  -6.78585e-16   3.43193e-16   0.992923      0.118764

Turns out, of course, that the singular directions are not very different from those indicated in

julia> m6.PCA
(item = 
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .
 P: unr       -0.05   1.0     .      .      .
 Q: deg       -0.36   0.38   1.0     .      .
 lQ: deg      -0.37   0.03   0.03   1.0     .
 lT: WD       -0.11   0.87   0.01   0.35   1.0

Normalized cumulative variances:
[0.4223, 0.6864, 0.9056, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5
 (Intercept)  -0.3    0.67  -0.01   0.67   0.07
 P: unr        0.6    0.38   0.22  -0.05  -0.67
 Q: deg        0.31  -0.34   0.7    0.46   0.28
 lQ: deg       0.31  -0.41  -0.62   0.55  -0.2
 lT: WD        0.6    0.34  -0.26  -0.15   0.66, subj = 
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .
 F: LF        -0.36   1.0     .      .      .      .
 P: unr       -0.35   0.89   1.0     .      .      .
 Q: deg       -0.41   0.45   0.73   1.0     .      .
 lQ: deg      -0.06   0.17   0.18   0.58   1.0     .
 lT: WD        0.26   0.1    0.02  -0.37  -0.51   1.0

Normalized cumulative variances:
[0.4768, 0.7357, 0.8811, 0.9429, 1.0, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6
 (Intercept)  -0.34  -0.02  0.84  -0.3   -0.28  -0.11
 F: LF         0.45   0.41  0.12  -0.42   0.43  -0.5
 P: unr        0.51   0.35  0.16  -0.16  -0.28   0.7
 Q: deg        0.52  -0.15  0.11   0.42  -0.55  -0.46
 lQ: deg       0.32  -0.51  0.43   0.28   0.57   0.22
 lT: WD       -0.21   0.65  0.24   0.67   0.18   0.0)

However, the component loadings from an SVD of λ can be applied to the results of ranef directly. I don't think the same is true for the loadings from the current PCA. We would need to undo the scaling in the covariance -> correlation conversion.