Is it possible to run an EFA in R that also accounts for non-independent observations (i.e., includes random effects)?

Suppose I have the following data, in which 20 participants each rated 10 items on five different dimensions. In reality I have much more of each, but I am just using this as an example.

set.seed(123)

library(data.table)
library(psych)

ppt <- rep(1:20, each = 10)
item <- rep(1:10, times = 20)
dim_a <- rnorm(200)
dim_b <- rnorm(200)
dim_c <- rnorm(200)
dim_d <- rnorm(200)
dim_e <- rnorm(200)

d <- as.data.table(cbind(ppt, item, dim_a, dim_b, dim_c, dim_d, dim_e))

I am interested in conducting a factor analysis on the five dimensions. I could do so with something like this. (In this example it suggests no factors be extracted because of the random data, but I proceed as though it suggested two.)

parallel <- fa.parallel(d[, c(3:7)], fm = 'pa', fa = 'fa')
efa2 <- fa(d[, c(3:7)], nfactors =2, rotate = "oblimin", fm = 'pa', max.iter = 100000)

However, this doesn't account for the fact that ratings were nested in participants, and in items. It seems like what I need is to add crossed random effects for subjects and items. After a lot of searching, it seems like this isn't possible in R. Am I right?

It may come down to your personal definition of what counts as an exploratory factor analysis (EFA) versus confirmatory factor analysis (CFA). I am of the opinion that your research question and rationale is more important than the label assigned a particular analytical technique.

All of that is to say that you can use what many view as a confirmatory technique to include random effects in a factor analysis. You will just have to use more judgment along the way. The package for the task I’d recommend is lavaan. You can view the specific walkthrough for a factor analysis with random effects here.

Here is the starting point of a worked out example using a data set that includes mood ratings captured via cell phone surveys sent to the same sample of participants throughout the day for a week.

Using SEM requires a bit more intervention by the analyst, but if you have even the slightest inkling about the factor structure it can be a powerful tool.

UPDATE: Taking the model from the second link

CFA.mod<-'
level: 1
W_PA=~Joy+Cheer+Enthus+Content+Relax+Calm
W_ANG=~Angry+Annoy
W_TRD=~Tired+Slug
level: 2
B_PA=~Joy+Cheer+Enthus+Content+Relax+Calm
B_ANG=~Angry+Annoy
B_TRD=~Tired+Slug
'

I have two levels - within-subjects and between-subjects that correspond to level: 1 and level: 2. I can add error covariances at either level of the model, and they do not have to remain entirely parallel as I make changes. lavInspect(my_fitted_model, "icc"), I can get estimates of the intra-class correlations for my observed variables. I can also use utilities such as semTools::reliability(my_fitted_model) to get different measures of each factor's reliability at both the within- and between-subjects level of the data, which is probably of more interest to you.

I do not think it is possible to specify more complex error covariance structures than those that involve estimating random intercepts by "cluster" (the lme4 analog of (1|ID), using lavaan, at least not at this time.

UPDATE 2: So lavaan is going to be very grumpy if there is not much between-subjects variance in your data. I don't see how the explicit example you provided will actually simulate random intercepts correctly.

But let's say you had a data with a single row per subject per observation and a single column per cluster variable + all other indicator variables.Something like the following:

N <- 200 # number of subjects
N_OBS <- 20 # Number of observations per participant
N_DIM <- 5 # Number of dimensions
N_VARS <- 10 # Number of variables per dimension

[additional simulation setup would go here]

id       a1       a2       a3       a4       a5       a6       a7       a8       a9      a10        b1        b2        b3        b4        b5        b6
1 200 1.552593 1.906698 1.819515 1.538699 1.590809 1.551324 1.634444 1.560588 1.710813 1.666928 0.7730447 0.8586777 0.4257731 0.8126416 0.6923401 0.8931680
2 200 1.594748 1.844632 1.815612 1.705720 1.854236 1.679682 1.689307 1.720757 1.601351 1.881109 0.8542609 0.8178276 0.9521182 0.9554177 0.7219753 0.7922130
3 200 1.707210 1.759805 1.663649 1.689524 1.500286 1.896003 1.693405 1.913407 1.640697 1.546836 0.9058836 0.8923759 0.8735811 0.7730722 0.7188808 1.1076694
4 200 1.673065 1.711226 1.594529 1.752306 1.633699 1.966551 1.445854 1.790901 1.639650 1.678080 0.9401823 0.5947513 0.7591943 0.8665099 0.8602549 0.8071486
5 200 1.848208 1.816846 1.709582 1.773056 1.854829 1.754225 1.801330 1.443547 1.850241 1.750587 0.8541575 0.5547727 0.7788956 0.7172802 0.8643167 0.4549975
6 200 1.646956 2.009243 1.733745 1.948446 1.630401 1.526548 1.626063 1.406556 1.703369 1.718094 0.5662642 0.8351231 0.5792455 0.7647101 0.9929229 0.9768771

where the variablesa1-a10 are indicators for your latent factor A, the variables b1-b10 are indicators for your latent factor B, and so on up to E.

You could then build a CFA model (assuming you have simulated the data correctly) with random intercepts as follows (note I am making use of the glue package to make the formula creation a little more compact here.

model_structure <- glue("
level: 1
W_A =~ {paste(paste0('a', 1:N_VARS), collapse='+')}
W_B =~ {paste(paste0('b', 1:N_VARS), collapse='+')}
W_C =~ {paste(paste0('c', 1:N_VARS), collapse='+')}
W_D =~ {paste(paste0('d', 1:N_VARS), collapse='+')}
W_E =~ {paste(paste0('e', 1:N_VARS), collapse='+')}
level: 2
B_A =~ {paste(paste0('a', 1:N_VARS), collapse='+')}
B_B =~ {paste(paste0('b', 1:N_VARS), collapse='+')}
B_C =~ {paste(paste0('c', 1:N_VARS), collapse='+')}
B_D =~ {paste(paste0('d', 1:N_VARS), collapse='+')}
B_E =~ {paste(paste0('e', 1:N_VARS), collapse='+')}
")

Then you would go ahead and fit the model:

fit <- cfa(model_structure, data=df, cluster="id")

and per lavaan's documentation - if you have convergence issues, you may have to adjust your call to run more iterations.

fit <- cfa(model_structure, data=df, cluster="id", optim.method = "em", em.iter.max = 20000, em.fx.tol = 1e-08, em.dx.tol = 1e-04)
• Thanks very much. I'm having a bit of trouble understanding how the random effects are specified in that example. I'm used to lme4. Can lavaan handle crossed random item and subject effects?
– Dave
Jun 22 '20 at 21:35
• So I think the short answer is that lavaan is only really capable right now of giving you (1|ID) sorts of random effects. So you can get intercepts for every variable decomposed into "between"- and "within"-subjects sources of variability, but you can't specify a more complete random effects covariance structure. I've updated my answer with some additional info using my own worked out example in the second link. Jun 30 '20 at 0:55
• Thank you so much! I think that specifying random intercepts would be good enough for me! I'm still a bit unclear how this is done though... Would you mind to write out how it might be done for my example data? Thank you very much!
– Dave
Jul 1 '20 at 16:04
• Added a little more, your specific example won't simulate random intercepts, and lavaan will not return a meaningful solution if I tried to apply a multilevel model to data created under those assumptions. Instead of gaming out a valid simulation data set, I went with your basic concept and wrote out some code for how you might approach such a problem "in the wild" - hope it is helpful enough to get you going. Jul 3 '20 at 5:02