Error Estimate in Summary is wrong, so is $Hessian

name matrix row col parameter estimate error estimate
1 A 1 6 0.39715247 163.52372
2 A 2 6 0.50366153 149.57542
3 A 3 6 0.57724183 141.16243
4 A 4 6 0.70277427 125.43848
5 A 5 6 0.79625063 53.94078
6 S 1 1 0.04081422 510.59245
7 S 2 2 0.03802001 510.18074
8 S 3 3 0.04082720 460.75238
9 S 4 4 0.03938708 428.60445
10 S 5 5 0.03628711 387.12985

and these values come from $hessian
$hessian

163.5237 -18.68289 -23.98918 -25.72917 -28.27548 52.367788 7.650826 11.152047 18.89175 -17.21338
0.0000 149.57542 -23.61426 -31.53574 -43.89935 -8.914412 33.044284 -6.495349 -18.81572 -24.23432
0.0000 0.00000 141.16243 -42.41420 -56.42054 -14.125920 -12.563846 31.168712 -31.56789 -36.47626
0.0000 0.00000 0.00000 125.43848 -99.01438 -10.837646 -13.872889 4.040232 45.16674 -50.46471
0.0000 0.00000 0.00000 0.00000 53.94078 -12.566096 -9.675148 -15.248109 -19.64663 18.52240
0.0000 0.00000 0.00000 0.00000 0.00000 510.592447 17.780684 7.882866 26.36950 27.23222
0.0000 0.00000 0.00000 0.00000 0.00000 0.000000 510.180737 11.295398 45.52152 30.39196
0.0000 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 460.752382 40.33177 63.59937
0.0000 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 0.000000 428.60445 43.38175
0.0000 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 0.000000 0.00000 387.12985

What is printed as $Hessian is the Cholesky decomposition of it, so what should really go in the Error Estimate column is

sqrt(diag(solve(t(factorresult@output$hessian) %*% factorresult@output$hessian )))

0.010674776 0.012604163 0.014016030 0.016688139 0.018605535 0.001966411 0.001975071 0.002205340 0.002347757 0.002583113

This jives reasonably well with bootstrap estimates (thx to the handy cov.wt function!):

mles<-function(dataset,wt){
manifests <- names(dataset)
latents <- c("G")
covwt<-cov.wt(dataset,wt)
mlevals<-mxRun(mxModel("One Factor", type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from=latents, to=manifests),
mxPath(from=manifests, arrows=2),
mxPath(from=latents, arrows=2,
free=F, values=1.0),
mxData(covwt$cov, type="cov",
numObs=500)))
return(as.vector(mlevals@output$estimate))}
boot(demoOneFactor,mles,R=1000)

Bootstrap Statistics :
original bias std. error
t1* 0.38841392 0.0084411481 0.008744156
t2* 0.49523971 0.0081500273 0.010768885
t3* 0.57578325 0.0008155712 0.012274162
t4* 0.69536404 0.0070740688 0.013864563
t5* 0.78976444 0.0059875654 0.014175315
t6* 0.03943445 0.0013193612 0.001660683
t7* 0.04021663 -0.0022253005 0.001440713
t8* 0.04123489 -0.0004473884 0.001870993
t9* 0.03662419 0.0027992703 0.002094885
t10* 0.03458172 0.0017000842 0.001910026

Accordingly, where $hessian is now in the output, we should print
t(factorresult@output$hessian) %*% factorresult@output$hessian
and
sqrt(diag(solve(t(factorresult@output$hessian) %*% factorresult@output$hessian )))
should go in the Error Estimate column

This may get messed up in the presence of non-linear constraints!