Nothing
splithalfT3 <-
function(X,n,m,p,r1,r2,r3,centopt,normopt,renormmode,wa_rel,wb_rel,wc_rel,addanal,conv,laba,labb,labc){
cat("This procedure performs a SPLIT-HALF analysis on X",fill=TRUE)
cat("NOTE:",fill=TRUE)
cat("In SPLIT-HALF analysis, the A-mode is taken as 'replication mode'",fill=TRUE)
cat("(which means that A-mode entities are considered a random sample",fill=TRUE)
cat("if this does not make sense, you should rearrange your data so that the A-mode is a replication mode)",fill=TRUE)
cat("The splitting into two halves can be done randomly (default), or into odd vs. even sequence numbers",fill=TRUE)
narg=nargs()
if (narg<16){
laba=paste("a",1:n,sep="")
labb=paste("b",1:m,sep="")
labc=paste("c",1:p,sep="")
}
X=as.matrix(X)
cat("If you prefer odd/even split, specify '1':",fill=TRUE)
nsplit=scan("",n=1)
if (length(nsplit)==0){
nsplit=0
}
if (nsplit==1){ #wi=[1:2:n 2:2:n]
check=0
wi=c(seq(1,n,2), seq(2,n,2))
cat("Splitting has been done into odd vs. even sequence numbers",fill=TRUE)
} else{
# create random splits
w=matrix(runif(n*1,0,1),n,1)
wi=ord(w)$a
cat("Splitting has been done randomly",fill=TRUE)
}
n1=ceiling(n/2)
n2=n-n1
X1=X[wi[1:n1],]
X2=X[wi[(n1+1):n],]
cat("You will now enter the split half procedure with the options specified so far",fill=TRUE)
cat("However, you may want to modify certain choices here (rather than rerunning the full Tucker3)",fill=TRUE)
cat("Do you want to modify certain choices? If so, specify '1':",fill=TRUE)
ccc=scan("",n=1)
if (length(ccc)==0){
ccc=0
}
if (ccc==1){
cat("Centering:",fill=TRUE)
cat(" 0 = none (default)", fill=TRUE)
cat(" 1= across A-mode",fill=TRUE)
cat(" 2= across B-mode",fill=TRUE)
cat(" 3= across C-mode",fill=TRUE)
cat(" 12= across A-mode and across B-mode",fill=TRUE)
cat(" 13= across A-mode and across C-mode",fill=TRUE)
cat(" 23= across B-mode and across C-mode",fill=TRUE)
cat(paste("Your centering choice was",centopt),fill=TRUE)
cat("If you want to change it, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
check=0
while(check==0){
cat("Specify centering option:",fill=TRUE)
centopt=scan("",n=1)
if (length(centopt)==0){
centopt=0
}
if((centopt!=0) & ((centopt!=1) & (centopt!=2) & (centopt!=3) & (centopt!=12) & (centopt!=13) & (centopt!=23))){
cat(" ",fill=TRUE)
cat("Error! Make a proper choice for centering the data")
cat(" ",fill=TRUE)
} else{
check=1
}
}
}
cat("Normalizing:",fill=TRUE)
cat(" 0= none (default)", fill=TRUE)
cat(" 1= within A-mode", fill=TRUE)
cat(" 2= within B-mode",fill=TRUE)
cat(" 3= within C-mode",fill=TRUE)
cat(paste("Your normalizing choice was",normopt),fill=TRUE)
cat("If you want to change it, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
check=0
while(check==0){
cat("Specify normalizing option",fill=TRUE)
normopt=scan("",n=1)
if (length(normopt)==0){
normopt=0
}
if ((normopt!=0) & ((normopt!=1) & (normopt!=2) & (normopt!=3))){
cat(" ",fill=TRUE)
cat("Error! Make a proper choice for normalizing the data")
cat(" ",fill=TRUE)
} else{
check=1
}
}
}
if (renormmode==1){
cat("Your choice was to renormalize A-mode",fill=TRUE)
}
if (renormmode==2){
cat("Your choice was to renormalize B-mode",fill=TRUE)
}
if (renormmode==3){
cat("Your choice was to renormalize C-mode",fill=TRUE)
} else{
if ((renormmode!=1) & (renormmode!=2) & (renormmode!=3)){
cat("Your choice was not to renormalize solution",fill=TRUE)
}
}
cat("If you want to change it, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
check2=0
while (check2==0){
cat("Which mode do you want to renormalize? Enter 1 (=A), 2 (=B) or 3 (=C):", fill=TRUE)
renormmode=scan("",n=1)
if (length(renormmode)==0){
renormmode=0
}
if ((renormmode!=1) &(renormmode!=2) & (renormmode!=3)){
cat(" ",fill=TRUE)
cat("Error! Make a proper choice for normalizing the data")
cat(" ",fill=TRUE)
} else{
check2=1
}
}
}
cat(paste("Numbers of components for A, B and C were: ",r1,", ",r2,", ",r3,sep=""),fill=TRUE)
cat("If you want to change them, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
cat("How many A-mode components do you want to use?",fill=TRUE)
check=0
while (check==0){
r1=scan("",n=1)
if (length(r1)==0){
r1=0
}
if ((r1==0) | ((floor(r1)-r1)!=0)){
cat(" ",fill=TRUE)
cat("Error! How many A-mode components do you want to use?",fill=TRUE)
cat(" ",fill=TRUE)
} else{
check=1
}
}
cat("How many B-mode components do you want to use?",fill=TRUE)
check=0
while (check==0){
r2=scan("",n=1)
if (length(r2)==0){
r2=0
}
if ((r2==0) | ((floor(r2)-r2)!=0)){
cat(" ",fill=TRUE)
cat("Error! How many B-mode components do you want to use?",fill=TRUE)
cat(" ",fill=TRUE)
} else{
check=1
}
}
cat("How many C-mode components do you want to use?",fill=TRUE)
check=0
while (check==0){
r3=scan("",n=1)
if (length(r3)==0){
r3=0
}
if ((r3==0) | ((floor(r3)-r3)!=0)){
cat(" ",fill=TRUE)
cat("Error! How many C-mode components do you want to use?",fill=TRUE)
cat(" ",fill=TRUE)
} else{
check=1
}
}
}
cat(paste("Your choice was to use a convergence criterion equal to",conv),fill=TRUE)
cat("If you want to change it, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
cat("Specify convergence criterion (default=1e-6)", fill=TRUE)
conv=scan("",n=1)
if (length(conv)==0){
conv=1e-6
}
}
cat(paste("Your choice was to use",addanal,"additional random starts in the analysis"),fill=TRUE)
cat("If you want to change it, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
check=0
cat("How many additional runs do you want to use?",fill=TRUE)
while (check==0){
addanal=scan("",n=1)
if (length(addanal)==0){
addanal=0
check=1
}
if ((floor(addanal)-addanal)!=0){
cat(" ",fill=TRUE)
cat("Error! How many additional runs do you want to use?",fill=TRUE)
cat(" ",fill=TRUE)
} else{
check=1
}
}
}
cat(paste("Relative simplicity rotation weights for A, B and C are: ",wa_rel,", ",wb_rel,", ",wc_rel,sep=""),fill=TRUE)
cat("If you want to change them, specify '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
cat("Specify relative weight for simplicity of A-mode (default=0):",fill=TRUE)
wa_rel=scan("",n=1)
cat("Specify relative weight for simplicity of B-mode (default=0):",fill=TRUE)
wb_rel=scan("",n=1)
cat("Specify relative weight for simplicity of C-mode (default=0):",fill=TRUE)
wc_rel=scan("",n=1)
if (length(wa_rel)==0){
wa_rel=0
}
if (length(wb_rel)==0){
wb_rel=0
}
if (length(wc_rel)==0){
wc_rel=0
}
}
}
cat("Analysis of FULL DATA",fill=TRUE)
# Full data analysis (reference)
# Full data Preprocessing
cc=centopt
if ((cc==1) | (cc==12) | (cc==13)){
X=cent3(X,n,m,p,1)
cat("X has been centered across A-mode",fill=TRUE)
}
if ((cc==2) | (cc==12) | (cc==23)){
X=cent3(X,n,m,p,2)
cat("X has been centered across B-mode",fill=TRUE)
}
if ((cc==3) | (cc==13) | (cc==23)){
X=cent3(X,n,m,p,3)
cat("X has been centered across C-mode",fill=TRUE)
} else{
cat("X has not been centered",fill=TRUE)
}
cc=normopt
if (cc==1){
X=norm3(X,n,m,p,1)
cat("X has been normalized within A-mode",fill=TRUE)
}
if (cc==2){
X=norm3(X,n,m,p,2)
cat("X has been normalized within B-mode",fill=TRUE)
}
if (cc==3){
X=norm3(X,n,m,p,3)
cat("X has been normalized within C-mode",fill=TRUE)
}
if ((cc!=1) & (cc!=2) & (cc!=3)){
cat("X has not been normalized",fill=TRUE)
}
# Full data analysis
Tuck3=T3func(X,n,m,p,r1,r2,r3,0,1e-6)
A=Tuck3$A
B=Tuck3$B
C=Tuck3$C
H=Tuck3$H
f=Tuck3$f
iter=Tuck3$iter
fp=Tuck3$fp
La=Tuck3$La
Lb=Tuck3$Lb
Lc=Tuck3$Lc
func=vector("numeric",length=1+addanal)
names(func)=paste("Start n.",1:(1+addanal),sep="")
func=Tuck3$fp
print(func)
for (run in 1:addanal){
cat(paste("Run no.",run+1,sep=" "),fill=TRUE)
Tuck3b=T3func(X,n,m,p,r1,r2,r3,1,1e-6)
func[run+1]=Tuck3b$fp
if (Tuck3b$fp>1.0001*Tuck3$fp){ # if fit more than .01% better is found, replace solution
A=Tuck3b$A
B=Tuck3b$B
C=Tuck3b$C
H=Tuck3b$H
f=Tuck3b$f
iter=Tuck3b$iter
fp=Tuck3b$fp
La=Tuck3b$La
Lb=Tuck3b$Lb
Lc=Tuck3b$Lc
}
}
if (addanal>=1){
cat("Fit % values from all runs:",fill=TRUE)
print(round(t(func), digits=2))
}
RNsol=renormsolT3(A,B,C,H,renormmode)
# Full data Rotation
VARM=varimcoco(RNsol$A,RNsol$B,RNsol$C,RNsol$H,wa_rel,wb_rel,wc_rel)
AS=VARM$AS
BT=VARM$BT
CU=VARM$CU
K=VARM$K
cat("Backnormalize A,B,C, and H",fill=TRUE)
if (renormmode==1){
Ds=diag(SUM(AS)$col^.5,nrow=r1)
AS=AS%*%solve(Ds)
K=Ds%*%K
}
if (renormmode==2){
Ds=diag(SUM(BT)$col^.5,nrow=r2)
BT=BT%*%solve(Ds)
K=K%*%kronecker(diag(r3),Ds)
}
if (renormmode==3){
Ds=diag(SUM(CU)$col^.5,nrow=r3)
CU=CU%*%solve(Ds)
K=K%*%kronecker(Ds,diag(r2))
}
Afull=AS
Bfull=BT
Cfull=CU
Kfull=K
# Split 1 data analysis (reference)
# Split 1 Preprocessing
cat("Analysis of SPLIT 1",fill=TRUE)
n_orig=n
n=n1
X_orig=X
X=X1
cc=centopt
if ((cc==1) | (cc==12) | (cc==13)){
X=cent3(X,n,m,p,1)
cat("X has been centered across A-mode",fill=TRUE)
}
if ((cc==2) | (cc==12) | (cc==23)){
X=cent3(X,n,m,p,2)
cat("X has been centered across B-mode",fill=TRUE)
}
if ((cc==3) | (cc==13) | (cc==23)){
X=cent3(X,n,m,p,3)
cat("X has been centered across C-mode",fill=TRUE)
}
cc=normopt
if (cc==1 ){
X=norm3(X,n,m,p,1)
cat("X has been normalized within A-mode",fill=TRUE)
}
if (cc==2){
X=norm3(X,n,m,p,2)
cat("X has been normalized within B-mode",fill=TRUE)
}
if (cc==3){
X=norm3(X,n,m,p,3)
cat("X has been normalized within C-mode",fill=TRUE)
}
if ((cc!=1) & (cc!=2) & (cc!=3)){
cat("X has not been normalized",fill=TRUE)
}
Xs1=X
# Split 1 analysis
Tuck3=T3func(X,n,m,p,r1,r2,r3,0,1e-6)
A=Tuck3$A
B=Tuck3$B
C=Tuck3$C
H=Tuck3$H
f=Tuck3$f
iter=Tuck3$iter
fp=Tuck3$fp
La=Tuck3$La
Lb=Tuck3$Lb
Lc=Tuck3$Lc
func=vector("numeric",length=1+addanal)
names(func)=paste("Start n.",1:(1+addanal),sep="")
func=Tuck3$fp
for (run in 1:addanal){
cat(paste("Run no.",run+1,sep=" "),fill=TRUE)
Tuck3b=T3func(X,n,m,p,r1,r2,r3,1,1e-6)
func[run+1]=Tuck3b$fp
if (Tuck3b$fp>1.0001*Tuck3$fp){ # if fit more than .01% better is found, replace solution
A=Tuck3b$A
B=Tuck3b$B
C=Tuck3b$C
H=Tuck3b$H
f=Tuck3b$f
iter=Tuck3b$iter
fp=Tuck3b$fp
La=Tuck3b$La
Lb=Tuck3b$Lb
Lc=Tuck3b$Lc
}
}
if (addanal>=1){
cat("Fit % values from all runs:",fill=TRUE)
print(round(t(func), digits=2))
}
cat("No Postprocessing used in splits, taken care off by reference rotation!",fill=TRUE)
As1=A
Bs1=B
Cs1=C
Ks1=H
# Split 2 data analysis (reference)
# Split 2 Preprocessing
cat("Analysis of SPLIT 2",fill=TRUE)
n=n2
X=X2
cc=centopt
if ((cc==1) | (cc==12) | (cc==13)){
X=cent3(X,n,m,p,1)
cat("X has been centered across A-mode",fill=TRUE)
}
if ((cc==2) | (cc==12) | (cc==23)){
X=cent3(X,n,m,p,2)
cat("X has been centered across B-mode",fill=TRUE)
}
if ((cc==3) | (cc==13) | (cc==23)){
X=cent3(X,n,m,p,3)
cat("X has been centered across C-mode",fill=TRUE)
}
cc=normopt
if (cc==1){
X=norm3(X,n,m,p,1)
cat("X has been normalized within A-mode",fill=TRUE)
}
if (cc==2){
X=norm3(X,n,m,p,2)
cat("X has been normalized within B-mode",fill=TRUE)
}
if (cc==3){
X=norm3(X,n,m,p,3)
cat("X has been normalized within C-mode",fill=TRUE)
}
if ((cc!=1) & (cc!=2) & (cc!=3)){
cat("X has not been normalized",fill=TRUE)
}
Xs2=X
# Split 2 analysis
Tuck3=T3func(X,n,m,p,r1,r2,r3,0,1e-6)
A=Tuck3$A
B=Tuck3$B
C=Tuck3$C
H=Tuck3$H
f=Tuck3$f
iter=Tuck3$iter
fp=Tuck3$fp
La=Tuck3$La
Lb=Tuck3$Lb
Lc=Tuck3$Lc
func=vector("numeric",length=1+addanal)
names(func)=paste("Start n.",1:(1+addanal),sep="")
func=Tuck3$fp
for (run in 1:addanal){
cat(paste("Run no.",run+1,sep=" "),fill=TRUE)
Tuck3b=T3func(X,n,m,p,r1,r2,r3,1,1e-6)
func[run+1]=Tuck3b$fp
if (Tuck3b$fp>1.0001*Tuck3$fp){ # if fit more than .01% better is found, replace solution
A=Tuck3b$A
B=Tuck3b$B
C=Tuck3b$C
H=Tuck3b$H
f=Tuck3b$f
iter=Tuck3b$iter
fp=Tuck3b$fp
La=Tuck3b$La
Lb=Tuck3b$Lb
Lc=Tuck3b$Lc
}
}
if (addanal>=1){
cat("Fit % values from all runs:",fill=TRUE)
print(round(t(func), digits=2))
}
cat("No Postprocessing used in splits, taken care off by reference rotation!",fill=TRUE)
As2=A
Bs2=B
Cs2=C
Ks2=H
n=n_orig
X=X_orig
# Compute stability indices:
Ss1=solve(t(As1)%*%As1)%*%t(As1)%*%Afull[wi[1:n1],]
Ss2=solve(t(As2)%*%As2)%*%t(As2)%*%Afull[wi[(n1+1):n],]
Ts1=solve(t(Bs1)%*%Bs1)%*%t(Bs1)%*%Bfull
Ts2=solve(t(Bs2)%*%Bs2)%*%t(Bs2)%*%Bfull
Us1=solve(t(Cs1)%*%Cs1)%*%t(Cs1)%*%Cfull
Us2=solve(t(Cs2)%*%Cs2)%*%t(Cs2)%*%Cfull
As1=As1%*%Ss1
As2=As2%*%Ss2
Bs1=Bs1%*%Ts1
Bs2=Bs2%*%Ts2
Cs1=Cs1%*%Us1
Cs2=Cs2%*%Us2
Ks1=solve(Ss1)%*%Ks1%*%solve(kronecker(t(Us1),t(Ts1)))
Ks2=solve(Ss2)%*%Ks2%*%solve(kronecker(t(Us2),t(Ts2)))
labCompA=paste("A",1:r1,sep="")
labCompB=paste("B",1:r2,sep="")
labCompC=paste("C",1:r3,sep="")
# make labels for columns of core
str=noquote(vector(mode="character",length=r3*r2))
i=1
for (k in 1:r3){
for (j in 1:r2){
str[i]=noquote(paste(" B",as.character(j),"xC",as.character(k),sep=""))
i=i+1
}
}
corelabelstr=str
cat("RESULTS stability analysis",fill=TRUE)
cat("Both splits are analyzed and rotated towards solution for full data",fill=TRUE)
cat("One can compare complete outcomes for different splits",fill=TRUE)
cat("or inspect correlation/congruence coefficients between corresponding columns of component matrices",fill=TRUE)
cat("Split-half solutions for B",fill=TRUE)
cat("B's next to each other, separated by column of 0's",fill=TRUE)
X=round(cbind(Bs1,0,Bs2),digits=2)
rownames(X)=labb
colnames(X)=c(labCompB,"-",labCompB)
print(X)
cat("Split-half solutions for C",fill=TRUE)
cat("C's next to each other, separated by column of 0's",fill=TRUE)
X=round(cbind(Cs1,0,Cs2),digits=2)
rownames(X)=labc
colnames(X)=c(labCompC,"-",labCompC)
print(X)
cat("Congruences for A in splits and in appropriate part of Afull",fill=TRUE)
X=round(cbind(diag(phi(Afull[wi[1:n1],],As1)),diag(phi(Afull[wi[(n1+1):n],],As2))),digits=2)
rownames(X)=labCompA
colnames(X)=c("SPL1","SPL2")
print(X)
cat("Correlations for A in splits and in appropriate part of Afull",fill=TRUE)
X=round(cbind(diag(phi(Cc(Afull[wi[1:n1],]),Cc(As1))),diag(phi(Cc(Afull[wi[(n1+1):n],]),Cc(As2)))),digits=2)
rownames(X)=labCompA
colnames(X)=c("SPL1","SPL2")
print(X)
cat("Congruence values for B-mode component matrix",fill=TRUE)
X=round(diag(phi(Bs1,Bs2)),digits=2)
names(X)=c(labCompB)
print(X)
cat("Congruence values for C-mode component matrix",fill=TRUE)
X=round(diag(phi(Cs1,Cs2)),digits=2)
names(X)=c(labCompC)
print(X)
cat("Relatively strong stability check of Core:",fill=TRUE)
cat("(simply comparing core matrices for two splits)",fill=TRUE)
cat("Core for split1",fill=TRUE)
Ks1=as.matrix(Ks1,digits=2)
rownames(Ks1)=labCompA
colnames(Ks1)=corelabelstr
print(round(Ks1,digits=2))
cat("Core for split2",fill=TRUE)
Ks2=as.matrix(Ks2,digits=2)
rownames(Ks2)=labCompA
colnames(Ks2)=corelabelstr
print(round(Ks2,digits=2))
cat("Weaker but Sufficient stability check of Core:",fill=TRUE)
cat("(computing cores in two splits, using full data solutions for A,B and C)",fill=TRUE)
A1=Afull[wi[1:n1],]
A2=Afull[wi[(n1+1):n],]
Hss1=solve(t(A1)%*%A1)%*%t(A1)%*%Xs1
Hss1=solve(t(Bfull)%*%Bfull)%*%t(Bfull)%*%permnew(Hss1,r1,m,p)
Hss1=solve(t(Cfull)%*%Cfull)%*%t(Cfull)%*%permnew(Hss1,r2,p,r1)
Hss1=permnew(Hss1,r3,r1,r2)
Hss2=solve(t(A2)%*%A2)%*%t(A2)%*%Xs2
Hss2=solve(t(Bfull)%*%Bfull)%*%t(Bfull)%*%permnew(Hss2,r1,m,p)
Hss2=solve(t(Cfull)%*%Cfull)%*%t(Cfull)%*%permnew(Hss2,r2,p,r1)
Hss2=permnew(Hss2,r3,r1,r2)
cat("Core for split1",fill=TRUE)
Hss1=as.matrix(Hss1)
rownames(Hss1)=labCompA
colnames(Hss1)=corelabelstr
print(round(Hss1,digits=2))
cat("Core for split2",fill=TRUE)
Hss2=as.matrix(Hss2)
rownames(Hss2)=labCompA
colnames(Hss2)=corelabelstr
print(round(Hss2,digits=2))
cat("Core for full data",fill=TRUE)
Kfull=as.matrix(Kfull)
rownames(Kfull)=labCompA
colnames(Kfull)=corelabelstr
print(round(Kfull,digits=2))
colnames(As1)=labCompA
colnames(As2)=labCompA
colnames(Afull)=labCompA
colnames(Bs1)=labCompB
colnames(Bs2)=labCompB
colnames(Bfull)=labCompB
colnames(Cs1)=labCompC
colnames(Cs2)=labCompC
colnames(Cfull)=labCompC
rownames(As1)=wi[1:n1]
rownames(As2)=wi[(n1+1):n]
rownames(Afull)=laba
rownames(Bs1)=labb
rownames(Bs2)=labb
rownames(Bfull)=labb
rownames(Cs1)=labc
rownames(Cs2)=labc
rownames(Cfull)=labc
out=list()
out$Afull=Afull
out$As1=As1
out$As2=As2
out$Bfull=Bfull
out$Bs1=Bs1
out$Bs2=Bs2
out$Cfull=Cfull
out$Cs1=Cs1
out$Cs2=Cs2
out$Kfull=Kfull
out$Ks1=Ks1
out$Ks2=Ks2
out$Kss1=Hss1
out$Kss2=Hss2
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.