Nothing
StatisBiplot <- function(X, InitTransform = "Standardize columns", dimens=2, SameVar=FALSE) {
mycall=match.call()
ContinuousDataTransform = c("Raw Data", "Substract the global mean", "Double centering", "Column centering", "Standardize columns", "Row centering",
"Standardize rows", "Divide by the column means and center", "Normalized residuals from independence")
if (is.numeric(InitialTransform))
InitialTransform = ContinuousDataTransform[InitialTransform]
ng = length(X) #Number of groups
StudyNames=names(X)
nri = matrix(0, ng, 1) #Number of rows in each group
nci = matrix(0, ng, 1) #Number of cols in each group
for (i in 1:ng) {
nri[i] = dim(X[[i]])[1]
nci[i] = dim(X[[i]])[2]
}
nr = nri[1]
if (sum(nri == nr) < ng)
stop("The number of rows must be the same in all ocassions")
if (SameVar){
nc = nci[1]
if (sum(nci == nc) < ng)
stop("The number of variables can not be the same in all ocassions (use SameVar=FALSE)")
}
# Extracting the names of the occasions
OccNames=names(X)
if (is.null(OccNames)) {
for (i in 1:ng) OccNames= c(OccNames, paste("Occasion_",i,sep=""))
}
# Initial transformation of data and calculation of statistics for the biplot
BiplotStatis=MultiTableStatistics(X)
BiplotStatis$call <- mycall
BiplotStatis$Type="Statis"
BiplotStatis$Title="Biplot induced by Statis"
BiplotStatis$Initial_Transformation=InitTransform
BiplotStatis$nrows=nr
BiplotStatis$ncols=sum(nci)
StatisRes=list()
StatisRes$Title="STATIS-ACT Biplot"
StatisRes$Type="STATIS-ACT"
StatisRes$NTables=ng
StatisRes$NRows=nr
if (SameVar)
StatisRes$NVars=nc
else
StatisRes$NVars=nci
StatisRes$RowLabels=rownames(X[[1]])
StatisRes$TableLabels=names(X)
if (SameVar) StatisRes$VarLabels=colnames(X[[1]])
X=MultiTableTransform(X, InitTransform = InitTransform)$X
#Calculation of the objects
Wt = list()
for (i in 1:ng) {
Wt[[i]] = X[[i]] %*% t(X[[i]])
}
# Calculation of the scalar products
P = matrix(0, ng, ng)
for (i in 1:ng) for (j in 1:ng) {
P[i, j] = sum(diag(Wt[[i]] %*% t(Wt[[j]])))
P[j, i] = P[i, j]
}
rownames(P)=StudyNames
colnames(P)=StudyNames
#RV: correlations among the occasions
RV = solve(diag(sqrt(diag(P)))) %*% P %*% solve(diag(sqrt(diag(P))))
rownames(RV)=StudyNames
colnames(RV)=StudyNames
# This information will be common to all the techniques, then a single
# function to plot is necassary
StatisRes$Data = X
StatisRes$Objects = Wt
StatisRes$ScalarProducts = P
StatisRes$RV = RV
Inter = svd(RV)
StatisRes$EigInter = Inter$d
names(StatisRes$EigInter)=paste("Dim",1:length(StudyNames))
StatisRes$InerInter=(Inter$d/sum(Inter$d))*100
names(StatisRes$InerInter)=paste("Dim",1:length(StudyNames))
StatisRes$InterStructure = -1 * Inter$u %*% diag(sqrt(Inter$d))
rownames(StatisRes$InterStructure)=StudyNames
colnames(StatisRes$InterStructure)=paste("Dim", 1:ng)
# Weigths for the Compromise
StatisRes$Weights = abs(Inter$u[, 1])/sum(abs(Inter$u[, 1]))
# Compromise
W=matrix(0,nr,nr)
for (i in 1:ng)
W=W + StatisRes$Weights[i]*Wt[[i]]
StatisRes$Compromise=W
# Euclidean Configuration of the compromise
Intra = svd(W)
r=sum(Intra$d>0.00001)
BiplotStatis$EigenValues=Intra$d[1:r]
BiplotStatis$Inertia=(Intra$d[1:r]/sum(diag(Intra$d[1:r])))*100
BiplotStatis$CumInertia = cumsum(BiplotStatis$Inertia)
BiplotStatis$alpha=1
BiplotStatis$Dimension=dimens
# Compromise-Consensus Coordinates and contributions
BiplotStatis$RowCoordinates = Intra$u[,1:dimens] %*% diag(sqrt(Intra$d[1:dimens]))
sf=apply((Intra$u[,1:r] %*% diag(sqrt(Intra$d[1:r])))^2,1,sum)
BiplotStatis$RowContributions=(diag(1/sf) %*% BiplotStatis$RowCoordinates^2)*100
rownames(BiplotStatis$RowCoordinates)=rownames(X[[1]])
colnames(BiplotStatis$RowCoordinates)=paste("Dim", 1:dimens)
rownames(BiplotStatis$RowContributions)=rownames(X[[1]])
colnames(BiplotStatis$RowContributions)=paste("Dim", 1:dimens)
# Trajectories for individuals (old kind)
trajin=list()
for (j in 1:ng)
trajin[[j]] = Wt[[j]] %*% Intra$u[,1:r] %*% diag(sqrt(1/Intra$d[1:r]))
StatisRes$TrajInd=list()
for (i in 1:nr){
Traj=NULL
for (j in 1:ng)
Traj=rbind(Traj , trajin[[j]][i,1:dimens])
rownames(Traj)=StudyNames
colnames(Traj)=paste("Dim", 1:dimens)
StatisRes$TrajInd[[i]]=Traj
}
names(StatisRes$TrajInd)=rownames(X[[1]])
# Squared cosines of the individual objects projected onto the compromise
contbt=list()
for (j in 1:ng){
sf=apply(trajin[[j]]^2,1,sum)
contbt[[j]] = (diag(1/sf) %*% trajin[[j]]^2)*100
}
StatisRes$ContribInd=list()
for (i in 1:nr){
Traj=NULL
for (j in 1:ng)
Traj=rbind(Traj , contbt[[j]][i,1:dimens])
rownames(Traj)=StudyNames
colnames(Traj)=paste("Dim", 1:dimens)
StatisRes$ContribInd[[i]]=Traj
}
names(StatisRes$ContribInd)=rownames(X[[1]])
# Coordinates and contributions of the variables on the biplot
# First we calculate the principal coordinates in order to calculate the contributions
trajvar=list()
for (j in 1:ng)
trajvar[[j]]=t(X[[j]]) %*% Intra$u[,1:r]
# Coordinates of the variables on the biplot
contvar=list()
for (j in 1:ng){
sf=apply(trajvar[[j]]^2,1,sum)
contvar[[j]] = (diag(1/sf) %*% trajvar[[j]]^2)*100
rownames(contvar[[j]])=paste(colnames(X[[j]]), "-", StudyNames[j], sep="")
}
# Standard coodinates to represent on the biplot (RMP-Biplot)
# Standard coordinates are needed to represent the biplot axes with graded scales
for (j in 1:ng){
trajvar[[j]]=t(X[[j]]) %*% Intra$u[,1:r] %*% diag(1/Intra$d[1:r])
rownames(trajvar[[j]])=paste(rownames(trajvar[[j]]), "-", StudyNames[j], sep="")
}
# This is just a reorganization of the column coordinates in order to represent them
# as a biplot using the standard procedure
BiplotStatis$ColCoordinates=NULL
for (j in 1:ng)
BiplotStatis$ColCoordinates=rbind(BiplotStatis$ColCoordinates, trajvar[[j]][,1:dimens])
colnames(BiplotStatis$ColCoordinates)<- paste("Dim ",1:dimens)
BiplotStatis$ColContributions=NULL
for (j in 1:ng)
BiplotStatis$ColContributions=rbind(BiplotStatis$ColContributions, contvar[[j]][,1:dimens])
colnames(BiplotStatis$ColContributions)<- paste("Dim ",1:dimens)
# Rescaling the biplot for better visualization. Used to be optional
# but now I do for every biplot.
# The row coordinates are multiplied for a constant and the coumn coordinates divided
# by the same constant. The scalar (inner product) in which biplot interpretation
# is based does not change, but the visuatization is much better.
# In a biplot with scales for each varibale, the rescaling does not affect
# the calculation of the scale marks.
# The rescaling is performed with all the dimensiond retained in the final solution
sca = sum(BiplotStatis$RowCoordinates^2)
scb = sum(BiplotStatis$ColCoordinates^2)
p=sum(nci)
sca = sca/nr
scb = scb/p
scf = sqrt(sqrt(scb/sca))
BiplotStatis$RowCoordinates = BiplotStatis$RowCoordinates * scf
BiplotStatis$ColCoordinates = BiplotStatis$ColCoordinates/scf
BiplotStatis$Scale_Factor=scf
for (i in 1:nr){
StatisRes$TrajInd[[i]]=StatisRes$TrajInd[[i]]*scf
}
BiplotStatis$ClusterType="us"
BiplotStatis$Clusters = as.factor(matrix(1,nrow(BiplotStatis$RowContributions), 1))
BiplotStatis$ClusterColors="blue"
BiplotStatis$ClusterNames="Cluster 1"
class(BiplotStatis) ="ContinuousBiplot"
# If the variables are the same for all the occasions or studies
# additional trajectories for the variables can be constructed
# selecting conveniently the var coordinates. This are original
# in our work, together with the simultaneous representation
# of rows and columns.
StatisRes$Biplot=BiplotStatis
StatisRes$SameIndiv = TRUE
StatisRes$SameVar = SameVar
if (SameVar){
p=nci[1]
StatisRes$TrajVar=list()
for (i in 1:p){
Traj=NULL
for (j in 1:ng)
Traj=rbind(Traj , trajvar[[j]][i,1:dimens]/scf)
rownames(Traj)=StudyNames
colnames(Traj)=paste("Dim", 1:dimens)
StatisRes$TrajVar[[i]]=Traj}
names(StatisRes$TrajVar)=colnames(X[[1]])
}
# Calculation of the Inertia of each study (or occasion) accounted for the STATIS solution.
# Comparing this with the inertia accounted by the biplot for each separate study
# we have an index of the goodness of the consensus (or compromise) in explaining the
# study.
StatisRes$Biplot$Structure = cor(StatisRes$Biplot$Non_Scaled_Data, StatisRes$Biplot$RowCoordinates)
class(StatisRes) = "StatisBiplot"
return(StatisRes)
}
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.