Nothing
pcamean <-
function(X,n,m,p,laba,labb,labc){
narg=nargs()
if (narg<5){
laba=paste("a",1:n,sep="")
labb=paste("b",1:m,sep="")
labc=paste("c",1:p,sep="")
}
X=as.matrix(X)
check=1
while (check==1){
cat("Over which mode do you want to compute means? Enter 1,2 or 3 (for A,B, or C)?",fill=TRUE)
aggregmode=scan("",n=1)
# aggregmode = mode over which means are computed
if (length(aggregmode)==0){
aggregmode=0
}
if (aggregmode==1){
check=0
y1=m
y2=p # B-mode x C-mode
check1=1
cat(paste("A", m,"by", p," matrix for modes B x C will be computed"),fill=TRUE)
while (check1==1)
{
cat("Which mode will serve as 'variables mode'? Enter 2 (=B) or 3 (=C):",fill=TRUE)
loadingsmode=scan("",n=1)
# loadingsmode = mode which is to be considered as the one
# for which you want to obtain simple structured loadings
# Note that the other mode then is the "scores" mode
# with scores normalized to unit sums of squares
if (length(loadingsmode)==0){
loadingsmode=0
}
if (loadingsmode==2){
check1=0
load='P'
cat("B will contain loadings, C will contain scores",fill=TRUE)
}
if (loadingsmode==3){
check1=0
load='Q'
cat("C will contain loadings, B will contain scores",fill=TRUE)
} else{
if ((loadingsmode!=2) & (loadingsmode!=3)){
cat(" ",fill=TRUE)
cat("Error! Select a proper number!",fill=TRUE)
cat(" ",fill=TRUE)
}
}
}
}
if (aggregmode==2){
check=0
y1=p
y2=n # C-mode x A-mode
X=permnew(X,n,m,p)
check1=1
while (check1==1){
cat(paste("A", p," by", n," matrix for modes C x A will be computed"),fill=TRUE)
cat("Which mode will serve as 'variables mode'? Enter 1 (=A) or 3 (=C):",fill=TRUE)
loadingsmode=scan("",n=1)
if (length(loadingsmode)==0){
loadingsmode=0
}
if (loadingsmode==3){
check1=0
load='P'
cat("C will contain loadings, A will contain scores",fill=TRUE)
}
if (loadingsmode==1){
check1=0
load='Q'
cat("A will contain loadings, C will contain scores",fill=TRUE)
} else{
if ((loadingsmode!=1) & (loadingsmode!=3)){
cat(" ",fill=TRUE)
cat("Error! Select a proper number!",fill=TRUE)
cat(" ",fill=TRUE)
}
}
}
}
if (aggregmode==3){
check=0
X=permnew(X,n,m,p)
X=permnew(X,m,p,n)
y1=n
y2=m # A-mode x B-mode
check1=1
while (check1==1){
cat(paste("A", n," by", m," matrix for modes A x B will be computed"),fill=TRUE)
cat("Which mode will serve as 'variables mode'? Enter 1 (=A) or 2 (=B):",fill=TRUE)
loadingsmode=scan("",n=1)
if (length(loadingsmode)==0){
loadingsmode=0
}
if (loadingsmode==1){
check1=0
load='P'
cat("A will contain loadings, B will contain scores",fill=TRUE)
}
if (loadingsmode==2){
check1=0
load='Q'
cat("B will contain loadings, A will contain scores",fill=TRUE)
} else{
if ((loadingsmode!=1) & (loadingsmode!=2)){
cat(" ",fill=TRUE)
cat("Error! Select a proper number!",fill=TRUE)
cat(" ",fill=TRUE)
}
}
}
}
if ((aggregmode!=1) & (aggregmode!=2) & (aggregmode!=3)){
cat(" ",fill=TRUE)
cat("Error! Select a proper number!",fill=TRUE)
cat(" ",fill=TRUE)
}
}
mm=t(SUM(X)$mc)
Y=matrix(mm,y1,y2) # P-mode x Q-mode
if (nrow(Y)==n){
rownames(Y)=laba
}
if (ncol(Y)==n){
colnames(Y)=laba
}
if (nrow(Y)==m){
rownames(Y)=labb
}
if (ncol(Y)==m){
colnames(Y)=labb
}
if (nrow(Y)==p){
rownames(Y)=labc
}
if (ncol(Y)==p){
colnames(Y)=labc
}
if (sum(Y^2)<1e-10){
cat("Means to analyze are exteremely small, probably data have been centered across the",fill=TRUE)
cat("same mode as over which means are computed, which makes no sense.",fill=TRUE)
cat("Therefore, PCA of means is not carried out!",fill=TRUE)
A1=NULL
A2=NULL
B1=NULL
B2=NULL
C1=NULL
C2=NULL
ev=NULL
} else{
cat("If you want to standardize 'variables', type '1':",fill=TRUE)
cc=scan("",n=1)
if (length(cc)==0){
cc=0
}
if (cc==1){
if (load=='Q'){
Y=nrm2(Cc(Y))
}
if (load=='P'){
Y=t(nrm2(Cc(t(Y))))
}
cat("'Variables' have been standardized",fill=TRUE)
} else{
cat("'Variables' have not been standardized",fill=TRUE)
}
}
if (sum(Y^2)>1e-10){
SVD=svd(Y) # P for rows Q for columns
ev=SVD$d^2
rr=length(ev)
cat("Eigenvalues and cumulative percentages (Scree plot in Figure)",fill=TRUE)
mat=matrix(,length(ev),2)
mat[,1]=ev
mat[,2]=cumsum(ev)/sum(ev)*100
labComp=paste("Comp.",1:rr,sep="")
names(ev)=labComp
rownames(mat)=labComp
colnames(mat)=c("Eigenvalue","Fit(%)")
print(round(mat,digits=2))
plot(t(c(1:rr)),ev,type="l",xlab="",ylab="")
rMax=min(nrow(Y),ncol(Y))
while(check==0){
cat("How many components do you want to maintain for rotation?",fill=TRUE)
r=scan("",n=1)
# check on the number of components
if (length(r)==0){
r=0
}
if ((r>1) & (r<=rMax) & ((floor(r)-r)==0)){
check=1
labComp=paste("Comp.",1:r,sep="")
if (load=='P'){
# varimax of P-mode loadings
NRMVAR=normvari(SVD$u[,1:r]%*%diag(SVD$d,nrow=rMax)[1:r,1:r])
Q1=SVD$v[,1:r]%*%NRMVAR$T
sg=diag(sign(colSums(NRMVAR$B)),nrow=ncol(NRMVAR$B)) # reflect if necessary
P1=NRMVAR$B%*%sg
Q1=Q1%*%sg
# Harris-Kaiser of P-mode loadings
NRMVAR2=normvari(SVD$u[,1:r])
Q2=SVD$v[,1:r]%*%diag(SVD$d,nrow=rMax)[1:r,1:r]%*%NRMVAR2$T
Ds=diag(colSums(Q2^2),nrow=ncol(Q2))
P2=NRMVAR2$B%*%Ds^.5
Q2=Q2%*%solve(Ds^.5)
sg=diag(sign(colSums(P2)),nrow=ncol(P2)) # reflect if necessary
P2=P2%*%sg
Q2=Q2%*%sg
cat("A1, B1 and/or C1 contain solution based on varimax rotation of loadings",fill=TRUE)
cat("A2, B2 and/or C2 contain solution based on oblique 'HKIC' orthomax rotation of loadings",fill=TRUE)
}
if (load=='Q'){
# varimax of Q-mode loadings
NRMVAR=normvari(SVD$v[,1:r]%*%diag(SVD$d,nrow=rMax)[1:r,1:r])
P1=SVD$u[,1:r]%*%NRMVAR$T
sg=diag(sign(colSums(NRMVAR$B)),nrow=ncol(NRMVAR$B)) # reflect if necessary
P1=P1%*%sg
Q1=NRMVAR$B%*%sg
# Harris-Kaiser of Q-mode loadings
NRMVAR2=normvari(SVD$v[,1:r])
P2=SVD$u[,1:r]%*%diag(SVD$d,nrow=rMax)[1:r,1:r]%*%NRMVAR2$T
Ds=diag(colSums(P2^2),nrow=ncol(P2))
Q2=NRMVAR2$B%*%Ds^.5
P2=P2%*%solve(Ds^.5)
sg=diag(sign(colSums(Q2)),nrow=ncol(Q2)) # reflect if necessary
P2=P2%*%sg
Q2=Q2%*%sg
cat("A1, B1 and/or C1 contain solution based on varimax rotation of loadings",fill=TRUE)
cat("A2, B2 and/or C2 contain solution based on oblique 'HKIC' orthomax rotation of loadings",fill=TRUE)
}
} else{
if (r==1){
labComp=c("Comp.1")
check=1
if (load=='P'){
P1=SVD$u[,1]*SVD$d[1]
sg=sign(sum(P1)) # reflect if necessary
P1=as.matrix(P1)*sg
Q1=as.matrix(SVD$v[,1])*sg
P2=P1
Q2=Q1
cat("A1, B1 and/or C1 contain principal component solution",fill=TRUE)
}
if (load=='Q'){
Q1=SVD$v[,1]*SVD$d[1]
sg=sign(sum(Q1)) # reflect if necessary
Q1=as.matrix(Q1)*sg
P1=as.matrix(SVD$u[,1])*sg
P2=P1
Q2=Q1
cat("A1=A2, B1=B2 and/or C1=C2 contain principal component solution",fill=TRUE)
}
} else{
cat(" ",fill=TRUE)
cat("Error! Select a proper number!",fill=TRUE)
cat(" ",fill=TRUE)
}
}
}
if (aggregmode==1 ){
A1=NULL
A2=NULL
B1=P1
B2=P2
C1=Q1
C2=Q2
}
if (aggregmode==2){
A1=Q1
A2=Q2
B1=NULL
B2=NULL
C1=P1
C2=P2
}
if (aggregmode==3){
A1=P1
A2=P2
B1=Q1
B2=Q2
C1=NULL
C2=NULL
}
if (identical(A1,NULL)==FALSE){
rownames(A1)=laba
rownames(A2)=laba
colnames(A1)=labComp
colnames(A2)=labComp
}
if (identical(B1,NULL)==FALSE){
rownames(B1)=labb
rownames(B2)=labb
colnames(B1)=labComp
colnames(B2)=labComp
}
if (identical(C1,NULL)==FALSE){
rownames(C1)=labc
rownames(C2)=labc
colnames(C1)=labComp
colnames(C2)=labComp
}
}
out=list()
out$Y=Y
out$ev=ev
out$A1=A1
out$A2=A2
out$B1=B1
out$B2=B2
out$C1=C1
out$C2=C2
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.