test.data <- Harman74.cor$cov
my.vss <- VSS(test.data) #suggests that 4 factor complexity two solution is optimal
print(my.vss[,1:12],digits =2)
VSS.plot(my.vss) #see graphic window for a plot
#produces this output
# dof chisq prob sqresid fit cfit.1 cfit.2 cfit.3 cfit.4 cfit.5 cfit.6 cfit.7
#1 252 4583 0.0e+00 17.2 0.79 0.79 0.00 0.00 0.00 0.00 0.00 0.00
#2 229 3105 0.0e+00 12.9 0.84 0.75 0.84 0.00 0.00 0.00 0.00 0.00
#4 186 1689 2.3e-240 8.0 0.90 0.66 0.87 0.90 0.90 0.00 0.00 0.00
#5 166 1398 9.3e-194 7.3 0.91 0.68 0.86 0.90 0.91 0.91 0.00 0.00
#6 147 1183 2.9e-161 6.5 0.92 0.53 0.83 0.88 0.91 0.92 0.92 0.00
#7 129 1002 5.8e-135 5.7 0.93 0.47 0.78 0.88 0.91 0.92 0.93 0.93
#8 112 803 5.3e-105 5.3 0.94 0.49 0.76 0.86 0.90 0.92 0.93 0.93
#compare the above solution to a "varimax" rotated solution which suggests 1 factor (g)
my.vss <- VSS(test.data,rotate="varimax") #suggests that 1 factor complexity one solution is optimal
print(my.vss[,1:14],digits =2)
VSS.plot(my.vss) #see graphic window for a plot
# dof chisq prob sqresid fit cfit.1 cfit.2 cfit.3 cfit.4 cfit.5 cfit.6 cfit.7 cfit.8 cresidual.1
#1 252 4583 0.0e+00 17.2 0.79 0.79 0.00 0.00 0.0 0.00 0.00 0.00 0.00 17
#2 229 3105 0.0e+00 12.9 0.84 0.54 0.84 0.00 0.0 0.00 0.00 0.00 0.00 38
#3 207 2193 0.0e+00 10.1 0.88 0.46 0.79 0.88 0.0 0.00 0.00 0.00 0.00 45
#4 186 1689 2.3e-240 8.0 0.90 0.42 0.73 0.87 0.9 0.00 0.00 0.00 0.00 48
#5 166 1398 9.3e-194 7.3 0.91 0.40 0.70 0.86 0.9 0.91 0.00 0.00 0.00 50
#6 147 1183 2.9e-161 6.5 0.92 0.39 0.69 0.86 0.9 0.92 0.92 0.00 0.00 51
#7 129 1002 5.8e-135 5.7 0.93 0.39 0.70 0.84 0.9 0.92 0.93 0.93 0.00 50
#8 112 803 5.3e-105 5.3 0.94 0.39 0.69 0.83 0.9 0.92 0.93 0.93 0.94 50
## The function is currently defined as
function (x,n=8,rotate="none",diagonal=FALSE,pc="pa",n.obs=1000,...) #apply the Very Simple Structure Criterion for up to n factors on data set x
#x is a data matrix
#n is the maximum number of factors to extract (default is 8)
#rotate is a string "none" or "varimax" for type of rotation (default is "none"
#diagonal is a boolean value for whether or not we should count the diagonal (default=FALSE)
# ... other parameters for factanal may be passed as well
#e.g., to do VSS on a covariance/correlation matrix with up to 8 factors and 3000 cases:
#VSS(covmat=msqcovar,n=8,rotate="none",n.obs=3000)
{ #start Function definition
#first some preliminary functions
#complexrow sweeps out all except the c largest loadings
#complexmat applies complexrow to the loading matrix
complexrow <- function(x,c) #sweep out all except c loadings
{ n=length(x) #how many columns in this row?
temp <- x #make a temporary copy of the row
x <- rep(0,n) #zero out x
for (j in 1:c)
{
locmax <- which.max(abs(temp)) #where is the maximum (absolute) value
x[locmax] <- sign(temp[locmax])*max(abs(temp)) #store it in x
temp[locmax] <- 0 #remove this value from the temp copy
}
return(x) #return the simplified (of complexity c) row
}
complexmat <- function(x,c) #do it for every row (could tapply somehow?)
{
nrows <- dim(x)[1]
ncols <- dim(x)[2]
for (i in 1:nrows)
{x[i,] <- complexrow(x[i,],c)} #simplify each row of the loading matrix
return(x)
}
#now do the main Very Simple Structure routine
complexfit <- array(0,dim=c(n,n)) #store these separately for complex fits
complexresid <- array(0,dim=c(n,n))
vss.df <- data.frame(dof=rep(0,n),chisq=0,prob=0,sqresid=0,fit=0) #keep the basic results here
if (dim(x)[1]!=dim(x)[2]) x <- cor(x,use="pairwise") # if given a rectangular
for (i in 1:n) #loop through 1 to the number of factors requested
{
if(!(pc=="pc")) { if ( pc=="pa") {
f <- factor.pa(x,i,rotate=rotate,...) #do a factor analysis with i factors and the rotations specified in the VSS call
if (i==1)
{original <- x #just find this stuff once
sqoriginal <- original*original #squared correlations
totaloriginal <- sum(sqoriginal) - diagonal*sum(diag(sqoriginal) ) #sum of squared correlations - the diagonal
}} else {
f <- factanal(x,i,rotation=rotate,covmat=x,n.obs=n.obs,...) #do a factor analysis with i factors and the rotations specified in the VSS call
if (i==1)
{original <- x #just find this stuff once
sqoriginal <- original*original #squared correlations
totaloriginal <- sum(sqoriginal) - diagonal*sum(diag(sqoriginal) ) #sum of squared correlations - the diagonal
}}
} else {f <- principal(x,i)
if (i==1)
{original <- x #the input to pc is a correlation matrix, so we don't need to find it again
sqoriginal <- original*original #squared correlations
totaloriginal <- sum(sqoriginal) - diagonal*sum(diag(sqoriginal) ) #sum of squared correlations - the diagonal
}
if((rotate=="varimax") & (i>1)) {f <- varimax(f$loadings)}
}
load <- as.matrix(f$loadings ) #the loading matrix
model <- load residual <- original-model #find the residual R* = R - FF'
sqresid <- residual*residual #square the residuals
totalresid <- sum(sqresid)- diagonal * sum(diag(sqresid) ) #sum squared residuals - the main diagonal
fit <- 1-totalresid/totaloriginal #fit is 1-sumsquared residuals/sumsquared original (of off diagonal elements
if ((pc=="mle")) {
vss.df[i,1] <- f$dof #degrees of freedom from the factor analysis
vss.df[i,2] <- f$STATISTIC #chi square from the factor analysis
vss.df[i,3] <- f$PVAL #probability value of this complete solution
}
vss.df[i,4] <- totalresid #residual given complete model
vss.df[i,5] <- fit #fit of complete model
#now do complexities -- how many factors account for each item
for (c in 1:i)
{
simpleload <- complexmat(load,c) #find the simple structure version of the loadings for complexity c
model <- simpleload residual <- original- model #R* = R - SS'
sqresid <- residual*residual
totalsimple <- sum(sqresid) -diagonal * sum(diag(sqresid)) #default is to not count the diagonal
simplefit <- 1-totalsimple/totaloriginal
complexresid[i,c] <-totalsimple
complexfit[i,c] <- simplefit
}
} #end of i loop for number of factors
vss.stats <- data.frame(vss.df,cfit=complexfit,cresidual=complexresid)
return(vss.stats)
} #end of VSS functionRun the code above in your browser using DataLab