Simulate data

library("MASS")
set.seed(42)
D <-3
P <- 10 
N <-1000

mu_theta <- rep(0,D) # the mean of eta
mu_epsilon <- rep(0,P) # the mean of epsilon
Phi <- diag(D)
Psi <- diag(P)
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3<-  c(0.00, 0.00, 0.85, 0.80, 0.00, 0.75, 0.75, 0.00, 0.80, 0.80)
L <-cbind(l1,l2,l3) # the loading matrix

Theta <-mvrnorm(N, mu_theta, Phi) # sample factor scores
Epsilon <-mvrnorm(N, mu_epsilon, Psi) # sample error vector
Y<-Theta%*%t(L)+Epsilon# generate observable data

Fit EFA

fa_model <- factanal(Y, D)
fa_model$loadings <- fa_model$loadings[, c(2, 3, 1)]
fa1 <- psych::fa(Y, nfactors = 3, rotate = "varimax")
cbind(fa_model$loadings, L, fa1$loadings[,c(2, 3, 1)])


jwmortensen/roster documentation built on May 30, 2019, 3:09 p.m.