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

Question

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.

Answer 1

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

# Note the 1* fixes the loading to 1

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')

# grabbing factor loadings (A) and factor covariance matrix (I)
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