Two models

We consider a one-factor model with 4 indicators and the Holzinger 3-factor model. For the one-factor simple model we generate a random sample:

library(lavaan)
## This is lavaan 0.6-7
## lavaan is BETA software! Please report any bugs.
#Simple factor model with random data
Sim.model <- "F=~x1+x2+x3+x4"
Sim.nobs <- 100
set.seed(1)
Sim.cov <- cov(simulateData(Sim.model, sample.nobs=Sim.nobs))

# Holzinger
HS.model <- " visual  =~ x1 + x2 + x3;textual =~ x4 + x5 + x6; speed   =~ x7 + x8 + x9 "
HS.cov <- cov(HolzingerSwineford1939[, paste0("x", 1:9)])
HS.nobs <- nrow(HolzingerSwineford1939)

Checking that the minimum in the GLS is the MLE

We put \(\Sigma(\hat \theta)^{-1}\) as \(V\) in GLS and minimize. It is confirmed that we regain the MLE.

# Simple model
fml <- sem(Sim.model, sample.cov=Sim.cov, sample.nobs=Sim.nobs)
sigmaimplied <- lavInspect(fml, "sigma.hat")
sigmaimpl_inv <- chol2inv(chol(sigmaimplied))
WLS.V <- lav_matrix_duplication_pre_post(kronecker(sigmaimpl_inv, sigmaimpl_inv))

#NB: need the sample covariance version with divisor n and not (n-1)
fgls <- sem(Sim.model,sample.cov=(Sim.nobs-1)*Sim.cov/Sim.nobs, sample.nobs=Sim.nobs, 
            estimator ="WLS", WLS.V = WLS.V)

mlest <- unname(coef(fml))
glsest<- unname(coef(fgls))
summary((mlest-glsest)/glsest)#identical..
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2.562e-07 -9.620e-08 -4.593e-08 -3.227e-08  5.941e-08  1.705e-07
# HS mdoel
fml <- sem(HS.model, sample.cov=HS.cov, sample.nobs=HS.nobs)
sigmaimplied <- lavInspect(fml, "sigma.hat")
sigmaimpl_inv <- chol2inv(chol(sigmaimplied))
WLS.V <- lav_matrix_duplication_pre_post(kronecker(sigmaimpl_inv, sigmaimpl_inv))

#NB: need the sample covariance version with divisor n and not (n-1)
fgls <- sem(HS.model,sample.cov=(HS.nobs-1)*HS.cov/HS.nobs, sample.nobs=HS.nobs, 
            estimator ="WLS", WLS.V = WLS.V)

mlest <- unname(coef(fml))
glsest<- unname(coef(fgls))
summary((mlest-glsest)/glsest)#identical...
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.524e-06 -3.464e-07 -6.245e-08 -1.691e-07  1.935e-07  6.973e-07

A high-level implementation of RLS based on lavaan functions

The following RLS code is an attempt to implement the RLS. It starts using the start=“Mplus” (default option) in lavaan.

According to help file: If “Mplus”, we use a similar scheme, but the factor loadings are estimated using the fabin3 estimator (tsls) per factor.

myrls1 <- function(model, S, nobs, max_iter=10, step_tol=1e-6, verbose=F){
  #start values a la lavaan 
  f <- sem(model, sample.cov=S, sample.nobs=nobs, do.fit=F)
  est.old <- par.df <- coef(f)
  i <- 1
  if (verbose)
    cat("\n Iteration ", i, unname(est.old), "\n")
  while(i < max_iter){
    sinv <- chol2inv(chol(lavInspect(f, "sigma.hat")))
    V <- 0.5*lav_matrix_duplication_pre_post(sinv %x% sinv)
    #update
    f <- sem(model, sample.cov=(nobs-1)*S/nobs, sample.nobs = nobs, se="none", 
             estimator="WLS", WLS.V = V)
    est.new <- coef(f)
    par.df <- rbind(par.df,est.new )
    rmse <- sqrt(mean((est.new - est.old)*(est.new - est.old)))
    if (verbose)
      cat("\n Iteration ", i, unname(est.new), "\n Rmse", rmse, "\n")
    if(rmse < step_tol)
      break
    i <- i+1
    est.old <- est.new
  }
  list(it = nrow(par.df), est=tail(par.df,1), it.df = par.df)
}

Running RLS on the simple one-factor model

Lets try the function and check whether RLS reaches the MLE for the simple model.

rlsest <- myrls1(Sim.model, Sim.cov, Sim.nobs, verbose=T)
## 
##  Iteration  1 0.9033813 1.024973 0.9125483 0.8406379 0.8779691 0.8972035 0.8514009 0.05 
## 
##  Iteration  1 0.9122521 1.039574 0.9261858 0.8789414 1.088146 0.9274464 1.014627 0.8022999 
##  Rmse 0.2827613 
## 
##  Iteration  2 0.9208464 1.0431 0.9375226 0.8887188 1.082394 0.9326856 1.005504 0.792909 
##  Rmse 0.008233801 
## 
##  Iteration  3 0.9208874 1.04314 0.9374205 0.8883194 1.083489 0.9315404 1.005982 0.7929615 
##  Rmse 0.0006035661 
## 
##  Iteration  4 0.9209628 1.043235 0.9375882 0.8884456 1.083465 0.9315597 1.005854 0.7928267 
##  Rmse 0.0001084867 
## 
##  Iteration  5 0.920964 1.043232 0.9375832 0.8884465 1.083481 0.931543 1.005855 0.7928297 
##  Rmse 8.32127e-06 
## 
##  Iteration  6 0.9209646 1.043235 0.9375867 0.8884487 1.083481 0.9315423 1.005853 0.7928269 
##  Rmse 2.193564e-06 
## 
##  Iteration  7 0.9209646 1.043235 0.9375864 0.8884488 1.083482 0.931542 1.005852 0.7928271 
##  Rmse 1.853705e-07
mlest <- coef(sem(Sim.model, sample.cov=Sim.cov, sample.nobs=Sim.nobs))
sqrt(mean((rlsest$est-mlest)^2))# yes, these agree
## [1] 7.932275e-08

This seems to work. We can also plot the iterations

library(reshape2)
library(ggplot2)
it.df <- data.frame(rlsest$it.df)
it.df$iteration <- 1:nrow(it.df)
library(reshape2)
#melt
m <- melt(it.df, id.vars="iteration")
library(ggplot2)
ggplot(m, aes(iteration, value, color=variable))+geom_line()+
  ggtitle("RLS work well in simple model")

Running RLS on Holzinger does not work?

rlsest <- myrls1(HS.model, HS.cov, HS.nobs, verbose=T)
## 
##  Iteration  1 0.7778315 1.107255 1.132894 0.9241826 1.225074 0.8544057 0.6791849 0.6908919 0.6374324 0.6753323 0.8298929 0.5981792 0.5915697 0.5109914 0.5075019 0.05 0.05 0.05 0 0 0 
## 
##  Iteration  1 0.5146489 0.6617405 1.070857 0.9413057 1.263045 1.638998 0.4613086 1.130812 0.864214 0.3588656 0.5133871 0.3198265 0.8926221 0.5909512 0.3533951 0.9072305 0.9931363 0.2531361 0.4187418 0.2376457 0.1490958 
##  Rmse 0.4082246 
## 
##  Iteration  2 0.5248951 0.705703 1.113051 0.9246252 1.133393 0.99229 0.5373452 1.145139 0.8526454 0.3741394 0.4421606 0.358756 0.803507 0.5081052 0.5416706 0.7145576 0.9532877 0.4218306 0.360285 0.2347466 0.1694364 
##  Rmse 0.165233 
## 
##  Iteration  3 0.5777643 0.7536782 1.109744 0.931888 1.263202 1.510748 0.5665621 1.127667 0.8363833 0.3713584 0.4484933 0.3541856 0.825685 0.5461965 0.4769818 0.8282346 0.9784268 0.249278 0.4192455 0.2594129 0.1508672 
##  Rmse 0.1283945 
## 
##  Iteration  4 0.5050088 0.6743521 1.113797 0.9187863 1.104491 0.8634221 0.5195844 1.141307 0.8549136 0.3705126 0.4438046 0.3585211 0.768045 0.468026 0.5604184 0.72531 0.9662647 0.4644334 0.3659063 0.2144678 0.1681582 
##  Rmse 0.1599649 
## 
##  Iteration  5 0.5855293 0.7637586 1.111174 0.9351223 1.295947 1.619136 0.5725493 1.125505 0.832913 0.3717907 0.4488599 0.3534792 0.8022594 0.5361053 0.4659274 0.8245062 0.9711742 0.2206314 0.4182098 0.2522626 0.1411216 
##  Rmse 0.1845242 
## 
##  Iteration  6 0.4869456 0.6534225 1.113642 0.9167126 1.085874 0.8222766 0.5083869 1.14309 0.8580487 0.3701577 0.4435005 0.3589822 0.7529042 0.4624614 0.5572561 0.7118003 0.9652715 0.4791442 0.3600027 0.2014105 0.1654417 
##  Rmse 0.1962511 
## 
##  Iteration  7 0.5891893 0.7681881 1.110724 0.9371876 1.31691 1.703757 0.5746619 1.124456 0.8313033 0.3718117 0.449563 0.3528338 0.7952146 0.5365491 0.4526996 0.8191042 0.9664793 0.2026114 0.4158034 0.2456209 0.1347261 
##  Rmse 0.215161 
## 
##  Iteration  8 0.4761315 0.6421271 1.113275 0.9156926 1.075899 0.8081821 0.5024409 1.14402 0.8594268 0.3699692 0.4433812 0.3592027 0.7460987 0.4627635 0.5522005 0.7043187 0.9647358 0.4852746 0.3567215 0.19657 0.1647037 
##  Rmse 0.2196734 
## 
##  Iteration  9 0.5911827 0.7703658 1.110374 0.938273 1.327493 1.74891 0.5756983 1.123826 0.8303977 0.3718116 0.4499382 0.3524962 0.7929725 0.5378327 0.4455727 0.8161105 0.9639626 0.1938711 0.4143588 0.2419558 0.1314834 
##  Rmse 0.2297285
mlest <- coef(sem(HS.model, sample.cov=HS.cov, sample.nobs=HS.nobs))
sqrt(mean((rlsest$est-mlest)^2))# no!!!
## [1] 0.158429
#plotting the iterations
it.df <- data.frame(rlsest$it.df)
it.df$iteration <- 1:nrow(it.df)
m <- melt(it.df, id.vars="iteration")
library(ggplot2)
ggplot(m, aes(iteration, value, color=variable))+geom_line()+
  ggtitle("RLS on Holzinger not converging")

Adding one residual covariance leads to convergence

HS.model <- paste(HS.model, "; x7~~x8")
rlsest <- myrls1(HS.model, HS.cov, HS.nobs, verbose=T)
## 
##  Iteration  1 0.7778315 1.107255 1.132894 0.9241826 1.225074 0.8544057 0 0.6791849 0.6908919 0.6374324 0.6753323 0.8298929 0.5981792 0.5915697 0.5109914 0.5075019 0.05 0.05 0.05 0 0 0 
## 
##  Iteration  1 0.5246734 0.6676951 1.07371 0.9422703 1.402786 2.650429 0.353608 0.4734869 1.125564 0.8626207 0.361085 0.510293 0.3201126 1.049613 0.77329 0.1130184 0.8948335 0.9907701 0.1283151 0.4157596 0.1760231 0.1038762 
##  Rmse 0.5417735 
## 
##  Iteration  2 0.5770663 0.7542846 1.113009 0.9275933 1.263242 2.520627 0.35044 0.5801948 1.129038 0.8318348 0.3751887 0.4417977 0.3581534 1.035533 0.7911494 0.1010683 0.7673286 0.9688788 0.1442705 0.3938423 0.18136 0.1005904 
##  Rmse 0.06222716 
## 
##  Iteration  3 0.5740774 0.7500386 1.114308 0.9258463 1.247245 2.514276 0.3526461 0.5736877 1.122697 0.8332825 0.3719925 0.4440919 0.3568807 1.036778 0.7941729 0.08912987 0.7837257 0.979546 0.1464641 0.399415 0.1842317 0.1024946 
##  Rmse 0.006802986 
## 
##  Iteration  4 0.5759629 0.7526447 1.114671 0.9262737 1.2456 2.515626 0.352653 0.5762044 1.122606 0.8319496 0.3722633 0.4436635 0.3570584 1.036465 0.7946786 0.08798099 0.7820312 0.978173 0.1464951 0.3993903 0.1842482 0.1020765 
##  Rmse 0.001171742 
## 
##  Iteration  5 0.575474 0.7520238 1.114878 0.9261415 1.244394 2.514189 0.3527015 0.5757321 1.122503 0.8321673 0.3722423 0.4436572 0.3570645 1.03647 0.7947909 0.08759952 0.7825839 0.9785262 0.1467142 0.3995507 0.1843897 0.1022122 
##  Rmse 0.0004857048 
## 
##  Iteration  6 0.5756858 0.7522638 1.114867 0.9261972 1.244602 2.514806 0.3527077 0.5758836 1.122497 0.8320952 0.3722489 0.443663 0.3570546 1.036466 0.7948124 0.08753944 0.7824833 0.978382 0.1466524 0.3995131 0.1843596 0.1021509 
##  Rmse 0.0001656284 
## 
##  Iteration  7 0.5756008 0.7521665 1.114878 0.9261777 1.24446 2.514527 0.3527066 0.5758304 1.1225 0.8321211 0.3722475 0.443659 0.3570588 1.036463 0.7948144 0.0875358 0.7825391 0.9784329 0.1466852 0.3995456 0.1843789 0.102174 
##  Rmse 7.627899e-05 
## 
##  Iteration  8 0.5756357 0.7522047 1.114878 0.9261857 1.244507 2.514645 0.3527083 0.5758483 1.122498 0.8321129 0.3722482 0.4436606 0.357057 1.036465 0.7948163 0.08752803 0.7825204 0.9784097 0.1466728 0.3995317 0.1843716 0.1021639 
##  Rmse 3.06683e-05 
## 
##  Iteration  9 0.5756211 0.7521886 1.114878 0.9261826 1.244486 2.514595 0.3527076 0.5758412 1.122499 0.8321163 0.3722479 0.44366 0.3570577 1.036464 0.7948158 0.08753057 0.7825292 0.9784189 0.1466781 0.3995387 0.184375 0.102168 
##  Rmse 1.30516e-05
mlest <- coef(sem(HS.model, sample.cov=HS.cov, sample.nobs=HS.nobs))
sqrt(mean((rlsest$est-mlest)^2))#yes
## [1] 4.583262e-06
#plotting the iterations
it.df <- data.frame(rlsest$it.df)
it.df$iteration <- 1:nrow(it.df)
m <- melt(it.df, id.vars="iteration")
library(ggplot2)
ggplot(m, aes(iteration, value, color=variable))+geom_line()+
  ggtitle("RLS on Holzinger with residual covariance IS converging")