procrustes {lazy.procrustes}

R Documentation

Procrustes Rotation with Prescribed Factor Correlations

Description

Procrustes Rotation with Prescribed Factor Correlations

Usage

procrustes(
  A,
  C,
  Phia = diag(ncol(A)),
  Phib = NULL,
  Tinit = NULL,
  maxiter = 500,
  eps = 1e-06,
  epsd = 1e-06,
  maxiter2 = 200,
  eps2 = 1e-06,
  epsd2 = 0.001,
  SQUAREM = 3,
  nSQUAREM = 1,
  minalpha = -999,
  maxalpha = -1,
  always = 0,
  reset1 = 0,
  reset2 = 1,
  print = 1
)

Arguments

A	The matrix to be rotated
C	The target matrix
Phia	The factor correlation matrix associated with A
Phib	The factor correlation of the rotated factors
Tinit	initial value of T matrix
maxiter	max # of iterations
eps	convergence criterion for rmse
epsd	convergence criterion for maximum absolute differences of Q.
maxiter2	max # of iterations of procotWKS
eps2	convergence criterion for rmse of procotWKS
epsd2	convergence criterion for maximum absolute differences of procotWKS
SQUAREM	= 3 : See the help of iSQUAREM in lazy.accel package.
nSQUAREM	= 1 : See the help of iSQUAREM in lazy.accel package.
minalpha	= -999 : See the help of iSQUAREM in lazy.accel package.
maxalpha	= -1 : See the help of iSQUAREM in lazy.accel package.
always	= 1 : See the help of iSQUAREM in lazy.accel package.
reset1	= 0 : See the help of iSQUAREM in lazy.accel package.
reset2	= 1 : See the help of iSQUAREM in lazy.accel package.
print	= 1 to print the result

Details

This function finds the factor rotation matrix T of the form:
g=T'f and B=A inv(T')
which minimizes the least squares criterion:
RSS = tr( (C - B)'(C - B) )
subject to corr(g)=Phib.
The rotation matrix T can be defined as:
where T = P inv(K) Q D R' and QQ'=Q'Q=I,
where corr(f)=Phia=P K2 P' and corr(g)=Phib=R D2 R'.
The missing elements of C matrix will be estimated so that they also minimize RSS.

Value

A list of B, T, Q, Cm, A, C, Phia, Phib, rmse,
where B is the rotated matrix, T is the rotation matrix,
Q is the orthogonal matrix which defines T, and
Cm is the target matrix with its missing elements replaced by LSE.

Examples

# Independent Cluster
seed <- 1701
set.seed(seed)

nvar <- 20
ndim0 <- 3
ps <- 0.1
df <- 500
phi <- 0.3
big=0.8

Lambda0 <- gendatafa_A( nvar, ndim0, large=big,small=1-big
                  , pc=0, sd=0 )$loadings
colnames(Lambda0) <- paste("f",1:ndim0,sep="")

Phi <- (1-phi)*diag(ndim0)+phi*matrix(1,ndim0,ndim0)
Sigma <- Lambda0%*%Phi%*%t(Lambda0)+ps*diag(nvar)
dS <- sqrt(diag(Sigma))
Sigma <- diag(1/dS)%*%Sigma%*%diag(1/dS)

S <- rWishart( 1, df, Sigma )
S <- S[,,1]/df
dS <- sqrt(diag(S))
S <- diag(1/dS)%*%S%*%diag(1/dS)

ndim <- ndim0


# temp <- lazy.fa::fa_hs( S, ndim=ndim, c="smc" )
# Lambda <- temp$Lambda
# psi <- temp$psic

temp <- eigen(S)
Lambda <- temp$vectors[,1:ndim]%*%diag(sqrt(temp$values[1:ndim]))
psi=diag(rep(mean(diag(S-Lambda%*%t(Lambda))),nvar))

Lambda00 <- Lambda0
Lambda00[Lambda00==big] <- 1
Lambda00[Lambda00==1-big] <- 0
Lambda000=Lambda00
Lambda000[Lambda00==1]=NA

phia <- 0; Phia <- (1-phia)*diag(ndim)+phia*matrix(1,ndim,ndim)
table <- NULL
for( p in seq(-0.45, 0.45, 0.05) ){
 Print(p)
 Phib <- (1-p)*diag(ndim)+p*matrix(1,ndim,ndim)
 res <- procrustes( Lambda, Lambda0, Phia=Phia, Phib, print=1, maxiter2=500 )
 table <- rbind(table,c(p,res$rmse))
 Print(res$B)
}
best=table[,2]==min(table[,2])
Print(table,best)

---