R: Using MLE to do a CFA by hand

Question

EDIT: To fix mistake in covariance matrix

Cross posted from here

Ok so I am trying to do a cfa by hand. I am doing the tutorial found here

The central equation is the model implied covariance matrix:

Sigma = LL^T+E

Where L is a matrix of loadings and E is a matrix of errors. The idea is to use maximum liklihood to get estimates for L and E which minimise the descrepency between the implied and the observed covariance matrix. The descrepency function to be minimized is:

F = log |Sigma| + tr(S*Sigma^-1) - log|S| - p

where sigma is the implied covariance matrix and S is the observed covariance matrix and p is the number of indicator items.

Ok so the observed covariance matrix is cut and pasted from the link above but put into an R cov matrix by:

dat<- '.804 .399 .500 .367 .451 .510 .487 .833 .433 .283 .372 .377 .621 .529 .805 .339 >.551 .543 .478 .362 .441 .733 .332 .341 .570 .461 .695 .438 .780 .556 .626 .455 .667 .438 >.693 .825'

dat<- strsplit(dat, ' ')

dat<-as.numeric(unlist(dat))

covar1<- matrix(dat, nrow=6, ncol=6, byrow=FALSE)

covar1<-t(covar1)

covar1<-as.matrix(forceSymmetric(covar1))

Starting values and function for ML is:

L1<- rep(.7,6)#starting values for loadings

E1<- diag(.3, 6, 6)#starting values for error

discrepency<- function(covar, L, E){
  #descrepency function 
  # I need to apply ML to this to find L and E
  sigma<- L%*%t(L)+E
  log(det(sigma))+
  sum(diag(covar%*%solve(sigma)))-
  log(det(covar))-
  nrow(covar)}

My question is how do I get the ML estimates for L and E. In other words how do I get estimates for L (column vector) and E (diagonal matrix) which minimize the descrepency value. I have tried stats4:::mle from base but with no success. Obviously I would like to use my starting values for L and E as noted above.

Many thanks!

EDIT: Here is where I have got to so far. I have tried the mle function as follows:

stats4:::mle(discrepency, fixed= list(covar=covar1), start = list(L=L1, E=E1))

The error I get is:

Error in optim(start, f, method = method, hessian = TRUE, ...) : (list) object cannot be coerced to type 'double'

I get the mle estimates to work by rewritting the descrepency function to take individual entries for each of the factor loading and error start values. Nevertheless while the estimates are in the ball park of the Mplus estimates they are still a little concerningly far away. — Philip Parker

Philip Parker Philip Parker · Accepted Answer · 2012-03-14T20:05:02

Success:

discrepency<- function(covar, L1,L2,L3,L4,L5,L6, E1,E2,E3,E4,E5,E6){
  #descrepency function 
  # I need to apply ML to this to find L and E
  sigma<- c(L1,L2,L3,L4,L5,L6)%*%t(c(L1,L2,L3,L4,L5,L6))+diag(c(E1,E2,E3,E4,E5,E6))
  log(det(sigma))+
  sum(diag(covar%*%solve(sigma)))-
  log(det(covar))-
  nrow(covar)}

With:

out<-stats4:::mle(discrepency, start = list(L1=.07, L2=.07, L3=.07, L4=.07, L5=.07,L6=.07,
      E1=.03,E2=.03,E3=.03,E4=.03,E5=.03,E6=.03), 
      fixed= list(covar=covar1), method="L-BFGS-B", 
      control=list(maxit=10000))

str(out)
out1<-round(matrix(out@coef, 6, 2),3)
colnames(out1)<- c('loadings', 'errors')

Gives MLE estimates that are very similar to MPlus and AMOS

R: Using MLE to do a CFA by hand

1 Answers