Nothing
#vss is just an alias to VSS to be consistent with naming conventions
#added the RMSEA, BIC, SABIC and complexity criteria 1/27/14
"VSS" <-
function (x,n=8,rotate="varimax",diagonal=FALSE,fm="minres",n.obs=NULL,plot=TRUE,title="Very Simple Structure",use="pairwise",cor="cor",...) #apply the Very Simple Structure Criterion for up to n factors on data set x
{vss(x=x,n=n,rotate=rotate,diagonal=diagonal,fm=fm,n.obs=n.obs,plot=plot,title=title,use=use,cor=cor,...) }
"vss" <-
function (x,n=8,rotate="varimax",diagonal=FALSE,fm="minres",n.obs=NULL,plot=TRUE,title="Very Simple Structure",use="pairwise",cor="cor",...) #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 "varimax"
#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)
{
cl <- match.call()
if (rotate=="oblimin") {if(!requireNamespace('GPArotation')) {stop("You must have GPArotation installed to use oblimin rotation")}}
old_rotate=rotate #used to remember which rotation to use
#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)
}
map <- function(x,n) {
nvar <- dim(x)[2]
min.partial <- rep(NA,n)
e <- eigen(x)
evect <- e$vectors
comp <- evect %*% diag(sqrt(e$values))
if( n >= nvar) {n1 <- nvar -1} else {n1 <- n}
for (i in 1:n1) {
c11.star <- x - comp[,1:i] %*% t(comp[,1:i])
d <- diag(1/sqrt(diag(c11.star)))
rstar <- d %*% c11.star %*% d
diag(rstar) <- 0
min.partial[i] <- sum(rstar * rstar) /(nvar*(nvar-1))
}
return(min.partial)
}
if(dim(x)[2] < n) n <- dim(x)[2]
#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=NA,prob=NA,sqresid=NA,fit=NA,RMSEA=NA,BIC=NA,SABIC=NA,complex=NA,eChisq=NA,SRMR=NA,eCRMS=NA,eBIC=NA) #keep the basic results here
if (dim(x)[1]!=dim(x)[2]) { n.obs <- dim(x)[1]
switch(cor,
cor = {x <- cor(x,use=use)},
cov = {x <- cov(x,use=use)
covar <- TRUE},
tet = {x <- tetrachoric(x)$rho},
tetrachoric = {x <- tetrachoric(x)$rho},
poly = {x <- polychoric(x)$rho},
polychoric = {x <- polychoric(x)$rho},
mixed = {x <- mixedCor(x,use=use)$rho},
Yuleb = {x <- YuleCor(x,,bonett=TRUE)$rho},
YuleQ = {x <- YuleCor(x,1)$rho},
YuleY = {x <- YuleCor(x,.5)$rho }
)
# x <- cor(x,use="pairwise") The case statement allows many types of correlations
} else {if(!is.matrix(x)) x <- as.matrix(x)}
# if given a rectangular
if(is.null(n.obs)) {message("n.obs was not specified and was arbitrarily set to 1000. This only affects the chi square values.")
n.obs <- 1000}
map.values <- map(x,n)
if (n > dim(x)[2]) {n <- dim(x)[2]} #in cases where there are very few variables
for (i in 1:n) #loop through 1 to the number of factors requested
{ PHI <- diag(i)
if(i<2) {(rotate="none")} else {rotate=old_rotate}
if(!(fm=="pc")) {
f <- fa(x,i,rotate=rotate,n.obs=n.obs,warnings=FALSE,fm=fm,scores="none",cor=cor,...) #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)
PHI <- diag(i)} else {
if(((rotate=="promax")| (rotate=="Promax") )& (i>1)) {f <- Promax(f$loadings)
PHI <- f$Phi} else {
if((rotate=="oblimin")& (i>1)) {f <- GPArotation::oblimin(f$loadings)
U <- f$Th
phi <- t(U) %*% U
PHI <- cov2cor(phi)
}}
}}
load <- as.matrix(f$loadings ) #the loading matrix
model <- load %*% PHI %*% t(load) #reproduce the correlation matrix by the factor law R= FF'
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 ((fm !="pc")) {
vss.df[i,"dof"] <- f$dof #degrees of freedom from the factor analysis
vss.df[i,"chisq"] <- f$STATISTIC #chi square from the factor analysis
vss.df[i,"prob"] <- f$PVAL #probability value of this complete solution\
vss.df[i,"eChisq"] <- f$chi
vss.df[i,"SRMR"] <- f$rms
vss.df[i,"eRMS"] <- f$rms
vss.df[i,"eCRMS"] <- f$crms
vss.df[i,"eBIC"] <- f$EBIC
if(!is.null(f$RMSEA)) {vss.df[i,"RMSEA"] <- f$RMSEA[1]} else {vss.df[i,"RMSEA"] <- NA}
if(!is.null(f$BIC)) {vss.df[i,"BIC"] <- f$BIC} else {vss.df[i,"BIC"] <- NA}
if(!is.null(f$SABIC)) {vss.df[i,"SABIC"] <- f$SABIC} else {vss.df[i,"SABIC"] <- NA}
if(!is.null(f$complexity)) {vss.df[i,"complex"] <- mean(f$complexity)} else {vss.df[i,"complex"] <- NA}
}
vss.df[i,"sqresid"] <- totalresid #residual given complete model
vss.df[i,"fit"] <- 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 %*% PHI %*% t(simpleload) #the model is now a simple structure version R ? SS'
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)
if (plot) VSS.plot(vss.stats,title=title)
vss.results <- list(title=title,map=map.values,cfit.1=complexfit[,1],cfit.2= complexfit[,2],vss.stats=vss.stats,call=cl)
class(vss.results) <- c("psych" ,"vss")
return(vss.results)
} #end of VSS function
"nfactors" <-
function(x,n=20,rotate="varimax",diagonal=FALSE,fm="minres",n.obs=NULL,title="Number of Factors",pch=16,use="pairwise",cor="cor",...) {
vs <- vss(x=x,n=n,rotate=rotate,diagonal=diagonal,fm=fm,n.obs=n.obs,plot=FALSE,title=title,use=use,cor=cor,...)
old.par <- par(no.readonly = TRUE) # save default, for resetting...
on.exit(par(old.par)) #and when we quit the function, restore to original values
op <- par(mfrow=c(2,2))
x <- vs$vss.stats
n = dim(x)
plot(x$cfit.1, ylim = c(0, 1), typ = "b", ylab = "Very Simple Structure Fit",
xlab = "Number of Factors",main="Very Simple Structure",pch=49)
lines(x$cfit.1)
x$cfit.2[1] <- NA
x$cfit.3[1] <- NA
x$cfit.3[2] <- NA
lines(x$cfit.2)
points(x$cfit.2,pch=50)
lines(x$cfit.3)
points(x$cfit.3,pch=51)
plot(vs$vss.stats[,"complex"],xlab="Number of factors",ylab="Complexity",typ="b",main="Complexity",pch=pch,...)
plot(vs$vss.stats[,"eBIC"],xlab="Number of factors",ylab="Empirical BIC",typ="b",main="Empirical BIC",pch=pch,...)
plot(vs$vss.stats[,"SRMR"],xlab="Number of factors",ylab="SRMR",typ="b",main="Root Mean Residual",pch=pch,...)
results <- list(title=title,map=vs$map,vss.stats=vs$vss.stats[,1:16],call=vs$call)
class(results) <- c("psych","vss")
return(results)
}
#added 2/17/23
#implementing an idea from an article by Larsen and Warne 2010
"eigenCi" <- function(x,n.iter=1000, use="pairwise", alpha=.05,plot=FALSE,root=FALSE) {
evo <- eigen(cor(x,use=use))$values
res <- list()
n.items <- NROW(x)
nvar <- NCOL(x)
for(i in 1:n.iter) {
ss <- sample(n.items,n.items,replace=TRUE)
R <- cor(x[ss,], use=use)
ev <- eigen(R)$values
res[[i]] <- c(ev)
}
res.mat <- matrix(unlist(res),ncol=nvar,byrow=TRUE)
#sd <- qnorm(alpha/2) * apply(res.mat,2, sd)
boot <- apply(res.mat,2,function(x) quantile(x,c(alpha/2,(1-alpha/2))))
norm <- qnorm(alpha/2)* sqrt((2 * evo^2)/n.items)
ci <- data.frame(ev=evo,nlow = evo+norm,nhigh=evo-norm, blow=boot[1,],bhigh =boot[2,])
if(plot) { if(root) {ci <- sqrt(ci)
matplot(ci,type="l",ylab="sqrt eigen values",main="Confidence intervals of eigen values")} else {
matplot(ci,type="l",ylab="eigen values",main="Confidence intervals of eigen values")}}
return(ci)
}
vssSelect <- function(keys,x,cor="cor", fm="minres",plot=FALSE) {
temp <- list()
for (i in 1:length(keys)) {
select <- selectFromKeys(keys[[i]])
VS <- vss(x[,select],cor=cor,fm=fm,plot=plot)
temp[[i]] <- list(first.max(VS$cfit.1),first.max(VS$cfit.2))
}
result <- as.data.frame(matrix(unlist(temp),ncol=2,byrow=TRUE))
colnames(result) <- c("VSS.1", "VSS.2")
rownames(result) <- names(keys)
return(result)
}
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.