#This is an R-script performing RFT (Retained-Components factor transformation), see Beauducel, A. & Spohn, F. (2014). Retained-Components Factor Transformation: Factor Loadings and Factor Score Predictors in the Column Space of Retained Components. Journal of Modern Applied Statistical Methods, 13(2). #Authors of the R-script: André Beauducel (University of Bonn) and Frank Spohn (University of Hamburg). # Please specify the directory containing your data to be analyzed path <- "C:/[...]" ## As one possibility: Import SPSS file: library(foreign) x <- as.matrix(read.spss(paste(path,"[SPSS_filename].sav",sep = ""), use.value.labels=FALSE, max.value.labels=Inf, to.data.frame=TRUE)) ## import as data set ## As an alternative: Read text file: # Please customize the file name containing your data matrix (ASCII-file with "." as decimal separator and tabulator as field separator). # x <- read.table(paste(path,"[text_filename].txt",sep = ""), header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE) eVec <- princomp(x, cor = TRUE)$loadings # compute eigenvectors with "princomp" (its "loadings" are essentially eigenvectors) eVec2 <- eigen(cor(x))$vectors # alternatively compute eigenvectors without "princomp" print(eVec, cutoff = 0) print(eVec2) # just to satisfy oneself that both methods render the same result (eVec2 is not used in further calculations) eVal <- princomp(x, cor = TRUE)$sdev**2 # compute eigenvalues with "princomp" eVal2 <- eigen(cor(x))$values # compute eigenvalues without "princomp" print(eVal) print(eVal2) # just to satisfy oneself that both methods render the same result (eVal2 is not used in further calculations) diag_eVal <- diag(eVal) # compute diagonal matrix with eigenvalues as elements A <- eVec%*%sqrt(diag_eVal) # compute loading matrix of all components of PCA print(A) # print PCA loadings to screen ## can be omitted q <- 3 # set number of retained components (please insert your own value here) M <- A[,1:q] # define M (loading matrix of retained components) N <- A[,(q+1):ncol(x)] # define N (loading matrix of not retained components) install.packages("psych") library(psych) ULFA <- fa(r=x, nfactors=q, rotate="none", fm="minres") # compute ULFA/Minres factor analysis (ULFA) MLFA <- factanal(x, factors=q, rotation="none") # compute maximum likelihood factor analysis (MLFA) FA <- ULFA # please set your preferred FA algorithm here (MLFA or ULFA) summary(FA) ## can be omitted Lambda <- loadings(FA) # extract loadings from FA print(Lambda,cutoff=0) # print unrotated FA loadings to screen ## can be omitted RFT <- Lambda-N%*%solve(t(N)%*%N)%*%t(N)%*%Lambda # compute unrotated RFT RFT <- round(RFT,digits=2) # specify the number of decimal places (you want to use) for output print(RFT,cutoff=0) # print unrotated RFT loadings to screen # Rotation -------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Subsequently different rotation methods can be chosen by deleting or inserting the comment character (#) at the beginning of each line of the respective paragraph. # Varimax: # Lambda_rot <- varimax(Lambda, normalize = TRUE, eps = 1e-5)$loadings # compute Varimax rotated FA loadings # RFT_rot <- varimax(RFT, normalize = TRUE, eps = 1e-5)$loadings # compute Varimax rotated RFT loadings # Promax: #promax_power <- 4 # define the power value for the Promax rotation #Lambda_rot <- promax(Lambda, m=promax_power)$loadings # compute Promax rotated FA loadings #RFT_rot <- promax(RFT, m=promax_power)$loadings # compute Promax rotated RFT loadings # Oblimin: install.packages("GPArotation") library(GPArotation) oblimin_kind <- 0 # define kind of Oblimin rotation(0=Quartimin, .5=Biquartimin, 1=Covarimin) # rotation of FA: obli_FA <- oblimin(Lambda, Tmat=diag(ncol(Lambda)), gam=oblimin_kind, normalize=TRUE, eps=1e-5, maxit=1000) Lambda_rot <- obli_FA$loadings # compute Oblimin rotated FA loadings (pattern matrix) # Phi from GPArotation package: Phi_Lambda_GPA <- obli_FA$Phi print(Phi_Lambda_GPA,digits=2,cutoff=0,sort=FALSE) # Phi "manually" computed: Phi_Lambda <- solve(t(Lambda_rot)%*%Lambda_rot)%*%t(Lambda_rot)%*%(Lambda%*%t(Lambda))%*%Lambda_rot%*%solve(t(Lambda_rot)%*%Lambda_rot) # rotation of RFT: obli_RFT <- oblimin(RFT, Tmat=diag(ncol(RFT)), gam=oblimin_kind, normalize=TRUE, eps=1e-5, maxit=1000); RFT_rot <- obli_RFT$loadings # compute Oblimin rotated RFT loadings (pattern matrix) # Phi from GPArotation package: Phi_RFT_GPA <- obli_RFT$Phi print(Phi_RFT_GPA,digits=2,cutoff=0,sort=FALSE) # Phi "manually" computed: Phi_RFT <- solve(t(RFT_rot)%*%RFT_rot)%*%t(RFT_rot)%*%(RFT%*%t(RFT))%*%RFT_rot%*%solve(t(RFT_rot)%*%RFT_rot) # compute structure matrix (identical with pattern matrix of orthogonal rotation): Lambda_strukt <- Lambda_rot%*%Phi_Lambda RFT_strukt <- RFT_rot%*%Phi_RFT print(Lambda_rot,digits=2,cutoff=0,sort=FALSE) # print rotated FA loadings to screen ## can also be omitted print(Phi_Lambda,digits=2,cutoff=0,sort=FALSE) print(Lambda_strukt,digits=2,cutoff=0,sort=FALSE) print(RFT_rot,digits=2,cutoff=0,sort=FALSE) # print rotated RFT to screen print(Phi_RFT,digits=2,cutoff=0,sort=FALSE) print(RFT_strukt,digits=2,cutoff=0,sort=FALSE) # -------------------------------------------------------------------------------------------------------------------------------------------------------------------- SS_Lambda <- rep(0,q) # create new vectors for sum of squares (SS) of loadings SS_Lambda_rot <- rep(0,q) SS_RFT <- rep(0,q) SS_RFT_rot <- rep(0,q) for (i in 1:q){SS_Lambda[i] <- sum(Lambda[,i]^2)} # compute SS of loadings for (i in 1:q){SS_Lambda_rot[i] <- sum(Lambda_rot[,i]^2)} for (i in 1:q){SS_RFT[i] <- sum(RFT[,i]^2)} # compute SS of loadings for (i in 1:q){SS_RFT_rot[i] <- sum(RFT_rot[,i]^2)} diff <- SS_RFT-SS_Lambda # compute differences between loadings from MLFA/ULFA and RFT names(diff) <- labels(Lambda)[[2]] print(diff) # print differences diff_rot <- SS_RFT_rot-SS_Lambda_rot # compute differences between loadings from rotated MLFA/ULFA and rotated RFT names(diff_rot) <- labels(Lambda)[[2]] print(diff_rot) # print differences write(t(Lambda), paste(path,"Lambda.txt"), sep="\t") ## save unrotated MLFA/ULFA to ASCII file with "." as decimal separator and tabulator as field separator write(t(Lambda_rot), paste(path,"Lambda_rot.txt"), sep="\t") ## save rotated MLFA/ULFA to file write(t(RFT), paste(path,"RFT.txt"), sep="\t") ## save RFT to file write(t(RFT_rot), paste(path,"RFT_rot.txt"), sep="\t") ## save rotated RFT to file #remove.packages("GPArotation") ## remove the package "GPArotation" used for Olimin rotation if you like #remove.packages("psych") ## remove the package "psych" used for ULFA if you like