# how can I perform confirmatory factor analysis when covariance matrix contains missing values?

I have dataset of test items which were administered in blocks. As a result, not all students answered each test item and there are some pairs of items for which no observations are shared, introducing missing values into the covariance matrix. As an example, suppose I have the following data. I'll provide an example in Python but your answer doesn't need to depend on a programming language.

import numpy as np
import pandas as pd
random.seed(23)

cov = np.array([[2,0.5],[0.5,1]])

s = np.random.multivariate_normal(mean=[0, 0],
cov=[[2,0.5],[0.5,1]], size=N)

m1 = 0.5 + 1*s[:, 0] + np.random.normal(scale=2, size=N)
m2 = 2 + 1*s[:, 1] + np.random.normal(scale=1.5, size=N)
m3 = -3 + 2*s[:, 0] -4*s[:, 1] + np.random.normal(scale=4, size=N)
m4 = -1 + 1.5*s[:, 0] + 2*s[:, 1] + np.random.normal(scale=0.5, size=N)
m5 = 3 - 1*s[:, 0] + 1*s[:, 1] + np.random.normal(scale=3, size=N)
m6 = 2 + 1.2*s[:, 0] + 5*s[:, 1] + np.random.normal(scale=1, size=N)

M = np.vstack([m1, m2, m3, m4, m5, m6]).T

M[:int(N/2), 4] = np.nan
M[int(N/2):, 5] = np.nan

print(pd.DataFrame(M).corr().values)


The resulting covariance matrix is

array([[ 1.  ,  0.12,  0.16,  0.47, -0.19,  0.32],
[ 0.12,  1.  , -0.3 ,  0.45,  0.06,  0.53],
[ 0.16, -0.3 ,  1.  , -0.16, -0.27, -0.39],
[ 0.47,  0.45, -0.16,  1.  , -0.12,  0.91],
[-0.19,  0.06, -0.27, -0.12,  1.  ,   nan],
[ 0.32,  0.53, -0.39,  0.91,   nan,  1.  ]])


which shows that I have missing values for the covariance on the last two measures. In my exact problem I have additional missing components.

I'd like to estimate a confirmatory factor analysis model using all of the available data, but I'm not sure how to do it with these missing values. I will be fixing a priori zeros (see Rubin and Thayer, 1982) in the loading matrix to secure identification, so I need the flexibility to be able to do this.

I'd be open to learning if this is possible in lavaan or factanal in R or using packages in Python. I'd also like to know the references for the methods used to address this problem. I'm thinking there must be some FIML approach. But most articles I've found address only the case where the data contain missing values but there at least a few observations that have nonmissing values for each pair of columns.

• FIML requires information from the individual observations, so you can't do it with just a covariance matrix. – Noah Apr 11 '19 at 22:21
• Do you have any references for a FIML approach in this context? And do you know if there is an existing package in R or Python that can handle this missing pattern? – jtorca Apr 11 '19 at 22:32
• I am less familiar with Python syntax. Can you explain why there would be missing values in your covariance matrix. Are you missing the raw values? Also you'll want non-standardized (i.e., not the correlation matrices) to feed into most SEM software. Otherwise you'll need to add back in the variance and (and sometimes mean vectors) depending on your model. If you have the raw vectors of your variables then FIML is fine and lavaan can handle missingness at random with the missing='fiml' argument added to the sem() function. – Matt Barstead Apr 12 '19 at 16:08
• I have a dataset where subjects were randomly administered different sets of questions and it just so happens there are pairs of questions for which no group took both questions. Thus, I can't compute a covariance matrix even pairwise. However, there should be sufficient information in there to recover all of the loadings. – jtorca Apr 12 '19 at 17:01

So I at least have the code, given the suggestion from @Matt Barstead. I used the lavaan package in R.

library(lavaan)
library(MASS)

set.seed(23)

N = 1000

cov = cbind(c(2,0.5),c(0.5,1))

s = mvrnorm(n=N, mu=c(0, 0), Sigma=cov)

m1 = 0.5 + 1*s[, 1] + rnorm(sd=2, n=N)
m2 = 2 + 1*s[, 2] + rnorm(sd=1.5, n=N)
m3 = -3 + 2*s[, 1] -4*s[, 2] + rnorm(sd=4, n=N)
m4 = -1 + 1.5*s[, 1] + 2*s[, 2] + rnorm(sd=0.5, n=N)
m5 = 3 - 1*s[, 1] + 1*s[, 2] + rnorm(sd=3, n=N)
m6 = 2 + 1.2*s[, 1] + 5*s[, 2] + rnorm(sd=1, n=N)

M = data.frame(cbind(m1, m2, m3, m4, m5, m6))

# creating missing data pattern that would create missing value in
# covariance matrix
M[1:N/2, 1] = NaN
M[(N/2):N, 2] = NaN

cfa.model <- '
f1 =~ 1*m1 + m3 + m4 + m5 + m6
f2 =~ 1*m2 + m3 + m4 + m5 + m6
f1 ~ f2
'

fcols = c('f1', 'f2')

K = length(fcols)
J = dim(M)[2]

cfa.fit = cfa(cfa.model, data=M, missing='fiml')

A = inspect(cfa.fit, what='est')\$lambda[c('m1', 'm2', 'm3', 'm4', 'm5', 'm6'), c('f1', 'f2')]
I = inspect(cfa.fit, what='cov.lv')


Result for A and I:

> A
f1         f2
m1  1.0000000  0.0000000
m2  0.0000000  1.0000000
m3  1.7334022 -3.9040606
m4  1.3276143  2.0064776
m5 -0.9336488  0.7713415
m6  0.9638958  5.1738190
> I
f1    f2
f1 2.224
f2 0.619 1.042


So I can recover the loadings I used to generate the data.

On a side note, I tried to use the sem package as well, and it worked fine when not I don't specify na.action=na.pass, but when I try to use FIML with na.action=na.pass, I get a warning message and I don't recover the parameters. Here is an incomplete snippet

cfa.fit = sem(cfa.model, data=M, na.action=na.pass)


But then I get

Warning message:
In sem.semmod(cfa.model, data = M, optimizer = optimizerOptim, na.action = na.pass) :
The following observed variables are in the input covariance or raw-moment matrix but do not appear in the model:
Intercept

• In what sense could you not get FIML to work? It works fine when I run that code. – Noah Apr 12 '19 at 20:13
• I might not know the proper argument when I try sem instead of lavaan. When using sem, I pass na.action=na.pass to try to use FIML, but then I get the error  The following observed variables are in the input covariance or raw-moment matrix but do not appear in the model: Intercept – jtorca Apr 12 '19 at 20:28