Log of Varimax Steps in R for HW7
=================================

> Lam1        
       [,1]       [,2]
[1,] 1.6480 -0.4018537
[2,] 2.6036 -1.0866706
[3,] 0.8122  0.9652192
[4,] 0.7853  1.1642072
[5,] 0.7575  1.0965249
[6,] 0.8911  1.0802817
    ### example of a two-column matrix of factor loadings

> varimax(Lam1)
$loadings
          [,1]      [,2]
[1,] 1.6427425 0.4228328
[2,] 2.8083618 0.2696104
[3,] 0.2610141 1.2341753
[4,] 0.1434510 1.3969597
[5,] 0.1508554 1.3241661
[6,] 0.2763260 1.3728480

$rotmat
           [,1]      [,2]
[1,]  0.8818138 0.4715976
[2,] -0.4715976 0.8818138

### Instead, make Lam0 (by normalizing the approximately 
orthogonal columns to have length 1) and apply varimax to it:

> Lam0 = Lam1 %*% diag(1/sqrt(diag(t(Lam1)%*%Lam1)))
          [,1]       [,2]
[1,] 0.4730076 -0.1640733
[2,] 0.7472831 -0.4436780
[3,] 0.2331170  0.3940904
[4,] 0.2253962  0.4753355
[5,] 0.2174170  0.4477014
[6,] 0.2557628  0.4410694

> varimax(Lam0)$load
             [,1]        [,2]
[1,]  0.496069222  0.06761338
[2,]  0.867152249 -0.05769909
[3,]  0.029658830  0.45691481
[4,] -0.013974705  0.52588206
[5,] -0.008591932  0.49762724
[6,]  0.028606899  0.50905647


### These steps indicate that we should first 
### create the orthonormal column matrices Lambda0,
### and then apply varimax to those and look for 
### components close to 1, 0.
### NB the result is still orthonormal.
###   But it seems logical to allow re-scaling of each 
###   column and then look for components which 
###   are <= 0.2 or >= 0.8 times the maximum component
###   for the column. 
### In this example, only 1 of 12 components is 
###   not <=0.2 or >= 0.8 times the max in column.

========================================================

## We next simulate a  p=5  and  q=2  factor-model 
##   dataset, with Lam = cbind(c(1,2,1,2,1),
##      c(-1,-1, .75, .75, .75))  and  psi0 = c(2,1,3,1,2)

## But note that I am simulating this time from a model in 
##   which there is NOT simple structure, even after rotation!

> Lam = cbind(c(1,2,1,2,1), c(-1,-1, .75, .75, .75))
  varimax(Lam)
...
Loadings:
     [,1]   [,2]
[1,]        -1.413
[2,]  0.790 -2.092
[3,]  1.243 -0.128
[4,]  1.978 -0.807
[5,]  1.243 -0.128

Compare this with:
> varimax(cbind(rep(1,6),rep(c(1,-1),3)))
...
Loadings:
     [,1] [,2]
[1,]  1    1
[2,]  1   -1
[3,]  1    1
[4,]  1   -1
[5,]  1    1
[6,]  1   -1

--------------------------------------------------
We continue with simulated example:

> SigMat = Lam %*% t(Lam) + diag(c(2,1,3,1,2))
  tmp = eigen(SigMat)
  SigRt = tmp$vec %*% diag(sqrt(tmp$val)) %*% t(tmp$vec)
  Yarr = array(rnorm(600*5), c(600,5)) %*% SigRt

> tmpfac = factanal(Yarr,2,rotation="none")
Uniquenesses:
[1] 0.304 0.284 0.644 0.127 0.555

Loadings:
     Factor1 Factor2
[1,]  0.498   0.669
[2,]  0.745   0.402
[3,]  0.520  -0.292
[4,]  0.909  -0.218
[5,]  0.567  -0.350

               Factor1 Factor2
SS loadings      2.221   0.864
Proportion Var   0.444   0.173
Cumulative Var   0.444   0.617

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 0.5 on 1 degree of freedom.
The p-value is 0.48

### Compare with results from direct LikProf maximization

> SY = var(Yarr)
  LikAlt = function(logps) LikProf(exp(logps),SY,2)[[1]]
  tmp = nlm(LikAlt, runif(5))
> tmp
$minimum
[1] 11.06982
$estimate
[1]  0.2026265  0.5573210  1.0617163 -0.2520151  0.5516301
$gradient
[1]  8.348877e-08  1.415756e-06 -1.273229e-06 -1.621814e-06  2.129852e-06
$code
[1] 1
$iterations
[1] 14

### So this does give convergence
> psiLim = exp(tmp$est)
  tmp = LikProf(psiLim,SY,2,T)
  psiLim = psiLim*attr(tmp,"sSqrHat")
> psiLim
[1] 1.2741407 1.8165999 3.0082601 0.8086658 1.8062912
### compares decently well with c(2,1,3,1,2), but not great !

> cbind(attr(tmp,"LamStHat"), Lam)    ### again:decent correspondence 
          [,1]       [,2] [,3]  [,4]   ### but not great
[1,] -0.9223185  1.2382688    1 -1.00
[2,] -1.4266413  0.7695351    2 -1.00
[3,] -0.6608765 -0.3705369    1  0.75
[4,] -2.6027306 -0.6234629    2  0.75
[5,] -0.7763759 -0.4795919    1  0.75

> auxLam = attr(tmp,"LamStHat")
  t(auxLam) %*% auxLam
### essentially diagonal
> auxfac = diag(t(auxLam) %*% auxLam)> auxfac
[1] 10.699701  2.881506
> auxLam %*% diag(-1/sqrt(auxfac))
           [,1]       [,2]
[1,] 0.2819649 -0.7294662    ### these are the fitted unit factor
[2,] 0.4361431 -0.4533344    ###  directions: very different from
[3,] 0.2020387  0.2182839    ###  the unrotated factors given by
[4,] 0.7956891  0.3672830    ###  the function "factanal" !?
[5,] 0.2373483  0.2825284

> Vfac = varimax(auxLam)
> Vfac
...
Loadings:
     [,1]   [,2]
[1,]         1.544
[2,] -0.692  1.466
[3,] -0.752
[4,] -2.465  1.043
[5,] -0.909

> round(c(.692,.752,.909)/2.465,3)
[1] 0.281 0.305 0.369
> round(c(1.466,1.043)/1.544,3)
[1] 0.949 0.676

### Here 3 of the loadings in the first column and 1 in 
###   the second column are in the "between" category 
###   (between .2 and .8 times the largest in column) .