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.



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) 
print(efa2$loadings,cutoff = 0.5)

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?


1 Answer 1


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

level: 1
level: 2

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]

> head(df)
   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)
  • $\begingroup$ 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? $\endgroup$
    – Dave
    Jun 22, 2020 at 21:35
  • $\begingroup$ 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. $\endgroup$ Jun 30, 2020 at 0:55
  • $\begingroup$ 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! $\endgroup$
    – Dave
    Jul 1, 2020 at 16:04
  • $\begingroup$ 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. $\endgroup$ Jul 3, 2020 at 5:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.