CODEWITHSUNDEEP
r
Philip Parker
Philip Parker

Reputation: 639

R: Using MLE to do a CFA by hand

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'

Upvotes: 2

Views: 575

Answers (1)

Philip Parker
Philip Parker

Reputation: 639

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

Upvotes: 0

Related Questions