Nothing
varimcoco <-
function(A,B,C,H,wa_rel,wb_rel,wc_rel,rot1,rot2,rot3,nanal){
conv=1e-6
narg=nargs()
if (narg<11){
nanal=5
}
if (narg<8){
rot1=1
rot2=1
rot3=1
nanal=5
}
n=nrow(A)
r1=ncol(A)
m=nrow(B)
r2=ncol(B)
p=nrow(C)
r3=ncol(C)
cat(paste("Results from ",nanal," random starts"),fill=TRUE)
# CHECK INPUT SPECIFICATION
sz=dim(H)
if ((sz[1]!=r1) | (sz[2]!=r2*r3)){
stop("Error: sizes of A,B and C incompatible with H!")
}
gam1=1
gam2=1
gam3=1
gama=1
gamb=1
gamc=1
# Compute Natural Weights for core rotation:
# ensure natural proportion for three core terms
w1=1/r2/r3
w2=1/r1/r3
w3=1/r1/r2
ww=(w1+w2+w3)
w1=w1/ww # normalizes weights to sum of 1
w2=w2/ww
w3=w3/ww
ssh=SUM(H)$ssq^2/r1/r2/r3 # mean sq of squared core elements
# makes contribution of qmax on uniform H equal to 1
w1=w1/ssh
w2=w2/ssh
w3=w3/ssh
v1=gam1*w1/(w1+w2+w3)
v2=gam2*w2/(w1+w2+w3)
v3=gam3*w3/(w1+w2+w3)
# Compute WEIGHTS for columnwise orthonormal component matrices
# relative weights times natural weights (these are "natural" given columnwise orthonormality!)
wa=n/r1
wb=m/r2
wc=p/r3
wa=wa_rel*wa
wb=wb_rel*wb
wc=wc_rel*wc
# START LOOP FOR .. RUNS OF ALGORITHM
Ssup=matrix(,r1*r1,nanal)
Tsup=matrix(,r2*r2,nanal)
Usup=matrix(,r3*r3,nanal)
func=c(1:nanal)
for (ana in 1:nanal){
if (rot1==1){
set.seed(ana^2,kind=NULL,normal.kind=NULL)
S=orth(matrix(runif(r1*r1,0,1),r1,r1)-.5)
} else{
S=diag(1,nrow=r1)
}
if (rot2==1){
set.seed(ana^2,kind=NULL,normal.kind=NULL)
T=orth(matrix(runif(r2*r2,0,1),r2,r2)-.5)
} else{
T=diag(1,nrow=r2)
}
if (rot3==1){
set.seed(ana^2,kind=NULL,normal.kind=NULL)
U=orth(matrix(runif(r3*r3,0,1),r3,r3)-.5)
} else{
U=diag(1,nrow=r3)
}
# Initialize Y, AS, BT and CU
# multiply by weights^.25
Y=(w1+w2+w3)^(.25)*S%*%H%*%kronecker(t(U),t(T))
AS=wa^(.25)*A%*%t(S)
BT=wb^(.25)*B%*%t(T)
CU=wc^(.25)*C%*%t(U)
# EVALUATE FUNCTION VALUE AT START
# three-way orthomax part
f1=sum((t(Y)*t(Y))^2)-v1/(r2*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r1,r2,r3) # Y of order r2 x r3r
f1=f1-v2/(r1*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r2,r3,r1) # Y of order r3 x r1r2
f1=f1-v3/(r1*r2)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r3,r1,r2) # Y of order r1 x r2r3
# orthomax for component matrices
L=AS
L=L*L
f2a=sum(L^2)-gama/n*sum((colSums(L))^2)
L=BT
L=L*L
f2b=sum(L^2)-gamb/m*sum((colSums(L))^2)
L=CU
L=L*L
f2c=sum(L^2)-gamc/p*sum((colSums(L))^2)
f2=f2a+f2b+f2c
f=f1+f2
# START OF ITERATIVE PART
iter=0
fold=f-2*conv*abs(f)
while ((abs(f-fold))>(conv*abs(f))){
iter=iter+1
fold=f
# UPDATING S
if (rot1==1){
# orthomax of Y' = (w1+w2+w3)^.25 H_F ' with gamma = v1
# + orthomax of AS = wa^.25A with gamma = gama
ORTHMAXs=orthmax2(t(Y),AS,v1,gama,conv)
AS=ORTHMAXs$B2
K=ORTHMAXs$B1
S=t(ORTHMAXs$T)%*%S
}
if (rot1==0){
K=t(Y)
}
Y=permnew(t(K[1:(r2*r3),]),r1,r2,r3) # Y of order r2 x r3r1
# UPDATING T
if (rot2==1){
# orthom of Y' = (w1+w2+w3)^.25 H_F ' with gamma=v2
#+ orthom of BT = wb^.25B with gamma=gamb
ORTHMAXt=orthmax2(t(Y),BT,v2,gamb,conv)
BT=ORTHMAXt$B2
K=ORTHMAXt$B1
T=t(ORTHMAXt$T)%*%T
}
if (rot2==0){
K=t(Y)
}
Y=permnew(t(K[1:(r1*r3),]),r2,r3,r1) # Y of order r3 x r1r2
# UPDATING U
if (rot3==1){
# orthom of Y' = (w1+w2+w3)^.25 H_F ' with gamma=v3
# + orthom of CU = wc^.25C with gamma=gamc
ORTHMAXu=orthmax2(t(Y),CU,v3,gamc,conv)
CU=ORTHMAXu$B2
K=ORTHMAXu$B1
U=t(ORTHMAXu$T)%*%U
}
if (rot3==0){
K=t(Y)
}
Y=permnew(t(K[1:(r1*r2),]),r3,r1,r2) # Y of order r1 x r2r3
# EVALUATE ORTHCOCO FUNCTION:
# three-way orthomax part
f1=sum((t(Y)*t(Y))^2)-v1/(r2*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r1,r2,r3) # Y of order r2 x r3r1
f1=f1-v2/(r1*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r2,r3,r1) # Y of order r3 x r1r2
f1=f1-v3/(r1*r2)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r3,r1,r2) # Y of order r1 x r2r3
# orthomax for component matrices
L=AS
L=L*L
f2a=sum(L^2)-gama/n*sum((colSums(L))^2)
L=BT
L=L*L
f2b=sum(L^2)-gamb/m*sum((colSums(L))^2)
L=CU
L=L*L
f2c=sum(L^2)-gamc/p*sum((colSums(L))^2)
f2=f2a+f2b+f2c
f=f1+f2
}
cat(paste("Run no. ",ana," f = ",round(f,digits=3),"(core: ",round(f1,digits=3),"; A: ",round(f2a,digits=3),"; B: ",round(f2b,digits=3),"; C: ",round(f2c,digits=3),"), ",iter," iters"),fill=TRUE)
func[ana]=f
Ssup[,ana]=S
Tsup[,ana]=T
Usup[,ana]=U
}
# FINISH LOOP FOR ANALYSES
# REPORT BEST SOLUTION
f=max(func)
maxi=which.max(func)
# Rotation matrices
S=matrix(Ssup[,maxi],nrow(S),ncol(S))
T=matrix(Tsup[,maxi],nrow(T),ncol(T))
U=matrix(Usup[,maxi],nrow(U),ncol(U))
cat(paste("Three-way orthomax function value for best solution is ",round(f,digits=4)),fill=TRUE)
# Rotated core:
AS=A%*%t(S)
BT=B%*%t(T)
CU=C%*%t(U)
# reflect to positive sums in A,B and C
sg=sign(colSums(AS))
for (i in 1:length(sg)){
if(sg[i]==0){
sg[i]=1
}
}
S=diag(sg,nrow=r1)%*%S
AS=AS%*%diag(sg,nrow=r1)
sg=sign(colSums(BT))
for (i in 1:length(sg)){
if (sg[i]==0){
sg[i]=1
}
}
T=diag(sg,nrow=r2)%*%T
BT=BT%*%diag(sg,nrow=r2)
sg=sign(colSums(CU))
for (i in 1:length(sg)){
if (sg[i]==0){
sg[i]=1
}
}
U=diag(sg,nrow=r3)%*%U
CU=CU%*%diag(sg,nrow=r3)
# compute core
K=S%*%H%*%kronecker(t(U),t(T))
# recompute final orthomax values
L=AS
L=L*L
f2a=(sum(L^2)-gama/n*sum((colSums(L))^2))*n/r1
L=BT
L=L*L
f2b=(sum(L^2)-gamb/m*sum((colSums(L))^2))*m/r2
L=CU
L=L*L
f2c=(sum(L^2)-gamc/p*sum((colSums(L))^2))*p/r3
# three-way orthomax part
Y=((w1+w2+w3)^.25)*K
f1=sum((t(Y)*t(Y))^2)-v1/(r2*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r1,r2,r3) # Y of order r2 x r3r1
f1=f1-v2/(r1*r3)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r2,r3,r1) # Y of order r3 x r1r2
f1=f1-v3/(r1*r2)*sum((colSums(t(Y)*t(Y)))^2)
Y=permnew(Y,r3,r1,r2) # Y of order r1 x r2r3
cat("Varimax values of core and AS, BT and CU (all based on natural weights)",fill=TRUE)
cat(" Core A B C ",fill=TRUE)
cat(paste(" ",round(f1,digits=3)," ",round(f2a,digits=3)," ",round(f2b,digits=3)," ",round(f2c,digits=3)),fill=TRUE)
cat("These simplicity values are based on 'natural' weights and therefore comparable across matrices",fill=TRUE)
cat("When multiplied by the relative weights, they give the contribution to the overall simplicity value",fill=TRUE)
cat("(they are I^2/p, J^2/q or K^2/r, resp., times the sum of the variances of squared values)",fill=TRUE)
out=list()
out$AS=AS
out$BT=BT
out$CU=CU
out$K=K
out$S=S
out$T=T
out$U=U
out$f=f
out$f1=f1
out$f2a=f2a
out$f2b=f2b
out$f2c=f2c
out$func=func
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.