R <- matrix(.5, 8, 8) diag(R) <- 1 # Given that Theta = 0 out <- eigen(R) eigenvectors <- out$vectors eigenvalues <- out$values eigenvalues <- diag(eigenvalues) eigenvectors %*% eigenvalues %*% t(eigenvectors) # eigenvectors -> factor loadings # eigenvalues -> factor variances # Assuming that all factors are independent. R2 <- matrix(.2, 8, 8) R2[1:4, 1:4] <- .6 R2[5:8, 5:8] <- .6 diag(R2) <- 1 eigen(R2) dat <- MASS::mvrnorm(1000000, rep(0, 8), R, empirical = TRUE) dat <- data.frame(dat) colnames(dat) <- paste0("x", 1:8) out <- lm(x1 ~ x2 + x3 + x4 + x5 + x6 + x7 + x8, data = dat) rsq1 <- summary(out)$r.squared dat <- MASS::mvrnorm(1000000, rep(0, 8), R2, empirical = TRUE) dat <- data.frame(dat) colnames(dat) <- paste0("x", 1:8) out <- lm(x1 ~ x2 + x3 + x4 + x5 + x6 + x7 + x8, data = dat) rsq2 <- summary(out)$r.squared C <- R C2 <- R2 diag(C) <- rsq1 diag(C2) <- rsq2 eigen(C) eigen(C2) vec1 <- eigen(C)$vectors[,1,drop = FALSE] val1 <- eigen(C)$values[1] newc <- vec1 %*% val1 %*% t(vec1) diag(newc) <- 1 newc # Model implied correlation matrix vec2 <- eigen(C2)$vectors[,1:2,drop = FALSE] val2 <- eigen(C2)$values[1:2] newc2 <- vec2 %*% diag(val2) %*% t(vec2) diag(newc2) <- 1 newc2 # Model implied correlation matrix # Tabachnick & Fidell rotmat <- matrix(c(cos(pi/4), -sin(pi/4), sin(pi/4), cos(pi/4)), 2, 2, byrow = TRUE) round(vec2 %*% solve(rotmat), 3) library(psych) library(lavaan) library(semTools) library(GPArotation) # Factanal # WISC-III set.seed(123321) R <- diag(13) R[2,1] <- .66 R[3,1:2] <- c(.57, .55) R[4,1:3] <- c(.70, .69, .54) R[5,1:4] <- c(.56, .59, .47, .64) R[6,1:5] <- c(.34, .34, .43, .35, .29) R[7,1:6] <- c(.47, .45, .39, .45, .38, .25) R[8,1:7] <- c(.21, .20, .27, .26, .25, .23, .18) R[9,1:8] <- c(.40, .39, .35, .40, .35, .20, .37, .28) R[10,1:9] <- c(.48, .49, .52, .46, .40, .32, .52, .27, .41) R[11,1:10] <- c(.41, .42, .39, .41, .34, .26, .49, .24, .37, .61) R[12,1:11] <- c(.35, .35, .41, .35, .34, .28, .33, .53, .36, .45, .38) R[13,1:12] <- c(.18, .18, .22, .17, .17, .14, .24, .15, .23, .31, .29, .24) R <- R + t(R) diag(R) <- 1 varnames <- c("Information", "Similarities", "Arithmetic", "Vocabulary", "Comprehension", "Digit Span", "Pic. Completion", "Coding Subscale", "Pic. Arrang.", "Block Design", "Object Assembly", "Symbol Search", "Mazes") colnames(R) <- rownames(R) <- varnames out1 <- princomp(covmat = R) summary(out1) loadings(out1) out1$sdev sum(out1$sdev^2) screeplot(out1) # Using maximum likelihood for factor extraction out2 <- factanal(covmat = R, factors = 3, n.obs = 200, rotation = "varimax") loadings(out2) names(out2) out2$loadings out2$uniqueness out2$correlation out2$method out3 <- factanal(covmat = R, factors = 3, n.obs = 200, rotation = "none") lambda3 <- out3$loadings library(GPArotation) GPFoblq(lambda3, method = "quartimin") GPForth(lambda3, method = "bifactor") out4 <- factanal(covmat = R, factors = 2, n.obs = 200, rotation = "none") lambda4 <- out4$loadings T4 <- matrix(0, 13, 2) T4[1:5, 1] <- 1 T4[6:13, 2] <- 1 GPFoblq(lambda4, method = "target", methodArgs = list(Target = T4)) ?rotations T42 <- matrix(0, 13, 2) T42[1:5, 1] <- 1 T42[6:13, 2] <- 1 W42 <- matrix(0, 13, 2) W42[1:5, 1] <- 1 W42[6:13, 2] <- 1 GPFoblq(lambda4, method = "pst", methodArgs = list(Target = T42, W = W42)) dat <- round(MASS::mvrnorm(300, rep(100, 13), lavaan::cor2cov(R, rep(15, 13)), empirical = TRUE)) dat <- data.frame(dat) names(dat) <- paste0("v", 1:13) out5 <- factanal(~v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12+v13 ,data = dat, factors = 3, rotation = "varimax", scores = "Bartlett") out5$scores ################## psych package colnames(dat) <- varnames out1 <- fa(dat, nfactors = 1, fm = "ml") out2 <- fa(dat, nfactors = 2, fm = "ml") loadings(out2) fa.diagram(out1) fa.diagram(out2) out1$STATISTIC out1$dof pchisq(out1$STATISTIC, out1$dof, lower.tail = FALSE) out1$RMSEA # Preacher et al. (2011) # AIC = chi^2 + J(J+1) - 2*df out1$STATISTIC + 13 * 14 - 2 * out1$dof pchisq(out1$STATISTIC - out2$STATISTIC, out1$dof - out2$dof, lower.tail = FALSE) fa.parallel(dat) out3 <- fa(dat, nfactors = 3, rotate = "bifactor") out4 <- fa(dat, nfactors = 4, rotate = "bifactor") out3a <- fa(dat, nfactors = 3, n.obs = 200, rotate = "bifactor", scores = "Anderson") factor.scores(dat, out3a, method = "Anderson") X <- factor.scores(dat, out2, method = "tenBerge") newdat <- dat[1, , drop = FALSE] predict(out3a, newdat, dat) # Find structure matrix out2$loadings %*% out2$Phi KMO(R) cortest.bartlett(R, n = 200) fa.sort(out2) fa.organize(out2, o = c(2, 1), i = 13:1) factor2cluster(out2) # only appropriate for very simple structure # Rhemtulla et al. (2012) out5 <- fa.poly(iqitems, nfactors = 3, rotate = "quartimin") fa.parallel(out5)