library(lavaan)## This is lavaan 0.6-7## lavaan is BETA software! Please report any bugs.# 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)Here is a rough but general implementation with step-size reduction
irwls2 <- function(model, S, nobs, max_iter=10, step_tol=1e-6, verbose=F){
f.new <- sem(model, sample.cov=S, sample.nobs=nobs, do.fit=F)
sinv <- chol2inv(chol(lavInspect(f.new, "sigma.hat")))
obj.new <- 10^6
est.new <- parameterEstimates(f.new)$est
cat("\n Startvalues \n")
cat(coef(f.new))
i <- 1
while(i < max_iter){
est.old <- est.new; obj.old <- obj.new; f.old <- f.new
# V has been updated
V <- 0.5*lav_matrix_duplication_pre_post(sinv %x% sinv)
f.new <- sem(model, sample.cov=S, sample.nobs = nobs, se="none",
estimator="WLS", WLS.V = V)
est.new <- parameterEstimates(f.new, remove.nonfree = T)# slow
obj.new <- fitmeasures(f.new, "fmin")
cat("\n Iteration", i,"\n ", coef(f.new), "\n objective", obj.new , "\n")
kk <- 1; alpha <- 1
while(obj.new > obj.old && kk < 7){
alpha <- 1/2^kk
cat(" !!!! Function value increases. Set alpha to", alpha, "\n")
newpars <- est.new
newpars$est <- alpha*est.new$est+(1-alpha)*coef(f.old)
cat(" new intermediate parameters for sigma_implied: ", newpars$est , "\n \n")
#calculate sigma implied from new intermediate parameters
f.new <- sem(model, sample.cov=S, sample.nobs = nobs, se="none",
estimator="WLS", WLS.V = V, start=newpars, do.fit=F)
sinv <- chol2inv(chol(lavInspect(f.new, "sigma.hat")))
#estimate model with new sigma implied
V <- 0.5*lav_matrix_duplication_pre_post(sinv %x% sinv)
f.new <- sem(model, sample.cov=S, sample.nobs = nobs, se="none",
estimator="WLS", WLS.V = V)
est.new <- parameterEstimates(f.new, remove.nonfree = T)# slow
obj.new <- fitmeasures(f.new, "fmin")
cat("which gives new estimates:",coef(f.new),
" \n and objective", obj.new , "\n \n ")
kk <- kk+1
}
if(kk==7){
cat("ERROR")
return(-999)
}
sinv <- chol2inv(chol(lavInspect(f.new, "sigma.hat")))
i <-i+1
}
return(est.new)
}rlssolution <-
irwls2(HS.model, HS.cov, HS.nobs)##
## Startvalues
## 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.5146491 0.6617407 1.070857 0.9413056 1.263047 1.638996 0.4628464 1.134582 0.8670949 0.3600619 0.5150971 0.3208933 0.895598 0.5929196 0.354575 0.9102544 0.9964464 0.2539798 0.4201376 0.2384378 0.1495928
## objective 0.2684643
##
## Iteration 2
## 0.5248952 0.7057031 1.113051 0.9246252 1.133393 0.9922897 0.5391364 1.148956 0.8554875 0.3753865 0.4436345 0.3599519 0.806185 0.5097988 0.5434764 0.7169401 0.9564655 0.4232369 0.3614863 0.2355293 0.1700012
## objective 0.1356554
##
## Iteration 3
## 0.577765 0.7536793 1.109746 0.931889 1.263206 1.510734 0.5684527 1.131424 0.8391703 0.3725968 0.4499878 0.3553662 0.8284369 0.5480141 0.4785768 0.8309933 0.9816863 0.2501108 0.4206418 0.2602782 0.1513705
## objective 0.1366088
## !!!! Function value increases. Set alpha to 0.5
## new intermediate parameters for sigma_implied: 0.5513301 0.7296912 1.111398 0.9282571 1.1983 1.251512 0.5537945 1.14019 0.8473289 0.3739917 0.4468112 0.357659 0.8173109 0.5289065 0.5110266 0.7739667 0.9690759 0.3366739 0.391064 0.2479037 0.1606858
##
## which gives new estimates: 0.5444949 0.7205679 1.110634 0.9250965 1.166704 1.050418 0.5446947 1.139711 0.8497183 0.371922 0.448142 0.3575083 0.8081248 0.4953692 0.5611594 0.7875003 0.9788764 0.398775 0.3963503 0.2522036 0.1735243
## and objective 0.1385296
##
## !!!! Function value increases. Set alpha to 0.25
## new intermediate parameters for sigma_implied: 0.5297951 0.7094193 1.112447 0.924743 1.141721 1.006822 0.540526 1.146644 0.8540452 0.3745204 0.4447614 0.359341 0.8066699 0.5061914 0.5478972 0.7345801 0.9620682 0.4171215 0.3702023 0.2396978 0.170882
##
## which gives new estimates: 0.5739973 0.7506829 1.110271 0.9310121 1.246317 1.424914 0.5659681 1.132592 0.8402332 0.3725937 0.4496245 0.3557019 0.8260706 0.5397289 0.4956257 0.8308615 0.9827929 0.2727884 0.4204203 0.2652458 0.1573354
## and objective 0.1364524
##
## !!!! Function value increases. Set alpha to 0.125
## new intermediate parameters for sigma_implied: 0.5310329 0.7113256 1.112704 0.9254235 1.147508 1.046368 0.5424904 1.14691 0.8535807 0.3750374 0.4443833 0.3594206 0.8086707 0.5135401 0.5374951 0.7311802 0.9597564 0.4044309 0.368853 0.2392438 0.168418
##
## which gives new estimates: 0.5713767 0.7483073 1.109819 0.9305206 1.237555 1.387522 0.56398 1.133376 0.8410808 0.3724505 0.4496778 0.355783 0.8294257 0.5394345 0.5016867 0.8280803 0.983084 0.2841346 0.4187589 0.2664343 0.1601404
## and objective 0.1377856
##
## !!!! Function value increases. Set alpha to 0.0625
## new intermediate parameters for sigma_implied: 0.5278003 0.7083658 1.112849 0.9249936 1.139903 1.016992 0.5406891 1.147982 0.8545871 0.375203 0.4440122 0.3596913 0.8076375 0.5116511 0.5408646 0.7238863 0.9581292 0.414543 0.3650658 0.2374608 0.1693849
##
## which gives new estimates: 0.574971 0.7513716 1.1098 0.9312702 1.25127 1.454013 0.5665131 1.132304 0.8399895 0.3725369 0.4498389 0.3555585 0.8293962 0.5445903 0.4890286 0.8301533 0.982448 0.2651739 0.4200304 0.2635233 0.1554962
## and objective 0.1371633
##
## !!!! Function value increases. Set alpha to 0.03125
## new intermediate parameters for sigma_implied: 0.5264601 0.7071302 1.112949 0.9248328 1.137076 1.006719 0.5399919 1.148435 0.8550032 0.3752975 0.4438284 0.3598146 0.8069103 0.5108861 0.5417749 0.720478 0.9572775 0.4182975 0.3633158 0.2364041 0.1695479
##
## which gives new estimates: 0.5762528 0.7524394 1.109759 0.9315555 1.256753 1.480004 0.567404 1.131899 0.8396146 0.3725615 0.4499135 0.3554678 0.8291282 0.546315 0.4842556 0.8305526 0.9821126 0.258145 0.4203182 0.2620893 0.1536198
## and objective 0.1369591
##
## !!!! Function value increases. Set alpha to 0.015625
## new intermediate parameters for sigma_implied: 0.5256976 0.7064333 1.113 0.9247335 1.13532 0.9999102 0.5395781 1.148689 0.8552395 0.3753424 0.4437326 0.3598818 0.8065435 0.5103694 0.5425511 0.7187153 0.9568662 0.4206574 0.3624055 0.2359443 0.1697453
##
## which gives new estimates: 0.5769932 0.7530478 1.109749 0.9317182 1.259908 1.495 0.5679167 1.131668 0.8393983 0.3725783 0.4499508 0.3554179 0.8288339 0.5471841 0.4814887 0.8307792 0.9819101 0.2541859 0.4204815 0.2612213 0.1525251
## and objective 0.1367995
##
## ERROR